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

Last change on this file since 167 was 167, checked in by sommeria, 10 years ago

small bug repair in phys and phys_polar: units

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