[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:


[1078]  5  % Data: output field structure


[574]  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


[1078]  16  %[Data,Data_1]=phys_polar(varargin)


[40]  17  %


 18  % OUTPUT:


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


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


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


[1071]  29  % Copyright 20082020, 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 


[1078]  46  function Data=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)


[574]  47  %


[1078]  48 


[933]  49  %% request input parameters


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


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


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


 53  num_lines= 2;


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


 55  if isfield(XmlData,'TransformInput')


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


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


 58  end


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


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


 61  end


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


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


 64  end


 65  end


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


[1078]  67  Data.TransformInput.PolarCentre=str2num(answer{1});


 68  Data.TransformInput.PolarReferenceRadius=str2num(answer{2});


 69  Data.TransformInput.PolarReferenceAngle=str2num(answer{3});


 70  % if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


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


 72  % end


[933]  73  return


 74  end


 75 


[1078]  76  %% default outputs


 77  Data=DataIn; %default output


 78  if isfield(Data,'CoordUnit')


 79  Data=rmfield(Data,'CoordUnit');


 80  end


 81  Data.ListVarName = {};


 82  Data.VarDimName={};


 83  Data.VarAttribute={};


 84  DataCell{1}=DataIn;


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


[1078]  86  DataCell{2}=[];%default


 87  checkpixel(1)=0;


 88  if isfield(DataCell{1},'CoorUnit')&& strcmp(DataCell{1}.CoorUnit,'px')


 89  checkpixel(1)=1;


 90  end


[40]  91  if nargin==2nargin==4


[1078]  92  if isfield(XmlData,'GeometryCalib') && ~isempty(XmlData.GeometryCalib)&& checkpixel(1)


[658]  93  Calib{1}=XmlData.GeometryCalib;


[40]  94  end


 95  Calib{2}=Calib{1};


 96  else


[1078]  97  Data.Txt='wrong input: need two or four structures';


[40]  98  end


[1078]  99  nbinput=1;


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


[1078]  101  checkpixel(2)=0;


 102  if isfield(DataCell{2},'CoorUnit')&& strcmp(DataCell{2}.CoorUnit,'px')


 103  checkpixel(2)=1;


 104  end


 105  DataCell{2}=DataIn_1;%default


 106  if isfield(XmlData_1,'GeometryCalib')&& ~isempty(XmlData_1.GeometryCalib) && checkpixel(2)


[658]  107  Calib{2}=XmlData_1.GeometryCalib;


[40]  108  end


[1078]  109  nbinput=2;


[40]  110  end


 111 


[1078]  112  %% parameters for polar coordinates (taken from the calibration data of the first field)


[40]  113  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


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


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


[1078]  117  angle_scale=180/pi;


 118  check_degree=1;%angle expressed in degrees by default


[933]  119  if isfield(XmlData,'TransformInput')


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


[1078]  121  if isequal(length(XmlData.TransformInput.PolarCentre),2)


[933]  122  origin_xy= XmlData.TransformInput.PolarCentre;


 123  end


[40]  124  end


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


 126  radius_offset=XmlData.TransformInput.PolarReferenceRadius;


 127  end


 128  if radius_offset > 0


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


[1078]  130  check_degree=0; %the output has the same unit asthe input


[933]  131  else


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


[1078]  133  check_degree=1;%angle expressed in degrees


[933]  134  end


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


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


 137  end


[40]  138  end


 139 


 140  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[1078]  141  %% get fields


 142  check_scalar=0;


 143  check_vector=0;


 144  nbvar=0;%counter for the number of output variables


 145  nbcoord=0;%counter for the number of variables for radial coordiantes (case of multiple field inputs)


 146  nbgrid=0;%counter for the number of gridded fields (all linearly interpolated on the same output polar grid)


 147  nbscattered=0;%counter of scattered fields


 148  radius_name='radius';


 149  theta_name='theta';


 150  U_r_name='U_r';


 151  U_theta_name='U_theta';


 152  for ifield=1:nbinput %1 or 2 input fields


 153  [CellInfo,NbDim,errormsg]=find_field_cells(DataCell{ifield});


 154  if ~isempty(errormsg)


 155  Data.Txt=['bad input to phys_polar: ' errormsg];


 156  return


[40]  157  end


[1078]  158  % end


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


[1078]  160  if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex)


 161  ZIndex=DataCell{ifield}.ZIndex;


[40]  162  else


 163  ZIndex=0;


 164  end


[1078]  165  for icell=1:numel(CellInfo)


 166  if NbDim(icell)==2


 167  % case of input field with scattered coordinates


 168  if strcmp(CellInfo{icell}.CoordType,'scattered')


 169  nbscattered=nbscattered+1;


 170  nbcoord=nbcoord+1;


 171  radius_name = rename_indexing(radius_name,Data.ListVarName);


 172  theta_name = rename_indexing(theta_name,Data.ListVarName);


 173  Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];


 174  dim_name = rename_indexing('nb_point',Data.VarDimName);


 175  Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];


 176  nbvar=nbvar+2;


 177  Data.VarAttribute{nbvar1}.Role='coord_x';


 178  check_unit=1;


 179  if isfield(DataCell{ifield},'CoordUnit')


 180  Data=rmfield(Data,'CoordUnit');


 181  Data.VarAttribute{nbvar1}.unit=DataCell{ifield}.CoordUnit;


 182  elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


 183  Data.VarAttribute{nbvar1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units


 184  else


 185  check_unit=0;


 186  end


 187  Data.VarAttribute{nbvar}.Role='coord_y';


 188  if check_degree


 189  Data.VarAttribute{nbvar}.unit='degree';


 190  elseif check_unit


 191  Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar1}.unit;


 192  end


 193 


 194  %transform u,v into polar coordinates


 195  X=DataCell{ifield}.(CellInfo{icell}.XName);


 196  Y=DataCell{ifield}.(CellInfo{icell}.YName);


 197  if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')


 198  UName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};


 199  VName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};


 200  if ~isempty(Calib{ifield})


 201  [X,Y,Z,DataCell{ifield}.(UName),DataCell{ifield}.(VName)]=...


 202  phys_XYUV(DataCell{ifield},Calib{ifield},ZIndex);


 203  end


 204  end


 205  [Theta,Radius] = cart2pol(Xorigin_xy(1),Yorigin_xy(2));


 206  Data.(radius_name)=Radiusradius_offset;


 207  Data.(theta_name)=Theta*angle_scaleangle_offset;


 208  if Z~=0


 209  Data.Z=Z;


 210  nbvar=nbvar+1;


 211  Data.ListVarName = [Data.ListVarName {'Z'}];


 212  Data.VarDimName=[Data.VarDimName {dim_name}];


 213  Data.VarAttribute{nbvar}.Role='coord_z';


 214  end


 215  if isfield(CellInfo{icell},'VarIndex_scalar')


 216  ScalarName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};


 217  ScalarName=rename_indexing(ScalarName,Data.ListVarName);


 218  Data.(ScalarName)=DataCell{ifield}.(ScalarName);


 219  nbvar=nbvar+1;


 220  Data.ListVarName = [Data.ListVarName {ScalarName}];


 221  Data.VarDimName=[Data.VarDimName {dim_name}];


 222  Data.VarAttribute{nbvar}.Role='scalar';


 223  end


 224  if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')


 225  U_r_name= rename_indexing(U_r_name,Data.ListVarName);


 226  U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);


 227  Data.(U_r_name)=DataCell{ifield}.(UName).*cos(Theta)+DataCell{ifield}.(VName).*sin(Theta);%radial velocity


 228  Data.(U_theta_name)=(DataCell{ifield}.(UName).*sin(Theta)+DataCell{ifield}.(VName).*cos(Theta));%./(Data.X)%+radius_ref);% azimuthal velocity component


 229  Data.ListVarName = [Data.ListVarName {U_r_name} {U_theta_name}];


 230  Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];


 231  Data.VarAttribute{nbvar+1}.Role='vector_x';


 232  Data.VarAttribute{nbvar+2}.Role='vector_y';


 233  nbvar=nbvar+2;


 234  end


 235  if isfield(CellInfo{icell},'VarIndex_errorflag')


 236  error_flag_name=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_errorflag};


 237  error_flag_newname= rename_indexing(error_flag_name,Data.ListVarName);


 238  Data.(error_flag_newname)=DataCell{ifield}.(error_flag_name);


 239  Data.ListVarName = [Data.ListVarName {error_flag_newname}];


 240  Data.VarDimName=[Data.VarDimName {dim_name}];


 241  nbvar=nbvar+1;


 242  Data.VarAttribute{nbvar}.Role='errorflag';


 243  end


 244 


 245  %caseof input fields on gridded coordinates (matrix)


 246  elseif strcmp(CellInfo{icell}.CoordType,'grid')


 247  if nbgrid==0% no gridded data yet, introduce the coordinate variables common to all gridded data


 248  nbcoord=nbcoord+1;%add new radial coordinates for the first gridded field


 249  radius_name = rename_indexing(radius_name,Data.ListVarName);


 250  theta_name = rename_indexing(theta_name,Data.ListVarName);


 251  Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];


 252  Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}];


 253  nbvar=nbvar+2;


 254  Data.VarAttribute{nbvar1}.Role='coord_x';


 255  Data.VarAttribute{nbvar}.Role='coord_y';


 256  check_unit=1;


 257  if isfield(DataCell{ifield},'CoordUnit')


 258  Data.VarAttribute{nbvar1}.unit=DataCell{ifield}.CoordUnit;


 259  elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


 260  Data.VarAttribute{nbvar1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units


 261  else


 262  check_unit=0;


 263  end


 264  if check_degree


 265  Data.VarAttribute{nbvar}.unit='degree';


 266  elseif check_unit


 267  Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar1}.unit;


 268  end


 269  end


 270  if isfield(CellInfo{icell},'VarIndex_scalar')


 271  nbgrid=nbgrid+1;


 272  nbvar=nbvar+1;


 273  Data.VarAttribute{nbvar}.Role='scalar';


 274  FieldName{nbgrid}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};


 275  A{nbgrid}=DataCell{ifield}.(FieldName{nbgrid});


 276  % Data.ListVarName=[Data.ListVarName {FieldName{nbgrid}}];


 277  % Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];


 278  nbpoint(nbgrid)=numel(A{nbgrid});


 279  check_scalar(nbgrid)=1;


 280  coord_x{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});


 281  coord_y{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});


 282  ZInd(nbgrid)=ZIndex;


 283  Calib_new{nbgrid}=Calib{ifield};


 284  end


 285  if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')


 286  FieldName{nbgrid+1}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};


 287  FieldName{nbgrid+2}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};


 288  A{nbgrid+1}=DataCell{ifield}.(FieldName{nbgrid+1});


 289  A{nbgrid+2}=DataCell{ifield}.(FieldName{nbgrid+2});


 290  % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}];


 291  %Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];


 292  Data.VarAttribute{nbvar+1}.Role='vector_x';


 293  Data.VarAttribute{nbvar+2}.Role='vector_y';


 294  nbpoint([nbgrid+1 nbgrid+2])=numel(A{nbgrid+1});


 295  check_vector(nbgrid+1)=1;


 296  check_vector(nbgrid+2)=1;


 297  coord_x{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});


 298  coord_y{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});


 299  coord_x{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});


 300  coord_y{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});


 301  ZInd(nbgrid+1)=ZIndex;


 302  ZInd(nbgrid+2)=ZIndex;


 303  Calib_new{nbgrid+1}=Calib{ifield};


 304  Calib_new{nbgrid+2}=Calib{ifield};


 305  nbgrid=nbgrid+2;


 306  nbvar=nbvar+2;


 307  end


 308  end


 309  end


 310  end


[40]  311  end


[567]  312 


[1078]  313  %% tranform cartesian to polar coordinates for gridded data


 314  if nbgrid~=0


 315  [A,Data.radius,Data.theta]=phys_Ima_polar(A,coord_x,coord_y,Calib_new,ZInd,origin_xy,radius_offset,angle_offset,angle_scale);


 316  for icell=1:numel(A)


 317  if icell<=numel(A)1 && check_vector(icell)==1 && check_vector(icell+1)==1 %transform u,v into polar coordiantes


 318  theta=Data.theta/angle_scaleangle_offset;


 319  [~,Theta]=meshgrid(Data.radius,theta);%grid in physical coordinates


 320  U_r_name= rename_indexing(U_r_name,Data.ListVarName);


 321  U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);


 322  Data.ListVarName=[Data.ListVarName {U_r_name,U_theta_name}];


 323  Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];


 324  Data.(U_r_name)=A{icell}.*cos(Theta)+A{icell+1}.*sin(Theta);%radial velocity


 325  Data.(U_theta_name)=(A{icell}.*sin(Theta)+A{icell+1}.*cos(Theta));%./(Data.X)%+radius_ref);% azimuthal velocity component


 326  elseif ~check_vector(icell)% for scalar fields


 327  FieldName{icell}= rename_indexing(FieldName{icell},Data.ListVarName);


 328  Data.ListVarName=[Data.ListVarName {FieldName{icell}}];


 329  Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];


 330  Data.(FieldName{icell})=A{icell};


 331  end


[40]  332  end


 333  end


 334 


[161]  335 


[40]  336  %


[1078]  337  % transform a single field into phys coordiantes


 338  function [X,Y,Z,U,V]=phys_XYUV(Data,Calib,ZIndex)


 339  %


 340  %% set default output


 341  %DataOut=Data;%default


 342  %DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters


 343  X=[];%default output


 344  Y=[];


 345  Z=0;


 346  U=[];


 347  V=[];


 348  %% transform X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)


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


 350  [X,Y,Z]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);


 351  Dt=1; %default


 352  if isfield(Data,'dt')&&~isempty(Data.dt)


 353  Dt=Data.dt;


[40]  354  end


[1078]  355  if isfield(Data,'Dt')&&~isempty(Data.Dt)


 356  Dt=Data.Dt;


[40]  357  end


[1078]  358  if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)


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


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


 361  U=(XOut_2XOut_1)/Dt;


 362  V=(YOut_2YOut_1)/Dt;


[40]  363  end


 364  end


 365 


 366  %%%%%%%%%%%%%%%%%%%%


[1078]  367  % tranform gridded field into polar coordiantes on a regular polar grid,


 368  % transform to phys coordiantes if requested by calibration input


 369  function [A_out,radius,theta]=phys_Ima_polar(A,coord_x,coord_y,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)


 370  rcorner=[];


 371  thetacorner=[];


[40]  372  npx=[];


 373  npy=[];


 374  for icell=1:length(A)


 375  siz=size(A{icell});


[1078]  376  npx(icell)=siz(2);


 377  npy(icell)=siz(1);


 378  x_edge=[linspace(coord_x{icell}(1),coord_x{icell}(end),npx(icell)) coord_x{icell}(end)*ones(1,npy(icell))...


 379  linspace(coord_x{icell}(end),coord_x{icell}(1),npx(icell)) coord_x{icell}(1)*ones(1,npy(icell))];%x coordinates of the image edge(four sides)


 380  y_edge=[coord_y{icell}(1)*ones(1,npx(icell)) linspace(coord_y{icell}(1),coord_y{icell}(end),npy(icell))...


 381  coord_y{icell}(end)*ones(1,npx(icell)) linspace(coord_y{icell}(end),coord_y{icell}(1),npy(icell))];%y coordinates of the image edge(four sides)


 382 


 383  % transform edges into phys coordinates if requested


 384  if ~isempty(CalibIn{icell})


 385  [x_edge,y_edge]=phys_XYZ(CalibIn{icell},x_edge,y_edge,ZIndex(icell));% physical coordinates of the image edge


[40]  386  end


[1078]  387 


 388  %transform the corner coordinates into polar ones


 389  x_edge=x_edgeorigin_xy(1);%shift to the origin of the polar coordinates


 390  y_edge=y_edgeorigin_xy(2);%shift to the origin of the polar coordinates


 391  [theta_edge,r_edge] = cart2pol(x_edge,y_edge);%theta and X are the polar coordinates angle and radius


 392  if (max(theta_edge)min(theta_edge))>pi %if the polar origin is inside the image


 393  r_edge=[0 max(r_edge)];


 394  theta_edge=[pi pi];


[40]  395  end


[1078]  396  rcorner=[rcorner r_edge];


 397  thetacorner=[thetacorner theta_edge];


[40]  398  end


[1078]  399  nbpoint=max(npx.*npy);


 400  Min_r=min(rcorner);


 401  Max_r=max(rcorner);


 402  Min_theta=min(thetacorner)*angle_scale;


 403  Max_theta=max(thetacorner)*angle_scale;


 404  Dr=round_uvmat((Max_rMin_r)/sqrt(nbpoint));


 405  Dtheta=round_uvmat((Max_thetaMin_theta)/sqrt(nbpoint));% get a simple mesh for the rescaled angle


 406  radius=Min_r:Dr:Max_r;% polar coordinates for projections


 407  theta=Min_theta:Dtheta:Max_theta;


 408  [Radius,Theta]=meshgrid(radius,theta/angle_scale);%grid in polar coordinates (angles in radians)


[40]  409  %transform X, Y in cartesian


[1078]  410  [X,Y] = pol2cart(Theta,Radius);% cartesian coordinates associated to the grid in polar coordinates


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


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


 413  radius=radiusradius_offset;


 414  theta=thetaangle_offset*angle_scale;


 415  [np_theta,np_r]=size(Radius);


 416 


 417  for icell=1:length(A)


 418  XIMA=X;


 419  YIMA=Y;


 420  if ~isempty(CalibIn{icell})%transform back to pixel if calibration parameters are introduced


 421  Z=0; %default


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


 423  if ZIndex(icell)==0


 424  ZIndex(icell)=1;


 425  end


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


 427  Z=SliceCoord(3); %to generalize for nonparallel planes


 428  if isfield(CalibIn{icell},'SliceAngle')


 429  norm_plane=angle2normal(CalibIn{icell}.SliceAngle);


 430  Z=Z(norm_plane(1)*(XSliceCoord(1))+norm_plane(2)*(YSliceCoord(2)))/norm_plane(3);


 431  end


 432  end


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


 434  end


 435  Dx=(coord_x{icell}(end)coord_x{icell}(1))/(npx(icell)1);


 436  Dy=(coord_y{icell}(end)coord_y{icell}(1))/(npy(icell)1);


 437  indx_ima=1+round((XIMAcoord_x{icell}(1))/Dx);%indices of the initial matrix close to the points of the new grid


 438  indy_ima=1+round((YIMAcoord_y{icell}(1))/Dy);


 439  Delta_x=1+(XIMAcoord_x{icell}(1))/Dxindx_ima;%


 440  Delta_y=1+(YIMAcoord_y{icell}(1))/Dyindy_ima;


 441  XIMA=reshape(indx_ima,1,[]);%indices reorganized in 'line'


 442  YIMA=reshape(indy_ima,1,[]);%indices reorganized in 'line'


 443  flagin=XIMA>=1 & XIMA<=npx(icell) & YIMA >=1 & YIMA<=npy(icell);%flagin=1 inside the original image


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


[1078]  445  checkuint8=isa(A{icell},'uint8');%check for image input with 8 bits


[1081]  446  checkuint16=isa(A{icell},'uint16');%check for image input with 16 bits


[1078]  447  A{icell}=double(A{icell});


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


[1078]  449  vec_A=reshape(A{icell}(:,:,1),1,[]);%put the original image in line


[164]  450  ind_in=find(flagin);


 451  ind_out=find(~flagin);


[1078]  452  ICOMB=((XIMA1)*npy(icell)+(npy(icell)+1YIMA));


[164]  453  ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A


 454  vec_B(ind_in)=vec_A(ICOMB);


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


[1078]  456  A_out{icell}=reshape(vec_B,np_theta,np_r);%new image in real coordinates


 457  DA_y=circshift(A_out{icell},1,1)A_out{icell};


 458  DA_y(end,:)=0;


 459  DA_x=circshift(A_out{icell},1,2)A_out{icell};


 460  DA_x(:,end)=0;


 461  A_out{icell}=A_out{icell}+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation


[164]  462  else


 463  for icolor=1:siz(3)


[1078]  464  vec_A=reshape(A{icell}(:,:,icolor),1,[]);%put the original image in line


 465  ind_in=find(flagin);


 466  ind_out=find(~flagin);


 467  ICOMB=((XIMA1)*npy(icell)+(npy(icell)+1YIMA));


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


 469  vec_B(ind_in)=vec_A(ICOMB);


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


 471  A_out{icell}(:,:,icolor)=reshape(vec_B,np_theta,np_r);%new image in real coordinates


 472  DA_y=circshift(A_out{icell}(:,:,icolor),1,1)A_out{icell}(:,:,icolor);


 473  DA_y(end,:)=0;


 474  DA_x=circshift(A_out{icell}(:,:,icolor),1,2)A_out{icell}(:,:,icolor);


 475  DA_x(:,end)=0;


 476  A_out{icell}(:,:,icolor)=A_out{icell}(:,:,icolor)+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation


[164]  477  end


 478  end


[1078]  479  if checkuint8


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


 481  elseif checkuint16


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


 483  end


[40]  484  end


 485 

