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

Last change on this file since 574 was 574, checked in by sommeria, 8 years ago

clean the transform fcts

File size: 11.9 KB
Line 
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%------------------------------------------------------------------------
29function DataOut=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)
30%------------------------------------------------------------------------
31Calib{1}=[];
32if 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};
41else
42    DataOut.Txt='wrong input: need two or four structures';
43end
44test_1=0;
45if 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
53end
54
55%parameters for polar coordinates (taken from the calibration data of the first field)
56%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
58if isfield(CalibData,'PolarCentre') && isnumeric(CalibData.PolarCentre)
59    if isequal(length(CalibData.PolarCentre),2);
60        origin_xy= CalibData.PolarCentre;
61    end
62end
63radius_offset=0;%reference radius used to offset the radial coordinate r
64angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
65if isfield(CalibData,'PolarReferenceRadius') && isnumeric(CalibData.PolarReferenceRadius)
66    radius_offset=CalibData.PolarReferenceRadius;
67end
68if radius_offset > 0
69    angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
70else
71    angle_scale=180/pi; %polar angle in degrees
72end
73if 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))
75end
76% new x coordinate = radius-radius_offset;
77% new y coordinate = theta*angle_scale-angle_offset
78
79%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80
81iscalar=0;
82%transform first field to cartesian phys coordiantes
83if  ~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
97end
98
99%transform second field (if exists) to cartesian phys coordiantes
100if 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
118end
119if 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
129end
130
131
132
133
134%------------------------------------------------
135function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
136
137DataOut=Data;
138% DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
139if isfield(Calib,'SliceCoord')
140    DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane
141end
142
143if 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
200end
201
202
203%%%%%%%%%%%%%%%%%%%%
204function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
205xcorner=[];
206ycorner=[];
207npx=[];
208npy=[];
209for 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];
234end
235Rangx(1)=min(xcorner);
236Rangx(2)=max(xcorner);
237Rangy(2)=min(ycorner);
238Rangy(1)=max(ycorner);
239% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
240npx=max(npx);
241npy=max(npy);
242x=linspace(Rangx(1),Rangx(2),npx);
243y=linspace(Rangy(1),Rangy(2),npy);
244[X,Y]=meshgrid(x,y);%grid in physical coordinates
245%transform X, Y in cartesian
246X=X+radius_offset;%
247Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius
248[X,Y] = pol2cart(Y,X);
249X=X+origin_xy(1);%shift to the origin of the polar coordinates
250Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
251for 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
278end
279
Note: See TracBrowser for help on using the repository browser.