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

Last change on this file since 730 was 695, checked in by sommeria, 10 years ago

phys corrected with linear interpolation of images

File size: 13.8 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%------------------------------------------------------------------------
4%%%%  Use the general syntax for transform fields %%%%
5% OUTPUT:
6% DataOut:   output field structure
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%------------------------------------------------------------------------
16function DataOut=phys(DataIn,XmlData,DataIn_1,XmlData_1)
17%------------------------------------------------------------------------
18
19% A FAIRE: 1- verifier si DataIn est une 'field structure'(.ListVarName'):
20% chercher ListVarAttribute, for each field (cell of variables):
21%   .CoordType: 'phys' or 'px'   (default==phys, no transform)
22%   .scale_factor: =dt (to transform displacement into velocity) default=1
23%   .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)
24%   (default='coord' if .Role='coord_x,_y...,
25%            'D_i' if '.Role='vector_x,...',
26%              'scalar', else (thenno change except scale factor)
27
28DataOut=[];
29DataOut_1=[]; %default second  output field
30if strcmp(DataIn,'*')
31    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
32        DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
33    end
34    return
35end
36
37%% analyse input and set default output
38DataOut=DataIn;%default first output field
39if nargin>=2 % nargin =nbre of input variables
40    if isfield(XmlData,'GeometryCalib')
41        Calib{1}=XmlData.GeometryCalib;
42    else
43        Calib{1}=[];
44    end
45    if nargin>=3  %two input fields
46        DataOut_1=DataIn_1;%default second output field
47        if nargin>=4 && isfield(XmlData_1,'GeometryCalib')
48            Calib{2}=XmlData_1.GeometryCalib;
49        else
50            Calib{2}=Calib{1};
51        end
52    end
53end
54
55%% get the z index defining the section plane
56if isfield(DataIn,'ZIndex')&&~isempty(DataIn.ZIndex)&&~isnan(DataIn.ZIndex)
57    ZIndex=DataIn.ZIndex;
58else
59    ZIndex=1;
60end
61
62%% transform first field
63iscalar=0;% counter of scalar fields
64if  ~isempty(Calib{1})
65    if ~isfield(Calib{1},'CalibrationType')||~isfield(Calib{1},'CoordUnit')
66        return %bad calib parameter input
67    end
68    if ~(isfield(DataIn,'CoordUnit')&& strcmp(DataIn.CoordUnit,'pixel'))
69        return % transform only fields in pixel coordinates
70    end
71    DataOut=phys_1(DataIn,Calib{1},ZIndex);% transform coordiantes and velocity components
72    %case of images or scalar: in case of two input fields, we need to project the transform  on the same regular grid
73    if isfield(DataIn,'A') && isfield(DataIn,'AX') && ~isempty(DataIn.AX) && isfield(DataIn,'AY')&&...
74                                           ~isempty(DataIn.AY) && length(DataIn.A)>1
75        iscalar=1;
76        A{1}=DataIn.A;
77    end
78end
79
80%% document the selected  plane position and angle if relevant
81if isfield(Calib{1},'SliceCoord')&&size(Calib{1}.SliceCoord,1)>=ZIndex
82    DataOut.PlaneCoord=Calib{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
83    if isfield(Calib{1},'SliceAngle') % transfer the slice rotation angles
84        if isequal(size(Calib{1}.SliceAngle,1),1)% case of a unique angle
85            DataOut.PlaneAngle=Calib{1}.SliceAngle;
86        else  % case of multiple planes with different angles: select the plane with index ZIndex
87            DataOut.PlaneAngle=Calib{1}.SliceAngle(ZIndex,:);
88        end
89    end
90end
91
92%% transform second field if relevant
93if ~isempty(DataOut_1)
94    if isfield(DataIn_1,'ZIndex') && ~isequal(DataIn_1.ZIndex,ZIndex)
95        DataOut_1.Txt='different plane indices for the two input fields';
96        return
97    end
98    if ~isfield(Calib{2},'CalibrationType')||~isfield(Calib{2},'CoordUnit')
99        return %bad calib parameter input
100    end
101    if ~(isfield(DataIn_1,'CoordUnit')&& strcmp(DataIn_1.CoordUnit,'pixel'))
102        return % transform only fields in pixel coordinates
103    end
104    DataOut_1=phys_1(DataOut_1,Calib{2},ZIndex);
105    if isfield(Calib{1},'SliceCoord')
106        if ~(isfield(Calib{2},'SliceCoord') && isequal(Calib{2}.SliceCoord,Calib{1}.SliceCoord))
107            DataOut_1.Txt='different plane positions for the two input fields';
108            return
109        end       
110        DataOut_1.PlaneCoord=DataOut.PlaneCoord;% same plane position for the two input fields
111        if isfield(Calib{1},'SliceAngle')
112            if ~(isfield(Calib{2},'SliceAngle') && isequal(Calib{2}.SliceAngle,Calib{1}.SliceAngle))
113                DataOut_1.Txt='different plane angles for the two input fields';
114                return
115            end
116            DataOut_1.PlaneAngle=DataOut.PlaneAngle; % same plane angle for the two input fields
117        end
118    end
119    if isfield(DataIn_1,'A')&&isfield(DataIn_1,'AX')&&~isempty(DataIn_1.AX) && isfield(DataIn_1,'AY')&&...
120            ~isempty(DataIn_1.AY)&&length(DataIn_1.A)>1
121        iscalar=iscalar+1;
122        Calib{iscalar}=Calib{2};
123        A{iscalar}=DataIn_1.A;
124    end
125end
126
127%% transform the scalar(s) or image(s)
128if iscalar~=0
129    [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima
130    if iscalar==1 && ~isempty(DataOut_1) % case for which only the second field is a scalar
131         DataOut_1.A=A{1};
132         DataOut_1.AX=AX;
133         DataOut_1.AY=AY;
134    else
135        DataOut.A=A{1};
136        DataOut.AX=AX;
137        DataOut.AY=AY;
138    end
139    if iscalar==2
140        DataOut_1.A=A{2};
141        DataOut_1.AX=AX;
142        DataOut_1.AY=AY;
143    end
144end
145
146% subtract fields
147if ~isempty(DataOut_1)
148DataOut=sub_field(DataOut,[],DataOut_1);
149end
150%------------------------------------------------
151%--- transform a single field
152function DataOut=phys_1(Data,Calib,ZIndex)
153%------------------------------------------------
154%% set default output
155DataOut=Data;%default
156DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters
157
158%% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
159if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
160  [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
161    Dt=1; %default
162    if isfield(Data,'dt')&&~isempty(Data.dt)
163        Dt=Data.dt;
164    end
165    if isfield(Data,'Dt')&&~isempty(Data.Dt)
166        Dt=Data.Dt;
167    end
168    if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
169        [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
170        [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
171        DataOut.U=(XOut_2-XOut_1)/Dt;
172        DataOut.V=(YOut_2-YOut_1)/Dt;
173    end
174%     if ~strcmp(Calib.CalibrationType,'rescale') && isfield(Data,'X_tps') && isfield(Data,'Y_tps')
175%         [DataOut.X_tps,DataOut.Y_tps]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
176%     end
177end
178
179%% suppress tps
180list_tps={'Coord_tps'  'U_tps'  'V_tps'  'SubRange'  'NbSites'};
181ind_remove=[];
182for ilist=1:numel(list_tps)
183    ind_tps=find(strcmp(list_tps{ilist},Data.ListVarName));
184    if ~isempty(ind_tps)
185        ind_remove=[ind_remove ind_tps];
186        DataOut=rmfield(DataOut,list_tps{ilist});
187    end
188end
189if isfield(DataOut,'VarAttribute') && numel(DataOut.VarAttribute)>=3 && isfield(DataOut.VarAttribute{3},'VarIndex_tps')
190    DataOut.VarAttribute{3}=rmfield(DataOut.VarAttribute{3},'VarIndex_tps');
191end
192if isfield(DataOut,'VarAttribute')&& numel(DataOut.VarAttribute)>=4 && isfield(DataOut.VarAttribute{4},'VarIndex_tps')
193    DataOut.VarAttribute{4}=rmfield(DataOut.VarAttribute{4},'VarIndex_tps');
194end
195if ~isempty(ind_remove)
196    DataOut.ListVarName(ind_remove)=[];
197    DataOut.VarDimName(ind_remove)=[];
198    DataOut.VarAttribute(ind_remove)=[];
199end
200   
201%% transform of spatial derivatives: TODO check the case with plane angles
202if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)
203    % estimate the Jacobian matrix DXpx/DXphys
204    for ip=1:length(Data.X)
205        [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
206        [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
207        [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
208        [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
209        %Jacobian matrix DXpphys/DXpx
210        DjXi(1,1)=(Xp1-Xm1);
211        DjXi(2,1)=(Yp1-Ym1);
212        DjXi(1,2)=(Xp2-Xm2);
213        DjXi(2,2)=(Yp2-Ym2);
214        DjUi(:,:)=Data.DjUi(ip,:,:);
215        DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
216        DataOut.DjUi(ip,:,:)=DjUi';
217    end
218    DataOut.DjUi =  DataOut.DjUi/Dt;   %     min(Data.DjUi(:,1,1))=DUDX
219end
220
221
222%%%%%%%%%%%%%%%%%%%%
223function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)
224xcorner=[];
225ycorner=[];
226npx=[];
227npy=[];
228dx=ones(1,numel(A));
229dy=ones(1,numel(A));
230for icell=1:numel(A)
231    siz=size(A{icell});
232    npx=[npx siz(2)];
233    npy=[npy siz(1)];
234    Calib=CalibIn{icell};
235    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
236    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
237    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
238    dx(icell)=(max(xcorner_new)-min(xcorner_new))/(siz(2)-1);
239    dy(icell)=(max(ycorner_new)-min(ycorner_new))/(siz(1)-1);
240    xcorner=[xcorner xcorner_new];
241    ycorner=[ycorner ycorner_new];
242end
243Rangx(1)=min(xcorner);
244Rangx(2)=max(xcorner);
245Rangy(2)=min(ycorner);
246Rangy(1)=max(ycorner);
247test_multi=(max(npx)~=min(npx)) || (max(npy)~=min(npy)); %different image lengths
248npX=1+round((Rangx(2)-Rangx(1))/min(dx));% nbre of pixels in the new image (use the finest resolution min(dx) in the set of images)
249npY=1+round((Rangy(1)-Rangy(2))/min(dy));
250x=linspace(Rangx(1),Rangx(2),npX);
251y=linspace(Rangy(1),Rangy(2),npY);
252[X,Y]=meshgrid(x,y);%grid in physical coordiantes
253%vec_B=[];
254A_out=cell(1,numel(A));
255for icell=1:length(A)
256    Calib=CalibIn{icell};
257    % rescaling of the image coordinates without change of the image array
258    if strcmp(Calib.CalibrationType,'rescale') && isequal(Calib,CalibIn{1})
259        A_out{icell}=A{icell};%no transform
260        Rangx=[0.5 npx-0.5];%image coordiantes of corners
261        Rangy=[npy-0.5 0.5];
262        [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
263        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
264    else         
265        % the image needs to be interpolated to the new coordinates
266        zphys=0; %default
267        if isfield(Calib,'SliceCoord') %.Z= index of plane
268           SliceCoord=Calib.SliceCoord(ZIndex,:);
269           zphys=SliceCoord(3); %to generalize for non-parallel planes
270           if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
271                H=Calib.InterfaceCoord(3);
272                if H>zphys
273                    zphys=H-(H-zphys)/Calib.RefractionIndex; %corrected z (virtual object)
274                end
275           end
276        end
277        xima=[0.5:npx-0.5];%image coordiantes of corners
278        yima=npy-0.5:-1:0.5;
279        [XIMA_init,YIMA_init]=meshgrid(xima,yima);%grid of initial image in px coordiantes
280        [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);% image coordinates for each point in the real
281        %[XPHYS_init,YPHYS_init]=phys_XYZ(Calib,XIMA_init,YIMA_init,ZIndex);
282                 testuint8=isa(A{icell},'uint8');
283         testuint16=isa(A{icell},'uint16');
284        if ndims(A{icell})==2 %(B/W images)
285        A_out{icell}=interp2(XIMA_init,YIMA_init,double(A{icell}),XIMA,YIMA);
286%         [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
287%         [XIMA_init,YIMA_init]=px_XYZ(CalibIn{icell},X,Y,zphys);% image coordinates for each point in the real space grid
288%         
289%         XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
290%         YIMA=reshape(round(YIMA),1,npX*npY);
291%         flagin=XIMA>=1 & XIMA<=npx(icell) & YIMA >=1 & YIMA<=npy(icell);%flagin=1 inside the original image
292
293%         if numel(siz)==2 %(B/W images)
294%             vec_A=reshape(A{icell},1,npx(icell)*npy(icell));%put the original image in line
295%             %ind_in=find(flagin);
296%             ind_out=find(~flagin);
297%             ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
298%             ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
299%             %vec_B(ind_in)=vec_A(ICOMB);
300%             vec_B(flagin)=vec_A(ICOMB);
301%             vec_B(~flagin)=zeros(size(ind_out));
302% %             vec_B(ind_out)=zeros(size(ind_out));
303%             A_out{icell}=reshape(vec_B,npY,npX);%new image in real coordinates
304         elseif ndims(A{icell})==3     
305             for icolor=1:size(A{icell},3)
306                 A{icell}=double(A{icell});
307                 A_out{icell}(:,:,icolor)=interp2(XIMA_init,YIMA_init,A{icell}(:,:,icolor),XIMA,YIMA);
308%                 vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
309%                % ind_in=find(flagin);
310%                 ind_out=find(~flagin);
311%                 ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
312%                 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
313%                 vec_B(flagin)=vec_A(ICOMB);
314%                 vec_B(~flagin)=zeros(size(ind_out));
315%                 A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
316             end
317         end
318        if testuint8
319            A_out{icell}=uint8(A_out{icell});
320        end
321        if testuint16
322            A_out{icell}=uint16(A_out{icell});
323        end     
324    end
325end
326
Note: See TracBrowser for help on using the repository browser.