# 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.mclc;clear all;% Generation of the rotation matrixthetaDeg = 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-axisrotation = [ ...[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 matrixD = input('Distance = ')perspectiveTrans = [                      0       1       0       0                      0       0       1       0                      1/D     0       0       0                     ] % projection-matrixpTilde = perspectiveTrans * rotation% Extracting the columns of the projection-matrixp1 = pTilde(:,1)p2 = pTilde(:,2)p3 = pTilde(:,3)p4 = pTilde(:,4)% Generation of a cube of ones of size 255x255x255for i = 1:255      for j = 1:255          for k = 1:255              f(i,j,k) = 1;          end      endenddeltaX = 1; deltaY = 1; deltaZ = 1;% Pre-multiply the columns of projection-matrixp1 = 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 = p0for 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      endend``

``>>[u,v] = projectionmatrix(thetaDeg,phiDeg,psiDeg,D);``
``function [u v] = projectionmatrix(thetaDeg,phiDeg,psiDeg,D)% Generation of the projection-matrix% projectionmatrix.m% Generation of the rotation matrixtheta = thetaDeg*(pi/180);phi = phiDeg*(pi/180);psi = psiDeg*(pi/180);% rotation matrix with respect to x-, y-, and z-axisrotation = [ 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 matrixperspectiveTrans = [0       1       0       0                  0       0       1       0                  1/D     0       0       0 ];% projection-matrixpTilde = perspectiveTrans * rotation;% Extracting the columns of the projection-matrixp1 = pTilde(:,1);p2 = pTilde(:,2);p3 = pTilde(:,3);p4 = pTilde(:,4);% Generation of a cube of ones of size 255x255x255f = 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-matrixp1 = p1 * deltaX;p2 = p2 * deltaY;p3 = p3 * deltaZ;% The pseudocode for the voxel-driven reprojection% Loop over the 3-D lattice pointsp0 = [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  endend``