source: trunk/src/transform_field/phys_polar.m @ 686

Last change on this file since 686 was 658, checked in by sommeria, 11 years ago

a few bugs corrected

File size: 11.8 KB
Line 
1%'phys_polar': transforms image (Unit='pixel') to polar (phys) coordinates using geometric calibration parameters
2
3%------------------------------------------------------------------------
4%%%%  Use the general syntax for transform fields %%%%
5% OUTPUT:
6% DataOut:   output field structure
7%      .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
8
9%INPUT:
10% DataIn:  first input field structure
11% XmlData: first input parameter structure,
12%        .GeometryCalib: substructure of the calibration parameters
13% DataIn_1: optional second input field structure
14% XmlData_1: optional second input parameter structure
15%         .GeometryCalib: substructure of the calibration parameters
16% transform image coordinates (px) to polar physical coordinates
17%[DataOut,DataOut_1]=phys_polar(varargin)
18%
19% OUTPUT:
20% DataOut: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
21% DataOut_1:  second data field (if two fields are in input)
22%
23%INPUT:
24% Data:  structure of input data (like UvData)
25% XmlData= structure containing the field .GeometryCalib with calibration parameters
26% Data_1:  second input field (not mandatory)
27% XmlData_1= calibration parameters for the second field
28%------------------------------------------------------------------------
29function DataOut=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)
30%------------------------------------------------------------------------
31Calib{1}=[];
32if nargin==2||nargin==4
33    DataOut=DataIn;%default
34    DataOut_1=[];%default
35    if isfield(XmlData,'GeometryCalib')
36        Calib{1}=XmlData.GeometryCalib;
37    end
38    Calib{2}=Calib{1};
39else
40    DataOut.Txt='wrong input: need two or four structures';
41end
42test_1=0;
43if nargin==4% case of two input fields
44    test_1=1;
45    DataOut_1=DataIn_1;%default
46    if isfield(XmlData_1,'GeometryCalib')
47        Calib{2}=XmlData_1.GeometryCalib;
48    end
49end
50
51%parameters for polar coordinates (taken from the calibration data of the first field)
52%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
54if isfield(XmlData,'PolarCentre') && isnumeric(XmlData.PolarCentre)
55    if isequal(length(XmlData.PolarCentre),2);
56        origin_xy= XmlData.PolarCentre;
57    end
58end
59radius_offset=0;%reference radius used to offset the radial coordinate r
60angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
61if isfield(XmlData,'PolarReferenceRadius') && isnumeric(XmlData.PolarReferenceRadius)
62    radius_offset=XmlData.PolarReferenceRadius;
63end
64if radius_offset > 0
65    angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
66else
67    angle_scale=180/pi; %polar angle in degrees
68end
69if isfield(XmlData,'PolarReferenceAngle') && isnumeric(XmlData.PolarReferenceAngle)
70    angle_offset=XmlData.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))
71end
72% new x coordinate = radius-radius_offset;
73% new y coordinate = theta*angle_scale-angle_offset
74
75%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76
77iscalar=0;
78%transform first field to cartesian phys coordiantes
79if  ~isempty(Calib{1})
80    DataOut=phys_1(Data,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);
81    %case of images or scalar
82    if isfield(Data,'A')&isfield(Data,'AX')&~isempty(Data.AX) & isfield(Data,'AY')&...
83                                           ~isempty(Data.AY)&length(Data.A)>1
84        iscalar=1;
85        A{1}=Data.A;
86    end
87    %transform of X,Y coordinates for vector fields
88    if isfield(Data,'ZIndex')&~isempty(Data.ZIndex)
89        ZIndex=Data.ZIndex;
90    else
91        ZIndex=0;
92    end
93end
94
95%transform second field (if exists) to cartesian phys coordiantes
96if test_1
97    DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale);
98    if isfield(Data_1,'A')&isfield(Data_1,'AX')&~isempty(Data_1.AX) & isfield(Data_1,'AY')&...
99                                       ~isempty(Data_1.AY)&length(Data_1.A)>1
100          iscalar=iscalar+1;
101          Calib{iscalar}=Calib{2};
102          A{iscalar}=Data_1.A;
103          if isfield(Data_1,'ZIndex')&~isequal(Data_1.ZIndex,ZIndex)
104              DataOut.Txt='inconsistent plane indexes in the two input fields';
105          end
106          if iscalar==1% case for which only the second field is a scalar
107               [A,AX,AY]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);
108               DataOut_1.A=A{1};
109               DataOut_1.AX=AX;
110               DataOut_1.AY=AY;
111               return
112          end
113    end
114end
115if iscalar~=0
116    [A,AX,AY]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%
117    DataOut.A=A{1};
118    DataOut.AX=AX;
119    DataOut.AY=AY;
120    if iscalar==2
121        DataOut_1.A=A{2};
122        DataOut_1.AX=AX;
123        DataOut_1.AY=AY;
124    end
125end
126
127
128
129
130%------------------------------------------------
131function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
132
133DataOut=Data;
134% DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
135if isfield(Calib,'SliceCoord')
136    DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane
137end
138
139if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
140    if isfield(Calib,'CoordUnit')
141        DataOut.CoordUnit=Calib.CoordUnit;
142    else
143        DataOut.CoordUnit='cm'; %default
144    end
145    DataOut.TimeUnit='s';
146    %transform of X,Y coordinates for vector fields
147    if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
148        Z=Data.ZIndex;
149    else
150        Z=0;
151    end
152    if isfield(Data,'X') &isfield(Data,'Y')&~isempty(Data.X) & ~isempty(Data.Y)
153        [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical
154        DataOut.X=DataOut.X-origin_xy(1);%origin of coordinates at the tank center
155        DataOut.Y=DataOut.Y-origin_xy(2);%origin of coordinates at the tank center
156        [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta  and X are the polar coordinates angle and radius
157          %shift and renormalize the polar coordinates
158        DataOut.X=DataOut.X-radius_offset;%
159        DataOut.Y=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
160        %transform velocity field if exists
161        if isfield(Data,'U')&isfield(Data,'V')&~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'dt')
162            if ~isempty(Data.dt)
163            [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
164            [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
165            UX=(XOut_2-XOut_1)/Data.dt;
166            VY=(YOut_2-YOut_1)/Data.dt;     
167            %transform u,v into polar coordiantes
168            DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity
169            DataOut.V=(-UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);%angular velocity calculated
170            %shift and renormalize the angular velocity
171            end
172        end
173        %transform of spatial derivatives
174        if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
175                && isfield(Data,'dt')
176            if ~isempty(Data.dt)
177                % estimate the Jacobian matrix DXpx/DXphys
178                for ip=1:length(Data.X)
179                    [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
180                    [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
181                    [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
182                    [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
183                    %Jacobian matrix DXpphys/DXpx
184                    DjXi(1,1)=(Xp1-Xm1);
185                    DjXi(2,1)=(Yp1-Ym1);
186                    DjXi(1,2)=(Xp2-Xm2);
187                    DjXi(2,2)=(Yp2-Ym2);
188                    DjUi(:,:)=Data.DjUi(ip,:,:);
189                    DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
190                    DataOut.DjUi(ip,:,:)=DjUi';
191                end
192                DataOut.DjUi =  DataOut.DjUi/Data.dt;   %     min(Data.DjUi(:,1,1))=DUDX
193            end
194        end
195    end
196end
197
198
199%%%%%%%%%%%%%%%%%%%%
200function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
201xcorner=[];
202ycorner=[];
203npx=[];
204npy=[];
205for icell=1:length(A)
206    siz=size(A{icell});
207    npx=[npx siz(2)];
208    npy=[npy siz(1)];
209    zphys=0; %default
210    if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
211       SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);
212       zphys=SliceCoord(3); %to generalize for non-parallel planes
213    end
214    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
215    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
216    [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates
217    %transform the corner coordinates into polar ones   
218    xcorner_new=xcorner_new-origin_xy(1);%shift to the origin of the polar coordinates
219    ycorner_new=ycorner_new-origin_xy(2);%shift to the origin of the polar coordinates       
220    [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta  and X are the polar coordinates angle and radius
221    if (max(theta)-min(theta))>pi   %if the polar origin is inside the image
222        xcorner_new=[0 max(xcorner_new)];
223        theta=[-pi pi];
224    end
225          %shift and renormalize the polar coordinates
226    xcorner_new=xcorner_new-radius_offset;%
227    ycorner_new=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
228    xcorner=[xcorner xcorner_new];
229    ycorner=[ycorner ycorner_new];
230end
231Rangx(1)=min(xcorner);
232Rangx(2)=max(xcorner);
233Rangy(2)=min(ycorner);
234Rangy(1)=max(ycorner);
235% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
236npx=max(npx);
237npy=max(npy);
238x=linspace(Rangx(1),Rangx(2),npx);
239y=linspace(Rangy(1),Rangy(2),npy);
240[X,Y]=meshgrid(x,y);%grid in physical coordinates
241%transform X, Y in cartesian
242X=X+radius_offset;%
243Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius
244[X,Y] = pol2cart(Y,X);
245X=X+origin_xy(1);%shift to the origin of the polar coordinates
246Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
247for icell=1:length(A)
248    siz=size(A{icell});
249    [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
250    XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
251    YIMA=reshape(round(YIMA),1,npx*npy);
252    flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
253    if numel(siz)==2 %(B/W images)
254        vec_A=reshape(A{icell}(:,:,1),1,npx*npy);%put the original image in line
255        ind_in=find(flagin);
256        ind_out=find(~flagin);
257        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
258        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
259        vec_B(ind_in)=vec_A(ICOMB);
260        vec_B(ind_out)=zeros(size(ind_out));
261        A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
262    else
263        for icolor=1:siz(3)
264                vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
265                ind_in=find(flagin);
266                ind_out=find(~flagin);
267                ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
268                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
269                vec_B(ind_in)=vec_A(ICOMB);
270                vec_B(ind_out)=zeros(size(ind_out));
271                A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
272        end
273    end
274end
275
Note: See TracBrowser for help on using the repository browser.