Changeset 1112 for trunk/src/transform_field

Timestamp:

Jan 26, 2022, 7:37:21 PM (3 years ago)

Author:

sommeria

Message:

set_slice separated from geometrey_calib

Location:

trunk/src/transform_field

Files:

• phys.m (modified) (10 diffs)
• phys_polar.m (modified) (13 diffs)

trunk/src/transform_field/phys.m

@@ -55 +55 @@
 DataOut=DataIn;%default first output field
 if nargin>=2 % nargin =nbre of input variables
+     Calib{1}=[];
     if isfield(XmlData,'GeometryCalib')
         Calib{1}=XmlData.GeometryCalib;
-    else
-        Calib{1}=[];
+    end
+    Slice{1}=Calib{1};
+    if isfield(XmlData,'Slice')
+        Slice{1}=XmlData.Slice;
     end
     if nargin>=3  %two input fields
         DataOut_1=DataIn_1;%default second output field
-        if nargin>=4 && isfield(XmlData_1,'GeometryCalib')
-            Calib{2}=XmlData_1.GeometryCalib;
-        else
-            Calib{2}=Calib{1}; 
+        Calib{2}=Calib{1};
+        if nargin>=4 
+            if isfield(XmlData_1,'GeometryCalib')
+                Calib{2}=XmlData_1.GeometryCalib;
+            end
+            Slice{2}=Calib{2};
+            if isfield(XmlData_1,'Slice')
+                Slice{2}=XmlData_1.Slice;
+            end
         end
     end
@@ -71 +79 @@
 
 %% get the z index defining the section plane
+ZIndex=1;
 if isfield(DataIn,'ZIndex')&&~isempty(DataIn.ZIndex)&&~isnan(DataIn.ZIndex)
     ZIndex=DataIn.ZIndex;
-else
-    ZIndex=1;
 end
 
@@ -82 +89 @@
 if  ~isempty(Calib{1})
     if isfield(Calib{1},'CalibrationType')&& isfield(Calib{1},'CoordUnit') && isfield(DataIn,'CoordUnit')&& strcmp(DataIn.CoordUnit,'pixel')   
-        DataOut=phys_1(DataIn,Calib{1},ZIndex);% transform coordinates and velocity components
+        DataOut=phys_1(DataIn,Calib{1},Slice{1},ZIndex);% transform coordinates and velocity components
         %case of images or scalar: in case of two input fields, we need to project the transform  on the same regular grid
         if isfield(DataIn,'A') && isfield(DataIn,'Coord_x') && ~isempty(DataIn.Coord_x) && isfield(DataIn,'Coord_y')&&...
@@ -94 +101 @@
 
 %% document the selected  plane position and angle if relevant
-if  checktransform && isfield(Calib{1},'SliceCoord')&&size(Calib{1}.SliceCoord,1)>=ZIndex
-    DataOut.PlaneCoord=Calib{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
-    if isfield(Calib{1},'SliceAngle') % transfer the slice rotation angles
-        if isequal(size(Calib{1}.SliceAngle,1),1)% case of a unique angle
-            DataOut.PlaneAngle=Calib{1}.SliceAngle;
+if  checktransform && isfield(Slice{1},'SliceCoord')&&size(Slice{1}.SliceCoord,1)>=ZIndex
+    DataOut.PlaneCoord=Slice{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex
+    if isfield(Slice{1},'SliceAngle') % transfer the slice rotation angles
+        if isequal(size(Slice{1}.SliceAngle,1),1)% case of a unique angle
+            DataOut.PlaneAngle=Slice{1}.SliceAngle;
         else  % case of multiple planes with different angles: select the plane with index ZIndex
-            DataOut.PlaneAngle=Calib{1}.SliceAngle(ZIndex,:);
+            DataOut.PlaneAngle=Slice{1}.SliceAngle(ZIndex,:);
         end
     end
@@ -113 +120 @@
     end
     if isfield(Calib{2},'CalibrationType')&&isfield(Calib{2},'CoordUnit') && isfield(DataIn_1,'CoordUnit')&& strcmp(DataIn_1.CoordUnit,'pixel')
-        DataOut_1=phys_1(DataOut_1,Calib{2},ZIndex);
-        if isfield(Calib{1},'SliceCoord')
-            if ~(isfield(Calib{2},'SliceCoord') && isequal(Calib{2}.SliceCoord,Calib{1}.SliceCoord))
+        DataOut_1=phys_1(DataOut_1,Calib{2},Slice{2},ZIndex);
+        if isfield(Slice{2},'SliceCoord')
+            if ~(isfield(Slice{2},'SliceCoord') && isequal(Slice{2}.SliceCoord,Slice{1}.SliceCoord))
                 DataOut_1.Txt='different plane positions for the two input fields';
                 return
             end
             DataOut_1.PlaneCoord=DataOut.PlaneCoord;% same plane position for the two input fields
-            if isfield(Calib{1},'SliceAngle')
-                if ~(isfield(Calib{2},'SliceAngle') && isequal(Calib{2}.SliceAngle,Calib{1}.SliceAngle))
+            if isfield(Slice{1},'SliceAngle')
+                if ~(isfield(Slice{2},'SliceAngle') && isequal(Slice{2}.SliceAngle,Slice{1}.SliceAngle))
                     DataOut_1.Txt='different plane angles for the two input fields';
                     return
@@ -131 +138 @@
                 ~isempty(DataIn_1.Coord_y)&&length(DataIn_1.A)>1
             iscalar=iscalar+1;
-            Calib{iscalar}=Calib{2};
+%             Calib{iscalar}=Calib{2};
             A{iscalar}=DataIn_1.A;
         end
@@ -163 +170 @@
 %------------------------------------------------
 %--- transform a single field
-function DataOut=phys_1(Data,Calib,ZIndex)
+function DataOut=phys_1(Data,Calib,Slice,ZIndex)
 %------------------------------------------------
 %% set default output
@@ -171 +178 @@
 %% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
 if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
-  [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
+  [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);
     Dt=1; %default
     if isfield(Data,'dt')&&~isempty(Data.dt)
@@ -180 +187 @@
     end
     if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
-        [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
-        [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
+        [XOut_1,YOut_1]=phys_XYZ(Calib,Slice,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
+        [XOut_2,YOut_2]=phys_XYZ(Calib,Slice,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
         DataOut.U=(XOut_2-XOut_1)/Dt;
         DataOut.V=(YOut_2-YOut_1)/Dt;
@@ -213 +220 @@
     % estimate the Jacobian matrix DXpx/DXphys
     for ip=1:length(Data.X)
-        [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
-        [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
-        [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
-        [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
+        [Xp1,Yp1]=phys_XYZ(Calib,Slice,Data.X(ip)+0.5,Data.Y(ip),ZIndex);
+        [Xm1,Ym1]=phys_XYZ(Calib,Slice,Data.X(ip)-0.5,Data.Y(ip),ZIndex);
+        [Xp2,Yp2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)+0.5,ZIndex);
+        [Xm2,Ym2]=phys_XYZ(Calib,Slice,Data.X(ip),Data.Y(ip)-0.5,ZIndex);
         %Jacobian matrix DXpphys/DXpx
         DjXi(1,1)=(Xp1-Xm1);

trunk/src/transform_field/phys_polar.m

@@ -78 +78 @@
 Data=DataIn; %default output
 if isfield(Data,'CoordUnit')
-Data=rmfield(Data,'CoordUnit');
+    Data=rmfield(Data,'CoordUnit');
 end
 Data.ListVarName = {};
@@ -85 +85 @@
 DataCell{1}=DataIn;
 Calib{1}=[];
+Slice{1}=[];
 DataCell{2}=[];%default
 checkpixel(1)=0;
-if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel') 
+if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel')
     checkpixel(1)=1;
 end
@@ -94 +95 @@
         Calib{1}=XmlData.GeometryCalib;
     end
+    Slice{1}=Calib{1};
+    if isfield(XmlData,'Slice')
+        Slice{1}=XmlData.Slice;
+    end
     Calib{2}=Calib{1};
+    Slice{2}=Slice{1};
 else
     Data.Txt='wrong input: need two or four structures';
@@ -101 +107 @@
 if nargin==4% case of two input fields
     checkpixel(2)=0;
-if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel') 
-    checkpixel(2)=1;
-end
+    if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel')
+        checkpixel(2)=1;
+    end
     DataCell{2}=DataIn_1;%default
     if isfield(XmlData_1,'GeometryCalib')&& ~isempty(XmlData_1.GeometryCalib) && checkpixel(2)
         Calib{2}=XmlData_1.GeometryCalib;
+    end
+    if isfield(XmlData_1,'Slice')
+        Slice{2}=XmlData_1.Slice;
     end
     nbinput=2;
@@ -218 +227 @@
                     if ~isempty(Calib{ifield})
                         [X,Y,Z,DataCell{ifield}.(UName),DataCell{ifield}.(VName)]=...
-                            phys_XYUV(DataCell{ifield},Calib{ifield},ZIndex);
+                            phys_XYUV(DataCell{ifield},Calib{ifield},Slice{ifield},ZIndex);
                     end
                 end
@@ -304 +313 @@
                     ZInd(nbgrid)=ZIndex;
                     Calib_new{nbgrid}=Calib{ifield};
+                    Slice_new{nbgrid}=Slice{ifield};
                 end
                 if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
@@ -325 +335 @@
                     Calib_new{nbgrid+1}=Calib{ifield};
                     Calib_new{nbgrid+2}=Calib{ifield};
+                    Slice_new{nbgrid+1}=Calib{ifield};
+                    Slice_new{nbgrid+2}=Calib{ifield};
                     nbgrid=nbgrid+2;
                     nbvar=nbvar+2;
@@ -335 +347 @@
 %% tranform cartesian to polar coordinates for gridded data
 if nbgrid~=0
-    [A,Data.radius,Data.theta]=phys_Ima_polar(A,coord_x,coord_y,Calib_new,ZInd,origin_xy,radius_offset,angle_offset,angle_scale);
+    [A,Data.radius,Data.theta]=phys_Ima_polar(A,coord_x,coord_y,Calib_new,Slice_new,ZInd,origin_xy,radius_offset,angle_offset,angle_scale);
     for icell=1:numel(A)
         if icell<=numel(A)-1 && check_vector(icell)==1 && check_vector(icell+1)==1   %transform u,v into polar coordinates
@@ -373 +385 @@
 %------------------------------------------------
 %--- transform a single field into phys coordiantes
-function [X,Y,Z,U,V]=phys_XYUV(Data,Calib,ZIndex)
+function [X,Y,Z,U,V]=phys_XYUV(Data,Calib,Slice,ZIndex)
 %------------------------------------------------
 %% set default output
@@ -385 +397 @@
 %% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
 if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
-    [X,Y,Z]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
+    [X,Y,Z]=phys_XYZ(Calib,Slice,Data.X,Data.Y,ZIndex);
     Dt=1; %default
     if isfield(Data,'dt')&&~isempty(Data.dt)
@@ -394 +406 @@
     end
     if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
-        [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
-        [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
+        [XOut_1,YOut_1]=phys_XYZ(Calib,Slice,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
+        [XOut_2,YOut_2]=phys_XYZ(Calib,Slice,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
         U=(XOut_2-XOut_1)/Dt;
         V=(YOut_2-YOut_1)/Dt;
@@ -404 +416 @@
 % tranform gridded field into polar coordiantes on a regular polar grid,
 % transform to phys coordiantes if requested by calibration input
-function [A_out,radius,theta]=phys_Ima_polar(A,coord_x,coord_y,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
+function [A_out,radius,theta]=phys_Ima_polar(A,coord_x,coord_y,CalibIn,SliceIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
 rcorner=[];
 thetacorner=[];
@@ -420 +432 @@
     % transform edges into phys coordinates if requested
     if ~isempty(CalibIn{icell})
-        [x_edge,y_edge]=phys_XYZ(CalibIn{icell},x_edge,y_edge,ZIndex(icell));% physical coordinates of the image edge
+        [x_edge,y_edge]=phys_XYZ(CalibIn{icell},SliceIn{icell},x_edge,y_edge,ZIndex(icell));% physical coordinates of the image edge
     end