3D Projective Viewing emphasizing Z-value

Ross Beveridge, October 30, 2018

This worksheet is a simplified and modified variant on the perspective projection pipeline notebook. This note book makes simplifications, namely placing the camera at the origin looking out the negative Z-axis, in order to concentrate attention on the peuedo-depth value, i.e. z-value, computation. 

from sage.plot.plot3d.base import SHOW_DEFAULTS
SHOW_DEFAULTS['aspect_ratio'] = (1,1,1)
SHOW_DEFAULTS['frame_aspect_ratio'] = (1,1,1)
SHOW_DEFAULTS['perspective_depth'] = false
%display latex
latex.matrix_delimiters(left='|', right='|')
latex.vector_delimiters(left='[', right=']')

Construct the house model.  This means a list of lists to define the vertices along with a 4 by 10 matrix of vertices in homogeneous coordinates wtih one vertex per column.  Finally, the faces are named in a dictionary that includes both a color to draw the face and a sequence of face vertex indices.

VVH1 = [[-8,-8,-30],[-8,2,-30],[0,8,-30],[8,2,-30],[8,-8,-30],[-8,-8,-54],[-8,2,-54],[0,8,-54],[8,2,-54],[8,-8,-54]]
MMH1 = matrix( len(VVH1), 3, VVH1).augment(matrix(ZZ, len(VVH1), 1, lambda x,y: 1)) 
MMH1 = MMH1.transpose()
faceV   = {
'houseFront' : ((0,1,2,3,4),'tomato'),     'houseBack' : ((5,6,7,8,9), 'green'),
'wallRight'   : ((0,1,6,5),  'brown'),      'wallLeft' : ((4,3,8,9),   'orange'), 
'roofRight'   : ((1,2,7,6),  'greenyellow'),'roofLeft' : ((3,2,7,8),   'springgreen'),
'floor'      : ((0,4,9,5),  'darkcyan')}
MMH1

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left|\begin{array}{rrrrrrrrrr} -8 & -8 & 0 & 8 & 8 & -8 & -8 & 0 & 8 & 8 \\ -8 & 2 & 8 & 2 & -8 & -8 & 2 & 8 & 2 & -8 \\ -30 & -30 & -30 & -30 & -30 & -54 & -54 & -54 & -54 & -54 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{array}\right|$$

Create a 3D view of the house that will facilitate your understanding of how to place camera relative to the house geometry.   

def gHouseFace(vis) : 
    vrts = [[MMH1[i,j] for i in range(3)] for j in vis[0]]
    return polygon3d(vrts,color=vis[1],alpha=0.7)
def gAxe(pt,c) : 
    p1 = [x * -1 for x in pt]
    return line3d([p1,pt],thickness=5,color=c)
def axes() : 
    bnd = 75
    return map(gAxe,[(bnd,0,0),(0,bnd,0),(0,0,bnd)],['red','green','blue'])
figcon1 = sum(map(gHouseFace,faceV.values()) + axes())
figcon1.show(aspect_ratio=(1,1,1))

bnd = 10
pt  = [(bnd,0,0),(0,bnd,0),(0,0,bnd)]
[x * -1 for x in pt[0]]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-10, 0, 0\right]$$

In this simplified version there will be no camera to world rotation and the camera Eye is placed at the origin. This will make relating pseudo-depth to the placement of the far clipping plane easier to understand going between the symbolic and computational example. Also notice that the house has been moved in the XY plane so as to be centered on the Z-axis

var('umin, umax, vmin, vmax, near, far')
casecams = {
   'sym' : {
       'eye'  : vector(SR, 3, (0, 0, 0)),
       'look' : vector(SR, 3, (0, 0, 30)),
       'up'   : vector(SR, 3, (0, 1,  0)),
       'bnds' : { 'left' : umin, 'right' : umax, 'bottom' : vmin, 'top' : vmax},
       'nefa' : { 'near' : near, 'far' : far}
   },
   'ex1' : {
       'eye'  : vector(SR, 3, (0, 0, 0)),
       'look' : vector(SR, 3, (0, 0, 30)),
       'up'   : vector(SR, 3, (0, 1,  0)),
       'bnds' : { 'left' :-20, 'right' : 20, 'bottom' : -20, 'top' : 20},
       'nefa' : { 'near' : -25, 'far' : -75}
   }
}
case = 'ex1'

What follows is a complete illustration of how to construct the transformation matrices - ultimately a single matrix - that accomplishes 3D project of a model in the world in to a camera's canonical view volume. 

Without the need to translate or rotation from world to camera, we are left with two parts (note, rotation and translation are shown, however they are both identity transforms).  

1. Apply prespective projection to place 3D points into canoncal view volume, introducing perspective effect
2. Normalize camera X, Y and Z so to match standard -1 to 1 for X and Y and 0 to 1 bounds on cannonical view volume

cam  = casecams[case]  # Choose a symbolic or numeric example
TM = matrix.identity(SR, 4); 
RM = matrix.identity(SR, 4); 
OM = matrix.identity(SR, 4);
OM[0,0] = 2 / (cam['bnds']['right'] - cam['bnds']['left'])
OM[1,1] = 2 / (cam['bnds']['top']   - cam['bnds']['bottom'])
OM[2,2] = 2 / (cam['nefa']['near']  - cam['nefa']['far'])
OM[0,3] = - ((cam['bnds']['right'] + cam['bnds']['left'])   / (cam['bnds']['right'] - cam['bnds']['left']))
OM[1,3] = - ((cam['bnds']['top']   + cam['bnds']['bottom']) / (cam['bnds']['top']   - cam['bnds']['bottom']))
OM[2,3] = - ((cam['nefa']['near']  + cam['nefa']['far'])    / (cam['nefa']['near']  - cam['nefa']['far']))
PM = matrix(SR, 4,4, 0)
PM[0,0] = cam['nefa']['near']  
PM[1,1] = cam['nefa']['near']  
PM[2,2] = cam['nefa']['near'] + cam['nefa']['far']
PM[2,3] = - cam['nefa']['near'] * cam['nefa']['far']
PM[3,2] = 1
WCM = OM * PM;
pretty_print(LatexExpr("{M_{WC} \: = \: M_{O} \: M_{P}}"))
pretty_print(LatexExpr("{M_{P} \: = \:}"), PM)
pretty_print(LatexExpr("{M_{O} \: = \:}"), OM)
pretty_print(LatexExpr("{M_{WC} \: = \:}"), WCM)

$$\newcommand{\Bold}[1]{\mathbf{#1}}{M_{WC} \: = \: M_{O} \: M_{P}}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}{M_{P} \: = \:} \left|\begin{array}{rrrr} -25 & 0 & 0 & 0 \\ 0 & -25 & 0 & 0 \\ 0 & 0 & -100 & -1875 \\ 0 & 0 & 1 & 0 \end{array}\right|$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}{M_{O} \: = \:} \left|\begin{array}{rrrr} \frac{1}{20} & 0 & 0 & 0 \\ 0 & \frac{1}{20} & 0 & 0 \\ 0 & 0 & \frac{1}{25} & 2 \\ 0 & 0 & 0 & 1 \end{array}\right|$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}{M_{WC} \: = \:} \left|\begin{array}{rrrr} -\frac{5}{4} & 0 & 0 & 0 \\ 0 & -\frac{5}{4} & 0 & 0 \\ 0 & 0 & -2 & -75 \\ 0 & 0 & 1 & 0 \end{array}\right|$$

To focus on the pseudo-depth calculation, consider a symbolic point on the optical axes with value z

var('z', 'pz')
zpt = matrix(SR,[0, 0, 'z', 1]).transpose()
Mpz = WCM * zpt
for i in range(4) :
    Mpz[i,0] = Mpz[i,0] / Mpz[3,0]
eq1 = pz == Mpz[2,0] 
if (case == 'sym') :
    pretty_print(Mpz, LatexExpr("{\:=\:}"), WCM, zpt)
    pretty_print(LatexExpr("{P_{cc} \: = \:}"), Mpz, LatexExpr("{\;\;} \hbox{and the z term only} {\;\;} "), eq1)
else :
    f(z) = Mpz[2,0]
    near = cam['nefa']['near']
    far  = cam['nefa']['far']
    range(near, far, (far-near)/min(100,(near-far)))
    v = [(z, f(z)) for z in range(near, far, (far-near)/min(1000,(near-far)))]
    show(points(v, rgbcolor=(0.2,0.6, 0.1), pointsize=30) )

Once the single matrix is constructed, provided a numerical case is being run, the code below transforms vertices from world coordinates into camera coordinates. 

def can_points(M) : 
    for j in range(M.ncols()):
       for i in [0..3]:
          M[i,j] = M[i,j]/M[3,j];
if (case != 'sym') :
    MMH2 = WCM * MMH1
    can_points(MMH2)
    MMH2 = MMH2.apply_map(RR)
    pretty_print(MMH2.n(digits=2))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left|\begin{array}{rrrrrrrrrr} -0.33 & -0.33 & 0.00 & 0.33 & 0.33 & -0.19 & -0.19 & 0.00 & 0.19 & 0.19 \\ -0.33 & 0.083 & 0.33 & 0.083 & -0.33 & -0.19 & 0.046 & 0.19 & 0.046 & -0.19 \\ 0.50 & 0.50 & 0.50 & 0.50 & 0.50 & -0.61 & -0.61 & -0.61 & -0.61 & -0.61 \\ 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{array}\right|$$

if (case != 'sym') :
    def gHouseFace(vis) : 
        vrts = [[MMH2[i,j] for i in range(3)] for j in vis[0]]
        return polygon3d(vrts,color=vis[1],alpha=0.7)
    def axe2(pt,c) : 
        return line3d([-pt,pt],thickness=5,color=c,opacity=0.1)
    def axes2() : 
        return map(axe2,[vector([1,0,0]),vector([0,1,0]),vector([0,0,1])],['red','green','blue'])
    figcon2 = sum(map(gHouseFace,faceV.values()) + axes2())
    figcon2.show(aspect_ratio=(1,1,1))