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
Related Question