Andrew Cooke | Contents | Latest | RSS | Twitter | Previous | Next


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.

Personal Projects

Lepl parser for Python.

Colorless Green.

Photography around Santiago.

SVG experiment.

Professional Portfolio

Calibration of seismometers.

Data access via web services.

Cache rewrite.

Extending OpenSSH.

C-ORM: docs, API.

Last 100 entries

[Link] Neat Python Exceptions; [Link] Fix for Windows 10 to Avoid Ads; [Link] Attacks on ZRTP; [Link] UK Jazz Invasion; [Review] Cuba; [Link] Aricle on Gender Reversal of US Presidential Debate; {OpenSuse] Fix for Network Offline in Updater Applet; [Link] Parkinson's Related to Gut Flora; Farellones Bike Park; [Meta] Tags; Update: Second Ride; Schwalbe Thunder Burt 2.1 v Continental X-King 2.4; Mountain Biking in Santiago; Books on Ethics; Security Fail from Command Driven Interface; Everything Old is New Again; Interesting Take on Trump's Lies; Chutney v6; References on Entropy; Amusing "Alexa.." broadcast; The Shame of Chile's Education System; Playing mp4 gifs in Firefox on Opensuses Leap 42.2; Concurrency at Microsoft; Globalisation: Uk -> Chile; OpenSuse 42.2 and Synaptics Touch-Pads; Even; Cherry Jam; Lebanese Writer Amin Maalouf; C++ - it's the language of the future; Learning From Trump; Chinese Writer Hu Fayun; And; Apricot Jam; Also; Excellent Article on USA Politics; Oh Metafilter; Prejudice Against The Rurals; Also, Zizek; Trump; Why Trump Won; Doxygen + Latex on CentOS 6; SMASH - Solve 5 Biggest Problems in Physics; Good article on racism, brexit, and social divides; Grandaddy are back!; Consciousness From Max Entropy; Democrats; Harvard Will Fix Black Poverty; Modelling Bicycle Wheels; Amusing Polling Outlier; If Labour keeps telling working class people...; Populism and Choice; Books on Defeat; Enrique Ferrari - Argentine Author; Transcript of German Scientists on Learning of Hiroshima; Calvert Journal; Owen Jones on Twitter; Possible Japanese Authors; Complex American Literature; Chutney v5; Weird Componentized Virus; Interesting Argentinian Author - Antonio Di Benedetto; Useful Thread on MetaPhysics; RAND on fighting online anarchy (2001); NSA Hacked; Very Good LRB Article on Brexit; Nussbaum on Anger; Tasting; Apple + Kiwi Jam; Hit Me; Sudoku - CSP + Chaos; Recycling Electronics In Santiago; Vector Displays in OpenGL; And Anti-Aliased; OpenGL - Render via Intermediate Texture; And Garmin Connect; Using Garmin Forerunner 230 With Linux; (Beating Dead Horse) StackOverflow; Current State of Justice in China; Axiom of Determinacy; Ewww; Fee Chaos Book; Course on Differential Geometry; Okay, but...; Sparse Matrices, Deep Learning; Sounds Bad; Applebaum Rape; Tomato Chutney v4; Have to add...; Culturally Liberal and Nothing More; Weird Finite / Infinite Result; Your diamond is a beaten up mess; Maths Books; Good Bike Route from Providencia / Las Condes to Panul; Iain Pears (Author of Complex Plots); Plum Jam; Excellent; More Recently; For a moment I forgot StackOverflow sucked; A Few Weeks On...; Chilean Book Recommendations; How To Write Shared Libraries

© 2006-2017 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.

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



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

My code for the orthographic projections is:

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

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

  def isometric(offset=(0,0)):
    s32 = sqrt(3) / 2
    return per_point(

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

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"

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(, v))
    eye = np.asarray(eye)
    screen = np.asarray(screen)
    normal = eye - screen
    normal2 =, 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 /, normal)
      p2 = eye + t * los # point in the plane
      return (, x),, y))
    return per_point(process)


Comment on this post