phys_polar

transform image coordinates (px) to physical coordinates
 then transform to polar coordinates: 
[DataOut,DataOut_1]=phys_polar(varargin)

 OUTPUT: 
 DataOut: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
 DataOut_1:  second data field (if two fields are in input)

INPUT:
 Data:  structure of input data (like UvData)
 CalibData= structure containing the field .GeometryCalib with calibration parameters
 Data_1:  second input field (not mandatory)
 CalibData_1= calibration parameters for the second field

0001 %transform image coordinates (px) to physical coordinates
0002 % then transform to polar coordinates:
0003 %[DataOut,DataOut_1]=phys_polar(varargin)
0004 %
0005 % OUTPUT:
0006 % DataOut: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
0007 % DataOut_1:  second data field (if two fields are in input)
0008 %
0009 %INPUT:
0010 % Data:  structure of input data (like UvData)
0011 % CalibData= structure containing the field .GeometryCalib with calibration parameters
0012 % Data_1:  second input field (not mandatory)
0013 % CalibData_1= calibration parameters for the second field
0014 
0015 function [DataOut,DataOut_1]=phys_polar(varargin)
0016 Calib{1}=[];
0017 if nargin==2||nargin==4
0018     Data=varargin{1};
0019     DataOut=Data;%default
0020     DataOut_1=[];%default
0021     CalibData=varargin{2};
0022     if isfield(CalibData,'GeometryCalib')
0023         Calib{1}=CalibData.GeometryCalib;
0024     end
0025     Calib{2}=Calib{1};
0026 else
0027     DataOut.Txt='wrong input: need two or four structures';
0028 end
0029 test_1=0;
0030 if nargin==4
0031     test_1=1;
0032     Data_1=varargin{3};
0033     DataOut_1=Data_1;%default
0034     CalibData_1=varargin{4};
0035     if isfield(CalibData_1,'GeometryCalib')
0036         Calib{2}=CalibData_1.GeometryCalib;
0037     end
0038 end
0039 
0040 %parameters for polar coordinates (taken from the calibration data of the first field)
0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0042 origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
0043 if isfield(Calib{1},'PolarCentre') && isnumeric(Calib{1}.PolarCentre)
0044     if isequal(length(Calib{1}.PolarCentre),2);
0045         origin_xy= Calib{1}.PolarCentre;
0046     end
0047 end
0048 radius_offset=0;%reference radius used to offset the radial coordinate r
0049 angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
0050 if isfield(Calib{1},'PolarReferenceRadius') && isnumeric(Calib{1}.PolarReferenceRadius)
0051     radius_offset=Calib{1}.PolarReferenceRadius;
0052 end
0053 if radius_offset > 0
0054     angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
0055 else
0056     angle_scale=180/pi; %polar angle in degrees
0057 end
0058 if isfield(Calib{1},'PolarReferenceAngle') && isnumeric(Calib{1}.PolarReferenceAngle)
0059     angle_offset=Calib{1}.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))
0060 end
0061 % new x coordinate = radius-radius_offset;
0062 % new y coordinate = theta*angle_scale-angle_offset
0063 
0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0065 
0066 iscalar=0;
0067 if  ~isempty(Calib{1})
0068     DataOut=phys_1(Data,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);
0069     %case of images or scalar
0070     if isfield(Data,'A')&isfield(Data,'AX')&~isempty(Data.AX) & isfield(Data,'AY')&...
0071                                            ~isempty(Data.AY)&length(Data.A)>1
0072         iscalar=1;
0073         A{1}=Data.A;
0074     end
0075     %transform of X,Y coordinates for vector fields
0076     if isfield(Data,'ZIndex')&~isempty(Data.ZIndex)
0077         ZIndex=Data.ZIndex;
0078     else
0079         ZIndex=0;
0080     end
0081 end
0082 
0083 if test_1
0084     DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale);
0085     if isfield(Data_1,'A')&isfield(Data_1,'AX')&~isempty(Data_1.AX) & isfield(Data_1,'AY')&...
0086                                        ~isempty(Data_1.AY)&length(Data_1.A)>1
0087           iscalar=iscalar+1;
0088           Calib{iscalar}=Calib{2};
0089           A{iscalar}=Data_1.A;
0090           if isfield(Data_1,'ZIndex')&~isequal(Data_1.ZIndex,ZIndex)
0091               DataOut.Txt='inconsistent plane indexes in the two input fields';
0092           end
0093           if iscalar==1% case for which only the second field is a scalar
0094                [A,AX,AY]=phys_Ima(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);
0095                DataOut_1.A=A{1};
0096                DataOut_1.AX=AX; 
0097                DataOut_1.AY=AY;
0098                return
0099           end
0100     end
0101 end
0102 if iscalar~=0
0103     [A,AX,AY]=phys_Ima(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%
0104     DataOut.A=A{1};
0105     DataOut.AX=AX; 
0106     DataOut.AY=AY;
0107     if iscalar==2
0108         DataOut_1.A=A{2};
0109         DataOut_1.AX=AX; 
0110         DataOut_1.AY=AY;
0111     end
0112 end
0113 
0114 %------------------------------------------------
0115 function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
0116 
0117 DataOut=Data;
0118 DataOut.CoordType='phys'; %put flag for physical coordinates
0119 if isfield(Calib,'CoordUnit')
0120     DataOut.CoordUnit=Calib.CoordUnit;
0121 else
0122     DataOut.CoordUnit='cm'; %default
0123 end
0124 DataOut.TimeUnit='s';
0125 %perform a geometry transform if Calib contains a field .GeometryCalib
0126 if isfield(Data,'CoordType') && isequal(Data.CoordType,'px') && ~isempty(Calib)
0127     if isfield(Data,'CoordUnit')
0128          DataOut=rmfield(DataOut,'CoordUnit');
0129     end
0130     %transform of X,Y coordinates for vector fields
0131     if isfield(Data,'ZIndex')&~isempty(Data.ZIndex)
0132         Z=Data.ZIndex;
0133     else
0134         Z=0;
0135     end
0136     if isfield(Data,'X') &isfield(Data,'Y')&~isempty(Data.X) & ~isempty(Data.Y)
0137         [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical
0138         DataOut.X=DataOut.X-origin_xy(1);%origin of coordinates at the tank center
0139         DataOut.Y=DataOut.Y-origin_xy(2);%origin of coordinates at the tank center
0140         [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta  and X are the polar coordinates angle and radius
0141           %shift and renormalize the polar coordinates
0142         DataOut.X=DataOut.X-radius_offset;%
0143         DataOut.Y=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
0144         %transform velocity field if exists
0145         if isfield(Data,'U')&isfield(Data,'V')&~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'dt') 
0146             if ~isempty(Data.dt)
0147             [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
0148             [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
0149             UX=(XOut_2-XOut_1)/Data.dt;
0150             VY=(YOut_2-YOut_1)/Data.dt;      
0151             %transform u,v into polar coordiantes
0152             DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity
0153             DataOut.V=(-UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);%angular velocity calculated
0154             %shift and renormalize the angular velocity
0155             end
0156         end
0157     end
0158 end
0159 
0160  
0161 %%%%%%%%%%%%%%%%%%%%
0162 function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
0163 xcorner=[];
0164 ycorner=[];
0165 npx=[];
0166 npy=[];
0167 
0168 for icell=1:length(A)
0169     siz=size(A{icell});
0170     npx=[npx siz(2)];
0171     npy=[npy siz(1)];
0172     zphys=0; %default
0173     if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
0174        SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);
0175        zphys=SliceCoord(3); %to generalize for non-parallel planes
0176     end
0177     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
0178     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
0179     [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates
0180     %transform the corner coordinates into polar ones
0181     xcorner_new=xcorner_new-origin_xy(1);%shift to the origin of the polar coordinates
0182     ycorner_new=ycorner_new-origin_xy(2);%shift to the origin of the polar coordinates
0183     [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta  and X are the polar coordinates angle and radius
0184     if (max(theta)-min(theta))>pi   %if the polar origin is inside the image
0185         xcorner_new=[0 max(xcorner_new)];
0186         theta=[-pi pi];
0187     end
0188           %shift and renormalize the polar coordinates
0189     xcorner_new=xcorner_new-radius_offset;%
0190     ycorner_new=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
0191     xcorner=[xcorner xcorner_new];
0192     ycorner=[ycorner ycorner_new];
0193 end
0194 Rangx(1)=min(xcorner);
0195 Rangx(2)=max(xcorner);
0196 Rangy(2)=min(ycorner);
0197 Rangy(1)=max(ycorner);
0198 % test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
0199 npx=max(npx);
0200 npy=max(npy);
0201 x=linspace(Rangx(1),Rangx(2),npx);
0202 y=linspace(Rangy(1),Rangy(2),npy);
0203 [X,Y]=meshgrid(x,y);%grid in physical coordinates
0204 %transform X, Y in cartesian
0205 X=X+radius_offset;%
0206 Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius
0207 [X,Y] = pol2cart(Y,X);
0208 X=X+origin_xy(1);%shift to the origin of the polar coordinates
0209 Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
0210 for icell=1:length(A) 
0211     [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
0212     XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
0213     YIMA=reshape(round(YIMA),1,npx*npy);
0214     flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
0215     vec_A=reshape(A{icell}(:,:,1),1,npx*npy);%put the original image in line
0216     ind_in=find(flagin);
0217     ind_out=find(~flagin);
0218     ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
0219     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
0220     vec_B(ind_in)=vec_A(ICOMB);
0221     vec_B(ind_out)=zeros(size(ind_out));
0222     A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
0223 end
0224 %Rangx=Rangx-radius_offset;
0225 
0226