source: trunk/src/phys_XYZ.m @ 659

Last change on this file since 659 was 242, checked in by sommeria, 14 years ago

phys_XYZ introduced as a main function, geometry_calib corrected to deal with reverse x axis in mode linear

File size: 4.0 KB
Line 
1%------------------------------------------------------------------------
2%'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
3% function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)
4%
5%OUTPUT:
6%
7%INPUT:
8%Z: index of plane
9function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex)
10%------------------------------------------------------------------------
11testangle=0;
12test_refraction=0;
13if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Zindex
14    if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,[0 0 0])
15        testangle=1;
16        om=norm(Calib.SliceAngle(Zindex,:));%norm of rotation angle in radians
17        OmAxis=Calib.SliceAngle(Zindex,:)/om; %unit vector marking the rotation axis
18        cos_om=cos(pi*om/180);
19        sin_om=sin(pi*om/180);
20        coeff=OmAxis(3)*(1-cos_om);
21        norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
22        norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
23        norm_plane(3)=OmAxis(3)*coeff+cos_om;
24        Z0=norm_plane*Calib.SliceCoord(Zindex,:)'/norm_plane(3);
25    else
26        Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte
27    end
28    Z0virt=Z0;
29    if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
30        H=Calib.InterfaceCoord(3);
31        if H>Z0
32            Z0virt=H-(H-Z0)/Calib.RefractionIndex; %corrected z (virtual object)
33            test_refraction=1;
34        end
35    end
36else
37    Z0=0;
38    Z0virt=0;
39end
40if ~exist('X','var')||~exist('Y','var')
41    Xphys=[];
42    Yphys=[];%default
43    return
44end
45%coordinate transform
46if ~isfield(Calib,'fx_fy')
47    Calib.fx_fy=[1 1];
48end
49if ~isfield(Calib,'Tx_Ty_Tz')
50    Calib.Tx_Ty_Tz=[0 0 1];
51end
52if ~isfield(Calib,'Cx_Cy')
53    Calib.Cx_Cy=[0 0];
54end
55if ~isfield(Calib,'kc')
56    Calib.kc=0;
57end
58if isfield(Calib,'R')
59    R=(Calib.R)';
60    if testangle
61        a=-norm_plane(1)/norm_plane(3);
62        b=-norm_plane(2)/norm_plane(3);
63        if test_refraction
64            a=a/Calib.RefractionIndex;
65            b=b/Calib.RefractionIndex;
66        end
67        R(1)=R(1)+a*R(3);
68        R(2)=R(2)+b*R(3);
69        R(4)=R(4)+a*R(6);
70        R(5)=R(5)+b*R(6);
71        R(7)=R(7)+a*R(9);
72        R(8)=R(8)+b*R(9);
73    end
74    Tx=Calib.Tx_Ty_Tz(1);
75    Ty=Calib.Tx_Ty_Tz(2);
76    Tz=Calib.Tx_Ty_Tz(3);
77    f=Calib.fx_fy(1);%dpy=1; sx=1
78    %dpx=Calib.fx_fy(2)/Calib.fx_fy(1);
79    Dx=R(5)*R(7)-R(4)*R(8);
80    Dy=R(1)*R(8)-R(2)*R(7);
81    D0=(R(2)*R(4)-R(1)*R(5));
82    Z11=R(6)*R(8)-R(5)*R(9);
83    Z12=R(2)*R(9)-R(3)*R(8);
84    Z21=R(4)*R(9)-R(6)*R(7);
85    Z22=R(3)*R(7)-R(1)*R(9);
86    Zx0=R(3)*R(5)-R(2)*R(6);
87    Zy0=R(1)*R(6)-R(3)*R(4);
88    A11=R(8)*Ty-R(5)*Tz+Z11*Z0virt;
89    A12=R(2)*Tz-R(8)*Tx+Z12*Z0virt;
90    A21=-R(7)*Ty+R(4)*Tz+Z21*Z0virt;
91    A22=-R(1)*Tz+R(7)*Tx+Z22*Z0virt;
92    %     X0=Calib.fx_fy(1)*(R(5)*Tx-R(2)*Ty+Zx0*Z0virt);
93    %     Y0=Calib.fx_fy(2)*(-R(4)*Tx+R(1)*Ty+Zy0*Z0virt);
94    X0=(R(5)*Tx-R(2)*Ty+Zx0*Z0virt);
95    Y0=(-R(4)*Tx+R(1)*Ty+Zy0*Z0virt);
96    %px to camera:
97    %     Xd=dpx*(X-Calib.Cx_Cy(1)); % sensor coordinates
98    %     Yd=(Y-Calib.Cx_Cy(2));
99    Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
100    Yd=(Y-Calib.Cx_Cy(2))/Calib.fx_fy(2);
101    dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd);%/(f*f); %distortion factor
102    Xu=Xd./dist_fact;%undistorted sensor coordinates
103    Yu=Yd./dist_fact;
104    denom=Dx*Xu+Dy*Yu+D0;
105    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
106    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
107    if testangle
108        Zphys=Z0+a*Xphys+b*Yphys;
109    else
110        Zphys=Z0;
111    end
112else
113    Xphys=-Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1);
114    Yphys=-Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2);
115end
116
117%'px_XYZ': transform phys coordinates to image coordinates (px)
118%
119% OUPUT:
120% X,Y: array of coordinates in the image cooresponding to the input physical positions
121%                    (origin at lower leftcorner, unit=pixel)
122
123% INPUT:
124% Calib: structure containing the calibration parameters (read from the ImaDoc .xml file)
125% Xphys, Yphys: array of x,y physical coordinates
126% [Z0]: corresponding array of z physical coordinates (0 by default)
Note: See TracBrowser for help on using the repository browser.