I have implemented a perspective projection algorithm according to chapter 6 Computer Graphics Principles and Practices (CGP&P) by Foley, van Dam, Feiner, Hughes (2nd edition). I have
N'per = M * Sper * Spar * T (-prp) * R * T (-vrp).
As I understand it, the final image should be in canonical form size of (-1,-1) to (1,1) and z in (0,-1). However, the final image X-Y dimensions (see Figure 1) do not seem correct. I'm mostly trying to understand how the final image size is determined. I have included the matlab code below. My frustum (f) is defined by eyepoint (EP) at a specified lat/lon that has been converted to ECEF; distances: near plane (nDist) = 300; view plane (vDist) = 900; and far plane (fDist) = 25000. A line of sight (LOS) vector created at the EP is the center of projection. The frustum correctly finds and returns the buildings that within it along the LOS. Field of View is (10 deg x 10 deg). Now I'm just trying to project those buildings onto a defined window so that I can "quantize" (paint?) the grid and indicate which building is located at which x,y pair in the view plane. Unfortunately, because the window is not returning at the indicated size, it makes the painting more difficult for me. And besides, I'd just like to know what I'm doing wrong to not end up with the correct dimensions.
Matlab code (no attempts at optimizations or anything. just brute-force implementation!
function iPersProj = getPersProj(bldg, bi, f, plotpersp, fPersPlot)
color = [rand rand rand];
face = eFaces.bottom;
iPersProjBtm = persproj(f, bldg, face);
face = eFaces.top;
iPersProjTop = persproj(f, bldg, face);
iPersProj = [iPersProjTop;iPersProjBtm];
hold on;
scatter3(iPersProjTop(:,1), ...
iPersProjTop(:,2), ...
iPersProjTop(:,3),'+','CData',color);
scatter3(iPersProjBtm(:,1), ...
iPersProjBtm(:,2), ...
iPersProjBtm(:,3),'o','CData',color);
pPersProj=[iPersProjTop;
iPersProjTop(1,:); ...
iPersProjBtm; ...
iPersProjBtm(1,:); ...
iPersProjBtm(2,:); ...
iPersProjTop(4,:); ...
iPersProjTop(3,:); ...
iPersProjBtm(3,:); ...
iPersProjBtm(4,:); ...
iPersProjTop(2,:); ...
iPersProjTop(1,:)];
line (pPersProj(:,1), pPersProj(:,2),'Color',color);
text (pPersProj(1,1), pPersProj(1,2), int2str(bi));
end
function proj = persproj(f, bldg, face)
vrp = f.vC; %center view plane
vpn = f.Z; % LOS for frustum
cop = -f.EP;
F = f.vDist - f.nDist;
B = f.vDist - f.fDist;
umin = -5;
vmin = -5;
umax = 5;
vmax = 5;
R = getrotation (f);
Tvrp = gettranslation(-vrp);
ed = R * Tvrp * [f.EP 1]'; %translate eyepoint to camera?
prp = [0 0 ed(3)];
sh = getsh(prp, umax, umin, vmax, vmin);
Tprp = gettranslation(-prp);
vrpp = -prp(3); %(sh * Tprp * [0;0;0;1]); %vrp-prime per CGP&P
zmin = -(vrpp + F)/(vrpp+B);
zmax = -(vrpp + B)/(vrpp+B);
zprj = -vrpp/(vrpp+B);
sper = getsper(vrpp, B, umax, umin, vmax, vmin);
M=[ 1 0 0 0; ...
0 1 0 0; ...
0 0 1/(1+zmin) -zmin/(1+zmin); ...
0 0 -1 0];
proj = zeros(4,4);
for i=1:4
Q=bldg.coords(i,:,face);
uvdw = M * sper * sh * Tprp * R * Tvrp * [Q';1];
proj (i,1) = uvdw(1);
proj (i,2) = uvdw(2);
proj (i,3) = uvdw(3);
end
end
function sper = getsper (vrpz, B, umax, umin, vmax, vmin)
dx=umax-umin;
dy=vmax-vmin;
sper=zeros(4,4);
sper(1,1) = 2*vrpz/(dx*(vrpz+B));
sper(2,2) = 2*vrpz/(dy*(vrpz+B));
sper(3,3) = -1/(vrpz+B);
sper(4,4) = 1;
end
function sh = getsh (prp, umax, umin, vmax, vmin)
sx=umax+umin;
sy=vmax+vmin;
cw = [sx/2 sy/2 0 1]';
dop = cw - [prp 1]';
shx = - dop(1)/dop(3);
shy = - dop(2)/dop(3);
sh=zeros(4,4);
sh(1,1) = 1;
sh(2,2) = 1;
sh(3,3) = 1;
sh(4,4) = 1;
sh(1,3) = shx;
sh(2,3) = shy;
end
function R = getrotation (f)
rz = f.Z;
rx=cross(f.Y, rz);
rx=rx/norm(rx);
ry=cross(rz,rx);
R=zeros(4,4);
R(1,1:3) = rx;
R(2,1:3) = ry;
R(3,1:3) = rz;
R(4,4) = 1;
end
function T = gettranslation(p)
T = zeros(4,4);
T(1:3,4) = p';
T(1,1) = 1;
T(2,2) = 1;
T(3,3) = 1;
T(4,4) = 1;
end
Figure 1: Prospective Projection but dimensions are not (-1,-1) to (1,1)1
