[172]  1  % transform image coordinates (px) to polar physical coordinates


[40]  2  %[DataOut,DataOut_1]=phys_polar(varargin)


 3  %


 4  % OUTPUT:


 5  % DataOut: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components


 6  % DataOut_1: second data field (if two fields are in input)


 7  %


 8  %INPUT:


 9  % Data: structure of input data (like UvData)


 10  % CalibData= structure containing the field .GeometryCalib with calibration parameters


 11  % Data_1: second input field (not mandatory)


 12  % CalibData_1= calibration parameters for the second field


 13 


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


 15  Calib{1}=[];


 16  if nargin==2nargin==4


 17  Data=varargin{1};


 18  DataOut=Data;%default


 19  DataOut_1=[];%default


 20  CalibData=varargin{2};


 21  if isfield(CalibData,'GeometryCalib')


 22  Calib{1}=CalibData.GeometryCalib;


 23  end


 24  Calib{2}=Calib{1};


 25  else


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


 27  end


 28  test_1=0;


[93]  29  if nargin==4% case of two input fields


[40]  30  test_1=1;


 31  Data_1=varargin{3};


 32  DataOut_1=Data_1;%default


 33  CalibData_1=varargin{4};


 34  if isfield(CalibData_1,'GeometryCalib')


 35  Calib{2}=CalibData_1.GeometryCalib;


 36  end


 37  end


 38 


 39  %parameters for polar coordinates (taken from the calibration data of the first field)


 40  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 41  origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates


 42  if isfield(Calib{1},'PolarCentre') && isnumeric(Calib{1}.PolarCentre)


 43  if isequal(length(Calib{1}.PolarCentre),2);


 44  origin_xy= Calib{1}.PolarCentre;


 45  end


 46  end


 47  radius_offset=0;%reference radius used to offset the radial coordinate r


 48  angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)


 49  if isfield(Calib{1},'PolarReferenceRadius') && isnumeric(Calib{1}.PolarReferenceRadius)


 50  radius_offset=Calib{1}.PolarReferenceRadius;


 51  end


 52  if radius_offset > 0


 53  angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius


 54  else


 55  angle_scale=180/pi; %polar angle in degrees


 56  end


 57  if isfield(Calib{1},'PolarReferenceAngle') && isnumeric(Calib{1}.PolarReferenceAngle)


 58  angle_offset=Calib{1}.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))


 59  end


 60  % new x coordinate = radiusradius_offset;


 61  % new y coordinate = theta*angle_scaleangle_offset


 62 


 63  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 64 


 65  iscalar=0;


[93]  66  %transform first field to cartesian phys coordiantes


[40]  67  if ~isempty(Calib{1})


 68  DataOut=phys_1(Data,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);


 69  %case of images or scalar


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


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


 72  iscalar=1;


 73  A{1}=Data.A;


 74  end


 75  %transform of X,Y coordinates for vector fields


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


 77  ZIndex=Data.ZIndex;


 78  else


 79  ZIndex=0;


 80  end


 81  end


[93]  82  %transform second field (if exists) to cartesian phys coordiantes


[40]  83  if test_1


 84  DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale);


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


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


 87  iscalar=iscalar+1;


 88  Calib{iscalar}=Calib{2};


 89  A{iscalar}=Data_1.A;


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


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


 92  end


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


[164]  94  [A,AX,AY]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);


[40]  95  DataOut_1.A=A{1};


 96  DataOut_1.AX=AX;


 97  DataOut_1.AY=AY;


 98  return


 99  end


 100  end


 101  end


 102  if iscalar~=0


[164]  103  [A,AX,AY]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%


[40]  104  DataOut.A=A{1};


 105  DataOut.AX=AX;


 106  DataOut.AY=AY;


 107  if iscalar==2


 108  DataOut_1.A=A{2};


 109  DataOut_1.AX=AX;


 110  DataOut_1.AY=AY;


 111  end


 112  end


 113 


[161]  114 


 115 


 116 


[40]  117  %


 118  function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)


 119 


 120  DataOut=Data;


[167]  121  % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates


[161]  122  if isfield(Calib,'SliceCoord')


 123  DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane


[40]  124  end


[161]  125 


 126  if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib)


 127  if isfield(Calib,'CoordUnit')


 128  DataOut.CoordUnit=Calib.CoordUnit;


 129  else


 130  DataOut.CoordUnit='cm'; %default


[40]  131  end


[161]  132  DataOut.TimeUnit='s';


[40]  133  %transform of X,Y coordinates for vector fields


[161]  134  if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)


[40]  135  Z=Data.ZIndex;


 136  else


 137  Z=0;


 138  end


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


 140  [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical


 141  DataOut.X=DataOut.Xorigin_xy(1);%origin of coordinates at the tank center


 142  DataOut.Y=DataOut.Yorigin_xy(2);%origin of coordinates at the tank center


 143  [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta and X are the polar coordinates angle and radius


 144  %shift and renormalize the polar coordinates


 145  DataOut.X=DataOut.Xradius_offset;%


 146  DataOut.Y=theta*angle_scaleangle_offset;% normalized angle: distance along reference radius


 147  %transform velocity field if exists


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


 149  if ~isempty(Data.dt)


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


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


 152  UX=(XOut_2XOut_1)/Data.dt;


 153  VY=(YOut_2YOut_1)/Data.dt;


 154  %transform u,v into polar coordiantes


 155  DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity


 156  DataOut.V=(UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);%angular velocity calculated


 157  %shift and renormalize the angular velocity


 158  end


 159  end


[93]  160  %transform of spatial derivatives


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


 162  && isfield(Data,'dt')


 163  if ~isempty(Data.dt)


 164  % estimate the Jacobian matrix DXpx/DXphys


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


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


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


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


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


 170  %Jacobian matrix DXpphys/DXpx


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


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


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


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


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


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


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


 178  end


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


 180  end


 181  end


[40]  182  end


 183  end


 184 


[164]  185 


[157]  186  %


 187  %'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters


 188  % function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)


 189  %


 190  %OUTPUT:


 191  %


 192  %INPUT:


 193  %Z: index of plane


 194  function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)


 195  %


 196  if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z


 197  Zindex=Z;


 198  Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE


 199  else


 200  Zphys=0;


 201  end


 202  if ~exist('X','var')~exist('Y','var')


 203  Xphys=[];


 204  Yphys=[];%default


 205  return


 206  end


 207  %coordinate transform


 208  if ~isfield(Calib,'fx_fy')


 209  Calib.fx_fy=[1 1];


 210  end


 211  if ~isfield(Calib,'Tx_Ty_Tz')


 212  Calib.Tx_Ty_Tz=[0 0 1];


 213  end


 214  if ~isfield(Calib,'Cx_Cy')


 215  Calib.Cx_Cy=[0 0];


 216  end


 217  if ~isfield(Calib,'kc')


 218  Calib.kc=0;


 219  end


 220  if isfield(Calib,'R')


 221  R=(Calib.R)';


 222  Tx=Calib.Tx_Ty_Tz(1);


 223  Ty=Calib.Tx_Ty_Tz(2);


 224  Tz=Calib.Tx_Ty_Tz(3);


 225  f=Calib.fx_fy(1);%dpy=1; sx=1


 226  dpx=Calib.fx_fy(2)/Calib.fx_fy(1);


 227  Dx=R(5)*R(7)R(4)*R(8);


 228  Dy=R(1)*R(8)R(2)*R(7);


 229  D0=f*(R(2)*R(4)R(1)*R(5));


 230  Z11=R(6)*R(8)R(5)*R(9);


 231  Z12=R(2)*R(9)R(3)*R(8);


 232  Z21=R(4)*R(9)R(6)*R(7);


 233  Z22=R(3)*R(7)R(1)*R(9);


 234  Zx0=R(3)*R(5)R(2)*R(6);


 235  Zy0=R(1)*R(6)R(3)*R(4);


 236  A11=R(8)*TyR(5)*Tz+Z11*Zphys;


 237  A12=R(2)*TzR(8)*Tx+Z12*Zphys;


 238  A21=R(7)*Ty+R(4)*Tz+Z21*Zphys;


 239  A22=R(1)*Tz+R(7)*Tx+Z22*Zphys;


 240  X0=f*(R(5)*TxR(2)*Ty+Zx0*Zphys);


 241  Y0=f*(R(4)*Tx+R(1)*Ty+Zy0*Zphys);


 242  %px to camera:


 243  Xd=dpx*(XCalib.Cx_Cy(1)); % sensor coordinates


 244  Yd=(YCalib.Cx_Cy(2));


 245  dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd)/(f*f); %distortion factor


 246  Xu=Xd./dist_fact;%undistorted sensor coordinates


 247  Yu=Yd./dist_fact;


 248  denom=Dx*Xu+Dy*Yu+D0;


 249  Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates


 250  Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;


 251  else


 252  Xphys=Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1);


 253  Yphys=Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2);


 254  end


[164]  255 


[40]  256  %%%%%%%%%%%%%%%%%%%%


[164]  257  function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)


[40]  258  xcorner=[];


 259  ycorner=[];


 260  npx=[];


 261  npy=[];


 262  for icell=1:length(A)


 263  siz=size(A{icell});


 264  npx=[npx siz(2)];


 265  npy=[npy siz(1)];


 266  zphys=0; %default


 267  if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane


 268  SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);


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


 270  end


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


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


 273  [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates


 274  %transform the corner coordinates into polar ones


 275  xcorner_new=xcorner_neworigin_xy(1);%shift to the origin of the polar coordinates


 276  ycorner_new=ycorner_neworigin_xy(2);%shift to the origin of the polar coordinates


 277  [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta and X are the polar coordinates angle and radius


 278  if (max(theta)min(theta))>pi %if the polar origin is inside the image


 279  xcorner_new=[0 max(xcorner_new)];


 280  theta=[pi pi];


 281  end


 282  %shift and renormalize the polar coordinates


 283  xcorner_new=xcorner_newradius_offset;%


 284  ycorner_new=theta*angle_scaleangle_offset;% normalized angle: distance along reference radius


 285  xcorner=[xcorner xcorner_new];


 286  ycorner=[ycorner ycorner_new];


 287  end


 288  Rangx(1)=min(xcorner);


 289  Rangx(2)=max(xcorner);


 290  Rangy(2)=min(ycorner);


 291  Rangy(1)=max(ycorner);


 292  % test_multi=(max(npx)~=min(npx))  (max(npy)~=min(npy));


 293  npx=max(npx);


 294  npy=max(npy);


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


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


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


 298  %transform X, Y in cartesian


 299  X=X+radius_offset;%


 300  Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius


 301  [X,Y] = pol2cart(Y,X);


 302  X=X+origin_xy(1);%shift to the origin of the polar coordinates


 303  Y=Y+origin_xy(2);%shift to the origin of the polar coordinates


 304  for icell=1:length(A)


[164]  305  siz=size(A{icell});


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


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


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


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


[164]  310  if numel(siz)==2 %(B/W images)


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


 312  ind_in=find(flagin);


 313  ind_out=find(~flagin);


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


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


 316  vec_B(ind_in)=vec_A(ICOMB);


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


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


 319  else


 320  for icolor=1:siz(3)


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


 322  ind_in=find(flagin);


 323  ind_out=find(~flagin);


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


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


 326  vec_B(ind_in)=vec_A(ICOMB);


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


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


 329  end


 330  end


[40]  331  end


 332 

