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

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