[862]  1  %'civ2vel_3C': combine velocity fields from two cameras to get three velocity components


[839]  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)


[863]  57  ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)


[839]  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)


[863]  64  ParamOut.ProjObject='on';%can use projection object(option 'off'/'on',


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


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


[863]  67  ParamOut.OutputSubDirMode='two'; % the two first input lines are used to define the output subfolder


[839]  68  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


[863]  69  %check the input files


[839]  70  first_j=[];


[863]  71  if size(Param.InputTable,1)<2


 72  msgbox_uvmat('WARNING',['two or three input file series are needed'])


 73  end


[839]  74  if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end


 75  PairString='';


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


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


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


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


 80  if ~exist(FirstFileName,'file')


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


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


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


 84  end


 85  return


 86  end


 87 


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


 89  ParamOut=[]; %default output


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


 91  checkrun=1;


 92  if ischar(Param)


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


 94  checkrun=0;


 95  end


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


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


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


 99 


[863]  100 


 101  %% root input file(s) name, type and index series


 102  RootPath=Param.InputTable(:,1);


 103  RootFile=Param.InputTable(:,3);


 104  SubDir=Param.InputTable(:,2);


 105  NomType=Param.InputTable(:,4);


 106  FileExt=Param.InputTable(:,5);


 107  hdisp=disp_uvmat('WAITING...','checking the file series',checkrun);


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


 109  if ~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  %%%%%%%%%%%%


 118  NbView=numel(i1_series);%number of input file series (lines in InputTable)


 119  NbField_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)


 120  NbField_i=size(i1_series{1},2); %nb of fields for the i index


 121  NbField=NbField_j*NbField_i; %total number of fields


 122 


[839]  123  %% define the directory for result file (with path=RootPath{1})


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


[863]  125  %


[839]  126  % if ~isfield(Param,'InputFields')


 127  % Param.InputFields.FieldName='';


 128  % end


 129 


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


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


 132  if size(time,1)>1


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


[863]  134  if diff_time>0


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


[863]  136  end


[839]  137  end


 138  if ~isempty(errormsg)


 139  disp_uvmat('WARNING',errormsg,checkrun)


 140  end


[863]  141  time=mean(time,1); %averaged time taken for the merged field


[839]  142  if isfield(XmlData{1},'GeometryCalib')


[863]  143  tsaiA=XmlData{1}.GeometryCalib;


 144  else


 145  disp_uvmat('ERROR','no geometric calibration available for image A',checkrun)


 146  return


 147  end


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


 149  tsaiB=XmlData{2}.GeometryCalib;


 150  else


 151  disp_uvmat('ERROR','no geometric calibration available for image B',checkrun)


 152  return


 153  end


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


 155 


[862]  156  %% grid of physical positions (given by projection plane)


 157  if ~Param.CheckObject


[863]  158  disp_uvmat('ERROR','a projection plane with interpolation is needed',checkrun)


 159  return


[862]  160  end


 161  ObjectData=Param.ProjObject;


[863]  162  xI=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);


 163  yI=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);


 164  [XI,YI]=meshgrid(xI,yI);


 165  U=zeros(size(XI,1),size(XI,2));


 166  V=zeros(size(XI,1),size(XI,2));


 167  W=zeros(size(XI,1),size(XI,2));


[862]  168 


[863]  169  %% MAIN LOOP ON FIELDS


 170  warning off


[880]  171 


 172  for index=1:NbField


[863]  173  update_waitbar(WaitbarHandle,index/NbField)


[839]  174  if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')


 175  disp('program stopped by user')


 176  return


 177  end


 178 


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


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


 181  timeread=zeros(1,NbView);


[864]  182 


 183  %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken


 184  %at the middle between to time step


 185  clear ZItemp


 186  ZItemp=zeros(size(XI,1),size(XI,2),2);


 187  for indextemp=index:index+1;


 188  if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys


 189 


[880]  190  [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(indextemp),'.nc']);


[864]  191 


[880]  192 


[864]  193  if exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector


 194  temp=find(Data{3}.Civ3_FF==0);


 195  Zphys=Data{3}.Zphys(temp);


 196  Yphys=Data{3}.Yphys(temp);


 197  Xphys=Data{3}.Xphys(temp);


 198  else


 199  Zphys=Data{3}.Zphys;


 200  Yphys=Data{3}.Yphys;


 201  Xphys=Data{3}.Xphys;


 202  end


 203 


 204 


 205 


 206  elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys


 207 


 208 


 209  %test if the seconde camera is the same for both folder


 210  for i=3:4


 211  indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1


 212  indline(i)=strfind(Param.InputTable{i,2},''); % indice of the "" is the folder name1


 213  camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)1);% extract the second camera name


 214  end


 215 


 216  if strcmp(camname{3},camname{4})==0


 217  disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)


 218  return


 219  end


 220 


[880]  221  [Data{3},tild,errormsg] = nc2struct([Param.InputTable{3,1},'/',Param.InputTable{3,2},'/',Param.InputTable{3,3},'_',int2str(indextemp),'.nc']);


[864]  222  if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector


 223  temp=find(Data{3}.Civ3_FF==0);


 224  Xmid3=Data{3}.Xmid(temp);


 225  Ymid3=Data{3}.Ymid(temp);


 226  U3=Data{3}.Uphys(temp);


 227  V3=Data{3}.Vphys(temp);


 228  else


 229  Xmid3=Data{3}.Xmid;


 230  Ymid3=Data{3}.Ymid;


 231  U3=Data{3}.Uphys;


 232  V3=Data{3}.Vphys;


 233  end


 234  %temporary gridd of merging the 2 stereos datas


 235  [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));


 236 


 237  %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera


 238  %(Dalsa 3)


 239  x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3(U3)/2,xq,yq);


 240  y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3(V3)/2,xq,yq);


 241 


 242 


[880]  243  [Data{4},tild,errormsg] = nc2struct([Param.InputTable{4,1},'/',Param.InputTable{4,2},'/',Param.InputTable{4,3},'_',int2str(indextemp),'.nc']);


[864]  244  if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector


 245  temp=find(Data{4}.Civ3_FF==0);


 246  Xmid4=Data{4}.Xmid(temp);


 247  Ymid4=Data{4}.Ymid(temp);


 248  U4=Data{4}.Uphys(temp);


 249  V4=Data{4}.Vphys(temp);


 250  else


 251  Xmid4=Data{4}.Xmid;


 252  Ymid4=Data{4}.Ymid;


 253  U4=Data{4}.Uphys;


 254  V4=Data{4}.Vphys;


 255  end


 256 


 257  %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera


 258  %(Dalsa 3)


 259  x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4(U4)/2,xq,yq);


 260  y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4(V4)/2,xq,yq);


 261 


 262  xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1);


 263  ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1);


 264  u=reshape(x4Qx3Q,length(xq(:,1)).*length(xq(1,:)),1);


 265  v=reshape(y4Qy3Q,length(yq(:,1)).*length(yq(1,:)),1);


 266 


 267 


 268  [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys


 269  %remove NaN


 270  tempNaN=isnan(Zphys);tempind=find(tempNaN==1);


 271  Zphys(tempind)=[];


 272  Xphys(tempind)=[];


 273  Yphys(tempind)=[];


 274  error(tempind)=[];


 275 


 276  end


 277 


 278 


 279 


 280 


 281  ZItemp(:,:,indextemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd


 282 


 283  end


 284  ZI=mean(ZItemp,3); %mean between two the two time step


[878]  285  Vtest=ZItemp(:,:,2)ZItemp(:,:,1);


[864]  286 


[863]  287  [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a


 288  [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b


 289 


[878]  290 


[863]  291  for iview=1:2


[839]  292  %% reading input file(s)


[863]  293  [Data{iview},tild,errormsg]=read_civdata(filecell{iview,index},{'vec(U,V)'},'*');


[839]  294  if ~isempty(errormsg)


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


 296  return


 297  end


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


[863]  299  if isfield(Data{iview},'Time')&& isequal(Data{iview}.Time,Data{1}.Time)


 300  Time=Data{iview}.Time;


 301  else


 302  disp_uvmat('ERROR','Time undefined or not synchronous',checkrun)


 303  return


[839]  304  end


[863]  305  if isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Data{1}.Dt)


 306  Dt=Data{iview}.Dt;


 307  else


 308  disp_uvmat('ERROR','Dt undefined or not synchronous',checkrun)


 309  return


[839]  310  end


 311  end


[864]  312  %remove wrong vector


 313  temp=find(Data{1}.FF==0);


 314  X1=Data{1}.X(temp);


 315  Y1=Data{1}.Y(temp);


 316  U1=Data{1}.U(temp);


 317  V1=Data{1}.V(temp);


 318 


 319  Ua=griddata(X1,Y1,U1,Xa,Ya);


 320  Va=griddata(X1,Y1,V1,Xa,Ya);


[863]  321 


[878]  322  [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X


 323  [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~


[864]  324 


[878]  325  %remove wrong vector


[864]  326  temp=find(Data{2}.FF==0);


 327  X2=Data{2}.X(temp);


 328  Y2=Data{2}.Y(temp);


 329  U2=Data{2}.U(temp);


 330  V2=Data{2}.V(temp);


 331  Ub=griddata(X2,Y2,U2,Xb,Yb);


 332  Vb=griddata(X2,Y2,V2,Xb,Yb);


[878]  333 


 334  [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X


 335  [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~


 336 


 337 


 338  % System to solve


[863]  339  S=ones(size(XI,1),size(XI,2),3);


 340  D=ones(size(XI,1),size(XI,2),3,3);


[878]  341 


[863]  342  S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;


 343  S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;


 344  S(:,:,3)=A(:,:,1,3).*Ua+A(:,:,2,3).*Va+B(:,:,1,3).*Ub+B(:,:,2,3).*Vb;


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


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


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


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


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


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


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


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


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


 354  for indj=1:size(XI,1)


 355  for indi=1:size(XI,2)


[878]  356  dxyz=(squeeze(D(indj,indi,:,:))*1000)\(squeeze(S(indj,indi,:))*1000); % solving...


[863]  357  U(indj,indi)=dxyz(1);


 358  V(indj,indi)=dxyz(2);


 359  W(indj,indi)=dxyz(3);


 360  end


 361  end


 362  Error=zeros(size(XI,1),size(XI,2),4);


 363  Error(:,:,1)=A(:,:,1,1).*U+A(:,:,1,2).*V+A(:,:,1,3).*WUa;


 364  Error(:,:,2)=A(:,:,2,1).*U+A(:,:,2,2).*V+A(:,:,2,3).*WVa;


 365  Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*WUb;


 366  Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*WVb;


 367 


[878]  368 


 369 


 370 


[839]  371  %% generating the name of the merged field


 372  i1=i1_series{1}(index);


 373  if ~isempty(i2_series{end})


 374  i2=i2_series{end}(index);


 375  else


 376  i2=i1;


 377  end


 378  j1=1;


 379  j2=1;


 380  if ~isempty(j1_series{1})


 381  j1=j1_series{1}(index);


 382  if ~isempty(j2_series{end})


 383  j2=j2_series{end}(index);


 384  else


 385  j2=j1;


 386  end


 387  end


[863]  388  OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_12',i1,i2,j1,j2);


 389 


[839]  390  %% recording the merged field


[863]  391  if index==1% initiate the structure at first index


 392  MergeData.ListGlobalAttribute={'Conventions','Time','Dt'};


[839]  393  MergeData.Conventions='uvmat';


[863]  394  MergeData.Time=Time;


 395  MergeData.Dt=Dt;


 396  MergeData.ListVarName={'coord_x','coord_y','Z','U','V','W','Error'};


 397  MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...


 398  {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};


 399  MergeData.coord_x=xI;


 400  MergeData.coord_y=yI;


 401  MergeData.Z=ZI;


 402  end


 403  MergeData.U=U/Dt;


 404  MergeData.V=V/Dt;


 405  MergeData.W=W/Dt;


[878]  406 


 407  mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;


 408  mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;


 409  MergeData.Error=(sqrt(mfx^2+mfy^2)/4).*sqrt(sum(Error.*Error,3));


[863]  410  errormsg=struct2nc(OutputFile,MergeData);%save result file


 411  if isempty(errormsg)


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


 413  else


 414  disp(errormsg)


 415  end


[839]  416  end


 417 


 418 


[878]  419  function [A]=get_coeff(Calib,X,Y,x,y,z) % compute A~ coefficients


[863]  420  R=(Calib.R)';%rotation matrix


 421  T_z=Calib.Tx_Ty_Tz(3);


 422  T=R(7)*x+R(8)*y+R(9)*z+T_z;


[878]  423 


[863]  424  A(:,:,1,1)=(R(1)R(7)*X)./T;


 425  A(:,:,1,2)=(R(2)R(8)*X)./T;


 426  A(:,:,1,3)=(R(3)R(9)*X)./T;


 427  A(:,:,2,1)=(R(4)R(7)*Y)./T;


 428  A(:,:,2,2)=(R(5)R(8)*Y)./T;


 429  A(:,:,2,3)=(R(6)R(9)*Y)./T;


[839]  430 


[878]  431  function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U


[839]  432 


[878]  433  X1d=XdUd/2;


 434  X2d=Xd+Ud/2;


 435  Y1d=YdVd/2;


 436  Y2d=Yd+Vd/2;


[839]  437 


[878]  438  X1=(X1dCalib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(2).*(X1dCalib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(2).*(Y1dCalib.Cx_Cy(2)).^2 ).^(1);


 439  X2=(X2dCalib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(2).*(X2dCalib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(2).*(Y2dCalib.Cx_Cy(2)).^2 ).^(1);


 440  Y1=(Y1dCalib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(2).*(X1dCalib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(2).*(Y1dCalib.Cx_Cy(2)).^2 ).^(1);


 441  Y2=(Y2dCalib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(2).*(X2dCalib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(2).*(Y2dCalib.Cx_Cy(2)).^2 ).^(1);


[839]  442 


[878]  443  U=X2X1;


 444  V=Y2Y1;


 445  X=X1+U/2;


 446  Y=Y1+V/2;


 447 


 448 


 449 


 450  function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) % get H from stereo data


[864]  451  z=0;


 452  error=0;


[839]  453 


 454 


[864]  455  %% first image


 456  Calib_A=XmlData{1}.GeometryCalib;


 457  R=(Calib_A.R)';


 458  x_a=xmid u/2;


[878]  459  y_a=ymid v/2;


[864]  460  z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);


 461  Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;


 462  Ya=(R(4)*x_a+R(5)*y_a+Calib_A.Tx_Ty_Tz(1,2))./z_a;


[839]  463 


[864]  464  A_1_1=R(1)R(7)*Xa;


 465  A_1_2=R(2)R(8)*Xa;


 466  A_1_3=R(3)R(9)*Xa;


 467  A_2_1=R(4)R(7)*Ya;


 468  A_2_2=R(5)R(8)*Ya;


 469  A_2_3=R(6)R(9)*Ya;


 470  Det=A_1_1.*A_2_2A_1_2.*A_2_1;


 471  Dxa=(A_1_2.*A_2_3A_2_2.*A_1_3)./Det;


 472  Dya=(A_2_1.*A_1_3A_1_1.*A_2_3)./Det;


 473 


 474  %% second image


[878]  475  %loading shift angle


 476 


[864]  477  Calib_B=XmlData{2}.GeometryCalib;


 478  R=(Calib_B.R)';


[878]  479 


 480 


[864]  481  x_b=xmid+ u/2;


 482  y_b=ymid+ v/2;


 483  z_b=R(7)*x_b+R(8)*y_b+Calib_B.Tx_Ty_Tz(1,3);


 484  Xb=(R(1)*x_b+R(2)*y_b+Calib_B.Tx_Ty_Tz(1,1))./z_b;


 485  Yb=(R(4)*x_b+R(5)*y_b+Calib_B.Tx_Ty_Tz(1,2))./z_b;


 486  B_1_1=R(1)R(7)*Xb;


 487  B_1_2=R(2)R(8)*Xb;


 488  B_1_3=R(3)R(9)*Xb;


 489  B_2_1=R(4)R(7)*Yb;


 490  B_2_2=R(5)R(8)*Yb;


 491  B_2_3=R(6)R(9)*Yb;


 492  Det=B_1_1.*B_2_2B_1_2.*B_2_1;


 493  Dxb=(B_1_2.*B_2_3B_2_2.*B_1_3)./Det;


 494  Dyb=(B_2_1.*B_1_3B_1_1.*B_2_3)./Det;


 495 


 496  %% result


 497  Den=(DxbDxa).*(DxbDxa)+(DybDya).*(DybDya);


[878]  498  error=((DybDya).*(u)(DxbDxa).*(v))./Den;


 499  % ex=error.*(DybDya);


 500  % ey=error.*(DxbDxa);


[864]  501 


[878]  502  % z1=u./(DxbDxa);


 503  % z2=v./(DybDya);


 504  z=((DxbDxa).*(u)+(DybDya).*(v))./Den;


 505 


[864]  506  xnew(1,:)=Dxa.*z+x_a;


 507  xnew(2,:)=Dxb.*z+x_b;


 508  ynew(1,:)=Dya.*z+y_a;


 509  ynew(2,:)=Dyb.*z+y_b;


[878]  510  Xphy=mean(xnew,1);


 511  Yphy=mean(ynew,1);


[864]  512 


 513 


 514 


[878]  515 

