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

Last change on this file since 1061 was 1061, checked in by g7moreau, 5 years ago
  • Update copyright to year 2019
File size: 15.0 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-2019, LEGI UMR 5519 / CNRS UGA 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%------------------------------------------------------------------------
48%% request input parameters
49if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
50    prompt = {'origin [x y] of polar coordinates';'reference radius';'reference angle(degrees)'};
51    dlg_title = 'set the parameters for the polar coordinates';
52    num_lines= 2;
53    def     = { '[0 0]';'0';'0'};
54    if isfield(XmlData,'TransformInput')
55        if isfield(XmlData.TransformInput,'PolarCentre')
56            def{1}=num2str(XmlData.TransformInput.PolarCentre);
57        end
58        if isfield(XmlData.TransformInput,'PolarReferenceRadius')
59            def{2}=num2str(XmlData.TransformInput.PolarReferenceRadius);
60        end
61        if isfield(XmlData.TransformInput,'PolarReferenceAngle')
62            def{3}=num2str(XmlData.TransformInput.PolarReferenceAngle);
63        end
64    end
65    answer = inputdlg(prompt,dlg_title,num_lines,def);
66    DataOut.TransformInput.PolarCentre=str2num(answer{1});
67    DataOut.TransformInput.PolarReferenceRadius=str2num(answer{2});
68    DataOut.TransformInput.PolarReferenceAngle=str2num(answer{3});
69    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
70        DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
71    end
72    return
73end
74
75Calib{1}=[];
76if nargin==2||nargin==4
77    DataOut=DataIn;%default
78    DataOut_1=[];%default
79    if isfield(XmlData,'GeometryCalib')
80        Calib{1}=XmlData.GeometryCalib;
81    end
82    Calib{2}=Calib{1};
83else
84    DataOut.Txt='wrong input: need two or four structures';
85end
86test_1=0;
87if nargin==4% case of two input fields
88    test_1=1;
89    DataOut_1=DataIn_1;%default
90    if isfield(XmlData_1,'GeometryCalib')
91        Calib{2}=XmlData_1.GeometryCalib;
92    end
93end
94
95%parameters for polar coordinates (taken from the calibration data of the first field)
96%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
97XmlData.PolarReferenceRadius=0;%450;
98XmlData.PolarReferenceAngle=0;%450*pi/2;
99origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
100radius_offset=0;%reference radius used to offset the radial coordinate r
101angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
102angle_scale=180/pi;
103if isfield(XmlData,'TransformInput')
104    if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre)
105        if isequal(length(XmlData.TransformInput.PolarCentre),2);
106            origin_xy= XmlData.TransformInput.PolarCentre;
107        end
108    end
109    if isfield(XmlData.TransformInput,'PolarReferenceRadius') && isnumeric(XmlData.TransformInput.PolarReferenceRadius)
110        radius_offset=XmlData.TransformInput.PolarReferenceRadius;
111    end
112    if radius_offset > 0
113        angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
114    else
115        angle_scale=180/pi; %polar angle in degrees
116    end
117    if isfield(XmlData.TransformInput,'PolarReferenceAngle') && isnumeric(XmlData.TransformInput.PolarReferenceAngle)
118        angle_offset=(pi/180)*XmlData.TransformInput.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))
119    end
120end
121
122%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123
124iscalar=0;
125%transform first field to cartesian phys coordiantes
126if  ~isempty(Calib{1})
127    DataOut=phys_1(DataIn,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);
128    %case of images or scalar
129    if isfield(DataIn,'A')&isfield(DataIn,'Coord_x')&~isempty(DataIn.Coord_x) & isfield(DataIn,'Coord_y')&...
130                                           ~isempty(DataIn.Coord_y)&length(DataIn.A)>1
131        iscalar=1;
132        A{1}=DataIn.A;
133    end
134    %transform of X,Y coordinates for vector fields
135    if isfield(DataIn,'ZIndex')&~isempty(DataIn.ZIndex)
136        ZIndex=DataIn.ZIndex;
137    else
138        ZIndex=0;
139    end
140end
141
142%transform second field (if exists) to cartesian phys coordiantes
143if test_1
144    DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale);
145    if isfield(Data_1,'A')&isfield(Data_1,'Coord_x')&~isempty(Data_1.Coord_x) & isfield(Data_1,'Coord_y')&...
146                                       ~isempty(Data_1.Coord_y)&length(Data_1.A)>1
147          iscalar=iscalar+1;
148          Calib{iscalar}=Calib{2};
149          A{iscalar}=Data_1.A;
150          if isfield(Data_1,'ZIndex')&~isequal(Data_1.ZIndex,ZIndex)
151              DataOut.Txt='inconsistent plane indexes in the two input fields';
152          end
153          if iscalar==1% case for which only the second field is a scalar
154               [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);
155               DataOut_1.A=A{1};
156               DataOut_1.Coord_x=Coord_x;
157               DataOut_1.Coord_y=Coord_y;
158               return
159          end
160    end
161end
162if iscalar~=0
163    [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%
164    DataOut.A=A{1};
165    DataOut.Coord_x=Coord_x;
166    DataOut.Coord_y=Coord_y;
167    if iscalar==2
168        DataOut_1.A=A{2};
169        DataOut_1.Coord_x=Coord_x;
170        DataOut_1.Coord_y=Coord_y;
171    end
172end
173
174
175
176
177%------------------------------------------------
178function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
179
180DataOut=Data;
181% DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
182if isfield(Calib,'SliceCoord')
183    DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane
184end
185
186if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
187    if isfield(Calib,'CoordUnit')
188        DataOut.CoordUnit=Calib.CoordUnit;
189    else
190        DataOut.CoordUnit='cm'; %default
191    end
192    DataOut.TimeUnit='s';
193    %transform of X,Y coordinates for vector fields
194    if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
195        Z=Data.ZIndex;
196    else
197        Z=0;
198    end
199    if isfield(Data,'X') &isfield(Data,'Y')&~isempty(Data.X) & ~isempty(Data.Y)
200        [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical
201        DataOut.X=DataOut.X-origin_xy(1);%origin of coordinates at the tank center
202        DataOut.Y=DataOut.Y-origin_xy(2);%origin of coordinates at the tank center
203        [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta  and X are the polar coordinates angle and radius
204          %shift and renormalize the polar coordinates
205        DataOut.X=DataOut.X-radius_offset;%shift the origin of radius, taken as the new X coordinate
206        DataOut.Y=(theta-angle_offset)*angle_scale;% normalized angle: distance along reference radius,taken as the new Y coordinate
207        %transform velocity field if exists
208        if isfield(Data,'U') & isfield(Data,'V') & ~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'Dt')
209            if ~isempty(Data.Dt)
210            [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);% X,Y positions of the vector origin in phys
211            [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);% X,Y positions of the vector end in phys
212            UX=(XOut_2-XOut_1)/Data.Dt;% phys velocity u component
213            VY=(YOut_2-YOut_1)/Data.Dt; % phys velocity v component     
214            %transform u,v into polar coordiantes
215            DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity
216            DataOut.V=(-UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);% azimuthal velocity component
217            %shift and renormalize the angular velocity
218            end
219        end
220        %transform of spatial derivatives
221        if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
222                && isfield(Data,'Dt')
223            if ~isempty(Data.Dt)
224                % estimate the Jacobian matrix DXpx/DXphys
225                for ip=1:length(Data.X)
226                    [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
227                    [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
228                    [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
229                    [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
230                    %Jacobian matrix DXpphys/DXpx
231                    DjXi(1,1)=(Xp1-Xm1);
232                    DjXi(2,1)=(Yp1-Ym1);
233                    DjXi(1,2)=(Xp2-Xm2);
234                    DjXi(2,2)=(Yp2-Ym2);
235                    DjUi(:,:)=Data.DjUi(ip,:,:);
236                    DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
237                    DataOut.DjUi(ip,:,:)=DjUi';
238                end
239                DataOut.DjUi =  DataOut.DjUi/Data.Dt;   %     min(Data.DjUi(:,1,1))=DUDX
240            end
241        end
242    end
243end
244
245
246%%%%%%%%%%%%%%%%%%%%
247function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
248xcorner=[];
249ycorner=[];
250npx=[];
251npy=[];
252NbPoints=20; % nbre of points used to determine the image edge
253for icell=1:length(A)
254    siz=size(A{icell});
255    npx=[npx siz(2)];
256    npy=[npy siz(1)];
257    zphys=0; %default
258    if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
259        if ZIndex==0
260            ZIndex=1;
261        end
262       SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);
263       zphys=SliceCoord(3); %to generalize for non-parallel planes
264    end
265%     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
266%     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
267      edge_x=linspace(0.5,siz(1)-0.5,NbPoints);     
268      edge_y=linspace(0.5,siz(2)-0.5,NbPoints);
269      xima=[edge_y (siz(2)-0.5)*ones(1,NbPoints) edge_y 0.5*ones(1,NbPoints)];%image coordinates of corners
270      yima=[0.5*ones(1,NbPoints) edge_x (siz(1)-0.5)*ones(1,NbPoints) edge_x];%image coordinates of corners
271    [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates
272    %transform the corner coordinates into polar ones   
273    xcorner_new=xcorner_new-origin_xy(1);%shift to the origin of the polar coordinates
274    ycorner_new=ycorner_new-origin_xy(2);%shift to the origin of the polar coordinates       
275    [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta  and X are the polar coordinates angle and radius
276    if (max(theta)-min(theta))>pi   %if the polar origin is inside the image
277        xcorner_new=[0 max(xcorner_new)];
278        theta=[-pi pi];
279    end
280          %shift and renormalize the polar coordinates
281    xcorner_new=xcorner_new-radius_offset;%
282    ycorner_new=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
283    xcorner=[xcorner xcorner_new];
284    ycorner=[ycorner ycorner_new];
285end
286Rangx(1)=min(xcorner);
287Rangx(2)=max(xcorner);
288Rangy(2)=min(ycorner);
289Rangy(1)=max(ycorner);
290% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
291npx=max(npx);
292npy=max(npy);
293x=linspace(Rangx(1),Rangx(2),npx);
294y=linspace(Rangy(1),Rangy(2),npy);
295[X,Y]=meshgrid(x,y);%grid in physical coordinates
296%transform X, Y in cartesian
297X=X+radius_offset;%
298Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius
299[X,Y] = pol2cart(Y,X);
300X=X+origin_xy(1);%shift to the origin of the polar coordinates
301Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
302for icell=1:length(A)
303    siz=size(A{icell});
304    [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
305    XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
306    YIMA=reshape(round(YIMA),1,npx*npy);
307    flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
308    if numel(siz)==2 %(B/W images)
309        vec_A=reshape(A{icell}(:,:,1),1,npx*npy);%put the original image in line
310        ind_in=find(flagin);
311        ind_out=find(~flagin);
312        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
313        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
314        vec_B(ind_in)=vec_A(ICOMB);
315        vec_B(ind_out)=zeros(size(ind_out));
316        A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
317    else
318        for icolor=1:siz(3)
319                vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
320                ind_in=find(flagin);
321                ind_out=find(~flagin);
322                ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
323                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
324                vec_B(ind_in)=vec_A(ICOMB);
325                vec_B(ind_out)=zeros(size(ind_out));
326                A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
327        end
328    end
329end
330
Note: See TracBrowser for help on using the repository browser.