How to calculate the plane at infinity? either in Matlab or Python

Question

• I tried to implement this function to compute the plane at infinity $n_\Pi$ from the equation: where $w^{-1}$ is the image of the absolute conic IAC, $M_j$ is the Projection Matrix of image $j$,

• but I'm not sure if this implementation is right, because I don't have a good idea about symbolic expressions and how they are solved.

• so I just mimicked other codes, but I'm still confused, because when I used this part to calculate the intrinsic matrix, the results of the intrinsic matrix weren't converging.
• can you please tell me if there are issues in my implementation, and how to fix it please?

The whole code, and data is here

Thanks in Advance

function P_infty = PlaneAtInfinity(K,P)

% Symbolic variables

X =  sym('X', 'real');
Y =  sym('Y', 'real');
Z =  sym('Z', 'real');
L2 = sym('L2','real');

N = [X; Y; Z];
W_inv = K*K';
% Quadric
Q = [W_inv       ,   (W_inv * N)    ;
    (N' * W_inv) , (N' * W_inv * N)];

% Autocalibration equation
Calibration = P * Q * P';
% func = cross(W_inv, Calibration);
% solution = Optimize(zeros(3,3),@func);
% solve linear equations
solution = solve(Calibration(1, 1) == (L2 * W_inv(1, 1)), ...
    Calibration(2, 2) == (L2 * W_inv(2, 2)), ...
    Calibration(3, 3) == (L2 * W_inv(3, 3)), ...
    Calibration(1, 3) == (L2 * W_inv(1, 3)));

P_infty = [double(solution.X(1));...
    double(solution.Y(1));...
    double(solution.Z(1))];
end

and the main script:

P = [4.35310807409356, 1.44651724424134, 36.0849770394682, 6.6274893804074;
22.9809618064503, 8.56379922408173, -2.94837061198918, 36.6987313462383;
-0.00438517470946907, -0.00160149495653558, 0.241477646472821, -0.00696180145320084];

K = [870.258271610840, 0, 279.100746891319;
0, 812.178628533806, 260.930050222247;
0, 0, 1];
PlaneAtInfinity(K,P)

Answer 1

The plane at infinity has to be calculated as following:

function plane = computePlaneAtInfinity(P, K)

%Input
%   P - Projection matrices
%   K - Approximate Values of Intrinsics
%
%Output
%   plane - coordinate of plane at infinity

% Compute the DIAC W^-1 
W_invert = K * K';

% Construct Symbolic Variables to Solve for Plane at Infinity
% X,Y,Z is the coordinate of plane at infinity
% XX is the scale
X = sym('X', 'real'); 
Y = sym('Y', 'real');
Z = sym('Z', 'real');
XX = sym('XX', 'real');

% Define Normal to Plane at Infinity
N = [X; Y; Z]; 

% Equation of Dual Absolute Quadric (DAQ)
Q = [W_invert, (W_invert * N); (N' * W_invert), (N' * W_invert * N)]; 

% Select Any One Projection Matrix
M = P(:, :, 2);

% Left hand side of the equation
LHS = M * Q * M';

% Solve for [X, Y, Z] considering the System of Linear Equations
% We need 4 equations for 4 variables X,Y,Z,XX
S = solve(LHS(1, 1) == (XX * W_invert(1, 1)), ...
          LHS(1, 2) == (XX * W_invert(1, 2)), ...
          LHS(1, 3) == (XX * W_invert(1, 3)), ...
          LHS(2, 2) == (XX * W_invert(2, 2)));

plane = [double(S.X(1)); double(S.Y(1)); double(S.Z(1))];

end
```