Changeset 1078 for trunk/src/transform_field/phys_polar.m
 Timestamp:
 Mar 30, 2020, 3:48:19 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/transform_field/phys_polar.m
r1071 r1078 3 3 %%%% Use the general syntax for transform fields %%%% 4 4 % OUTPUT: 5 % Data Out: output field structure5 % Data: output field structure 6 6 % .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components 7 7 % … … 14 14 % .GeometryCalib: substructure of the calibration parameters 15 15 % transform image coordinates (px) to polar physical coordinates 16 %[Data Out,DataOut_1]=phys_polar(varargin)16 %[Data,Data_1]=phys_polar(varargin) 17 17 % 18 18 % OUTPUT: 19 % Data Out: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components20 % Data Out_1: second data field (if two fields are in input)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) 21 21 % 22 22 %INPUT: … … 44 44 %======================================================================= 45 45 46 function Data Out=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)46 function Data=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1) 47 47 % 48 48 49 %% request input parameters 49 50 if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) … … 64 65 end 65 66 answer = inputdlg(prompt,dlg_title,num_lines,def); 66 Data Out.TransformInput.PolarCentre=str2num(answer{1});67 Data Out.TransformInput.PolarReferenceRadius=str2num(answer{2});68 Data Out.TransformInput.PolarReferenceAngle=str2num(answer{3});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 units71 end67 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 72 73 return 73 74 end 74 75 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; 75 85 Calib{1}=[]; 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 76 91 if nargin==2nargin==4 77 DataOut=DataIn;%default 78 DataOut_1=[];%default 79 if isfield(XmlData,'GeometryCalib') 92 if isfield(XmlData,'GeometryCalib') && ~isempty(XmlData.GeometryCalib)&& checkpixel(1) 80 93 Calib{1}=XmlData.GeometryCalib; 81 94 end 82 95 Calib{2}=Calib{1}; 83 96 else 84 Data Out.Txt='wrong input: need two or four structures';85 end 86 test_1=0;97 Data.Txt='wrong input: need two or four structures'; 98 end 99 nbinput=1; 87 100 if nargin==4% case of two input fields 88 test_1=1; 89 DataOut_1=DataIn_1;%default 90 if isfield(XmlData_1,'GeometryCalib') 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) 91 107 Calib{2}=XmlData_1.GeometryCalib; 92 108 end 93 end 94 95 %parameters for polar coordinates (taken from the calibration data of the first field) 109 nbinput=2; 110 end 111 112 %% parameters for polar coordinates (taken from the calibration data of the first field) 96 113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 97 XmlData.PolarReferenceRadius=0;%450;98 XmlData.PolarReferenceAngle=0;%450*pi/2;99 114 origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates 100 115 radius_offset=0;%reference radius used to offset the radial coordinate r 101 116 angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default) 102 angle_scale=180/pi; 117 angle_scale=180/pi; 118 check_degree=1;%angle expressed in degrees by default 103 119 if isfield(XmlData,'TransformInput') 104 120 if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre) 105 if isequal(length(XmlData.TransformInput.PolarCentre),2) ;121 if isequal(length(XmlData.TransformInput.PolarCentre),2) 106 122 origin_xy= XmlData.TransformInput.PolarCentre; 107 123 end … … 112 128 if radius_offset > 0 113 129 angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius 130 check_degree=0; %the output has the same unit asthe input 114 131 else 115 132 angle_scale=180/pi; %polar angle in degrees 133 check_degree=1;%angle expressed in degrees 116 134 end 117 135 if isfield(XmlData.TransformInput,'PolarReferenceAngle') && isnumeric(XmlData.TransformInput.PolarReferenceAngle) … … 121 139 122 140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 123 124 iscalar=0; 125 %transform first field to cartesian phys coordiantes 126 if ~isempty(Calib{1}) 127 DataOut=phys_1(DataIn,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale); 128 %case of images or scalar 129 if isfield(DataIn,'A')&isfield(DataIn,'Coord_x')&~isempty(DataIn.Coord_x) & isfield(DataIn,'Coord_y')&... 130 ~isempty(DataIn.Coord_y)&length(DataIn.A)>1 131 iscalar=1; 132 A{1}=DataIn.A; 133 end 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 157 end 158 % end 134 159 %transform of X,Y coordinates for vector fields 135 if isfield(Data In,'ZIndex')&~isempty(DataIn.ZIndex)136 ZIndex=Data In.ZIndex;160 if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex) 161 ZIndex=DataCell{ifield}.ZIndex; 137 162 else 138 163 ZIndex=0; 139 164 end 140 end 141 142 %transform second field (if exists) to cartesian phys coordiantes 143 if test_1 144 DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale); 145 if isfield(Data_1,'A')&isfield(Data_1,'Coord_x')&~isempty(Data_1.Coord_x) & isfield(Data_1,'Coord_y')&... 146 ~isempty(Data_1.Coord_y)&length(Data_1.A)>1 147 iscalar=iscalar+1; 148 Calib{iscalar}=Calib{2}; 149 A{iscalar}=Data_1.A; 150 if isfield(Data_1,'ZIndex')&~isequal(Data_1.ZIndex,ZIndex) 151 DataOut.Txt='inconsistent plane indexes in the two input fields'; 152 end 153 if iscalar==1% case for which only the second field is a scalar 154 [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale); 155 DataOut_1.A=A{1}; 156 DataOut_1.Coord_x=Coord_x; 157 DataOut_1.Coord_y=Coord_y; 158 return 159 end 160 end 161 end 162 if iscalar~=0 163 [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);% 164 DataOut.A=A{1}; 165 DataOut.Coord_x=Coord_x; 166 DataOut.Coord_y=Coord_y; 167 if iscalar==2 168 DataOut_1.A=A{2}; 169 DataOut_1.Coord_x=Coord_x; 170 DataOut_1.Coord_y=Coord_y; 171 end 172 end 173 174 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 311 end 312 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 332 end 333 end 175 334 176 335 177 336 % 178 function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale) 179 180 DataOut=Data; 181 % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates 182 if isfield(Calib,'SliceCoord') 183 DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane 184 end 185 186 if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib) 187 if isfield(Calib,'CoordUnit') 188 DataOut.CoordUnit=Calib.CoordUnit; 189 else 190 DataOut.CoordUnit='cm'; %default 191 end 192 DataOut.TimeUnit='s'; 193 %transform of X,Y coordinates for vector fields 194 if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 195 Z=Data.ZIndex; 196 else 197 Z=0; 198 end 199 if isfield(Data,'X') &isfield(Data,'Y')&~isempty(Data.X) & ~isempty(Data.Y) 200 [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical 201 DataOut.X=DataOut.Xorigin_xy(1);%origin of coordinates at the tank center 202 DataOut.Y=DataOut.Yorigin_xy(2);%origin of coordinates at the tank center 203 [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta and X are the polar coordinates angle and radius 204 %shift and renormalize the polar coordinates 205 DataOut.X=DataOut.Xradius_offset;%shift the origin of radius, taken as the new X coordinate 206 DataOut.Y=(thetaangle_offset)*angle_scale;% normalized angle: distance along reference radius,taken as the new Y coordinate 207 %transform velocity field if exists 208 if isfield(Data,'U') & isfield(Data,'V') & ~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'Dt') 209 if ~isempty(Data.Dt) 210 [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 211 [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 212 UX=(XOut_2XOut_1)/Data.Dt;% phys velocity u component 213 VY=(YOut_2YOut_1)/Data.Dt; % phys velocity v component 214 %transform u,v into polar coordiantes 215 DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity 216 DataOut.V=(UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);% azimuthal velocity component 217 %shift and renormalize the angular velocity 218 end 219 end 220 %transform of spatial derivatives 221 if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)... 222 && isfield(Data,'Dt') 223 if ~isempty(Data.Dt) 224 % estimate the Jacobian matrix DXpx/DXphys 225 for ip=1:length(Data.X) 226 [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z); 227 [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)0.5,Data.Y(ip),Z); 228 [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z); 229 [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)0.5,Z); 230 %Jacobian matrix DXpphys/DXpx 231 DjXi(1,1)=(Xp1Xm1); 232 DjXi(2,1)=(Yp1Ym1); 233 DjXi(1,2)=(Xp2Xm2); 234 DjXi(2,2)=(Yp2Ym2); 235 DjUi(:,:)=Data.DjUi(ip,:,:); 236 DjUi=(DjXi*DjUi')/DjXi;% =J1*M*J , curvature effects (derivatives of J) neglected 237 DataOut.DjUi(ip,:,:)=DjUi'; 238 end 239 DataOut.DjUi = DataOut.DjUi/Data.Dt; % min(Data.DjUi(:,1,1))=DUDX 240 end 241 end 242 end 243 end 244 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; 354 end 355 if isfield(Data,'Dt')&&~isempty(Data.Dt) 356 Dt=Data.Dt; 357 end 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; 363 end 364 end 245 365 246 366 %%%%%%%%%%%%%%%%%%%% 247 function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale) 248 xcorner=[]; 249 ycorner=[]; 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=[]; 250 372 npx=[]; 251 373 npy=[]; 252 NbPoints=20; % nbre of points used to determine the image edge253 374 for icell=1:length(A) 254 375 siz=size(A{icell}); 255 npx=[npx siz(2)]; 256 npy=[npy siz(1)]; 257 zphys=0; %default 258 if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane 259 if ZIndex==0 260 ZIndex=1; 261 end 262 SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:); 263 zphys=SliceCoord(3); %to generalize for nonparallel planes 264 end 265 % xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordinates of corners 266 % yima=[0.5 0.5 siz(1)0.5 siz(1)0.5]; 267 edge_x=linspace(0.5,siz(1)0.5,NbPoints); 268 edge_y=linspace(0.5,siz(2)0.5,NbPoints); 269 xima=[edge_y (siz(2)0.5)*ones(1,NbPoints) edge_y 0.5*ones(1,NbPoints)];%image coordinates of corners 270 yima=[0.5*ones(1,NbPoints) edge_x (siz(1)0.5)*ones(1,NbPoints) edge_x];%image coordinates of corners 271 [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates 272 %transform the corner coordinates into polar ones 273 xcorner_new=xcorner_neworigin_xy(1);%shift to the origin of the polar coordinates 274 ycorner_new=ycorner_neworigin_xy(2);%shift to the origin of the polar coordinates 275 [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta and X are the polar coordinates angle and radius 276 if (max(theta)min(theta))>pi %if the polar origin is inside the image 277 xcorner_new=[0 max(xcorner_new)]; 278 theta=[pi pi]; 279 end 280 %shift and renormalize the polar coordinates 281 xcorner_new=xcorner_newradius_offset;% 282 ycorner_new=theta*angle_scaleangle_offset;% normalized angle: distance along reference radius 283 xcorner=[xcorner xcorner_new]; 284 ycorner=[ycorner ycorner_new]; 285 end 286 Rangx(1)=min(xcorner); 287 Rangx(2)=max(xcorner); 288 Rangy(2)=min(ycorner); 289 Rangy(1)=max(ycorner); 290 % test_multi=(max(npx)~=min(npx))  (max(npy)~=min(npy)); 291 npx=max(npx); 292 npy=max(npy); 293 x=linspace(Rangx(1),Rangx(2),npx); 294 y=linspace(Rangy(1),Rangy(2),npy); 295 [X,Y]=meshgrid(x,y);%grid in physical coordinates 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 386 end 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]; 395 end 396 rcorner=[rcorner r_edge]; 397 thetacorner=[thetacorner theta_edge]; 398 end 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) 296 409 %transform X, Y in cartesian 297 X=X+radius_offset;% 298 Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius 299 [X,Y] = pol2cart(Y,X); 300 X=X+origin_xy(1);%shift to the origin of the polar coordinates 301 Y=Y+origin_xy(2);%shift to the origin of the polar coordinates 302 for icell=1:length(A) 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 303 444 siz=size(A{icell}); 304 [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid 305 XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line' 306 YIMA=reshape(round(YIMA),1,npx*npy); 307 flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image 445 checkuint8=isa(A{icell},'uint8');%check for image input with 8 bits 446 checkuint16=isa(A{icell},'uint8');%check for image input with 16 bits 447 A{icell}=double(A{icell}); 308 448 if numel(siz)==2 %(B/W images) 309 vec_A=reshape(A{icell}(:,:,1),1, npx*npy);%put the original image in line449 vec_A=reshape(A{icell}(:,:,1),1,[]);%put the original image in line 310 450 ind_in=find(flagin); 311 451 ind_out=find(~flagin); 312 ICOMB=((XIMA1)*npy +(npy+1YIMA));452 ICOMB=((XIMA1)*npy(icell)+(npy(icell)+1YIMA)); 313 453 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 314 454 vec_B(ind_in)=vec_A(ICOMB); 315 455 vec_B(ind_out)=zeros(size(ind_out)); 316 A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates 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 317 462 else 318 463 for icolor=1:siz(3) 319 vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line 320 ind_in=find(flagin); 321 ind_out=find(~flagin); 322 ICOMB=((XIMA1)*npy+(npy+1YIMA)); 323 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 324 vec_B(ind_in)=vec_A(ICOMB); 325 vec_B(ind_out)=zeros(size(ind_out)); 326 A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates 327 end 328 end 329 end 330 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 477 end 478 end 479 if checkuint8 480 A_out{icell}=uint8(A_out{icell}); 481 elseif checkuint16 482 A_out{icell}=uint16(A_out{icell}); 483 end 484 end 485
Note: See TracChangeset
for help on using the changeset viewer.