[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 2008-2020, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
|
---|
[810] | 30 | % http://www.legi.grenoble-inp.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==2||nargin==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{nbvar-1}.Role='coord_x';
|
---|
| 178 | check_unit=1;
|
---|
| 179 | if isfield(DataCell{ifield},'CoordUnit')
|
---|
| 180 | Data=rmfield(Data,'CoordUnit');
|
---|
| 181 | Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
|
---|
| 182 | elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
|
---|
| 183 | Data.VarAttribute{nbvar-1}.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{nbvar-1}.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(X-origin_xy(1),Y-origin_xy(2));
|
---|
| 206 | Data.(radius_name)=Radius-radius_offset;
|
---|
| 207 | Data.(theta_name)=Theta*angle_scale-angle_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{nbvar-1}.Role='coord_x';
|
---|
| 255 | Data.VarAttribute{nbvar}.Role='coord_y';
|
---|
| 256 | check_unit=1;
|
---|
| 257 | if isfield(DataCell{ifield},'CoordUnit')
|
---|
| 258 | Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
|
---|
| 259 | elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
|
---|
| 260 | Data.VarAttribute{nbvar-1}.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{nbvar-1}.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_scale-angle_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.X-Data.U/2,Data.Y-Data.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_2-XOut_1)/Dt;
|
---|
| 362 | V=(YOut_2-YOut_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_edge-origin_xy(1);%shift to the origin of the polar coordinates
|
---|
| 390 | y_edge=y_edge-origin_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_r-Min_r)/sqrt(nbpoint));
|
---|
| 405 | Dtheta=round_uvmat((Max_theta-Min_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=radius-radius_offset;
|
---|
| 414 | theta=theta-angle_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 non-parallel planes
|
---|
| 428 | if isfield(CalibIn{icell},'SliceAngle')
|
---|
| 429 | norm_plane=angle2normal(CalibIn{icell}.SliceAngle);
|
---|
| 430 | Z=Z-(norm_plane(1)*(X-SliceCoord(1))+norm_plane(2)*(Y-SliceCoord(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((XIMA-coord_x{icell}(1))/Dx);%indices of the initial matrix close to the points of the new grid
|
---|
| 438 | indy_ima=1+round((YIMA-coord_y{icell}(1))/Dy);
|
---|
| 439 | Delta_x=1+(XIMA-coord_x{icell}(1))/Dx-indx_ima;%
|
---|
| 440 | Delta_y=1+(YIMA-coord_y{icell}(1))/Dy-indy_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=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
|
---|
[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=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
|
---|
| 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 |
|
---|