Index: /trunk/src/calc_field.m =================================================================== --- /trunk/src/calc_field.m (revision 514) +++ /trunk/src/calc_field.m (revision 515) @@ -38,12 +38,12 @@ %% list of field options implemented -FieldOptions={'velocity';...%image correlation corresponding to a vel vector - 'ima_cor';...%image correlation corresponding to a vel vector - 'norm_vel';...%norm of the velocity - 'vort';...%vorticity - 'div';...%divergence - 'strain';...%rate of strain - 'u';... %u velocity component - 'v';... %v velocity component +FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector + 'C';...%image correlation corresponding to a vel vector + 'norm(U,V)';...%norm of the velocity + 'curl(U,V)';...%vorticity + 'div(U,V)';...%divergence + 'strain(U,V)';...%rate of strain + 'U';... %u velocity component + 'V';... %v velocity component 'w';... %w velocity component 'w_normal';... %w velocity component normal to the plane @@ -77,8 +77,4 @@ end FieldList=FieldList(check_calc==1); -% if isempty(FieldList) -% DataOut=DataIn; -% return -% end if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X)) nbcoord=3; @@ -91,18 +87,35 @@ units_cell={}; +%% reproduce global attributes +DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; +for ilist=1:numel(DataOut.ListGlobalAttribute) + DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist}); +end + %% interpolation with new civ data -if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der) - %reproduce global attributes - DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; - for ilist=1:numel(DataOut.ListGlobalAttribute) - DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist}); - end - +[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn); +icell_tps=[]; +for icell=1:numel(CellVarIndex) + VarType=VarTypeCell{icell}; + if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps) + icell_tps=[icell_tps icell]; + end +end + +%if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der) +if ~isempty(icell_tps) %create a default grid if needed if ~exist('Coord_interp','var') - XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates - YMax=max(max(DataIn.SubRange(2,:,:))); - XMin=min(min(DataIn.SubRange(1,:,:))); - YMin=min(min(DataIn.SubRange(2,:,:))); + for ind=1:numel(icell_tps) + SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps}); + XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates + YMax(ind)=max(max(SubRange(2,:,:))); + XMin(ind)=min(min(SubRange(1,:,:))); + YMin(ind)=min(min(SubRange(2,:,:))); + end + XMax=max(XMax); + YMax=max(YMax); + XMin=min(XMin); + YMin=min(YMin); if ~isfield(DataIn,'Mesh') DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps)); @@ -119,11 +132,7 @@ coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to visualisation coord_y=YMin:DataIn.Mesh:YMax; -% npx=length(coord_x); -% npy=length(coord_y); DataOut.coord_x=[XMin XMax]; DataOut.coord_y=[YMin YMax]; [XI,YI]=meshgrid(coord_x,coord_y); -% XI=reshape(XI,[],1); -% YI=reshape(YI,[],1); Coord_interp=cat(3,XI,YI);%[XI YI]; end @@ -131,5 +140,5 @@ npy=size(Coord_interp,1); Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3)); -% npy=length(coord_y); + %initialise output nb_sites=size(Coord_interp,1); @@ -152,73 +161,76 @@ % interpolate data in each subdomain - for isub=1:NbSubDomain - nbvec_sub=DataIn.NbSites(isub); - check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)'); - ind_sel=find(sum(check_range,2)==nb_coord); - %rho smoothing parameter - % epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites - % ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs - nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap) - if check_grid - EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources' - end - if check_der - [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources' - end - ListVar={}; - for ilist=1:length(FieldList) - var_count=numel(ListVar); - switch FieldList{ilist} - case 'velocity' - ListVar=[ListVar {'U', 'V'}]; - VarAttribute{var_count+1}.Role='vector_x'; - VarAttribute{var_count+2}.Role='vector_y'; - DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); - DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); - case 'u' - ListVar=[ListVar {'u'}]; - VarAttribute{var_count+1}.Role='scalar'; - DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); - case 'v' - ListVar=[ListVar {'v'}]; - VarAttribute{var_count+1}.Role='scalar'; - DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); - case 'norm_vel' - ListVar=[ListVar {'norm_vel'}]; - VarAttribute{var_count+1}.Role='scalar'; - U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub); - V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub); - DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V); - case 'vort' - ListVar=[ListVar {'vort'}]; - VarAttribute{var_count+1}.Role='scalar'; - DataOut.vort(ind_sel)=DataOut.vort(ind_sel)-EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub); - case 'div' - ListVar=[ListVar {'div'}]; - VarAttribute{var_count+1}.Role='scalar'; - DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub); - case 'strain' - ListVar=[ListVar {'strain'}]; - VarAttribute{var_count+1}.Role='scalar'; - DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub); + for icell=icell_tps + for isub=1:NbSubDomain + nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps}); + SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps}); + check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)'); + ind_sel=find(sum(check_range,2)==nb_coord); + nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap) + Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps}); + if check_grid + EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources' end - end - end - DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang - nbval(nbval==0)=1;% to avoid division by zero for averaging - if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed - DataOut.ListVarName=[DataOut.ListVarName {'FF'}]; - DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}]; - DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag'; - end - DataOut.ListVarName=[DataOut.ListVarName ListVar]; - for ifield=1:numel(ListVar) - VarDimName{ifield}={'coord_y','coord_x'}; - DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval; - DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx); - end - DataOut.FF=reshape(DataOut.FF,npy,npx); - DataOut.VarDimName=[DataOut.VarDimName VarDimName]; - DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute]; + if check_der + [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources' + end + ListVar={}; + U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps}); + V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps}); + for ilist=1:length(FieldList) + var_count=numel(ListVar); + switch FieldList{ilist} + case 'velocity' + ListVar=[ListVar {'U', 'V'}]; + VarAttribute{var_count+1}.Role='vector_x'; + VarAttribute{var_count+2}.Role='vector_y'; + DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub); + DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub); + case 'u' + ListVar=[ListVar {'u'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub); + case 'v' + ListVar=[ListVar {'v'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub); + case 'norm_vel' + ListVar=[ListVar {'norm_vel'}]; + VarAttribute{var_count+1}.Role='scalar'; + U=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub); + V=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub); + DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V); + case 'vort' + ListVar=[ListVar {'vort'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.vort(ind_sel)=DataOut.vort(ind_sel)-EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub); + case 'div' + ListVar=[ListVar {'div'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub); + case 'strain' + ListVar=[ListVar {'strain'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub); + end + end + end + DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang + nbval(nbval==0)=1;% to avoid division by zero for averaging + if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed + DataOut.ListVarName=[DataOut.ListVarName {'FF'}]; + DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}]; + DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag'; + end + DataOut.ListVarName=[DataOut.ListVarName ListVar]; + for ifield=1:numel(ListVar) + VarDimName{ifield}={'coord_y','coord_x'}; + DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval; + DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx); + end + DataOut.FF=reshape(DataOut.FF,npy,npx); + DataOut.VarDimName=[DataOut.VarDimName VarDimName]; + DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute]; + end else Index: /trunk/src/calc_field_interp.m =================================================================== --- /trunk/src/calc_field_interp.m (revision 515) +++ /trunk/src/calc_field_interp.m (revision 515) @@ -0,0 +1,88 @@ + +%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them +%--------------------------------------------------------------------- +% [DataOut,VarAttribute,errormsg]=calc_field_interp(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp) +% +% OUTPUT: +% DataOut: structure representing the output fields +% +% INPUT: +% Coord_tps: +% NbSites +% SubRange +% FieldVar +% Operation: cell array representing the list of operations (eg div, rot..) +% Coord_interp: coordiantes of sites on which the fields need to be calculated + +function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,FieldVar,Operation,XI,YI) + +%% nbre of subdomains +% if ndims(Coord_interp)==3 +% nb_coord=size(Coord_interp,3); +% npx=size(Coord_interp,2); +% npy=size(Coord_interp,1); +% nb_sites=npx*npy; +% Coord_interp=reshape(Coord_interp,nb_sites,nb_coord); +% else +% nb_coord=size(Coord_interp,2); +% nb_sites=size(Coord_interp,1); +% end +VarVal=[]; +ListVarName={}; +VarAttribute={}; +errormsg=''; +check_u=0; +check_v=0; +for ilist=1:length(Operation) + switch Operation{ilist} + case {'U'} + check_u=1; + case {'V'} + check_v=1; + case {'vec(U,V)','norm(U,V)'} + check_u=1; + check_v=1; + end +end +if check_u + F_u = TriScatteredInterp(Coord,FieldVar(:,1),'linear'); +end +if check_v + F_v = TriScatteredInterp(Coord,FieldVar(:,2),'linear'); +end +for ilist=1:length(Operation) + nbvar=numel(ListVarName); + switch Operation{ilist} + case 'vec(U,V)' + VarVal{nbvar+1}=F_u(XI,YI); + VarVal{nbvar+2}=F_v(XI,YI); + ListVarName{nbvar+1}='U'; + ListVarName{nbvar+2}='V'; + VarAttribute{nbvar+1}.Role='vector_x'; + VarAttribute{nbvar+2}.Role='vector_y'; + case 'U' + VarVal{nbvar+1}=F_u(XI,YI); + ListVarName{nbvar+1}='U'; + VarAttribute{nbvar+1}.Role='scalar'; + case 'V' + VarVal{nbvar+1}=F_v(XI,YI); + ListVarName{nbvar+1}='V'; + VarAttribute{nbvar+1}.Role='scalar'; + case 'norm(U,V)' + VarVal{nbvar+1}=sqrt(F_u(XI,YI).*F_u(XI,YI)+F_v(XI,YI).*F_v(XI,YI)); + ListVarName{nbvar+1}='norm(U,V)'; + VarAttribute{nbvar+1}.Role='scalar'; + end +end +nbvar=numel(ListVarName); +ListVarName{nbvar+1}='FF'; +VarVal{nbvar+1}=isnan(VarVal{nbvar}); +VarAttribute{nbvar+1}.Role='errorflag'; + +% Attr_FF.Role='errorflag'; +% VarAttribute=[VarAttribute {Attr_FF}]; + + + + + Index: /trunk/src/calc_field_tps.m =================================================================== --- /trunk/src/calc_field_tps.m (revision 515) +++ /trunk/src/calc_field_tps.m (revision 515) @@ -0,0 +1,128 @@ + +%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them +%--------------------------------------------------------------------- +% [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp) +% +% OUTPUT: +% DataOut: structure representing the output fields +% +% INPUT: +% Coord_tps: +% NbSites +% SubRange +% FieldVar +% Operation: cell array representing the list of operations (eg div, rot..) +% Coord_interp: coordiantes of sites on which the fields need to be calculated + +function [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp) + +%list of defined scalars to display in menus (in addition to 'ima_cor'). +% a type is associated to each scalar: +% 'discrete': related to the individual velocity vectors, not interpolated by patch +% 'vel': calculated from velocity components, continuous field (interpolated with velocity) +% 'der': needs spatial derivatives +% 'var': the scalar name corresponds to a field name in the netcdf files +% a specific variable name for civ1 and civ2 fields are also associated, if +% the scalar is calculated from other fields, as explicited below +errormsg=''; + +%% nbre of subdomains +if ndims(Coord_interp)==3 + nb_coord=size(Coord_interp,3); + npx=size(Coord_interp,2); + npy=size(Coord_interp,1); + nb_sites=npx*npy; + Coord_interp=reshape(Coord_interp,nb_sites,nb_coord); +else + nb_coord=size(Coord_interp,2); + nb_sites=size(Coord_interp,1); +end +NbSubDomain=size(Coord_tps,3); +nbval=zeros(nb_sites,1); + +%% list of operations +check_grid=0; +check_der=0; +for ilist=1:length(Operation) + OperationType=regexprep(Operation{ilist},'(.+',''); + switch OperationType + case 'vec' + DataOut.U=zeros(nb_sites,1); + DataOut.V=zeros(nb_sites,1); + VarAttribute{1}.Role='vector_x'; + VarAttribute{2}.Role='vector_y'; + case {'U','V','norm'} + check_grid=1; + DataOut.(OperationType)=zeros(nb_sites,1); + VarAttribute{1}.Role='scalar'; + case {'curl','div','strain'} + check_der=1; + DataOut.(OperationType)=zeros(nb_sites,1); + VarAttribute{1}.Role='scalar'; + end +end +Attr_FF.Role='errorflag'; +VarAttribute=[VarAttribute {Attr_FF}]; + +%% loop on subdomains +for isub=1:NbSubDomain + nbvec_sub=NbSites(isub); + check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)'); + ind_sel=find(sum(check_range,2)==nb_coord); + nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap) + if check_grid + EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources' + end + if check_der + [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources' + end + ListVar={}; + for ilist=1:length(Operation) + var_count=numel(ListVar); + switch Operation{ilist} + case 'vec(U,V)' + ListVar=[ListVar {'U', 'V'}]; + VarAttribute{var_count+1}.Role='vector_x'; + VarAttribute{var_count+2}.Role='vector_y'; + DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1); + DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2); + case 'U' + ListVar=[ListVar {'U'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1); + case 'V' + ListVar=[ListVar {'V'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2); + case 'norm(U,V)' + ListVar=[ListVar {'norm'}]; + VarAttribute{var_count+1}.Role='scalar'; + U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1); + V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2); + DataOut.norm(ind_sel)=sqrt(U.*U+V.*V); + case 'curl(U,V)' + ListVar=[ListVar {'curl'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2); + case 'div(U,V)' + ListVar=[ListVar {'div'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2); + case 'strain(U,V)' + ListVar=[ListVar {'strain'}]; + VarAttribute{var_count+1}.Role='scalar'; + DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2); + end + end +end +DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang +nbval(nbval==0)=1;% to avoid division by zero for averaging +ListFieldOut=fieldnames(DataOut); +for ifield=1:numel(ListFieldOut) + DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval; + DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx); +end + + + + Index: /trunk/src/calc_tps.m =================================================================== --- /trunk/src/calc_tps.m (revision 514) +++ /trunk/src/calc_tps.m (revision 515) @@ -1,3 +1,3 @@ -function DataOut=calc_tps(DataIn) +function DataOut=calc_tps(DataIn,checkall) DataOut=DataIn;%default SubDomain=1000; %default, estimated nbre of vectors in a subdomain used for tps @@ -5,17 +5,63 @@ SubDomain=DataIn.SubDomain;% end -[DataOut.SubRange,DataOut.NbSites,DataOut.Coord_tps,DataOut.U_tps,DataOut.V_tps] =... - filter_tps([DataIn.X(DataIn.FF==0) DataIn.Y(DataIn.FF==0)],DataIn.U(DataIn.FF==0),DataIn.V(DataIn.FF==0),[],SubDomain,0); -nbvar=numel(DataIn.ListVarName); -DataOut.ListVarName=[DataIn.ListVarName {'SubRange','NbSites','Coord_tps','U_tps','V_tps'}]; -DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},... - {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}]; -DataOut.VarAttribute{nbvar+3}.Role='coord_tps'; -DataOut.VarAttribute{nbvar+4}.Role='vector_x'; -DataOut.VarAttribute{nbvar+5}.Role='vector_y'; -if isfield(DataOut,'ListDimName')%cleaning - DataOut=rmfield(DataOut,'ListDimName'); +[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn); +nbtps=0; +for icell=1:numel(CellVarIndex); + VarType=VarTypeCell{icell}; + if NbDimVec(icell)>=2 && ~isempty(VarType.coord_x) + nbtps=nbtps+1; + X=DataIn.(DataIn.ListVarName{VarType.coord_x}); + Y=DataIn.(DataIn.ListVarName{VarType.coord_y}); + if ~isempty(VarType.vector_x)&&~isempty(VarType.vector_y) + Attr=DataIn.VarAttribute{VarType.vector_x}; + if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'FieldRequest')&&strcmp(Attr.FieldRequest,'derivatives'))) + U=DataIn.(DataIn.ListVarName{VarType.vector_x}); + V=DataIn.(DataIn.ListVarName{VarType.vector_y}); + else + continue + end + end + if ~isempty(VarType.errorflag) + FF=DataIn.(DataIn.ListVarName{VarType.errorflag}); + X=X(FF==0); + Y=Y(FF==0); + U=U(FF==0); + V=V(FF==0); + end + if nbtps==1 + ListNewVar={'SubRange','NbSites','Coord_tps','U_tps','V_tps'}; + ListNewDim={'nb_tps','nb_subdomain'}; + DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},... + {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}]; + else + ListNewVar={['SubRange_' num2str(nbtps-1)],['NbSites_' num2str(nbtps-1)],['Coord_tps_' num2str(nbtps-1)],['U_tps_' num2str(nbtps-1)] ,['V_tps_' num2str(nbtps-1)]}; + ListNewDim={['nb_tps_' num2str(nbtps-1)],['nb_subdomain_' num2str(nbtps-1)]}; + DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds',ListNewDim{2}},ListNewDim(2),... + {ListNewDim{1},'nb_coord',ListNewDim{2}},ListNewDim,ListNewDim}]; + end + DataOut.ListVarName=[DataIn.ListVarName ListNewVar]; + + [DataOut.(ListNewVar{1}),DataOut.(ListNewVar{2}),DataOut.(ListNewVar{3}),DataOut.(ListNewVar{4}),DataOut.(ListNewVar{5})] =... + filter_tps([X Y],U,V,[],SubDomain,0); + nbvar=numel(DataIn.ListVarName); + + DataOut.VarAttribute{nbvar+3}.Role='coord_tps'; + DataOut.VarAttribute{nbvar+4}=DataIn.VarAttribute{VarType.vector_x};%reproduce attributes of velocity + DataOut.VarAttribute{nbvar+4}.Role='vector_x_tps'; + DataIn.VarAttribute{VarType.vector_x}.VarIndex_tps=nbvar+4;% indicte the correspondance with initial data + DataOut.VarAttribute{nbvar+5}=DataIn.VarAttribute{VarType.vector_y};%reproduce attributes of velocity + DataOut.VarAttribute{nbvar+5}.Role='vector_y_tps'; +% if isfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest') +% DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest'); +% end +% if isfield(DataOut.VarAttribute{VarType.vector_x},'Operation') +% DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'Operation'); +% end + if isfield(DataOut,'ListDimName')%cleaning'FieldRequest' + DataOut=rmfield(DataOut,'ListDimName'); + end + if isfield(DataOut,'DimValue')%cleaning + DataOut=rmfield(DataOut,'DimValue'); + end + end end -if isfield(DataOut,'DimValue')%cleaning - DataOut=rmfield(DataOut,'DimValue'); -end Index: /trunk/src/check_files.m =================================================================== --- /trunk/src/check_files.m (revision 514) +++ /trunk/src/check_files.m (revision 515) @@ -32,5 +32,6 @@ svn_info.status=[]; list_fct={... - 'calc_field';...% defines fields (velocity, vort, div...) from civx data and calculate them + 'calc_field_interp';...% defines fields (velocity, vort, div...) from civx data and calculate them + 'calc_field_tps';...% defines fields (velocity, vort, div...) and calculate them 'cell2tab';... %transform a Matlab cell in a character array suitable for display in a table 'check_field_structure';...% check the validity of the field struture representation consistant with the netcdf format @@ -64,5 +65,4 @@ 'get_file_series';...% determine the list of file names and file indices for functions called by 'series'. 'get_file_type';...% determine info about a file (image, multimage, civdata,...) . - 'griddata_uvmat';...%make 2D linear interpolation using griddata, with input appropriate for both Matlab 6.5 and 7 'hist_update';...% update of a current global histogram by inclusion of a new field 'imadoc2struct';...%convert the image documentation file ImaDoc into a Matlab structure Index: /trunk/src/fill_GUI.m =================================================================== --- /trunk/src/fill_GUI.m (revision 514) +++ /trunk/src/fill_GUI.m (revision 515) @@ -31,6 +31,6 @@ hh=[]; input_data=Param.(fields{ifield}); - display(fields{ifield}) - display(input_data) +% display(fields{ifield}) +% display(input_data) check_done=0; if isfield(handles,fields{ifield}) @@ -68,6 +68,9 @@ end case 'edit' + input_string=''; if isnumeric(input_data) + if numel(input_data)>0 input_string=num2str(input_data(ibox)); + end else input_string=input_data; Index: /trunk/src/find_field_cells.m =================================================================== --- /trunk/src/find_field_cells.m (revision 514) +++ /trunk/src/find_field_cells.m (revision 515) @@ -19,4 +19,8 @@ % .scalar: scalar field (default) % .coord: vector of indices of coordinate variables corresponding to matrix dimensions +% +% .FieldRequest= 'interp', 'derivatives' indicate whether interpolation or derivatives (tps) is needed to calculate the requested field +% .FieldNames = cell of fields to calculate from the fied cell +% % errormsg: error message % @@ -49,57 +53,8 @@ % GNU General Public License (file UVMAT/COPYING.txt) for more details.for input in proj_field and plot_field % group the variables into 'fields' with common dimensions -%------------------------------------------------------------------------ -% function [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data) -% -% OUTPUT: -% CellVaxIndex: cell whose elements are arrays of indices in the list data.ListVarName -% CellvarIndex{i} represents a set of variables with the same dimensions -% NbDim: array with the length of CellVarIndex, giving its space dimension -% CellVarType: cell array of structures with fields -% .coord_x, y, z: indices (in .ListVarname) of variables representing unstructured coordinates x, y, z -% .vector_x,_y,_z: indices of variables giving the vector components x, y, z -% .warnflag: index of warnflag -% .errorflag: index of error flag -% .ancillary: indices of ancillary variables -% .image : B/W image, (behaves like scalar) -% .color : color image, the last index, which is not a coordinate variable, represent the 3 color components rgb -% .discrete: like scalar, but set of data points without continuity, represented as dots in a usual plot, instead of continuous lines otherwise -% .scalar: scalar field (default) -% .coord: vector of indices of coordinate variables corresponding to matrix dimensions -% errormsg: error message -% -% INPUT: -% Data: structure representing fields, output of check_field_structure -% .ListVarName: cell array listing the names (cahr strings) of the variables -% .VarDimName: cell array of cells containing the set of dimension names for each variable of .ListVarName -% .VarAttribute: cell array of structures containing the variable attributes: -% .VarAttribute{ilist}.key=value, where ilist is the variable -% index, key is the name of the attribute, value its value (char string or number) -% -% HELP: -% to get the dimensions of arrays common to the field #icell -% VarIndex=CellVarIndex{icell}; % list of variable indices -% DimIndex=Data.VarDimIndex{VarIndex(1)} % list of dimensions for each variable in the cell #icell -% -%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -% Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org. -%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA -% This file is part of the toolbox UVMAT. -% -% UVMAT is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% UVMAT is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License (file UVMAT/COPYING.txt) for more details. -%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA function [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data) CellVarIndex={}; -NbDim=[]; CellVarType=[]; errormsg=[]; @@ -108,5 +63,4 @@ NbDim=[]; -ivardim=0; VarDimIndex=[]; VarDimName={}; @@ -116,6 +70,9 @@ end -%% role of variables +%% role of variables and list of requested operations Role=num2cell(blanks(nbvar));%initialize a cell array of nbvar blanks +FieldRequest=regexprep(Role,' ',''); % fieldRequest set to '' by default +Operation=cell(size(Role)); % fieldRequest set to {} by default +CheckSub=zeros(size(Role));% =1 for fields to substract Role=regexprep(Role,' ','scalar'); % Role set to 'scalar' by default if isfield(Data,'VarAttribute') @@ -124,8 +81,18 @@ Role{ivar}=Data.VarAttribute{ivar}.Role; end + if isfield(Data.VarAttribute{ivar},'FieldRequest') + FieldRequest{ivar}=Data.VarAttribute{ivar}.FieldRequest; + end + if isfield(Data.VarAttribute{ivar},'Operation') + Operation{ivar}=Data.VarAttribute{ivar}.Operation; + end + if isfield(Data.VarAttribute{ivar},'CheckSub') + CheckSub(ivar)=Data.VarAttribute{ivar}.CheckSub; + end end end %% loop on the list of variables, group them by common dimensions +CellVarType=cell(1,length(CellVarIndex)); for ivar=1:nbvar if ischar(Data.VarDimName{ivar}) @@ -146,15 +113,16 @@ icell=icell+1; CellVarIndex{icell}=ivar;%put the current variable index in the new cell - NbDim(icell)=numel(DimCell);%default - end - - %look for dimension variables -% if numel(DimCell)==1% if the variable has a single dimension -% if strcmp(DimCell{1},Data.ListVarName{ivar}) %|| strcmp(Role{ivar},'dimvar') -% ivardim=ivardim+1; -% VarDimIndex(ivardim)=ivar;%index of the variable -% VarDimName{ivardim}=DimCell{1};%name of the dimension -% end -% end + NbDim(icell)=numel(DimCell);%default + CellVarType{icell}=[]; + end + if ~isempty(FieldRequest{ivar}) + CellVarType{icell}.FieldRequest=FieldRequest{ivar}; + end + if ~isempty(Operation{ivar}) + CellVarType{icell}.Operation=Operation{ivar}; + end + if CheckSub(ivar) + CellVarType{icell}.CheckSub=1; + end end @@ -163,7 +131,9 @@ ind_dim_var_cell=find(checksinglecell); %CoordType(ind_dim_var_cell)='dim_var';% to be used in output +%VarDimIndex=cell(size(ind_dim_var_cell)); +VarDimName=cell(size(ind_dim_var_cell)); for icoord=1:numel(ind_dim_var_cell) -VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)}; -VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1}; + VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)}; + VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1}; end @@ -201,5 +171,4 @@ index_remove=[]; -CellVarType=cell(1,length(CellVarIndex)); for icell=1:length(CellVarIndex) if checksinglecell(icell) @@ -262,10 +231,4 @@ coord(idim)=VarDimIndex(ind_coord); end -% for ivardim=1:numel(VarDimName) -% if strcmp(VarDimName{ivardim},DimCell{idim}) -% coord(idim)=VarDimIndex(ivardim); -% break -% end -% end end NbDim(icell)=numel(find(coord)); Index: unk/src/griddata_uvmat.m =================================================================== --- /trunk/src/griddata_uvmat.m (revision 514) +++ (revision ) @@ -1,28 +1,0 @@ -%'griddata_uvmat': function griddata_uvmat(vec2_X,vec2_Y,vec2_U,vec_X,vec_Y,'linear') -%adapt the input of the matlab function griddata to the appropriate version of Matlab -function ZI = griddata_uvmat(X,Y,Z,XI,YI) -% if ~exist('rho','var')|| isequal(rho,0) -ZI=griddata(X,Y,Z,XI,YI,'linear'); -% txt=ver('MATLAB'); -% Release=txt.Release; -% relnumb=str2num(Release(3:4)); -% if relnumb >= 20 -% ZI=griddata(double(X),double(Y),double(Z),double(XI),double(YI),'linear',{'QJ'}); -% elseif relnumb >=14 -% ZI=griddata(X,Y,Z,XI,YI,'linear',{'QJ'}); -% else -% ZI=griddata(X,Y,Z,XI,YI,'linear'); -% end -% else %smooth with thin plate spline -% [ZI,Z_diff]=patch_uvmat(X,Y,Z,XI,YI,rho); -% diff_norm=mean(Z_diff.*Z_diff) -% ind_good=find(abs(Z_diff)<5*diff_norm); -% nb_remove=numel(Z_diff)-numel(ind_good) -% if nb_remove>0 -% X=X(ind_good); -% Y=Y(ind_good); -% Z=Z(ind_good); -% [ZI,Z_diff]=patch_uvmat(X,Y,Z,XI,YI,rho); -% diff_norm_new=mean(Z_diff.*Z_diff) -% end -% end Index: /trunk/src/plot_field.m =================================================================== --- /trunk/src/plot_field.m (revision 514) +++ /trunk/src/plot_field.m (revision 515) @@ -122,5 +122,5 @@ return end -index_2D=find(NbDim==2,2);%find 2D fields (at most 2) +index_2D=find(NbDim==2);%find 2D fields index_3D=find(NbDim>2,1); if ~isempty(index_3D) @@ -560,5 +560,5 @@ if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected if test_vec - errormsg='error in plot_field: attempt to plot two vector fields'; + errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with mode interp'; return else Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 514) +++ /trunk/src/proj_field.m (revision 515) @@ -84,49 +84,19 @@ ProjData=[]; -%% case of no projection (object is used only as graph display) -if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside')) +%% in the absence of object Type or projection mode, or object coordinaes, the output is empty +if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode') return end - -%% check coincidence of coordinate units -if isfield(FieldData,'CoordUnit') && isfield(ObjectData,'CoordUnit')&&~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit) - errormsg='inconsistent coord units for field and projection object'; +% case of no projection (object is used only as graph display) +if isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside') return end - -%% in the absence of object Type or projection mode, or object coordinaes, the input field is just tranfered without change -if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode') - ProjData=FieldData; - return -end -if ~isfield(ObjectData,'Coord') +if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord) if strcmp(ObjectData.Type,'plane') ObjectData.Coord=[0 0 0];%default else - ProjData=FieldData; return end end - -%% OBSOLETE -if isfield(ObjectData,'XMax') && ~isempty(ObjectData.XMax) - ObjectData.RangeX(1)=ObjectData.XMax; -end -if isfield(ObjectData,'XMin') && ~isempty(ObjectData.XMin) - ObjectData.RangeX(2)=ObjectData.XMin; -end -if isfield(ObjectData,'YMax') && ~isempty(ObjectData.YMax) - ObjectData.RangeY(1)=ObjectData.YMax; -end -if isfield(ObjectData,'YMin') && ~isempty(ObjectData.YMin) - ObjectData.RangeY(2)=ObjectData.YMin; -end -if isfield(ObjectData,'ZMax') && ~isempty(ObjectData.ZMax) - ObjectData.RangeZ(1)=ObjectData.ZMax; -end -if isfield(ObjectData,'ZMin') && ~isempty(ObjectData.ZMin) - ObjectData.RangeZ(2)=ObjectData.ZMin; -end -%%%%%%%%%% %% apply projection depending on the object type @@ -401,5 +371,5 @@ testproj(VarType.color)=1; VarIndex=VarIndex(find(testproj(VarIndex)));%select only the projected variables - if testX %case of unstructured coordinates + if check_unstructured%case of unstructured coordinates eval(['nbpoint=numel(FieldData.' FieldData.ListVarName{VarIndex(1)} ');']) for ivar=[VarIndex VarType.coord_x VarType.coord_y VarType.errorflag] @@ -547,6 +517,5 @@ ProjData.NbDim=1; %initialisation of the input parameters and defaultoutput -ProjMode='projection';%direct projection on the line by default -if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; +ProjMode=ObjectData.ProjMode; % ProjAngle=90; %90 degrees projection by default @@ -609,5 +578,5 @@ continue end - testX=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates + check_unstructured=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates test_tps=~isempty(VarType.coord_tps); testU=~isempty(VarType.vector_x) && ~isempty(VarType.vector_y);% test for vectors @@ -631,5 +600,5 @@ end % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y - if testX + if check_unstructured if ~isequal(ProjMode,'interp') if width==0 @@ -664,5 +633,5 @@ %%%%%%% % A FAIRE CALCULER MEAN DES QUANTITES %%%%%% %case of unstructured coordinates - if testX + if check_unstructured for ip=1:siz_line(1)-1 %Loop on the segments of the polyline linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip)); @@ -917,15 +886,4 @@ %----------------------------------------------------------------- -%% initialisation of the input parameters of the projection plane -ProjMode='projection';%direct projection by default -if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; - -%% axis origin -if isempty(ObjectData.Coord) - ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default - ObjectData.Coord(1,2)=0; - ObjectData.Coord(1,3)=0; -end - %% rotation angles PlaneAngle=[0 0 0]; @@ -951,6 +909,6 @@ %% mesh sizes DX and DY -DX=0; -DY=0; %default +DX=FieldData.Mesh; +DY=FieldData.Mesh; %default if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) DX=abs(ObjectData.DX);%mesh of interpolation points @@ -959,5 +917,5 @@ DY=abs(ObjectData.DY);%mesh of interpolation points end -if ~strcmp(ProjMode,'projection') && (DX==0||DY==0) +if ~strcmp(ObjectData.ProjMode,'projection') && (DX==0||DY==0) errormsg='DX or DY missing'; display(errormsg) @@ -970,4 +928,8 @@ testYMin=0; testYMax=0; +XMin=FieldData.XMin;%default +XMax=FieldData.XMax;%default +YMin=FieldData.YMin;%default +YMax=FieldData.YMax;%default if isfield(ObjectData,'RangeX') XMin=min(ObjectData.RangeX); @@ -1003,9 +965,6 @@ ListIndex={}; -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions -%----------------------------------------------------------------- -idimvar=0; - +% CellVarIndex=cells of variable index arrays [CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(FieldData); if ~isempty(errormsg) @@ -1014,4 +973,32 @@ end +%% projection modes +check_grid=0; +ProjMode=cell(size(VarTypeCell)); +for icell=1:numel(VarTypeCell)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates + ProjMode{icell}=ObjectData.ProjMode; + if isfield(VarTypeCell{icell},'FieldRequest') + switch VarTypeCell{icell}.FieldRequest + case 'interp' + ProjMode{icell}='interp'; + case 'derivatives' + ProjMode{icell}='filter'; + end + end + if strcmp(ProjMode{icell},'interp')||strcmp(ProjMode{icell},'filter') + check_grid=1; + end +end + +%% define the new coordinates in case of interpolation on a grid +if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates + ProjData.ListVarName={'coord_y','coord_x'}; + ProjData.VarDimName={'coord_y','coord_x'}; + ProjData.VarAttribute={[],[]}; + ProjData.coord_y=[YMin YMax]; + ProjData.coord_x=[XMin XMax]; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS % CellVarIndex=cells of variable index arrays @@ -1020,4 +1007,6 @@ nbcoord=0;%number of added coordinate variables brought by projection nbvar=0; +vector_x_proj=[]; +vector_y_proj=[]; for icell=1:length(CellVarIndex) NbDim=NbDimVec(icell); @@ -1027,7 +1016,4 @@ VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName VarType=VarTypeCell{icell}; -% ivar_X=VarType.coord_x; -% ivar_Y=VarType.coord_y; -% ivar_Z=VarType.coord_z; ivar_U=VarType.vector_x; ivar_V=VarType.vector_y; @@ -1037,9 +1023,4 @@ end ivar_W=VarType.vector_z; -% ivar_Anc=VarType.ancillary; -% test_anc=zeros(size(VarIndex)); -% test_anc(ivar_Anc)=ones(size(ivar_Anc)); -% ivar_F=VarType.warnflag; -% ivar_FF=VarType.errorflag; %type of coordinates @@ -1058,5 +1039,5 @@ end coord_z=0;%default - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch CoordType @@ -1064,14 +1045,14 @@ %% case of input fields with unstructured coordinates case 'unstructured' - if strcmp(ObjectData.ProjMode,'filter') - continue %skip for filter (needs tps) - end - coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x}); - coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y}); + if strcmp(ProjMode{icell},'filter') + continue %skip for filter (needs tps field cell) + end + coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});% initial x coordinates + coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});% initial y coordinates if ~isempty(VarType.coord_z) - coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z}); + coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});% initial z coordinates end - % translate initial coordinates + % translate initial coordinates to account for the new origin coord_x=coord_x-ObjectData.Coord(1,1); coord_y=coord_y-ObjectData.Coord(1,2); @@ -1085,5 +1066,4 @@ fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane indcut=find(abs(fieldZ) <= width); - size(indcut) for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; @@ -1096,5 +1076,5 @@ end - %rotate coordinates if needed: + %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane Psi=PlaneAngle(1); Theta=PlaneAngle(2); @@ -1104,6 +1084,5 @@ coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); coord_Y=coord_Y+coord_z *sin(Theta); - coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER - + coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); else @@ -1112,5 +1091,5 @@ end - %restriction to the range of x and y if imposed + %restriction to the range of X and Y if imposed by the projection object testin=ones(size(coord_X)); %default testbound=0; @@ -1139,5 +1118,5 @@ for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; - eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) + FieldData.(VarName)=FieldData.(VarName)(indcut); end coord_X=coord_X(indcut); @@ -1149,6 +1128,7 @@ % different cases of projection - switch ObjectData.ProjMode - case 'projection' + switch ProjMode{icell} + case 'projection' + nbvar=numel(ProjData.ListVarName); for ivar=VarIndex %transfer variables to the projection plane VarName=FieldData.ListVarName{ivar}; @@ -1172,20 +1152,6 @@ coord_x_proj=XMin:DX:XMax; coord_y_proj=YMin:DY:YMax; - DimCell={'coord_y','coord_x'}; - ProjData.ListVarName={'coord_y','coord_x'}; - ProjData.VarDimName={'coord_y','coord_x'}; - nbcoord=2; - ProjData.coord_y=[YMin YMax]; - ProjData.coord_x=[XMin XMax]; -% if isempty(ivar_X), ivar_X=0; end; -% if isempty(ivar_Y), ivar_Y=0; end; -% if isempty(ivar_Z), ivar_Z=0; end; -% if isempty(ivar_U), ivar_U=0; end; -% if isempty(ivar_V), ivar_V=0; end; -% if isempty(ivar_W), ivar_W=0; end; -% if isempty(ivar_F), ivar_F=0; end; -% if isempty(ivar_FF), ivar_FF=0; end; - % if ~isequal(ivar_FF,0) - if ~isempty(VarType.errorflag) + [XI,YI]=meshgrid(coord_x_proj,coord_y_proj); + if ~isempty(VarType.errorflag) VarName_FF=FieldData.ListVarName{VarType.errorflag}; indsel=find(FieldData.(VarName_FF)==0); @@ -1193,57 +1159,55 @@ coord_Y=coord_Y(indsel); end - - FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); testFF=0; - %for ivar=VarIndex % loop on field variables to project - for ivar=[VarType.scalar VarType.vector_x VarType.vector_y] - VarName=FieldData.ListVarName{ivar}; -% if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||... -% ( numel(test_anc)>=ivar && test_anc(ivar)==1)) - ivar_new=ivar_new+1; - ProjData.ListVarName=[ProjData.ListVarName {VarName}]; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; - if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar - ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; - end - %if ~isequal(ivar_FF,0) - if ~isempty(VarType.errorflag) - FieldData.(VarName)=FieldData.(VarName)(indsel); - end - ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj'); - varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj)); - FFlag= isnan(varline); %detect undefined values NaN - indnan=find(FFlag); - if~isempty(indnan) - varline(indnan)=zeros(size(indnan)); - ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj)); - FF(indnan)=ones(size(indnan)); - testFF=1; - end - if ivar==ivar_U - ivar_U=ivar_new; - end - if ivar==ivar_V - ivar_V=ivar_new; - end - if ivar==ivar_W - ivar_W=ivar_new; - end -% end - end - if testFF - ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); - ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; - ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; - end -% case 'filter'%interpolate data on a regular grid -% errormsg='tps required for filter option' - - end - - %% case of tps interpolation (applies only in filter mode) + nbvar=numel(ProjData.ListVarName); + if isfield(VarType,'vector_x')&&isfield(VarType,'vector_y')&&~isempty(VarType.vector_x) + VarName_x=FieldData.ListVarName{VarType.vector_x}; + VarName_y=FieldData.ListVarName{VarType.vector_y}; + if ~isempty(VarType.errorflag) + FieldData.(VarName_x)=FieldData.(VarName_x)(indsel); + FieldData.(VarName_y)=FieldData.(VarName_y)(indsel); + end + FieldVar=cat(2,FieldData.(VarName_x),FieldData.(VarName_y)); + if ~isfield(VarType,'CheckSub') || ~VarType.CheckSub + vector_x_proj=numel(ProjData.ListVarName)+1; + vector_y_proj=numel(ProjData.ListVarName)+2; + end + end + if ~isempty(VarType.scalar) + VarName_scalar=FieldData.ListVarName{VarType.scalar}; + if ~isempty(VarType.errorflag) + FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel); + end + FieldVar=FieldData.(VarName_scalar); + end + [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldVar,VarType.Operation,XI,YI); + if isfield(VarType,'CheckSub') && VarType.CheckSub && ~isempty(vector_x_proj) + ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1}; + ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2}; + else + VarDimName=cell(size(ListFieldProj)); + for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise + ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+',''); + if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName))) + ListFieldProj{ilist}=[ListFieldProj{ilist} '_1']; + end + ProjData.(ListFieldProj{ilist})=VarVal{ilist}; + VarDimName{ilist}={'coord_y','coord_x'}; + end + ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; + ProjData.VarDimName=[ProjData.VarDimName VarDimName]; + ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; + end + end + + %% case of tps interpolation (applies only in filter mode and for spatial derivatives) case 'tps' - if strcmp(ObjectData.ProjMode,'filter') + if strcmp(ProjMode{icell},'filter') + Coord=FieldData.(FieldData.ListVarName{VarType.coord_tps}); + NbSites=FieldData.(FieldData.ListVarName{VarType.nbsites_tps}); + SubRange=FieldData.(FieldData.ListVarName{VarType.subrange_tps}); + if isfield(VarType,'vector_x_tps')&&isfield(VarType,'vector_y_tps') + FieldVar=cat(3,FieldData.(FieldData.ListVarName{VarType.vector_x_tps}),FieldData.(FieldData.ListVarName{VarType.vector_y_tps})); + end coord_x_proj=XMin:DX:XMax; coord_y_proj=YMin:DY:YMax; @@ -1253,17 +1217,22 @@ XI=XI+ObjectData.Coord(1,1); YI=YI+ObjectData.Coord(1,2); - DataOut=calc_field(FieldData.FieldList,FieldData,cat(3,XI,YI)); - ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName]; - ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName]; - ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute]; - for ilist=1:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise - VarName=DataOut.ListVarName{ilist}; -% ProjData.(VarName)=DataOut.(VarName); - if ilist>=3 - ProjData.(VarName)=reshape(DataOut.(VarName),np_y,np_x); - end - end - ProjData.coord_x=[XMin XMax]; - ProjData.coord_y=[YMin YMax]; + [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbSites,SubRange,FieldVar,VarType.Operation,cat(3,XI,YI)); + ListFieldProj=(fieldnames(DataOut))'; + VarDimName=cell(size(ListFieldProj)); + for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise + VarName=ListFieldProj{ilist}; + ProjData.(VarName)=DataOut.(VarName); + VarDimName{ilist}={'coord_y','coord_x'}; + end +% if ~isfield(ProjData,'coord_x') +% ProjData.coord_x=[XMin XMax]; +% ProjData.coord_y=[YMin YMax]; +% ListFieldProj=[{'coord_x','coord_y'} ListFieldProj']; +% VarDimName=[{'coord_x','coord_y'} VarDimName]; +% VarAttribute=[{[],[]} VarAttribute]; +% end + ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; + ProjData.VarDimName=[ProjData.VarDimName VarDimName]; + ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; end @@ -1272,5 +1241,5 @@ VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions - eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions + DimValue=size(FieldData.(VarName));%input matrix dimensions DimValue(DimValue==1)=[];%remove singleton dimensions NbDim=numel(DimValue);%update number of space dimensions @@ -1405,7 +1374,15 @@ end % case with no interpolation - if isequal(ProjMode,'projection') && (~testangle || test90y || test90x) + if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x) if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax - ProjData=FieldData;% no change by projection + ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)]; + ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)]; + ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)]; + ProjData.(AYProjName)=FieldData.(AYName); + ProjData.(AXProjName)=FieldData.(AXName); + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + ProjData.(VarName)=FieldData.(VarName);% no change by projection + end else indY=NbDim-1; @@ -1489,5 +1466,5 @@ YIMA=reshape(round(YIMA),1,npX*npY); flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image - if isequal(ObjectData.ProjMode,'filter') + if isequal(ProjMode{icell},'filter') npx_filter=ceil(abs(DX/DAX)); npy_filter=ceil(abs(DY/DAY)); @@ -1581,4 +1558,27 @@ end +% %prepare substraction in case of two input fields +% SubData.ListVarName={}; +% SubData.VarDimName={}; +% SubData.VarAttribute={}; +% check_remove=zeros(size(ProjData.ListVarName)); +% for iproj=1:numel(ProjData.VarAttribute) +% if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1) +% VarName=ProjData.ListVarName{iproj}; +% SubData.ListVarName=[SubData.ListVarName {VarName}]; +% SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}]; +% SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}]; +% SubData.(VarName)=ProjData.(VarName); +% check_remove(iproj)=1; +% end +% end +% if ~isempty(find(check_remove)) +% ind_remove=find(check_remove); +% ProjData.ListVarName(ind_remove)=[]; +% ProjData.VarDimName(ind_remove)=[]; +% ProjData.VarAttribute(ind_remove)=[]; +% ProjData=sub_field(ProjData,[],SubData); +% end + %----------------------------------------------------------------- %projection in a volume @@ -1586,8 +1586,4 @@ %----------------------------------------------------------------- ProjData=FieldData;%default output - -%% initialisation of the input parameters of the projection plane -ProjMode='projection';%direct projection by default -if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; %% axis origin Index: /trunk/src/read_civdata.m =================================================================== --- /trunk/src/read_civdata.m (revision 514) +++ /trunk/src/read_civdata.m (revision 515) @@ -10,5 +10,5 @@ % .NbCoord: number of vector components % .NbDim: number of dimensions (=2 or 3) -% .dt: time interval for the corresponding image pair +% .Dt: time interval for the corresponding image pair % .Time: absolute time (average of the initial image pair) % .CivStage: =0, @@ -63,7 +63,19 @@ end errormsg=''; +if ischar(FieldNames), FieldNames={FieldNames}; end; +FieldRequest='vec'; +for ilist=1:length(FieldNames) + if ~isempty(FieldNames{ilist}) + switch FieldNames{ilist} + case{'U','V','norm(U,V)'} + FieldRequest='interp'; + case {'curl(U,V)','div(U,V)','strain(U,V)'} + FieldRequest='derivatives'; + end + end +end %% reading data -[varlist,role,units,VelTypeOut]=varcivx_generator(FieldNames,VelType,CivStage); +[varlist,role,VelTypeOut]=varcivx_generator(FieldRequest,VelType,CivStage); if isempty(varlist) erromsg=['error in read_civdata: unknow velocity type ' VelType]; @@ -97,10 +109,31 @@ end Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'Dt','Time'}]; +ivar_U_tps=[]; var_ind=find(vardetect); for ivar=1:numel(var_ind) Field.VarAttribute{ivar}.Role=role{var_ind(ivar)}; - Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)}; + Field.VarAttribute{ivar}.FieldRequest=FieldRequest; + if strcmp(role{var_ind(ivar)},'vector_x') + Field.VarAttribute{ivar}.Operation=FieldNames; + ivar_U=ivar; + end + if strcmp(role{var_ind(ivar)},'vector_x_tps') + Field.VarAttribute{ivar}.Operation=FieldNames; + ivar_U_tps=ivar; + end +% Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)}; Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel end +if ~isempty(ivar_U_tps) + Field.VarAttribute{ivar_U}.VarIndex_tps=ivar_U_tps; +end +% if strcmp(FieldRequest,'derivatives')% fields will be calculated from the tps +% Field.VarAttribute{ivar_U_tps}.FieldRequest=FieldRequest; +% Field.VarAttribute{ivar_U_tps}.Operation=FieldNames; +% else% fields will be calculated from the initial fields +% Field.VarAttribute{ivar_U}.FieldRequest=FieldRequest;% +% Field.VarAttribute{ivar_U}.Operation=FieldNames; +% end + Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'NbCoord','NbDim','TimeUnit','CoordUnit'}]; % %% update list of global attributes @@ -121,29 +154,15 @@ % if vel_type=[] or'*', a priority choice is done, civ2 considered better than civ1 ) -function [var,role,units,vel_type_out,errormsg]=varcivx_generator(FieldNames,vel_type,CivStage) +function [var,role,vel_type_out,errormsg]=varcivx_generator(FieldRequest,vel_type,CivStage) %% default input values if ~exist('vel_type','var'),vel_type='';end; if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed -if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar -if ischar(FieldNames), FieldNames={FieldNames}; end; errormsg=''; %% select the priority order for automatic vel_type selection -testder=0; -testpatch=0; -for ilist=1:length(FieldNames) - if ~isempty(FieldNames{ilist}) - switch FieldNames{ilist} - case{'u','v'} - testpatch=1; - case {'vort','div','strain'} - testder=1; - end - end -end -if strcmp(vel_type,'civ2') && testder +if strcmp(vel_type,'civ2') && strcmp(FieldRequest,'derivatives') vel_type='filter2'; -elseif strcmp(vel_type,'civ1') && testder +elseif strcmp(vel_type,'civ1') && strcmp(FieldRequest,'derivatives') vel_type='filter1'; end @@ -153,5 +172,5 @@ vel_type='filter2'; case {4,5}% civ2 available but not filter2 - if testder% derivatives needed + if strcmp(FieldRequest,'derivatives')% derivatives needed vel_type='filter1'; else @@ -171,5 +190,5 @@ 'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'}; role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'}; - units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''}; + % units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''}; case 'filter1' var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';... @@ -178,10 +197,10 @@ role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',... 'vector_y_tps','vector_z_tps','ancillary','ancillary'}; - units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''}; + % units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''}; case 'civ2' var={'X','Y','Z','U','V','W','C','F','FF';... 'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF'}; role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'}; - units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''}; + % units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''}; case 'filter2' var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';... @@ -190,5 +209,8 @@ role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',... 'vector_y_tps','vector_z_tps','ancillary','ancillary'}; - units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''}; + % units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''}; +end +if ~strcmp(FieldRequest,'derivatives') + var=var(:,1:9);%suppress tps if not needed end vel_type_out=vel_type; Index: /trunk/src/sub_field.m =================================================================== --- /trunk/src/sub_field.m (revision 514) +++ /trunk/src/sub_field.m (revision 515) @@ -7,5 +7,5 @@ % when scalar and vectors are combined, the fields are just merged in a single matlab structure for common visualisation %----------------------------------------------------------------------- -% function SubData=sub_field(Field,Field_1) +% function SubData=sub_field(Field,XmlData,Field_1) % % OUPUT: @@ -16,8 +16,17 @@ % Field_1:matlab structure representing the second field -function [SubData,errormsg]=sub_field(Field,Field_1) +function SubData=sub_field(Field,XmlData,Field_1) +SubData=[]; +if strcmp(Field,'*') + return +end +if nargin<3 + SubData=Field; + return +end %% global attributes -test_attr=0; +SubData.ListGlobalAttribute={};%default +%transfer global attributes of Field if isfield(Field,'ListGlobalAttribute') SubData.ListGlobalAttribute=Field.ListGlobalAttribute; @@ -26,26 +35,16 @@ SubData.(AttrName)=Field.(AttrName); end - test_attr=1; end +%transfer global attributes of Field_1 if isfield(Field_1,'ListGlobalAttribute') for ilist=1:numel(Field_1.ListGlobalAttribute) - test_1=1; AttrName=Field_1.ListGlobalAttribute{ilist}; - if test_attr - for i_prev=1:numel(Field.ListGlobalAttribute) - if isequal(Field.ListGlobalAttribute{i_prev},AttrName) - test_1=0; %attribute already written - eval(['Val=Field.' AttrName ';']) - eval(['Val_1=Field_1.' AttrName ';']) - if isequal(Val,Val_1) - break% data already written - else - eval(['SubData.' AttrName '_1=Field_1.' AttrName ';']) - end - end - end + AttrNameNew=AttrName; + while ~isempty(find(strcmp(AttrNameNew,SubData.ListGlobalAttribute)))&&~isequal(Field_1.(AttrNameNew),Field.(AttrNameNew)) + AttrNameNew=[AttrNameNew '_1']; end - if test_1 - eval(['SubData.' AttrName '=Field_1.' AttrName ';']) + if ~isfield(Field,AttrName) || ~isequal(Field_1.(AttrName),Field.(AttrName)) + SubData.ListGlobalAttribute=[SubData.ListGlobalAttribute {AttrNameNew}]; + SubData.(AttrNameNew)=Field_1.(AttrName); end end @@ -53,4 +52,5 @@ %% variables +%reproduce variables of the first field and list its dimensions SubData.ListVarName=Field.ListVarName; SubData.VarDimName=Field.VarDimName; @@ -58,339 +58,124 @@ SubData.VarAttribute=Field.VarAttribute; end -%reproduce the first field Field by default -for ivar=1:numel(Field.ListVarName) - VarName=Field.ListVarName{ivar}; - SubData.(VarName)=Field.(VarName); +ListDimName={}; +for ilist=1:numel(Field.ListVarName) + VarName=Field.ListVarName{ilist}; + SubData.(VarName)=Field.(VarName); + SubData.VarAttribute{ilist}.CheckSub=0; + DimCell=Field.VarDimName{ilist}; + if ischar(DimCell) + DimCell={DimCell}; + end + for idim=1:numel(DimCell) + if isempty(find(strcmp(DimCell{idim},ListDimName))) + ListDimName=[ListDimName DimCell(idim)]; + end + end end -%% check the two input fields -[CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_cells(Field); -if ~isempty(errormsg) - errormsg=['invalid first input to sub_field:' errormsg]; - return -end -[CellVarIndex_1,NbDim_1,VarTypeCell_1,errormsg]=find_field_cells(Field_1); -if ~isempty(errormsg) - errormsg=['invalid second input to sub_field:' errormsg]; - return -end -iselect=find(NbDim==2); -for icell=iselect - if ~isempty(VarTypeCell{icell}.coord_tps) - NbDim(icell)=0; +%% field request +FieldRequest=cell(size(Field.ListVarName)); +for ilist=1:numel(Field.VarAttribute) + if isfield(Field.VarAttribute{ilist},'FieldRequest') + FieldRequest{ilist}=Field.VarAttribute{ilist}.FieldRequest; end end -iselect=find(NbDim==2); -if ~isequal(numel(iselect),1) - errormsg='invalid first input to sub_field: it must contain a single 2D field cell'; - return -end -iselect_1=find(NbDim_1==2); -for icell=iselect_1 - if ~isempty(VarTypeCell_1{icell}.coord_tps) - NbDim_1(icell)=0; +FieldRequest_1=cell(size(Field_1.ListVarName)); +for ilist=1:numel(Field_1.VarAttribute) + if isfield(Field_1.VarAttribute{ilist},'FieldRequest') + FieldRequest_1{ilist}=Field_1.VarAttribute{ilist}.FieldRequest; end end -iselect_1=find(NbDim_1==2); -if ~isequal(numel(iselect_1),1) - errormsg='invalid second input to sub_field: it must contain a single 2D field cell'; - return -end -VarType=VarTypeCell{iselect}; -VarType_1=VarTypeCell_1{iselect_1}; -testX=~isempty(VarType.coord_x)&& ~isempty(VarType.coord_y);%unstructured coordiantes -testX_1=~isempty(VarType_1.coord_x)&& ~isempty(VarType_1.coord_y);%unstructured coordiantes -testU=~isempty(VarType.vector_x)&& ~isempty(VarType.vector_y);%vector field -testU_1=~isempty(VarType_1.vector_x)&& ~isempty(VarType_1.vector_y);%vector field -testfalse_1=~isempty(VarType_1.errorflag); -ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields -if numel(ivar_C)>1 - errormsg='too many scalar fields in the first input of sub_field.m'; - return -end -ivar_C_1=[VarType_1.scalar VarType_1.image VarType_1.color VarType_1.ancillary]; %defines index (indices) for the scalar or ancillary fields -if numel(ivar_C_1)>1 - errormsg='too many scalar fields in the second input of sub_field.m'; - return + +%% rename the dimensions of the second field if identical to those of the first +for ilist=1:numel(Field_1.VarDimName) + DimCell=Field_1.VarDimName{ilist}; + if ischar(DimCell) + DimCell={DimCell}; + end + for idim=1:numel(DimCell) + ind_dim=find(strcmp(DimCell{idim},ListDimName)); + if ~isempty(ind_dim) + if ischar(Field_1.VarDimName{ilist}) + Field_1.VarDimName{ilist}=Field_1.VarDimName(ilist); + end + Field_1.VarDimName{ilist}{idim}=[ListDimName{ind_dim} '_1']; + end + end end -%% substract two vector fields or two scalars -if (testU && testU_1) || (~testU && ~testU_1) - %check coincidence in positions - %unstructured coordinates for the first field - if testX - XName=Field.ListVarName{VarType.coord_x}; - YName=Field.ListVarName{VarType.coord_y}; - vec_X=Field.(XName); - vec_Y=Field.(YName); - nbpoints=numel(vec_X); - vec_X=reshape(vec_X,nbpoints,1); - vec_Y=reshape(vec_Y,nbpoints,1); - if testX_1 %unstructured coordinates for the second field - X_1_Name=Field_1.ListVarName{VarType_1.coord_x}; - Y_1_Name=Field_1.ListVarName{VarType_1.coord_y}; - eval(['vec_X_1=Field_1.' X_1_Name ';']) - eval(['vec_Y_1=Field_1.' Y_1_Name ';']) +%look for coordinates common to Field in Field_1 +ind_remove=zeros(size(Field_1.ListVarName)); +for ilist=1:numel(Field_1.ListVarName) + if ischar(Field_1.VarDimName{ilist})||numel(Field_1.VarDimName{ilist})==1 + OldDim=Field_1.VarDimName{ilist}; + if ischar(OldDim) + OldDim=Field_1.VarDimName(ilist); + end + VarVal=Field_1.(Field_1.ListVarName{ilist}); + for i1=1:numel(Field.ListVarName) + if (isempty(FieldRequest{i1})&&isempty(FieldRequest_1{ilist})||strcmp(FieldRequest{i1},FieldRequest_1{ilist})) && isequal(Field.(Field.ListVarName{i1}),VarVal) + ind_remove(ilist)=1; + NewDim=Field.VarDimName{i1}; + if ischar(NewDim) + NewDim={NewDim}; + end + Field_1.VarDimName=regexprep_r(Field_1.VarDimName,OldDim{1},NewDim{1}); + end + end + end +end +Field_1.ListVarName(find(ind_remove))=[];%removes these redondent coordinates +Field_1.VarDimName(find(ind_remove))=[]; +Field_1.VarAttribute(find(ind_remove))=[]; - else %structured coordinates for the second field - y_1_Name=Field_1.ListVarName{VarType_1.coord(1)}; - x_1_Name=Field_1.ListVarName{VarType_1.coord(2)}; - eval(['y_1=Field_1.' y_1_Name ';']) - eval(['x_1=Field_1.' x_1_Name ';']) - if isequal(numel(x_1),2) - x_1=linspace(x_1(1),x_1(2),nbpoints_x_1); - end - if isequal(numel(y_1),2) - y_1=linspace(y_1(1),y_1(2),nbpoints_y_1); - end - [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1); - end - vec_X_1=reshape(vec_X_1,[],1); - vec_Y_1=reshape(vec_Y_1,[],1); - if testfalse_1 - FFName_1=Field_1.ListVarName{VarType_1.errorflag}; - eval(['vec_FF_1=Field_1.' FFName_1 ';']) - vec_FF_1=reshape(vec_FF_1,[],1); - indsel=find(~vec_FF_1); - vec_X_1=vec_X_1(indsel); - vec_Y_1=vec_Y_1(indsel); - end - if testU % vector fields - U_1_Name=Field_1.ListVarName{VarType_1.vector_x}; - V_1_Name=Field_1.ListVarName{VarType_1.vector_y}; - eval(['vec_U_1=Field_1.' U_1_Name ';']) - eval(['vec_V_1=Field_1.' V_1_Name ';']) - nbpoints_x_1=size(vec_U_1,2); - nbpoints_y_1=size(vec_U_1,1); - vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1); - vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1); - if testfalse_1 - vec_U_1=vec_U_1(indsel); - vec_V_1=vec_V_1(indsel); - end - else %(~testU && ~testU_1) - A_1_Name=Field_1.ListVarName{ivar_C_1}; - eval(['vec_A_1=Field_1.' A_1_Name ';']) - nbpoints_x_1=size(vec_A_1,2); - nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2) - vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1); - if testfalse_1 - vec_A_1=vec_A_1(indsel); - end - end - - if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same - if testU - vec_U_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_U_1,vec_X,vec_Y); %interpolate vectors in the second field - vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y); %interpolate vectors in the second field - else - vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,double(vec_A_1),vec_X,vec_Y); %interpolate vectors in the second field - end - end - if testU - UName=Field.ListVarName{VarType.vector_x}; - VName=Field.ListVarName{VarType.vector_y}; - eval(['vec_U=Field.' UName ';']) - eval(['vec_V=Field.' VName ';']) - vec_U=reshape(vec_U,numel(vec_U),1); - vec_V=reshape(vec_V,numel(vec_V),1); - eval(['SubData.' UName '=vec_U-vec_U_1;']) - eval(['SubData.' VName '=vec_V-vec_V_1;']) - else - AName=Field.ListVarName{ivar_C}; - % size(Field.vort) - eval(['SubData.' AName '=Field.' AName '-vec_A_1;']) - end - else %structured coordiantes - XName=Field.ListVarName{VarType.coord(2)}; - YName=Field.ListVarName{VarType.coord(1)}; - eval(['x=Field.' XName ';']) - eval(['y=Field.' YName ';']) - if testX_1 %unstructured coordinates for the second field - errormsg='the second input scalar is not on a regular grid: comparison option not implemented'; - return - else - XName_1=Field.ListVarName{VarType_1.coord(2)}; - YName_1=Field.ListVarName{VarType_1.coord(1)}; - eval(['x_1=Field_1.' XName_1 ';']) - eval(['y_1=Field_1.' YName_1 ';']) - end - if testU % vector fields - UName=Field.ListVarName{VarType.vector_x}; - VName=Field.ListVarName{VarType.vector_y}; - U_1_Name=Field_1.ListVarName{VarType_1.vector_x}; - V_1_Name=Field_1.ListVarName{VarType_1.vector_y}; - eval(['U_1=Field_1.' U_1_Name ';']) - eval(['V_1=Field_1.' V_1_Name ';']) - if ~isequal(x_1,x)||~isequal(y_1,y) - [X_1,Y_1]=meshgrid(x_1,y_1); - U_1 =interp2(X_1,Y_1,U_1,x,y'); - V_1 =interp2(X_1,Y_1,V_1,x,y'); - end - eval(['SubData.' UName '=Field.' UName '-U_1;']) - eval(['SubData.' VName '=Field.' VName '-V_1;']) - else - AName=Field.ListVarName{ivar_C}; - A_1_Name=Field_1.ListVarName{ivar_C_1}; - eval(['A_1=double(Field_1.' A_1_Name ');']) - if ~isequal(x_1,x)||~isequal(y_1,y) - [X_1,Y_1]=meshgrid(x_1,y_1); - A_1 =interp2(X_1,Y_1,A_1,x,y'); - end - eval(['SubData.' AName '=double(Field.' AName ')-A_1;']) - end - end +%append the other variables of the second field, modifying their name if needed +for ilist=1:numel(Field_1.ListVarName) + VarName=Field_1.ListVarName{ilist}; + ind_prev=find(strcmp(VarName,Field.ListVarName)); + if isempty(ind_prev)% variable name does not exist in Field + VarNameNew=VarName; + else % variable name exists in Field + VarNameNew=[VarName '_1']; + end + SubData.ListVarName=[SubData.ListVarName {VarNameNew}]; + SubData.VarDimName=[SubData.VarDimName Field_1.VarDimName(ilist)]; + SubData.(VarNameNew)=Field_1.(VarName); + SubData.VarAttribute=[SubData.VarAttribute Field_1.VarAttribute(ilist)]; + SubData.VarAttribute{end}.CheckSub=1;% mark that the field needs to be substracted end -% merge a vector field and a scalar as second input -if testU && ~testU_1 - AName_1=Field_1.ListVarName{ivar_C_1}; - if isfield(Field_1,'VarAttribute') & numel(Field_1.VarAttribute)>=ivar_C_1 - AAttr=Field_1.VarAttribute{ivar_C_1} ; - else - AAttr=[]; +%append the other variables of the second field, modifying their name if needed + +%% substrat fields when possible +[CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(SubData); +ind_remove=zeros(size(SubData.ListVarName)); +ivar=[]; +ivar_1=[]; +for icell=1:numel(CellVarIndex) + if ~isempty(CellVarIndex{icell}) + if isfield(CellVarType{icell},'scalar') && numel(CellVarType{icell}.scalar)==2 && SubData.VarAttribute{CellVarType{icell}.scalar(2)}.CheckSub; + ivar=[ivar CellVarType{icell}.scalar(1)]; + ivar_1=[ivar_1 CellVarType{icell}.scalar(2)]; + end + if isfield(CellVarType{icell},'vector_x') && numel(CellVarType{icell}.vector_x)==2 && SubData.VarAttribute{CellVarType{icell}.vector_x(2)}.CheckSub; + ivar=[ivar CellVarType{icell}.vector_x(1)]; + ivar_1=[ivar_1 CellVarType{icell}.vector_x(2)]; + end + if isfield(CellVarType{icell},'vector_y') && numel(CellVarType{icell}.vector_y)==2 && SubData.VarAttribute{CellVarType{icell}.vector_y(2)}.CheckSub; + ivar=[ivar CellVarType{icell}.vector_y(1)]; + ivar_1=[ivar_1 CellVarType{icell}.vector_y(2)]; + end end - if testX_1 %unstructured coordinate - XName_1=Field_1.ListVarName{VarType_1.coord_x}; - YName_1=Field_1.ListVarName{VarType_1.coord_y}; - DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]); - if isfield(Field_1,'VarAttribute') - if numel(Field_1.VarAttribute)>=VarType_1.coord_x - XAttr=Field_1.VarAttribute{VarType_1.coord_x}; - else - XAttr=[]; - end - if numel(Field_1.VarAttribute)>=VarType_1.coord_y - YAttr=Field_1.VarAttribute{VarType_1.coord_y}; - else - YAttr=[]; - end - SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr}]; - end - else - XName_1=Field_1.ListVarName{VarType_1.coord(2)}; - YName_1=Field_1.ListVarName{VarType_1.coord(1)}; -% DimCell=[{YName_1} {XName_1}]; - if isfield(Field_1,'VarAttribute') - if numel(Field_1.VarAttribute)>=VarType_1.coord(2) - XAttr=Field_1.VarAttribute{VarType_1.coord(2)} ; - else - XAttr=[]; - end - if numel(Field_1.VarAttribute)>=VarType_1.coord(1) - YAttr=Field_1.VarAttribute{VarType_1.coord(1)} ; - else - YAttr=[]; - end - SubData.VarAttribute=[SubData.VarAttribute {YAttr} {XAttr}]; - end - end - %look for previously used variable names - XName_1_1=XName_1;%default - YName_1_1=YName_1;%default - AName_1_1=AName_1;%default - for iprev=1:numel(SubData.ListVarName) - switch SubData.ListVarName{iprev} - case XName_1 - XName_1_1=[XName_1 '_1']; - case YName_1 - YName_1_1=[YName_1 '_1']; - case AName_1 - AName_1_1=[AName_1 '_1']; - end - end - if ~testX_1% if the second field has structured coordinates - ilist_1=find(strcmp(AName_1,Field_1.ListVarName)); - DimCell_1=Field_1.VarDimName{ilist_1}; - if numel(DimCell_1)>2 - DimCell={XName_1_1,YName_1_1, [{YName_1_1,XName_1_1} DimCell_1(3:end)]}; - else - DimCell={XName_1_1,YName_1_1, {YName_1_1,XName_1_1}}; - end - else - DimCell=[DimCell Field_1.VarDimName(ivar_C_1)]; - end - SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}]; - % check that a different dimension name is used for the two fields - if testX_1% if the second field has unstructured coordinates - for icell=1:numel(DimCell) - if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1}) - DimCell{icell}{1}=[DimCell{icell}{1} '_1']; - end - end - end - SubData.VarDimName=[SubData.VarDimName DimCell]; - if isfield(Field_1,'VarAttribute') - SubData.VarAttribute=[SubData.VarAttribute {AAttr}]; - end - eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';']) - eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';']) - eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';']) end - -%merge a scalar as the first input and a vector field as second input -if ~testU && testU_1 - UName_1=Field_1.ListVarName{VarType_1.vector_x}; - VName_1=Field_1.ListVarName{VarType_1.vector_y}; - UAttr=Field_1.VarAttribute{VarType_1.vector_x}; - VAttr=Field_1.VarAttribute{VarType_1.vector_y}; - if testX_1 %unstructured coordinate for the second field - XName_1=Field_1.ListVarName{VarType_1.coord_x}; - YName_1=Field_1.ListVarName{VarType_1.coord_y}; - - XAttr=Field_1.VarAttribute{VarType_1.coord_x}; - YAttr=Field_1.VarAttribute{VarType_1.coord_y}; -% SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}]; - DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]); - else - XName_1=Field_1.ListVarName{VarType_1.coord(2)}; - YName_1=Field_1.ListVarName{VarType_1.coord(1)}; - if numel(Field_1.VarAttribute)>=VarType_1.coord(2) - XAttr=Field_1.VarAttribute{VarType_1.coord(2)}; - else - XAttr=[]; - end - if numel(Field_1.VarAttribute)>=VarType_1.coord(1) - YAttr=Field_1.VarAttribute{VarType_1.coord(1)}; - else - YAttr=[]; - end - end - %check for the existence of the same variable name - XName_1_1=XName_1; %default - YName_1_1=YName_1; %default - UName_1_1=UName_1; %default - VName_1_1=VName_1; %default - for iprev=1:numel(SubData.ListVarName) - switch SubData.ListVarName{iprev} - case XName_1 - XName_1_1=[XName_1 '_1']; - case YName_1 - YName_1_1=[YName_1 '_1']; - case UName_1 - UName_1_1=[UName_1 '_1']; - case VName_1 - VName_1_1=[VName_1 '_1']; - end - end +for imod=1:numel(ivar) + VarName=SubData.ListVarName{ivar(imod)}; + VarName_1=SubData.ListVarName{ivar_1(imod)}; + SubData.(VarName)=SubData.(VarName)-SubData.(VarName_1); + ind_remove(ivar_1(imod))=1; +end +SubData.ListVarName(find(ind_remove))=[]; +SubData.VarDimName(find(ind_remove))=[]; +SubData.VarAttribute(find(ind_remove))=[]; + - if ~testX_1 - DimCell=[{XName_1_1} {YName_1_1}]; - end - SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}]; - DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])]; - SubData.VarDimName=[SubData.VarDimName DimCell]; - if isfield(SubData,'VarAttribute') - if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName)) - for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)-4 - SubData.VarAttribute{ivar}=[]; - end - end - SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}]; - end - SubData.(XName_1_1)=Field_1.(XName_1); - SubData.(YName_1_1)=Field_1.(YName_1); - SubData.(UName_1_1)=Field_1.(UName_1); - SubData.(VName_1_1)=Field_1.(VName_1); -end - Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 514) +++ /trunk/src/uvmat.m (revision 515) @@ -70,16 +70,12 @@ % (calc_tps.m) calculate tps coefficients (for filter projection or spatial derivatives). % | -% (calc_field.m) calculate field -% | % UvData.Field-------------->histogram % _____________|____________ % | | -% proj_field.m proj_field.m project the field on the projection objects +% proj_field.m proj_field.m project the field on the projection objects (use calc_field.m) % | | % UvData.PlotAxes ViewData.PlotAxes (on view_field) % | | % plot_field.m (uvmat) plot_field.m (view_field) plot the projected fields -% -% rmq: calc_field can be performed instead at the level of proj_field when needed % % @@ -209,5 +205,5 @@ %% TRANSFORM menu: builtin fcts -transform_menu={'';'phys';'px';'phys_polar'}; +transform_menu={'';'sub_field';'phys';'phys_polar'}; UvData.OpenParam.NbBuiltin=numel(transform_menu); %number of functions path_list=(num2cell(blanks(UvData.OpenParam.NbBuiltin)))';%initialize a cell array of nbvar blanks @@ -299,22 +295,22 @@ %% plot input field if exists -if testinputfield - %delete drawn objects - hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects - for iobj=1:length(hother) - delete_object(hother(iobj)) - end - if isempty(inputfile) - errormsg=refresh_field(handles,[],[],[],[],[],[],{Field}); - set(handles.MenuTools,'Enable','on') - set(handles.OBJECT_txt,'Visible','on') - set(handles.edit_object,'Visible','on') -% set(handles.ListObject_1,'Visible','on') - set(handles.frame_object,'Visible','on') - if ~isempty(errormsg) - msgbox_uvmat('ERROR',errormsg) - end - end -end +% if testinputfield +% %delete drawn objects +% hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects +% for iobj=1:length(hother) +% delete_object(hother(iobj)) +% end +% if isempty(inputfile) +% errormsg=refresh_field(handles,[],[],[],[],[],[],{Field}); +% set(handles.MenuTools,'Enable','on') +% set(handles.OBJECT_txt,'Visible','on') +% set(handles.edit_object,'Visible','on') +% % set(handles.ListObject_1,'Visible','on') +% set(handles.frame_object,'Visible','on') +% if ~isempty(errormsg) +% msgbox_uvmat('ERROR',errormsg) +% end +% end +% end set_vec_col_bar(handles) %update the display of color code for vectors @@ -679,9 +675,24 @@ set(handles.j1,'String',num2stra(j1,NomType)); set(handles.j2,'String',num2stra(j2,NomType)); - else %read the current field index if the second file series is opened - i1=str2num(get(handles.i1,'String')); - i2=str2num(get(handles.i2,'String')); - j1=stra2num(get(handles.j1,'String')); - j2=stra2num(get(handles.j2,'String')); + else %read the current field index to synchronise with the first series + i1_s=str2num(get(handles.i1,'String')); + i2_0=str2num(get(handles.i2,'String')); + if ~isempty(i2_0) + i2_s=i2_0; + else + i2_s=i2; + end + j1_0=stra2num(get(handles.j1,'String')); + if ~isempty(j1_0) + j1_s=j1_0; + else + j1_s=j1; + end + j2_0=stra2num(get(handles.j2,'String')); + if ~isempty(j2_0) + j2_s=j2_0; + else + j2_s=j2; + end end @@ -705,11 +716,12 @@ end %updtate the indices of the second field series to correspond to the newly opened one - FileName_1=fullfile_uvmat(Input.RootPath_1,Input.SubDir_1,Input.RootFile_1,Input.FileExt_1,Input.NomType_1,i1,i2,j1,j2); + FileName_1=fullfile_uvmat(Input.RootPath_1,Input.SubDir_1,Input.RootFile_1,Input.FileExt_1,Input.NomType_1,i1_s,i2_s,j1_s,j2_s); if exist(FileName_1,'file') + FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1_s,i2_s,j1_s,j2_s); + else FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1,i2,j1,j2); - set(handles.FileIndex_1,'String',FileIndex_1) - else msgbox_uvmat('WARNING','unable to synchronise the indices of the two series') end + set(handles.FileIndex_1,'String',FileIndex_1) end @@ -757,4 +769,5 @@ % --- Update information about a new field series (indices to scan, timing, % calibration from an xml file, then refresh current plots + function update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileType,VideoObject,index) %------------------------------------------------------------------------ @@ -943,5 +956,5 @@ end if ~get(handles.CheckFixLimits,'Value') - set(handles.transform_fct,'Value',2); % phys transform by default if fixedLimits is off + set(handles.transform_fct,'Value',3); % phys transform by default if fixedLimits is off end if isfield(GeometryCalib,'SliceCoord') @@ -1014,9 +1027,9 @@ end end -[ref_j,ref_i]=find(i1_series); +[ref_j,ref_i]=find(squeeze(i1_series(1,:,:))); if ~isempty(j1_series) state_j='on'; if index==1 - if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first vzlue + if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first value scan_option='j'; %scan j indext end @@ -1057,6 +1070,9 @@ %% apply the effect of the transform fct and view the field +transform=get(handles.path_transform,'UserData'); +if index==2 && (~isa(transform,'function_handle')||nargin(transform)<3) + set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields +end transform_fct_Callback([],[],handles) -%run0_Callback([],[], handles); %view field mask_test=get(handles.CheckMask,'value'); if mask_test @@ -1106,39 +1122,63 @@ function i1_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ -set(handles.i1,'BackgroundColor',[0.7 0.7 0.7]) +update_ij(handles,1) + +%------------------------------------------------------------------------ +function i2_Callback(hObject, eventdata, handles) +%------------------------------------------------------------------------ +update_ij(handles,2) + +%------------------------------------------------------------------------ +function j1_Callback(hObject, eventdata, handles) +%------------------------------------------------------------------------ +update_ij(handles,3) + +%------------------------------------------------------------------------ +function j2_Callback(hObject, eventdata, handles) +%------------------------------------------------------------------------ +update_ij(handles,4) + +%------------------------------------------------------------------------ +%--- update the index display after action on edit boxes i1, i2, j1 or j2 +function update_ij(handles,index_rank) NomType=get(handles.NomType,'String'); -i1=stra2num(get(handles.i1,'String')); -i2=stra2num(get(handles.i2,'String')); -j1=stra2num(get(handles.j1,'String')); -j2=stra2num(get(handles.j2,'String')); -indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2); -%indices=name_generator('',num1,num_a,'',NomType,1,num2,num_b,''); +indices=get(handles.FileIndex,'String'); +[tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(indices);% the indices for the second series taken from FileIndex +switch index_rank + case 1 + indices=fullfile_uvmat('','','','',NomType,stra2num(get(handles.i1,'String')),i2,j1,j2); + set(handles.i1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation + case 2 + indices=fullfile_uvmat('','','','',NomType,i1,stra2num(get(handles.i2,'String')),j1,j2); + set(handles.i2,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation + case 3 + indices=fullfile_uvmat('','','','',NomType,i1,i2,stra2num(get(handles.j1,'String')),j2); + set(handles.j1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation + case 4 + indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,stra2num(get(handles.j2,'String'))); + set(handles.j2,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation +end set(handles.FileIndex,'String',indices) -set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7]) -if get(handles.SubField,'Value')==1 +set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation +% update the second index if relevant +if strcmp(get(handles.FileIndex_1,'Visible'),'on') NomType_1=get(handles.NomType_1,'String'); - indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2); - %indices=name_generator('',num1,num_a,'',NomType_1,1,num2,num_b,''); - set(handles.FileIndex_1,'String',indices) - set(handles.FileIndex_1,'BackgroundColor',[0.7 0.7 0.7]) -end - -%------------------------------------------------------------------------ -function i2_Callback(hObject, eventdata, handles) -%------------------------------------------------------------------------ -set(handles.i2,'BackgroundColor',[0.7 0.7 0.7]) -i1_Callback(hObject, eventdata, handles) - -%------------------------------------------------------------------------ -function j1_Callback(hObject, eventdata, handles) -%------------------------------------------------------------------------ -set(handles.j1,'BackgroundColor',[0.7 0.7 0.7]) -i1_Callback(hObject, eventdata, handles) - -%------------------------------------------------------------------------ -function j2_Callback(hObject, eventdata, handles) -%------------------------------------------------------------------------ -set(handles.j2,'BackgroundColor',[0.7 0.7 0.7]) -i1_Callback(hObject, eventdata, handles) + indices_1=get(handles.FileIndex_1,'String'); + [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(indices_1);% the indices for the second series taken from FileIndex_1 + switch index_rank + case 1 + indices_1=fullfile_uvmat('','','','',NomType_1,stra2num(get(handles.i1,'String')),i2_1,j1_1,j2_1); + case 2 + indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,stra2num(get(handles.i2,'String')),j1_1,j2_1); + case 3 + indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,stra2num(get(handles.j1,'String')),j2_1); + case 4 + indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,j1_1,stra2num(get(handles.j2,'String'))); + end + set(handles.FileIndex_1,'String',indices_1) + set(handles.FileIndex_1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation +end + +%------------------------------------------------------------------------ %------------------------------------------------------------------------ @@ -1523,5 +1563,5 @@ if sub_value % a second input file has been entered [InputFile.RootPath_1,InputFile.SubDir_1,InputFile.RootFile_1,InputFile.FileIndex_1,InputFile.FileExt_1,InputFile.NomType_1]=read_file_boxes_1(handles); - [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(InputFile.FileIndex_1); + [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(InputFile.FileIndex_1);% the indices for the second series taken from FileIndex_1 else filename_1=[]; @@ -1529,7 +1569,5 @@ %% increment (or decrement) the field indices and update the input filename(s) -CheckSearch=0; if isempty(increment) - CheckSearch=1;% search for the next available file set(handles.CheckFixPair,'Value',0) end @@ -1545,4 +1583,5 @@ i2_1=i2_1+increment; end + else % case of scanning along index j (burst numbers) j1=j1+increment; @@ -1604,9 +1643,9 @@ if ~isempty(errormsg),return,end siz=size(UvData.i1_series{1}); - ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1); + ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1); i1_subseries=UvData.i1_series{1}(ref_indices); ref_indices=ref_indices(i1_subseries>0); if isempty(ref_indices)% case of pairs (free index i) - ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2); + ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2); i1_subseries=UvData.i1_series{1}(ref_indices); ref_indices=ref_indices(i1_subseries>0); @@ -1623,25 +1662,69 @@ if ~isempty(UvData.j2_series{1}) j2=UvData.j2_series{1}(ref_indices(end)); - end - % case of a second file series (TODO:revise) - if numel(UvData.i1_series)>=2 - i1_subseries=UvData.i1_series{2}(ref_i+1,ref_j+1,:); - i1_subseries=i1_subseries(i1_subseries>0); - i1_1=i1_subseries(end); + end + + % case of a second file series + if sub_value + ref_i_1=i1_1; + if ~isempty(i2_1) + ref_i_1=floor((i1_1+i2_1)/2);% current reference index i + end + ref_j_1=1; + if ~isempty(j1_1) + ref_j_1=j1_1; + if ~isempty(j2_1) + ref_j_1=floor((j1_1+j2_1)/2);% current reference index j + end + end + if ~isempty(increment) + if get(handles.scan_i,'Value')==1% case of scanning along index i + ref_i_1=ref_i_1+increment;% increment the current reference index i + else % case of scanning along index j (burst numbers) + ref_j_1=ref_j_1+increment;% increment the current reference index j if scan_j option is used + end + else % free increment, synchronise the ref indices with the first series + ref_i_1=ref_i; + ref_j_1=ref_j; + end + if numel(UvData.i1_series)==1 + UvData.i1_series{2}=UvData.i1_series{1}; + UvData.j1_series{2}=UvData.j1_series{1}; + UvData.i2_series{2}=UvData.i2_series{1}; + UvData.j2_series{2}=UvData.j2_series{1}; + end + if ref_i_1<0 + errormsg='minimum i index reached'; + elseif ref_j_1<0 + errormsg='minimum j index reached'; + elseif ref_i_1+1>size(UvData.i1_series{2},3) + errormsg='maximum i index reached for the second series (reload the input file to update the index bound)'; + elseif ref_j_1+1>size(UvData.i1_series{2},2) + errormsg='maximum j index reached for the second series(reload the input file to update the index bound)'; + end + if ~isempty(errormsg),return,end + siz=size(UvData.i1_series{2}); + ref_indices=ref_i_1*siz(1)*siz(2)+ref_j_1*siz(1)+1:ref_i_1*siz(1)*siz(2)+(ref_j_1+1)*siz(1); + i1_subseries=UvData.i1_series{2}(ref_indices); + ref_indices=ref_indices(i1_subseries>0); + if isempty(ref_indices)% case of pairs (free index i) + ref_indices=ref_i_1*siz(1)*siz(2)+1:(ref_i_1+1)*siz(1)*siz(2); + i1_subseries=UvData.i1_series{2}(ref_indices); + ref_indices=ref_indices(i1_subseries>0); + end + i1_1=UvData.i1_series{2}(ref_indices(end)); if ~isempty(UvData.i2_series{2}) - i2_subseries=UvData.i2_series{2}(ref_i+1,ref_j+1,:); - i2_subseries=i2_subseries(i2_subseries>0); - i2_1=i2_subseries(end); + i2_1=UvData.i2_series{2}(ref_indices(end)); end if ~isempty(UvData.j1_series{2}) - j1_subseries=UvData.j1_series{2}(ref_i+1,ref_j+1,:); - j1_subseries=j1_subseries(j1_subseries>0); - j1_1=j1_subseries(end); + j1_1=UvData.j1_series{2}(ref_indices(end)); end if ~isempty(UvData.j2_series{2}) - j2_subseries=UvData.j2_series{2}(ref_i+1,ref_j+1,:); - j2_subseries=j2_subseries(j2_subseries>0); - j2_1=j2_subseries(end); - end + j2_1=UvData.j2_series{1}(ref_indices(end)); + end + else% the second series (if needed) is the same file as the first + i1_1=i1; + i2_1=i2; + j1_1=j1; + j2_1=j2; end end @@ -1652,5 +1735,5 @@ %% refresh plots -errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2); +errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2,i1_1,i2_1,j1_1,j2_1); %% update the index counters if the index move is successfull @@ -1846,10 +1929,11 @@ set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow drawnow -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); -filename=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; -filename_1=[];%default +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); +filename=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt]; +filename_1='';%default +FileIndex_1=''; if get(handles.SubField,'Value') - [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles); - filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1]; + [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles); + filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1]; end num_i1=stra2num(get(handles.i1,'String')); @@ -1857,6 +1941,7 @@ num_j1=stra2num(get(handles.j1,'String')); num_j2=stra2num(get(handles.j2,'String')); - -errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2); +[tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(FileIndex_1); + +errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1); if ~isempty(errormsg) @@ -1883,5 +1968,5 @@ % Field: structure describing an optional input field (then replace the input file) -function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,Field) +function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1) %------------------------------------------------------------------------ @@ -1923,4 +2008,15 @@ if length(masknumber)>=z_index set(handles.masklevel,'Value',z_index) +end + +%% test for need of tps +check_proj_tps=0; +if (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx')) + for iobj=1:numel(UvData.Object) + if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter') + check_proj_tps=1; + break + end + end end @@ -1977,4 +2073,8 @@ ParamIn.VelType=VelType; ParamIn.GUIName='get_field'; +end +check_tps=0; +if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamIn.FieldName,'velocity')&&~strcmp(ParamIn.FieldName,'get_field...') + check_tps=1;%tps needed to get the requested field end [Field{1},ParamOut,errormsg] = read_field(FileName,UvData.FileType{1},ParamIn,frame_index); @@ -2030,13 +2130,13 @@ ParamIn_1=UvData.MovieObject{2}; if ~strcmp(NomType_1,'*') - frame_index_1=num_j1;%frame index for movies or multimage + frame_index_1=j1_1;%frame index for movies or multimage else - frame_index_1=num_i1; + frame_index_1=i1_1; end case 'multimage' if ~strcmp(NomType_1,'*') - frame_index_1=num_j1;%frame index for movies or multimage + frame_index_1=j1_1;%frame index for movies or multimage else - frame_index_1=num_i1; + frame_index_1=i1_1; end case 'vol' %TODO: update @@ -2055,8 +2155,6 @@ end test_keepdata_1=0;% test for keeping the previous stored data if the input files are unchanged - if ~isequal(NomType_1,'*')%in case of a series of files (not avi movie) - if isfield(UvData,'FileName_1')%&& isfield(UvData,'VelType_1') && isfield(UvData,'FieldName_1') - test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1); - end + if ~isequal(NomType_1,'*')&& isfield(UvData,'FileName_1') + test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1); end if test_keepdata_1 @@ -2068,9 +2166,15 @@ ParamIn_1.VelType=VelType_1; ParamIn_1.GUIName='get_field_1'; - end + end [Field{2},ParamOut_1,errormsg] = read_field(Name,UvData.FileType{2},ParamIn_1,frame_index_1); if ~isempty(errormsg) errormsg=['error in reading ' FieldName_1 ' in ' FileName_1 ': ' errormsg]; return + end + if isstruct(ParamOut_1)&&~strcmp(ParamOut_1.FieldName,'get_field...')&& (strcmp(UvData.FileType{2},'civdata')||strcmp(UvData.FileType{2},'civx'))... + &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...') + if ~check_proj_tps + % Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2}); + end end end @@ -2130,5 +2234,5 @@ test_veltype_1=1; set(handles.VelType_1,'Visible','on') - menu=set_veltype_display(ParamOut_1.CivStage); + menu=set_veltype_display(ParamOut_1.CivStage,UvData.FileType{2}); index_menu=strcmp(ParamOut_1.VelType,menu); set(handles.VelType_1,'Value',1+find(index_menu,1)) @@ -2219,18 +2323,17 @@ end if numel(UvData.XmlData)==2 - [tild,tild,tild,num_i1,num_i2,num_j1,num_j2]=fileparts_uvmat(['xx' get(handles.FileIndex_1,'String') get(handles.FileExt_1,'String')]); - if isempty(num_i2) - num_i2=num_i1; - end - if isempty(num_j1) - num_j1=1; - end - if isempty(num_j2) - num_j2=num_j1; + if isempty(i2_1) + i2_1=num_i1; + end + if isempty(j1_1) + j1_1=1; + end + if isempty(j2_1) + j2_1=j1_1; end siz=size(UvData.XmlData{2}.Time); - if ~isempty(num_i1) && siz(1)>=max(num_i1,num_i2) && siz(2)>=max(num_j1,num_j2) - abstime_1=(UvData.XmlData{2}.Time(num_i1,num_j1)+UvData.XmlData{2}.Time(num_i2,num_j2))/2;%overset the time read from files - Field{2}.Dt=(UvData.XmlData{2}.Time(num_i2,num_j2)-UvData.XmlData{2}.Time(num_i1,num_j1)); + if ~isempty(i1_1) && siz(1)>=max(i1_1,i2_1) && siz(2)>=max(j1_1,j2_1) + abstime_1=(UvData.XmlData{2}.Time(i1_1,j1_1)+UvData.XmlData{2}.Time(i2_1,j2_1))/2;%overset the time read from files + Field{2}.Dt=(UvData.XmlData{2}.Time(i2_1,j2_1)-UvData.XmlData{2}.Time(i1_1,j1_1)); end end @@ -2264,5 +2367,7 @@ %% apply coordinate transform or other user fct transform=get(handles.path_transform,'UserData'); -if ~isempty(transform) +if isempty(transform) + UvData.Field=Field{1}; +else XmlData=[];%default XmlData_1=[];%default @@ -2294,39 +2399,6 @@ end -%% check whether tps is needed, then calculate tps coefficients if needed -check_proj_tps=0; -if (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx')) - for iobj=1:numel(UvData.Object) - if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter') - check_proj_tps=1; - break - end - end -end -check_tps=0; -if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamOut.FieldName,'velocity')&&~strcmp(ParamOut.FieldName,'get_field...') - check_tps=1;%tps needed to get the requested field -end -if (check_tps ||check_proj_tps)&&~isfield(UvData.Field,'Coord_tps') - UvData.Field=calc_tps(UvData.Field); -end -UvData.Field.FieldList=[{ParamOut.FieldName} {ParamOut.ColorVar}]; - -%% calculate scalar -if isstruct(ParamOut)&&~strcmp(ParamOut.FieldName,'get_field...')&& (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))... - &&~strcmp(ParamOut.FieldName,'velocity') && ~strcmp(ParamOut.FieldName,'get_field...') - if ~check_proj_tps - UvData.Field=calc_field([{ParamOut.FieldName} {ParamOut.ColorVar}],UvData.Field); - end -end -if isstruct(ParamOut_1)&& numel(Field)==2 && ~strcmp(ParamOut_1.FieldName,'get_field...')&& ~test_keepdata_1 && (strcmp(UvData.FileType{2},'civdata')||strcmp(UvData.FileType{2},'civx'))... - &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...') - if check_proj_tps - Field{2}.FieldList=[{ParamOut_1.FieldName} {ParamOut_1.ColorVar}]; - else - Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2}); - end -end - +%% calculate tps coefficients if needed +UvData.Field=calc_tps(UvData.Field,check_proj_tps); %% analyse input field @@ -2942,5 +3014,12 @@ end set(handles.uvmat,'UserData',UvData); - run0_Callback(hObject, eventdata, handles); %run + transform_fct_list=get(handles.transform_fct,'String'); + transform_fct=transform_fct_list(get(handles.transform_fct,'Value')); + if strcmp(transform_fct,'sub_field') + set(handles.transform_fct,'Value',1)%suppress the sub_field transform + transform_fct_Callback(hObject, eventdata, handles); + else + run0_Callback(hObject, eventdata, handles) + end else MenuBrowse_1_Callback(hObject, eventdata, handles) @@ -3114,6 +3193,4 @@ FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1]; [tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(get(handles.FileIndex,'String')); -% set(handles.FileIndex_1,'Visible','on') -% set(handles.FileExt_1,'Visible','on') switch field_1 case 'get_field...' @@ -3128,10 +3205,7 @@ set(hhget_field.list_fig,'Value',1) set(hhget_field.list_fig,'String',{'uvmat'}) -% set(handles.transform_fct,'Value',1)% no transform by default -% set(handles.path_transform,'String','') if check_new UvData.FileType{2}=UvData.FileType{1}; set(handles.FileIndex_1,'String',get(handles.FileIndex,'String')) -% set(handles.FileExt_1,'String',get(handles.FileExt,'String')) set(handles.uvmat,'UserData',UvData) end @@ -3239,51 +3313,60 @@ %------------------------------------------------------------------------ -% --- Executes on button press in VelType_1. +% --- Executes on choice selection in VelType_1. function VelType_1_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ set(handles.FixVelType,'Value',1)% the velocity type is now imposed by the GUI (not automatic) UvData=get(handles.uvmat,'UserData'); -%refresh field with a second FileName=first file name -set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow +set(handles.run0,'BackgroundColor',[1 1 0])%paint run0 button in yellow to indicate its activation drawnow -InputFile=read_GUI(handles.InputFile); -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); -FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; - +InputFile=read_GUI(handles.InputFile);% read the input file parameters +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); +[RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles); +FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];% name of the first input file + +check_refresh=0; if isempty(InputFile.VelType_1) - FileName_1='';% we plot the current field without the second field + FileName_1='';% we plot the first input field without the second field set(handles.SubField,'Value',0) - SubField_Callback(hObject, eventdata, handles) + SubField_Callback(hObject, eventdata, handles)% activate SubField_Callback and refresh current display, removing the second field elseif get(handles.SubField,'Value')% if subfield is already 'on' - [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles); - FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1]; -else - FileName_1=FileName;% we compare two fields in the same file by default + FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];% name of the second input file + check_refresh=1;%will refresh the current display +else% we introduce the same file (with a different field) for the second series + FileName_1=FileName;% we compare two fields in the same file UvData.FileType{2}=UvData.FileType{1}; + UvData.XmlData{2}= UvData.XmlData{1}; set(handles.SubField,'Value',1) -end -if isfield(UvData,'Field_1') - UvData=rmfield(UvData,'Field_1');% removes the stored second field if it exists -end -UvData.FileName_1='';% desactivate the use of a constant second file -set(handles.uvmat,'UserData',UvData) -num_i1=stra2num(get(handles.i1,'String')); -num_i2=stra2num(get(handles.i2,'String')); -num_j1=stra2num(get(handles.j1,'String')); -num_j2=stra2num(get(handles.j2,'String')); - -errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2); - -if ~isempty(errormsg) - msgbox_uvmat('ERROR',errormsg); -else - set(handles.i1,'BackgroundColor',[1 1 1]) - set(handles.i2,'BackgroundColor',[1 1 1]) - set(handles.j1,'BackgroundColor',[1 1 1]) - set(handles.j2,'BackgroundColor',[1 1 1]) - set(handles.FileIndex,'BackgroundColor',[1 1 1]) - set(handles.FileIndex_1,'BackgroundColor',[1 1 1]) -end -set(handles.run0,'BackgroundColor',[1 0 0]) + transform=get(handles.path_transform,'UserData'); + if (~isa(transform,'function_handle')||nargin(transform)<3) + set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields + transform_fct_Callback(hObject, eventdata, handles)% activate transform_fct_Callback and refresh current display + else + check_refresh=1; + end +end + +% refresh the current display if it has not been done previously +if check_refresh + UvData.FileName_1='';% desactivate the use of a constant second file + set(handles.uvmat,'UserData',UvData) + num_i1=stra2num(get(handles.i1,'String')); + num_i2=stra2num(get(handles.i2,'String')); + num_j1=stra2num(get(handles.j1,'String')); + num_j2=stra2num(get(handles.j2,'String')); + [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(['xx' FileIndex_1]); + errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1); + if ~isempty(errormsg) + msgbox_uvmat('ERROR',errormsg); + else + set(handles.i1,'BackgroundColor',[1 1 1]) + set(handles.i2,'BackgroundColor',[1 1 1]) + set(handles.j1,'BackgroundColor',[1 1 1]) + set(handles.j2,'BackgroundColor',[1 1 1]) + set(handles.FileIndex,'BackgroundColor',[1 1 1]) + set(handles.FileIndex_1,'BackgroundColor',[1 1 1]) + end + set(handles.run0,'BackgroundColor',[1 0 0]) +end %----------------------------------------------- @@ -3605,7 +3688,10 @@ %% refresh the current plot -run0_Callback(hObject, eventdata, handles) - - +if isempty(list_path{ichoice}) || nargin(transform_handle)<3 + set(handles.SubField,'Value',0) + SubField_Callback(hObject, eventdata, handles) +else + run0_Callback(hObject, eventdata, handles) +end %------------------------------------------------ @@ -3961,8 +4047,7 @@ %replot all the objects within the new projected field for IndexObj=1:numel(list_str) - IndexObj - hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat + hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat; if isempty(hobject) || ~ishandle(hobject) - hobject=handles.PlotAxes + hobject=handles.PlotAxes; end if isequal(IndexObj,get(handles.ListObject,'Value')) @@ -4151,23 +4236,23 @@ data.Type='plane'; end - if isfield(UvData,'Field') - Field=UvData.Field; - if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh) - data.RangeX=[UvData.Field.XMin UvData.Field.XMax]; - if strcmp(data.Type,'line')||strcmp(data.Type,'polyline') - data.RangeY=UvData.Field.Mesh; - else - data.RangeY=[UvData.Field.YMin UvData.Field.YMax]; - end - data.DX=UvData.Field.Mesh; - data.DY=UvData.Field.Mesh; - end - if isfield(Field,'NbDim')&& isequal(Field.NbDim,3) - data.Coord=[0 0 0]; %default - end - if isfield(Field,'CoordUnit') - data.CoordUnit=Field.CoordUnit; - end - end +% if isfield(UvData,'Field') +% Field=UvData.Field; +% if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh) +% data.RangeX=[UvData.Field.XMin UvData.Field.XMax]; +% if strcmp(data.Type,'line')||strcmp(data.Type,'polyline') +% data.RangeY=UvData.Field.Mesh; +% else +% data.RangeY=[UvData.Field.YMin UvData.Field.YMax]; +% end +% data.DX=UvData.Field.Mesh; +% data.DY=UvData.Field.Mesh; +% end +% if isfield(Field,'NbDim')&& isequal(Field.NbDim,3) +% data.Coord=[0 0 0]; %default +% end +% if isfield(Field,'CoordUnit') +% data.CoordUnit=Field.CoordUnit; +% end +% end hset_object=set_object(data,[],ZBounds); hhset_object=guidata(hset_object); @@ -4297,5 +4382,5 @@ UvData=get(huvmat,'UserData'); %[xx,xx,FileBase]=read_file_boxes(handles); -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); FileBase=fullfile(RootPath,RootFile); %read the current input file name @@ -4587,6 +4672,6 @@ pos(1)=pos(1)+pos(3)-0.311+0.04; %0.311= width of the geometry_calib interface (units relative to the srcreen) pos(2)=pos(2)-0.02; -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); -FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); +FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt]; set(handles.view_xml,'Backgroundcolor',[1 1 0])%indicate the reading of the current xml file by geometry_calib if isfield(UvData.OpenParam,'CalOrigin') @@ -4720,6 +4805,6 @@ %prepare display of the set_grid GUI -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); -FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); +FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt]; CoordList=get(handles.transform_fct,'String'); val=get(handles.transform_fct,'Value'); @@ -4746,9 +4831,9 @@ %------------------------------------------------------------------------ series; %first display of the GUI to fill waiting time -[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); -param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; +[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); +param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt]; if isequal(get(handles.SubField,'Value'),1) - [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles); - FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1]; + [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles); + FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1]; if ~isequal(FileName_1,param.FileName) param.FileName_1=FileName_1; @@ -4788,6 +4873,6 @@ function MenuPIV_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ - [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles); - FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt]; + [RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles); + FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt]; civ(FileName);% interface de civ(not in the uvmat file)