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

Last change on this file since 497 was 172, checked in by sommeria, 14 years ago

cleaning of documentation

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