Projection of a cube

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

[Sean de Wolski]

Walter, you can edit others' posts? That's awesome!

Answer 1

[Sean de Wolski]

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

[Srinidhi Murthy]

Thank you very much, Sean.

Yes, the 'f' is never used. My logic itself is wrong here. I have been thinking that I am creating a voxelized cube and then projecting it. I was finding it very hard to imagine. I think now I have a clear idea.

[Srinidhi Murthy]

Hi Sean. I have a problem in understanding the bilinear interpolation. I am confused with updating the projection values over the neighborhood. Can you please suggest me something here? Thank you.

[Sean de Wolski]

Have you looked at MATLAB's stock interp3/interpn? What about John's interpne:

http://www.mathworks.com/matlabcentral/fileexchange/10522-interpne

Can you please be more specific as to what you mean by "updating projection values over a neighborhood."