Ignore:
Timestamp:
Apr 25, 2012, 2:16:56 PM (12 years ago)
Author:
sommeria
Message:

bugs corrected and cleaning in phys

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/transform_field/phys.m

    r250 r396  
    44
    55% OUTPUT:
    6 % DataOut:   structure representing the modified field
    7 % DataOut_1: structure representing the second modified field
     6% DataOut:   structure representing the first field in phys coordinates
     7% DataOut_1: structure representing the second  field in phys coordinates
    88
    99%INPUT:
    … …  
    1414%       .CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px'
    1515% 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.
    1617
    1718function [DataOut,DataOut_1]=phys(varargin)
    … …  
    2425%            'D_i' if '.Role='vector_x,...',
    2526%              'scalar', else (thenno change except scale factor)
    26 Calib{1}=[];
    27 if 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};
    36 else
    37     DataOut.Txt='wrong input: need two or four structures';
    38 end
    39 test_1=0;
    40 if 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
    48 end
    49 iscalar=0;
    50 if  ~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
    58 end
    59 %transform of X,Y coordinates for vector fields
    60 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
    61     ZIndex=Data.ZIndex;
     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}=CalibData_1.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;
    6249else
    6350    ZIndex=1;
    6451end
    65 if 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
    83 end
     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(Data_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)
    84119if iscalar~=0
    85120    [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;
     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
    89130    if iscalar==2
    90131        DataOut_1.A=A{2};
    … …  
    94135end
    95136
    96 % DataOut.VarDimName{2}
    97 % DataOut.VarDimName{3}
    98 % DataOut.VarDimName{4}
    99 % DataOut.VarDimName{5}
    100 % DataOut.VarDimName{6}
    101 % DataOut.VarDimName{7}
    102 % DataOut.VarAttribute{1}
    103 % DataOut.VarAttribute{2}
    104 % DataOut.VarAttribute{3}
    105 % DataOut.VarAttribute{4}
    106 % DataOut.VarAttribute{5}
    107 % DataOut.VarAttribute{6}
    108 % DataOut.VarAttribute{7}
    109137%------------------------------------------------
    110 function DataOut=phys_1(Data,Calib)
    111 % for icell=1:length(Data)
    112 
     138%--- transform a single field
     139function DataOut=phys_1(Data,Calib,ZIndex)
     140%------------------------------------------------
     141%% set default output
    113142DataOut=Data;%default
    114 % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
    115 if isfield(Calib,'SliceCoord') && isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
    116     DataOut.PlaneCoord=Calib.SliceCoord(Data.ZIndex,:);% transfer the slice position
    117     if isfield(Calib,'SliceAngle') % transfer the slice rotation angles
    118         if isequal(size(Calib.SliceAngle,1),1)
    119              DataOut.PlaneAngle=Calib.SliceAngle;
    120         else
    121             DataOut.PlaneAngle=Calib.SliceAngle(Data.ZIndex,:);
    122         end
    123     end
    124 end
    125 % The transform ACTS ONLY IF .CoordType='px'and Calib defined
    126 if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
    127     if isfield(Calib,'CoordUnit')
    128         DataOut.CoordUnit=Calib.CoordUnit;
    129     else
    130         DataOut.CoordUnit='cm'; %default
    131     end
    132     DataOut.TimeUnit='s';
    133     %transform of X,Y coordinates for vector fields
    134     test_z=0;
    135     if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
    136         Z=Data.ZIndex;
    137         test_z=1;
    138     else
    139         Z=0;
    140     end
    141     if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
    142         [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z);
    143         if test_z
    144              DataOut.ListVarName=[DataOut.ListVarName(1:2) {'Z'} DataOut.ListVarName(3:end)];
    145              DataOut.VarDimName=[DataOut.VarDimName(1:2) DataOut.VarDimName(1) DataOut.VarDimName(3:end)];
    146              ZAttribute{1}.Role='coord_z';
    147              DataOut.VarAttribute=[DataOut.VarAttribute(1:2) ZAttribute DataOut.VarAttribute(3:end)];
    148         end
    149         if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt')
    150             if ~isempty(Data.dt)
    151             [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);
    152             [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);
    153             DataOut.U=(XOut_2-XOut_1)/Data.dt;
    154             DataOut.V=(YOut_2-YOut_1)/Data.dt;
    155             end
    156         end
    157     end
    158     %transform of an image or scalar: done in phys_ima
    159      
    160     %transform of spatial derivatives
    161     if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
    162           && isfield(Data,'dt')   
    163         if ~isempty(Data.dt)
    164             % estimate the Jacobian matrix DXpx/DXphys
    165             for ip=1:length(Data.X)
    166                 [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
    167                 [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
    168                 [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
    169                 [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
    170             %Jacobian matrix DXpphys/DXpx
    171                DjXi(1,1)=(Xp1-Xm1);
    172                DjXi(2,1)=(Yp1-Ym1);
    173                DjXi(1,2)=(Xp2-Xm2);
    174                DjXi(2,2)=(Yp2-Ym2);
    175                DjUi(:,:)=Data.DjUi(ip,:,:);
    176                DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
    177                DataOut.DjUi(ip,:,:)=DjUi';
    178             end
    179             DataOut.DjUi =  DataOut.DjUi/Data.dt;   %     min(Data.DjUi(:,1,1))=DUDX                         
    180         end
    181     end
    182 end
     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%% transform of spatial derivatives: TODO check the case with plane angles
     167if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
     168      && isfield(Data,'dt')   
     169    if ~isempty(Data.dt)
     170        % estimate the Jacobian matrix DXpx/DXphys
     171        for ip=1:length(Data.X)
     172            [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
     173            [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
     174            [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
     175            [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
     176        %Jacobian matrix DXpphys/DXpx
     177           DjXi(1,1)=(Xp1-Xm1);
     178           DjXi(2,1)=(Yp1-Ym1);
     179           DjXi(1,2)=(Xp2-Xm2);
     180           DjXi(2,2)=(Yp2-Ym2);
     181           DjUi(:,:)=Data.DjUi(ip,:,:);
     182           DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
     183           DataOut.DjUi(ip,:,:)=DjUi';
     184        end
     185        DataOut.DjUi =  DataOut.DjUi/Dt;   %     min(Data.DjUi(:,1,1))=DUDX                         
     186    end
     187end
     188
    183189
    184190%%%%%%%%%%%%%%%%%%%%
    … …  
    217223for icell=1:length(A)
    218224    Calib=CalibIn{icell};
    219     if isfield(Calib,'R') || isfield(Calib,'kc')|| test_multi ||~isequal(Calib,CalibIn{1})% the image needs to be interpolated to the new coordinates
     225    % rescaling of the image coordinates without change of the image array
     226    if strcmp(Calib.CalibrationType,'rescale') && isequal(Calib,CalibIn{1})
     227        A_out{icell}=A{icell};%no transform
     228        Rangx=[0.5 npx-0.5];%image coordiantes of corners
     229        Rangy=[npy-0.5 0.5];
     230        [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
     231        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
     232    else         
     233        % the image needs to be interpolated to the new coordinates
    220234        zphys=0; %default
    221235        if isfield(Calib,'SliceCoord') %.Z= index of plane
    … …  
    263277        if testuint16
    264278            A_out{icell}=uint16(A_out{icell});
    265         end
    266     else%     
    267         A_out{icell}=A{icell};%no transform
    268         Rangx=[0.5 npx-0.5];%image coordiantes of corners
    269         Rangy=[npy-0.5 0.5];
    270         [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
    271         [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
    272     end
    273 end
    274 
     279        end     
     280    end
     281end
     282
Note: See TracChangeset for help on using the changeset viewer.