Changeset 1108


Ignore:
Timestamp:
Jan 7, 2022, 10:55:01 AM (2 years ago)
Author:
sommeria
Message:

phys_XYZ modified for angular scan

Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/angle2normal.m

    r1080 r1108  
    11%calculate the components of the unit vector norm_plane normal to the plane
    22%defined by the rotation vector PlaneAngle (in degree)
    3 % this gives the equation of the plane as norm_plane(1)x + norm_plane(2)y + norm_plane(2)z = 0
     3% this gives the equation of the plane as norm_plane(1)x + norm_plane(2)y +norm_plane(3)z = cte
    44
    55function norm_plane=rotate(PlaneAngle)
  • trunk/src/phys_XYZ.m

    r1107 r1108  
    44%
    55%OUTPUT:
    6 %
     6% Xphys,Yphys,Zphys: vector of phys coordinates corresponding to the input vector of image coordinates
    77%INPUT:
    8 %Z: index of plane
     8% Calib: Matlab structure containing the calibration parameters (pinhole camera model, see
     9% http://servforge.legi.grenoble-inp.fr/projects/soft-uvmat/wiki/UvmatHelp#GeometryCalib) and the
     10%    parameters describing the illumination plane(s)
     11%    .Tx_Ty_Tz: translation (3 phys coordinates) defining the origine of the camera frame
     12%    .R : rotation matrix from phys to camera frame
     13%    .fx_fy: focal length along each direction of the image
     14% X, Y: vectors of X and Y image coordinates
     15% ZIndex: index defining the current illumination plane in a volume scan
    916
    1017%=======================================================================
     
    2835function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex)
    2936%------------------------------------------------------------------------
    30 testangle=0;
    31 test_refraction=0;
     37testangle=0;% =1 if the illumination plane is tilted with respect to the horizontal plane Xphys Yphys
     38test_refraction=0;% =1 if the considered points are viewed through an horizontal interface (located at z=Calib.InterfaceCoord(3)')
    3239Zphys=0; %default output
    3340if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&size(Calib.SliceCoord,1)>=Zindex
    3441    if isfield(Calib, 'SliceAngle') && size(Calib.SliceAngle,1)>=Zindex && ~isequal(Calib.SliceAngle(Zindex,:),[0 0 0])
    3542        testangle=1;
    36         norm_plane=angle2normal(Calib.SliceAngle(Zindex,:));
     43        norm_plane=angle2normal(Calib.SliceAngle(Zindex,:));% coordinates of the unit vector normal to the current illumination plane
    3744    end
    3845    Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte
    3946    Z0virt=Z0;
    4047    if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
    41         H=Calib.InterfaceCoord(3);
     48        H=Calib.InterfaceCoord(3);% z position of the water surface
    4249        if H>Z0
    4350            Z0virt=H-(H-Z0)/Calib.RefractionIndex; %corrected z (virtual object)
     
    6976if isfield(Calib,'R')
    7077    R=(Calib.R)';
     78    c=Z0virt;
    7179    if testangle
     80        % equation of the illumination plane: z=ax+by+c
    7281        a=-norm_plane(1)/norm_plane(3);
    7382        b=-norm_plane(2)/norm_plane(3);
    7483        if test_refraction
    75             a=a/Calib.RefractionIndex;
    76             b=b/Calib.RefractionIndex;
     84            avirt=a/Calib.RefractionIndex;
     85            bvirt=b/Calib.RefractionIndex;
    7786        end
    78         R(1)=R(1)+a*R(3);
    79         R(2)=R(2)+b*R(3);
    80         R(4)=R(4)+a*R(6);
    81         R(5)=R(5)+b*R(6);
    82         R(7)=R(7)+a*R(9);
    83         R(8)=R(8)+b*R(9);
     87        c=Z0virt-avirt*Calib.SliceCoord(Zindex,1)-bvirt*Calib.SliceCoord(Zindex,2);% Z0 = (virtual) z coordinate on the rotation axis (assumed horizontal)
     88                               % c=z coordinate at (x,y)=(0,0)
     89        R(1)=R(1)+avirt*R(3);
     90        R(2)=R(2)+bvirt*R(3);
     91        R(4)=R(4)+avirt*R(6);
     92        R(5)=R(5)+bvirt*R(6);
     93        R(7)=R(7)+avirt*R(9);
     94        R(8)=R(8)+bvirt*R(9);         
    8495    end
    8596    Tx=Calib.Tx_Ty_Tz(1);
     
    93104    Z21=R(4)*R(9)-R(6)*R(7);
    94105    Z22=R(3)*R(7)-R(1)*R(9);
    95     Zx0=R(3)*R(5)-R(2)*R(6);
    96     Zy0=R(1)*R(6)-R(3)*R(4);
    97     A11=R(8)*Ty-R(5)*Tz+Z11*Z0virt;
    98     A12=R(2)*Tz-R(8)*Tx+Z12*Z0virt;
    99     A21=-R(7)*Ty+R(4)*Tz+Z21*Z0virt;
    100     A22=-R(1)*Tz+R(7)*Tx+Z22*Z0virt;
    101     X0=(R(5)*Tx-R(2)*Ty+Zx0*Z0virt);
    102     Y0=(-R(4)*Tx+R(1)*Ty+Zy0*Z0virt);
     106     Zx0=R(3)*R(5)-R(2)*R(6);
     107     Zy0=R(1)*R(6)-R(3)*R(4);
     108    B11=R(8)*Ty-R(5)*Tz+Z11*c;
     109    B12=R(2)*Tz-R(8)*Tx+Z12*c;
     110    B21=-R(7)*Ty+R(4)*Tz+Z21*c;
     111    B22=-R(1)*Tz+R(7)*Tx+Z22*c;
     112    X0=(R(5)*Tx-R(2)*Ty+Zx0*c);
     113    Y0=(-R(4)*Tx+R(1)*Ty+Zy0*c);
    103114    %px to camera:
    104115    Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
     
    116127    end
    117128    denom=Dx*Xu+Dy*Yu+D0;
    118     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
    119     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
     129    Xphys=(B11.*Xu+B12.*Yu+X0)./denom;%world coordinates
     130    Yphys=(B21.*Xu+B22.*Yu+Y0)./denom;
    120131    if testangle
    121132        Zphys=Z0+a*Xphys+b*Yphys;
  • trunk/src/series.m

    r1107 r1108  
    10571057if ~isempty(VideoObject)% case of movies
    10581058    imainfo=get(VideoObject);
     1059    if isfield(imainfo,'NumFrames')
     1060        imainfo.NumberOfFrames=imainfo.NumFrames;
     1061    end
    10591062    if isempty(j1_series) % frame index along i
    10601063        Time=zeros(imainfo.NumberOfFrames+1,2);
     
    17071710    OutputDir='';
    17081711    answer='';
    1709     if isfield(Param,'OutputSubDir')% possibly update the output dir if it already exists
     1712    if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')% possibly update the output dir if it already exists
    17101713        PathOut=get(handles.OutputPath,'String');
    17111714        if ~exist(PathOut,'dir') % test if  the dir  already exist
     
    17341737            end
    17351738        end
     1739
    17361740        SubDirOut=[Param.OutputSubDir Param.OutputDirExt];
    17371741        SubDirOutNew=SubDirOut;
     
    17821786        OutputDir=Param.ActionInput.LogPath;
    17831787    end
    1784     if isfield(Param,'OutputSubDir') 
    1785     set(handles.OutputSubDir,'String',Param.OutputSubDir)
    1786     set(handles.OutputDirExt,'String',Param.OutputDirExt)
    1787     drawnow
     1788    if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')
     1789        set(handles.OutputSubDir,'String',Param.OutputSubDir)
     1790        set(handles.OutputDirExt,'String',Param.OutputDirExt)
     1791        drawnow
    17881792    end
    17891793    if get(handles.Replicate,'Value')
     
    25512555            end
    25522556            Coord_z=get(handles.Coord_z,'String');
    2553             if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))
     2557            if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))REFRESH
    25542558                FieldList={};
    25552559                set(handles.Coord_z,'String','')
  • trunk/src/uvmat.m

    r1107 r1108  
    13741374Angle_1=linspace(SliceData.SliceAngle_1(1),SliceData.SliceAngle_1(2),SliceData.NbSlice);
    13751375Angle_2=linspace(SliceData.SliceAngle_2(1),SliceData.SliceAngle_2(2),SliceData.NbSlice);
    1376 GeometryCalib.SliceAngle(:,1)=Angle_1';%rotation around x axis (to generalise)
    1377 GeometryCalib.SliceAngle(:,2)=Angle_2';%rotation around y axis (to generalise)
     1376GeometryCalib.SliceAngle(:,1)=Angle_1';%rotation angle around x axis
     1377GeometryCalib.SliceAngle(:,2)=Angle_2';%rotation angle around y axis
    13781378GeometryCalib.SliceAngle(:,3)=0;
    13791379if SliceData.CheckRefraction
     
    14351435            else
    14361436                if check_update
    1437                     display([XmlName ' updated with calibration parameters'])
     1437                    display([XmlName ' updated with slice positions'])
    14381438                else
    1439                     display([XmlName ' created with calibration parameters'])
     1439                    display([XmlName ' created with slice positions'])
    14401440                end
    14411441            end
Note: See TracChangeset for help on using the changeset viewer.