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 endenddeltaX = 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 = 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

## Best Answer