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

Last change on this file since 494 was 494, checked in by sommeria, 12 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

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