Home > . > phys.m

# phys

## PURPOSE

'phys': transforms image (px) to real world (phys) coordinates using geometric calibration parameters

## SYNOPSIS

function [DataOut,DataOut_1]=phys(varargin)

## DESCRIPTION

    'phys': transforms image (px) to real world (phys) coordinates using geometric calibration parameters
     OUTPUT:
     DataOut:   structure of modified data (like UvData)

    INPUT:
     Data:  structure of input data (like UvData)
           with fields A (image or scalar matrix), AX, AY
           X,Y,U,V, DjUi
     Data.CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px'
     Calib: structure containing calibration parameters or a subtree Calib.GeometryCalib =calibration data (tsai parameters)

## CROSS-REFERENCE INFORMATION

This function calls:

- phys_XYZ 'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
- px_XYZ 'px_XYZ': transform phys coordinates to image coordinates (px)

This function is called by:

## SUBFUNCTIONS

- function DataOut=phys_1(Data,Calib)
- function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)

## SOURCE CODE

    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