source: trunk/src/read_civxdata.m @ 239

Last change on this file since 239 was 236, checked in by sommeria, 13 years ago

correct Matlab PIV, remove call to image tool box. Improve menu of uvmat VelType? (replacement of buttons)

File size: 11.0 KB
Line 
1%'read_civxdata': reads civx data from netcdf files
2%------------------------------------------------------------------
3%
4% function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)
5%
6% OUTPUT:
7% Field: structure representing the selected field, containing
8%            .Txt: (char string) error message if any
9%            .ListGlobalAttribute: list of global attributes containing:
10%                    .NbCoord: number of vector components
11%                    .NbDim: number of dimensions (=2 or 3)
12%                    .dt: time interval for the corresponding image pair
13%                    .Time: absolute time (average of the initial image pair)
14%                    .CivStage: =0,
15%                       =1, civ1 has been performed only
16%                       =2, fix1 has been performed
17%                       =3, pacth1 has been performed
18%                       =4, civ2 has been performed
19%                       =5, fix2 has been performed
20%                       =6, pacth2 has been performed
21%                     .CoordUnit: 'pixel'
22%             .ListVarName: {'X'  'Y'  'U'  'V'  'F'  'FF'}
23%                   .X, .Y, .Z: set of vector coordinates
24%                    .U,.V,.W: corresponding set of vector components
25%                    .F: warning flags
26%                    .FF: error flag, =0 for good vectors
27%                     .C: scalar associated with velocity (used for vector colors)
28%                    .DijU; matrix of spatial derivatives (DijU(1,1,:)=DUDX,
29%                        DijU(1,2,:)=DUDY, Dij(2,1,:)=DVDX, DijU(2,2,:)=DVDY
30%
31% VelTypeOut: velocity type corresponding to the selected field: ='civ1','interp1','interp2','civ2'....
32%
33% INPUT:
34% filename: file name (string).
35% FieldNames =cell of field names to get, which can contain the strings:
36%             'ima_cor': image correlation, vec_c or vec2_C
37%             'vort','div','strain': requires velocity derivatives DUDX...
38%             'error': error estimate (vec_E or vec2_E)
39%             
40% VelType : character string indicating the types of velocity fields to read ('civ1','civ2'...)
41%            if vel_type=[] or'*', a  priority choice, given by vel_type_out{1,2}, is done depending
42%            if vel_type='filter'; a structured field is sought (filter2 in priority, then filter1)
43
44
45% FUNCTIONS called:
46% 'varcivx_generator':, sets the names of vaiables to read in the netcdf file
47% 'nc2struct': reads a netcdf file
48
49function [Field,VelTypeOut,errormsg]=read_civxdata(filename,FieldNames,VelType)
50errormsg='';
51DataTest=nc2struct(filename,'ListGlobalAttribute','Conventions');
52if isequal(DataTest.Conventions,'uvmat/civdata')%test for new civ format
53     [Field,VelTypeOut]=read_civdata(filename,FieldNames,VelType);
54     return
55end
56   
57%% default input
58if ~exist('VelType','var')
59    VelType=[];
60end
61if isequal(VelType,'*')
62    VelType=[];
63end
64if ~exist('FieldNames','var')
65    FieldNames=[]; %default
66end
67
68%% reading data
69VelTypeOut=VelType;%default
70[var,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);%determine the names of constants and variables to read
71[Field,vardetect,ichoice]=nc2struct(filename,var);%read the variables in the netcdf file
72if isfield(Field,'Txt')
73    errormsg=Field.Txt;
74    return
75end
76if vardetect(1)==0
77     errormsg=[ 'requested field not available in ' filename '/' VelType];
78     return
79end
80var_ind=find(vardetect);
81for ivar=1:length(var_ind)
82    Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
83    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
84end
85VelTypeOut=VelType;
86if ~isempty(ichoice)
87    VelTypeOut=vel_type_out_cell{ichoice};
88end
89
90%% adjust for Djui:
91if isfield(Field,'DjUi')
92    Field.ListVarName(end-2:end)=[];
93    Field.ListVarName{end}='DjUi';
94    Field.VarDimName(end-2:end)=[];
95    Field.VarAttribute(end-2:end)=[];
96end
97
98%% renaming for standard conventions
99Field.NbCoord=Field.nb_coord;
100Field.NbDim=Field.nb_dim;
101
102%% CivStage
103if isfield(Field,'patch2')&& isequal(Field.patch2,1)
104    Field.CivStage=6;
105elseif isfield(Field,'fix2')&& isequal(Field.fix2,1)
106    Field.CivStage=5;
107elseif isfield(Field,'civ2')&& isequal(Field.civ2,1)
108    Field.CivStage=4;
109elseif isfield(Field,'patch')&& isequal(Field.patch,1)
110    Field.CivStage=3;
111elseif isfield(Field,'fix')&& isequal(Field.fix,1)
112    Field.CivStage=2;
113else
114    Field.CivStage=1;
115end
116
117%determine the appropriate constant for time and dt for the PIV pair
118test_civ1=isequal(VelTypeOut,'civ1')||isequal(VelTypeOut,'interp1')||isequal(VelTypeOut,'filter1');
119test_civ2=isequal(VelTypeOut,'civ2')||isequal(VelTypeOut,'interp2')||isequal(VelTypeOut,'filter2');
120Field.Time=0; %default
121Field.TimeUnit='s';
122if test_civ1
123    if isfield(Field,'absolut_time_T0')
124        Field.Time=double(Field.absolut_time_T0);
125        Field.dt=double(Field.dt);
126    else
127       errormsg='the input file is not civx';
128       Field.CivStage=0;
129       Field.dt=0;
130    end
131elseif test_civ2
132    Field.Time=double(Field.absolut_time_T0_2);
133    Field.dt=double(Field.dt2);
134else
135    errormsg='the input file is not civx';
136    Field.CivStage=0;
137    Field.dt=0;
138end
139
140
141
142%% rescale fields to pixel coordinates
143if isfield(Field,'pixcmx')
144    Field.pixcmx=double(Field.pixcmx);
145    Field.pixcmy=double(Field.pixcmy);
146    Field.U=Field.U*Field.pixcmx;
147    Field.V=Field.V*Field.pixcmy;
148    Field.X=Field.X*Field.pixcmx;
149    Field.Y=Field.Y*Field.pixcmy;
150end
151if ~isequal(Field.dt,0)
152    Field.U=Field.U*Field.dt;%translate in px displacement
153    Field.V=Field.V*Field.dt;
154    if isfield(Field,'DjUi')
155       Field.DjUi(:,1,1)=Field.dt*Field.DjUi(:,1,1);
156       Field.DjUi(:,2,2)=Field.dt*Field.DjUi(:,2,2);
157       Field.DjUi(:,1,2)=(Field.pixcmy/Field.pixcmx)*Field.dt*Field.DjUi(:,1,2);
158       Field.DjUi(:,2,1)=(Field.pixcmx/Field.pixcmy)*Field.dt*Field.DjUi(:,2,1);
159    end
160end
161
162%% update list of global attributes
163List=Field.ListGlobalAttribute;
164ind_remove=[];
165for ilist=1:length(List)
166    switch(List{ilist})
167        case {'patch2','fix2','civ2','patch','fix','dt2','absolut_time_T0','absolut_time_T0_2','nb_coord','nb_dim','pixcmx','pixcmy'}
168            ind_remove=[ind_remove ilist];
169            Field=rmfield(Field,List{ilist});
170    end
171end
172List(ind_remove)=[];
173Field.ListGlobalAttribute=[{'NbCoord'},{'NbDim'} List {'Time','TimeUnit','CivStage','CoordUnit'}];
174Field.CoordUnit='pixel';
175
176
177%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
178% [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
179%INPUT:
180% FieldNames =cell of field names to get, which can contain the strings:
181%             'ima_cor': image correlation, vec_c or vec2_C
182%             'vort','div','strain': requires velocity derivatives DUDX...
183%             'error': error estimate (vec_E or vec2_E)
184%             
185% vel_type: character string indicating the types of velocity fields to read ('civ1','civ2'...)
186%            if vel_type=[] or'*', a  priority choice, given by vel_type_out{1,2}, is done depending
187%            if vel_type='filter'; a structured field is sought (filter2 in priority, then filter1)
188
189function [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
190
191%% default input values
192if ~exist('vel_type','var'),vel_type=[];end;
193if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
194if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar
195if ischar(FieldNames), FieldNames={FieldNames}; end;
196
197%% select the priority order for automatic vel_type selection
198testder=0;
199for ilist=1:length(FieldNames)
200    if ~isempty(FieldNames{ilist})
201    switch FieldNames{ilist}
202        case {'vort','div','strain'}
203            testder=1;
204    end
205    end
206end     
207if isempty(vel_type) || isequal(vel_type,'*') %undefined velocity type (civ1,civ2...)
208    if testder
209         vel_type_out{1}='filter2'; %priority to filter2 for scalar reading, filter1 as second
210        vel_type_out{2}='filter1';
211    else
212        vel_type_out{1}='civ2'; %priority to civ2 for vector reading, civ1 as second priority     
213        vel_type_out{2}='civ1';
214    end
215elseif isequal(vel_type,'filter')
216        vel_type_out{1}='filter2'; %priority to filter2 for scalar reading, filter1 as second
217        vel_type_out{2}='filter1';
218        if ~testder
219            vel_type_out{3}='civ1';%civ1 as third priority if derivatives are not needed
220        end
221elseif testder
222    test_civ1=isequal(vel_type,'civ1')||isequal(vel_type,'interp1')||isequal(vel_type,'filter1');
223    if test_civ1
224        vel_type_out{1}='filter1'; %switch to filter for reading spatial derivatives
225    else
226        vel_type_out{1}='filter2';
227    end
228else   
229    vel_type_out{1}=vel_type;%imposed velocity field
230end
231vel_type_out=vel_type_out';
232
233%% determine names of netcdf variables to read
234var={'X','Y','Z','U','V','W','C','F','FF'};
235role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
236units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
237if testder
238    var=[var {'DjUi(:,1,1)','DjUi(:,1,2)','DjUi(:,2,1)','DjUi(:,2,2)'}];
239    role=[role {'tensor','tensor','tensor','tensor'}];
240    units=[units {'pixel','pixel','pixel','pixel'}];
241end
242for ilist=1:length(vel_type_out)
243    var=[var;varname1(vel_type_out{ilist},FieldNames)];
244end
245
246%------------------------------------------------------------------------ 
247%--- determine  var names to read
248function varin=varname1(vel_type,FieldNames)
249%------------------------------------------------------------------------
250testder=0;
251C1='';
252C2='';
253for ilist=1:length(FieldNames)
254    if ~isempty(FieldNames{ilist})
255    switch FieldNames{ilist}
256        case 'ima_cor' %image correlation corresponding to a vel vector
257            C1='vec_C';
258            C2='vec2_C';
259        case 'error'
260            C1='vec_E';
261            C2='vec2_E';
262        case {'vort','div','strain'}
263            testder=1;
264    end
265    end
266end     
267switch vel_type
268    case 'civ1'
269        varin={'vec_X','vec_Y','vec_Z','vec_U','vec_V','vec_W',C1,'vec_F','vec_FixFlag'};
270    case 'interp1'
271        varin={'vec_patch_X','vec_patch_Y','','vec_patch0_U','vec_patch0_V','','','',''};
272    case 'filter1'
273        varin={'vec_patch_X','vec_patch_Y','','vec_patch_U','vec_patch_V','','','',''};
274    case 'civ2'
275        varin={'vec2_X','vec2_Y','vec2_Z','vec2_U','vec2_V','vec2_W',C2,'vec2_F','vec2_FixFlag'};
276    case 'interp2'
277        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch0_U','vec2_patch0_V','vec2_patch0_W','','',''};
278    case 'filter2'
279        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch_U','vec2_patch_V','vec2_patch0_W','','',''};
280end
281if testder
282     switch vel_type
283        case 'filter1'
284            varin=[varin {'vec_patch_DUDX','vec_patch_DVDX','vec_patch_DUDY','vec_patch_DVDY'}];
285        case 'filter2'
286            varin=[varin {'vec2_patch_DUDX','vec2_patch_DVDX','vec2_patch_DUDY','vec2_patch_DVDY'}];
287    end   
288end
Note: See TracBrowser for help on using the repository browser.