I see that you need help in computing projection matrix given that you have corresponding 2D and 3D points.
To compute projection matrix, you can follow the given steps:
- First, the Direct Linear Transformation (DLT) method is used to compute the projection matrix (P).
- This involves setting up a system of linear equations derived from the correspondence between 3D world points and 2D image points, which is then solved using Singular Value Decomposition (SVD).
- Once the projection matrix is obtained, it is decomposed into the intrinsic matrix (K), rotation matrix (R), and translation vector (t).
- This decomposition is achieved using RQ decomposition, where the matrix (M), a 3x3 submatrix of (P), is decomposed into (K) and (R) by performing a QR decomposition on its transpose and then rearranging the results.
- The intrinsic matrix (K) is normalized to ensure that its bottom-right element is 1, and the translation vector (t) is extracted from the last column of (P).
Kindly refer to the code below for detailed implementation of the above steps:
function cameraCalibration(worldPoints, imagePoints)
% Ensure the points are in homogeneous coordinates
if size(worldPoints, 2) ~= 3 || size(imagePoints, 2) ~= 2
error('worldPoints should be Nx3 and imagePoints should be Nx2');
end
% Number of points
n = size(worldPoints, 1);
% Construct the matrix A
A = [];
for i = 1:n
X = worldPoints(i, 1);
Y = worldPoints(i, 2);
Z = worldPoints(i, 3);
u = imagePoints(i, 1);
v = imagePoints(i, 2);
A = [A;
X, Y, Z, 1, 0, 0, 0, 0, -u*X, -u*Y, -u*Z, -u;
0, 0, 0, 0, X, Y, Z, 1, -v*X, -v*Y, -v*Z, -v];
end
% Solve for p using SVD
[~, ~, V] = svd(A);
P = reshape(V(:, end), [4, 3])'; % Projection matrix
% Extract the 3x3 submatrix from P
M = P(:, 1:3);
% Perform RQ decomposition
[K, R] = rq(M);
% Ensure K has positive diagonal entries
T = diag(sign(diag(K)));
K = K * T;
R = T * R;
% Compute the translation vector
t = K \ P(:, 4);
% Normalize K
K = K / K(3, 3);
% Display the results
disp('Projection Matrix (P):');
disp(P);
disp('Intrinsic Matrix (K):');
disp(K);
disp('Rotation Matrix (R):');
disp(R);
disp('Translation Vector (t):');
disp(t);
end
function [R, Q] = rq(M)
% RQ decomposition using QR decomposition of the transposed matrix
[Q, R] = qr(flipud(M)');
R = flipud(R');
Q = flipud(Q');
Q = Q(:, end:-1:1);
R = R(end:-1:1, :);
end
% Example usage
worldPoints = [0, 0, 0; 1, 0, 0; 1, 1, 0; 0, 1, 0]; % Example 3D points
imagePoints = [100, 100; 200, 100; 200, 200; 100, 200]; % Corresponding 2D points
cameraCalibration(worldPoints, imagePoints);
In the process above the intrinsic matrix is used as an alias for Calibration matrix.
I hope this helps!
