[2]  1  %'phys': transforms image (px) to real world (phys) coordinates using geometric calibration parameters


 2  % OUTPUT:


 3  % DataOut: structure representing the modified field


 4  % DataOut_1: structure representing the second modified field


 5 


 6  %INPUT:


 7  % Data: structure of input data


 8  % with fields .A (image or scalar matrix), AX, AY


 9  % .X,.Y,.U,.V, .DjUi


 10  % .ZIndex: index of plane in multilevel case


 11  % Data.CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px'


 12  % Calib: structure containing calibration parameters or a subtree Calib.GeometryCalib =calibration data (tsai parameters)


 13 


 14  function [DataOut,DataOut_1]=phys(varargin)


 15  % A FAIRE: 1 verifier si DataIn est une 'field structure'(.ListVarName'):


 16  % chercher ListVarAttribute, for each field (cell of variables):


 17  % .CoordType: 'phys' or 'px' (default==phys, no transform)


 18  % .scale_factor: =dt (to transform displacement into velocity) default=1


 19  % .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)


 20  % (default='coord' if .Role='coord_x,_y...,


 21  % 'D_i' if '.Role='vector_x,...',


 22  % 'scalar', else (thenno change except scale factor)


 23  Calib{1}=[];


 24  if nargin==2nargin==4 % nargin =nbre of input variables


 25  Data=varargin{1};


 26  DataOut=Data;%default


 27  DataOut_1=[];%default


 28  CalibData=varargin{2};


 29  if isfield(CalibData,'GeometryCalib')


 30  Calib{1}=CalibData.GeometryCalib;


 31  end


 32  Calib{2}=Calib{1};


 33  else


 34  DataOut.Txt='wrong input: need two or four structures';


 35  end


 36  test_1=0;


 37  if nargin==4


 38  test_1=1;


 39  Data_1=varargin{3};


 40  DataOut_1=Data_1;%default


 41  CalibData_1=varargin{4};


 42  if isfield(CalibData_1,'GeometryCalib')


 43  Calib{2}=CalibData_1.GeometryCalib;


 44  end


 45  end


 46  iscalar=0;


 47  if ~isempty(Calib{1})


 48  DataOut=phys_1(Data,Calib{1});


 49  %case of images or scalar: in case of two input fields, we need to project the transform of on the same regular grid


 50  if isfield(Data,'A') && isfield(Data,'AX') && ~isempty(Data.AX) && isfield(Data,'AY')&&...


 51  ~isempty(Data.AY) && length(Data.A)>1


 52  iscalar=1;


 53  A{1}=Data.A;


 54  end


 55  end


 56  %transform of X,Y coordinates for vector fields


 57  if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)


 58  ZIndex=Data.ZIndex;


 59  else


 60  ZIndex=0;


 61  end


 62  if test_1


 63  DataOut_1=phys_1(Data_1,Calib{2});


 64  if isfield(Data_1,'A')&&isfield(Data_1,'AX')&&~isempty(Data_1.AX) && isfield(Data_1,'AY')&&...


 65  ~isempty(Data_1.AY)&&length(Data_1.A)>1


 66  iscalar=iscalar+1;


 67  Calib{iscalar}=Calib{2};


 68  A{iscalar}=Data_1.A;


 69  if isfield(Data_1,'ZIndex') && ~isequal(Data_1.ZIndex,ZIndex)


 70  DataOut.Txt='inconsistent plane indexes in the two input fields';


 71  end


 72  if iscalar==1% case for which only the second field is a scalar


 73  [A,AX,AY]=phys_Ima(A,Calib,ZIndex);


 74  DataOut_1.A=A{1};


 75  DataOut_1.AX=AX;


 76  DataOut_1.AY=AY;


 77  return


 78  end


 79  end


 80  end


 81  if iscalar~=0


 82  [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima


 83  DataOut.A=A{1};


 84  DataOut.AX=AX;


 85  DataOut.AY=AY;


 86  if iscalar==2


 87  DataOut_1.A=A{2};


 88  DataOut_1.AX=AX;


 89  DataOut_1.AY=AY;


 90  end


 91  end


 92 


 93  %


 94  function DataOut=phys_1(Data,Calib)


 95  % for icell=1:length(Data)


 96 


 97  DataOut=Data;%default


 98  DataOut.CoordType='phys'; %put flag for physical coordinates


 99  % The transform ACTS ONLY IF .CoordType='px'and Calib defined


 100  if isfield(Data,'CoordType')&& isequal(Data.CoordType,'px')&& ~isempty(Calib)


 101  if isfield(Calib,'CoordUnit')


 102  DataOut.CoordUnit=Calib.CoordUnit;


 103  else


 104  DataOut.CoordUnit='cm'; %default


 105  % elseif isfield(DataOut,'CoordUnit')


 106  % DataOut=rmfield(DataOut,'CoordUnit');


 107  end


 108  DataOut.TimeUnit='s';


 109  %transform of X,Y coordinates for vector fields


 110  if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)


 111  Z=Data.ZIndex;


 112  else


 113  Z=0;


 114  end


 115  if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)


 116  [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z);


 117  if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt')


 118  if ~isempty(Data.dt)


 119  [XOut_1,YOut_1]=phys_XYZ(Calib,Data.XData.U/2,Data.YData.V/2,Z);


 120  [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);


 121  DataOut.U=(XOut_2XOut_1)/Data.dt;


 122  DataOut.V=(YOut_2YOut_1)/Data.dt;


 123  end


 124  end


 125  end


 126  %transform of an image or scalar: done in phys_ima


 127 


 128  %transform of spatial derivatives


 129  if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...


 130  && isfield(Data,'dt')


 131  if ~isempty(Data.dt)


 132  % estimate the Jacobian matrix DXpx/DXphys


 133  for ip=1:length(Data.X)


 134  [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);


 135  [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)0.5,Data.Y(ip),Z);


 136  [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);


 137  [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)0.5,Z);


 138  %Jacobian matrix DXpphys/DXpx


 139  DjXi(1,1)=(Xp1Xm1);


 140  DjXi(2,1)=(Yp1Ym1);


 141  DjXi(1,2)=(Xp2Xm2);


 142  DjXi(2,2)=(Yp2Ym2);


 143  DjUi(:,:)=Data.DjUi(ip,:,:);


 144  DjUi=(DjXi*DjUi')/DjXi;% =J1*M*J , curvature effects (derivatives of J) neglected


 145  DataOut.DjUi(ip,:,:)=DjUi';


 146  end


 147  DataOut.DjUi = DataOut.DjUi/Data.dt; % min(Data.DjUi(:,1,1))=DUDX


 148  end


 149  end


 150  end


 151 


 152 


 153  %%%%%%%%%%%%%%%%%%%%


 154  function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)


 155  xcorner=[];


 156  ycorner=[];


 157  npx=[];


 158  npy=[];


 159  for icell=1:length(A)


 160  siz=size(A{icell});


 161  npx=[npx siz(2)];


 162  npy=[npy siz(1)];


 163  Calib=CalibIn{icell};


 164  xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordiantes of corners


 165  yima=[0.5 0.5 siz(1)0.5 siz(1)0.5];


 166  [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates


 167  xcorner=[xcorner xcorner_new];


 168  ycorner=[ycorner ycorner_new];


 169  end


 170  Rangx(1)=min(xcorner);


 171  Rangx(2)=max(xcorner);


 172  Rangy(2)=min(ycorner);


 173  Rangy(1)=max(ycorner);


 174  test_multi=(max(npx)~=min(npx))  (max(npy)~=min(npy));


 175  npx=max(npx);


 176  npy=max(npy);


 177  x=linspace(Rangx(1),Rangx(2),npx);


 178  y=linspace(Rangy(1),Rangy(2),npy);


 179  [X,Y]=meshgrid(x,y);%grid in physical coordiantes


 180  vec_B=[];


 181  A_out={};


 182  for icell=1:length(A)


 183  Calib=CalibIn{icell};


 184  if (isfield(Calib,'R') && ~isequal(Calib.R(2,1),0) && ~isequal(Calib.R(1,2),0)) ...


 185  ((isfield(Calib,'kappa1')&& ~isequal(Calib.kappa1,0)))  test_multi  ~isequal(Calib,CalibIn{1})


 186  zphys=0; %default


 187  if isfield(Calib,'SliceCoord') %.Z= index of plane


 188  SliceCoord=Calib.SliceCoord(ZIndex,:);


 189  zphys=SliceCoord(3); %to generalize for nonparallel planes


 190  end


 191  [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid


 192  XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'


 193  YIMA=reshape(round(YIMA),1,npx*npy);


 194  flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image


 195  testuint8=isa(A{icell},'uint8');


 196  testuint16=isa(A{icell},'uint16');


 197  if numel(siz)==2 %(B/W images)


 198  vec_A=reshape(A{icell},1,npx*npy);%put the original image in line


 199  ind_in=find(flagin);


 200  ind_out=find(~flagin);


 201  ICOMB=((XIMA1)*npy+(npy+1YIMA));


 202  ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A


 203  vec_B(ind_in)=vec_A(ICOMB);


 204  vec_B(ind_out)=zeros(size(ind_out));


 205  A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates


 206  elseif numel(siz)==3


 207  for icolor=1:siz(3)


 208  vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line


 209  ind_in=find(flagin);


 210  ind_out=find(~flagin);


 211  ICOMB=((XIMA1)*npy+(npy+1YIMA));


 212  ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A


 213  vec_B(ind_in)=vec_A(ICOMB);


 214  vec_B(ind_out)=zeros(size(ind_out));


 215  A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates


 216  end


 217  end


 218  if testuint8


 219  A_out{icell}=uint8(A_out{icell});


 220  end


 221  if testuint16


 222  A_out{icell}=uint16(A_out{icell});


 223  end


 224  else%


 225 


 226  A_out{icell}=A{icell};%no transform


 227  Rangx=[0.5 npx0.5];%image coordiantes of corners


 228  Rangy=[npy0.5 0.5];


 229  [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[ZIndex ZIndex]);%case of translations without rotation and quadratic deformation


 230  [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[ZIndex ZIndex]);


 231  end


 232  end

