source: trunk/src/read_civxdata.m @ 106

Last change on this file since 106 was 93, checked in by sommeria, 14 years ago

FFT: improved to deal with NaN data
merge_proj: corrected to reproduce dt if unique
uvmat.fig: tooltip corrected
plot_field: bug for isocontour corrected,
im_filter: cleaning
phys_polar: spatial derivative included (still to check)
set_obeject.fig: minor correction
struct2nc: comments improved
uvmat: button NB implemented
read_civxdata: error message improved

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