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

Last change on this file since 1127 was 1127, checked in by g7moreau, 10 months ago

Update Joel email

File size: 27.3 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% Data:   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%[Data,Data_1]=phys_polar(varargin)
17%
18% OUTPUT:
19% Data: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
20% Data_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-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
30%   http://www.legi.grenoble-inp.fr
31%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.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 Data=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)
47%------------------------------------------------------------------------
48
49%% request input parameters
50if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
51    prompt = {'origin [x y] of polar coordinates';'reference radius';'reference angle(degrees)';'angle direction and switch x y(+/-)'};
52    dlg_title = 'set the parameters for the polar coordinates';
53    num_lines= 2;
54    def     = { '[0 0]';'';'0';'+'};
55    if isfield(XmlData,'TransformInput')
56        if isfield(XmlData.TransformInput,'PolarCentre')
57            def{1}=num2str(XmlData.TransformInput.PolarCentre);
58        end
59        if isfield(XmlData.TransformInput,'PolarReferenceRadius')
60            def{2}=num2str(XmlData.TransformInput.PolarReferenceRadius);
61        end
62        if isfield(XmlData.TransformInput,'PolarReferenceAngle')
63            def{3}=num2str(XmlData.TransformInput.PolarReferenceAngle);
64        end
65        if isfield(XmlData.TransformInput,'PolarAngleDirection')
66            def{4}=XmlData.TransformInput.PolarAngleDirection;
67        end
68    end
69    answer = inputdlg(prompt,dlg_title,num_lines,def);
70    Data.TransformInput.PolarCentre=str2num(answer{1});
71    Data.TransformInput.PolarReferenceRadius=str2num(answer{2});
72    Data.TransformInput.PolarReferenceAngle=str2num(answer{3});
73    Data.TransformInput.PolarAngleDirection=answer{4};
74    return
75end
76
77%% default outputs
78Data=DataIn; %default output
79if isfield(Data,'CoordUnit')
80    Data=rmfield(Data,'CoordUnit');
81end
82Data.ListVarName = {};
83Data.VarDimName={};
84Data.VarAttribute={};
85DataCell{1}=DataIn;
86Calib{1}=[];
87Slice{1}=[];
88DataCell{2}=[];%default
89checkpixel(1)=0;
90if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel')
91    checkpixel(1)=1;
92end
93if nargin==2||nargin==4
94    if isfield(XmlData,'GeometryCalib') && ~isempty(XmlData.GeometryCalib)&& checkpixel(1)
95        Calib{1}=XmlData.GeometryCalib;
96    end
97    Slice{1}=Calib{1};
98    if isfield(XmlData,'Slice')
99        Slice{1}=XmlData.Slice;
100    end
101    Calib{2}=Calib{1};
102    Slice{2}=Slice{1};
103else
104    Data.Txt='wrong input: need two or four structures';
105end
106nbinput=1;
107if nargin==4% case of two input fields
108    checkpixel(2)=0;
109    if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel')
110        checkpixel(2)=1;
111    end
112    DataCell{2}=DataIn_1;%default
113    if isfield(XmlData_1,'GeometryCalib')&& ~isempty(XmlData_1.GeometryCalib) && checkpixel(2)
114        Calib{2}=XmlData_1.GeometryCalib;
115    end
116    if isfield(XmlData_1,'Slice')
117        Slice{2}=XmlData_1.Slice;
118    end
119    nbinput=2;
120end
121
122%% parameters for polar coordinates (taken from the calibration data of the first field)
123%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
125radius_offset=0;%reference radius used to offset the radial coordinate r
126angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
127angle_scale=180/pi;
128check_reverse=false;
129check_degree=1;%angle expressed in degrees by default
130if isfield(XmlData,'TransformInput')
131    if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre)
132        if isequal(length(XmlData.TransformInput.PolarCentre),2)
133            origin_xy= XmlData.TransformInput.PolarCentre;
134        end
135    end
136    if isfield(XmlData.TransformInput,'PolarReferenceRadius') && ~isempty(XmlData.TransformInput.PolarReferenceRadius)
137        radius_offset=XmlData.TransformInput.PolarReferenceRadius;
138    end
139    if radius_offset > 0
140        angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
141        check_degree=0; %the output has the same unit as the input
142    else
143        angle_scale=180/pi; %polar angle in degrees
144        check_degree=1;%angle expressed in degrees
145    end
146    if isfield(XmlData.TransformInput,'PolarReferenceAngle') && isnumeric(XmlData.TransformInput.PolarReferenceAngle)
147        angle_offset=(pi/180)*XmlData.TransformInput.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))
148    end
149    check_reverse=isfield(XmlData.TransformInput,'PolarAngleDirection')&& strcmp(XmlData.TransformInput.PolarAngleDirection,'-');
150end
151
152%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153%% get fields
154
155nbvar=0;%counter for the number of output variables
156nbcoord=0;%counter for the number of variablecheck_degrees for radial coordiantes (case of multiple field inputs)
157nbgrid=0;%counter for the number of gridded fields (all linearly interpolated on the same output polar grid)
158nbscattered=0;%counter of scattered fields
159radius_name='radius';
160theta_name='theta';
161U_r_name='U_r';
162U_theta_name='U_theta';
163for ifield=1:nbinput %1 or 2 input fields
164    [CellInfo,NbDim,errormsg]=find_field_cells(DataCell{ifield});
165    if ~isempty(errormsg)
166        Data.Txt=['bad input to phys_polar: ' errormsg];
167        return
168    end
169    %transform of X,Y coordinates for vector fields
170    if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex)
171        ZIndex=DataCell{ifield}.ZIndex;
172    else
173        ZIndex=0;
174    end
175    check_scalar=zeros(1,numel(CellInfo));
176    check_vector=zeros(1,numel(CellInfo));
177    for icell=1:numel(CellInfo)
178        if NbDim(icell)==2
179            % case of input field with scattered coordinates
180            if strcmp(CellInfo{icell}.CoordType,'scattered')
181                nbscattered=nbscattered+1;
182                nbcoord=nbcoord+1;
183                radius_name = rename_indexing(radius_name,Data.ListVarName);
184                theta_name = rename_indexing(theta_name,Data.ListVarName);
185                Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];
186                dim_name = rename_indexing('nb_point',Data.VarDimName);
187                Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];
188                nbvar=nbvar+2;
189                Data.VarAttribute{nbvar-1}.Role='coord_x';
190                check_unit=1;
191                %unit of output field
192                if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
193                    radius_unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
194                elseif isfield(DataCell{ifield},'CoordUnit')
195                    radius_unit=DataCell{ifield}.CoordUnit;
196                else
197                    radius_unit='';
198                end
199                Data.VarAttribute{nbvar-1}.units=radius_unit;
200                if check_degree
201                     Data.VarAttribute{nbvar}.units='degree';
202                else %case of a reference radius
203                    Data.VarAttribute{nbvar}.units=radius_unit;
204                    Data.CoordUnit=radius_unit;
205                end
206%                 if isfield(DataCell{ifield},'CoordUnit')
207%                     Data=rmfield(Data,'CoordUnit');
208%                     Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
209%                 elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
210%                     Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
211%                 else
212%                     check_unit=0;
213%                 end
214                Data.VarAttribute{nbvar}.Role='coord_y';
215%                 if check_degree
216%                 Data.VarAttribute{nbvar}.units='degree';
217%                 elseif check_unit
218%                     Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units;
219%                 end
220 
221                %transform u,v into polar coordinates
222                X=DataCell{ifield}.(CellInfo{icell}.XName);
223                Y=DataCell{ifield}.(CellInfo{icell}.YName);
224                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
225                    UName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};
226                    VName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};
227                    if ~isempty(Calib{ifield})
228                        [X,Y,Z,DataCell{ifield}.(UName),DataCell{ifield}.(VName)]=...
229                            phys_XYUV(DataCell{ifield},Calib{ifield},Slice{ifield},ZIndex);
230                    end
231                end
232                [Theta,Radius] = cart2pol(X-origin_xy(1),Y-origin_xy(2));
233                Data.(radius_name)=Radius-radius_offset;
234                Data.(theta_name)=Theta*angle_scale-angle_offset;
235                if Z~=0
236                    Data.Z=Z;
237                    nbvar=nbvar+1;
238                    Data.ListVarName = [Data.ListVarName {'Z'}];
239                    Data.VarDimName=[Data.VarDimName {dim_name}];
240                    Data.VarAttribute{nbvar}.Role='coord_z';
241                end
242                if isfield(CellInfo{icell},'VarIndex_scalar')
243                    ScalarName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};
244                    ScalarName=rename_indexing(ScalarName,Data.ListVarName);
245                    Data.(ScalarName)=DataCell{ifield}.(ScalarName);
246                    nbvar=nbvar+1;
247                    Data.ListVarName = [Data.ListVarName {ScalarName}];
248                    Data.VarDimName=[Data.VarDimName {dim_name}];
249                    Data.VarAttribute{nbvar}.Role='scalar';
250                end
251                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
252                    U_r_name= rename_indexing(U_r_name,Data.ListVarName);
253                    U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);
254                    Data.(U_r_name)=DataCell{ifield}.(UName).*cos(Theta)+DataCell{ifield}.(VName).*sin(Theta);%radial velocity
255                    Data.(U_theta_name)=(-DataCell{ifield}.(UName).*sin(Theta)+DataCell{ifield}.(VName).*cos(Theta));%./(Data.X)%+radius_ref);% azimuthal velocity component
256                    Data.ListVarName = [Data.ListVarName {U_r_name} {U_theta_name}];
257                    Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];
258                    Data.VarAttribute{nbvar+1}.Role='vector_x';
259                    Data.VarAttribute{nbvar+2}.Role='vector_y';
260                    nbvar=nbvar+2;
261                end
262                if isfield(CellInfo{icell},'VarIndex_errorflag')
263                    error_flag_name=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_errorflag};
264                    error_flag_newname= rename_indexing(error_flag_name,Data.ListVarName);
265                    Data.(error_flag_newname)=DataCell{ifield}.(error_flag_name);
266                    Data.ListVarName = [Data.ListVarName {error_flag_newname}];
267                    Data.VarDimName=[Data.VarDimName {dim_name}];
268                    nbvar=nbvar+1;
269                    Data.VarAttribute{nbvar}.Role='errorflag';
270                end
271               
272           %caseof input fields on gridded coordinates (matrix)
273            elseif strcmp(CellInfo{icell}.CoordType,'grid')
274                if nbgrid==0% no gridded data yet, introduce the coordinate variables common to all gridded data
275                    nbcoord=nbcoord+1;%add new radial coordinates for the first gridded field
276                    radius_name = rename_indexing(radius_name,Data.ListVarName);% add an index to the name, or increment an existing index,
277                    theta_name = rename_indexing(theta_name,Data.ListVarName);% if the proposed Name already exists in the list
278                    Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];%add polar coordinates to the list of variables
279                    Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}];
280                    nbvar=nbvar+2;
281                    if check_reverse
282                        Data.VarAttribute{nbvar-1}.Role='coord_y';
283                        Data.VarAttribute{nbvar}.Role='coord_x';
284                    else
285                        Data.VarAttribute{nbvar-1}.Role='coord_x';
286                        Data.VarAttribute{nbvar}.Role='coord_y';
287                    end
288                    check_unit=1;
289
290                    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
291                        Data.VarAttribute{nbvar-1}.units=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
292                    elseif isfield(DataCell{ifield},'CoordUnit')
293                        Data.VarAttribute{nbvar-1}.units=DataCell{ifield}.CoordUnit;%radius in coord units
294                    else
295                        check_unit=0;
296                    end
297                    if check_degree
298                        Data.VarAttribute{nbvar}.units='degree';%angle in degree
299                    elseif check_unit
300                        Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units;% angle in coord unit (normalised by reference radiuss)
301                    end
302                end
303                if isfield(CellInfo{icell},'VarIndex_scalar')
304                    nbgrid=nbgrid+1;
305                    nbvar=nbvar+1;
306                    Data.VarAttribute{nbvar}.Role='scalar';
307                    FieldName{nbgrid}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};
308                    A{nbgrid}=DataCell{ifield}.(FieldName{nbgrid});
309                    nbpoint(nbgrid)=numel(A{nbgrid});
310                    check_scalar(nbgrid)=1;
311                    coord_x{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
312                    coord_y{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
313                    ZInd(nbgrid)=ZIndex;
314                    Calib_new{nbgrid}=Calib{ifield};
315                    Slice_new{nbgrid}=Slice{ifield};
316                end
317                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
318                    FieldName{nbgrid+1}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};
319                    FieldName{nbgrid+2}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};
320                    A{nbgrid+1}=DataCell{ifield}.(FieldName{nbgrid+1});
321                    A{nbgrid+2}=DataCell{ifield}.(FieldName{nbgrid+2});
322                    % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}];
323                    %Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];
324                    Data.VarAttribute{nbvar+1}.Role='vector_x';
325                    Data.VarAttribute{nbvar+2}.Role='vector_y';
326                    nbpoint([nbgrid+1 nbgrid+2])=numel(A{nbgrid+1});
327                    check_vector(nbgrid+1)=1;
328                    check_vector(nbgrid+2)=1;
329                    coord_x{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
330                    coord_y{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
331                    coord_x{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
332                    coord_y{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
333                    ZInd(nbgrid+1)=ZIndex;
334                    ZInd(nbgrid+2)=ZIndex;
335                    Calib_new{nbgrid+1}=Calib{ifield};
336                    Calib_new{nbgrid+2}=Calib{ifield};
337                    Slice_new{nbgrid+1}=Calib{ifield};
338                    Slice_new{nbgrid+2}=Calib{ifield};
339                    nbgrid=nbgrid+2;
340                    nbvar=nbvar+2;
341                end
342            end
343        end
344    end
345end
346
347%% tranform cartesian to polar coordinates for gridded data
348if nbgrid~=0
349    [A,Data.radius,Data.theta]=phys_Ima_polar(A,coord_x,coord_y,Calib_new,Slice_new,ZInd,origin_xy,radius_offset,angle_offset,angle_scale);
350    for icell=1:numel(A)
351        if icell<=numel(A)-1 && check_vector(icell)==1 && check_vector(icell+1)==1   %transform u,v into polar coordinates
352            theta=Data.theta/angle_scale-angle_offset;
353            [~,Theta]=meshgrid(Data.radius,theta);%grid in physical coordinates
354            U_r_name= rename_indexing(U_r_name,Data.ListVarName);
355            U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);       
356                Data.(U_r_name)=A{icell}.*cos(Theta)+A{icell+1}.*sin(Theta);%radial velocity
357                Data.(U_theta_name)=(-A{icell}.*sin(Theta)+A{icell+1}.*cos(Theta));% azimuthal velocity component
358            if check_reverse
359                Data.(U_theta_name)=(Data.(U_theta_name))';
360                Data.(U_r_name)=Data.(U_r_name)';
361                Data.ListVarName=[Data.ListVarName {U_theta_name,U_r_name}];
362                Data.VarDimName=[Data.VarDimName {{radius_name,theta_name}} {{radius_name,theta_name}}];
363            else
364                Data.ListVarName=[Data.ListVarName {U_r_name,U_theta_name}];
365                Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];
366            end
367        elseif ~check_vector(icell)% for scalar fields
368            FieldName{icell}= rename_indexing(FieldName{icell},Data.ListVarName);
369            Data.ListVarName=[Data.ListVarName FieldName(icell)];       
370            if check_reverse
371                Data.(FieldName{icell})=A{icell}';
372                Data.VarDimName=[Data.VarDimName {{radius_name,theta_name}}];
373            else
374                Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];
375                Data.(FieldName{icell})=A{icell};
376            end
377        end
378    end
379end
380if check_reverse
381    Data.(theta_name)=-Data.(theta_name);
382end
383
384
385%------------------------------------------------
386%--- transform a single field into phys coordiantes
387function [X,Y,Z,U,V]=phys_XYUV(Data,Calib,Slice,ZIndex)
388%------------------------------------------------
389%% set default output
390%DataOut=Data;%default
391%DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters
392X=[];%default output
393Y=[];
394Z=0;
395U=[];
396V=[];
397%% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
398if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
399    [X,Y,Z]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);
400    Dt=1; %default
401    if isfield(Data,'dt')&&~isempty(Data.dt)
402        Dt=Data.dt;
403    end
404    if isfield(Data,'Dt')&&~isempty(Data.Dt)
405        Dt=Data.Dt;
406    end
407    if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
408        [XOut_1,YOut_1]=phys_XYZ(Calib,Slice,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
409        [XOut_2,YOut_2]=phys_XYZ(Calib,Slice,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
410        U=(XOut_2-XOut_1)/Dt;
411        V=(YOut_2-YOut_1)/Dt;
412    end
413end
414
415%%%%%%%%%%%%%%%%%%%%
416% tranform gridded field into polar coordiantes on a regular polar grid,
417% transform to phys coordiantes if requested by calibration input
418function [A_out,radius,theta]=phys_Ima_polar(A,coord_x,coord_y,CalibIn,SliceIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
419rcorner=[];
420thetacorner=[];
421npx=[];
422npy=[];
423for icell=1:length(A)
424    siz=size(A{icell});
425    npx(icell)=siz(2);
426    npy(icell)=siz(1);
427    x_edge=[linspace(coord_x{icell}(1),coord_x{icell}(end),npx(icell)) coord_x{icell}(end)*ones(1,npy(icell))...
428        linspace(coord_x{icell}(end),coord_x{icell}(1),npx(icell)) coord_x{icell}(1)*ones(1,npy(icell))];%x coordinates of the image edge(four sides)
429    y_edge=[coord_y{icell}(1)*ones(1,npx(icell)) linspace(coord_y{icell}(1),coord_y{icell}(end),npy(icell))...
430        coord_y{icell}(end)*ones(1,npx(icell)) linspace(coord_y{icell}(end),coord_y{icell}(1),npy(icell))];%y coordinates of the image edge(four sides)
431   
432    % transform edges into phys coordinates if requested
433    if ~isempty(CalibIn{icell})
434        [x_edge,y_edge]=phys_XYZ(CalibIn{icell},SliceIn{icell},x_edge,y_edge,ZIndex(icell));% physical coordinates of the image edge
435    end
436   
437    %transform the corner coordinates into polar ones
438    x_edge=x_edge-origin_xy(1);%shift to the origin of the polar coordinates
439    y_edge=y_edge-origin_xy(2);%shift to the origin of the polar coordinates
440    [theta_edge,r_edge] = cart2pol(x_edge,y_edge);%theta  and X are the polar coordinates angle and radius
441    if (max(theta_edge)-min(theta_edge))>pi   %if the polar origin is inside the image
442        r_edge=[0 max(r_edge)];
443        theta_edge=[-pi pi];
444    end
445    rcorner=[rcorner r_edge];
446    thetacorner=[thetacorner theta_edge];
447end
448nbpoint=max(npx.*npy);
449Min_r=min(rcorner);
450Max_r=max(rcorner);
451Min_theta=min(thetacorner)*angle_scale;
452Max_theta=max(thetacorner)*angle_scale;
453Dr=round_uvmat((Max_r-Min_r)/sqrt(nbpoint));
454Dtheta=round_uvmat((Max_theta-Min_theta)/sqrt(nbpoint));% get a simple mesh for the rescaled angle
455radius=Min_r:Dr:Max_r;% polar coordinates for projections
456theta=Min_theta:Dtheta:Max_theta;
457%theta=Max_theta:-Dtheta:Min_theta;
458[Radius,Theta]=meshgrid(radius,theta/angle_scale);%grid in polar coordinates (angles in radians)
459%transform X, Y in cartesian
460[X,Y] = pol2cart(Theta,Radius);% cartesian coordinates associated to the grid in polar coordinates
461X=X+origin_xy(1);%shift to the origin of the polar coordinates
462Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
463radius=radius-radius_offset;
464theta=theta-angle_offset*angle_scale;
465[np_theta,np_r]=size(Radius);
466
467for icell=1:length(A)
468    XIMA=X;
469    YIMA=Y;
470    if ~isempty(CalibIn{icell})%transform back to pixel if calibration parameters are introduced
471        Z=0; %default
472        if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
473            if ZIndex(icell)==0
474                ZIndex(icell)=1;
475            end
476            SliceCoord=CalibIn{icell}.SliceCoord(ZIndex(icell),:);
477            Z=SliceCoord(3); %to generalize for non-parallel planes
478            if isfield(CalibIn{icell},'SliceAngle')
479            norm_plane=angle2normal(CalibIn{icell}.SliceAngle);
480            Z=Z-(norm_plane(1)*(X-SliceCoord(1))+norm_plane(2)*(Y-SliceCoord(2)))/norm_plane(3);
481            end
482        end
483        [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,Z);%corresponding image indices for each point in the real space grid
484    end
485    Dx=(coord_x{icell}(end)-coord_x{icell}(1))/(npx(icell)-1);
486    Dy=(coord_y{icell}(end)-coord_y{icell}(1))/(npy(icell)-1);
487    indx_ima=1+round((XIMA-coord_x{icell}(1))/Dx);%indices of the initial matrix close to the points of the new grid
488    %indy_ima=1+round((YIMA-coord_y{icell}(1))/Dy);
489    indy_ima=1+round((coord_y{icell}(end)-YIMA)/Dy);
490     Delta_x=1+(XIMA-coord_x{icell}(1))/Dx-indx_ima;%error in the index discretisation
491     Delta_y=1+(coord_y{icell}(end)-YIMA)/Dy-indy_ima;
492    XIMA=reshape(indx_ima,1,[]);%indices reorganized in 'line'
493    YIMA=reshape(indy_ima,1,[]);%indices reorganized in 'line'
494    flagin=XIMA>=1 & XIMA<=npx(icell) & YIMA >=1 & YIMA<=npy(icell);%flagin=1 inside the original image
495    siz=size(A{icell});
496    checkuint8=isa(A{icell},'uint8');%check for image input with 8 bits
497    checkuint16=isa(A{icell},'uint16');%check for image input with 16 bits
498    A{icell}=double(A{icell});
499    if numel(siz)==2 %(B/W images)
500        vec_A=reshape(A{icell}(:,:,1),1,[]);%put the original image in line
501        ind_in=find(flagin);
502        ind_out=find(~flagin);
503        ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));% indices in vec_A
504        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
505        vec_B(ind_in)=vec_A(ICOMB);
506        vec_B(ind_out)=zeros(size(ind_out));
507        A_out{icell}=reshape(vec_B,np_theta,np_r);%new image in real coordinates
508        DA_y=circshift(A_out{icell},-1,1)-A_out{icell};% derivative
509        DA_y(end,:)=0;
510        DA_x=circshift(A_out{icell},-1,2)-A_out{icell};
511        DA_x(:,end)=0;
512        A_out{icell}=A_out{icell}+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation
513    else
514        for icolor=1:siz(3)
515            vec_A=reshape(A{icell}(:,:,icolor),1,[]);%put the original image in line
516            ind_in=find(flagin);
517            ind_out=find(~flagin);
518            ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
519            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
520            vec_B(ind_in)=vec_A(ICOMB);
521            vec_B(ind_out)=zeros(size(ind_out));
522            A_out{icell}(:,:,icolor)=reshape(vec_B,np_theta,np_r);%new image in real coordinates
523            DA_y=circshift(A_out{icell}(:,:,icolor),-1,1)-A_out{icell}(:,:,icolor);
524            DA_y(end,:)=0;
525            DA_x=circshift(A_out{icell}(:,:,icolor),-1,2)-A_out{icell}(:,:,icolor);
526            DA_x(:,end)=0;
527            A_out{icell}(:,:,icolor)=A_out{icell}(:,:,icolor)+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation
528        end
529    end
530    if checkuint8
531        A_out{icell}=uint8(A_out{icell});
532    elseif checkuint16
533        A_out{icell}=uint16(A_out{icell});
534    end
535end
536
Note: See TracBrowser for help on using the repository browser.