Ignore:
Timestamp:
Jan 26, 2022, 7:37:21 PM (2 years ago)
Author:
sommeria
Message:

set_slice separated from geometrey_calib

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/transform_field/phys.m

    r1107 r1112  
    5555DataOut=DataIn;%default first output field
    5656if nargin>=2 % nargin =nbre of input variables
     57     Calib{1}=[];
    5758    if isfield(XmlData,'GeometryCalib')
    5859        Calib{1}=XmlData.GeometryCalib;
    59     else
    60         Calib{1}=[];
     60    end
     61    Slice{1}=Calib{1};
     62    if isfield(XmlData,'Slice')
     63        Slice{1}=XmlData.Slice;
    6164    end
    6265    if nargin>=3  %two input fields
    6366        DataOut_1=DataIn_1;%default second output field
    64         if nargin>=4 && isfield(XmlData_1,'GeometryCalib')
    65             Calib{2}=XmlData_1.GeometryCalib;
    66         else
    67             Calib{2}=Calib{1};
     67        Calib{2}=Calib{1};
     68        if nargin>=4
     69            if isfield(XmlData_1,'GeometryCalib')
     70                Calib{2}=XmlData_1.GeometryCalib;
     71            end
     72            Slice{2}=Calib{2};
     73            if isfield(XmlData_1,'Slice')
     74                Slice{2}=XmlData_1.Slice;
     75            end
    6876        end
    6977    end
     
    7179
    7280%% get the z index defining the section plane
     81ZIndex=1;
    7382if isfield(DataIn,'ZIndex')&&~isempty(DataIn.ZIndex)&&~isnan(DataIn.ZIndex)
    7483    ZIndex=DataIn.ZIndex;
    75 else
    76     ZIndex=1;
    7784end
    7885
     
    8289if  ~isempty(Calib{1})
    8390    if isfield(Calib{1},'CalibrationType')&& isfield(Calib{1},'CoordUnit') && isfield(DataIn,'CoordUnit')&& strcmp(DataIn.CoordUnit,'pixel')   
    84         DataOut=phys_1(DataIn,Calib{1},ZIndex);% transform coordinates and velocity components
     91        DataOut=phys_1(DataIn,Calib{1},Slice{1},ZIndex);% transform coordinates and velocity components
    8592        %case of images or scalar: in case of two input fields, we need to project the transform  on the same regular grid
    8693        if isfield(DataIn,'A') && isfield(DataIn,'Coord_x') && ~isempty(DataIn.Coord_x) && isfield(DataIn,'Coord_y')&&...
     
    94101
    95102%% document the selected  plane position and angle if relevant
    96 if  checktransform && isfield(Calib{1},'SliceCoord')&&size(Calib{1}.SliceCoord,1)>=ZIndex
    97     DataOut.PlaneCoord=Calib{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
    98     if isfield(Calib{1},'SliceAngle') % transfer the slice rotation angles
    99         if isequal(size(Calib{1}.SliceAngle,1),1)% case of a unique angle
    100             DataOut.PlaneAngle=Calib{1}.SliceAngle;
     103if  checktransform && isfield(Slice{1},'SliceCoord')&&size(Slice{1}.SliceCoord,1)>=ZIndex
     104    DataOut.PlaneCoord=Slice{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
     105    if isfield(Slice{1},'SliceAngle') % transfer the slice rotation angles
     106        if isequal(size(Slice{1}.SliceAngle,1),1)% case of a unique angle
     107            DataOut.PlaneAngle=Slice{1}.SliceAngle;
    101108        else  % case of multiple planes with different angles: select the plane with index ZIndex
    102             DataOut.PlaneAngle=Calib{1}.SliceAngle(ZIndex,:);
     109            DataOut.PlaneAngle=Slice{1}.SliceAngle(ZIndex,:);
    103110        end
    104111    end
     
    113120    end
    114121    if isfield(Calib{2},'CalibrationType')&&isfield(Calib{2},'CoordUnit') && isfield(DataIn_1,'CoordUnit')&& strcmp(DataIn_1.CoordUnit,'pixel')
    115         DataOut_1=phys_1(DataOut_1,Calib{2},ZIndex);
    116         if isfield(Calib{1},'SliceCoord')
    117             if ~(isfield(Calib{2},'SliceCoord') && isequal(Calib{2}.SliceCoord,Calib{1}.SliceCoord))
     122        DataOut_1=phys_1(DataOut_1,Calib{2},Slice{2},ZIndex);
     123        if isfield(Slice{2},'SliceCoord')
     124            if ~(isfield(Slice{2},'SliceCoord') && isequal(Slice{2}.SliceCoord,Slice{1}.SliceCoord))
    118125                DataOut_1.Txt='different plane positions for the two input fields';
    119126                return
    120127            end
    121128            DataOut_1.PlaneCoord=DataOut.PlaneCoord;% same plane position for the two input fields
    122             if isfield(Calib{1},'SliceAngle')
    123                 if ~(isfield(Calib{2},'SliceAngle') && isequal(Calib{2}.SliceAngle,Calib{1}.SliceAngle))
     129            if isfield(Slice{1},'SliceAngle')
     130                if ~(isfield(Slice{2},'SliceAngle') && isequal(Slice{2}.SliceAngle,Slice{1}.SliceAngle))
    124131                    DataOut_1.Txt='different plane angles for the two input fields';
    125132                    return
     
    131138                ~isempty(DataIn_1.Coord_y)&&length(DataIn_1.A)>1
    132139            iscalar=iscalar+1;
    133             Calib{iscalar}=Calib{2};
     140%             Calib{iscalar}=Calib{2};
    134141            A{iscalar}=DataIn_1.A;
    135142        end
     
    163170%------------------------------------------------
    164171%--- transform a single field
    165 function DataOut=phys_1(Data,Calib,ZIndex)
     172function DataOut=phys_1(Data,Calib,Slice,ZIndex)
    166173%------------------------------------------------
    167174%% set default output
     
    171178%% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
    172179if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
    173   [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
     180  [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);
    174181    Dt=1; %default
    175182    if isfield(Data,'dt')&&~isempty(Data.dt)
     
    180187    end
    181188    if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
    182         [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
    183         [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
     189        [XOut_1,YOut_1]=phys_XYZ(Calib,Slice,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
     190        [XOut_2,YOut_2]=phys_XYZ(Calib,Slice,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
    184191        DataOut.U=(XOut_2-XOut_1)/Dt;
    185192        DataOut.V=(YOut_2-YOut_1)/Dt;
     
    213220    % estimate the Jacobian matrix DXpx/DXphys
    214221    for ip=1:length(Data.X)
    215         [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
    216         [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
    217         [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
    218         [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
     222        [Xp1,Yp1]=phys_XYZ(Calib,Slice,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
     223        [Xm1,Ym1]=phys_XYZ(Calib,Slice,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
     224        [Xp2,Yp2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
     225        [Xm2,Ym2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
    219226        %Jacobian matrix DXpphys/DXpx
    220227        DjXi(1,1)=(Xp1-Xm1);
Note: See TracChangeset for help on using the changeset viewer.