[574]  1  %'phys_polar': transforms image (Unit='pixel') to polar (phys) coordinates using geometric calibration parameters


 2 


 3  %


 4  %%%% Use the general syntax for transform fields %%%%


 5  % OUTPUT:


 6  % DataOut: output field structure


 7  % .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components


 8 


 9  %INPUT:


 10  % DataIn: first input field structure


 11  % XmlData: first input parameter structure,


 12  % .GeometryCalib: substructure of the calibration parameters


 13  % DataIn_1: optional second input field structure


 14  % XmlData_1: optional second input parameter structure


 15  % .GeometryCalib: substructure of the calibration parameters


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


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


 18  %


 19  % OUTPUT:


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


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


 22  %


 23  %INPUT:


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


[658]  25  % XmlData= structure containing the field .GeometryCalib with calibration parameters


[40]  26  % Data_1: second input field (not mandatory)


[658]  27  % XmlData_1= calibration parameters for the second field


[574]  28  %


 29  function DataOut=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)


 30  %


[40]  31  Calib{1}=[];


 32  if nargin==2nargin==4


[658]  33  DataOut=DataIn;%default


[40]  34  DataOut_1=[];%default


[658]  35  if isfield(XmlData,'GeometryCalib')


 36  Calib{1}=XmlData.GeometryCalib;


[40]  37  end


 38  Calib{2}=Calib{1};


 39  else


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


 41  end


 42  test_1=0;


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


[40]  44  test_1=1;


[658]  45  DataOut_1=DataIn_1;%default


 46  if isfield(XmlData_1,'GeometryCalib')


 47  Calib{2}=XmlData_1.GeometryCalib;


[40]  48  end


 49  end


 50 


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


 52  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


[658]  54  if isfield(XmlData,'PolarCentre') && isnumeric(XmlData.PolarCentre)


 55  if isequal(length(XmlData.PolarCentre),2);


 56  origin_xy= XmlData.PolarCentre;


[40]  57  end


 58  end


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


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


[658]  61  if isfield(XmlData,'PolarReferenceRadius') && isnumeric(XmlData.PolarReferenceRadius)


 62  radius_offset=XmlData.PolarReferenceRadius;


[40]  63  end


 64  if radius_offset > 0


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


 66  else


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


 68  end


[658]  69  if isfield(XmlData,'PolarReferenceAngle') && isnumeric(XmlData.PolarReferenceAngle)


 70  angle_offset=XmlData.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))


[40]  71  end


 72  % new x coordinate = radiusradius_offset;


 73  % new y coordinate = theta*angle_scaleangle_offset


 74 


 75  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 76 


 77  iscalar=0;


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


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


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


 81  %case of images or scalar


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


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


 84  iscalar=1;


 85  A{1}=Data.A;


 86  end


 87  %transform of X,Y coordinates for vector fields


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


 89  ZIndex=Data.ZIndex;


 90  else


 91  ZIndex=0;


 92  end


 93  end


[567]  94 


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


[40]  96  if test_1


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


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


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


 100  iscalar=iscalar+1;


 101  Calib{iscalar}=Calib{2};


 102  A{iscalar}=Data_1.A;


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


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


 105  end


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


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


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


 109  DataOut_1.AX=AX;


 110  DataOut_1.AY=AY;


 111  return


 112  end


 113  end


 114  end


 115  if iscalar~=0


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


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


 118  DataOut.AX=AX;


 119  DataOut.AY=AY;


 120  if iscalar==2


 121  DataOut_1.A=A{2};


 122  DataOut_1.AX=AX;


 123  DataOut_1.AY=AY;


 124  end


 125  end


 126 


[161]  127 


 128 


 129 


[40]  130  %


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


 132 


 133  DataOut=Data;


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


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


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


[40]  137  end


[161]  138 


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


 140  if isfield(Calib,'CoordUnit')


 141  DataOut.CoordUnit=Calib.CoordUnit;


 142  else


 143  DataOut.CoordUnit='cm'; %default


[40]  144  end


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


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


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


[40]  148  Z=Data.ZIndex;


 149  else


 150  Z=0;


 151  end


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


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


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


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


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


 157  %shift and renormalize the polar coordinates


 158  DataOut.X=DataOut.Xradius_offset;%


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


 160  %transform velocity field if exists


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


 162  if ~isempty(Data.dt)


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


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


 165  UX=(XOut_2XOut_1)/Data.dt;


 166  VY=(YOut_2YOut_1)/Data.dt;


 167  %transform u,v into polar coordiantes


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


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


 170  %shift and renormalize the angular velocity


 171  end


 172  end


[93]  173  %transform of spatial derivatives


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


 175  && isfield(Data,'dt')


 176  if ~isempty(Data.dt)


 177  % estimate the Jacobian matrix DXpx/DXphys


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


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


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


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


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


 183  %Jacobian matrix DXpphys/DXpx


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


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


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


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


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


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


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


 191  end


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


 193  end


 194  end


[40]  195  end


 196  end


 197 


[164]  198 


[40]  199  %%%%%%%%%%%%%%%%%%%%


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


[40]  201  xcorner=[];


 202  ycorner=[];


 203  npx=[];


 204  npy=[];


 205  for icell=1:length(A)


 206  siz=size(A{icell});


 207  npx=[npx siz(2)];


 208  npy=[npy siz(1)];


 209  zphys=0; %default


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


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


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


 213  end


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


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


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


 217  %transform the corner coordinates into polar ones


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


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


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


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


 222  xcorner_new=[0 max(xcorner_new)];


 223  theta=[pi pi];


 224  end


 225  %shift and renormalize the polar coordinates


 226  xcorner_new=xcorner_newradius_offset;%


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


 228  xcorner=[xcorner xcorner_new];


 229  ycorner=[ycorner ycorner_new];


 230  end


 231  Rangx(1)=min(xcorner);


 232  Rangx(2)=max(xcorner);


 233  Rangy(2)=min(ycorner);


 234  Rangy(1)=max(ycorner);


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


 236  npx=max(npx);


 237  npy=max(npy);


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


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


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


 241  %transform X, Y in cartesian


 242  X=X+radius_offset;%


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


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


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


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


 247  for icell=1:length(A)


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


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


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


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


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


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


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


 255  ind_in=find(flagin);


 256  ind_out=find(~flagin);


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


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


 259  vec_B(ind_in)=vec_A(ICOMB);


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


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


 262  else


 263  for icolor=1:siz(3)


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


 265  ind_in=find(flagin);


 266  ind_out=find(~flagin);


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


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


 269  vec_B(ind_in)=vec_A(ICOMB);


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


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


 272  end


 273  end


[40]  274  end


 275 

