Ignore:
Timestamp:
May 13, 2024, 9:49:09 PM (6 months ago)
Author:
sommeria
Message:

false flags modified in civ_series with new conventions, low pass filter adde to transform_field

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/civ2vel_3C.m

    r1141 r1144  
    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)
    60     ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     60    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    6161    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    6262    ParamOut.FieldTransform = 'off';%use the phys  transform function without choice
     
    108108hdisp=disp_uvmat('WAITING...','checking the file series',checkrun);
    109109[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
    110 if ~isempty(hdisp),delete(hdisp),end;
     110if ~isempty(hdisp),delete(hdisp),end
    111111%%%%%%%%%%%%
    112112% The cell array filecell is the list of input file names, while
     
    124124%% define the directory for result file (with path=RootPath{1})
    125125OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    126 %
    127 % if ~isfield(Param,'InputFields')
    128 %     Param.InputFields.FieldName='';
    129 % end
    130126
    131127%% calibration data and timing: read the ImaDoc files
     
    175171    CheckOverwrite=Param.CheckOverwrite;
    176172end
     173
    177174for index=1:NbField
    178175   
    179 
    180    
    181       %% generating the name of the merged field
     176    %% generating the name of the merged field
    182177    i1=i1_series{1}(index);
    183178    if ~isempty(i2_series{end})
     
    198193    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2);
    199194   
    200     %%
    201    
    202    
    203195    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    204196        disp('program stopped by user')
     
    206198    end
    207199   
    208      if (~CheckOverwrite && exist(OutputFile,'file')) 
    209             disp('existing output file already exists, skip to next field')
    210             continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    211      end   
    212      
     200    if (~CheckOverwrite && exist(OutputFile,'file'))
     201        disp('existing output file already exists, skip to next field')
     202        continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     203    end
     204   
    213205    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    214206    Data=cell(1,NbView);%initiate the set Data
     
    217209    %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken
    218210    %at the middle between to time step
    219    clear ZItemp
    220    ZItemp=zeros(size(XI,1),size(XI,2),2);
    221    
    222    if index==1
     211    clear ZItemp
     212    ZItemp=zeros(size(XI,1),size(XI,2),2);
     213    CheckZ=0;
     214   
     215    if index==1
    223216        first_img=i1_series{1,1}(1,1); %id of the first image of the series
    224    end
    225      
    226      idtemp=0;
    227  for indextemp=index:index+1;
    228      idtemp=idtemp+1;
    229     if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys
    230      
     217    end
     218   
     219    idtemp=0;
     220    for indextemp=index:index+1
     221        idtemp=idtemp+1;
     222        if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys
     223           
     224            [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
     225           
     226            if  exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector
     227                temp=find(Data{3}.Civ3_FF==0);
     228                Zphys=Data{3}.Zphys(temp);
     229                Yphys=Data{3}.Yphys(temp);
     230                Xphys=Data{3}.Xphys(temp);
     231            else
     232                Zphys=Data{3}.Zphys;
     233                Yphys=Data{3}.Yphys;
     234                Xphys=Data{3}.Xphys;
     235            end
     236           
     237        elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys
     238           
     239            %test if the seconde camera is the same for both folder
     240            for i=3:4
     241                indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1
     242                indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1
     243                camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name
     244            end
     245           
     246            if strcmp(camname{3},camname{4})==0
     247                disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)
     248                return
     249            end
     250           
     251            [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
     252           
     253            if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector
     254                temp=find(Data{3}.Civ3_FF==0);
     255                Xmid3=Data{3}.Xmid(temp);
     256                Ymid3=Data{3}.Ymid(temp);
     257                U3=Data{3}.Uphys(temp);
     258                V3=Data{3}.Vphys(temp);
     259            else
     260                Xmid3=Data{3}.Xmid;
     261                Ymid3=Data{3}.Ymid;
     262                U3=Data{3}.Uphys;
     263                V3=Data{3}.Vphys;
     264            end
     265            %temporary gridd of merging the 2 stereos datas
     266            [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)-min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)-min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2));
     267           
     268            %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera
     269            %(Dalsa 3)
     270            x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq);
     271            y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq);
     272           
     273            [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp-1),'.nc']);
     274            if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector
     275                temp=find(Data{4}.Civ3_FF==0);
     276                Xmid4=Data{4}.Xmid(temp);
     277                Ymid4=Data{4}.Ymid(temp);
     278                U4=Data{4}.Uphys(temp);
     279                V4=Data{4}.Vphys(temp);
     280            else
     281                Xmid4=Data{4}.Xmid;
     282                Ymid4=Data{4}.Ymid;
     283                U4=Data{4}.Uphys;
     284                V4=Data{4}.Vphys;
     285            end
     286           
     287            %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera
     288            %(Dalsa 3)
     289            x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq);
     290            y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq);
     291           
     292            xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1);
     293            ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1);
     294            u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1);
     295            v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1);
     296           
     297           
     298            [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys
     299            %remove NaN
     300            tempNaN=isnan(Zphys);tempind=find(tempNaN==1);
     301            Zphys(tempind)=[];
     302            Xphys(tempind)=[];
     303            Yphys(tempind)=[];       
     304        end
    231305       
    232        
    233         [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
    234        
    235         if  exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector
    236             temp=find(Data{3}.Civ3_FF==0);
    237             Zphys=Data{3}.Zphys(temp);
    238             Yphys=Data{3}.Yphys(temp);
    239             Xphys=Data{3}.Xphys(temp);
    240         else
    241             Zphys=Data{3}.Zphys;
    242             Yphys=Data{3}.Yphys;
    243             Xphys=Data{3}.Xphys;
    244         end
    245        
    246        
    247        
    248     elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys
    249        
    250        
    251         %test if the seconde camera is the same for both folder
    252         for i=3:4
    253         indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1
    254         indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1
    255         camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name
    256         end
    257        
    258         if strcmp(camname{3},camname{4})==0
    259             disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)
    260             return
    261         end
    262        
    263    
    264        
    265         [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
    266    
    267         if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector
    268             temp=find(Data{3}.Civ3_FF==0);
    269             Xmid3=Data{3}.Xmid(temp);
    270             Ymid3=Data{3}.Ymid(temp);
    271             U3=Data{3}.Uphys(temp);
    272             V3=Data{3}.Vphys(temp);
    273         else
    274             Xmid3=Data{3}.Xmid;
    275             Ymid3=Data{3}.Ymid;
    276             U3=Data{3}.Uphys;
    277             V3=Data{3}.Vphys;
    278         end
    279         %temporary gridd of merging the 2 stereos datas
    280         [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)-min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)-min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2));
    281        
    282         %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera
    283         %(Dalsa 3)
    284         x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq);
    285         y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq);
    286        
    287        
    288 
    289          [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp-1),'.nc']);
    290         if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector
    291             temp=find(Data{4}.Civ3_FF==0);
    292             Xmid4=Data{4}.Xmid(temp);
    293             Ymid4=Data{4}.Ymid(temp);
    294             U4=Data{4}.Uphys(temp);
    295             V4=Data{4}.Vphys(temp);
    296         else
    297             Xmid4=Data{4}.Xmid;
    298             Ymid4=Data{4}.Ymid;
    299             U4=Data{4}.Uphys;
    300             V4=Data{4}.Vphys;
    301         end
    302        
    303         %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera
    304         %(Dalsa 3)
    305         x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq);
    306         y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq);
    307        
    308         xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1);
    309         ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1);
    310         u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1);
    311         v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1);
    312        
    313        
    314         [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys
    315         %remove NaN
    316         tempNaN=isnan(Zphys);tempind=find(tempNaN==1);
    317         Zphys(tempind)=[];
    318         Xphys(tempind)=[];
    319         Yphys(tempind)=[];
    320         error(tempind)=[];
    321          
    322     end
    323    
    324             if NbView>2   
    325        ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd
    326             end
    327    
    328 end
    329     ZI=mean(ZItemp,3); %mean between two the two time step
    330     Vtest=ZItemp(:,:,2)-ZItemp(:,:,1);
     306        if NbView>2
     307            ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd
     308            CheckZ=1;
     309        end     
     310    end
     311    ZI=mean(ZItemp,3); %mean between the two time step
    331312   
    332313    [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,[],XI,YI,ZI);% set of image coordinates on view a
    333314    [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,[],XI,YI,ZI);% set of image coordinates on view b
    334315   
    335    
     316    
    336317    for iview=1:2
    337318        %% reading input file(s)
     
    342323        end
    343324        % get the time defined in the current file if not already defined from the xml file
    344         if isfield(Data{iview},'Time')&& isequal(Data{iview}.Time,Data{1}.Time)
     325        if isfield(Data{iview},'Time')&& (Data{iview}.Time-Data{1}.Time)<0.0001
    345326            Time=Data{iview}.Time;
    346327        else
     
    355336        end
    356337    end
    357     %remove wrong vector 
     338    %remove wrong vector
    358339    if isfield(Data{1},'FF')
    359340        temp=find(Data{1}.FF==0);
     
    373354    [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~
    374355   
    375     %remove wrong vector 
     356    %remove wrong vector
    376357    if isfield(Data{2},'FF')
    377358        temp=find(Data{2}.FF==0);
     
    396377    S=ones(size(XI,1),size(XI,2),3);
    397378    D=ones(size(XI,1),size(XI,2),3,3);
    398 
     379   
    399380    S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
    400381    S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
     
    416397            W(indj,indi)=dxyz(3);
    417398        end
    418     end   
     399    end
    419400    Error=zeros(size(XI,1),size(XI,2),4);
    420401    Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua;
     
    424405   
    425406   
    426 
    427    
    428  
    429    
    430407    %% recording the merged field
    431408    if index==1% initiate the structure at first index
    432         MergeData.ListGlobalAttribute={'Conventions','Time','Dt'};
     409        MergeData.ListGlobalAttribute={'Conventions','Time','Dt','CoordUnit'};
    433410        MergeData.Conventions='uvmat';
    434         MergeData.Time=Time;
    435         MergeData.Dt=Dt;
    436         MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'};
     411        if isfield (XmlData{1}.GeometryCalib,'CoordUnit') && isfield (XmlData{2}.GeometryCalib,'CoordUnit') && strcmp(XmlData{1}.GeometryCalib.CoordUnit, XmlData{2}.GeometryCalib.CoordUnit)
     412        MergeData.CoordUnit=XmlData{1}.GeometryCalib.CoordUnit;
     413        else
     414            disp_uvmat('ERROR','inconsistent coord units in the two input velocity series',checkrun)
     415            return
     416        end
     417        MergeData.ListVarName={'coord_x','coord_y','U','V','W','Error'};
    437418        MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
    438                 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     419            {'coord_y','coord_x'},{'coord_y','coord_x'}};
     420        MergeData.VarAttribute{6}.unit='pixel'; %error estimate expressed in pixel
     421        if CheckZ
     422            MergeData.ListVarName=[MergeData.ListVarName {'Z'}];
     423            MergeData.VarDimName=[MergeData.ListVarName {'coord_y','coord_x'}];
     424            MergeData.Z=ZI;
     425        end
    439426        MergeData.coord_x=xI;
    440427        MergeData.coord_y=yI;
    441428    end
     429    MergeData.Time=Time;
     430    MergeData.Dt=Dt;
    442431    MergeData.U=U/Dt;
    443432    MergeData.V=V/Dt;
    444433    MergeData.W=W/Dt;
    445     MergeData.Z=ZI;
    446    
    447 %     mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
    448 %     mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
    449     MergeData.Error=0.5*sqrt(sum(Error.^2,3));
     434   
     435   
     436    mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
     437    mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
     438    MergeData.Error=0.25*(mfx+mfy)*sqrt(sum(Error.^2,3));
    450439    errormsg=struct2nc(OutputFile,MergeData);%save result file
    451440    if isempty(errormsg)
Note: See TracChangeset for help on using the changeset viewer.