source: trunk/src/phys.m @ 36

Last change on this file since 36 was 19, checked in by gostiaux, 14 years ago

the private files have been moved down to the root folder

File size: 9.3 KB
Line 
1%'phys': transforms image (px) to real world (phys) coordinates using geometric calibration parameters
2% OUTPUT:
3% DataOut:   structure representing the modified field
4% DataOut_1: structure representing the second modified field
5
6%INPUT:
7% Data:  structure of input data
8%       with fields .A (image or scalar matrix), AX, AY
9%       .X,.Y,.U,.V, .DjUi
10%       .ZIndex: index of plane in multilevel case
11% Data.CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px'
12% Calib: structure containing calibration parameters or a subtree Calib.GeometryCalib =calibration data (tsai parameters)
13
14function [DataOut,DataOut_1]=phys(varargin)
15% A FAIRE: 1- verifier si DataIn est une 'field structure'(.ListVarName'):
16% chercher ListVarAttribute, for each field (cell of variables):
17%   .CoordType: 'phys' or 'px'   (default==phys, no transform)
18%   .scale_factor: =dt (to transform displacement into velocity) default=1
19%   .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)
20%   (default='coord' if .Role='coord_x,_y...,
21%            'D_i' if '.Role='vector_x,...',
22%              'scalar', else (thenno change except scale factor)
23Calib{1}=[];
24if nargin==2||nargin==4 % nargin =nbre of input variables
25    Data=varargin{1};
26    DataOut=Data;%default
27    DataOut_1=[];%default
28    CalibData=varargin{2};
29    if isfield(CalibData,'GeometryCalib')
30        Calib{1}=CalibData.GeometryCalib;
31    end
32    Calib{2}=Calib{1};
33else
34    DataOut.Txt='wrong input: need two or four structures';
35end
36test_1=0;
37if nargin==4
38    test_1=1;
39    Data_1=varargin{3};
40    DataOut_1=Data_1;%default
41    CalibData_1=varargin{4};
42    if isfield(CalibData_1,'GeometryCalib')
43        Calib{2}=CalibData_1.GeometryCalib;
44    end
45end
46iscalar=0;
47if  ~isempty(Calib{1})
48    DataOut=phys_1(Data,Calib{1});
49    %case of images or scalar: in case of two input fields, we need to project the transform of on the same regular grid
50    if isfield(Data,'A') && isfield(Data,'AX') && ~isempty(Data.AX) && isfield(Data,'AY')&&...
51                                           ~isempty(Data.AY) && length(Data.A)>1
52        iscalar=1;
53        A{1}=Data.A;
54    end
55end
56%transform of X,Y coordinates for vector fields
57if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)
58    ZIndex=Data.ZIndex;
59else
60    ZIndex=0;
61end
62if test_1
63    DataOut_1=phys_1(Data_1,Calib{2});
64    if isfield(Data_1,'A')&&isfield(Data_1,'AX')&&~isempty(Data_1.AX) && isfield(Data_1,'AY')&&...
65                                       ~isempty(Data_1.AY)&&length(Data_1.A)>1
66          iscalar=iscalar+1;
67          Calib{iscalar}=Calib{2};
68          A{iscalar}=Data_1.A;
69          if isfield(Data_1,'ZIndex') && ~isequal(Data_1.ZIndex,ZIndex)
70              DataOut.Txt='inconsistent plane indexes in the two input fields';
71          end
72          if iscalar==1% case for which only the second field is a scalar
73               [A,AX,AY]=phys_Ima(A,Calib,ZIndex);
74               DataOut_1.A=A{1};
75               DataOut_1.AX=AX;
76               DataOut_1.AY=AY;
77               return
78          end
79    end
80end
81if iscalar~=0
82    [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima
83    DataOut.A=A{1};
84    DataOut.AX=AX;
85    DataOut.AY=AY;
86    if iscalar==2
87        DataOut_1.A=A{2};
88        DataOut_1.AX=AX;
89        DataOut_1.AY=AY;
90    end
91end
92
93%------------------------------------------------
94function DataOut=phys_1(Data,Calib)
95% for icell=1:length(Data)
96
97DataOut=Data;%default
98DataOut.CoordType='phys'; %put flag for physical coordinates
99% The transform ACTS ONLY IF .CoordType='px'and Calib defined
100if isfield(Data,'CoordType')&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
101    if isfield(Calib,'CoordUnit')
102        DataOut.CoordUnit=Calib.CoordUnit;
103    else
104        DataOut.CoordUnit='cm'; %default
105%     elseif isfield(DataOut,'CoordUnit')
106%         DataOut=rmfield(DataOut,'CoordUnit');
107    end
108    DataOut.TimeUnit='s';
109    %transform of X,Y coordinates for vector fields
110    if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)
111        Z=Data.ZIndex;
112    else
113        Z=0;
114    end
115    if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
116        [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z);
117        if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt')
118            if ~isempty(Data.dt)
119            [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
120            [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
121            DataOut.U=(XOut_2-XOut_1)/Data.dt;
122            DataOut.V=(YOut_2-YOut_1)/Data.dt;
123            end
124        end
125    end
126    %transform of an image or scalar: done in phys_ima
127     
128    %transform of spatial derivatives
129    if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
130          && isfield(Data,'dt')   
131        if ~isempty(Data.dt)
132            % estimate the Jacobian matrix DXpx/DXphys
133            for ip=1:length(Data.X)
134                [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
135                [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
136                [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
137                [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
138            %Jacobian matrix DXpphys/DXpx
139               DjXi(1,1)=(Xp1-Xm1);
140               DjXi(2,1)=(Yp1-Ym1);
141               DjXi(1,2)=(Xp2-Xm2);
142               DjXi(2,2)=(Yp2-Ym2);
143               DjUi(:,:)=Data.DjUi(ip,:,:);
144               DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
145               DataOut.DjUi(ip,:,:)=DjUi';
146            end
147            DataOut.DjUi =  DataOut.DjUi/Data.dt;   %     min(Data.DjUi(:,1,1))=DUDX                         
148        end
149    end
150end
151
152
153%%%%%%%%%%%%%%%%%%%%
154function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)
155xcorner=[];
156ycorner=[];
157npx=[];
158npy=[];
159for icell=1:length(A)
160    siz=size(A{icell});
161    npx=[npx siz(2)];
162    npy=[npy siz(1)];
163    Calib=CalibIn{icell};
164    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
165    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
166    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
167    xcorner=[xcorner xcorner_new];
168    ycorner=[ycorner ycorner_new];
169end
170Rangx(1)=min(xcorner);
171Rangx(2)=max(xcorner);
172Rangy(2)=min(ycorner);
173Rangy(1)=max(ycorner);
174test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
175npx=max(npx);
176npy=max(npy);
177x=linspace(Rangx(1),Rangx(2),npx);
178y=linspace(Rangy(1),Rangy(2),npy);
179[X,Y]=meshgrid(x,y);%grid in physical coordiantes
180vec_B=[];
181A_out={};
182for icell=1:length(A)
183    Calib=CalibIn{icell};
184    if (isfield(Calib,'R') && ~isequal(Calib.R(2,1),0) && ~isequal(Calib.R(1,2),0)) ||...
185        ((isfield(Calib,'kappa1')&& ~isequal(Calib.kappa1,0))) || test_multi || ~isequal(Calib,CalibIn{1})
186        zphys=0; %default
187        if isfield(Calib,'SliceCoord') %.Z= index of plane
188           SliceCoord=Calib.SliceCoord(ZIndex,:);
189           zphys=SliceCoord(3); %to generalize for non-parallel planes
190        end
191        [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
192        XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
193        YIMA=reshape(round(YIMA),1,npx*npy);
194        flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
195        testuint8=isa(A{icell},'uint8');
196        testuint16=isa(A{icell},'uint16');
197        if numel(siz)==2 %(B/W images)
198            vec_A=reshape(A{icell},1,npx*npy);%put the original image in line
199            ind_in=find(flagin);
200            ind_out=find(~flagin);
201            ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
202            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
203            vec_B(ind_in)=vec_A(ICOMB);
204            vec_B(ind_out)=zeros(size(ind_out));
205            A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
206        elseif numel(siz)==3     
207            for icolor=1:siz(3)
208                vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
209                ind_in=find(flagin);
210                ind_out=find(~flagin);
211                ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
212                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
213                vec_B(ind_in)=vec_A(ICOMB);
214                vec_B(ind_out)=zeros(size(ind_out));
215                A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
216            end
217        end
218        if testuint8
219            A_out{icell}=uint8(A_out{icell});
220        end
221        if testuint16
222            A_out{icell}=uint16(A_out{icell});
223        end
224    else%
225       
226        A_out{icell}=A{icell};%no transform
227        Rangx=[0.5 npx-0.5];%image coordiantes of corners
228        Rangy=[npy-0.5 0.5];
229        [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[ZIndex ZIndex]);%case of translations without rotation and quadratic deformation
230        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[ZIndex ZIndex]);
231    end
232end
Note: See TracBrowser for help on using the repository browser.