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

Last change on this file since 40 was 40, checked in by sommeria, 11 years ago

-relabel_i_j added to the svn (relabel an image series with two indices)
-Menu projection Object added to uvmat, rationalizes the call of projection objects (to improve yet)
-get_plot_handles: modified in relationwith the Menu projection object
-mouse_down: corrected in relation with previous change in field transforms
-transform_field functions added to the svn
-set_object and keyboard_callback: cleaning

File size: 11.6 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 modified field
7% DataOut_1: structure representing the second modified field
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
17function [DataOut,DataOut_1]=phys(varargin)
18% A FAIRE: 1- verifier si DataIn est une 'field structure'(.ListVarName'):
19% chercher ListVarAttribute, for each field (cell of variables):
20%   .CoordType: 'phys' or 'px'   (default==phys, no transform)
21%   .scale_factor: =dt (to transform displacement into velocity) default=1
22%   .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)
23%   (default='coord' if .Role='coord_x,_y...,
24%            'D_i' if '.Role='vector_x,...',
25%              'scalar', else (thenno change except scale factor)
26Calib{1}=[];
27if nargin==2||nargin==4 % nargin =nbre of input variables
28    Data=varargin{1};
29    DataOut=Data;%default
30    DataOut_1=[];%default
31    CalibData=varargin{2};
32    if isfield(CalibData,'GeometryCalib')
33        Calib{1}=CalibData.GeometryCalib;
34    end
35    Calib{2}=Calib{1};
36else
37    DataOut.Txt='wrong input: need two or four structures';
38end
39test_1=0;
40if nargin==4
41    test_1=1;
42    Data_1=varargin{3};
43    DataOut_1=Data_1;%default
44    CalibData_1=varargin{4};
45    if isfield(CalibData_1,'GeometryCalib')
46        Calib{2}=CalibData_1.GeometryCalib;
47    end
48end
49iscalar=0;
50if  ~isempty(Calib{1})
51    DataOut=phys_1(Data,Calib{1});
52    %case of images or scalar: in case of two input fields, we need to project the transform of on the same regular grid
53    if isfield(Data,'A') && isfield(Data,'AX') && ~isempty(Data.AX) && isfield(Data,'AY')&&...
54                                           ~isempty(Data.AY) && length(Data.A)>1
55        iscalar=1;
56        A{1}=Data.A;
57    end
58end
59%transform of X,Y coordinates for vector fields
60if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)
61    ZIndex=Data.ZIndex;
62else
63    ZIndex=0;
64end
65if test_1
66    DataOut_1=phys_1(Data_1,Calib{2});
67    if isfield(Data_1,'A')&&isfield(Data_1,'AX')&&~isempty(Data_1.AX) && isfield(Data_1,'AY')&&...
68                                       ~isempty(Data_1.AY)&&length(Data_1.A)>1
69          iscalar=iscalar+1;
70          Calib{iscalar}=Calib{2};
71          A{iscalar}=Data_1.A;
72          if isfield(Data_1,'ZIndex') && ~isequal(Data_1.ZIndex,ZIndex)
73              DataOut.Txt='inconsistent plane indexes in the two input fields';
74          end
75          if iscalar==1% case for which only the second field is a scalar
76               [A,AX,AY]=phys_Ima(A,Calib,ZIndex);
77               DataOut_1.A=A{1};
78               DataOut_1.AX=AX;
79               DataOut_1.AY=AY;
80               return
81          end
82    end
83end
84if iscalar~=0
85    [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima
86    DataOut.A=A{1};
87    DataOut.AX=AX;
88    DataOut.AY=AY;
89    if iscalar==2
90        DataOut_1.A=A{2};
91        DataOut_1.AX=AX;
92        DataOut_1.AY=AY;
93    end
94end
95
96%------------------------------------------------
97function DataOut=phys_1(Data,Calib)
98% for icell=1:length(Data)
99
100DataOut=Data;%default
101DataOut.CoordType='phys'; %put flag for physical coordinates
102% The transform ACTS ONLY IF .CoordType='px'and Calib defined
103if isfield(Data,'CoordType')&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
104    if isfield(Calib,'CoordUnit')
105        DataOut.CoordUnit=Calib.CoordUnit;
106    else
107        DataOut.CoordUnit='cm'; %default
108%     elseif isfield(DataOut,'CoordUnit')
109%         DataOut=rmfield(DataOut,'CoordUnit');
110    end
111    DataOut.TimeUnit='s';
112    %transform of X,Y coordinates for vector fields
113    if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)
114        Z=Data.ZIndex;
115    else
116        Z=0;
117    end
118    if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
119        [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z);
120        if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt')
121            if ~isempty(Data.dt)
122            [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
123            [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
124            DataOut.U=(XOut_2-XOut_1)/Data.dt;
125            DataOut.V=(YOut_2-YOut_1)/Data.dt;
126            end
127        end
128    end
129    %transform of an image or scalar: done in phys_ima
130     
131    %transform of spatial derivatives
132    if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
133          && isfield(Data,'dt')   
134        if ~isempty(Data.dt)
135            % estimate the Jacobian matrix DXpx/DXphys
136            for ip=1:length(Data.X)
137                [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
138                [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
139                [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
140                [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
141            %Jacobian matrix DXpphys/DXpx
142               DjXi(1,1)=(Xp1-Xm1);
143               DjXi(2,1)=(Yp1-Ym1);
144               DjXi(1,2)=(Xp2-Xm2);
145               DjXi(2,2)=(Yp2-Ym2);
146               DjUi(:,:)=Data.DjUi(ip,:,:);
147               DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
148               DataOut.DjUi(ip,:,:)=DjUi';
149            end
150            DataOut.DjUi =  DataOut.DjUi/Data.dt;   %     min(Data.DjUi(:,1,1))=DUDX                         
151        end
152    end
153end
154
155
156%%%%%%%%%%%%%%%%%%%%
157function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)
158xcorner=[];
159ycorner=[];
160npx=[];
161npy=[];
162for icell=1:length(A)
163    siz=size(A{icell});
164    npx=[npx siz(2)];
165    npy=[npy siz(1)];
166    Calib=CalibIn{icell};
167    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
168    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
169    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
170    xcorner=[xcorner xcorner_new];
171    ycorner=[ycorner ycorner_new];
172end
173Rangx(1)=min(xcorner);
174Rangx(2)=max(xcorner);
175Rangy(2)=min(ycorner);
176Rangy(1)=max(ycorner);
177test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
178npx=max(npx);
179npy=max(npy);
180x=linspace(Rangx(1),Rangx(2),npx);
181y=linspace(Rangy(1),Rangy(2),npy);
182[X,Y]=meshgrid(x,y);%grid in physical coordiantes
183vec_B=[];
184A_out={};
185for icell=1:length(A)
186    Calib=CalibIn{icell};
187    if (isfield(Calib,'R') && ~isequal(Calib.R(2,1),0) && ~isequal(Calib.R(1,2),0)) ||...
188        ((isfield(Calib,'kappa1')&& ~isequal(Calib.kappa1,0))) || test_multi || ~isequal(Calib,CalibIn{1})
189        zphys=0; %default
190        if isfield(Calib,'SliceCoord') %.Z= index of plane
191           SliceCoord=Calib.SliceCoord(ZIndex,:);
192           zphys=SliceCoord(3); %to generalize for non-parallel planes
193        end
194        [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
195        XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
196        YIMA=reshape(round(YIMA),1,npx*npy);
197        flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
198        testuint8=isa(A{icell},'uint8');
199        testuint16=isa(A{icell},'uint16');
200        if numel(siz)==2 %(B/W images)
201            vec_A=reshape(A{icell},1,npx*npy);%put the original image in line
202            ind_in=find(flagin);
203            ind_out=find(~flagin);
204            ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
205            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
206            vec_B(ind_in)=vec_A(ICOMB);
207            vec_B(ind_out)=zeros(size(ind_out));
208            A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
209        elseif numel(siz)==3     
210            for icolor=1:siz(3)
211                vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
212                ind_in=find(flagin);
213                ind_out=find(~flagin);
214                ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
215                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
216                vec_B(ind_in)=vec_A(ICOMB);
217                vec_B(ind_out)=zeros(size(ind_out));
218                A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
219            end
220        end
221        if testuint8
222            A_out{icell}=uint8(A_out{icell});
223        end
224        if testuint16
225            A_out{icell}=uint16(A_out{icell});
226        end
227    else%
228       
229        A_out{icell}=A{icell};%no transform
230        Rangx=[0.5 npx-0.5];%image coordiantes of corners
231        Rangy=[npy-0.5 0.5];
232        [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[ZIndex ZIndex]);%case of translations without rotation and quadratic deformation
233        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[ZIndex ZIndex]);
234    end
235end
236
237%'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
238% function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)
239%
240%OUTPUT:
241%
242%INPUT:
243%Z: index of plane
244function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
245if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
246    Zindex=Z;
247    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
248else
249%     if exist('Z','var')
250%         Zphys=Z;
251%     else
252        Zphys=0;
253%     end
254end
255if ~exist('X','var')||~exist('Y','var')
256    Xphys=[];
257    Yphys=[];%default
258    return
259end
260Xphys=X;%default
261Yphys=Y;
262%image transform
263if isfield(Calib,'R')
264    R=(Calib.R)';
265    Dx=R(5)*R(7)-R(4)*R(8);
266    Dy=R(1)*R(8)-R(2)*R(7);
267    D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
268    Z11=R(6)*R(8)-R(5)*R(9);
269    Z12=R(2)*R(9)-R(3)*R(8); 
270    Z21=R(4)*R(9)-R(6)*R(7);
271    Z22=R(3)*R(7)-R(1)*R(9);
272    Zx0=R(3)*R(5)-R(2)*R(6);
273    Zy0=R(1)*R(6)-R(3)*R(4);
274    A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
275    A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
276    A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
277    A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
278    X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
279    Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
280        %px to camera:
281    Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
282    Yd=Calib.dpy*(Y-Calib.Cy);
283    dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
284    Xu=dist_fact.*Xd;%undistorted sensor coordinates
285    Yu=dist_fact.*Yd;
286    denom=Dx*Xu+Dy*Yu+D0;
287    % denom2=denom.*denom;
288    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
289    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
290end
291
292%'px_XYZ': transform phys coordinates to image coordinates (px)
293%
294% OUPUT:
295% X,Y: array of coordinates in the image cooresponding to the input physical positions
296%                    (origin at lower leftcorner, unit=pixel)
297
298% INPUT:
299% Calib: structure containing the calibration parameters (read from the ImaDoc .xml file)
300% Xphys, Yphys: array of x,y physical coordinates
301% [Zphys]: corresponding array of z physical coordinates (0 by default)
302
303
304
305
Note: See TracBrowser for help on using the repository browser.