Description of phys

0001 %'phys': transforms image (px) to real world (phys) coordinates using geometric calibration parameters
0002 % OUTPUT:
0003 % DataOut:   structure of modified data (like UvData)
0004 %
0005 %INPUT:
0006 % Data:  structure of input data (like UvData)
0007 %       with fields A (image or scalar matrix), AX, AY
0008 %       X,Y,U,V, DjUi
0009 % Data.CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px'
0010 % Calib: structure containing calibration parameters or a subtree Calib.GeometryCalib =calibration data (tsai parameters)
0011 
0012 function [DataOut,DataOut_1]=phys(varargin)
0013 % A FAIRE: 1- verifier si DataIn est une 'field structure'(.ListVarName'):
0014 % chercher ListVarAttribute, for each field (cell of variables):
0015 %   .CoordType: 'phys' or 'px'   (default==phys, no transform)
0016 %   .scale_factor: =dt (to transform displacement into velocity) default=1
0017 %   .covariance: 'scalar', 'coord', 'D_i': covariant (like velocity), 'D^i': contravariant (like gradient), 'D^jD_i' (like strain tensor)
0018 %   (default='coord' if .Role='coord_x,_y...,
0019 %            'D_i' if '.Role='vector_x,...',
0020 %              'scalar', else (thenno change except scale factor)
0021 Calib{1}=[];
0022 if nargin==2||nargin==4 % nargin =nbre of input variables
0023     Data=varargin{1};
0024     DataOut=Data;%default
0025     DataOut_1=[];%default
0026     CalibData=varargin{2};
0027     if isfield(CalibData,'GeometryCalib')
0028         Calib{1}=CalibData.GeometryCalib;
0029     end
0030     Calib{2}=Calib{1};
0031 else
0032     DataOut.Txt='wrong input: need two or four structures';
0033 end
0034 test_1=0;
0035 if nargin==4
0036     test_1=1;
0037     Data_1=varargin{3};
0038     DataOut_1=Data_1;%default
0039     CalibData_1=varargin{4};
0040     if isfield(CalibData_1,'GeometryCalib')
0041         Calib{2}=CalibData_1.GeometryCalib;
0042     end
0043 end
0044 iscalar=0;
0045 if  ~isempty(Calib{1})
0046     DataOut=phys_1(Data,Calib{1});
0047     %case of images or scalar
0048     if isfield(Data,'A') && isfield(Data,'AX') && ~isempty(Data.AX) && isfield(Data,'AY')&&...
0049                                            ~isempty(Data.AY) && length(Data.A)>1
0050         iscalar=1;
0051         A{1}=Data.A;
0052     end
0053 end
0054 %transform of X,Y coordinates for vector fields
0055 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)
0056     ZIndex=Data.ZIndex;
0057 else
0058     ZIndex=0;
0059 end
0060 if test_1
0061     DataOut_1=phys_1(Data_1,Calib{2});
0062     if isfield(Data_1,'A')&&isfield(Data_1,'AX')&&~isempty(Data_1.AX) && isfield(Data_1,'AY')&&...
0063                                        ~isempty(Data_1.AY)&&length(Data_1.A)>1
0064           iscalar=iscalar+1;
0065           Calib{iscalar}=Calib{2};
0066           A{iscalar}=Data_1.A;
0067           if isfield(Data_1,'ZIndex') && ~isequal(Data_1.ZIndex,ZIndex)
0068               DataOut.Txt='inconsistent plane indexes in the two input fields';
0069           end
0070           if iscalar==1% case for which only the second field is a scalar
0071                [A,AX,AY]=phys_Ima(A,Calib,ZIndex);
0072                DataOut_1.A=A{1};
0073                DataOut_1.AX=AX; 
0074                DataOut_1.AY=AY;
0075                return
0076           end
0077     end
0078 end
0079 if iscalar~=0
0080     [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima
0081     DataOut.A=A{1};
0082     DataOut.AX=AX; 
0083     DataOut.AY=AY;
0084     if iscalar==2
0085         DataOut_1.A=A{2};
0086         DataOut_1.AX=AX; 
0087         DataOut_1.AY=AY;
0088     end
0089 end
0090 
0091 %------------------------------------------------
0092 function DataOut=phys_1(Data,Calib)
0093 % for icell=1:length(Data)
0094 
0095 DataOut=Data;%default
0096 DataOut.CoordType='phys'; %put flag for physical coordinates
0097 % The transform ACTS ONLY IF .CoordType='px'and Calib defined
0098 if isfield(Data,'CoordType')&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
0099     if isfield(Calib,'CoordUnit')
0100         DataOut.CoordUnit=Calib.CoordUnit;
0101     else
0102         DataOut.CoordUnit='cm'; %default
0103 %     elseif isfield(DataOut,'CoordUnit')
0104 %         DataOut=rmfield(DataOut,'CoordUnit');
0105     end
0106     DataOut.TimeUnit='s';
0107     %transform of X,Y coordinates for vector fields
0108     if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)
0109         Z=Data.ZIndex;
0110     else
0111         Z=0;
0112     end
0113     if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
0114         [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); 
0115         if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt') 
0116             if ~isempty(Data.dt)
0117             [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
0118             [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
0119             DataOut.U=(XOut_2-XOut_1)/Data.dt;
0120             DataOut.V=(YOut_2-YOut_1)/Data.dt;
0121             end
0122         end
0123     end
0124     %transform of an image or scalar: done in phys_ima
0125       
0126     %transform of spatial derivatives
0127     if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
0128           && isfield(Data,'dt')    
0129         if ~isempty(Data.dt)
0130             % estimate the Jacobian matrix DXpx/DXphys
0131             for ip=1:length(Data.X) 
0132                 [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
0133                 [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
0134                 [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
0135                 [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z); 
0136             %Jacobian matrix DXpphys/DXpx
0137                DjXi(1,1)=(Xp1-Xm1);
0138                DjXi(2,1)=(Yp1-Ym1);
0139                DjXi(1,2)=(Xp2-Xm2);
0140                DjXi(2,2)=(Yp2-Ym2);
0141                DjUi(:,:)=Data.DjUi(ip,:,:);
0142                DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
0143                DataOut.DjUi(ip,:,:)=DjUi';
0144             end
0145             DataOut.DjUi =  DataOut.DjUi/Data.dt;   %     min(Data.DjUi(:,1,1))=DUDX
0146         end
0147     end
0148 end
0149 
0150 
0151 %%%%%%%%%%%%%%%%%%%%
0152 function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)
0153 xcorner=[];
0154 ycorner=[];
0155 npx=[];
0156 npy=[];
0157 for icell=1:length(A)
0158     siz=size(A{icell});
0159     npx=[npx siz(2)];
0160     npy=[npy siz(1)];
0161     Calib=CalibIn{icell};
0162     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
0163     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
0164     [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
0165     xcorner=[xcorner xcorner_new];
0166     ycorner=[ycorner ycorner_new];
0167 end
0168 Rangx(1)=min(xcorner);
0169 Rangx(2)=max(xcorner);
0170 Rangy(2)=min(ycorner);
0171 Rangy(1)=max(ycorner);
0172 test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy)); 
0173 npx=max(npx);
0174 npy=max(npy);
0175 x=linspace(Rangx(1),Rangx(2),npx);
0176 y=linspace(Rangy(1),Rangy(2),npy);
0177 [X,Y]=meshgrid(x,y);%grid in physical coordiantes
0178 vec_B=[];
0179 A_out={};
0180 for icell=1:length(A) 
0181     Calib=CalibIn{icell};
0182     if (isfield(Calib,'R') && ~isequal(Calib.R(2,1),0) && ~isequal(Calib.R(1,2),0)) ||...
0183         ((isfield(Calib,'kappa1')&& ~isequal(Calib.kappa1,0))) || test_multi || ~isequal(Calib,CalibIn{1})
0184         zphys=0; %default
0185         if isfield(Calib,'SliceCoord') %.Z= index of plane
0186            SliceCoord=Calib.SliceCoord(ZIndex,:);
0187            zphys=SliceCoord(3); %to generalize for non-parallel planes
0188         end
0189         [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
0190         XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
0191         YIMA=reshape(round(YIMA),1,npx*npy);
0192         flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
0193         if numel(siz)==2 %(B/W images)
0194             vec_A=reshape(A{icell},1,npx*npy);%put the original image in line
0195             ind_in=find(flagin);
0196             ind_out=find(~flagin);
0197             ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
0198             ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
0199             vec_B(ind_in)=vec_A(ICOMB);
0200             vec_B(ind_out)=zeros(size(ind_out));
0201             A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
0202         elseif numel(siz)==3     
0203             for icolor=1:siz(3)
0204                 vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
0205                 ind_in=find(flagin);
0206                 ind_out=find(~flagin);
0207                 ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
0208                 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
0209                 vec_B(ind_in)=vec_A(ICOMB);
0210                 vec_B(ind_out)=zeros(size(ind_out));
0211                 A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
0212             end
0213         end
0214     else%
0215         
0216         A_out{icell}=A{icell};%no transform
0217         Rangx=[0.5 npx-0.5];%image coordiantes of corners
0218         Rangy=[npy-0.5 0.5];
0219         [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[ZIndex ZIndex]);%case of translations without rotation and quadratic deformation
0220         [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[ZIndex ZIndex]);
0221     end
0222 end
phys

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
