Ignore:
Timestamp:
Jun 22, 2016, 1:36:50 PM (5 years ago)
Author:
sommeria
Message:

update calib modfied + various updates

File:
1 edited

Legend:

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

    r927 r954  
    1 %'civ2vel_3C': combine velocity fields from two cameras to get three velocity components
     1%'civ2vel_3C': combine velocity fields from two cameras to get the three velocity components
     2% used with the GUI 'series':
     3%           first input line =raw PIV camera 1 (image coordinates)
     4%           second input line=raw PIV camera 2 (image coordinates)
     5% Three modes:
     6%   1) no additional input: measurements assumed in the reference plane (laser sheet)
     7%   2) measurement surface obtained by stereoscopic comparison of the images of the two cameras.
     8%           third input line =surface z(x,y) given by correlation between camera 1 and 2 (expressed in phys apparent coordinates)
     9%   3)  surface z(x,y) given by adding the displacements of each camera with a third intermediate camera (#3) used to reduce the
     10% to reduce the angle for stereoscopic view.
     11%           third input line =correlation between camera 1 and 3 (expressed in phys apparent coordinates)
     12%           fourth input line =corelation between camera 2 and 3 (expressed in phys apparent coordinates)
     13%               
    214%------------------------------------------------------------------------
    315% function ParamOut=civ2vel_3C(Param)
     
    719%
    820%INPUT:
     21%
    922% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
    1023% In batch mode, Param is the name of the corresponding xml file containing the same information
     
    1528% see the current structure Param)
    1629%    .InputTable: cell of input file names, (several lines for multiple input)
    17 %                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
     30%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}.
     31%                3 or 4 lines used as described above
    1832%    .OutputSubDir: name of the subdirectory for data outputs
    1933%    .OutputDirExt: directory extension for data outputs
     
    2236%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
    2337%             .RUN =0 for GUI input, =1 for function activation
    24 %             .RunMode='local','background', 'cluster': type of function  use
    25 %             
     38%             .RunMode='local','background', 'cluster': type of function  use         
    2639%    .IndexRange: set the file or frame indices on which the action must be performed
    2740%    .InputFields: sub structure describing the input fields withfields
     
    5265
    5366function ParamOut=civ2vel_3C(Param)
    54 disp('test')
     67
    5568%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed
    5669if isstruct(Param) && isequal(Param.Action.RUN,0)
     
    124137%% define the directory for result file (with path=RootPath{1})
    125138OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    126 %
    127 % if ~isfield(Param,'InputFields')
    128 %     Param.InputFields.FieldName='';
    129 % end
    130139
    131140%% calibration data and timing: read the ImaDoc files
     
    160169    return
    161170end
    162 ObjectData=Param.ProjObject;
     171ObjectData=Param.ProjObject;% Object attached to the GUI series
    163172xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
    164173yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
     
    168177W=zeros(size(XI,1),size(XI,2));
    169178
    170 %% MAIN LOOP ON FIELDS
     179%%%% MAIN LOOP ON FIELDS %%%%%
    171180warning off
    172181
     
    175184    CheckOverwrite=Param.CheckOverwrite;
    176185end
    177 for index=1:NbField
    178    
    179     update_waitbar(WaitbarHandle,index/NbField)
    180    
    181    
    182    
    183    
    184       %% generating the name of the merged field
    185     i1=i1_series{1}(index);
     186for field_index=1:NbField
     187   
     188    update_waitbar(WaitbarHandle,field_index/NbField)% waitbar to visualise progress (in RUN mode)
     189   
     190    %% generating the name of the output file for the merged field
     191    i1=i1_series{1}(field_index);
    186192    if ~isempty(i2_series{end})
    187         i2=i2_series{end}(index);
     193        i2=i2_series{end}(field_index);
    188194    else
    189195        i2=i1;
     
    192198    j2=1;
    193199    if ~isempty(j1_series{1})
    194         j1=j1_series{1}(index);
     200        j1=j1_series{1}(field_index);
    195201        if ~isempty(j2_series{end})
    196             j2=j2_series{end}(index);
     202            j2=j2_series{end}(field_index);
    197203        else
    198204            j2=j1;
     
    201207    OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1-2',i1,i2,j1,j2);
    202208   
    203     %%
    204    
    205    
     209    %% program stop if requested on the GUI
    206210    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    207211        disp('program stopped by user')
     
    209213    end
    210214   
    211      if (~CheckOverwrite && exist(OutputFile,'file')) 
    212             disp('existing output file already exists, skip to next field')
    213             continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    214      end   
    215      
     215    %% check existence of the output file
     216    if (~CheckOverwrite && exist(OutputFile,'file'))
     217        disp('existing output file already exists, skip to next field')
     218        continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     219    end
     220   
    216221    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    217222    Data=cell(1,NbView);%initiate the set Data
    218     timeread=zeros(1,NbView);
    219223   
    220224    %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken
    221225    %at the middle between to time step
    222    clear ZItemp
    223    ZItemp=zeros(size(XI,1),size(XI,2),2);
    224    
    225    if index==1
     226    ZItemp=zeros(size(XI,1),size(XI,2),2);
     227    if field_index==1
    226228        first_img=i1_series{1,1}(1,1); %id of the first image of the series
    227    end
    228      
    229      idtemp=0;
    230  for indextemp=index:index+1;
    231      idtemp=idtemp+1;
    232     if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys
    233      
    234        
    235        
    236         [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
    237        
    238         if  exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector
    239             temp=find(Data{3}.Civ3_FF==0);
    240             Zphys=Data{3}.Zphys(temp);
    241             Yphys=Data{3}.Yphys(temp);
    242             Xphys=Data{3}.Xphys(temp);
    243         else
    244             Zphys=Data{3}.Zphys;
    245             Yphys=Data{3}.Yphys;
    246             Xphys=Data{3}.Xphys;
     229    end
     230   
     231    idtemp=0;
     232    % get the surface shape corresponding to the PIV measurements
     233    for indextemp=field_index:field_index+1;%TODO: generalise to field index intervals>1 for PIV
     234        idtemp=idtemp+1;
     235        InputFile_3=fullfile(Param.InputTable{3,1},Param.InputTable{3,2},[Param.InputTable{3,3} '_' int2str(first_img+indextemp-1) '.nc']);
     236        if NbView==3 % if there is only 1 stereo folder (2 cameras only), extract directly Xphys,Yphys and Zphys 
     237            % Data{1}: =raw PIV camera 1 only
     238            % Data{2}: =raw PIV camera 2 only
     239            % Data{3}: =correlation between camera 1 and 2
     240            [Data{3},tild,errormsg] = nc2struct(InputFile_3);%read input file       
     241            if  exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector
     242                temp=find(Data{3}.Civ3_FF==0);
     243                Zphys=Data{3}.Zphys(temp);
     244                Yphys=Data{3}.Yphys(temp);
     245                Xphys=Data{3}.Xphys(temp);
     246            else
     247                Zphys=Data{3}.Zphys;
     248                Yphys=Data{3}.Yphys;
     249                Xphys=Data{3}.Xphys;
     250            end
     251           
     252        elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys
     253            % Data{1}: =raw PIV camera 1 only (left)
     254            % Data{2}: =raw PIV camera 2 only (right) (no PIV done with middle camera)
     255            % Data{3}: =corelation between camera 1 and 3 (left and middle)
     256            % Data{4}: =corelation between camera 2 and 3 (right and middle)
     257            % test if the second camera (3) is the same for both folders
     258            for i=3:4
     259                indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1
     260                indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1
     261                camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name
     262            end       
     263            if strcmp(camname{3},camname{4})==0
     264                disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)
     265                return
     266            end     
     267            [Data{3},tild,errormsg] = nc2struct(InputFile_3);       
     268            if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector
     269                temp=find(Data{3}.Civ3_FF==0);
     270                Xmid3=Data{3}.Xmid(temp);
     271                Ymid3=Data{3}.Ymid(temp);
     272                U3=Data{3}.Uphys(temp);
     273                V3=Data{3}.Vphys(temp);
     274            else
     275                Xmid3=Data{3}.Xmid;
     276                Ymid3=Data{3}.Ymid;
     277                U3=Data{3}.Uphys;
     278                V3=Data{3}.Vphys;
     279            end
     280            %temporary grid of merging the 2 stereos data
     281            [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));
     282           
     283            %1st folder : interpolate the first camera points on the second (common) camera
     284            %(Dalsa 3)
     285            x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq);
     286            y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq);
     287           
     288            InputFile_4=fullfile(Param.InputTable{4,1},Param.InputTable{4,2},[Param.InputTable{4,3} '_' int2str(first_img+indextemp-1) '.nc']);
     289            [Data{4},tild,errormsg] = nc2struct(InputFile_4);
     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            %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera
     303            %(Dalsa 3)
     304            x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq);
     305            y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq);
     306           
     307            %add the displacements of the two camera pairs
     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            % get the surface z(x,y) from the combined displacement
     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           
    247322        end
    248323       
     324        ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen grid
    249325       
    250        
    251     elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys
    252        
    253        
    254         %test if the seconde camera is the same for both folder
    255         for i=3:4
    256         indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1
    257         indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1
    258         camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name
    259         end
    260        
    261         if strcmp(camname{3},camname{4})==0
    262             disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)
    263             return
    264         end
    265        
    266    
    267        
    268         [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(first_img+indextemp-1),'.nc']);
    269    
    270         if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector
    271             temp=find(Data{3}.Civ3_FF==0);
    272             Xmid3=Data{3}.Xmid(temp);
    273             Ymid3=Data{3}.Ymid(temp);
    274             U3=Data{3}.Uphys(temp);
    275             V3=Data{3}.Vphys(temp);
    276         else
    277             Xmid3=Data{3}.Xmid;
    278             Ymid3=Data{3}.Ymid;
    279             U3=Data{3}.Uphys;
    280             V3=Data{3}.Vphys;
    281         end
    282         %temporary gridd of merging the 2 stereos datas
    283         [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));
    284        
    285         %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera
    286         %(Dalsa 3)
    287         x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq);
    288         y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq);
    289        
    290        
    291 
    292          [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(first_img+indextemp-1),'.nc']);
    293         if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector
    294             temp=find(Data{4}.Civ3_FF==0);
    295             Xmid4=Data{4}.Xmid(temp);
    296             Ymid4=Data{4}.Ymid(temp);
    297             U4=Data{4}.Uphys(temp);
    298             V4=Data{4}.Vphys(temp);
    299         else
    300             Xmid4=Data{4}.Xmid;
    301             Ymid4=Data{4}.Ymid;
    302             U4=Data{4}.Uphys;
    303             V4=Data{4}.Vphys;
    304         end
    305        
    306         %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera
    307         %(Dalsa 3)
    308         x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq);
    309         y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq);
    310        
    311         xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1);
    312         ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1);
    313         u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1);
    314         v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1);
    315        
    316        
    317         [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys
    318         %remove NaN
    319         tempNaN=isnan(Zphys);tempind=find(tempNaN==1);
    320         Zphys(tempind)=[];
    321         Xphys(tempind)=[];
    322         Yphys(tempind)=[];
    323         error(tempind)=[];
    324          
    325     end
    326    
    327    
    328    
    329    
    330        ZItemp(:,:,idtemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd
    331    
    332 end
    333     ZI=mean(ZItemp,3); %mean between two the two time step
     326    end
     327    ZI=mean(ZItemp,3); %mean between two the two times used for surface measurement
    334328    Vtest=ZItemp(:,:,2)-ZItemp(:,:,1);
    335329   
     
    337331    [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b
    338332   
    339    
     333    
    340334    for iview=1:2
    341         %% reading input file(s)
    342         [Data{iview},tild,errormsg]=read_civdata(filecell{iview,index},{'vec(U,V)'},'*');
     335        %% reading PIV input file(s)
     336        [Data{iview},tild,errormsg]=read_civdata(filecell{iview,field_index},{'vec(U,V)'},'*');
    343337        if ~isempty(errormsg)
    344338            disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)
     
    359353        end
    360354    end
    361     %remove wrong vector
     355    %remove false vectors
    362356    temp=find(Data{1}.FF==0);
    363357    X1=Data{1}.X(temp);
     
    369363    Va=griddata(X1,Y1,V1,Xa,Ya);
    370364   
    371 %     [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X
     365    [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X
    372366    [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~
    373367   
    374     %remove wrong vector
     368    %remove false vectors
    375369    temp=find(Data{2}.FF==0);
    376370    X2=Data{2}.X(temp);
     
    380374    Ub=griddata(X2,Y2,U2,Xb,Yb);
    381375    Vb=griddata(X2,Y2,V2,Xb,Yb);
    382 
    383 %     [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X
     376   
     377    [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X
    384378    [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~
    385    
     379    
    386380   
    387381    % System to solve
    388382    S=ones(size(XI,1),size(XI,2),3);
    389383    D=ones(size(XI,1),size(XI,2),3,3);
    390 
     384   
    391385    S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
    392386    S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
     
    408402            W(indj,indi)=dxyz(3);
    409403        end
    410     end   
     404    end
    411405    Error=zeros(size(XI,1),size(XI,2),4);
    412406    Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*W-Ua;
     
    415409    Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb;
    416410   
    417    
    418 
    419    
    420  
    421    
    422411    %% recording the merged field
    423     if index==1% initiate the structure at first index
     412    if field_index==1% initiate the structure at first index
    424413        MergeData.ListGlobalAttribute={'Conventions','Time','Dt'};
    425414        MergeData.Conventions='uvmat';
     
    428417        MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'};
    429418        MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
    430                 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     419            {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
    431420        MergeData.coord_x=xI;
    432421        MergeData.coord_y=yI;
     
    437426    MergeData.Z=ZI;
    438427   
    439 %     mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
    440 %     mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
     428    %     mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
     429    %     mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
    441430    MergeData.Error=0.5*sqrt(sum(Error.^2,3));
    442431    errormsg=struct2nc(OutputFile,MergeData);%save result file
     
    461450A(:,:,2,3)=(R(6)-R(9)*Y)./T;
    462451
    463 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X  and Ud to U
    464 
     452function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd)
     453% convert image coordinates to view angles, after removal of  quadratic distorsion
     454% input in pixel, output in radians
    465455X1d=Xd-Ud/2;
    466456X2d=Xd+Ud/2;
     
    483473z=0;
    484474error=0;
    485 
    486475
    487476%% first image
Note: See TracChangeset for help on using the changeset viewer.