[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  %=======================================================================


[924]  29  % Copyright 20082016, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[810]  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  %


[933]  48  %% request input parameters


 49  if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)


 50  prompt = {'origin [x y] of polar coordinates';'reference radius';'reference angle(degrees)'};


 51  dlg_title = 'set the parameters for the polar coordinates';


 52  num_lines= 2;


 53  def = { '[0 0]';'0';'0'};


 54  if isfield(XmlData,'TransformInput')


 55  if isfield(XmlData.TransformInput,'PolarCentre')


 56  def{1}=num2str(XmlData.TransformInput.PolarCentre);


 57  end


 58  if isfield(XmlData.TransformInput,'PolarReferenceRadius')


 59  def{2}=num2str(XmlData.TransformInput.PolarReferenceRadius);


 60  end


 61  if isfield(XmlData.TransformInput,'PolarReferenceAngle')


 62  def{3}=num2str(XmlData.TransformInput.PolarReferenceAngle);


 63  end


 64  end


 65  answer = inputdlg(prompt,dlg_title,num_lines,def);


 66  DataOut.TransformInput.PolarCentre=str2num(answer{1});


 67  DataOut.TransformInput.PolarReferenceRadius=str2num(answer{2});


 68  DataOut.TransformInput.PolarReferenceAngle=str2num(answer{3});


[935]  69  if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


 70  DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units


 71  end


[933]  72  return


 73  end


 74 


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


 76  if nargin==2nargin==4


[658]  77  DataOut=DataIn;%default


[40]  78  DataOut_1=[];%default


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


 80  Calib{1}=XmlData.GeometryCalib;


[40]  81  end


 82  Calib{2}=Calib{1};


 83  else


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


 85  end


 86  test_1=0;


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


[40]  88  test_1=1;


[658]  89  DataOut_1=DataIn_1;%default


 90  if isfield(XmlData_1,'GeometryCalib')


 91  Calib{2}=XmlData_1.GeometryCalib;


[40]  92  end


 93  end


 94 


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


 96  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[933]  97  XmlData.PolarReferenceRadius=450;


 98  XmlData.PolarReferenceAngle=450*pi/2;


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


[933]  100  radius_offset=0;%reference radius used to offset the radial coordinate r


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


 102  if isfield(XmlData,'TransformInput')


 103  if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre)


 104  if isequal(length(XmlData.TransformInput.PolarCentre),2);


 105  origin_xy= XmlData.TransformInput.PolarCentre;


 106  end


[40]  107  end


[933]  108  if isfield(XmlData.TransformInput,'PolarReferenceRadius') && isnumeric(XmlData.TransformInput.PolarReferenceRadius)


 109  radius_offset=XmlData.TransformInput.PolarReferenceRadius;


 110  end


 111  if radius_offset > 0


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


 113  else


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


 115  end


 116  if isfield(XmlData.TransformInput,'PolarReferenceAngle') && isnumeric(XmlData.TransformInput.PolarReferenceAngle)


 117  angle_offset=(pi/180)*XmlData.TransformInput.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))


 118  end


[40]  119  end


 120 


 121  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 122 


 123  iscalar=0;


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


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


[880]  126  DataOut=phys_1(DataIn,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);


[40]  127  %case of images or scalar


[880]  128  if isfield(DataIn,'A')&isfield(DataIn,'Coord_x')&~isempty(DataIn.Coord_x) & isfield(DataIn,'Coord_y')&...


 129  ~isempty(DataIn.Coord_y)&length(DataIn.A)>1


[40]  130  iscalar=1;


[880]  131  A{1}=DataIn.A;


[40]  132  end


 133  %transform of X,Y coordinates for vector fields


[880]  134  if isfield(DataIn,'ZIndex')&~isempty(DataIn.ZIndex)


 135  ZIndex=DataIn.ZIndex;


[40]  136  else


 137  ZIndex=0;


 138  end


 139  end


[567]  140 


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


[40]  142  if test_1


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


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


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


[40]  146  iscalar=iscalar+1;


 147  Calib{iscalar}=Calib{2};


 148  A{iscalar}=Data_1.A;


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


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


 151  end


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


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


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


[782]  155  DataOut_1.Coord_x=Coord_x;


 156  DataOut_1.Coord_y=Coord_y;


[40]  157  return


 158  end


 159  end


 160  end


 161  if iscalar~=0


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


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


[782]  164  DataOut.Coord_x=Coord_x;


 165  DataOut.Coord_y=Coord_y;


[40]  166  if iscalar==2


 167  DataOut_1.A=A{2};


[782]  168  DataOut_1.Coord_x=Coord_x;


 169  DataOut_1.Coord_y=Coord_y;


[40]  170  end


 171  end


 172 


[161]  173 


 174 


 175 


[40]  176  %


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


 178 


 179  DataOut=Data;


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


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


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


[40]  183  end


[161]  184 


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


 186  if isfield(Calib,'CoordUnit')


 187  DataOut.CoordUnit=Calib.CoordUnit;


 188  else


 189  DataOut.CoordUnit='cm'; %default


[40]  190  end


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


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


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


[40]  194  Z=Data.ZIndex;


 195  else


 196  Z=0;


 197  end


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


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


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


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


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


 203  %shift and renormalize the polar coordinates


[933]  204  DataOut.X=DataOut.Xradius_offset;%shift the origin of radius, taken as the new X coordinate


 205  DataOut.Y=(thetaangle_offset)*angle_scale;% normalized angle: distance along reference radius,taken as the new Y coordinate


[40]  206  %transform velocity field if exists


[933]  207  if isfield(Data,'U') & isfield(Data,'V') & ~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'Dt')


 208  if ~isempty(Data.Dt)


 209  [XOut_1,YOut_1]=phys_XYZ(Calib,Data.XData.U/2,Data.YData.V/2,Z);% X,Y positions of the vector origin in phys


 210  [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);% X,Y positions of the vector end in phys


 211  UX=(XOut_2XOut_1)/Data.Dt;% phys velocity u component


 212  VY=(YOut_2YOut_1)/Data.Dt; % phys velocity v component


[40]  213  %transform u,v into polar coordiantes


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


[933]  215  DataOut.V=(UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);% azimuthal velocity component


[40]  216  %shift and renormalize the angular velocity


 217  end


 218  end


[93]  219  %transform of spatial derivatives


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


[933]  221  && isfield(Data,'Dt')


 222  if ~isempty(Data.Dt)


[93]  223  % estimate the Jacobian matrix DXpx/DXphys


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


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


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


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


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


 229  %Jacobian matrix DXpphys/DXpx


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


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


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


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


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


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


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


 237  end


[933]  238  DataOut.DjUi = DataOut.DjUi/Data.Dt; % min(Data.DjUi(:,1,1))=DUDX


[93]  239  end


 240  end


[40]  241  end


 242  end


 243 


[164]  244 


[40]  245  %%%%%%%%%%%%%%%%%%%%


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


[40]  247  xcorner=[];


 248  ycorner=[];


 249  npx=[];


 250  npy=[];


 251  for icell=1:length(A)


 252  siz=size(A{icell});


 253  npx=[npx siz(2)];


 254  npy=[npy siz(1)];


 255  zphys=0; %default


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


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


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


 259  end


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


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


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


 263  %transform the corner coordinates into polar ones


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


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


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


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


 268  xcorner_new=[0 max(xcorner_new)];


 269  theta=[pi pi];


 270  end


 271  %shift and renormalize the polar coordinates


 272  xcorner_new=xcorner_newradius_offset;%


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


 274  xcorner=[xcorner xcorner_new];


 275  ycorner=[ycorner ycorner_new];


 276  end


 277  Rangx(1)=min(xcorner);


 278  Rangx(2)=max(xcorner);


 279  Rangy(2)=min(ycorner);


 280  Rangy(1)=max(ycorner);


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


 282  npx=max(npx);


 283  npy=max(npy);


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


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


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


 287  %transform X, Y in cartesian


 288  X=X+radius_offset;%


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


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


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


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


 293  for icell=1:length(A)


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


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


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


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


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


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


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


 301  ind_in=find(flagin);


 302  ind_out=find(~flagin);


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


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


 305  vec_B(ind_in)=vec_A(ICOMB);


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


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


 308  else


 309  for icolor=1:siz(3)


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


 311  ind_in=find(flagin);


 312  ind_out=find(~flagin);


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


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


 315  vec_B(ind_in)=vec_A(ICOMB);


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


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


 318  end


 319  end


[40]  320  end


 321 

