Home > . > phys_polar.m

# phys_polar

## PURPOSE

transform image coordinates (px) to physical coordinates

## SYNOPSIS

function [DataOut,DataOut_1]=phys_polar(varargin)

## DESCRIPTION

    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

## CROSS-REFERENCE INFORMATION

This function calls:

- phys_XYZ 'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
- px_XYZ 'px_XYZ': transform phys coordinates to image coordinates (px)

This function is called by:

## SUBFUNCTIONS

- function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
- function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)

## SOURCE CODE

    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