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

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.

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

Crypto AG DID work with NSA / GCHQ; UNUMS (Universal Number Format); MOOCs (Massive Open Online Courses); Interesting Looking Game; Euler's Theorem for Polynomials; Weeks 3-6; Reddit Comment; Differential Cryptanalysis For Dummies; Japanese Graphic Design; Books To Be Re-Read; And Today I Learned Bugs Need Clear Examples; Factoring a 67 bit prime in your head; Islamic Geometric Art; Useful Julia Backtraces from Tasks; Nothing, however, is lost with less discomfort than that which, when lost, cannot be missed; Article on Didion; Cost of Living by City; British Slavery; Derrida on Metaphor; African SciFi; Traits in Julia; Alternative Japanese Lit; Pulic Key as Address (Snow); Why Information Grows; The Blindness Of The Chilean Elite; Some Victoriagate Links; This Is Why I Left StackOverflow; New TLS Implementation; Maths for Physicists; How I Am 8; 1000 Word Philosophy; Cyberpunk Reading List; Detailed Discussion of Message Dispatch in ParserCombinator Library for Julia; FizzBuzz in Julia w Dependent Types; kokko - Design Shop in Osaka; Summary of Greece, Currently; LLVM and GPUs; See Also; Schoolgirl Groyps (Maths); Japanese Lit; Another Example - Modular Arithmetic; Music from United; Python 2 and 3 compatible alternative.; Read Agatha Christie for the Plot; A Constructive Look at TempleOS; Music Thread w Many Recommendations; Fixed Version; A Useful Julia Macro To Define Equality And Hash; k3b cdrom access, OpenSuse 13.1; Week 2; From outside, the UK looks less than stellar; Huge Fonts in VirtualBox; Keen - Complex Emergencies; The Fallen of World War II; Some Spanish Fiction; Calling C From Fortran 95; Bjork DJ Set; Z3 Example With Python; Week 1; Useful Guide To Starting With IJulia; UK Election + Media; Review: Reinventing Organizations; Inline Assembly With Julia / LLVM; Against the definition of types; Dumb Crypto Paper; The Search For Quasi-Periodicity...; Is There An Alternative To Processing?; CARDIAC (CARDboard Illustrative Aid to Computation); The Bolivian Case Against Chile At The Hague; Clear, Cogent Economic Arguments For Immigration; A Program To Say If I Am Working; Decent Cards For Ill People; New Photo; Luksic And Barrick Gold; President Bachelet's Speech; Baltimore Primer; libxml2 Parsing Stream; configure.ac Recipe For Library Path; The Davalos Affair For Idiots; Not The Onion: Google Fireside Chat w Kissinger; Bicycle Wheels, Inertia, and Energy; Another Tax Fraud; Google's Borg; A Verion That Redirects To Local HTTP Server; Spanish Accents For Idiots; Aluminium Cans; Advice on Spray Painting; Female View of Online Chat From a Male; UX Reading List; S4 Subgroups - Geometric Interpretation; Fucking Email; The SQM Affair For Idiots; Using Kolmogorov Complexity; Oblique Strategies in bash; Curses Tools; Markov Chain Monte Carlo Without all the Bullshit; Email Para Matias Godoy Mercado; The Penta Affair For Idiots; Example Code To Create numpy Array in C; Good Article on Bias in Graphic Design (NYTimes); Do You Backup github?

© 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

Comment on this post