[839]  1  %'civ2vel_3C': combine velocity fields from two camerasto get three velocity components


 2  %


 3  % function ParamOut=civ2vel_3C(Param)


 4  %


 5  %OUTPUT


 6  % ParamOut: sets options in the GUI series.fig needed for the function


 7  %


 8  %INPUT:


 9  % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.


 10  % In batch mode, Param is the name of the corresponding xml file containing the same information


 11  % when Param.Action.RUN=0 (as activated when the current Action is selected


 12  % in series), the function ouput paramOut set the activation of the needed GUI elements


 13  %


 14  % Param contains the elements:(use the menu bar command 'export/GUI config' in series to


 15  % see the current structure Param)


 16  % .InputTable: cell of input file names, (several lines for multiple input)


 17  % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}


 18  % .OutputSubDir: name of the subdirectory for data outputs


 19  % .OutputDirExt: directory extension for data outputs


 20  % .Action: .ActionName: name of the current activated function


 21  % .ActionPath: path of the current activated function


 22  % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct


 23  % .RUN =0 for GUI input, =1 for function activation


 24  % .RunMode='local','background', 'cluster': type of function use


 25  %


 26  % .IndexRange: set the file or frame indices on which the action must be performed


 27  % .InputFields: sub structure describing the input fields withfields


 28  % .FieldName: name(s) of the field


 29  % .VelType: velocity type


 30  % .FieldName_1: name of the second field in case of two input series


 31  % .VelType_1: velocity type of the second field in case of two input series


 32  % .Coord_y: name of y coordinate variable


 33  % .Coord_x: name of x coordinate variable'


 34 


 35  %=======================================================================


 36  % Copyright 20082014, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


 37  % http://www.legi.grenobleinp.fr


 38  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 39  %


 40  % This file is part of the toolbox UVMAT.


 41  %


 42  % UVMAT is free software; you can redistribute it and/or modify


 43  % it under the terms of the GNU General Public License as published


 44  % by the Free Software Foundation; either version 2 of the license,


 45  % or (at your option) any later version.


 46  %


 47  % UVMAT is distributed in the hope that it will be useful,


 48  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 49  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 50  % GNU General Public License (see LICENSE.txt) for more details.


 51  %=======================================================================


 52 


 53  function ParamOut=civ2vel_3C(Param)


 54 


 55  %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed


 56  if 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)


 58  ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)


 59  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)


 61  ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)


 62  ParamOut.FieldTransform = 'off';%use the phys transform function without choice


 63  %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',


 65  ParamOut.Mask='off';%can use mask option (option 'off'/'on', 'off' by default)


 66  ParamOut.OutputDirExt='.vel3C';%set the output dir extension


 67  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  first_j=[];


 70  if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end


 71  PairString='';


 72  if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end


 73  [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);


 74  FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...


 75  Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);


 76  if ~exist(FirstFileName,'file')


 77  msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])


 78  elseif isequal(size(Param.InputTable,1),1) && ~isfield(Param,'ProjObject')


 79  msgbox_uvmat('WARNING','You may need a projection object of type plane for merge_proj')


 80  end


 81  return


 82  end


 83 


 84  %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%


 85  ParamOut=[]; %default output


 86  %% read input parameters from an xml file if input is a file name (batch mode)


 87  checkrun=1;


 88  if ischar(Param)


 89  Param=xml2struct(Param);% read Param as input file (batch case)


 90  checkrun=0;


 91  end


 92  hseries=findobj(allchild(0),'Tag','series');


 93  RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series


 94  WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series


 95 


 96  %% define the directory for result file (with path=RootPath{1})


 97  OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files


 98  %


 99  % if ~isfield(Param,'InputFields')


 100  % Param.InputFields.FieldName='';


 101  % end


 102 


 103  %% calibration data and timing: read the ImaDoc files


 104  [XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);


 105  if size(time,1)>1


 106  diff_time=max(max(diff(time)));


 107  if diff_time>0


 108  disp_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time) ': the mean time is chosen in result'],checkrun)


 109  end


 110  end


 111  if ~isempty(errormsg)


 112  disp_uvmat('WARNING',errormsg,checkrun)


 113  end


 114  time=mean(time,1); %averaged time taken for the merged field


 115  if isfield(XmlData{1},'GeometryCalib')


 116  tsaiA=XmlData{1}.GeometryCalib;


 117  else


 118  msgbox_uvmat('ERROR','no geometric calibration available for image A')


 119  return


 120  end


 121  if isfield(XmlData{2},'GeometryCalib')


 122  tsaiB=XmlData{2}.GeometryCalib;


 123  else


 124  msgbox_uvmat('ERROR','no geometric calibration available for image B')


 125  return


 126  end


 127  [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);


 128 


 129  %% MAIN LOOP ON FIELDS


 130  for index=1:NbField


 131  update_waitbar(WaitbarHandle,index/NbField)


 132  if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')


 133  disp('program stopped by user')


 134  return


 135  end


 136 


 137  %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%


 138  Data=cell(1,NbView);%initiate the set Data


 139  timeread=zeros(1,NbView);


 140  for iview=1:NbView


 141  %% reading input file(s)


 142  [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));


 143  if ~isempty(errormsg)


 144  disp_uvmat('ERROR',['ERROR in civ2vel_3C/read_field/' errormsg],checkrun)


 145  return


 146  end


 147  % get the time defined in the current file if not already defined from the xml file


 148  if ~isempty(time) && isfield(Data{iview},'Time')


 149  timeread(iview)=Data{iview}.Time;


 150  end


 151  if ~isempty(NbSlice_calib)


 152  Data{iview}.ZIndex=mod(i1_series{iview}(index)1,NbSlice_calib{iview})+1;%Zindex for phys transform


 153  end


 154 


 155  %% transform the input field (e.g; phys) if requested (no transform involving two input fields)


 156  Data{iview}=phys(Data{iview},XmlData{iview});


 157 


 158  %% projection on object (gridded plane)


 159  % if Param.CheckObject


 160  [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);


 161  if ~isempty(errormsg)


 162  disp_uvmat('ERROR',['ERROR in merge_proge/proj_field: ' errormsg],checkrun)


 163  return


 164  end


 165 


 166  end


 167  %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%%


 168 


 169  %% merge the NbView fields


 170  [MergeData,errormsg]=merge_field(Data);


 171  if ~isempty(errormsg)


 172  disp_uvmat('ERROR',errormsg,checkrun);


 173  return


 174  end


 175 


 176  %% time of the merged field: take the average of the different views


 177  if ~isempty(time)


 178  timeread=time(index);


 179  elseif ~isempty(find(timeread))% time defined from ImaDoc


 180  timeread=mean(timeread(timeread~=0));% take average over times form the files (when defined)


 181  else


 182  timeread=index;% take time=file index


 183  end


 184 


 185  %% generating the name of the merged field


 186  i1=i1_series{1}(index);


 187  if ~isempty(i2_series{end})


 188  i2=i2_series{end}(index);


 189  else


 190  i2=i1;


 191  end


 192  j1=1;


 193  j2=1;


 194  if ~isempty(j1_series{1})


 195  j1=j1_series{1}(index);


 196  if ~isempty(j2_series{end})


 197  j2=j2_series{end}(index);


 198  else


 199  j2=j1;


 200  end


 201  end


 202  OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomType{1},i1,i2,j1,j2);


 203 


 204  %% recording the merged field


 205 


 206  MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim'};


 207  MergeData.Conventions='uvmat';


 208  MergeData.nb_coord=2;


 209  MergeData.nb_dim=2;


 210  dt=[];


 211  if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)


 212  dt=Data{1}.dt;


 213  end


 214  for iview =2:numel(Data)


 215  if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))


 216  dt=[];%dt not the same for all fields


 217  end


 218  end


 219  if ~isempty(timeread)


 220  MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];


 221  MergeData.Time=timeread;


 222  end


 223  if ~isempty(dt)


 224  MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'dt'}];


 225  MergeData.dt=dt;


 226  end


 227  error=struct2nc(OutputFile,MergeData);%save result file


 228  if isempty(error)


 229  disp(['output file ' OutputFile ' written'])


 230  else


 231  disp(error)


 232  end


 233  end


 234 


 235  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 236  % %read the velocity fields


 237  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 238  %read field A


 239  [Field,VelTypeOut]=read_civxdata(file_A,[],vel_type);


 240  %removes false vectors


 241  if isfield(Field,'FF')


 242  Field.X=Field.X(find(Field.FF==0));


 243  Field.Y=Field.Y(find(Field.FF==0));


 244  Field.U=Field.U(find(Field.FF==0));


 245  Field.V=Field.V(find(Field.FF==0));


 246  end


 247  %interpolate on the grid common to both images in phys coordinates


 248  dXa= griddata_uvmat(Field.X,Field.Y,Field.U,XimaA,YimaA);


 249  dYa= griddata_uvmat(Field.X,Field.Y,Field.V,XimaA,YimaA);


 250  dt=Field.dt;


 251  time=Field.Time;


 252 


 253  %read field B


 254  [Field,VelTypeOut]=read_civxdata(file_B,[],vel_type);


 255  if ~isequal(Field.dt,dt)


 256  msgbox_uvmat('ERROR','different time intervals for the two velocity fields ')


 257  return


 258  end


 259  if ~isequal(Field.Time,time)


 260  msgbox_uvmat('ERROR','different times for the two velocity fields ')


 261  return


 262  end


 263  %removes false vectors


 264  if isfield(Field,'FF')


 265  Field.X=Field.X(find(Field.FF==0));


 266  Field.Y=Field.Y(find(Field.FF==0));


 267  Field.U=Field.U(find(Field.FF==0));


 268  Field.V=Field.V(find(Field.FF==0));


 269  end


 270  %interpolate on XimaB


 271  dXb=griddata_uvmat(Field.X,Field.Y,Field.U,XimaB,YimaB);


 272  dYb=griddata_uvmat(Field.X,Field.Y,Field.V,XimaB,YimaB);


 273  %eliminate NotaNumber


 274  ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));


 275  dXa=dXa(ind_Nan);


 276  dYa=dYa(ind_Nan);


 277  dXb=dXb(ind_Nan);


 278  dYb=dYb(ind_Nan);


 279  grid_phys1(:,1)=grid_real_x(ind_Nan);


 280  grid_phys1(:,2)=grid_real_y(ind_Nan);


 281 


 282  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 283  %compute the differential coefficients of the geometric calibration


 284  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 285  [A11,A12,A13,A21,A22,A23]=pxcm_tsai(tsaiA,grid_phys1);


 286  [B11,B12,B13,B21,B22,B23]=pxcm_tsai(tsaiB,grid_phys1);


 287 


 288  C1=A11.*A22A12.*A21;


 289  C2=A13.*A22A12.*A23;


 290  C3=A13.*A21A11.*A23;


 291  D1=B11.*B22B12.*B21;


 292  D2=B13.*B22B12.*B23;


 293  D3=B13.*B21B11.*B23;


 294  A1=(A22.*D1.*(C1.*D3C3.*D1)+A21.*D1.*(C2.*D1C1.*D2));


 295  A2=(A12.*D1.*(C3.*D1C1.*D3)+A11.*D1.*(C1.*D2C2.*D1));


 296  B1=(B22.*C1.*(C3.*D1C1.*D3)+B21.*C1.*(C1.*D2C2.*D1));


 297  B2=(B12.*C1.*(C1.*D3C3.*D1)+B11.*C1.*(C2.*D1C1.*D2));


 298  Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);


 299 


 300  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 301  %Projection for compatible displacements


 302  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 303  Ua=dXaLambda.*A1;


 304  Va=dYaLambda.*A2;


 305  Ub=dXbLambda.*B1;


 306  Vb=dYbLambda.*B2;


 307 


 308  %%%%%%%%%%%%%%%%%%%%%%%%%%%%


 309  %Calculations of displacements and error


 310  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 311  U=(A22.*D2.*UaA12.*D2.*VaB22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2C2.*D1);


 312  V=(A21.*D3.*UaA11.*D3.*VaB21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1C1.*D3);


 313  W=(A22.*D1.*UaA12.*D1.*VaB22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1C1.*D2);


 314  W1=(A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*UbB11.*C1.*Vb)./(C1.*D3C3.*D1);


 315 


 316  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));


 317 


 318  ind_error=(find(error<thresh_patch));


 319  U=U(ind_error);


 320  V=V(ind_error);


 321  W=W(ind_error);%correction for water interface


 322  error=error(ind_error);


 323 


 324  %create nc grid file


 325  Result.ListGlobalAttribute={'nb_coord','nb_dim','constant_pixcm','absolut_time_T0','hart','dt','civ'};


 326  Result.nb_coord=3;%grid file, no velocity


 327  Result.nb_dim=2;


 328  Result.constant_pixcm=0;%no linear correspondance with images


 329  Result.absolut_time_T0=time;%absolute time of the field


 330  Result.hart=0;


 331  Result.dt=dt;%time interval for image correlation (put by default)


 332  % cte.title='grid';


 333  Result.civ=0;%not a civ file (no direct correspondance with an image)


 334  % Result.ListDimName={'nb_vectors'}


 335  % Result.DimValue=length(U);


 336  Result.ListVarName={'vec_X','vec_Y','vec_U','vec_V','vec_W','vec_E'};


 337  Result.VarDimName={'nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors','nb_vectors'}


 338  Result.vec_X= grid_phys1(ind_error,1);


 339  Result.vec_Y= grid_phys1(ind_error,2);


 340  Result.vec_U=U/dt;


 341  Result.vec_V=V/dt;


 342  Result.vec_W=W/dt;


 343  Result.vec_E=error;


 344  % error=write_netcdf(file_st,cte,fieldlabels,grid_phys);


 345  error=struct2nc(file_st,Result);


 346  display([file_st ' written'])


 347 


 348 


 349 


 350  %'pxcm_tsai': find differentials of the Tsai calibration


 351  function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)


 352  R=(a.R)';


 353 


 354  x=var_phys(:,1);


 355  y=var_phys(:,2);


 356 


 357  if isfield(a,'PlanePos')


 358  prompt={'Plane 1 Index','Plane 2 Index'};


 359  Rep=inputdlg(prompt,'Target displacement test');


 360  Z1=str2double(Rep(1));


 361  Z2=str2double(Rep(2));


 362  z=(a.PlanePos(Z2,3)+a.PlanePos(Z1,3))/2


 363  else


 364  z=0;


 365  end


 366 


 367  %transform coeff for differentiels


 368  a.C11=R(1)*R(8)R(2)*R(7);


 369  a.C12=R(2)*R(7)R(1)*R(8);


 370  a.C21=R(4)*R(8)R(5)*R(7);


 371  a.C22=R(5)*R(7)R(4)*R(8);


 372  a.C1x=R(3)*R(7)R(9)*R(1);


 373  a.C1y=R(3)*R(8)R(9)*R(2);


 374  a.C2x=R(6)*R(7)R(9)*R(4);


 375  a.C2y=R(6)*R(8)R(9)*R(5);


 376 


 377  % %dependence in x,y


 378  % denom=(R(7)*x+R(8)*y+R(9)*z+a.Tz).*(R(7)*x+R(8)*y+R(9)*z+a.Tz);


 379  % A11=(a.f*a.sx*(a.C11*ya.C1x*z+R(1)*a.TzR(7)*a.Tx)./denom)/a.dpx;


 380  % A12=(a.f*a.sx*(a.C12*xa.C1y*z+R(2)*a.TzR(8)*a.Tx)./denom)/a.dpx;


 381  % A21=(a.f*a.sx*(a.C21*ya.C2x*z+R(4)*a.TzR(7)*a.Ty)./denom)/a.dpy;


 382  % A22=(a.f*(a.C22*xa.C2y*z+R(5)*a.TzR(8)*a.Ty)./denom)/a.dpy;


 383  % A13=(a.f*(a.C1x*x+a.C1y*y+R(3)*a.TzR(9)*a.Tx)./denom)/a.dpx;


 384  % A23=(a.f*(a.C2x*x+a.C2y*y+R(6)*a.TzR(9)*a.Ty)./denom)/a.dpy;


 385 


 386  %dependence in x,y


 387  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));


 388  A11=(a.fx_fy(1)*(a.C11*ya.C1x*z+R(1)*a.Tx_Ty_Tz(3)R(7)*a.Tx_Ty_Tz(1))./denom);


 389  A12=(a.fx_fy(1)*(a.C12*xa.C1y*z+R(2)*a.Tx_Ty_Tz(3)R(8)*a.Tx_Ty_Tz(1))./denom);


 390  A21=(a.fx_fy(1)*(a.C21*ya.C2x*z+R(4)*a.Tx_Ty_Tz(3)R(7)*a.Tx_Ty_Tz(2))./denom);


 391  A22=(a.fx_fy(2)*(a.C22*xa.C2y*z+R(5)*a.Tx_Ty_Tz(3)R(8)*a.Tx_Ty_Tz(2))./denom);


 392  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);


 393  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);


 394 


 395 


 396 


 397 


 398 


 399 


 400 


 401 

