Assignment   a3_perspective
By           Daniel B. Carr (Extended Feb 2004, Revised Mar 2007)
Concepts:    Homogeneous coordinates 
             Projection matrix and perspective
             Rotation 
             Horizon and vanish points
             Necker Cube, illusions, infinity points 
             Constructing a caricature of a building
             Shape channel and depth cues
             Matching lightness
             
Due:         Plots from 1.4, 2.5, 3.2-3.3, 3.4, 5.2 
             Color intensities from 5.2    

    
1. Projection

It took artists a long time to figure out perspective!
Now we can project 3D line drawings with ease.
   
1.1 A projection function and a projection matrix

Assume that we have already transformed our drawing 
to suitable units such as centimeters or inches.  

The equation for projecting to a plane z=0 with (0,0,d)
as the center of projection is (in right handed coordinates)
     
xp = d*x/(d-z)    
yp = d*y/(d-z)   
   
Here d is the viewing distance in inches.  The vectors x, y and
z give the 3-D coordinates of points to be projected.  The points
can be endpoints of line segments, vertices of triangles and so on. 
The resulting vectors, xp and yp, are the projected coordinates in
the plane z=0.
	
Different implementations for calculating projected coordinates have
different merits. A simple alternative for R is the function below.


##Run   create a projection function for a 3 column matrix: x,y,z 

proj2D = function(mat,d = 20){ #
  xp = d*mat[,1]/(d-mat[,3])
  yp = d*mat[,2]/(d-mat[,3])
  return(cbind(xp,yp))
}
## End

The computer graphics community often uses homogeneous
coordinates and multiplication by 4 x 4 transformation
matrix to handle rotation, scaling, translation, and projection.

The first step (after scaling the data) is to convert
Euclidean coordinates to natural homogeneous coordinates.
In other words augment the matrix with x, y, and z with a
fourth column of 1's 

The next step is to post multiple by the 4 x 4 projection
matrix given below. Here the notation [x y z 1] means
a matrix with the suggested column vectors.  

   [x' y' z' n'] = [x y z 1] %*% [ d 0 0  0 
                                   0 d 0  0
                                   0 0 0 -1
                                   0 0 0  d ]
   	  
Implementations using the transpose of the matrix
on each side of the equation are conceptually equivalent.
Remember that t(A %*% B) = t(B) %*% t(A) where t() is the
Transpose operator and %*% is matrix multiplication.  

Note that the 4 x 4 transformation matrix contains
no data coordinates and only one parameter, the viewing
distance d.  One matrix will transform all the coordinates.     
 
We see that x' = d*x
            y' = d*y
            z' = 0
            n' = d-z

Let x',y',z', and n' be the resultant column vectors whose 
values taken as rows give the natural homogeneous coordinates
of points to plot.  

The last step is to convert these coordinates back to 
corresponding Euclidean coordinates have a last coordinate of 1.  
Thus the Euclidean coordinates in the plane with zp=0 are given by 

   xp = x'/n' = d*x/(d-z)
   yp = y'/n' = d*y/(d-z) 

where / means division by corresponding elements.
This is the desired results from above


The remain parts 
   zp = z'/n' = 0/(d-z)
   np =  1    = (d-z)/(d-z)
are not needed.  
  

Yes there are many wasted operations such as multiplying
zero and adding zero, but the hardware speed more than
compensates for this. 

##Run:  define a project matrix for post multiplication 

projmat = function(d=24){
             mat = matrix(
                c(d, 0, 0,0,
                  0, d, 0,0,
                  0, 0, 0,-1,
	     	      0, 0, 0,d),
                  ncol=4,byrow=T)
                return(mat)
             }
##End

1.2 Railroad tracks

In a perspective view, the rails of a straight railroad track
converge as the tracks get far away.     
      
A primitive line drawing will illustrate this.  The rails are just
two lines and the ties are line segments connected and orthogonal to
the rails.  (This may look like gaps in a concrete sidewalk.) 

Let x coordinates for the left and right rails be at -1.5 and 1.5 feet.
Let y coordinates be -6 feet as if our eyes are 6 feet above the track.

(My youngest son, who has grown recently, would be sure to suggest
that I, a short person, should use a smaller number. [:-).  If we want
to be more specific, we could use our height minus 1/2 have our head
height. Our eyes are located about in the middle of our heads.    

Let our position be a z=20 in front of the project plane.  Let z
coordinates for the ties be in steps of 10 feet starting 10 feet
behind the plane z=0.  These coordinates will have negative values
in a right hand coordinate system.  We will use  -10*(1:n) with n=20.   

1.3 Track construction

The construction here is a bit tricky.  It will draw all n 
railroad ties with one lines() command.  It will draw both
of the rails with one lines command.  

The basic trick is to stack line segments separated by NA's
in x and y vectors.  An NA in either the x or y coordinate
signals the end of a poly-line.  We will stick the NA's in the x's
and the corresponding value for y will not matter.   

For x coordinates, each railroad tie will have three values c(1,-1,NA).
We will repeat this n times to draw n ties.  All the y's will be
6 feet deep.  We need 3*n of them to match the x's.  The z's will
have three -10 values followed by 3 -20 values and so on.  It is easy
to do this a rep() as indicated below.  Finally we need to add 1's for 
natural homogeneous coordinates.   

## Run

n = 20    # number of railroad ties  
x = rep(c(-1,1,NA),n)        
y = rep(-6,3*n)        # flat land
z =  rep(-10*(1:n),rep(3,n))
one = rep(1,3*n)
hmat = cbind(x,y,z,one)  # matrix of homogenous coordinates

## End

1.4 Projection with natural homogeneous coordinates and plotting

## Run
# projection matrix with our feet 20 feet in front of z=0 
projM = projmat(20) 
	
newhmat = hmat %*% projM          # project
xyp = newhmat[,1:2]/newhmat[,4]  # Transform x and y to Euclidean coordinates
xp = xyp[,1]                     # extra x
yp = xyp[,2]                     # extra y
 
windows()
plot(xp,yp,type='n',
     main="Very Thin Railroad Tracks",xlab='',ylab='')
box(col=1)

lines(xp,yp)  # ties
    
# extra points from the first and last railroad ties
last = length(xp)
subsleft=c(1,last-2)  
trackx = c(xp[subsleft],NA,xp[subsleft+1])
tracky = c(yp[subsleft],NA,yp[subsleft+1])
lines(trackx,tracky,lwd=2)

##End

# Yes, we could have combined xp with trackx and yp with tracky and
# and used one lines() command. 
 
Look at the plot.  Note that the rails appear to be converging at a point.
This point is called a vanishing point. It is located on the horizon.   
The y screen coordinate of vanishing point and the horizon is y=0.
This is amazingly (or obvious depending on your knowledge) the same as
as the y coordinate for our eyes in the construction.

   
2.  Rotation of rectangular parallelepipeds

This example does not use homogeneous coordinates.

2.1 The function below produces a rotation matrix
for post multiplication: dat %*% rotate3d()
where dat has n rows and three columns.  
                            
## Run
rotate3d = function(aboutx=20,abouty=25,aboutz=0){
   thetax = aboutx*pi/180  # convert to radians
   thetay = abouty*pi/180
   thetaz = aboutz*pi/180

   cx = cos(thetax); sx = sin(thetax)
   cy = cos(thetay); sy = sin(thetay)
   cz = cos(thetaz); sz = sin(thetaz)

   matXfixed = matrix(c(
          1,  0,  0,
          0, cx, sx,
          0,-sx, cx),ncol=3,byrow=T)

   matYfixed = matrix(c(
          cy,  0, -sy,
           0,  1,   0,
          sy,  0,  cy),ncol=3,byrow=T)

   matZfixed = matrix(c(
	    cz, sz, 0,
         -sz, cz, 0,
           0,  0, 1),ncol=3,byrow=T)

# note matrix multiplication is not commutative 
   return(matXfixed%*%matYfixed%*%matZfixed)
}

##End

2.2 Which way does the about rotation go about the axes

## Run test
# Does +x rotate about z into +y
x = c(1,0,0)
tmp = x %*% rotate3d(aboutx=0,abouty=0,aboutz=90)
round(tmp)

# Does +y rotate about x into +z
y = c(0,1,0)
tmp = y %*% rotate3d(aboutx=90,abouty=0,aboutz=0)
round(tmp)

# Does +z rotate about y into +x
z = c(0,0,1)
tmp = z %*% rotate3d(aboutx=0,abouty=90,aboutz=0)
round(tmp)

## End

Note 1: There are sign choices that determine which way a rotation goes
about an axis.  If rotate3d() does not meet an  established convention
at least we know what it does. 

Note 2:  Matrix multiplication is not commutative so the
order of rotating about axes makes a difference.  The order is
hard coded in the above function.  A more general function would
allow specification of the rotation order.  


2.3 Produce of vertices for drawing a centered
    rectangular parallelepiped (RP) with the coordinate axes
    orthogonal to the faces.

## Run:  A function to vertices
RPvertices = function(x=1,y=1,z=1){
                x=abs(x);y=abs(y);z=abs(z)
                mat = expand.grid(list(x=c(-x,x),y=c(-y,y),z=c(-z,z)))
                return(as.matrix(mat))
             }
          
rp = RPvertices()
rp

## End


2.4 Rotate with modest angles about x and y, 
    Project vertices and plot the row numbers
    at the locations


## Run
    
mat = rp %*% rotate3d(10,10,0) # Rotate

d = 15 # Set viewing distance

# quick project d*x/(d-z),  d*y/(d-z)
xp = d*mat[,1]/(d-mat[,3])
yp = d*mat[,2]/(d-mat[,3])

plot(xp,yp,type='n')
text(xp,yp,lab=1:length(xp))
##End

# Now add some edges by using vertex subscripts.
# Note that subscripting produces NA's when subscripts are NA's

#Run  
subs = c(1,2,4,3,1,NA,5,6,8,7,5,1,NA,3,7)
lines(xp[subs],yp[subs])
#END

2.5  Start fixing the 2.4 plot: Step 1

Turn off the axes (axes=F is a plot argument)
Use empty character strings for axis labels and main label.  
Add the subscripts to draw all remaining edges
Add a blue dot with pch= 19 at the project location of the 
     (-1,-1,-1) points so we know the lower left back point.

## Run Solution
plot(xp,yp,type='n',axes=F,xlab='',ylab='',main='')
#text(xp,yp,lab=1:length(xp))

subs = c(1,2,4,3,1,NA,5,6,8,7,5,1,NA,3,7,??)
lines(xp[subs],yp[subs])
points(xp[1],yp[1],pch=19,col="blue",cex=2)
#END

3. Make the projected cube that looks like the constructed cube

A projected cube may not look like the intended cube for two major reasons.
 
First the linear transformation of data coordinates to screen coordinates
depends on 
1) the actual the window size: The R windows device has problemes  
2) the internal margin space plus other reserved space
3) the computation that associates data limits with the plotting space  

The linear transformation may not treat the x and y axes equally.  
If axes are treated equally the cube be still be wrong size  

Second, our actual viewing distance may be different than the constuction
viewing distance leading to distorting in the z direction.       

3.1 Construct a cube with a 3 inch edge

## Run:  Create vertices of a 3-inch edge cube and project

cubeV = RPvertices(1.5,1.5,1.5) # V = Vertices
cubePV = proj2D(cubeV,d=24)         # PV = projected vertices

## End

3.2 Produce a suitable plot frame for projection and rotation

We will be viewing a cube with different rotations about the origin.
After rotation the diagonal from the cube center to a vertex could
lie along the x axis.  The length is 1.5*sqrt(3) = 2.598. 

Create a windows with a plot frame from [-2.7, 2.7] for both
x and y axes that actually has these measurements in inches.  
That is, one unit on each axis should measure 1 inch.   

The good enough solution for my screen resolution was

## Run

windows(width=6.3,height=7.3)
par(xaxs='i',yaxs='i',mai=c(.5,.5,.5,.5))
lims = c(-2.7,2.7)
plot(lims,lims,type='n',xlab='',ylab='')

## End

Your screen size and resolution is likely different than mine. 
Experiment with the window width and height until is pretty close
to the correct size.  If the cube is too big to work, 
you may use a smaller one.  
  
The axis style 'i' in Splus uses data limits as the plotting
frame limits and generates nice internal tick locations. I
am not positive that R uses the exact data limits.   

3.3  Plot the non-rotated cube with red points on the front square
     Vertices and smaller blue points at the back square vertices

## Run
lines(cubePV[subs,]) 
points(cubePV[5:8,1],cubePV[5:8,2],pch=19,col="red",cex=1.7)
points(cubePV[1:4,1],cubePV[1:4,2],pch=19,col="blue",cex=1.1) 

## End

Such lines drawing of a cube produce an optical illusion.  
This was first published in 1832 by Swiss crystallographer
Louis Albert Necker and the cubes are called Necker cubes. 

The illusion is that there are two stable interpretations and 
our minds can switch between them.  The red vertex square can appear
to be in front or in back.  The larger red dots should favor seeing
that square in front, but I still see it the other way. 

Our resolution of such depth ambiguity can be local.  The artist
MC Escher has exploited this to produce drawing in which water
appears to flow continuously down hill but the stream is forms
a connected loop. 

Note that there is only one infinity point this rotation. The point is
is at the center where parallel lines would converge.       

3.4  Three infinity points with two emphasized

In a rotated view of the cube the parallel edges, when extended in the 
correct direction, would meet at an infinity point.  There are
three sets of parallel edges and hence three infinity points.

A more extreme perspective make the convergence of parallel edges less
A chose rotation can two of the infinite point more obvious a the example of the the other infinity point.

## Run

windows(width=6.3,height=7.3)
par(xaxs='i',yaxs='i',mai=c(.5,.5,.5,.5))
lims = c(-2.7,2.7)
plot(lims,lims,type='n',xlab='',ylab='')

cubePRV = proj2D(cubeV%*%rotate3d(aboutx=35,abouty=6.5),d=5.5)

lines(cubePRV[subs,])     # assume subs is still above from you solution
                          # to 2.5. 

# extract x and y cooridinates for use in arrows
xPRV=cubePRV[,1]
yPRV=cubePRV[,2]

lines(cubePRV[c(3,4),],lwd=2)

# top blue arrows
subs1 = c(5,6)
subs2 = c(1,2)
arrows(xPRV[subs1],yPRV[subs1],xPRV[subs2],yPRV[subs2],col='blue',
       lwd=2, length=.1)

# bottom blue arrows
subs1 = c(7,8)
subs2 = c(3,4)
arrows(xPRV[subs1],yPRV[subs1],xPRV[subs2],yPRV[subs2],col='blue',
       length=.2,lwd=3)

# back downward red arrows
subs1 = c(3,4)
subs2 = c(1,2)
arrows(xPRV[subs1],yPRV[subs1],xPRV[subs2],yPRV[subs2],col='red',
        length=.1,lwd=2)

# lines with width cues for the intended view
lines(cubePRV[c(7,8),],lwd=5)
lines(cubePRV[c(5,6),],lwd=3)

# front downward red arrows
subs1 = c(7,8)
subs2 = c(5,6)
arrows(xPRV[subs1],yPRV[subs1],xPRV[subs2],yPRV[subs2],col='red',lwd=3)

## End

4. Rectangles in 3D with a constant x, or y, or z.

The following function produces rectangles with a normal that is 
parallel to a coordinate axis. The specification gives
two opposite vertices of the rectangle in 3D.  One coordinate
value of the two triples will have same value.  For example
(2,2,0) and (-3,-3,0) and the same depth, z=0.  The x and y
coordinates of opposite vertices are (2,2) and (-3,3).  

## Run
 
rect2d = function(v1,v2){
	 dat = expand.grid(list(x=c(v1[1],v2[1]),
	                        y=c(v1[2],v2[2]),
	                        z=c(v1[3],v2[3])))
	if(v1[3]==v2[3]){
	  subs = c(1,2,4,3)
	} else if(v1[2]==v2[2]){
	  subs = c(1,2,6,5)
   } else if(v1[1]==v2[1]){
	  subs = c(1,3,7,5)
	} else subs = c(1,2,4,3,1,
	              5,6,8,7,5,
	              1,3,7,8,
	              4,2,6)
	return(rbind(dat[subs,],c(NA,NA,NA)))
}

## End

5.  Lightness matching and perspective processing through the gray
    level channel. 

Research has produced various theories about human color processing.

Friedhoff and Benzon 1991 (Visualization, the second computer revolution,
W. H. Freemn and Company, New York) 

describe three processing channels.  One channel they call the
shape channel responds to gray level intensities, yields high resolution
shapes and processes depth cues.  Monocular or "pictoral" depth cues
include linear perspective, interposition or oculation, shadows, detail
perspective and aerial perspective.

(A more current resource such as Vision Science by Palmer may
provide more current knowledge and theories.)     

At a conference I once saw a remarkable demonstration supporting
the notion that linear perspective is processed through a gray
level channel.  It was a two color perspective line drawing shown using
a slide projector. There was one color for "background" and one for the
"foreground" lines.  The prespective draw hard to look at and 
appeared flat.  However it was clearly 3D when shown in black and white.   
The colors for the flat view had been carefully matched in gray level. 
The material below reflects efforts to produce a similar example.  
So far nothing I have not found has been as dramatic an example.

The phenomena seem fragile. It could be that slides provide much better color
control than working with monitors, or at least old monitors.  We don't
have control of the light purity from the red, green and blue phosphers
and we have limited control of the rgb intensities.
My monitor is better than in past so may there is still hope. 
At the very least, this exercise should help fix the concept in mind.

How can we match gray levels? 
The YIQ color transformation was develop to prevent oversaturated color
from bleeding across scan lines in television. As I remember more bits were
used to encode y, the gray level (or lightness) at the expense of
other color coordinates I and Q  The Y is computed from RGB
intensities on a [0,1] scale by  

Y = .299R + .587G + .114B

This should get us in the right ball park.  

5.1 Creating caricature building for a perspective view.

## Run

xmin = 8
xmax = 17
ymin = -1
ymax = 9
zmin = -5
zmax = 5

#______________front__________________
v1 = c(xmin,ymin,zmax)
v2 = c(xmax,ymax,zmax)
matf = rect2d(v1,v2)

# front door
v1 = c(xmin+3.5,ymin,zmax)
v2 = c(xmin+5.5,ymin+3,zmax)
matfd = rect2d(v1,v2)

# front windows
v1 = c(xmin+1.5,ymin+5,zmax)
v2 = c(xmin+2.5,ymin+6,zmax)
matf1 = rect2d(v1,v2)

v1 = c(xmin+4,ymin+5,zmax)
v2 = c(xmin+5,ymin+6,zmax)
matf2 = rect2d(v1,v2)

v1 = c(xmin+6.5,ymin+5,zmax)
v2 = c(xmin+7.5,ymin+6,zmax)
matf3 = rect2d(v1,v2)

v1 = c(xmin+1.5,ymin+8,zmax)
v2 = c(xmin+2.5,ymin+9,zmax)
matf4 = rect2d(v1,v2)

v1 = c(xmin+4,ymin+8,zmax)
v2 = c(xmin+5,ymin+9,zmax)
matf5 = rect2d(v1,v2)

v1 = c(xmin+6.5,ymin+8,zmax)
v2 = c(xmin+7.5,ymin+9,zmax)
matf6 = rect2d(v1,v2)

# _________________left face____________________

v1 = c(xmin,ymin,zmin)
v2 = c(xmin,ymax,zmax)
matl = rect2d(v1,v2)

# left window1
v1 = c(xmin,ymin+8,-4)
v2 = c(xmin,ymin+9,-3)
matl1 = rect2d(v1,v2)

# left window2
v1 = c(xmin,ymin+8,-2)
v2 = c(xmin,ymin+9,-1)
matl2 = rect2d(v1,v2)

# left window3
v1 = c(xmin,ymin+8,0)
v2 = c(xmin,ymin+9,1)
matl3 = rect2d(v1,v2)

# left window4
v1 = c(xmin,ymin+8,2)
v2 = c(xmin,ymin+9,3)
matl4 = rect2d(v1,v2)

# left window5
v1 = c(xmin,ymin+5,-4)
v2 = c(xmin,ymin+6,-3)
matl5 = rect2d(v1,v2)

# left window6
v1 = c(xmin,ymin+5,-2)
v2 = c(xmin,ymin+6,-1)
matl6 = rect2d(v1,v2)

# left window7
v1 = c(xmin,ymin+5,0)
v2 = c(xmin,ymin+6,1)
matl7 = rect2d(v1,v2)

# left window8
v1 = c(xmin,ymin+5,2)
v2 = c(xmin,ymin+6,3)
matl8 = rect2d(v1,v2)


# left window9
v1 = c(xmin,ymin+2,-4)
v2 = c(xmin,ymin+3,-3)
matl9 = rect2d(v1,v2)


# left window10
v1 = c(xmin,ymin+2,-2)
v2 = c(xmin,ymin+3,-1)
matl10 = rect2d(v1,v2)

# left window11
v1 = c(xmin,ymin+2,0)
v2 = c(xmin,ymin+3,1)
matl11 = rect2d(v1,v2)

# left window12
v1 = c(xmin,ymin+2,2)
v2 = c(xmin,ymin+3,3)
matl12 = rect2d(v1,v2)

bigmat = rbind(matf,matfd,matf1,matf2,matf3,matf4,matf5,matf6,
		matl,matl1,matl2,matl3,matl4,
		matl5,matl6,matl7,matl8,
		matl9,matl10,matl11,matl12)

pts = proj2D(bigmat)

## End   

5.2  Attempt to pick a different set of colors for foreground
and background the are distinguishable in hue, but that are so
close in gray level that the perspective does work. 
Below in 5.3 are some of my commented out attempts. 

Do not assume that the gray level calculation below is precise.  It 
is only provided as a rough guide.   

## Run

# good old black and white show perspective quite well
# modify the rgb intensities below
rC = c(0,1)
gC = c(0,1)
bC = c(0,1)

pen= rgb(rC,gC,bC)
.299*rC+.587*gC+.114*bC # gray

windows()
par(mai=c(.3,.3,.3,.3),bg=pen[2])
plot(c(6,22),c(ymin-.5,ymax+3),type='n',axes=F,xlab='',ylab='')
polygon(pts,density=0,col=pen[1],lwd=2)

## End


5.3  Some attempts 

#1 Not too good 
rC = c(1,0)
gC = c(0,.51)
bC = c(1,1)

#2 A little better
rC = c(.85,.3)
gC = c(.3,.57)
bC = c(1,1)


#3 More like 1
rC = c(1,.3)
gC = c(.3,.65)
bC = c(1,1)

#4 It is hard to see the lines
rC = c(.65,.3)
gC = c(.3,.475)
bC = c(1,1)

#5 Hard to see here too
rC = c(.60,.0)
gC = c(.0,.34)
bC = c(1,1)

#6 I an guessing the slide I saw was red and green
# Getting rid of blue here does makes the red stand out
# more obvious 
rC = c(.60,.0)
gC = c(.0,.34)
bC = c(0,0)


#7 Sticking in some blue seem to help 
rC = c(.60,.0)
gC = c(.0,.44)
bC = c(.4,.2)