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