# C[omp]ute

Welcome to my blog, which was once a mailing list of the same name and is still generated by mail. Please reply via the "comment" links.

Always interested in offers/projects/new ideas. Eclectic experience in fields like: numerical computing; Python web; Java enterprise; functional languages; GPGPU; SQL databases; etc. Based in Santiago, Chile; telecommute worldwide. CV; email.

© 2006-2015 Andrew Cooke (site) / post authors (content).

## Cube - Series of Images for Laser Printer

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

Date: Fri, 26 Dec 2008 17:29:04 -0300 (CLST)

I've been messing around with Python (using NumPy for the calculations and
PyScript for the final Postscript rendering), generating a series of
images based on a cube that is "wrapped" by one or more lines.

Initially points on the surface of a cube are generated in 3D - this is
implemented using a Python generator.  Then they are filtered by a
"pipeline" of transforms that filter, project to 2D, and scale/rotate.

The initial generator takes as parameters further generators, one for each
line).  So, for example, the length of each line might be given by
"constant(10)", or "repeated_uniform(1, 5)" (repeated sampling from a
uniform distribution from 1 to 5), or [1, 2, 3] (lengths for 3 lines).

The transforms are not necessarily generators - they are functions that
are composed.  This is because some transforms require the entire set of
points to function correctly.  But in practice most, those that operate on
individual points, are implemented as generators.

http://www.acooke.org/cube.html

I'm wondering about printing some out on "nice" paper and seeing if I can
get a gallery to stock them.  Is that crazy?  I have no idea...

Andrew

### 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