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 00 0 1 01/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:255for j = 1:255for k = 1:255f(i,j,k) = 1;endendenddeltaX = 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 = p0for i = 1:255v1 = v1 + p1;v2 = v1;for j = 1:255v2 = v2 + p2;v3 = v2;for k = 1:255v3 = v3 + p3u = v3(1)/v3(3)v = v3(2)/v3(3)update the projection over the neighbourhood of p = [u v]' using bilinear interpolationendendend