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


[1126]  29  % Copyright 20082024, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[810]  30  % http://www.legi.grenobleinp.fr


[1127]  31  % Joel.Sommeria  Joel.Sommeria (A) univgrenoblealpes.fr


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


[1084]  51  prompt = {'origin [x y] of polar coordinates';'reference radius';'reference angle(degrees)';'angle direction and switch x y(+/)'};


[933]  52  dlg_title = 'set the parameters for the polar coordinates';


 53  num_lines= 2;


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


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


[1084]  65  if isfield(XmlData.TransformInput,'PolarAngleDirection')


 66  def{4}=XmlData.TransformInput.PolarAngleDirection;


 67  end


[933]  68  end


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


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


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


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


[1084]  73  Data.TransformInput.PolarAngleDirection=answer{4};


[933]  74  return


 75  end


 76 


[1078]  77  %% default outputs


 78  Data=DataIn; %default output


 79  if isfield(Data,'CoordUnit')


[1112]  80  Data=rmfield(Data,'CoordUnit');


[1078]  81  end


 82  Data.ListVarName = {};


 83  Data.VarDimName={};


 84  Data.VarAttribute={};


 85  DataCell{1}=DataIn;


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


[1112]  87  Slice{1}=[];


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


 89  checkpixel(1)=0;


[1112]  90  if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel')


[1078]  91  checkpixel(1)=1;


 92  end


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


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


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


[40]  96  end


[1112]  97  Slice{1}=Calib{1};


 98  if isfield(XmlData,'Slice')


 99  Slice{1}=XmlData.Slice;


 100  end


[40]  101  Calib{2}=Calib{1};


[1112]  102  Slice{2}=Slice{1};


[40]  103  else


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


[40]  105  end


[1078]  106  nbinput=1;


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


[1078]  108  checkpixel(2)=0;


[1112]  109  if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel')


 110  checkpixel(2)=1;


 111  end


[1078]  112  DataCell{2}=DataIn_1;%default


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


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


[40]  115  end


[1112]  116  if isfield(XmlData_1,'Slice')


 117  Slice{2}=XmlData_1.Slice;


 118  end


[1078]  119  nbinput=2;


[40]  120  end


 121 


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


[40]  123  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


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


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


[1078]  127  angle_scale=180/pi;


[1104]  128  check_reverse=false;


[1078]  129  check_degree=1;%angle expressed in degrees by default


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


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


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


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


 134  end


[40]  135  end


[1094]  136  if isfield(XmlData.TransformInput,'PolarReferenceRadius') && ~isempty(XmlData.TransformInput.PolarReferenceRadius)


[933]  137  radius_offset=XmlData.TransformInput.PolarReferenceRadius;


 138  end


 139  if radius_offset > 0


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


[1094]  141  check_degree=0; %the output has the same unit as the input


[933]  142  else


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


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


[933]  145  end


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


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


 148  end


[1084]  149  check_reverse=isfield(XmlData.TransformInput,'PolarAngleDirection')&& strcmp(XmlData.TransformInput.PolarAngleDirection,'');


[40]  150  end


 151 


 152  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[1078]  153  %% get fields


[1084]  154 


[1078]  155  nbvar=0;%counter for the number of output variables


[1094]  156  nbcoord=0;%counter for the number of variablecheck_degrees for radial coordiantes (case of multiple field inputs)


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


 158  nbscattered=0;%counter of scattered fields


 159  radius_name='radius';


 160  theta_name='theta';


 161  U_r_name='U_r';


 162  U_theta_name='U_theta';


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


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


 165  if ~isempty(errormsg)


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


 167  return


[40]  168  end


 169  %transform of X,Y coordinates for vector fields


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


 171  ZIndex=DataCell{ifield}.ZIndex;


[40]  172  else


 173  ZIndex=0;


 174  end


[1084]  175  check_scalar=zeros(1,numel(CellInfo));


 176  check_vector=zeros(1,numel(CellInfo));


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


 178  if NbDim(icell)==2


 179  % case of input field with scattered coordinates


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


 181  nbscattered=nbscattered+1;


 182  nbcoord=nbcoord+1;


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


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


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


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


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


 188  nbvar=nbvar+2;


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


 190  check_unit=1;


[1094]  191  %unit of output field


 192  if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


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


 194  elseif isfield(DataCell{ifield},'CoordUnit')


 195  radius_unit=DataCell{ifield}.CoordUnit;


[1078]  196  else


[1094]  197  radius_unit='';


[1078]  198  end


[1094]  199  Data.VarAttribute{nbvar1}.units=radius_unit;


[1078]  200  if check_degree


[1094]  201  Data.VarAttribute{nbvar}.units='degree';


 202  else %case of a reference radius


 203  Data.VarAttribute{nbvar}.units=radius_unit;


 204  Data.CoordUnit=radius_unit;


[1078]  205  end


[1094]  206  % if isfield(DataCell{ifield},'CoordUnit')


 207  % Data=rmfield(Data,'CoordUnit');


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


 209  % elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


 210  % 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


 211  % else


 212  % check_unit=0;


 213  % end


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


 215  % if check_degree


 216  % Data.VarAttribute{nbvar}.units='degree';


 217  % elseif check_unit


 218  % Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar1}.units;


 219  % end


[1078]  220 


 221  %transform u,v into polar coordinates


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


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


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


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


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


 227  if ~isempty(Calib{ifield})


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


[1112]  229  phys_XYUV(DataCell{ifield},Calib{ifield},Slice{ifield},ZIndex);


[1078]  230  end


 231  end


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


 233  Data.(radius_name)=Radiusradius_offset;


 234  Data.(theta_name)=Theta*angle_scaleangle_offset;


 235  if Z~=0


 236  Data.Z=Z;


 237  nbvar=nbvar+1;


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


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


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


 241  end


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


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


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


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


 246  nbvar=nbvar+1;


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


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


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


 250  end


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


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


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


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


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


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


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


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


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


 260  nbvar=nbvar+2;


 261  end


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


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


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


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


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


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


 268  nbvar=nbvar+1;


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


 270  end


 271 


 272  %caseof input fields on gridded coordinates (matrix)


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


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


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


[1094]  276  radius_name = rename_indexing(radius_name,Data.ListVarName);% add an index to the name, or increment an existing index,


 277  theta_name = rename_indexing(theta_name,Data.ListVarName);% if the proposed Name already exists in the list


 278  Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];%add polar coordinates to the list of variables


[1078]  279  Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}];


 280  nbvar=nbvar+2;


[1084]  281  if check_reverse


[1094]  282  Data.VarAttribute{nbvar1}.Role='coord_y';


 283  Data.VarAttribute{nbvar}.Role='coord_x';


[1084]  284  else


[1094]  285  Data.VarAttribute{nbvar1}.Role='coord_x';


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


[1084]  287  end


[1078]  288  check_unit=1;


[1094]  289 


 290  if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')


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


 292  elseif isfield(DataCell{ifield},'CoordUnit')


 293  Data.VarAttribute{nbvar1}.units=DataCell{ifield}.CoordUnit;%radius in coord units


[1078]  294  else


 295  check_unit=0;


 296  end


 297  if check_degree


[1094]  298  Data.VarAttribute{nbvar}.units='degree';%angle in degree


[1078]  299  elseif check_unit


[1094]  300  Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar1}.units;% angle in coord unit (normalised by reference radiuss)


[1078]  301  end


 302  end


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


 304  nbgrid=nbgrid+1;


 305  nbvar=nbvar+1;


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


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


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


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


 310  check_scalar(nbgrid)=1;


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


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


 313  ZInd(nbgrid)=ZIndex;


 314  Calib_new{nbgrid}=Calib{ifield};


[1112]  315  Slice_new{nbgrid}=Slice{ifield};


[1078]  316  end


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


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


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


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


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


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


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


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


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


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


 327  check_vector(nbgrid+1)=1;


 328  check_vector(nbgrid+2)=1;


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


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


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


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


 333  ZInd(nbgrid+1)=ZIndex;


 334  ZInd(nbgrid+2)=ZIndex;


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


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


[1112]  337  Slice_new{nbgrid+1}=Calib{ifield};


 338  Slice_new{nbgrid+2}=Calib{ifield};


[1078]  339  nbgrid=nbgrid+2;


 340  nbvar=nbvar+2;


 341  end


 342  end


 343  end


 344  end


[40]  345  end


[567]  346 


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


 348  if nbgrid~=0


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


[1078]  350  for icell=1:numel(A)


[1084]  351  if icell<=numel(A)1 && check_vector(icell)==1 && check_vector(icell+1)==1 %transform u,v into polar coordinates


[1078]  352  theta=Data.theta/angle_scaleangle_offset;


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


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


[1084]  355  U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);


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


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


 358  if check_reverse


 359  Data.(U_theta_name)=(Data.(U_theta_name))';


 360  Data.(U_r_name)=Data.(U_r_name)';


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


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


 363  else


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


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


 366  end


[1078]  367  elseif ~check_vector(icell)% for scalar fields


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


[1084]  369  Data.ListVarName=[Data.ListVarName FieldName(icell)];


 370  if check_reverse


 371  Data.(FieldName{icell})=A{icell}';


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


 373  else


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


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


 376  end


[1078]  377  end


[40]  378  end


 379  end


[1084]  380  if check_reverse


 381  Data.(theta_name)=Data.(theta_name);


 382  end


[40]  383 


[161]  384 


[40]  385  %


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


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


[1078]  388  %


 389  %% set default output


 390  %DataOut=Data;%default


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


 392  X=[];%default output


 393  Y=[];


 394  Z=0;


 395  U=[];


 396  V=[];


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


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


[1112]  399  [X,Y,Z]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);


[1078]  400  Dt=1; %default


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


 402  Dt=Data.dt;


[40]  403  end


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


 405  Dt=Data.Dt;


[40]  406  end


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


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


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


[1078]  410  U=(XOut_2XOut_1)/Dt;


 411  V=(YOut_2YOut_1)/Dt;


[40]  412  end


 413  end


 414 


 415  %%%%%%%%%%%%%%%%%%%%


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


 417  % transform to phys coordiantes if requested by calibration input


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


[1078]  419  rcorner=[];


 420  thetacorner=[];


[40]  421  npx=[];


 422  npy=[];


 423  for icell=1:length(A)


 424  siz=size(A{icell});


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


 426  npy(icell)=siz(1);


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


 428  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)


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


 430  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)


 431 


 432  % transform edges into phys coordinates if requested


 433  if ~isempty(CalibIn{icell})


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


[40]  435  end


[1078]  436 


 437  %transform the corner coordinates into polar ones


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


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


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


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


 442  r_edge=[0 max(r_edge)];


 443  theta_edge=[pi pi];


[40]  444  end


[1078]  445  rcorner=[rcorner r_edge];


 446  thetacorner=[thetacorner theta_edge];


[40]  447  end


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


 449  Min_r=min(rcorner);


 450  Max_r=max(rcorner);


 451  Min_theta=min(thetacorner)*angle_scale;


 452  Max_theta=max(thetacorner)*angle_scale;


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


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


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


 456  theta=Min_theta:Dtheta:Max_theta;


[1084]  457  %theta=Max_theta:Dtheta:Min_theta;


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


[40]  459  %transform X, Y in cartesian


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


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


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


 463  radius=radiusradius_offset;


 464  theta=thetaangle_offset*angle_scale;


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


 466 


 467  for icell=1:length(A)


 468  XIMA=X;


 469  YIMA=Y;


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


 471  Z=0; %default


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


 473  if ZIndex(icell)==0


 474  ZIndex(icell)=1;


 475  end


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


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


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


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


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


 481  end


 482  end


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


 484  end


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


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


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


[1084]  488  %indy_ima=1+round((YIMAcoord_y{icell}(1))/Dy);


 489  indy_ima=1+round((coord_y{icell}(end)YIMA)/Dy);


 490  Delta_x=1+(XIMAcoord_x{icell}(1))/Dxindx_ima;%error in the index discretisation


 491  Delta_y=1+(coord_y{icell}(end)YIMA)/Dyindy_ima;


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


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


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


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


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


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


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


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


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


[164]  501  ind_in=find(flagin);


 502  ind_out=find(~flagin);


[1084]  503  ICOMB=((XIMA1)*npy(icell)+(npy(icell)+1YIMA));% indices in vec_A


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


 505  vec_B(ind_in)=vec_A(ICOMB);


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


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


[1084]  508  DA_y=circshift(A_out{icell},1,1)A_out{icell};% derivative


[1078]  509  DA_y(end,:)=0;


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


 511  DA_x(:,end)=0;


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


[164]  513  else


 514  for icolor=1:siz(3)


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


 516  ind_in=find(flagin);


 517  ind_out=find(~flagin);


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


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


 520  vec_B(ind_in)=vec_A(ICOMB);


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


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


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


 524  DA_y(end,:)=0;


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


 526  DA_x(:,end)=0;


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


[164]  528  end


 529  end


[1078]  530  if checkuint8


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


 532  elseif checkuint16


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


 534  end


[40]  535  end


 536 

