Changeset 515


Ignore:
Timestamp:
Aug 15, 2012, 11:36:12 PM (9 years ago)
Author:
sommeria
Message:

improvement of calc-field and combination of two fields

Location:
trunk/src
Files:
2 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/calc_field.m

    r501 r515  
    3838
    3939%% list of field options implemented
    40 FieldOptions={'velocity';...%image correlation corresponding to a vel vector
    41     'ima_cor';...%image correlation corresponding to a vel vector
    42     'norm_vel';...%norm of the velocity
    43     'vort';...%vorticity
    44     'div';...%divergence
    45     'strain';...%rate of strain
    46     'u';... %u velocity component
    47     'v';... %v velocity component
     40FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector
     41    'C';...%image correlation corresponding to a vel vector
     42    'norm(U,V)';...%norm of the velocity
     43    'curl(U,V)';...%vorticity
     44    'div(U,V)';...%divergence
     45    'strain(U,V)';...%rate of strain
     46    'U';... %u velocity component
     47    'V';... %v velocity component
    4848    'w';... %w velocity component
    4949    'w_normal';... %w velocity component normal to the plane
     
    7777end
    7878FieldList=FieldList(check_calc==1);
    79 % if isempty(FieldList)
    80 %     DataOut=DataIn;
    81 %     return
    82 % end
    8379if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
    8480    nbcoord=3;
     
    9187units_cell={};
    9288
     89%% reproduce global attributes
     90DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
     91for ilist=1:numel(DataOut.ListGlobalAttribute)
     92    DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
     93end
     94
    9395%% interpolation with new civ data
    94 if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
    95     %reproduce global attributes
    96     DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
    97     for ilist=1:numel(DataOut.ListGlobalAttribute)
    98         DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
    99     end
    100 
     96[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
     97icell_tps=[];
     98for icell=1:numel(CellVarIndex)
     99    VarType=VarTypeCell{icell};
     100    if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps)     
     101        icell_tps=[icell_tps icell];
     102    end
     103end
     104
     105%if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
     106if ~isempty(icell_tps)
    101107    %create a default grid if needed
    102108    if  ~exist('Coord_interp','var')
    103             XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
    104             YMax=max(max(DataIn.SubRange(2,:,:)));
    105             XMin=min(min(DataIn.SubRange(1,:,:)));
    106             YMin=min(min(DataIn.SubRange(2,:,:)));
     109        for ind=1:numel(icell_tps)
     110            SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps});
     111            XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates
     112            YMax(ind)=max(max(SubRange(2,:,:)));
     113            XMin(ind)=min(min(SubRange(1,:,:)));
     114            YMin(ind)=min(min(SubRange(2,:,:)));
     115        end
     116        XMax=max(XMax);
     117        YMax=max(YMax);
     118        XMin=min(XMin);
     119        YMin=min(YMin);
    107120        if ~isfield(DataIn,'Mesh')
    108121            DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps));
     
    119132        coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to  visualisation
    120133        coord_y=YMin:DataIn.Mesh:YMax;
    121 %         npx=length(coord_x);
    122 %         npy=length(coord_y);
    123134        DataOut.coord_x=[XMin XMax];
    124135        DataOut.coord_y=[YMin YMax];
    125136        [XI,YI]=meshgrid(coord_x,coord_y);
    126 %         XI=reshape(XI,[],1);
    127 %         YI=reshape(YI,[],1);
    128137        Coord_interp=cat(3,XI,YI);%[XI YI];
    129138    end
     
    131140    npy=size(Coord_interp,1);
    132141    Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3));
    133 %         npy=length(coord_y);
     142
    134143    %initialise output
    135144    nb_sites=size(Coord_interp,1);
     
    152161   
    153162    % interpolate data in each subdomain
    154     for isub=1:NbSubDomain
    155         nbvec_sub=DataIn.NbSites(isub);
    156         check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
    157         ind_sel=find(sum(check_range,2)==nb_coord);
    158         %rho smoothing parameter
    159         %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
    160         %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
    161         nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
    162         if check_grid
    163             EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
    164         end
    165         if check_der
    166             [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'
    167         end
    168         ListVar={};
    169         for ilist=1:length(FieldList)
    170             var_count=numel(ListVar);
    171             switch FieldList{ilist}
    172                 case 'velocity'
    173                     ListVar=[ListVar {'U', 'V'}];
    174                     VarAttribute{var_count+1}.Role='vector_x';
    175                     VarAttribute{var_count+2}.Role='vector_y';
    176                     DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    177                     DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    178                 case 'u'
    179                     ListVar=[ListVar {'u'}];
    180                     VarAttribute{var_count+1}.Role='scalar';
    181                     DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    182                 case 'v'
    183                     ListVar=[ListVar {'v'}];
    184                     VarAttribute{var_count+1}.Role='scalar';
    185                     DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    186                 case 'norm_vel'
    187                     ListVar=[ListVar {'norm_vel'}];
    188                     VarAttribute{var_count+1}.Role='scalar';
    189                     U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    190                     V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    191                     DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
    192                 case 'vort'
    193                     ListVar=[ListVar {'vort'}];
    194                     VarAttribute{var_count+1}.Role='scalar';
    195                     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);
    196                 case 'div'
    197                     ListVar=[ListVar {'div'}];
    198                     VarAttribute{var_count+1}.Role='scalar';
    199                     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);
    200                 case 'strain'
    201                     ListVar=[ListVar {'strain'}];
    202                     VarAttribute{var_count+1}.Role='scalar';
    203                     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);
     163    for icell=icell_tps
     164        for isub=1:NbSubDomain
     165            nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps});
     166            SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps});
     167            check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
     168            ind_sel=find(sum(check_range,2)==nb_coord);
     169            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     170            Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps});
     171            if check_grid
     172                EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
    204173            end
    205         end
    206     end
    207     DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
    208     nbval(nbval==0)=1;% to avoid division by zero for averaging
    209     if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
    210         DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
    211         DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
    212         DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
    213     end
    214     DataOut.ListVarName=[DataOut.ListVarName ListVar];
    215     for ifield=1:numel(ListVar)
    216         VarDimName{ifield}={'coord_y','coord_x'};
    217         DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
    218         DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
    219     end
    220     DataOut.FF=reshape(DataOut.FF,npy,npx);
    221     DataOut.VarDimName=[DataOut.VarDimName VarDimName];     
    222     DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
     174            if check_der
     175                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
     176            end
     177            ListVar={};
     178            U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps});
     179            V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps});
     180            for ilist=1:length(FieldList)
     181                var_count=numel(ListVar);
     182                switch FieldList{ilist}
     183                    case 'velocity'
     184                        ListVar=[ListVar {'U', 'V'}];
     185                        VarAttribute{var_count+1}.Role='vector_x';
     186                        VarAttribute{var_count+2}.Role='vector_y';
     187                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     188                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     189                    case 'u'
     190                        ListVar=[ListVar {'u'}];
     191                        VarAttribute{var_count+1}.Role='scalar';
     192                        DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     193                    case 'v'
     194                        ListVar=[ListVar {'v'}];
     195                        VarAttribute{var_count+1}.Role='scalar';
     196                        DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     197                    case 'norm_vel'
     198                        ListVar=[ListVar {'norm_vel'}];
     199                        VarAttribute{var_count+1}.Role='scalar';
     200                        U=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
     201                        V=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
     202                        DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
     203                    case 'vort'
     204                        ListVar=[ListVar {'vort'}];
     205                        VarAttribute{var_count+1}.Role='scalar';
     206                        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);
     207                    case 'div'
     208                        ListVar=[ListVar {'div'}];
     209                        VarAttribute{var_count+1}.Role='scalar';
     210                        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);
     211                    case 'strain'
     212                        ListVar=[ListVar {'strain'}];
     213                        VarAttribute{var_count+1}.Role='scalar';
     214                        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);
     215                end
     216            end
     217        end
     218        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
     219        nbval(nbval==0)=1;% to avoid division by zero for averaging
     220        if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
     221            DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
     222            DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
     223            DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
     224        end
     225        DataOut.ListVarName=[DataOut.ListVarName ListVar];
     226        for ifield=1:numel(ListVar)
     227            VarDimName{ifield}={'coord_y','coord_x'};
     228            DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
     229            DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
     230        end
     231        DataOut.FF=reshape(DataOut.FF,npy,npx);
     232        DataOut.VarDimName=[DataOut.VarDimName VarDimName];
     233        DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
     234    end
    223235else
    224236
  • trunk/src/calc_tps.m

    r508 r515  
    1 function DataOut=calc_tps(DataIn)     
     1function DataOut=calc_tps(DataIn,checkall)     
    22DataOut=DataIn;%default
    33SubDomain=1000; %default, estimated nbre of vectors in a subdomain used for tps
     
    55    SubDomain=DataIn.SubDomain;%
    66end
    7 [DataOut.SubRange,DataOut.NbSites,DataOut.Coord_tps,DataOut.U_tps,DataOut.V_tps] =...
    8     filter_tps([DataIn.X(DataIn.FF==0) DataIn.Y(DataIn.FF==0)],DataIn.U(DataIn.FF==0),DataIn.V(DataIn.FF==0),[],SubDomain,0);
    9 nbvar=numel(DataIn.ListVarName);
    10 DataOut.ListVarName=[DataIn.ListVarName {'SubRange','NbSites','Coord_tps','U_tps','V_tps'}];
    11 DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
    12     {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
    13 DataOut.VarAttribute{nbvar+3}.Role='coord_tps';
    14 DataOut.VarAttribute{nbvar+4}.Role='vector_x';
    15 DataOut.VarAttribute{nbvar+5}.Role='vector_y';
    16 if isfield(DataOut,'ListDimName')%cleaning
    17     DataOut=rmfield(DataOut,'ListDimName');
     7[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
     8nbtps=0;
     9for icell=1:numel(CellVarIndex);
     10    VarType=VarTypeCell{icell};
     11    if NbDimVec(icell)>=2 && ~isempty(VarType.coord_x)
     12        nbtps=nbtps+1;
     13        X=DataIn.(DataIn.ListVarName{VarType.coord_x});
     14        Y=DataIn.(DataIn.ListVarName{VarType.coord_y});
     15        if ~isempty(VarType.vector_x)&&~isempty(VarType.vector_y)
     16            Attr=DataIn.VarAttribute{VarType.vector_x};
     17            if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'FieldRequest')&&strcmp(Attr.FieldRequest,'derivatives')))               
     18                U=DataIn.(DataIn.ListVarName{VarType.vector_x});
     19                V=DataIn.(DataIn.ListVarName{VarType.vector_y});
     20            else
     21                continue
     22            end
     23        end
     24        if ~isempty(VarType.errorflag)
     25            FF=DataIn.(DataIn.ListVarName{VarType.errorflag});
     26            X=X(FF==0);
     27            Y=Y(FF==0);
     28            U=U(FF==0);
     29            V=V(FF==0);
     30        end
     31        if nbtps==1
     32            ListNewVar={'SubRange','NbSites','Coord_tps','U_tps','V_tps'};
     33            ListNewDim={'nb_tps','nb_subdomain'};
     34            DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
     35                {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
     36        else
     37            ListNewVar={['SubRange_' num2str(nbtps-1)],['NbSites_' num2str(nbtps-1)],['Coord_tps_' num2str(nbtps-1)],['U_tps_' num2str(nbtps-1)] ,['V_tps_' num2str(nbtps-1)]};
     38            ListNewDim={['nb_tps_' num2str(nbtps-1)],['nb_subdomain_' num2str(nbtps-1)]};
     39            DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds',ListNewDim{2}},ListNewDim(2),...
     40                {ListNewDim{1},'nb_coord',ListNewDim{2}},ListNewDim,ListNewDim}];
     41        end
     42        DataOut.ListVarName=[DataIn.ListVarName ListNewVar];
     43       
     44        [DataOut.(ListNewVar{1}),DataOut.(ListNewVar{2}),DataOut.(ListNewVar{3}),DataOut.(ListNewVar{4}),DataOut.(ListNewVar{5})] =...
     45            filter_tps([X Y],U,V,[],SubDomain,0);
     46        nbvar=numel(DataIn.ListVarName);
     47       
     48        DataOut.VarAttribute{nbvar+3}.Role='coord_tps';
     49        DataOut.VarAttribute{nbvar+4}=DataIn.VarAttribute{VarType.vector_x};%reproduce attributes of velocity
     50         DataOut.VarAttribute{nbvar+4}.Role='vector_x_tps';
     51         DataIn.VarAttribute{VarType.vector_x}.VarIndex_tps=nbvar+4;% indicte the correspondance with initial data
     52        DataOut.VarAttribute{nbvar+5}=DataIn.VarAttribute{VarType.vector_y};%reproduce attributes of velocity
     53         DataOut.VarAttribute{nbvar+5}.Role='vector_y_tps';
     54%          if isfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest')
     55%              DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest');
     56%          end
     57%          if isfield(DataOut.VarAttribute{VarType.vector_x},'Operation')
     58%              DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'Operation');
     59%          end
     60        if isfield(DataOut,'ListDimName')%cleaning'FieldRequest'
     61            DataOut=rmfield(DataOut,'ListDimName');
     62        end
     63        if isfield(DataOut,'DimValue')%cleaning
     64            DataOut=rmfield(DataOut,'DimValue');
     65        end
     66    end
    1867end
    19 if isfield(DataOut,'DimValue')%cleaning
    20     DataOut=rmfield(DataOut,'DimValue');
    21 end
  • trunk/src/check_files.m

    r512 r515  
    3232svn_info.status=[];
    3333list_fct={...
    34     'calc_field';...% defines fields (velocity, vort, div...) from civx data and calculate them
     34    'calc_field_interp';...% defines fields (velocity, vort, div...) from civx data and calculate them
     35    'calc_field_tps';...% defines fields (velocity, vort, div...) and calculate them
    3536    'cell2tab';... %transform a Matlab cell in a character array suitable for display in a table
    3637    'check_field_structure';...% check the validity of the field struture representation consistant with the netcdf format
     
    6465    'get_file_series';...% determine the list of file names and file indices for functions called by 'series'.
    6566    'get_file_type';...% determine info about a file (image, multimage, civdata,...) .
    66     'griddata_uvmat';...%make 2D linear interpolation using griddata, with input appropriate for both Matlab 6.5 and 7
    6767    'hist_update';...%  update of a current global histogram by inclusion of a new field
    6868    'imadoc2struct';...%convert the image documentation file ImaDoc into a Matlab structure
  • trunk/src/fill_GUI.m

    r500 r515  
    3131        hh=[];
    3232        input_data=Param.(fields{ifield});
    33                     display(fields{ifield})
    34                     display(input_data)
     33%                     display(fields{ifield})
     34%                     display(input_data)
    3535        check_done=0;
    3636        if isfield(handles,fields{ifield})
     
    6868                    end
    6969                case 'edit'
     70                    input_string='';
    7071                    if isnumeric(input_data)
     72                        if numel(input_data)>0
    7173                        input_string=num2str(input_data(ibox));
     74                        end
    7275                    else
    7376                        input_string=input_data;
  • trunk/src/find_field_cells.m

    r514 r515  
    1919%      .scalar: scalar field (default)
    2020%      .coord: vector of indices of coordinate variables corresponding to matrix dimensions
     21%
     22%      .FieldRequest= 'interp', 'derivatives' indicate whether interpolation  or derivatives (tps) is needed to calculate the requested field
     23%      .FieldNames = cell of fields to calculate from the fied cell
     24%
    2125% errormsg: error message
    2226%   
     
    4953%     GNU General Public License (file UVMAT/COPYING.txt) for more details.for input in proj_field and plot_field
    5054%    group the variables  into 'fields' with common dimensions
    51 %------------------------------------------------------------------------
    52 % function  [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
    53 %
    54 % OUTPUT:
    55 % CellVaxIndex: cell whose elements are arrays of indices in the list data.ListVarName 
    56 %              CellvarIndex{i} represents a set of variables with the same dimensions
    57 % NbDim: array with the length of CellVarIndex, giving its  space dimension
    58 % CellVarType: cell array of structures with fields
    59 %      .coord_x, y, z: indices (in .ListVarname) of variables representing  unstructured coordinates x, y, z
    60 %      .vector_x,_y,_z: indices of variables giving the vector components x, y, z
    61 %      .warnflag: index of warnflag
    62 %      .errorflag: index of error flag
    63 %      .ancillary: indices of ancillary variables
    64 %      .image   : B/W image, (behaves like scalar)
    65 %      .color : color image, the last index, which is not a coordinate variable, represent the 3 color components rgb
    66 %      .discrete: like scalar, but set of data points without continuity, represented as dots in a usual plot, instead of continuous lines otherwise
    67 %      .scalar: scalar field (default)
    68 %      .coord: vector of indices of coordinate variables corresponding to matrix dimensions
    69 % errormsg: error message
    70 %   
    71 % INPUT:
    72 % Data: structure representing fields, output of check_field_structure
    73 %            .ListVarName: cell array listing the names (cahr strings) of the variables
    74 %            .VarDimName: cell array of cells containing the set of dimension names for each variable of .ListVarName
    75 %            .VarAttribute: cell array of structures containing the variable attributes:
    76 %                     .VarAttribute{ilist}.key=value, where ilist is the variable
    77 %                     index, key is the name of the attribute, value its value (char string or number)
    78 %
    79 % HELP:
    80 % to get the dimensions of arrays common to the field #icell
    81 %         VarIndex=CellVarIndex{icell}; % list of variable indices
    82 %         DimIndex=Data.VarDimIndex{VarIndex(1)} % list of dimensions for each variable in the cell #icell
    83 %
    84 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    85 %  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
    86 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    87 %     This file is part of the toolbox UVMAT.
    88 %
    89 %     UVMAT is free software; you can redistribute it and/or modify
    90 %     it under the terms of the GNU General Public License as published by
    91 %     the Free Software Foundation; either version 2 of the License, or
    92 %     (at your option) any later version.
    93 %
    94 %     UVMAT is distributed in the hope that it will be useful,
    95 %     but WITHOUT ANY WARRANTY; without even the implied warranty of
    96 %     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    97 %     GNU General Public License (file UVMAT/COPYING.txt) for more details.
    98 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    9955
    10056function [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
    10157CellVarIndex={};
    10258
    103 NbDim=[];
    10459CellVarType=[];
    10560errormsg=[];
     
    10863
    10964NbDim=[];
    110 ivardim=0;
    11165VarDimIndex=[];
    11266VarDimName={};
     
    11670end
    11771
    118 %% role of variables
     72%% role of variables and list of requested operations
    11973Role=num2cell(blanks(nbvar));%initialize a cell array of nbvar blanks
     74FieldRequest=regexprep(Role,' ',''); % fieldRequest set to '' by default
     75Operation=cell(size(Role)); % fieldRequest set to {} by default
     76CheckSub=zeros(size(Role));% =1 for fields to substract
    12077Role=regexprep(Role,' ','scalar'); % Role set to 'scalar' by default
    12178if isfield(Data,'VarAttribute')
     
    12481            Role{ivar}=Data.VarAttribute{ivar}.Role;
    12582        end
     83        if isfield(Data.VarAttribute{ivar},'FieldRequest')
     84            FieldRequest{ivar}=Data.VarAttribute{ivar}.FieldRequest;
     85        end
     86        if isfield(Data.VarAttribute{ivar},'Operation')
     87            Operation{ivar}=Data.VarAttribute{ivar}.Operation;
     88        end
     89        if isfield(Data.VarAttribute{ivar},'CheckSub')
     90            CheckSub(ivar)=Data.VarAttribute{ivar}.CheckSub;
     91        end
    12692    end
    12793end
    12894
    12995%% loop on the list of variables, group them by common dimensions
     96CellVarType=cell(1,length(CellVarIndex));
    13097for ivar=1:nbvar
    13198    if ischar(Data.VarDimName{ivar})
     
    146113        icell=icell+1;
    147114        CellVarIndex{icell}=ivar;%put the current variable index in the new cell
    148         NbDim(icell)=numel(DimCell);%default
    149     end
    150    
    151     %look for dimension variables
    152 %     if numel(DimCell)==1% if the variable has a single dimension
    153 %         if strcmp(DimCell{1},Data.ListVarName{ivar}) %|| strcmp(Role{ivar},'dimvar')
    154 %             ivardim=ivardim+1;
    155 %             VarDimIndex(ivardim)=ivar;%index of the variable
    156 %             VarDimName{ivardim}=DimCell{1};%name of the dimension
    157 %         end
    158 %     end
     115        NbDim(icell)=numel(DimCell);%default   
     116        CellVarType{icell}=[];
     117    end
     118    if ~isempty(FieldRequest{ivar})
     119       CellVarType{icell}.FieldRequest=FieldRequest{ivar};
     120    end
     121    if ~isempty(Operation{ivar})
     122       CellVarType{icell}.Operation=Operation{ivar};
     123    end
     124    if CheckSub(ivar)
     125    CellVarType{icell}.CheckSub=1;
     126    end
    159127end
    160128
     
    163131ind_dim_var_cell=find(checksinglecell);
    164132%CoordType(ind_dim_var_cell)='dim_var';% to be used in output
     133%VarDimIndex=cell(size(ind_dim_var_cell));
     134VarDimName=cell(size(ind_dim_var_cell));
    165135for icoord=1:numel(ind_dim_var_cell)
    166 VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)};
    167 VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1};
     136    VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)};
     137    VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1};
    168138end
    169139
     
    201171
    202172index_remove=[];
    203 CellVarType=cell(1,length(CellVarIndex));
    204173for icell=1:length(CellVarIndex)
    205174    if checksinglecell(icell)
     
    262231                    coord(idim)=VarDimIndex(ind_coord);
    263232                end
    264 %                 for ivardim=1:numel(VarDimName)
    265 %                     if strcmp(VarDimName{ivardim},DimCell{idim})
    266 %                         coord(idim)=VarDimIndex(ivardim);
    267 %                         break
    268 %                     end
    269 %                 end
    270233            end
    271234            NbDim(icell)=numel(find(coord));
  • trunk/src/plot_field.m

    r512 r515  
    122122    return
    123123end
    124 index_2D=find(NbDim==2,2);%find 2D fields (at most 2)
     124index_2D=find(NbDim==2);%find 2D fields
    125125index_3D=find(NbDim>2,1);
    126126if ~isempty(index_3D)
     
    560560    if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected
    561561        if test_vec
    562             errormsg='error in plot_field: attempt to plot two vector fields';
     562            errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with mode interp';
    563563            return
    564564        else
  • trunk/src/proj_field.m

    r512 r515  
    8484ProjData=[];
    8585
    86 %% case of no projection (object is used only as graph display)
    87 if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
     86%% in the absence of object Type or projection mode, or object coordinaes, the output is empty
     87if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
    8888    return
    8989end
    90 
    91 %% check coincidence of coordinate units
    92 if isfield(FieldData,'CoordUnit') && isfield(ObjectData,'CoordUnit')&&~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit)
    93     errormsg='inconsistent coord units for field and projection object';
     90% case of no projection (object is used only as graph display)
     91if isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside')
    9492    return
    9593end
    96    
    97 %% in the absence of object Type or projection mode, or object coordinaes, the input field is just tranfered without change
    98 if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
    99     ProjData=FieldData;
    100     return
    101 end
    102 if ~isfield(ObjectData,'Coord')
     94if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord)
    10395    if strcmp(ObjectData.Type,'plane')
    10496        ObjectData.Coord=[0 0 0];%default
    10597    else
    106         ProjData=FieldData;
    10798        return
    10899    end
    109100end
    110        
    111 %% OBSOLETE
    112 if isfield(ObjectData,'XMax') && ~isempty(ObjectData.XMax)
    113     ObjectData.RangeX(1)=ObjectData.XMax;
    114 end
    115 if isfield(ObjectData,'XMin') && ~isempty(ObjectData.XMin)
    116     ObjectData.RangeX(2)=ObjectData.XMin;
    117 end
    118 if isfield(ObjectData,'YMax') && ~isempty(ObjectData.YMax)
    119     ObjectData.RangeY(1)=ObjectData.YMax;
    120 end
    121 if isfield(ObjectData,'YMin') && ~isempty(ObjectData.YMin)
    122     ObjectData.RangeY(2)=ObjectData.YMin;
    123 end
    124 if isfield(ObjectData,'ZMax') && ~isempty(ObjectData.ZMax)
    125     ObjectData.RangeZ(1)=ObjectData.ZMax;
    126 end
    127 if isfield(ObjectData,'ZMin') && ~isempty(ObjectData.ZMin)
    128     ObjectData.RangeZ(2)=ObjectData.ZMin;
    129 end
    130 %%%%%%%%%%
    131101
    132102%% apply projection depending on the object type
     
    401371    testproj(VarType.color)=1;
    402372    VarIndex=VarIndex(find(testproj(VarIndex)));%select only the projected variables
    403     if testX %case of unstructured coordinates
     373    if check_unstructured%case of unstructured coordinates
    404374         eval(['nbpoint=numel(FieldData.' FieldData.ListVarName{VarIndex(1)} ');'])
    405375         for ivar=[VarIndex VarType.coord_x VarType.coord_y VarType.errorflag]
     
    547517ProjData.NbDim=1;
    548518%initialisation of the input parameters and defaultoutput
    549 ProjMode='projection';%direct projection on the line by default
    550 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
     519ProjMode=ObjectData.ProjMode;
    551520% ProjAngle=90; %90 degrees projection by default
    552521
     
    609578        continue
    610579    end
    611     testX=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
     580    check_unstructured=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
    612581    test_tps=~isempty(VarType.coord_tps);
    613582    testU=~isempty(VarType.vector_x) && ~isempty(VarType.vector_y);% test for vectors
     
    631600    end   
    632601    % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
    633     if testX
     602    if check_unstructured
    634603        if  ~isequal(ProjMode,'interp')
    635604            if width==0
     
    664633  %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    665634   %case of unstructured coordinates
    666     if testX   
     635    if check_unstructured   
    667636        for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
    668637            linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip)); 
     
    917886%-----------------------------------------------------------------
    918887
    919 %% initialisation of the input parameters of the projection plane
    920 ProjMode='projection';%direct projection by default
    921 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    922 
    923 %% axis origin
    924 if isempty(ObjectData.Coord)
    925     ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
    926     ObjectData.Coord(1,2)=0;
    927     ObjectData.Coord(1,3)=0;
    928 end
    929 
    930888%% rotation angles
    931889PlaneAngle=[0 0 0];
     
    951909
    952910%% mesh sizes DX and DY
    953 DX=0;
    954 DY=0; %default
     911DX=FieldData.Mesh;
     912DY=FieldData.Mesh; %default
    955913if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
    956914     DX=abs(ObjectData.DX);%mesh of interpolation points
     
    959917     DY=abs(ObjectData.DY);%mesh of interpolation points
    960918end
    961 if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0)
     919if  ~strcmp(ObjectData.ProjMode,'projection') && (DX==0||DY==0)
    962920        errormsg='DX or DY missing';
    963921        display(errormsg)
     
    970928testYMin=0;
    971929testYMax=0;
     930XMin=FieldData.XMin;%default
     931XMax=FieldData.XMax;%default
     932YMin=FieldData.YMin;%default
     933YMax=FieldData.YMax;%default
    972934if isfield(ObjectData,'RangeX')
    973935        XMin=min(ObjectData.RangeX);
     
    1003965ListIndex={};
    1004966
    1005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1006967%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
    1007 %-----------------------------------------------------------------
    1008 idimvar=0;
    1009 
     968% CellVarIndex=cells of variable index arrays
    1010969[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(FieldData);
    1011970if ~isempty(errormsg)
     
    1014973end
    1015974
     975%% projection modes
     976check_grid=0;
     977ProjMode=cell(size(VarTypeCell));
     978for icell=1:numel(VarTypeCell)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
     979    ProjMode{icell}=ObjectData.ProjMode;
     980    if isfield(VarTypeCell{icell},'FieldRequest')
     981        switch VarTypeCell{icell}.FieldRequest
     982            case 'interp'
     983                ProjMode{icell}='interp';
     984            case 'derivatives'
     985                ProjMode{icell}='filter';
     986        end
     987    end
     988    if strcmp(ProjMode{icell},'interp')||strcmp(ProjMode{icell},'filter')
     989        check_grid=1;
     990    end
     991end
     992
     993%% define the new coordinates in case of interpolation on a grid
     994if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates
     995    ProjData.ListVarName={'coord_y','coord_x'};
     996    ProjData.VarDimName={'coord_y','coord_x'}; 
     997    ProjData.VarAttribute={[],[]};
     998    ProjData.coord_y=[YMin YMax];
     999    ProjData.coord_x=[XMin XMax];
     1000end
     1001
     1002%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    10161003% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    10171004% CellVarIndex=cells of variable index arrays
     
    10201007nbcoord=0;%number of added coordinate variables brought by projection
    10211008nbvar=0;
     1009vector_x_proj=[];
     1010vector_y_proj=[];
    10221011for icell=1:length(CellVarIndex)
    10231012    NbDim=NbDimVec(icell);
     
    10271016    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
    10281017    VarType=VarTypeCell{icell};
    1029 %     ivar_X=VarType.coord_x;
    1030 %     ivar_Y=VarType.coord_y;
    1031 %     ivar_Z=VarType.coord_z;
    10321018    ivar_U=VarType.vector_x;
    10331019    ivar_V=VarType.vector_y;
     
    10371023    end
    10381024    ivar_W=VarType.vector_z;
    1039 %     ivar_Anc=VarType.ancillary;
    1040 %     test_anc=zeros(size(VarIndex));
    1041 %     test_anc(ivar_Anc)=ones(size(ivar_Anc));
    1042 %     ivar_F=VarType.warnflag;
    1043 %     ivar_FF=VarType.errorflag;
    10441025
    10451026    %type of coordinates
     
    10581039    end
    10591040    coord_z=0;%default
    1060    
     1041
    10611042    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    10621043    switch CoordType
     
    10641045        %% case of input fields with unstructured coordinates
    10651046        case 'unstructured'
    1066             if strcmp(ObjectData.ProjMode,'filter')
    1067                 continue %skip for filter (needs tps)
    1068             end
    1069             coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});
    1070             coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});
     1047            if strcmp(ProjMode{icell},'filter')
     1048                continue %skip for filter (needs tps field cell)
     1049            end
     1050            coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});% initial x coordinates
     1051            coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});% initial y coordinates
    10711052            if ~isempty(VarType.coord_z)
    1072                 coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});
     1053                coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});% initial z coordinates
    10731054            end
    10741055           
    1075             % translate  initial coordinates
     1056            % translate  initial coordinates to account for the new origin
    10761057            coord_x=coord_x-ObjectData.Coord(1,1);
    10771058            coord_y=coord_y-ObjectData.Coord(1,2);
     
    10851066                fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
    10861067                indcut=find(abs(fieldZ) <= width);
    1087                 size(indcut)
    10881068                for ivar=VarIndex
    10891069                    VarName=FieldData.ListVarName{ivar};
     
    10961076            end
    10971077           
    1098             %rotate coordinates if needed:
     1078            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    10991079            Psi=PlaneAngle(1);
    11001080            Theta=PlaneAngle(2);
     
    11041084                coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    11051085                coord_Y=coord_Y+coord_z *sin(Theta);
    1106                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1107                
     1086                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER               
    11081087                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    11091088            else
     
    11121091            end
    11131092           
    1114             %restriction to the range of x and y if imposed
     1093            %restriction to the range of X and Y if imposed by the projection object
    11151094            testin=ones(size(coord_X)); %default
    11161095            testbound=0;
     
    11391118                for ivar=VarIndex
    11401119                    VarName=FieldData.ListVarName{ivar};
    1141                     eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1120                    FieldData.(VarName)=FieldData.(VarName)(indcut);
    11421121                end
    11431122                coord_X=coord_X(indcut);
     
    11491128           
    11501129            % different cases of projection
    1151             switch ObjectData.ProjMode
    1152                 case 'projection'
     1130            switch ProjMode{icell}
     1131                case 'projection' 
     1132                    nbvar=numel(ProjData.ListVarName);
    11531133                    for ivar=VarIndex %transfer variables to the projection plane
    11541134                        VarName=FieldData.ListVarName{ivar};
     
    11721152                    coord_x_proj=XMin:DX:XMax;
    11731153                    coord_y_proj=YMin:DY:YMax;
    1174                     DimCell={'coord_y','coord_x'};
    1175                     ProjData.ListVarName={'coord_y','coord_x'};
    1176                     ProjData.VarDimName={'coord_y','coord_x'};
    1177                     nbcoord=2;
    1178                     ProjData.coord_y=[YMin YMax];
    1179                     ProjData.coord_x=[XMin XMax];
    1180 %                     if isempty(ivar_X), ivar_X=0; end;
    1181 %                     if isempty(ivar_Y), ivar_Y=0; end;
    1182 %                     if isempty(ivar_Z), ivar_Z=0; end;
    1183 %                     if isempty(ivar_U), ivar_U=0; end;
    1184 %                     if isempty(ivar_V), ivar_V=0; end;
    1185 %                     if isempty(ivar_W), ivar_W=0; end;
    1186 %                     if isempty(ivar_F), ivar_F=0; end;
    1187 %                     if isempty(ivar_FF), ivar_FF=0; end;
    1188                   %  if ~isequal(ivar_FF,0)
    1189                    if ~isempty(VarType.errorflag)
     1154                    [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);
     1155                    if ~isempty(VarType.errorflag)
    11901156                        VarName_FF=FieldData.ListVarName{VarType.errorflag};
    11911157                        indsel=find(FieldData.(VarName_FF)==0);
     
    11931159                        coord_Y=coord_Y(indsel);
    11941160                    end
    1195                    
    1196                     FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
    11971161                    testFF=0;
    1198                     %for ivar=VarIndex % loop on field variables to project
    1199                     for ivar=[VarType.scalar VarType.vector_x VarType.vector_y]
    1200                         VarName=FieldData.ListVarName{ivar};
    1201 %                         if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||...
    1202 %                                 ( numel(test_anc)>=ivar && test_anc(ivar)==1))
    1203                             ivar_new=ivar_new+1;
    1204                             ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    1205                             ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1206                             if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1207                                 ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
    1208                             end
    1209                             %if  ~isequal(ivar_FF,0)
    1210                             if ~isempty(VarType.errorflag)
    1211                                 FieldData.(VarName)=FieldData.(VarName)(indsel);
    1212                             end
    1213                             ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj');
    1214                             varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
    1215                             FFlag= isnan(varline); %detect undefined values NaN
    1216                             indnan=find(FFlag);
    1217                             if~isempty(indnan)
    1218                                 varline(indnan)=zeros(size(indnan));
    1219                                 ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj));
    1220                                 FF(indnan)=ones(size(indnan));
    1221                                 testFF=1;
    1222                             end
    1223                             if ivar==ivar_U
    1224                                 ivar_U=ivar_new;
    1225                             end
    1226                             if ivar==ivar_V
    1227                                 ivar_V=ivar_new;
    1228                             end
    1229                             if ivar==ivar_W
    1230                                 ivar_W=ivar_new;
    1231                             end
    1232 %                         end
    1233                     end
    1234                     if testFF
    1235                         ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    1236                         ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1237                         ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1238                         ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    1239                     end
    1240 %                     case 'filter'%interpolate data on a regular grid
    1241 %                         errormsg='tps required for filter option'
    1242                        
    1243             end
    1244            
    1245             %% case of tps interpolation (applies only in filter mode)
     1162                    nbvar=numel(ProjData.ListVarName);
     1163                    if isfield(VarType,'vector_x')&&isfield(VarType,'vector_y')&&~isempty(VarType.vector_x)
     1164                        VarName_x=FieldData.ListVarName{VarType.vector_x};
     1165                        VarName_y=FieldData.ListVarName{VarType.vector_y};
     1166                        if ~isempty(VarType.errorflag)
     1167                            FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
     1168                            FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
     1169                        end
     1170                        FieldVar=cat(2,FieldData.(VarName_x),FieldData.(VarName_y));
     1171                        if ~isfield(VarType,'CheckSub') || ~VarType.CheckSub
     1172                            vector_x_proj=numel(ProjData.ListVarName)+1;
     1173                            vector_y_proj=numel(ProjData.ListVarName)+2;
     1174                        end
     1175                    end
     1176                    if ~isempty(VarType.scalar)
     1177                        VarName_scalar=FieldData.ListVarName{VarType.scalar};
     1178                        if ~isempty(VarType.errorflag)
     1179                            FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
     1180                        end
     1181                        FieldVar=FieldData.(VarName_scalar);
     1182                    end
     1183                    [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldVar,VarType.Operation,XI,YI);
     1184                    if isfield(VarType,'CheckSub') && VarType.CheckSub && ~isempty(vector_x_proj)
     1185                        ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1};
     1186                        ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2};
     1187                    else
     1188                        VarDimName=cell(size(ListFieldProj));
     1189                        for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1190                            ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+','');
     1191                            if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName)))
     1192                                ListFieldProj{ilist}=[ListFieldProj{ilist} '_1'];
     1193                            end                       
     1194                            ProjData.(ListFieldProj{ilist})=VarVal{ilist};
     1195                            VarDimName{ilist}={'coord_y','coord_x'};
     1196                        end
     1197                        ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
     1198                        ProjData.VarDimName=[ProjData.VarDimName VarDimName];
     1199                        ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
     1200                    end
     1201            end
     1202
     1203            %% case of tps interpolation (applies only in filter mode and for spatial derivatives)
    12461204        case 'tps'
    1247             if strcmp(ObjectData.ProjMode,'filter')
     1205            if strcmp(ProjMode{icell},'filter')
     1206                Coord=FieldData.(FieldData.ListVarName{VarType.coord_tps});
     1207                NbSites=FieldData.(FieldData.ListVarName{VarType.nbsites_tps});
     1208                SubRange=FieldData.(FieldData.ListVarName{VarType.subrange_tps});
     1209                if isfield(VarType,'vector_x_tps')&&isfield(VarType,'vector_y_tps')
     1210                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{VarType.vector_x_tps}),FieldData.(FieldData.ListVarName{VarType.vector_y_tps}));
     1211                end
    12481212                coord_x_proj=XMin:DX:XMax;
    12491213                coord_y_proj=YMin:DY:YMax;
     
    12531217                XI=XI+ObjectData.Coord(1,1);
    12541218                YI=YI+ObjectData.Coord(1,2);
    1255                 DataOut=calc_field(FieldData.FieldList,FieldData,cat(3,XI,YI));
    1256                 ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
    1257                 ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
    1258                 ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];   
    1259                 for ilist=1:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    1260                     VarName=DataOut.ListVarName{ilist};
    1261 %                     ProjData.(VarName)=DataOut.(VarName);
    1262                     if ilist>=3
    1263                     ProjData.(VarName)=reshape(DataOut.(VarName),np_y,np_x);
    1264                     end
    1265                 end
    1266                 ProjData.coord_x=[XMin XMax];
    1267                 ProjData.coord_y=[YMin YMax];
     1219                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbSites,SubRange,FieldVar,VarType.Operation,cat(3,XI,YI));   
     1220                ListFieldProj=(fieldnames(DataOut))';
     1221                VarDimName=cell(size(ListFieldProj));
     1222                for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1223                    VarName=ListFieldProj{ilist};
     1224                    ProjData.(VarName)=DataOut.(VarName);
     1225                    VarDimName{ilist}={'coord_y','coord_x'};
     1226                end
     1227%                 if ~isfield(ProjData,'coord_x')
     1228%                 ProjData.coord_x=[XMin XMax];
     1229%                 ProjData.coord_y=[YMin YMax];
     1230%                 ListFieldProj=[{'coord_x','coord_y'} ListFieldProj'];
     1231%                 VarDimName=[{'coord_x','coord_y'} VarDimName];
     1232%                 VarAttribute=[{[],[]} VarAttribute];
     1233%                 end
     1234                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
     1235                ProjData.VarDimName=[ProjData.VarDimName VarDimName];
     1236                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    12681237            end
    12691238           
     
    12721241           
    12731242            VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
    1274             eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
     1243            DimValue=size(FieldData.(VarName));%input matrix dimensions
    12751244            DimValue(DimValue==1)=[];%remove singleton dimensions
    12761245            NbDim=numel(DimValue);%update number of space dimensions
     
    14051374            end
    14061375            % case with no  interpolation
    1407             if isequal(ProjMode,'projection') && (~testangle || test90y || test90x)
     1376            if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x)
    14081377                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax
    1409                     ProjData=FieldData;% no change by projection
     1378                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)];
     1379                    ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)]; 
     1380                    ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)];
     1381                    ProjData.(AYProjName)=FieldData.(AYName);
     1382                    ProjData.(AXProjName)=FieldData.(AXName);
     1383                    for ivar=VarIndex
     1384                        VarName=FieldData.ListVarName{ivar};
     1385                        ProjData.(VarName)=FieldData.(VarName);% no change by projection
     1386                    end
    14101387                else
    14111388                    indY=NbDim-1;
     
    14891466                    YIMA=reshape(round(YIMA),1,npX*npY);
    14901467                    flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    1491                     if isequal(ObjectData.ProjMode,'filter')
     1468                    if isequal(ProjMode{icell},'filter')
    14921469                        npx_filter=ceil(abs(DX/DAX));
    14931470                        npy_filter=ceil(abs(DY/DAY));
     
    15811558end
    15821559
     1560% %prepare substraction in case of two input fields
     1561% SubData.ListVarName={};
     1562% SubData.VarDimName={};
     1563% SubData.VarAttribute={};
     1564% check_remove=zeros(size(ProjData.ListVarName));
     1565% for iproj=1:numel(ProjData.VarAttribute)
     1566%     if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1)
     1567%         VarName=ProjData.ListVarName{iproj};
     1568%         SubData.ListVarName=[SubData.ListVarName {VarName}];
     1569%         SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}];
     1570%         SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}];
     1571%         SubData.(VarName)=ProjData.(VarName);
     1572%         check_remove(iproj)=1;       
     1573%     end
     1574% end
     1575% if ~isempty(find(check_remove))
     1576%     ind_remove=find(check_remove);
     1577%     ProjData.ListVarName(ind_remove)=[];
     1578%     ProjData.VarDimName(ind_remove)=[];
     1579%     ProjData.VarAttribute(ind_remove)=[];
     1580%     ProjData=sub_field(ProjData,[],SubData);
     1581% end   
     1582
    15831583%-----------------------------------------------------------------
    15841584%projection in a volume
     
    15861586%-----------------------------------------------------------------
    15871587ProjData=FieldData;%default output
    1588 
    1589 %% initialisation of the input parameters of the projection plane
    1590 ProjMode='projection';%direct projection by default
    1591 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    15921588
    15931589%% axis origin
  • trunk/src/read_civdata.m

    r497 r515  
    1010%                    .NbCoord: number of vector components
    1111%                    .NbDim: number of dimensions (=2 or 3)
    12 %                    .dt: time interval for the corresponding image pair
     12%                    .Dt: time interval for the corresponding image pair
    1313%                    .Time: absolute time (average of the initial image pair)
    1414%                    .CivStage: =0,
     
    6363end
    6464errormsg='';
     65if ischar(FieldNames), FieldNames={FieldNames}; end;
     66FieldRequest='vec';
     67for ilist=1:length(FieldNames)
     68    if ~isempty(FieldNames{ilist})
     69        switch FieldNames{ilist}
     70            case{'U','V','norm(U,V)'}
     71                FieldRequest='interp';
     72            case {'curl(U,V)','div(U,V)','strain(U,V)'}
     73                FieldRequest='derivatives';
     74        end
     75    end
     76end
    6577
    6678%% reading data
    67 [varlist,role,units,VelTypeOut]=varcivx_generator(FieldNames,VelType,CivStage);
     79[varlist,role,VelTypeOut]=varcivx_generator(FieldRequest,VelType,CivStage);
    6880if isempty(varlist)
    6981    erromsg=['error in read_civdata: unknow velocity type ' VelType];
     
    97109end
    98110Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'Dt','Time'}];
     111ivar_U_tps=[];
    99112var_ind=find(vardetect);
    100113for ivar=1:numel(var_ind)
    101114    Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
    102     Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)};
     115    Field.VarAttribute{ivar}.FieldRequest=FieldRequest;
     116    if strcmp(role{var_ind(ivar)},'vector_x')
     117        Field.VarAttribute{ivar}.Operation=FieldNames;
     118        ivar_U=ivar;
     119    end
     120    if strcmp(role{var_ind(ivar)},'vector_x_tps')
     121        Field.VarAttribute{ivar}.Operation=FieldNames;
     122        ivar_U_tps=ivar;
     123    end
     124%     Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)};
    103125    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
    104126end
     127if ~isempty(ivar_U_tps)
     128    Field.VarAttribute{ivar_U}.VarIndex_tps=ivar_U_tps;
     129end
     130% if strcmp(FieldRequest,'derivatives')% fields will be calculated from the tps
     131%     Field.VarAttribute{ivar_U_tps}.FieldRequest=FieldRequest;
     132%     Field.VarAttribute{ivar_U_tps}.Operation=FieldNames;
     133% else% fields will be calculated from the initial fields
     134%     Field.VarAttribute{ivar_U}.FieldRequest=FieldRequest;%
     135%     Field.VarAttribute{ivar_U}.Operation=FieldNames;
     136% end
     137
    105138Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'NbCoord','NbDim','TimeUnit','CoordUnit'}];
    106139% %% update list of global attributes
     
    121154%            if vel_type=[] or'*', a  priority choice is done, civ2 considered better than civ1 )
    122155
    123 function [var,role,units,vel_type_out,errormsg]=varcivx_generator(FieldNames,vel_type,CivStage)
     156function [var,role,vel_type_out,errormsg]=varcivx_generator(FieldRequest,vel_type,CivStage)
    124157
    125158%% default input values
    126159if ~exist('vel_type','var'),vel_type='';end;
    127160if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
    128 if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar
    129 if ischar(FieldNames), FieldNames={FieldNames}; end;
    130161errormsg='';
    131162
    132163%% select the priority order for automatic vel_type selection
    133 testder=0;
    134 testpatch=0;
    135 for ilist=1:length(FieldNames)
    136     if ~isempty(FieldNames{ilist})
    137         switch FieldNames{ilist}
    138             case{'u','v'}
    139                 testpatch=1;
    140             case {'vort','div','strain'}
    141                 testder=1;
    142         end
    143     end
    144 end
    145 if strcmp(vel_type,'civ2') && testder
     164if strcmp(vel_type,'civ2') && strcmp(FieldRequest,'derivatives')
    146165    vel_type='filter2';
    147 elseif strcmp(vel_type,'civ1') && testder
     166elseif strcmp(vel_type,'civ1') && strcmp(FieldRequest,'derivatives')
    148167    vel_type='filter1';
    149168end
     
    153172            vel_type='filter2';
    154173        case {4,5}% civ2 available but not filter2
    155             if testder% derivatives needed
     174            if strcmp(FieldRequest,'derivatives')% derivatives needed
    156175                vel_type='filter1';
    157176            else
     
    171190            'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
    172191        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
    173         units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
     192    %    units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
    174193    case 'filter1'
    175194        var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';...
     
    178197        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
    179198            'vector_y_tps','vector_z_tps','ancillary','ancillary'};
    180         units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
     199     %   units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
    181200    case 'civ2'
    182201        var={'X','Y','Z','U','V','W','C','F','FF';...
    183202            'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF'};
    184203        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
    185         units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
     204      %  units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
    186205    case 'filter2'
    187206        var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';...
     
    190209        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
    191210            'vector_y_tps','vector_z_tps','ancillary','ancillary'};
    192         units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
     211       % units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
     212end
     213if ~strcmp(FieldRequest,'derivatives')
     214    var=var(:,1:9);%suppress tps if not needed
    193215end
    194216vel_type_out=vel_type;
  • trunk/src/sub_field.m

    r512 r515  
    77% when scalar and vectors are combined, the fields are just merged in a single matlab structure for common visualisation
    88%-----------------------------------------------------------------------
    9 % function SubData=sub_field(Field,Field_1)
     9% function SubData=sub_field(Field,XmlData,Field_1)
    1010%
    1111% OUPUT:
     
    1616% Field_1:matlab structure representing the second field
    1717
    18 function [SubData,errormsg]=sub_field(Field,Field_1)
     18function SubData=sub_field(Field,XmlData,Field_1)
    1919
     20SubData=[];
     21if strcmp(Field,'*')
     22    return
     23end
     24if nargin<3
     25    SubData=Field;
     26    return
     27end
    2028%% global attributes
    21 test_attr=0;
     29SubData.ListGlobalAttribute={};%default
     30%transfer global attributes of Field
    2231if isfield(Field,'ListGlobalAttribute')
    2332    SubData.ListGlobalAttribute=Field.ListGlobalAttribute;
     
    2635        SubData.(AttrName)=Field.(AttrName);
    2736    end
    28     test_attr=1;
    2937end
     38%transfer global attributes of Field_1
    3039if isfield(Field_1,'ListGlobalAttribute')
    3140    for ilist=1:numel(Field_1.ListGlobalAttribute)
    32         test_1=1;
    3341        AttrName=Field_1.ListGlobalAttribute{ilist};
    34         if test_attr
    35             for i_prev=1:numel(Field.ListGlobalAttribute)
    36                 if isequal(Field.ListGlobalAttribute{i_prev},AttrName)
    37                     test_1=0; %attribute already written
    38                     eval(['Val=Field.' AttrName ';'])                 
    39                     eval(['Val_1=Field_1.' AttrName ';'])
    40                     if isequal(Val,Val_1)           
    41                         break% data already written
    42                     else
    43                         eval(['SubData.' AttrName '_1=Field_1.' AttrName ';'])
    44                     end
    45                 end
    46             end
     42        AttrNameNew=AttrName;
     43        while ~isempty(find(strcmp(AttrNameNew,SubData.ListGlobalAttribute)))&&~isequal(Field_1.(AttrNameNew),Field.(AttrNameNew))
     44            AttrNameNew=[AttrNameNew '_1'];
    4745        end
    48         if test_1
    49             eval(['SubData.' AttrName '=Field_1.' AttrName ';'])
     46        if ~isfield(Field,AttrName) || ~isequal(Field_1.(AttrName),Field.(AttrName))
     47        SubData.ListGlobalAttribute=[SubData.ListGlobalAttribute {AttrNameNew}];
     48        SubData.(AttrNameNew)=Field_1.(AttrName);
    5049        end
    5150    end
     
    5352
    5453%% variables
     54%reproduce variables of the first field and list its dimensions
    5555SubData.ListVarName=Field.ListVarName;
    5656SubData.VarDimName=Field.VarDimName;
     
    5858    SubData.VarAttribute=Field.VarAttribute;
    5959end
    60 %reproduce the first field Field by default
    61 for ivar=1:numel(Field.ListVarName)
    62    VarName=Field.ListVarName{ivar};
    63    SubData.(VarName)=Field.(VarName);
     60ListDimName={};
     61for ilist=1:numel(Field.ListVarName)
     62    VarName=Field.ListVarName{ilist};
     63    SubData.(VarName)=Field.(VarName);
     64    SubData.VarAttribute{ilist}.CheckSub=0;
     65    DimCell=Field.VarDimName{ilist};
     66    if ischar(DimCell)
     67        DimCell={DimCell};
     68    end
     69    for idim=1:numel(DimCell)
     70        if isempty(find(strcmp(DimCell{idim},ListDimName)))
     71            ListDimName=[ListDimName DimCell(idim)];
     72        end
     73    end
    6474end
    6575
    66 %% check the two input fields     
    67 [CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_cells(Field);
    68 if ~isempty(errormsg)
    69     errormsg=['invalid  first input to sub_field:' errormsg];
    70     return
    71 end
    72 [CellVarIndex_1,NbDim_1,VarTypeCell_1,errormsg]=find_field_cells(Field_1);
    73 if ~isempty(errormsg)
    74     errormsg=['invalid second input to sub_field:' errormsg];
    75     return
    76 end
    77 iselect=find(NbDim==2);
    78 for icell=iselect
    79     if ~isempty(VarTypeCell{icell}.coord_tps)
    80         NbDim(icell)=0;
     76%% field request
     77FieldRequest=cell(size(Field.ListVarName));
     78for ilist=1:numel(Field.VarAttribute)
     79    if isfield(Field.VarAttribute{ilist},'FieldRequest')
     80        FieldRequest{ilist}=Field.VarAttribute{ilist}.FieldRequest;
    8181    end
    8282end
    83 iselect=find(NbDim==2);
    84 if ~isequal(numel(iselect),1)
    85     errormsg='invalid  first input to sub_field: it must  contain a single 2D field cell';
    86     return
    87 end
    88 iselect_1=find(NbDim_1==2);
    89 for icell=iselect_1
    90     if ~isempty(VarTypeCell_1{icell}.coord_tps)
    91         NbDim_1(icell)=0;
     83FieldRequest_1=cell(size(Field_1.ListVarName));
     84for ilist=1:numel(Field_1.VarAttribute)
     85    if isfield(Field_1.VarAttribute{ilist},'FieldRequest')
     86        FieldRequest_1{ilist}=Field_1.VarAttribute{ilist}.FieldRequest;
    9287    end
    9388end
    94 iselect_1=find(NbDim_1==2);
    95 if ~isequal(numel(iselect_1),1)
    96     errormsg='invalid  second input to sub_field: it must  contain a single 2D field cell';
    97     return
    98 end
    99 VarType=VarTypeCell{iselect};
    100 VarType_1=VarTypeCell_1{iselect_1};
    101 testX=~isempty(VarType.coord_x)&& ~isempty(VarType.coord_y);%unstructured coordiantes
    102 testX_1=~isempty(VarType_1.coord_x)&& ~isempty(VarType_1.coord_y);%unstructured coordiantes
    103 testU=~isempty(VarType.vector_x)&& ~isempty(VarType.vector_y);%vector field
    104 testU_1=~isempty(VarType_1.vector_x)&& ~isempty(VarType_1.vector_y);%vector field
    105 testfalse_1=~isempty(VarType_1.errorflag);
    106 ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields
    107 if numel(ivar_C)>1
    108     errormsg='too many scalar fields in the first input of sub_field.m';
    109     return
    110 end
    111 ivar_C_1=[VarType_1.scalar VarType_1.image VarType_1.color VarType_1.ancillary]; %defines index (indices) for the scalar or ancillary fields
    112 if numel(ivar_C_1)>1
    113     errormsg='too many scalar fields in the second input of sub_field.m';
    114     return
     89
     90%% rename the dimensions of the second field if identical to those of the first
     91for ilist=1:numel(Field_1.VarDimName)
     92    DimCell=Field_1.VarDimName{ilist};
     93    if ischar(DimCell)
     94        DimCell={DimCell};
     95    end
     96    for idim=1:numel(DimCell)
     97        ind_dim=find(strcmp(DimCell{idim},ListDimName));
     98        if ~isempty(ind_dim)
     99            if ischar(Field_1.VarDimName{ilist})
     100                Field_1.VarDimName{ilist}=Field_1.VarDimName(ilist);
     101            end
     102            Field_1.VarDimName{ilist}{idim}=[ListDimName{ind_dim} '_1'];
     103        end
     104    end
    115105end
    116106
    117 %% substract two vector fields or two scalars
    118 if (testU && testU_1) || (~testU && ~testU_1)
    119    %check coincidence in positions
    120    %unstructured coordinates for the first field
    121    if testX 
    122        XName=Field.ListVarName{VarType.coord_x};
    123        YName=Field.ListVarName{VarType.coord_y};
    124        vec_X=Field.(XName);
    125        vec_Y=Field.(YName);
    126        nbpoints=numel(vec_X);
    127        vec_X=reshape(vec_X,nbpoints,1);
    128        vec_Y=reshape(vec_Y,nbpoints,1);
    129        if testX_1 %unstructured coordinates for the second field
    130             X_1_Name=Field_1.ListVarName{VarType_1.coord_x};
    131             Y_1_Name=Field_1.ListVarName{VarType_1.coord_y};
    132             eval(['vec_X_1=Field_1.' X_1_Name ';'])
    133             eval(['vec_Y_1=Field_1.' Y_1_Name ';'])
     107%look for coordinates common to Field in Field_1
     108ind_remove=zeros(size(Field_1.ListVarName));
     109for ilist=1:numel(Field_1.ListVarName)
     110    if ischar(Field_1.VarDimName{ilist})||numel(Field_1.VarDimName{ilist})==1
     111        OldDim=Field_1.VarDimName{ilist};
     112        if ischar(OldDim)
     113            OldDim=Field_1.VarDimName(ilist);
     114        end
     115        VarVal=Field_1.(Field_1.ListVarName{ilist});
     116        for i1=1:numel(Field.ListVarName)
     117            if (isempty(FieldRequest{i1})&&isempty(FieldRequest_1{ilist})||strcmp(FieldRequest{i1},FieldRequest_1{ilist})) && isequal(Field.(Field.ListVarName{i1}),VarVal)
     118               ind_remove(ilist)=1;
     119               NewDim=Field.VarDimName{i1};
     120               if ischar(NewDim)
     121                   NewDim={NewDim};
     122               end
     123               Field_1.VarDimName=regexprep_r(Field_1.VarDimName,OldDim{1},NewDim{1});
     124            end
     125        end
     126    end
     127end
     128Field_1.ListVarName(find(ind_remove))=[];%removes these redondent coordinates
     129Field_1.VarDimName(find(ind_remove))=[];
     130Field_1.VarAttribute(find(ind_remove))=[];
    134131
    135        else   %structured coordinates for the second field
    136            y_1_Name=Field_1.ListVarName{VarType_1.coord(1)};
    137            x_1_Name=Field_1.ListVarName{VarType_1.coord(2)};
    138            eval(['y_1=Field_1.' y_1_Name ';'])
    139            eval(['x_1=Field_1.' x_1_Name ';'])
    140            if isequal(numel(x_1),2) 
    141                x_1=linspace(x_1(1),x_1(2),nbpoints_x_1);
    142            end
    143            if isequal(numel(y_1),2) 
    144                y_1=linspace(y_1(1),y_1(2),nbpoints_y_1);
    145            end
    146            [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);
    147        end
    148        vec_X_1=reshape(vec_X_1,[],1);
    149        vec_Y_1=reshape(vec_Y_1,[],1);
    150        if testfalse_1
    151            FFName_1=Field_1.ListVarName{VarType_1.errorflag};         
    152            eval(['vec_FF_1=Field_1.' FFName_1 ';'])
    153            vec_FF_1=reshape(vec_FF_1,[],1);
    154            indsel=find(~vec_FF_1);
    155            vec_X_1=vec_X_1(indsel);
    156            vec_Y_1=vec_Y_1(indsel);
    157        end
    158        if testU % vector fields
    159             U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
    160             V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
    161             eval(['vec_U_1=Field_1.' U_1_Name ';'])
    162             eval(['vec_V_1=Field_1.' V_1_Name ';'])
    163             nbpoints_x_1=size(vec_U_1,2);
    164             nbpoints_y_1=size(vec_U_1,1);
    165             vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1);
    166             vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1);
    167             if testfalse_1
    168                 vec_U_1=vec_U_1(indsel);
    169                 vec_V_1=vec_V_1(indsel);
    170             end           
    171        else %(~testU && ~testU_1)
    172            A_1_Name=Field_1.ListVarName{ivar_C_1};
    173            eval(['vec_A_1=Field_1.' A_1_Name ';'])   
    174            nbpoints_x_1=size(vec_A_1,2);
    175            nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2)
    176            vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1);
    177            if testfalse_1
    178                 vec_A_1=vec_A_1(indsel);
    179            end
    180        end
    181 
    182        if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same
    183            if testU
    184                vec_U_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_U_1,vec_X,vec_Y);  %interpolate vectors in the second field
    185                vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y);  %interpolate vectors in the second field   
    186            else
    187                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
    188            end
    189        end
    190        if testU
    191            UName=Field.ListVarName{VarType.vector_x};
    192            VName=Field.ListVarName{VarType.vector_y}; 
    193            eval(['vec_U=Field.' UName ';'])
    194            eval(['vec_V=Field.' VName ';'])       
    195            vec_U=reshape(vec_U,numel(vec_U),1);
    196            vec_V=reshape(vec_V,numel(vec_V),1);
    197            eval(['SubData.' UName '=vec_U-vec_U_1;'])
    198            eval(['SubData.' VName '=vec_V-vec_V_1;'])
    199        else
    200            AName=Field.ListVarName{ivar_C};
    201           % size(Field.vort)
    202            eval(['SubData.' AName '=Field.' AName '-vec_A_1;'])
    203        end
    204    else  %structured coordiantes
    205        XName=Field.ListVarName{VarType.coord(2)};
    206        YName=Field.ListVarName{VarType.coord(1)};
    207        eval(['x=Field.' XName ';'])
    208        eval(['y=Field.' YName ';'])
    209        if testX_1 %unstructured coordinates for the second field
    210            errormsg='the second input scalar is not on a regular grid: comparison option not implemented';
    211            return
    212        else
    213            XName_1=Field.ListVarName{VarType_1.coord(2)};
    214            YName_1=Field.ListVarName{VarType_1.coord(1)};
    215            eval(['x_1=Field_1.' XName_1 ';'])
    216            eval(['y_1=Field_1.' YName_1 ';'])
    217        end
    218        if testU % vector fields
    219            UName=Field.ListVarName{VarType.vector_x};
    220            VName=Field.ListVarName{VarType.vector_y};
    221            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
    222            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
    223            eval(['U_1=Field_1.' U_1_Name ';'])
    224            eval(['V_1=Field_1.' V_1_Name ';'])
    225            if ~isequal(x_1,x)||~isequal(y_1,y)
    226                 [X_1,Y_1]=meshgrid(x_1,y_1);
    227                 U_1 =interp2(X_1,Y_1,U_1,x,y');
    228                 V_1 =interp2(X_1,Y_1,V_1,x,y');
    229            end
    230            eval(['SubData.' UName '=Field.' UName '-U_1;'])
    231            eval(['SubData.' VName '=Field.' VName '-V_1;'])
    232        else
    233            AName=Field.ListVarName{ivar_C};
    234            A_1_Name=Field_1.ListVarName{ivar_C_1};
    235            eval(['A_1=double(Field_1.' A_1_Name ');'])
    236            if ~isequal(x_1,x)||~isequal(y_1,y)
    237                 [X_1,Y_1]=meshgrid(x_1,y_1);
    238                 A_1 =interp2(X_1,Y_1,A_1,x,y');
    239            end
    240            eval(['SubData.' AName '=double(Field.' AName ')-A_1;'])
    241        end
    242    end
     132%append the other variables of the second field, modifying their name if needed
     133for ilist=1:numel(Field_1.ListVarName)
     134    VarName=Field_1.ListVarName{ilist};
     135    ind_prev=find(strcmp(VarName,Field.ListVarName));
     136    if isempty(ind_prev)% variable name does not exist in Field
     137        VarNameNew=VarName;
     138    else  % variable name exists in Field     
     139            VarNameNew=[VarName '_1'];   
     140    end
     141        SubData.ListVarName=[SubData.ListVarName {VarNameNew}];
     142        SubData.VarDimName=[SubData.VarDimName Field_1.VarDimName(ilist)];
     143        SubData.(VarNameNew)=Field_1.(VarName);
     144        SubData.VarAttribute=[SubData.VarAttribute Field_1.VarAttribute(ilist)];
     145        SubData.VarAttribute{end}.CheckSub=1;% mark that the field needs to be substracted
    243146end
    244147
    245 % merge a vector field and a scalar as second input
    246 if testU && ~testU_1
    247     AName_1=Field_1.ListVarName{ivar_C_1};
    248     if isfield(Field_1,'VarAttribute') & numel(Field_1.VarAttribute)>=ivar_C_1
    249         AAttr=Field_1.VarAttribute{ivar_C_1} ;
    250     else
    251         AAttr=[];
     148%append the other variables of the second field, modifying their name if needed
     149
     150%% substrat fields when possible
     151[CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(SubData);
     152ind_remove=zeros(size(SubData.ListVarName));
     153ivar=[];
     154ivar_1=[];
     155for icell=1:numel(CellVarIndex)
     156    if ~isempty(CellVarIndex{icell})
     157        if isfield(CellVarType{icell},'scalar') && numel(CellVarType{icell}.scalar)==2 && SubData.VarAttribute{CellVarType{icell}.scalar(2)}.CheckSub;
     158            ivar=[ivar CellVarType{icell}.scalar(1)];
     159            ivar_1=[ivar_1 CellVarType{icell}.scalar(2)];
     160        end
     161        if isfield(CellVarType{icell},'vector_x') && numel(CellVarType{icell}.vector_x)==2 && SubData.VarAttribute{CellVarType{icell}.vector_x(2)}.CheckSub;
     162            ivar=[ivar CellVarType{icell}.vector_x(1)];
     163            ivar_1=[ivar_1 CellVarType{icell}.vector_x(2)];
     164        end
     165        if isfield(CellVarType{icell},'vector_y') && numel(CellVarType{icell}.vector_y)==2 && SubData.VarAttribute{CellVarType{icell}.vector_y(2)}.CheckSub;
     166            ivar=[ivar CellVarType{icell}.vector_y(1)];
     167            ivar_1=[ivar_1 CellVarType{icell}.vector_y(2)];
     168        end
    252169    end
    253     if testX_1 %unstructured coordinate
    254        XName_1=Field_1.ListVarName{VarType_1.coord_x};
    255        YName_1=Field_1.ListVarName{VarType_1.coord_y};
    256        DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
    257        if isfield(Field_1,'VarAttribute')
    258            if numel(Field_1.VarAttribute)>=VarType_1.coord_x
    259                 XAttr=Field_1.VarAttribute{VarType_1.coord_x};
    260            else
    261                 XAttr=[];
    262            end
    263            if numel(Field_1.VarAttribute)>=VarType_1.coord_y
    264                YAttr=Field_1.VarAttribute{VarType_1.coord_y};
    265            else
    266                YAttr=[];
    267            end
    268            SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr}];
    269        end
    270     else
    271        XName_1=Field_1.ListVarName{VarType_1.coord(2)};
    272        YName_1=Field_1.ListVarName{VarType_1.coord(1)};
    273 %        DimCell=[{YName_1} {XName_1}];
    274        if isfield(Field_1,'VarAttribute')
    275            if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
    276                 XAttr=Field_1.VarAttribute{VarType_1.coord(2)} ;
    277            else
    278                 XAttr=[];
    279            end
    280            if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
    281                YAttr=Field_1.VarAttribute{VarType_1.coord(1)} ;
    282            else
    283                YAttr=[];
    284            end
    285            SubData.VarAttribute=[SubData.VarAttribute {YAttr} {XAttr}];
    286        end
    287     end 
    288     %look for previously used variable names
    289     XName_1_1=XName_1;%default
    290     YName_1_1=YName_1;%default
    291     AName_1_1=AName_1;%default
    292     for iprev=1:numel(SubData.ListVarName)
    293         switch SubData.ListVarName{iprev}
    294             case XName_1
    295                 XName_1_1=[XName_1 '_1'];
    296             case YName_1
    297                 YName_1_1=[YName_1 '_1'];
    298             case AName_1
    299                 AName_1_1=[AName_1 '_1'];
    300         end
    301     end     
    302     if ~testX_1% if the second field has structured coordinates
    303         ilist_1=find(strcmp(AName_1,Field_1.ListVarName));
    304         DimCell_1=Field_1.VarDimName{ilist_1};
    305         if numel(DimCell_1)>2
    306              DimCell={XName_1_1,YName_1_1, [{YName_1_1,XName_1_1} DimCell_1(3:end)]};
    307         else
    308           DimCell={XName_1_1,YName_1_1, {YName_1_1,XName_1_1}};
    309         end
    310     else
    311         DimCell=[DimCell Field_1.VarDimName(ivar_C_1)];
    312     end
    313     SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}];
    314     % check that a  different dimension name is used for the two fields
    315      if testX_1% if the second field has unstructured coordinates
    316         for icell=1:numel(DimCell)
    317             if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1})
    318                 DimCell{icell}{1}=[DimCell{icell}{1} '_1'];
    319             end
    320         end
    321      end
    322     SubData.VarDimName=[SubData.VarDimName DimCell];
    323     if isfield(Field_1,'VarAttribute')
    324         SubData.VarAttribute=[SubData.VarAttribute {AAttr}];
    325     end
    326     eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';'])
    327     eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';'])
    328     eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';'])
    329170end
    330 
    331 %merge a scalar as the first input and a vector field as second input
    332 if ~testU && testU_1
    333     UName_1=Field_1.ListVarName{VarType_1.vector_x};
    334     VName_1=Field_1.ListVarName{VarType_1.vector_y};
    335     UAttr=Field_1.VarAttribute{VarType_1.vector_x};
    336     VAttr=Field_1.VarAttribute{VarType_1.vector_y};
    337     if testX_1 %unstructured coordinate for the second field
    338        XName_1=Field_1.ListVarName{VarType_1.coord_x};
    339        YName_1=Field_1.ListVarName{VarType_1.coord_y};
    340        
    341        XAttr=Field_1.VarAttribute{VarType_1.coord_x};
    342        YAttr=Field_1.VarAttribute{VarType_1.coord_y};
    343 %        SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}];
    344        DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
    345     else
    346        XName_1=Field_1.ListVarName{VarType_1.coord(2)};
    347        YName_1=Field_1.ListVarName{VarType_1.coord(1)};
    348        if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
    349            XAttr=Field_1.VarAttribute{VarType_1.coord(2)};
    350        else
    351            XAttr=[];
    352        end
    353        if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
    354            YAttr=Field_1.VarAttribute{VarType_1.coord(1)};
    355        else
    356            YAttr=[];
    357        end     
    358     end 
    359     %check for the existence of the same  variable name
    360     XName_1_1=XName_1; %default
    361     YName_1_1=YName_1; %default
    362     UName_1_1=UName_1; %default
    363     VName_1_1=VName_1; %default
    364     for iprev=1:numel(SubData.ListVarName)
    365         switch SubData.ListVarName{iprev}
    366             case XName_1
    367                 XName_1_1=[XName_1 '_1'];
    368             case YName_1
    369                 YName_1_1=[YName_1 '_1'];
    370             case UName_1
    371                 UName_1_1=[UName_1 '_1'];
    372             case VName_1
    373                 VName_1_1=[VName_1 '_1'];
    374         end
    375     end   
     171for imod=1:numel(ivar)
     172        VarName=SubData.ListVarName{ivar(imod)};
     173        VarName_1=SubData.ListVarName{ivar_1(imod)};
     174        SubData.(VarName)=SubData.(VarName)-SubData.(VarName_1);
     175        ind_remove(ivar_1(imod))=1;
     176end
     177SubData.ListVarName(find(ind_remove))=[];
     178SubData.VarDimName(find(ind_remove))=[];
     179SubData.VarAttribute(find(ind_remove))=[];
     180       
    376181   
    377     if ~testX_1
    378           DimCell=[{XName_1_1} {YName_1_1}];
    379     end
    380     SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}];
    381     DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])];
    382     SubData.VarDimName=[SubData.VarDimName DimCell];
    383     if isfield(SubData,'VarAttribute')
    384         if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName))
    385             for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)-4
    386                 SubData.VarAttribute{ivar}=[];
    387             end
    388         end
    389         SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}];
    390     end
    391     SubData.(XName_1_1)=Field_1.(XName_1);
    392     SubData.(YName_1_1)=Field_1.(YName_1);
    393     SubData.(UName_1_1)=Field_1.(UName_1);
    394     SubData.(VName_1_1)=Field_1.(VName_1);
    395 end
    396  
  • trunk/src/uvmat.m

    r512 r515  
    7070%                        (calc_tps.m)               calculate tps coefficients (for filter projection or spatial derivatives).
    7171%                            |
    72 %                       (calc_field.m)               calculate field
    73 %                            |
    7472%                       UvData.Field-------------->histogram
    7573%               _____________|____________
    7674%              |                          |                   
    77 %        proj_field.m               proj_field.m       project the field on the projection objects              
     75%        proj_field.m               proj_field.m       project the field on the projection objects (use calc_field.m)           
    7876%              |                          |
    7977%         UvData.PlotAxes          ViewData.PlotAxes (on view_field)
    8078%              |                          |
    8179%       plot_field.m (uvmat)       plot_field.m (view_field)      plot the projected fields
    82 %
    83 % rmq: calc_field can be performed instead at the level of proj_field when needed
    8480%
    8581%
     
    209205
    210206%% TRANSFORM menu: builtin fcts
    211 transform_menu={'';'phys';'px';'phys_polar'};
     207transform_menu={'';'sub_field';'phys';'phys_polar'};
    212208UvData.OpenParam.NbBuiltin=numel(transform_menu); %number of functions
    213209path_list=(num2cell(blanks(UvData.OpenParam.NbBuiltin)))';%initialize a cell array of nbvar blanks
     
    299295
    300296%% plot input field if exists
    301 if testinputfield
    302     %delete drawn objects
    303     hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects
    304     for iobj=1:length(hother)
    305         delete_object(hother(iobj))
    306     end 
    307     if isempty(inputfile)
    308         errormsg=refresh_field(handles,[],[],[],[],[],[],{Field});
    309         set(handles.MenuTools,'Enable','on')
    310         set(handles.OBJECT_txt,'Visible','on')
    311         set(handles.edit_object,'Visible','on')
    312 %         set(handles.ListObject_1,'Visible','on')
    313         set(handles.frame_object,'Visible','on')
    314         if ~isempty(errormsg)
    315             msgbox_uvmat('ERROR',errormsg)
    316         end
    317     end
    318 end
     297% if testinputfield
     298%     %delete drawn objects
     299%     hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects
     300%     for iobj=1:length(hother)
     301%         delete_object(hother(iobj))
     302%     end 
     303%     if isempty(inputfile)
     304%         errormsg=refresh_field(handles,[],[],[],[],[],[],{Field});
     305%         set(handles.MenuTools,'Enable','on')
     306%         set(handles.OBJECT_txt,'Visible','on')
     307%         set(handles.edit_object,'Visible','on')
     308% %         set(handles.ListObject_1,'Visible','on')
     309%         set(handles.frame_object,'Visible','on')
     310%         if ~isempty(errormsg)
     311%             msgbox_uvmat('ERROR',errormsg)
     312%         end
     313%     end
     314% end
    319315
    320316set_vec_col_bar(handles) %update the display of color code for vectors
     
    679675            set(handles.j1,'String',num2stra(j1,NomType));
    680676            set(handles.j2,'String',num2stra(j2,NomType));
    681         else %read the current field index if the second file series is opened
    682             i1=str2num(get(handles.i1,'String'));
    683             i2=str2num(get(handles.i2,'String'));
    684             j1=stra2num(get(handles.j1,'String'));
    685             j2=stra2num(get(handles.j2,'String'));
     677        else %read the current field index to synchronise with the first series
     678            i1_s=str2num(get(handles.i1,'String'));
     679            i2_0=str2num(get(handles.i2,'String'));
     680            if ~isempty(i2_0)
     681                i2_s=i2_0;
     682            else
     683               i2_s=i2;
     684            end
     685            j1_0=stra2num(get(handles.j1,'String'));
     686            if ~isempty(j1_0)
     687                j1_s=j1_0;
     688            else
     689                j1_s=j1;
     690            end
     691            j2_0=stra2num(get(handles.j2,'String'));
     692            if ~isempty(j2_0)
     693                j2_s=j2_0;
     694            else
     695                j2_s=j2;
     696            end
    686697        end
    687698       
     
    705716            end
    706717            %updtate the indices of the second field series to correspond to the newly opened one
    707             FileName_1=fullfile_uvmat(Input.RootPath_1,Input.SubDir_1,Input.RootFile_1,Input.FileExt_1,Input.NomType_1,i1,i2,j1,j2);
     718            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);
    708719            if exist(FileName_1,'file')
     720                FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1_s,i2_s,j1_s,j2_s);
     721            else
    709722                FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1,i2,j1,j2);
    710                 set(handles.FileIndex_1,'String',FileIndex_1)
    711             else
    712723                msgbox_uvmat('WARNING','unable to synchronise the indices of the two series')
    713724            end
     725            set(handles.FileIndex_1,'String',FileIndex_1)
    714726        end
    715727       
     
    757769% --- Update information about a new field series (indices to scan, timing,
    758770%     calibration from an xml file, then refresh current plots
     771
    759772function update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileType,VideoObject,index)
    760773%------------------------------------------------------------------------
     
    943956        end
    944957        if ~get(handles.CheckFixLimits,'Value')
    945             set(handles.transform_fct,'Value',2); % phys transform by default if fixedLimits is off
     958            set(handles.transform_fct,'Value',3); % phys transform by default if fixedLimits is off
    946959        end
    947960        if isfield(GeometryCalib,'SliceCoord')           
     
    10141027    end
    10151028end
    1016 [ref_j,ref_i]=find(i1_series);
     1029[ref_j,ref_i]=find(squeeze(i1_series(1,:,:)));
    10171030if ~isempty(j1_series)
    10181031        state_j='on';
    10191032        if index==1
    1020             if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first vzlue
     1033            if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first value
    10211034                scan_option='j'; %scan j indext               
    10221035            end
     
    10571070
    10581071%% apply the effect of the transform fct and view the field 
     1072transform=get(handles.path_transform,'UserData');
     1073if index==2 && (~isa(transform,'function_handle')||nargin(transform)<3)
     1074    set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields
     1075end
    10591076transform_fct_Callback([],[],handles)
    1060 %run0_Callback([],[], handles); %view field
    10611077mask_test=get(handles.CheckMask,'value');
    10621078if mask_test
     
    11061122function i1_Callback(hObject, eventdata, handles)
    11071123%------------------------------------------------------------------------
    1108 set(handles.i1,'BackgroundColor',[0.7 0.7 0.7])
     1124update_ij(handles,1)
     1125
     1126%------------------------------------------------------------------------
     1127function i2_Callback(hObject, eventdata, handles)
     1128%------------------------------------------------------------------------
     1129update_ij(handles,2)
     1130
     1131%------------------------------------------------------------------------
     1132function j1_Callback(hObject, eventdata, handles)
     1133%------------------------------------------------------------------------
     1134update_ij(handles,3)
     1135
     1136%------------------------------------------------------------------------
     1137function j2_Callback(hObject, eventdata, handles)
     1138%------------------------------------------------------------------------
     1139update_ij(handles,4)
     1140
     1141%------------------------------------------------------------------------
     1142%--- update the index display after action on edit boxes i1, i2, j1 or j2
     1143function update_ij(handles,index_rank)
    11091144NomType=get(handles.NomType,'String');
    1110 i1=stra2num(get(handles.i1,'String'));
    1111 i2=stra2num(get(handles.i2,'String'));
    1112 j1=stra2num(get(handles.j1,'String'));
    1113 j2=stra2num(get(handles.j2,'String'));
    1114 indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2);
    1115 %indices=name_generator('',num1,num_a,'',NomType,1,num2,num_b,'');
     1145indices=get(handles.FileIndex,'String');
     1146[tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(indices);% the indices for the second series taken from FileIndex
     1147switch index_rank
     1148    case 1
     1149        indices=fullfile_uvmat('','','','',NomType,stra2num(get(handles.i1,'String')),i2,j1,j2);
     1150        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
     1151    case 2
     1152        indices=fullfile_uvmat('','','','',NomType,i1,stra2num(get(handles.i2,'String')),j1,j2);
     1153        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
     1154    case 3
     1155        indices=fullfile_uvmat('','','','',NomType,i1,i2,stra2num(get(handles.j1,'String')),j2);
     1156        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
     1157    case 4
     1158        indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,stra2num(get(handles.j2,'String')));
     1159        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
     1160end
    11161161set(handles.FileIndex,'String',indices)
    1117 set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7])
    1118 if get(handles.SubField,'Value')==1
     1162set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
     1163% update the second index if relevant
     1164if strcmp(get(handles.FileIndex_1,'Visible'),'on')
    11191165    NomType_1=get(handles.NomType_1,'String');
    1120     indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2);
    1121      %indices=name_generator('',num1,num_a,'',NomType_1,1,num2,num_b,'');
    1122      set(handles.FileIndex_1,'String',indices)
    1123      set(handles.FileIndex_1,'BackgroundColor',[0.7 0.7 0.7])
    1124 end
    1125 
    1126 %------------------------------------------------------------------------
    1127 function i2_Callback(hObject, eventdata, handles)
    1128 %------------------------------------------------------------------------
    1129 set(handles.i2,'BackgroundColor',[0.7 0.7 0.7])
    1130 i1_Callback(hObject, eventdata, handles)
    1131 
    1132 %------------------------------------------------------------------------
    1133 function j1_Callback(hObject, eventdata, handles)
    1134 %------------------------------------------------------------------------
    1135 set(handles.j1,'BackgroundColor',[0.7 0.7 0.7])
    1136 i1_Callback(hObject, eventdata, handles)
    1137 
    1138 %------------------------------------------------------------------------
    1139 function j2_Callback(hObject, eventdata, handles)
    1140 %------------------------------------------------------------------------
    1141 set(handles.j2,'BackgroundColor',[0.7 0.7 0.7])
    1142 i1_Callback(hObject, eventdata, handles)
     1166    indices_1=get(handles.FileIndex_1,'String');
     1167    [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
     1168    switch index_rank
     1169        case 1
     1170            indices_1=fullfile_uvmat('','','','',NomType_1,stra2num(get(handles.i1,'String')),i2_1,j1_1,j2_1);
     1171        case 2
     1172            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,stra2num(get(handles.i2,'String')),j1_1,j2_1);
     1173        case 3
     1174            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,stra2num(get(handles.j1,'String')),j2_1);
     1175        case 4
     1176            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,j1_1,stra2num(get(handles.j2,'String')));
     1177    end
     1178    set(handles.FileIndex_1,'String',indices_1)
     1179    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
     1180end
     1181   
     1182%------------------------------------------------------------------------
    11431183
    11441184%------------------------------------------------------------------------
     
    15231563if sub_value % a second input file has been entered
    15241564     [InputFile.RootPath_1,InputFile.SubDir_1,InputFile.RootFile_1,InputFile.FileIndex_1,InputFile.FileExt_1,InputFile.NomType_1]=read_file_boxes_1(handles);   
    1525     [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(InputFile.FileIndex_1);
     1565    [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
    15261566else
    15271567    filename_1=[];
     
    15291569
    15301570%% increment (or decrement) the field indices and update the input filename(s)
    1531 CheckSearch=0;
    15321571if isempty(increment)
    1533     CheckSearch=1;% search for the next available file
    15341572    set(handles.CheckFixPair,'Value',0)
    15351573end
     
    15451583            i2_1=i2_1+increment;
    15461584        end
     1585       
    15471586    else % case of scanning along index j (burst numbers)
    15481587        j1=j1+increment;
     
    16041643    if ~isempty(errormsg),return,end
    16051644    siz=size(UvData.i1_series{1});
    1606     ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1);   
     1645    ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1);
    16071646    i1_subseries=UvData.i1_series{1}(ref_indices);
    16081647    ref_indices=ref_indices(i1_subseries>0);
    16091648    if isempty(ref_indices)% case of pairs (free index i)
    1610         ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2); 
     1649        ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2);
    16111650        i1_subseries=UvData.i1_series{1}(ref_indices);
    16121651        ref_indices=ref_indices(i1_subseries>0);
     
    16231662    if ~isempty(UvData.j2_series{1})
    16241663        j2=UvData.j2_series{1}(ref_indices(end));
    1625     end 
    1626     % case of a second file series (TODO:revise)
    1627     if numel(UvData.i1_series)>=2
    1628         i1_subseries=UvData.i1_series{2}(ref_i+1,ref_j+1,:);
    1629         i1_subseries=i1_subseries(i1_subseries>0);
    1630         i1_1=i1_subseries(end);
     1664    end
     1665   
     1666    % case of a second file series
     1667    if sub_value
     1668        ref_i_1=i1_1;
     1669        if ~isempty(i2_1)
     1670            ref_i_1=floor((i1_1+i2_1)/2);% current reference index i
     1671        end
     1672        ref_j_1=1;
     1673        if ~isempty(j1_1)
     1674            ref_j_1=j1_1;
     1675            if ~isempty(j2_1)
     1676                ref_j_1=floor((j1_1+j2_1)/2);% current reference index j
     1677            end
     1678        end
     1679        if ~isempty(increment)
     1680            if get(handles.scan_i,'Value')==1% case of scanning along index i
     1681                ref_i_1=ref_i_1+increment;% increment the current reference index i
     1682            else % case of scanning along index j (burst numbers)
     1683                ref_j_1=ref_j_1+increment;% increment the current reference index j if scan_j option is used
     1684            end
     1685        else % free increment, synchronise the ref indices with the first series
     1686            ref_i_1=ref_i;
     1687            ref_j_1=ref_j;
     1688        end
     1689        if numel(UvData.i1_series)==1
     1690            UvData.i1_series{2}=UvData.i1_series{1};
     1691            UvData.j1_series{2}=UvData.j1_series{1};
     1692            UvData.i2_series{2}=UvData.i2_series{1};
     1693            UvData.j2_series{2}=UvData.j2_series{1};
     1694        end
     1695        if ref_i_1<0
     1696            errormsg='minimum i index reached';
     1697        elseif ref_j_1<0
     1698            errormsg='minimum j index reached';
     1699        elseif ref_i_1+1>size(UvData.i1_series{2},3)
     1700            errormsg='maximum i index reached for the second series (reload the input file to update the index bound)';
     1701        elseif ref_j_1+1>size(UvData.i1_series{2},2)
     1702            errormsg='maximum j index reached for the second series(reload the input file to update the index bound)';
     1703        end
     1704        if ~isempty(errormsg),return,end
     1705        siz=size(UvData.i1_series{2});
     1706        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);
     1707        i1_subseries=UvData.i1_series{2}(ref_indices);
     1708        ref_indices=ref_indices(i1_subseries>0);
     1709        if isempty(ref_indices)% case of pairs (free index i)
     1710            ref_indices=ref_i_1*siz(1)*siz(2)+1:(ref_i_1+1)*siz(1)*siz(2);
     1711            i1_subseries=UvData.i1_series{2}(ref_indices);
     1712            ref_indices=ref_indices(i1_subseries>0);
     1713        end
     1714        i1_1=UvData.i1_series{2}(ref_indices(end));
    16311715        if ~isempty(UvData.i2_series{2})
    1632             i2_subseries=UvData.i2_series{2}(ref_i+1,ref_j+1,:);
    1633             i2_subseries=i2_subseries(i2_subseries>0);
    1634             i2_1=i2_subseries(end);
     1716            i2_1=UvData.i2_series{2}(ref_indices(end));
    16351717        end
    16361718        if ~isempty(UvData.j1_series{2})
    1637             j1_subseries=UvData.j1_series{2}(ref_i+1,ref_j+1,:);
    1638             j1_subseries=j1_subseries(j1_subseries>0);
    1639             j1_1=j1_subseries(end);
     1719            j1_1=UvData.j1_series{2}(ref_indices(end));
    16401720        end
    16411721        if ~isempty(UvData.j2_series{2})
    1642             j2_subseries=UvData.j2_series{2}(ref_i+1,ref_j+1,:);
    1643             j2_subseries=j2_subseries(j2_subseries>0);
    1644             j2_1=j2_subseries(end);
    1645         end
     1722            j2_1=UvData.j2_series{1}(ref_indices(end));
     1723        end   
     1724    else% the second series (if needed) is the same file as the first
     1725        i1_1=i1;
     1726        i2_1=i2;
     1727        j1_1=j1;
     1728        j2_1=j2;
    16461729    end
    16471730end
     
    16521735
    16531736%% refresh plots
    1654 errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2);
     1737errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2,i1_1,i2_1,j1_1,j2_1);
    16551738
    16561739%% update the index counters if the index move is successfull
     
    18461929set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
    18471930drawnow
    1848 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    1849 filename=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
    1850 filename_1=[];%default
     1931[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     1932filename=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
     1933filename_1='';%default
     1934FileIndex_1='';
    18511935if get(handles.SubField,'Value')
    1852     [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
    1853     filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
     1936    [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
     1937    filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
    18541938end
    18551939num_i1=stra2num(get(handles.i1,'String'));
     
    18571941num_j1=stra2num(get(handles.j1,'String'));
    18581942num_j2=stra2num(get(handles.j2,'String'));
    1859 
    1860 errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2);
     1943[tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(FileIndex_1);
     1944
     1945errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1);
    18611946
    18621947if ~isempty(errormsg)
     
    18831968% Field: structure describing an optional input field (then replace the input file)
    18841969
    1885 function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,Field)
     1970function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1)
    18861971%------------------------------------------------------------------------
    18871972
     
    19232008if length(masknumber)>=z_index
    19242009    set(handles.masklevel,'Value',z_index)
     2010end
     2011
     2012%% test for need of tps
     2013check_proj_tps=0;
     2014if  (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))
     2015    for iobj=1:numel(UvData.Object)
     2016        if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter')
     2017            check_proj_tps=1;
     2018            break
     2019        end
     2020    end
    19252021end
    19262022
     
    19772073ParamIn.VelType=VelType;
    19782074ParamIn.GUIName='get_field';
     2075end
     2076check_tps=0;         
     2077if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamIn.FieldName,'velocity')&&~strcmp(ParamIn.FieldName,'get_field...')
     2078       check_tps=1;%tps needed to get the requested field
    19792079end
    19802080[Field{1},ParamOut,errormsg] = read_field(FileName,UvData.FileType{1},ParamIn,frame_index);
     
    20302130            ParamIn_1=UvData.MovieObject{2};
    20312131                        if ~strcmp(NomType_1,'*')
    2032                 frame_index_1=num_j1;%frame index for movies or multimage
     2132                frame_index_1=j1_1;%frame index for movies or multimage
    20332133            else
    2034                 frame_index_1=num_i1;
     2134                frame_index_1=i1_1;
    20352135            end 
    20362136         case 'multimage'
    20372137            if ~strcmp(NomType_1,'*')
    2038                 frame_index_1=num_j1;%frame index for movies or multimage
     2138                frame_index_1=j1_1;%frame index for movies or multimage
    20392139            else
    2040                 frame_index_1=num_i1;
     2140                frame_index_1=i1_1;
    20412141            end   
    20422142        case 'vol' %TODO: update
     
    20552155    end
    20562156    test_keepdata_1=0;% test for keeping the previous stored data if the input files are unchanged
    2057     if ~isequal(NomType_1,'*')%in case of a series of files (not avi movie)
    2058         if isfield(UvData,'FileName_1')%&& isfield(UvData,'VelType_1') && isfield(UvData,'FieldName_1')
    2059             test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1);
    2060         end
     2157    if ~isequal(NomType_1,'*')&& isfield(UvData,'FileName_1')
     2158           test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1);
    20612159    end
    20622160    if test_keepdata_1
     
    20682166        ParamIn_1.VelType=VelType_1;
    20692167        ParamIn_1.GUIName='get_field_1';
    2070         end
     2168        end 
    20712169        [Field{2},ParamOut_1,errormsg] = read_field(Name,UvData.FileType{2},ParamIn_1,frame_index_1);
    20722170        if ~isempty(errormsg)
    20732171            errormsg=['error in reading ' FieldName_1 ' in ' FileName_1 ': ' errormsg];
    20742172            return
     2173        end
     2174        if isstruct(ParamOut_1)&&~strcmp(ParamOut_1.FieldName,'get_field...')&& (strcmp(UvData.FileType{2},'civdata')||strcmp(UvData.FileType{2},'civx'))...
     2175                &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...')
     2176            if ~check_proj_tps
     2177             %   Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
     2178            end
    20752179        end
    20762180    end
     
    21302234        test_veltype_1=1;
    21312235        set(handles.VelType_1,'Visible','on')
    2132         menu=set_veltype_display(ParamOut_1.CivStage);
     2236        menu=set_veltype_display(ParamOut_1.CivStage,UvData.FileType{2});
    21332237        index_menu=strcmp(ParamOut_1.VelType,menu);
    21342238        set(handles.VelType_1,'Value',1+find(index_menu,1))
     
    22192323    end
    22202324    if numel(UvData.XmlData)==2
    2221         [tild,tild,tild,num_i1,num_i2,num_j1,num_j2]=fileparts_uvmat(['xx' get(handles.FileIndex_1,'String') get(handles.FileExt_1,'String')]);
    2222         if isempty(num_i2)
    2223             num_i2=num_i1;
    2224         end
    2225         if isempty(num_j1)
    2226             num_j1=1;
    2227         end
    2228         if isempty(num_j2)
    2229             num_j2=num_j1;
     2325        if isempty(i2_1)
     2326            i2_1=num_i1;
     2327        end
     2328        if isempty(j1_1)
     2329            j1_1=1;
     2330        end
     2331        if isempty(j2_1)
     2332            j2_1=j1_1;
    22302333        end
    22312334        siz=size(UvData.XmlData{2}.Time);
    2232         if ~isempty(num_i1) && siz(1)>=max(num_i1,num_i2) && siz(2)>=max(num_j1,num_j2)
    2233             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
    2234             Field{2}.Dt=(UvData.XmlData{2}.Time(num_i2,num_j2)-UvData.XmlData{2}.Time(num_i1,num_j1));
     2335        if ~isempty(i1_1) && siz(1)>=max(i1_1,i2_1) && siz(2)>=max(j1_1,j2_1)
     2336            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
     2337            Field{2}.Dt=(UvData.XmlData{2}.Time(i2_1,j2_1)-UvData.XmlData{2}.Time(i1_1,j1_1));
    22352338        end
    22362339    end
     
    22642367%% apply coordinate transform or other user fct
    22652368transform=get(handles.path_transform,'UserData');
    2266 if ~isempty(transform)
     2369if isempty(transform)
     2370    UvData.Field=Field{1};
     2371else
    22672372    XmlData=[];%default
    22682373    XmlData_1=[];%default
     
    22942399end
    22952400
    2296 %% check whether tps is needed, then calculate tps coefficients if needed
    2297 check_proj_tps=0;
    2298 if  (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))
    2299     for iobj=1:numel(UvData.Object)
    2300         if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter')
    2301             check_proj_tps=1;
    2302             break
    2303         end
    2304     end
    2305 end
    2306 check_tps=0;         
    2307 if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamOut.FieldName,'velocity')&&~strcmp(ParamOut.FieldName,'get_field...')
    2308        check_tps=1;%tps needed to get the requested field
    2309 end
    2310 if (check_tps ||check_proj_tps)&&~isfield(UvData.Field,'Coord_tps')
    2311     UvData.Field=calc_tps(UvData.Field);
    2312 end
    2313 UvData.Field.FieldList=[{ParamOut.FieldName} {ParamOut.ColorVar}];
    2314 
    2315 %% calculate scalar
    2316 if isstruct(ParamOut)&&~strcmp(ParamOut.FieldName,'get_field...')&& (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))...
    2317          &&~strcmp(ParamOut.FieldName,'velocity') && ~strcmp(ParamOut.FieldName,'get_field...')
    2318     if ~check_proj_tps
    2319         UvData.Field=calc_field([{ParamOut.FieldName} {ParamOut.ColorVar}],UvData.Field);
    2320     end
    2321 end
    2322 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'))...
    2323         &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...')
    2324     if check_proj_tps
    2325         Field{2}.FieldList=[{ParamOut_1.FieldName} {ParamOut_1.ColorVar}];
    2326     else
    2327      Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
    2328     end
    2329 end
    2330 
     2401%% calculate tps coefficients if needed
     2402UvData.Field=calc_tps(UvData.Field,check_proj_tps);
    23312403
    23322404%% analyse input field
     
    29423014    end
    29433015    set(handles.uvmat,'UserData',UvData);
    2944     run0_Callback(hObject, eventdata, handles); %run
     3016    transform_fct_list=get(handles.transform_fct,'String');
     3017    transform_fct=transform_fct_list(get(handles.transform_fct,'Value'));
     3018    if strcmp(transform_fct,'sub_field')
     3019        set(handles.transform_fct,'Value',1)%suppress the sub_field transform
     3020        transform_fct_Callback(hObject, eventdata, handles);
     3021    else
     3022        run0_Callback(hObject, eventdata, handles)
     3023    end 
    29453024else
    29463025    MenuBrowse_1_Callback(hObject, eventdata, handles)
     
    31143193FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
    31153194[tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(get(handles.FileIndex,'String'));
    3116 % set(handles.FileIndex_1,'Visible','on')
    3117 % set(handles.FileExt_1,'Visible','on')
    31183195switch field_1
    31193196    case 'get_field...'
     
    31283205        set(hhget_field.list_fig,'Value',1)
    31293206        set(hhget_field.list_fig,'String',{'uvmat'})
    3130 %         set(handles.transform_fct,'Value',1)% no transform by default
    3131 %         set(handles.path_transform,'String','')
    31323207        if check_new
    31333208            UvData.FileType{2}=UvData.FileType{1};
    31343209            set(handles.FileIndex_1,'String',get(handles.FileIndex,'String'))
    3135 %             set(handles.FileExt_1,'String',get(handles.FileExt,'String'))
    31363210              set(handles.uvmat,'UserData',UvData)
    31373211        end
     
    32393313
    32403314%------------------------------------------------------------------------
    3241 % --- Executes on button press in VelType_1.
     3315% --- Executes on choice selection in VelType_1.
    32423316function VelType_1_Callback(hObject, eventdata, handles)
    32433317%------------------------------------------------------------------------
    32443318set(handles.FixVelType,'Value',1)% the velocity type is now imposed by the GUI (not automatic)
    32453319UvData=get(handles.uvmat,'UserData');
    3246 %refresh field with a second FileName=first file name
    3247 set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
     3320set(handles.run0,'BackgroundColor',[1 1 0])%paint run0 button in yellow to indicate its activation
    32483321drawnow   
    3249 InputFile=read_GUI(handles.InputFile);
    3250 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    3251 FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
    3252 
     3322InputFile=read_GUI(handles.InputFile);% read the input file parameters
     3323[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     3324[RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
     3325FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];% name of the first input file
     3326
     3327check_refresh=0;
    32533328if isempty(InputFile.VelType_1)
    3254         FileName_1='';% we plot the current field without the second field
     3329        FileName_1='';% we plot the first input field without the second field
    32553330        set(handles.SubField,'Value',0)
    3256         SubField_Callback(hObject, eventdata, handles)
     3331        SubField_Callback(hObject, eventdata, handles)% activate SubField_Callback and refresh current display, removing the second field
    32573332elseif get(handles.SubField,'Value')% if subfield is already 'on'
    3258       [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
    3259      FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
    3260 else
    3261      FileName_1=FileName;% we compare two fields in the same file by default
     3333     FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];% name of the second input file
     3334     check_refresh=1;%will refresh the current display
     3335else% we introduce the same file (with a different field) for the second series
     3336     FileName_1=FileName;% we compare two fields in the same file
    32623337     UvData.FileType{2}=UvData.FileType{1};
     3338     UvData.XmlData{2}= UvData.XmlData{1};
    32633339     set(handles.SubField,'Value',1)
    3264 end
    3265 if isfield(UvData,'Field_1')
    3266     UvData=rmfield(UvData,'Field_1');% removes the stored second field if it exists
    3267 end
    3268 UvData.FileName_1='';% desactivate the use of a constant second file
    3269 set(handles.uvmat,'UserData',UvData)
    3270 num_i1=stra2num(get(handles.i1,'String'));
    3271 num_i2=stra2num(get(handles.i2,'String'));
    3272 num_j1=stra2num(get(handles.j1,'String'));
    3273 num_j2=stra2num(get(handles.j2,'String'));
    3274 
    3275 errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2);
    3276 
    3277 if ~isempty(errormsg)
    3278     msgbox_uvmat('ERROR',errormsg);
    3279 else
    3280     set(handles.i1,'BackgroundColor',[1 1 1])
    3281     set(handles.i2,'BackgroundColor',[1 1 1])
    3282     set(handles.j1,'BackgroundColor',[1 1 1])
    3283     set(handles.j2,'BackgroundColor',[1 1 1])
    3284     set(handles.FileIndex,'BackgroundColor',[1 1 1])
    3285     set(handles.FileIndex_1,'BackgroundColor',[1 1 1])
    3286 end
    3287 set(handles.run0,'BackgroundColor',[1 0 0])
     3340     transform=get(handles.path_transform,'UserData');
     3341     if (~isa(transform,'function_handle')||nargin(transform)<3)
     3342        set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields
     3343        transform_fct_Callback(hObject, eventdata, handles)% activate transform_fct_Callback and refresh current display
     3344     else
     3345         check_refresh=1;
     3346     end 
     3347end
     3348
     3349% refresh the current display if it has not been done previously
     3350if check_refresh
     3351    UvData.FileName_1='';% desactivate the use of a constant second file
     3352    set(handles.uvmat,'UserData',UvData)
     3353    num_i1=stra2num(get(handles.i1,'String'));
     3354    num_i2=stra2num(get(handles.i2,'String'));
     3355    num_j1=stra2num(get(handles.j1,'String'));
     3356    num_j2=stra2num(get(handles.j2,'String'));
     3357    [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(['xx' FileIndex_1]);
     3358    errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1);
     3359    if ~isempty(errormsg)
     3360        msgbox_uvmat('ERROR',errormsg);
     3361    else
     3362        set(handles.i1,'BackgroundColor',[1 1 1])
     3363        set(handles.i2,'BackgroundColor',[1 1 1])
     3364        set(handles.j1,'BackgroundColor',[1 1 1])
     3365        set(handles.j2,'BackgroundColor',[1 1 1])
     3366        set(handles.FileIndex,'BackgroundColor',[1 1 1])
     3367        set(handles.FileIndex_1,'BackgroundColor',[1 1 1])
     3368    end
     3369    set(handles.run0,'BackgroundColor',[1 0 0])
     3370end
    32883371
    32893372%-----------------------------------------------
     
    36053688
    36063689%% refresh the current plot
    3607 run0_Callback(hObject, eventdata, handles)
    3608 
    3609 
     3690if isempty(list_path{ichoice}) || nargin(transform_handle)<3
     3691    set(handles.SubField,'Value',0)
     3692    SubField_Callback(hObject, eventdata, handles)
     3693else
     3694    run0_Callback(hObject, eventdata, handles)
     3695end
    36103696
    36113697%------------------------------------------------
     
    39614047%replot all the objects within the new projected field
    39624048for IndexObj=1:numel(list_str)
    3963     IndexObj
    3964         hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat
     4049        hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat;
    39654050        if isempty(hobject) || ~ishandle(hobject)
    3966             hobject=handles.PlotAxes
     4051            hobject=handles.PlotAxes;
    39674052        end
    39684053        if isequal(IndexObj,get(handles.ListObject,'Value'))
     
    41514236        data.Type='plane';
    41524237    end
    4153     if isfield(UvData,'Field')
    4154         Field=UvData.Field;
    4155         if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh)
    4156             data.RangeX=[UvData.Field.XMin UvData.Field.XMax];
    4157             if strcmp(data.Type,'line')||strcmp(data.Type,'polyline')
    4158                 data.RangeY=UvData.Field.Mesh;
    4159             else
    4160                 data.RangeY=[UvData.Field.YMin UvData.Field.YMax];
    4161             end
    4162             data.DX=UvData.Field.Mesh;
    4163             data.DY=UvData.Field.Mesh;
    4164         end
    4165         if isfield(Field,'NbDim')&& isequal(Field.NbDim,3)
    4166             data.Coord=[0 0 0]; %default
    4167         end
    4168         if isfield(Field,'CoordUnit')
    4169             data.CoordUnit=Field.CoordUnit;
    4170         end
    4171     end
     4238%     if isfield(UvData,'Field')
     4239%         Field=UvData.Field;
     4240%         if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh)
     4241%             data.RangeX=[UvData.Field.XMin UvData.Field.XMax];
     4242%             if strcmp(data.Type,'line')||strcmp(data.Type,'polyline')
     4243%                 data.RangeY=UvData.Field.Mesh;
     4244%             else
     4245%                 data.RangeY=[UvData.Field.YMin UvData.Field.YMax];
     4246%             end
     4247%             data.DX=UvData.Field.Mesh;
     4248%             data.DY=UvData.Field.Mesh;
     4249%         end
     4250%         if isfield(Field,'NbDim')&& isequal(Field.NbDim,3)
     4251%             data.Coord=[0 0 0]; %default
     4252%         end
     4253%         if isfield(Field,'CoordUnit')
     4254%             data.CoordUnit=Field.CoordUnit;
     4255%         end
     4256%     end
    41724257    hset_object=set_object(data,[],ZBounds);
    41734258    hhset_object=guidata(hset_object);
     
    42974382UvData=get(huvmat,'UserData');
    42984383%[xx,xx,FileBase]=read_file_boxes(handles);
    4299 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
     4384[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
    43004385FileBase=fullfile(RootPath,RootFile);
    43014386 %read the current input file name
     
    45874672pos(1)=pos(1)+pos(3)-0.311+0.04; %0.311= width of the geometry_calib interface (units relative to the srcreen)
    45884673pos(2)=pos(2)-0.02;
    4589 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    4590 FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
     4674[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     4675FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
    45914676set(handles.view_xml,'Backgroundcolor',[1 1 0])%indicate the reading of the current xml file by geometry_calib
    45924677if isfield(UvData.OpenParam,'CalOrigin')
     
    47204805
    47214806%prepare display of the set_grid GUI
    4722 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    4723 FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
     4807[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     4808FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
    47244809CoordList=get(handles.transform_fct,'String');
    47254810val=get(handles.transform_fct,'Value');
     
    47464831%------------------------------------------------------------------------
    47474832series; %first display of the GUI to fill waiting time
    4748 [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    4749 param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
     4833[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     4834param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
    47504835if isequal(get(handles.SubField,'Value'),1)
    4751     [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
    4752     FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
     4836    [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
     4837    FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
    47534838    if ~isequal(FileName_1,param.FileName)
    47544839        param.FileName_1=FileName_1;
     
    47884873function MenuPIV_Callback(hObject, eventdata, handles)
    47894874%------------------------------------------------------------------------
    4790  [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
    4791  FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
     4875 [RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     4876 FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
    47924877civ(FileName);% interface de civ(not in the uvmat file)
    47934878
Note: See TracChangeset for help on using the changeset viewer.