C[omp]ute

Projections

From: "andrew cooke" <andrew@...>

Date: Sat, 27 Dec 2008 08:43:40 -0300 (CLST)

The conversion from 3D (spatial model) to 2D points (paper plot) is worth
recording.  This is a summary of what I learnt from googling around - it
may not be correct (the info on the web about projective geometry isn't
completely consistent).

There are two common kinds of projections - orthographic and perspective.

Orthographic projections are calculated using an affine transformation - 
a linear calculation using a matrix multiplication (with a shift that can
be implemented using an additional vector addition, or by using an extra
dimension with "1" extending the vectors).

Orthographic projections can be further divided into three groups.  First,
those perpendicular to "principal" planes (like plans), which are not very
interesting here.  Second, axonometric projections which show a "skew"
view.  These are equivalent to a perspective projection "from infinity"
(so parallel lines remain parallel in projection) and the most common is
probably the isometric view.  Third, oblique projections where one plane
is shown "face on" while others are skewed.  The most common here is
"cavalier perspective" with an angle of 45 degrees and no scaling of
lengths.

My code for the orthographic projections is:

  def _orthographic(transform, offset):
    t = np.asarray(transform)
    o = np.asarray(offset)
    return lambda p: np.dot(t, p) + o

  def cavalier(offset=(0,0)):
    rsq2 = 1 / sqrt(2)
    return per_point(
      _orthographic(
	transform=((-1,rsq2,0), (0,-rsq2,1)),
	offset=offset))

  def isometric(offset=(0,0)):
    s32 = sqrt(3) / 2
    return per_point(
      _orthographic(
	transform=((-s32,s32,0),(-1./2,-1./2,1)),
	offset=offset))

where "per_point" is the function to convert from a pointwise transform to
one for the whole image (image is a list of paths; paths are lists of
points).


I couldn't find a simple description  of how to calculate a perspective
projection, but it's simple to derive by hand if you construct the
projection by taking the figure generated by rays from eye to object that
intersect a plane (the "image plane").

I imagine the image plane as a sheet hanging in space, through which I can
vaguely see the object that is to be projected; the perspective projection
is just tracing the outline seen.

For an "undistorted" view the image plane is perpendicular to the line of
sight to the "centre" of the object (ie the point on the object that will
be at the centre of the image).  Also, vertical lines in the object will
be vertical in the image.  This is sufficient to define the geometry if we
place the eye at E (capital letters for vectors) and the centre of the
image at C (so the image plane contains C and is perpendicular to EC).

I will use, for example, EC as shorthand for C-E - the vector from E to C
(the reversal of order between "EC" and "C-E" is a real pain).

There are two coordinate systems in use here, confusingly: Z is the
spatial vertical direction, while X and Y are the horizontal and vertical
directions in the image plane.  From the discussion above
("..undistorted...") we want Y to appear (from E) parallel to Z.  We can
use this to define X and Y in terms of C, E and Z.

First, let us define unit(V) = V / |V| (where |V| is the length of V, an
arbitrary vector) and N = unit(CE) (the normal to the image plane,
pointing towards E).

Then X = unit(Z^N) since it is perpendicular to N (in the image plane) and
perpendicular to the vertical.  I am using "^" for the "cross" or "vector"
product.

Similarly, Y = N^X

Now we can do the projection.  The eye looks towards a point P on the
object.  The line EP intersects the image plane at Q (there is no need for
Q to be between E and P - it may lie on the extended line, in which case
we get weird "fisheye" effects).

So the projection consists in finding Q for each P on the object and then
resolving Q along X and Y.

How to find Q?

We know that Q lies on EP, so for some parameter t, Q = E + t EP (a simple
parametrised line).  We also know that Q is in the image line and so is
perpendicular to N: 0 = N.Q (where "." is the "dot" or "scalar" product).

We can therefore solve for t:
Q = E + t EP
0 = N.Q
0 = N.(E + t EP)
-(N.E) = t (N.EP)
t = -(N.E) / (N.EP)

So Q = E - EP (N.E) / (N.EP)

The final projection to 2D coords is then just the dot product "along" X
and Y: (x, y) = (Q.X, Q.Y).

My code is a bit ugly, but seems to work:

  def perspective(eye, screen):
    def unit(v):
      return v / sqrt(np.dot(v, v))
    eye = np.asarray(eye)
    screen = np.asarray(screen)
    normal = eye - screen
    normal2 = np.dot(normal, normal)
    n = normal / sqrt(normal2)
    # screen axes
    x = unit(np.cross(Z, n))
    y = np.cross(n, x)
    def process(p1):
      los = p1 - eye
      t = -normal2 / np.dot(los, normal)
      p2 = eye + t * los # point in the plane
      return (np.dot(p2, x), np.dot(p2, y))
    return per_point(process)

Andrew