MATLAB: Projection of a cube

cubeprojection matrixvoxel-driven

I am learning and trying to implement the voxel-driven reprojection using projection matrices. I have derived the perspective transformation matrix and constructed the projection-matrix.
I converted the perspective transformation matrix into program and could generate the projection matrix. But I am stuck with the reprojection part. I have the pseudocode for the voxel-driven reprojection. I have tried a million times, but was unable to program it properly. Can somebody suggest or give me a tip on where I am going wrong.
Here is the code:
% EDITORIAL NOTE : I reformatted Srinidhi's rotation matrix as the width limitations made it too messy to follow easily — Walter
% Generation of the projection-matrix
% projectionmatrix.m
clc;
clear all;
% Generation of the rotation matrix
thetaDeg = input('theta = ')
phiDeg = input('phi = ')
psiDeg = input('psi = ')
theta = thetaDeg*(pi/180);
phi = phiDeg*(pi/180);
psi = psiDeg*(pi/180);
% rotation matrix with respect to x-, y-, and z-axis
rotation = [ ...
[cos(theta)*cos(phi) (cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) (cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) 0]; ...
[sin(phi)*cos(theta) (sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) (sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) 0]; ...
[-sin(theta) cos(theta)*sin(psi) cos(theta)*cos(psi) 0]; ...
[0 0 0 1]];
% Generation of the 3x4 perspective transformation matrix
D = input('Distance = ')
perspectiveTrans = [
0 1 0 0
0 0 1 0
1/D 0 0 0
]
% projection-matrix
pTilde = perspectiveTrans * rotation
% Extracting the columns of the projection-matrix
p1 = pTilde(:,1)
p2 = pTilde(:,2)
p3 = pTilde(:,3)
p4 = pTilde(:,4)
% Generation of a cube of ones of size 255x255x255
for i = 1:255
for j = 1:255
for k = 1:255
f(i,j,k) = 1;
end
end
end
deltaX = 1; deltaY = 1; deltaZ = 1;
% Pre-multiply the columns of projection-matrix
p1 = p1 * deltaX;
p2 = p2 * deltaY;
p3 = p3 * deltaZ;
The pseudocode for the voxel-driven reprojection
% Loop over the 3-D lattice points
p0 = [0 0 0]';
v1 = p0
for i = 1:255
v1 = v1 + p1;
v2 = v1;
for j = 1:255
v2 = v2 + p2;
v3 = v2;
for k = 1:255
v3 = v3 + p3
u = v3(1)/v3(3)
v = v3(2)/v3(3)
update the projection over the neighbourhood of p = [u v]' using bilinear interpolation
end
end
end

Best Answer

  • Okay; I've rewritten this with a few pointers and ideas. First off; NEVER use clear all. It does all sorts of bad things. Second, use a function (entering input values 4 times every time you run this must get kind of old). The other comments are contained in the below.
    Call this on the command line with:
    >>[u,v] = projectionmatrix(thetaDeg,phiDeg,psiDeg,D);
    Set breakpoints to debug it as it runs.
    function [u v] = projectionmatrix(thetaDeg,phiDeg,psiDeg,D)
    % Generation of the projection-matrix
    % projectionmatrix.m
    % Generation of the rotation matrix
    theta = thetaDeg*(pi/180);
    phi = phiDeg*(pi/180);
    psi = psiDeg*(pi/180);
    % rotation matrix with respect to x-, y-, and z-axis
    rotation = [ cos(theta)*cos(phi) cos(phi)*sin(theta)*sin(psi)-sin(phi)*cos(psi) cos(phi)*sin(theta)*cos(psi)+sin(phi)*sin(psi) 0
    sin(phi)*cos(theta) sin(phi)*sin(theta)*sin(psi)+cos(phi)*cos(psi) sin(phi)*sin(theta)*cos(psi)-cos(phi)*sin(psi) 0
    -sin(theta) cos(theta)*sin(psi) cos(theta)*cos(psi) 0
    0 0 0 1];
    % Generation of the 3x4 perspective transformation matrix
    perspectiveTrans = [0 1 0 0
    0 0 1 0
    1/D 0 0 0 ];
    % projection-matrix
    pTilde = perspectiveTrans * rotation;
    % Extracting the columns of the projection-matrix
    p1 = pTilde(:,1);
    p2 = pTilde(:,2);
    p3 = pTilde(:,3);
    p4 = pTilde(:,4);
    % Generation of a cube of ones of size 255x255x255
    f = ones(255,255,255); %one liner; faster better!
    %f is never used, why did you make it?
    deltaX = 1; deltaY = 1; deltaZ = 1;
    % Pre-multiply the columns of projection-matrix
    p1 = p1 * deltaX;
    p2 = p2 * deltaY;
    p3 = p3 * deltaZ;
    % The pseudocode for the voxel-driven reprojection
    % Loop over the 3-D lattice points
    p0 = [0 0 0]';
    v1 = p0;
    u = zeros(255,255,255);
    v = u;
    for i = 1:255
    v1 = v1 + p1;
    v2 = v1;
    for j = 1:255
    v2 = v2 + p2;
    v3 = v2;
    for k = 1:255
    v3 = v3 + p3;
    u(i,j,k) = v3(1)/v3(3);%Store; don't overwrite!
    v(i,j,k) = v3(2)/v3(3);
    %update the projection over the neighbourhood of p = [u v]' using bilinear interpolation
    end
    end
    end