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


 2  %


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


 4  % OUTPUT:


 5  % DataOut: output field structure


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


[810]  7  %


[574]  8  %INPUT:


 9  % DataIn: first input field structure


 10  % XmlData: first input parameter structure,


 11  % .GeometryCalib: substructure of the calibration parameters


 12  % DataIn_1: optional second input field structure


 13  % XmlData_1: optional second input parameter structure


 14  % .GeometryCalib: substructure of the calibration parameters


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


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


 17  %


 18  % OUTPUT:


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


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


 21  %


 22  %INPUT:


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


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


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


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


[810]  27 


 28  %=======================================================================


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


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


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


 32  %


 33  % This file is part of the toolbox UVMAT.


 34  %


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


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


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


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


 39  %


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


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


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


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


 44  %=======================================================================


 45 


[574]  46  function DataOut=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)


 47  %


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


 49  if nargin==2nargin==4


[658]  50  DataOut=DataIn;%default


[40]  51  DataOut_1=[];%default


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


 53  Calib{1}=XmlData.GeometryCalib;


[40]  54  end


 55  Calib{2}=Calib{1};


 56  else


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


 58  end


 59  test_1=0;


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


[40]  61  test_1=1;


[658]  62  DataOut_1=DataIn_1;%default


 63  if isfield(XmlData_1,'GeometryCalib')


 64  Calib{2}=XmlData_1.GeometryCalib;


[40]  65  end


 66  end


 67 


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


 69  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


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


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


 73  origin_xy= XmlData.PolarCentre;


[40]  74  end


 75  end


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


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


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


 79  radius_offset=XmlData.PolarReferenceRadius;


[40]  80  end


 81  if radius_offset > 0


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


 83  else


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


 85  end


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


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


[40]  88  end


 89  % new x coordinate = radiusradius_offset;


 90  % new y coordinate = theta*angle_scaleangle_offset


 91 


 92  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 93 


 94  iscalar=0;


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


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


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


 98  %case of images or scalar


[782]  99  if isfield(Data,'A')&isfield(Data,'Coord_x')&~isempty(Data.Coord_x) & isfield(Data,'Coord_y')&...


 100  ~isempty(Data.Coord_y)&length(Data.A)>1


[40]  101  iscalar=1;


 102  A{1}=Data.A;


 103  end


 104  %transform of X,Y coordinates for vector fields


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


 106  ZIndex=Data.ZIndex;


 107  else


 108  ZIndex=0;


 109  end


 110  end


[567]  111 


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


[40]  113  if test_1


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


[782]  115  if isfield(Data_1,'A')&isfield(Data_1,'Coord_x')&~isempty(Data_1.Coord_x) & isfield(Data_1,'Coord_y')&...


 116  ~isempty(Data_1.Coord_y)&length(Data_1.A)>1


[40]  117  iscalar=iscalar+1;


 118  Calib{iscalar}=Calib{2};


 119  A{iscalar}=Data_1.A;


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


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


 122  end


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


[782]  124  [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);


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


[782]  126  DataOut_1.Coord_x=Coord_x;


 127  DataOut_1.Coord_y=Coord_y;


[40]  128  return


 129  end


 130  end


 131  end


 132  if iscalar~=0


[782]  133  [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%


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


[782]  135  DataOut.Coord_x=Coord_x;


 136  DataOut.Coord_y=Coord_y;


[40]  137  if iscalar==2


 138  DataOut_1.A=A{2};


[782]  139  DataOut_1.Coord_x=Coord_x;


 140  DataOut_1.Coord_y=Coord_y;


[40]  141  end


 142  end


 143 


[161]  144 


 145 


 146 


[40]  147  %


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


 149 


 150  DataOut=Data;


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


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


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


[40]  154  end


[161]  155 


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


 157  if isfield(Calib,'CoordUnit')


 158  DataOut.CoordUnit=Calib.CoordUnit;


 159  else


 160  DataOut.CoordUnit='cm'; %default


[40]  161  end


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


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


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


[40]  165  Z=Data.ZIndex;


 166  else


 167  Z=0;


 168  end


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


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


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


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


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


 174  %shift and renormalize the polar coordinates


 175  DataOut.X=DataOut.Xradius_offset;%


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


 177  %transform velocity field if exists


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


 179  if ~isempty(Data.dt)


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


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


 182  UX=(XOut_2XOut_1)/Data.dt;


 183  VY=(YOut_2YOut_1)/Data.dt;


 184  %transform u,v into polar coordiantes


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


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


 187  %shift and renormalize the angular velocity


 188  end


 189  end


[93]  190  %transform of spatial derivatives


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


 192  && isfield(Data,'dt')


 193  if ~isempty(Data.dt)


 194  % estimate the Jacobian matrix DXpx/DXphys


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


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


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


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


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


 200  %Jacobian matrix DXpphys/DXpx


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


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


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


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


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


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


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


 208  end


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


 210  end


 211  end


[40]  212  end


 213  end


 214 


[164]  215 


[40]  216  %%%%%%%%%%%%%%%%%%%%


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


[40]  218  xcorner=[];


 219  ycorner=[];


 220  npx=[];


 221  npy=[];


 222  for icell=1:length(A)


 223  siz=size(A{icell});


 224  npx=[npx siz(2)];


 225  npy=[npy siz(1)];


 226  zphys=0; %default


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


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


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


 230  end


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


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


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


 234  %transform the corner coordinates into polar ones


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


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


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


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


 239  xcorner_new=[0 max(xcorner_new)];


 240  theta=[pi pi];


 241  end


 242  %shift and renormalize the polar coordinates


 243  xcorner_new=xcorner_newradius_offset;%


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


 245  xcorner=[xcorner xcorner_new];


 246  ycorner=[ycorner ycorner_new];


 247  end


 248  Rangx(1)=min(xcorner);


 249  Rangx(2)=max(xcorner);


 250  Rangy(2)=min(ycorner);


 251  Rangy(1)=max(ycorner);


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


 253  npx=max(npx);


 254  npy=max(npy);


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


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


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


 258  %transform X, Y in cartesian


 259  X=X+radius_offset;%


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


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


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


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


 264  for icell=1:length(A)


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


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


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


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


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


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


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


 272  ind_in=find(flagin);


 273  ind_out=find(~flagin);


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


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


 276  vec_B(ind_in)=vec_A(ICOMB);


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


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


 279  else


 280  for icolor=1:siz(3)


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


 282  ind_in=find(flagin);


 283  ind_out=find(~flagin);


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


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


 286  vec_B(ind_in)=vec_A(ICOMB);


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


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


 289  end


 290  end


[40]  291  end


 292 

