source: trunk/src/transform_field/phys.m @ 1160

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

Update Joel email

File size: 10.0 KB
Line 
1%'phys': transforms image (Unit='pixel') to real world (phys) coordinates using geometric calibration parameters.  It acts if the input field contains the tag 'CoordTUnit' with value 'pixel'
2%------------------------------------------------------------------------
3%%%%  Use the general syntax for transform fields %%%%
4% OUTPUT:
5% DataOut:   output field structure
6%
7%INPUT:
8% DataIn:  first input field structure
9% XmlData: first input parameter structure,
10%        .GeometryCalib: substructure of the calibration parameters
11% DataIn_1: optional second input field structure
12% XmlData_1: optional second input parameter structure
13%         .GeometryCalib: substructure of the calibration parameters
14
15%=======================================================================
16% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
17%   http://www.legi.grenoble-inp.fr
18%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
19%
20%     This file is part of the toolbox UVMAT.
21%
22%     UVMAT is free software; you can redistribute it and/or modify
23%     it under the terms of the GNU General Public License as published
24%     by the Free Software Foundation; either version 2 of the license,
25%     or (at your option) any later version.
26%
27%     UVMAT is distributed in the hope that it will be useful,
28%     but WITHOUT ANY WARRANTY; without even the implied warranty of
29%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30%     GNU General Public License (see LICENSE.txt) for more details.
31%=======================================================================
32
33function DataOut=phys(DataIn,XmlData,DataIn_1,XmlData_1)
34%------------------------------------------------------------------------
35
36% A FAIRE: 1- verifier si DataIn est une 'field structure'(.ListVarName'):
37% chercher ListVarAttribute, for each field (cell of variables):
38%   .CoordType: 'phys' or 'px'   (default==phys, no transform)
39%   .scale_factor: =dt (to transform displacement into velocity) default=1
40%   .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)
41%   (default='coord' if .Role='coord_x,_y...,
42%            'D_i' if '.Role='vector_x,...',
43%              'scalar', else (thenno change except scale factor)
44
45DataOut=[];
46DataOut_1=[]; %default second  output field
47if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
48    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
49        DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
50    end
51    return
52end
53
54%% analyse input and set default output
55DataOut=DataIn;%default first output field
56if nargin>=2 % nargin =nbre of input variables
57     Calib{1}=[];
58    if isfield(XmlData,'GeometryCalib')
59        Calib{1}=XmlData.GeometryCalib;
60    end
61    Slice{1}=Calib{1};
62    if isfield(XmlData,'Slice')
63        Slice{1}=XmlData.Slice;
64    end
65    if nargin>=3  %two input fields
66        DataOut_1=DataIn_1;%default second output field
67        Calib{2}=Calib{1};
68        if nargin>=4
69            if isfield(XmlData_1,'GeometryCalib')
70                Calib{2}=XmlData_1.GeometryCalib;
71            end
72            Slice{2}=Calib{2};
73            if isfield(XmlData_1,'Slice')
74                Slice{2}=XmlData_1.Slice;
75            end
76        end
77    end
78end
79
80%% get the z index defining the section plane
81ZIndex=1;
82if isfield(DataIn,'ZIndex')&&~isempty(DataIn.ZIndex)&&~isnan(DataIn.ZIndex)
83    ZIndex=DataIn.ZIndex;
84end
85
86%% transform first field
87iscalar=0;% counter of scalar fields
88checktransform=0;
89if  ~isempty(Calib{1})
90    if isfield(Calib{1},'CalibrationType')&& isfield(Calib{1},'CoordUnit') && isfield(DataIn,'CoordUnit')&& strcmp(DataIn.CoordUnit,'pixel')   
91        DataOut=phys_1(DataIn,Calib{1},Slice{1},ZIndex);% transform coordinates and velocity components
92        %case of images or scalar: in case of two input fields, we need to project the transform  on the same regular grid
93        if isfield(DataIn,'A') && isfield(DataIn,'Coord_x') && ~isempty(DataIn.Coord_x) && isfield(DataIn,'Coord_y')&&...
94                ~isempty(DataIn.Coord_y) && length(DataIn.A)>1
95            iscalar=1;
96            A{1}=DataIn.A;
97        end
98        checktransform=1;
99    end
100end
101
102%% document the selected  plane position and angle if relevant
103if  checktransform && isfield(Slice{1},'SliceCoord')&&size(Slice{1}.SliceCoord,1)>=ZIndex
104    DataOut.PlaneCoord=Slice{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
105    if isfield(Slice{1},'SliceAngle') % transfer the slice rotation angles
106        if isequal(size(Slice{1}.SliceAngle,1),1)% case of a unique angle
107            DataOut.PlaneAngle=Slice{1}.SliceAngle;
108        else  % case of multiple planes with different angles: select the plane with index ZIndex
109            DataOut.PlaneAngle=Slice{1}.SliceAngle(ZIndex,:);
110        end
111    end
112end
113
114%% transform second field if relevant
115checktransform_1=0;
116if ~isempty(DataOut_1)
117    if isfield(DataIn_1,'ZIndex') && ~isequal(DataIn_1.ZIndex,ZIndex)
118        DataOut_1.Txt='different plane indices for the two input fields';
119        return
120    end
121    if isfield(Calib{2},'CalibrationType')&&isfield(Calib{2},'CoordUnit') && isfield(DataIn_1,'CoordUnit')&& strcmp(DataIn_1.CoordUnit,'pixel')
122        DataOut_1=phys_1(DataOut_1,Calib{2},Slice{2},ZIndex);
123        if isfield(Slice{2},'SliceCoord')
124            if ~(isfield(Slice{2},'SliceCoord') && isequal(Slice{2}.SliceCoord,Slice{1}.SliceCoord))
125                DataOut_1.Txt='different plane positions for the two input fields';
126                return
127            end
128            DataOut_1.PlaneCoord=DataOut.PlaneCoord;% same plane position for the two input fields
129            if isfield(Slice{1},'SliceAngle')
130                if ~(isfield(Slice{2},'SliceAngle') && isequal(Slice{2}.SliceAngle,Slice{1}.SliceAngle))
131                    DataOut_1.Txt='different plane angles for the two input fields';
132                    return
133                end
134                DataOut_1.PlaneAngle=DataOut.PlaneAngle; % same plane angle for the two input fields
135            end
136        end
137        if isfield(DataIn_1,'A')&&isfield(DataIn_1,'Coord_x')&&~isempty(DataIn_1.Coord_x) && isfield(DataIn_1,'Coord_y')&&...
138                ~isempty(DataIn_1.Coord_y)&&length(DataIn_1.A)>1
139            iscalar=iscalar+1;
140%             Calib{iscalar}=Calib{2};
141            A{iscalar}=DataIn_1.A;
142        end
143        checktransform_1=1;
144    end
145end
146
147%% transform the scalar(s) or image(s)
148if checktransform && iscalar~=0
149    [A,Coord_x,Coord_y]=phys_ima(A,XmlData,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima
150    if iscalar==1 && ~isempty(DataOut_1) % case for which only the second field is a scalar
151         DataOut_1.A=A{1};
152         DataOut_1.Coord_x=Coord_x;
153         DataOut_1.Coord_y=Coord_y;
154    else
155        DataOut.A=A{1};
156        DataOut.Coord_x=Coord_x;
157        DataOut.Coord_y=Coord_y;
158    end
159    if iscalar==2
160        DataOut_1.A=A{2};
161        DataOut_1.Coord_x=Coord_x;
162        DataOut_1.Coord_y=Coord_y;
163    end
164end
165
166% subtract fields
167if ~isempty(DataOut_1)
168    DataOut=sub_field(DataOut,[],DataOut_1);
169end
170%------------------------------------------------
171%--- transform a single field
172function DataOut=phys_1(Data,Calib,Slice,ZIndex)
173%------------------------------------------------
174%% set default output
175DataOut=Data;%default
176DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters
177
178%% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
179if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
180  [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);
181    Dt=1; %default
182    if isfield(Data,'dt')&&~isempty(Data.dt)
183        Dt=Data.dt;
184    end
185    if isfield(Data,'Dt')&&~isempty(Data.Dt)
186        Dt=Data.Dt;
187    end
188    if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
189        [XOut_1,YOut_1]=phys_XYZ(Calib,Slice,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
190        [XOut_2,YOut_2]=phys_XYZ(Calib,Slice,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
191        DataOut.U=(XOut_2-XOut_1)/Dt;
192        DataOut.V=(YOut_2-YOut_1)/Dt;
193    end
194end
195
196%% suppress tps
197list_tps={'Coord_tps'  'U_tps'  'V_tps'  'SubRange'  'NbSites'};
198ind_remove=[];
199for ilist=1:numel(list_tps)
200    ind_tps=find(strcmp(list_tps{ilist},Data.ListVarName));
201    if ~isempty(ind_tps)
202        ind_remove=[ind_remove ind_tps];
203        DataOut=rmfield(DataOut,list_tps{ilist});
204    end
205end
206if isfield(DataOut,'VarAttribute') && numel(DataOut.VarAttribute)>=3 && isfield(DataOut.VarAttribute{3},'VarIndex_tps')
207    DataOut.VarAttribute{3}=rmfield(DataOut.VarAttribute{3},'VarIndex_tps');
208end
209if isfield(DataOut,'VarAttribute')&& numel(DataOut.VarAttribute)>=4 && isfield(DataOut.VarAttribute{4},'VarIndex_tps')
210    DataOut.VarAttribute{4}=rmfield(DataOut.VarAttribute{4},'VarIndex_tps');
211end
212if ~isempty(ind_remove)
213    DataOut.ListVarName(ind_remove)=[];
214    DataOut.VarDimName(ind_remove)=[];
215    DataOut.VarAttribute(ind_remove)=[];
216end
217   
218%% transform of spatial derivatives: TODO check the case with plane angles
219if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)
220    % estimate the Jacobian matrix DXpx/DXphys
221    for ip=1:length(Data.X)
222        [Xp1,Yp1]=phys_XYZ(Calib,Slice,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
223        [Xm1,Ym1]=phys_XYZ(Calib,Slice,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
224        [Xp2,Yp2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
225        [Xm2,Ym2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
226        %Jacobian matrix DXpphys/DXpx
227        DjXi(1,1)=(Xp1-Xm1);
228        DjXi(2,1)=(Yp1-Ym1);
229        DjXi(1,2)=(Xp2-Xm2);
230        DjXi(2,2)=(Yp2-Ym2);
231        DjUi(:,:)=Data.DjUi(ip,:,:);
232        DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
233        DataOut.DjUi(ip,:,:)=DjUi';
234    end
235    DataOut.DjUi =  DataOut.DjUi/Dt;   %     min(Data.DjUi(:,1,1))=DUDX
236end
237
238
239%%%%%%%%%%%%%%%%%%%%
240
Note: See TracBrowser for help on using the repository browser.