Changeset 1108 for trunk

Timestamp:

Jan 7, 2022, 10:55:01 AM (3 years ago)

Author:

sommeria

Message:

phys_XYZ modified for angular scan

Location:

trunk/src

Files:

• angle2normal.m (modified) (1 diff)
• phys_XYZ.m (modified) (5 diffs)
• series.m (modified) (5 diffs)
• uvmat.m (modified) (2 diffs)

trunk/src/angle2normal.m

@@ -1 +1 @@
 %calculate the components of the unit vector norm_plane normal to the plane
 %defined by the rotation vector PlaneAngle (in degree) 
-% this gives the equation of the plane as norm_plane(1)x + norm_plane(2)y + norm_plane(2)z = 0
+% this gives the equation of the plane as norm_plane(1)x + norm_plane(2)y +norm_plane(3)z = cte
 
 function norm_plane=rotate(PlaneAngle)

trunk/src/phys_XYZ.m

@@ -4 +4 @@
 %
 %OUTPUT:
-%
+% Xphys,Yphys,Zphys: vector of phys coordinates corresponding to the input vector of image coordinates
 %INPUT:
-%Z: index of plane
+% Calib: Matlab structure containing the calibration parameters (pinhole camera model, see 
+% http://servforge.legi.grenoble-inp.fr/projects/soft-uvmat/wiki/UvmatHelp#GeometryCalib) and the
+%    parameters describing the illumination plane(s)
+%    .Tx_Ty_Tz: translation (3 phys coordinates) defining the origine of the camera frame
+%    .R : rotation matrix from phys to camera frame
+%    .fx_fy: focal length along each direction of the image
+% X, Y: vectors of X and Y image coordinates
+% ZIndex: index defining the current illumination plane in a volume scan
 
 %=======================================================================
@@ -28 +35 @@
 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex)
 %------------------------------------------------------------------------
-testangle=0;
-test_refraction=0;
+testangle=0;% =1 if the illumination plane is tilted with respect to the horizontal plane Xphys Yphys
+test_refraction=0;% =1 if the considered points are viewed through an horizontal interface (located at z=Calib.InterfaceCoord(3)') 
 Zphys=0; %default output
 if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&size(Calib.SliceCoord,1)>=Zindex
     if isfield(Calib, 'SliceAngle') && size(Calib.SliceAngle,1)>=Zindex && ~isequal(Calib.SliceAngle(Zindex,:),[0 0 0])
         testangle=1;
-        norm_plane=angle2normal(Calib.SliceAngle(Zindex,:));
+        norm_plane=angle2normal(Calib.SliceAngle(Zindex,:));% coordinates of the unit vector normal to the current illumination plane
     end
     Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte
     Z0virt=Z0;
     if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
-        H=Calib.InterfaceCoord(3);
+        H=Calib.InterfaceCoord(3);% z position of the water surface
         if H>Z0
             Z0virt=H-(H-Z0)/Calib.RefractionIndex; %corrected z (virtual object)
@@ -69 +76 @@
 if isfield(Calib,'R')
     R=(Calib.R)';
+    c=Z0virt;
     if testangle
+        % equation of the illumination plane: z=ax+by+c
         a=-norm_plane(1)/norm_plane(3);
         b=-norm_plane(2)/norm_plane(3);
         if test_refraction
-            a=a/Calib.RefractionIndex;
-            b=b/Calib.RefractionIndex;
+            avirt=a/Calib.RefractionIndex;
+            bvirt=b/Calib.RefractionIndex;
         end
-        R(1)=R(1)+a*R(3);
-        R(2)=R(2)+b*R(3);
-        R(4)=R(4)+a*R(6);
-        R(5)=R(5)+b*R(6);
-        R(7)=R(7)+a*R(9);
-        R(8)=R(8)+b*R(9);
+        c=Z0virt-avirt*Calib.SliceCoord(Zindex,1)-bvirt*Calib.SliceCoord(Zindex,2);% Z0 = (virtual) z coordinate on the rotation axis (assumed horizontal)
+                               % c=z coordinate at (x,y)=(0,0)
+        R(1)=R(1)+avirt*R(3);
+        R(2)=R(2)+bvirt*R(3);
+        R(4)=R(4)+avirt*R(6);
+        R(5)=R(5)+bvirt*R(6);
+        R(7)=R(7)+avirt*R(9);
+        R(8)=R(8)+bvirt*R(9);          
     end
     Tx=Calib.Tx_Ty_Tz(1);
@@ -93 +104 @@
     Z21=R(4)*R(9)-R(6)*R(7);
     Z22=R(3)*R(7)-R(1)*R(9);
-    Zx0=R(3)*R(5)-R(2)*R(6);
-    Zy0=R(1)*R(6)-R(3)*R(4);
-    A11=R(8)*Ty-R(5)*Tz+Z11*Z0virt;
-    A12=R(2)*Tz-R(8)*Tx+Z12*Z0virt;
-    A21=-R(7)*Ty+R(4)*Tz+Z21*Z0virt;
-    A22=-R(1)*Tz+R(7)*Tx+Z22*Z0virt;
-    X0=(R(5)*Tx-R(2)*Ty+Zx0*Z0virt);
-    Y0=(-R(4)*Tx+R(1)*Ty+Zy0*Z0virt);
+     Zx0=R(3)*R(5)-R(2)*R(6);
+     Zy0=R(1)*R(6)-R(3)*R(4);
+    B11=R(8)*Ty-R(5)*Tz+Z11*c;
+    B12=R(2)*Tz-R(8)*Tx+Z12*c;
+    B21=-R(7)*Ty+R(4)*Tz+Z21*c;
+    B22=-R(1)*Tz+R(7)*Tx+Z22*c;
+    X0=(R(5)*Tx-R(2)*Ty+Zx0*c);
+    Y0=(-R(4)*Tx+R(1)*Ty+Zy0*c);
     %px to camera:
     Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
@@ -116 +127 @@
     end
     denom=Dx*Xu+Dy*Yu+D0;
-    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
-    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
+    Xphys=(B11.*Xu+B12.*Yu+X0)./denom;%world coordinates
+    Yphys=(B21.*Xu+B22.*Yu+Y0)./denom;
     if testangle
         Zphys=Z0+a*Xphys+b*Yphys;

trunk/src/series.m

@@ -1057 +1057 @@
 if ~isempty(VideoObject)% case of movies
     imainfo=get(VideoObject);
+    if isfield(imainfo,'NumFrames')
+        imainfo.NumberOfFrames=imainfo.NumFrames;
+    end
     if isempty(j1_series) % frame index along i
         Time=zeros(imainfo.NumberOfFrames+1,2);
@@ -1707 +1710 @@
     OutputDir='';
     answer='';
-    if isfield(Param,'OutputSubDir')% possibly update the output dir if it already exists
+    if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')% possibly update the output dir if it already exists
         PathOut=get(handles.OutputPath,'String');
         if ~exist(PathOut,'dir') % test if  the dir  already exist
@@ -1734 +1737 @@
             end
         end
+
         SubDirOut=[Param.OutputSubDir Param.OutputDirExt];
         SubDirOutNew=SubDirOut;
@@ -1782 +1786 @@
         OutputDir=Param.ActionInput.LogPath;
     end
-    if isfield(Param,'OutputSubDir') 
-    set(handles.OutputSubDir,'String',Param.OutputSubDir)
-    set(handles.OutputDirExt,'String',Param.OutputDirExt)
-    drawnow
+    if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')
+        set(handles.OutputSubDir,'String',Param.OutputSubDir)
+        set(handles.OutputDirExt,'String',Param.OutputDirExt)
+        drawnow
     end
     if get(handles.Replicate,'Value')
@@ -2551 +2555 @@
             end
             Coord_z=get(handles.Coord_z,'String');
-            if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))
+            if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))REFRESH
                 FieldList={};
                 set(handles.Coord_z,'String','')

trunk/src/uvmat.m

@@ -1374 +1374 @@
 Angle_1=linspace(SliceData.SliceAngle_1(1),SliceData.SliceAngle_1(2),SliceData.NbSlice);
 Angle_2=linspace(SliceData.SliceAngle_2(1),SliceData.SliceAngle_2(2),SliceData.NbSlice);
-GeometryCalib.SliceAngle(:,1)=Angle_1';%rotation around x axis (to generalise)
-GeometryCalib.SliceAngle(:,2)=Angle_2';%rotation around y axis (to generalise)
+GeometryCalib.SliceAngle(:,1)=Angle_1';%rotation angle around x axis 
+GeometryCalib.SliceAngle(:,2)=Angle_2';%rotation angle around y axis 
 GeometryCalib.SliceAngle(:,3)=0;
 if SliceData.CheckRefraction
@@ -1435 +1435 @@
             else
                 if check_update
-                    display([XmlName ' updated with calibration parameters'])
+                    display([XmlName ' updated with slice positions'])
                 else
-                    display([XmlName ' created with calibration parameters'])
+                    display([XmlName ' created with slice positions'])
                 end
             end