Changeset 863


Ignore:
Timestamp:
Feb 3, 2015, 7:59:58 PM (10 years ago)
Author:
sommeria
Message:

civ2vel3C_introduced

Location:
trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/calc_field_interp.m

    r809 r863  
    5050    Operator{ilist}='';%default empty operator (vec, norm,...)
    5151    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
    52     if isempty(r) % the field name is a variable itself
     52    if isempty(r) % no operator: the field name is a variable itself
    5353        ivar=find(strcmp(FieldName{ilist},Data.ListVarName));
    54         if isempty(ivar)
     54        if isempty(ivar)% the requested variable does not exist
    5555            check_skipped(ilist)=1; %variable not found
    56         elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
     56        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));% the variable exists and has not been already selected
    5757            if isfield(Data.VarAttribute{ivar},'Role') &&...
    58                 (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
    59                 check_interp(ilist)=0; % ancillary variable, not interpolated     
     58                    (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
     59                check_interp(ilist)=0; % ancillary variable, not interpolated ?????
     60                check_skipped(ilist)=1; %variable not used
     61            else
     62                InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables
    6063            end
    61             InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
    6264        end
    6365    else
    … …  
    161163
    162164%% put an error flag to indicate NaN data
    163 if exist('XI','var')&&~isempty(VarVal)
    164     nbvar=numel(VarVal);
    165     ListVarName{nbvar+1}='FF';
    166     VarVal{nbvar+1}=isnan(VarVal{nbvar});
    167     VarAttribute{nbvar+1}.Role='errorflag';
    168 end
     165% if exist('XI','var')&&~isempty(VarVal)
     166%     nbvar=numel(VarVal);
     167%     ListVarName{nbvar+1}='FF';
     168%     VarVal{nbvar+1}=isnan(VarVal{nbvar});
     169%     VarAttribute{nbvar+1}.Role='errorflag';
     170% end
    169171
    170172

  • trunk/src/check_files.m

    r856 r863  
    9595    'rotate_points';...%'rotate_points': associated with GUI rotate_points.fig to introduce (2D) rotation parameters
    9696    'rotate_points.fig';...
    97     'RUN_STLIN';...% combine 2 displacement fields for stereo PIV
    9897    'series';...% master function for analysis field series, with interface 'series.fig'
    9998    'series.fig';...% interface for 'series'

  • trunk/src/get_file_series.m

    r809 r863  
    6060            Param.IndexRange.PairString={Param.IndexRange.PairString};
    6161        end
     62        if ~isempty(Param.IndexRange.PairString{iview,1})
    6263        r=regexp(Param.IndexRange.PairString{iview,1},'(?<mode>(Di=)|(Dj=)) -*(?<num1>\d+)\|(?<num2>\d+)','names');%look for mode=Dj or Di
    6364        if isempty(r)
    6465            r=regexp(Param.IndexRange.PairString{iview,1},'(?<num1>\d+)(?<mode>-)(?<num2>\d+)','names');%look for burst pairs
     66        end
    6567        end
    6668        % TODO case of free pairs:

  • trunk/src/proj_field.m

    r854 r863  
    485485end
    486486
    487 
    488487%-----------------------------------------------------------------
    489488%project on a line
    … …  
    696695                end
    697696            elseif isequal(ProjMode,'interp_lin')  %filtering %linear interpolation:
     697                if ~isequal(errorflag,0)
     698                    VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
     699                    indsel=find(FieldData.(VarName_FF)==0);
     700                    coord_x=coord_x(indsel);
     701                    coord_y=coord_y(indsel);
     702                    for ivar=1:numel(CellInfo{icell}.VarIndex)
     703                        VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};
     704                        FieldData.(VarName)=FieldData.(VarName)(indsel);
     705                    end
     706                end                   
    698707                [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI);
    699708                ProjData.X=Xproj;

  • trunk/src/series.m

    r857 r863  
    22652265            end
    22662266        end
    2267     case 'first'
     2267    case 'one'
    22682268        OutputDirVisible='on';
    22692269        SubDirOut=InputTable{1,2}; %use the first subdir name (+OutputDirExt) as output  subdirectory
     2270    case 'two'
     2271        OutputDirVisible='on';   
     2272        SubDir=InputTable(1:2,2); %set of subdirectories
     2273        SubDirOut=SubDir{1};
     2274        if numel(SubDir)>1
     2275                SubDirOut=[SubDirOut '-' regexprep(SubDir{2},'^/','')];
     2276        end
    22702277    case 'last'
    22712278        OutputDirVisible='on';

  • trunk/src/series/civ2vel_3C.m

    r862 r863  
    5555%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed
    5656if isstruct(Param) && isequal(Param.Action.RUN,0)
    57     ParamOut.AllowInputSort='on';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
     57    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    5858    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
    5959    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
    … …  
    6262    ParamOut.FieldTransform = 'off';%use the phys  transform function without choice
    6363    %ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only)
    64     ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
     64    ParamOut.ProjObject='on';%can use projection object(option 'off'/'on',
    6565    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
    6666    ParamOut.OutputDirExt='.vel3C';%set the output dir extension
     67    ParamOut.OutputSubDirMode='two'; % the two first input lines are used to define the output subfolder
    6768    ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
    68       %check the input files
     69    %check the input files
    6970    first_j=[];
     71    if size(Param.InputTable,1)<2
     72        msgbox_uvmat('WARNING',['two or three input file series are needed'])
     73    end
    7074    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
    7175    PairString='';
    … …  
    9498WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    9599
     100
     101%% root input file(s) name, type and index series
     102RootPath=Param.InputTable(:,1);
     103RootFile=Param.InputTable(:,3);
     104SubDir=Param.InputTable(:,2);
     105NomType=Param.InputTable(:,4);
     106FileExt=Param.InputTable(:,5);
     107hdisp=disp_uvmat('WAITING...','checking the file series',checkrun);
     108[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     109if ~isempty(hdisp),delete(hdisp),end;
     110%%%%%%%%%%%%
     111% The cell array filecell is the list of input file names, while
     112% filecell{iview,fileindex}:
     113%        iview: line in the table corresponding to a given file series
     114%        fileindex: file index within  the file series,
     115% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
     116% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
     117%%%%%%%%%%%%
     118NbView=numel(i1_series);%number of input file series (lines in InputTable)
     119NbField_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
     120NbField_i=size(i1_series{1},2); %nb of fields for the i index
     121NbField=NbField_j*NbField_i; %total number of fields
     122
    96123%% define the directory for result file (with path=RootPath{1})
    97124OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    98 % 
     125%
    99126% if ~isfield(Param,'InputFields')
    100127%     Param.InputFields.FieldName='';
    … …  
    105132if size(time,1)>1
    106133    diff_time=max(max(diff(time)));
    107     if diff_time>0 
     134    if diff_time>0
    108135        disp_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time) ': the mean time is chosen in result'],checkrun)
    109     end   
     136    end
    110137end
    111138if ~isempty(errormsg)
    112139    disp_uvmat('WARNING',errormsg,checkrun)
    113140end
    114 time=mean(time,1); %averaged time taken for the merged field         
     141time=mean(time,1); %averaged time taken for the merged field
    115142if isfield(XmlData{1},'GeometryCalib')
    116      tsaiA=XmlData{1}.GeometryCalib;
    117  else
    118      disp_uvmat('ERROR','no geometric calibration available for image A',checkrun)
    119      return
    120  end
    121  if isfield(XmlData{2},'GeometryCalib')
    122      tsaiB=XmlData{2}.GeometryCalib;
    123  else
    124      disp_uvmat('ERROR','no geometric calibration available for image B',checkrun)
    125      return
    126  end
     143    tsaiA=XmlData{1}.GeometryCalib;
     144else
     145    disp_uvmat('ERROR','no geometric calibration available for image A',checkrun)
     146    return
     147end
     148if isfield(XmlData{2},'GeometryCalib')
     149    tsaiB=XmlData{2}.GeometryCalib;
     150else
     151    disp_uvmat('ERROR','no geometric calibration available for image B',checkrun)
     152    return
     153end
    127154[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
    128155
    129156%% grid of physical positions (given by projection plane)
    130157if ~Param.CheckObject
    131          disp_uvmat('ERROR','a projection plane with interpolation is needed',checkrun)
    132      return
     158    disp_uvmat('ERROR','a projection plane with interpolation is needed',checkrun)
     159    return
    133160end
    134161ObjectData=Param.ProjObject;
    135 
    136 [x,y]=meshgrid(ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2),ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2));
    137 z=zeros(size(x));
    138 [Field,ParamOut,errormsg] = read_field(FileName,FileType,ParamIn,num)
    139 %camera coordinates:initialisation 2 cameras
    140 NbCamera=2;
    141 X=zeros(NbCamera,size(x,1),size(y,1));
    142 Y=zeros(NbCamera,size(x,1),size(y,1));
    143 for icamera=1:NbCamera
    144 %camera coordinates
    145 Calib=XmlData{1}.GeometryCalib;
    146     xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx_Ty_Tz(1);
    147     yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Tx_Ty_Tz(2);
    148     zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tx_Ty_Tz(3);
    149    
    150     %undistorted image coordinates
    151     X(icamera,:,:)=xc./zc;
    152     Y(icamera,:,:)=yc./zc;
    153 end
    154 
    155 %% case of fixed planes (for moving interface, coeff need to be calculated for each field
    156 if isequal(size(InputTable),2)
    157 [S,D]=get_coeff(XmlData,X,Y);
    158 end
    159 
    160  %% MAIN LOOP ON FIELDS
     162xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
     163yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
     164[XI,YI]=meshgrid(xI,yI);
     165U=zeros(size(XI,1),size(XI,2));
     166V=zeros(size(XI,1),size(XI,2));
     167W=zeros(size(XI,1),size(XI,2));
     168
     169%% MAIN LOOP ON FIELDS
     170warning off
    161171for index=1:NbField
    162         update_waitbar(WaitbarHandle,index/NbField)
     172    update_waitbar(WaitbarHandle,index/NbField)
    163173    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    164174        disp('program stopped by user')
    … …  
    169179    Data=cell(1,NbView);%initiate the set Data
    170180    timeread=zeros(1,NbView);
    171     for iview=1:NbView
     181    [Data{3},tild,errormsg] = nc2struct(filecell{3,index});
     182    ZI=griddata(Data{3}.Xphys,Data{3}.Yphys,Data{3}.Zphys,XI,YI);
     183    [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a
     184    [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b
     185   
     186    %trouver z
     187    % trouver les coordonnées px sur chaque image.
     188    %A=
     189    for iview=1:2
    172190        %% reading input file(s)
    173         [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));
     191        [Data{iview},tild,errormsg]=read_civdata(filecell{iview,index},{'vec(U,V)'},'*');
    174192        if ~isempty(errormsg)
    175193            disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)
    … …  
    177195        end
    178196        % get the time defined in the current file if not already defined from the xml file
    179         if ~isempty(time) && isfield(Data{iview},'Time')
    180             timeread(iview)=Data{iview}.Time;
    181         end
    182         if ~isempty(NbSlice_calib)
    183             Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
    184         end
    185        
    186         %% transform the input field (e.g; phys) if requested (no transform involving two input fields)
    187             %camera coordinates
    188         Data{iview}=phys(Data{iview},XmlData{iview});
    189        
    190         %% projection on object (gridded plane)
    191 %         if Param.CheckObject
    192             [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);
    193             if ~isempty(errormsg)
    194                 disp_uvmat('ERROR',['ERROR in merge_proge/proj_field: ' errormsg],checkrun)
    195                 return
    196             end
    197 
    198     end
    199     %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%%
    200 
    201     %% merge the NbView fields
    202     [MergeData,errormsg]=merge_field(Data);
    203     if ~isempty(errormsg)
    204         disp_uvmat('ERROR',errormsg,checkrun);
    205         return
    206     end
    207 
    208     %% time of the merged field: take the average of the different views
    209     if ~isempty(time)
    210         timeread=time(index);   
    211     elseif ~isempty(find(timeread))% time defined from ImaDoc
    212         timeread=mean(timeread(timeread~=0));% take average over times form the files (when defined)
    213     else
    214         timeread=index;% take time=file index
    215     end
    216 
     197        if isfield(Data{iview},'Time')&& isequal(Data{iview}.Time,Data{1}.Time)
     198            Time=Data{iview}.Time;
     199        else
     200            disp_uvmat('ERROR','Time undefined or not synchronous',checkrun)
     201            return
     202        end
     203        if isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Data{1}.Dt)
     204            Dt=Data{iview}.Dt;
     205        else
     206            disp_uvmat('ERROR','Dt undefined or not synchronous',checkrun)
     207            return
     208        end
     209    end
     210    Ua=griddata(Data{1}.X,Data{1}.Y,Data{1}.U,Xa,Ya);
     211    Va=griddata(Data{1}.X,Data{1}.Y,Data{1}.V,Xa,Ya);
     212    A=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,YI,YI,ZI);
     213   
     214    Ub=griddata(Data{2}.X,Data{2}.Y,Data{1}.U,Xb,Yb);
     215    Vb=griddata(Data{2}.X,Data{2}.Y,Data{2}.V,Xb,Yb);
     216    B=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,YI,YI,ZI);
     217    S=ones(size(XI,1),size(XI,2),3);
     218    D=ones(size(XI,1),size(XI,2),3,3);
     219    S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
     220    S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
     221    S(:,:,3)=A(:,:,1,3).*Ua+A(:,:,2,3).*Va+B(:,:,1,3).*Ub+B(:,:,2,3).*Vb;
     222    D(:,:,1,1)=A(:,:,1,1).*A(:,:,1,1)+A(:,:,2,1).*A(:,:,2,1)+B(:,:,1,1).*B(:,:,1,1)+B(:,:,2,1).*B(:,:,2,1);
     223    D(:,:,1,2)=A(:,:,1,1).*A(:,:,1,2)+A(:,:,2,1).*A(:,:,2,2)+B(:,:,1,1).*B(:,:,1,2)+B(:,:,2,1).*B(:,:,2,2);
     224    D(:,:,1,3)=A(:,:,1,1).*A(:,:,1,3)+A(:,:,2,1).*A(:,:,2,3)+B(:,:,1,1).*B(:,:,1,3)+B(:,:,2,1).*B(:,:,2,3);
     225    D(:,:,2,1)=A(:,:,1,2).*A(:,:,1,1)+A(:,:,2,2).*A(:,:,2,1)+B(:,:,1,2).*B(:,:,1,1)+B(:,:,2,2).*B(:,:,2,1);
     226    D(:,:,2,2)=A(:,:,1,2).*A(:,:,1,2)+A(:,:,2,2).*A(:,:,2,2)+B(:,:,1,2).*B(:,:,1,2)+B(:,:,2,2).*B(:,:,2,2);
     227    D(:,:,2,3)=A(:,:,1,2).*A(:,:,1,3)+A(:,:,2,2).*A(:,:,2,3)+B(:,:,1,2).*B(:,:,1,3)+B(:,:,2,2).*B(:,:,2,3);
     228    D(:,:,3,1)=A(:,:,1,3).*A(:,:,1,1)+A(:,:,2,3).*A(:,:,2,1)+B(:,:,1,3).*B(:,:,1,1)+B(:,:,2,3).*B(:,:,2,1);
     229    D(:,:,3,2)=A(:,:,1,3).*A(:,:,1,2)+A(:,:,2,3).*A(:,:,2,2)+B(:,:,1,3).*B(:,:,1,2)+B(:,:,2,3).*B(:,:,2,2);
     230    D(:,:,3,3)=A(:,:,1,3).*A(:,:,1,3)+A(:,:,2,3).*A(:,:,2,3)+B(:,:,1,3).*B(:,:,1,3)+B(:,:,2,3).*B(:,:,2,3);
     231    for indj=1:size(XI,1)
     232        for indi=1:size(XI,2)
     233            dxyz=squeeze(S(indj,indi,:))\squeeze(D(indj,indi,:,:));
     234            U(indj,indi)=dxyz(1);
     235            V(indj,indi)=dxyz(2);
     236            W(indj,indi)=dxyz(3);
     237        end
     238    end   
     239    Error=zeros(size(XI,1),size(XI,2),4);
     240    Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua;
     241    Error(:,:,2)=A(:,:,2,1).*U+A(:,:,2,2).*V+A(:,:,2,3).*W-Va;
     242    Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*W-Ub;
     243    Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb;
     244   
    217245    %% generating the name of the merged field
    218246    i1=i1_series{1}(index);
    … …  
    232260        end
    233261    end
    234     OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomType{1},i1,i2,j1,j2);
    235 
     262    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2);
     263   
    236264    %% recording the merged field
    237    
    238         MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim'};
     265    if index==1% initiate the structure at first index
     266        MergeData.ListGlobalAttribute={'Conventions','Time','Dt'};
    239267        MergeData.Conventions='uvmat';
    240         MergeData.nb_coord=2;
    241         MergeData.nb_dim=2;
    242         dt=[];
    243         if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)
    244             dt=Data{1}.dt;
    245         end
    246         for iview =2:numel(Data)
    247             if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))
    248                 dt=[];%dt not the same for all fields
    249             end
    250         end
    251         if ~isempty(timeread)
    252             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
    253             MergeData.Time=timeread;
    254         end
    255         if ~isempty(dt)
    256             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'dt'}];
    257             MergeData.dt=dt;
    258         end
    259         error=struct2nc(OutputFile,MergeData);%save result file
    260         if isempty(error)
    261             disp(['output file ' OutputFile ' written'])
    262         else
    263             disp(error)
    264         end
    265 end
    266 
    267 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    268 % %read the velocity fields
    269 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    270 %read field A
    271 [Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);
    272 %removes false vectors
    273 if isfield(Field,'FF')
    274     Field.X=Field.X(find(Field.FF==0));
    275     Field.Y=Field.Y(find(Field.FF==0));
    276     Field.U=Field.U(find(Field.FF==0));
    277     Field.V=Field.V(find(Field.FF==0));
    278 end
    279 %interpolate on the grid common to both images in phys coordinates
    280 dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);
    281 dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);
    282 dt=Field.dt;
    283 time=Field.Time;
    284 
    285 %read field B
    286 [Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);
    287 if ~isequal(Field.dt,dt)
    288     msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')
    289      return
    290 end
    291 if ~isequal(Field.Time,time)
    292     msgbox_uvmat('ERROR','different times for the two velocity fields ')
    293      return
    294 end
    295 %removes false vectors
    296 if isfield(Field,'FF')
    297 Field.X=Field.X(find(Field.FF==0));
    298 Field.Y=Field.Y(find(Field.FF==0));
    299 Field.U=Field.U(find(Field.FF==0));
    300 Field.V=Field.V(find(Field.FF==0));
    301 end
    302 %interpolate on XimaB
    303 dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);
    304 dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);
    305 %eliminate Not-a-Number
    306 ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
    307 dXa=dXa(ind_Nan);
    308 dYa=dYa(ind_Nan);
    309 dXb=dXb(ind_Nan);
    310 dYb=dYb(ind_Nan);
    311 grid_phys1(:,1)=grid_real_x(ind_Nan);
    312 grid_phys1(:,2)=grid_real_y(ind_Nan);
    313  
    314 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    315 %compute the differential coefficients of the geometric calibration
    316 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    317 [A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);
    318 [B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);
    319 
    320 C1=A11.*A22-A12.*A21;
    321 C2=A13.*A22-A12.*A23;
    322 C3=A13.*A21-A11.*A23;
    323 D1=B11.*B22-B12.*B21;
    324 D2=B13.*B22-B12.*B23;
    325 D3=B13.*B21-B11.*B23;
    326 A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
    327 A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
    328 B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
    329 B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
    330 Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
    331 
    332 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    333 %Projection for compatible displacements
    334 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    335 Ua=dXa-Lambda.*A1;
    336 Va=dYa-Lambda.*A2;
    337 Ub=dXb-Lambda.*B1;
    338 Vb=dYb-Lambda.*B2;
    339 
    340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    341 %Calculations of displacements and error
    342 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    343 U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
    344 V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
    345 W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
    346 W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
    347 
    348 error=sqrt((A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb).*(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2));
    349 
    350 ind_error=(find(error<thresh_patch));
    351 U=U(ind_error);
    352 V=V(ind_error);
    353 W=W(ind_error);%correction for water interface
    354 error=error(ind_error);
    355 
    356 %create nc grid file
    357 Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};
    358 Result.nb_coord=3;%grid file, no velocity
    359 Result.nb_dim=2;
    360 Result.constant_pixcm=0;%no linear correspondance with images
    361 Result.absolut_time_T0=time;%absolute time of the field
    362 Result.hart=0;
    363 Result.dt=dt;%time interval for image correlation (put  by default)
    364 % cte.title='grid';
    365 Result.civ=0;%not a civ file (no direct correspondance with an image)
    366 % Result.ListDimName={'nb_vectors'}
    367 % Result.DimValue=length(U);
    368 Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};
    369 Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}
    370 Result.vec_X= grid_phys1(ind_error,1);
    371 Result.vec_Y= grid_phys1(ind_error,2);
    372 Result.vec_U=U/dt;
    373 Result.vec_V=V/dt;
    374 Result.vec_W=W/dt;
    375 Result.vec_E=error;
    376 % error=write_netcdf(file_st,cte,fieldlabels,grid_phys);
    377 error=struct2nc(file_st,Result);
    378 display([file_st ' written'])
    379 
    380 function [S,D]=get_coeff(XmlData,X,Y);
    381 Calib_a=XmlData{1}.GeometryCalib;
    382 R=(Calib_a.R)';%rotation matrix
    383 A(1,1)=R(1)-R(7)*X;
    384 A(1,2)=R(2)-R(8)*X;
    385 A(1,3)=R(3)-R(9)*X;
    386 A(2,1)=R(4)-R(7)*Y;
    387 A(2,2)=R(5)-R(8)*Y;
    388 A(2,3)=R(6)-R(9)*Y;
    389 
    390 %'pxcm_tsai': find differentials of the Tsai calibration
    391 function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)
    392 R=(a.R)';
    393 
    394 x=var_phys(:,1);
    395 y=var_phys(:,2);
    396 
    397 if isfield(a,'PlanePos')
    398     prompt={'Plane 1 Index','Plane 2 Index'};
    399     Rep=inputdlg(prompt,'Target displacement test');
    400     Z1=str2double(Rep(1));
    401     Z2=str2double(Rep(2));
    402     z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2
    403 else
    404     z=0;
    405 end
    406 
    407 %transform coeff for differentiels
    408 a.C11=R(1)*R(8)-R(2)*R(7);
    409 a.C12=R(2)*R(7)-R(1)*R(8);
    410 a.C21=R(4)*R(8)-R(5)*R(7);
    411 a.C22=R(5)*R(7)-R(4)*R(8);
    412 a.C1x=R(3)*R(7)-R(9)*R(1);
    413 a.C1y=R(3)*R(8)-R(9)*R(2);
    414 a.C2x=R(6)*R(7)-R(9)*R(4);
    415 a.C2y=R(6)*R(8)-R(9)*R(5);
    416 
    417 % %dependence in x,y
    418 % denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);
    419 % A11=(a.f*a.sx*(a.C11*y-a.C1x*z+R(1)*a.Tz-R(7)*a.Tx)./denom)/a.dpx;
    420 % A12=(a.f*a.sx*(a.C12*x-a.C1y*z+R(2)*a.Tz-R(8)*a.Tx)./denom)/a.dpx;
    421 % A21=(a.f*a.sx*(a.C21*y-a.C2x*z+R(4)*a.Tz-R(7)*a.Ty)./denom)/a.dpy;
    422 % A22=(a.f*(a.C22*x-a.C2y*z+R(5)*a.Tz-R(8)*a.Ty)./denom)/a.dpy;
    423 % A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.Tz-R(9)*a.Tx)./denom)/a.dpx;
    424 % A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.Tz-R(9)*a.Ty)./denom)/a.dpy;
    425 
    426 %dependence in x,y
    427 denom=(R(7)*x+R(8)*y+R(9)*z+a.Tx_Ty_Tz(3)).*(R(7)*x+R(8)*y+R(9)*z+a.Tx_Ty_Tz(3));
    428 A11=(a.fx_fy(1)*(a.C11*y-a.C1x*z+R(1)*a.Tx_Ty_Tz(3)-R(7)*a.Tx_Ty_Tz(1))./denom);
    429 A12=(a.fx_fy(1)*(a.C12*x-a.C1y*z+R(2)*a.Tx_Ty_Tz(3)-R(8)*a.Tx_Ty_Tz(1))./denom);
    430 A21=(a.fx_fy(1)*(a.C21*y-a.C2x*z+R(4)*a.Tx_Ty_Tz(3)-R(7)*a.Tx_Ty_Tz(2))./denom);
    431 A22=(a.fx_fy(2)*(a.C22*x-a.C2y*z+R(5)*a.Tx_Ty_Tz(3)-R(8)*a.Tx_Ty_Tz(2))./denom);
    432 A13=(a.fx_fy(2)*(a.C1x*x+a.C1y*y+R(3)*a.Tx_Ty_Tz(3)-R(9)*a.Tx_Ty_Tz(1))./denom);
    433 A23=(a.fx_fy(2)*(a.C2x*x+a.C2y*y+R(6)*a.Tx_Ty_Tz(3)-R(9)*a.Tx_Ty_Tz(2))./denom);
    434 
    435 
    436 
    437 
    438 
    439 
    440 
    441 
     268        MergeData.Time=Time;
     269        MergeData.Dt=Dt;
     270        MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'};
     271        MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
     272                {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     273        MergeData.coord_x=xI;
     274        MergeData.coord_y=yI;
     275        MergeData.Z=ZI;
     276    end
     277    MergeData.U=U/Dt;
     278    MergeData.V=V/Dt;
     279    MergeData.W=W/Dt;
     280    MergeData.Error=sqrt(sum(Error.*Error,3));
     281    errormsg=struct2nc(OutputFile,MergeData);%save result file
     282    if isempty(errormsg)
     283        disp(['output file ' OutputFile ' written'])
     284    else
     285        disp(errormsg)
     286    end
     287end
     288
     289
     290function A=get_coeff(Calib,X,Y,x,y,z)
     291R=(Calib.R)';%rotation matrix
     292T_z=Calib.Tx_Ty_Tz(3);
     293T=R(7)*x+R(8)*y+R(9)*z+T_z;
     294A(:,:,1,1)=(R(1)-R(7)*X)./T;
     295A(:,:,1,2)=(R(2)-R(8)*X)./T;
     296A(:,:,1,3)=(R(3)-R(9)*X)./T;
     297A(:,:,2,1)=(R(4)-R(7)*Y)./T;
     298A(:,:,2,2)=(R(5)-R(8)*Y)./T;
     299A(:,:,2,3)=(R(6)-R(9)*Y)./T;
     300
     301
     302
     303
     304
     305
     306
Note: See TracChangeset for help on using the changeset viewer.