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

Last change on this file since 1081 was 1081, checked in by sommeria, 4 years ago

phys_polar corrected

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