cormullion | juliabloggers.com

Re-posted from: https://cormullion.github.io/blog/2019/04/03/pentachoron.html

This page contains quite a few animated GIFs. I apologise in advance if your browser or network connection doesn’t like them. It’s much easier to add a GIF than it is to embed a video into a GitHub pages blog, unfortunately.

Pentachoron

image label

I encountered the word “pentachoron” recently for the first time, and I thought it had to be something cool. Turns out it is, so I thought I’d have a go at drawing one…

Unfortunately, it’s a 4-dimensional object, which is going to tax my limited 2D skills. But I like being out of my depth.

First steps

The first job is to define a couple of types to help organize all the various points.

Here’s an immutable structure to hold the coordinates of a point in four dimensions. In four dimensions, the axes are usually called x, y, z, and w.

struct Point4D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end

Point4D(a::Array{Float64, 1}) = Point4D(a...)

Base.size(pt::Point4D) = (4, )

Base.getindex(pt::Point4D, i) = [pt.x, pt.y, pt.z, pt.w][i]

By adopting the AbstractArray as our stepmother type, we can pick up many useful behaviors. size() and getindex() can be taught how to handle Point4Ds. If you browse the official documentation for a while you’ll find out how to effectively tap in to many desirable pre-defined abilities that our new type can inherit.

While we’re here, let’s do three-dimensional points as well:

struct Point3D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
end

Base.size(pt::Point3D) = (3, )

Conversion

The primary task we have to address is how to convert a 4D point into a 2D point. Let’s start with the easier task: how to convert a 3D point into a 2D point, ie how can we draw a 3D shape on a flat surface?

Consider a simple cube. The front face and the back face could have the same X and Y coordinates, and vary only by their Z values.

(This is just a diagram, it’s not really 3D…)

image label

using Luxor

@svg begin
    fontface("Menlo")
    fontsize(8)
    setblend(blend(
        boxtopcenter(BoundingBox()),
        boxmiddlecenter(BoundingBox()),
        "skyblue",
        "white"))
    box(boxtopleft(BoundingBox()),
        boxmiddleright(BoundingBox()), :fill)

    setblend(blend(
        boxmiddlecenter(BoundingBox()),
        boxbottomcenter(BoundingBox()),
        "grey95",
        "grey45"
        ))
    box(boxmiddleleft(BoundingBox()),
        boxbottomright(BoundingBox()), :fill)

    sethue("black")
    setline(2)
    bx1 = box(O, 250, 250, vertices=true)
    poly(bx1, :stroke, close=true)
    label.(["-1/1/1", "-1/-1/1", "1/-1/1", "1/1/1"],
        slope.(O, bx1), bx1)

    setline(1)
    bx2 = box(O, 150, 150, vertices=true)
    poly(bx2, :stroke, close=true)
    label.(["-1/1/0", "-1/-1/0", "1/-1/0", "1/1/0"],
        slope.(O, bx2), bx2, offset=-45)

    map((x, y) -> line(x, y, :stroke), bx1, bx2)
end 500 500 "images/pentachoron/cube.svg"

So the idea is to project the cube from 3D to 2D by keeping the first two values, and multiplying or modifying them by the third value. We’ll make a convert() function to do this:

function convert(Point, pt3::Point3D)
    k = 1/(K - pt3.z)
    return Point(pt3.x * k, pt3.y * k)
end

K is just a constant which provides a consistent value for depth:

const K = 4.0

Testing this quickly in Luxor.jl gives promising results:

@svg begin
    cube = Point3D[
        Point3D(-1, -1,  1),
        Point3D(-1,  1,  1),
        Point3D( 1, -1,  1),
        Point3D( 1,  1,  1),
        Point3D(-1, -1, -1),
        Point3D(-1,  1, -1),
        Point3D( 1, -1, -1),
        Point3D( 1,  1, -1),
    ]
    circle.(convert.(Point, cube) * 300, 5, :fill)
end 500 500 "images/pentachoron/basic3_to_2.png"

image label

It’s a simple type of perspective projection.

Using the same principle, let’s make a method for converting a 4D point.

function convert(Point3D, pt4::Point4D)
    k = 1/(K - pt4.w)
    return Point3D(pt4.x * k, pt4.y * k, pt4.z * k)
end

Let’s hope it works.

We can combine these into a single utility function called flatten() that takes a list of 4D points and double-maps them into a list of 2D points suitable for drawing.

function flatten(shape4)
    return map(pt3 -> convert(Point, pt3), map(pt4 -> convert(Point3D, pt4), shape4))
end

Going up a level

To test this, we’ll define the vertices of a unit pentachoron:

const n = -1/√5

const pentachoron = [Point4D(vertex...) for vertex in [
    [ 1.0,  1.0,  1.0, n],
    [ 1.0, -1.0, -1.0, n],
    [-1.0,  1.0, -1.0, n],
    [-1.0, -1.0,  1.0, n],
    [ 0.0,  0.0,  0.0, n + √5]]]

(According to Wikipedia, other names for the pentachoron include:

Regular 5-cell
C5
pentatope
pentahedroid
tetrahedral pyramid
4-simplex

I think pentachoron is the coolest, though.)

Here’s a list of “faces”, with each face defined by three of the vertices:

const pentachoronfaces = [
    [1, 2, 3],
    [1, 2, 4],
    [1, 2, 5],
    [1, 3, 4],
    [1, 3, 5],
    [1, 4, 5],
    [2, 3, 4],
    [2, 3, 5],
    [2, 4, 5],
    [3, 4, 5]]

A quick test:

@svg begin
    setopacity(0.2)
    pentachoron2D = flatten(pentachoron)
    for (n, face) in enumerate(pentachoronfaces)
        randomhue()
        poly(1500 * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end 600 250 "images/pentachoron/firstdraw.svg"

image label

This isn’t very interesting, although I think it’s correct. To see a more appealing display, let’s make it dance…

Swings and roundabouts

You usually rotate 2D points about a 1D point. You usually rotate 3D points about a 2D line (often one of the XYZ axes). So logically you’d rotate 4D points with reference to a 3D plane. We’ll define some matrices that do 4D rotations relative to a plane defined by two of the X, Y, Z, and W axes.

As I’m thinking about it, the XY plane is typically the plane of the drawing surface. If you think of the XY plane as a computer screen straight in front of you, the XZ plane is parallel to your desk or floor, and the YZ plane is like a wall beside your desk on your right or left side.

But what about the XW plane? And the YW and ZW planes, for that matter? This is the mystery of 4D shapes: we can’t see these planes, we can only imagine their existence by watching shapes moving through and around them.

Here’s a function that generates a suitable rotation matrix for rotating a 4D point in the XY plane:

function XY(θ)
    [cos(θ)   -sin(θ)  0        0;
    sin(θ)    cos(θ)   0        0;
    0         0        1        0;
    0         0        0        1]
end

and here’s another one that rotates a 4D point in the XW plane:

function XW(θ)
    [cos(θ)   0        0        -sin(θ);
    0         1        0        0;
    0         0        1        0;
    sin(θ)    0        0        cos(θ)]
end

Oh, it looks like we’re having a matrix party, so we might as well do the others while we’re here:

function XZ(θ)
    [cos(θ)   0        -sin(θ)  0;
    0         1        0        0;
    sin(θ)    0        cos(θ)   0;
    0         0        0        1]
end


function YZ(θ)
    [1        0        0        0;
    0         cos(θ)   -sin(θ)  0;
    0         sin(θ)   cos(θ)   0;
    0         0        0        1]
end

function YW(θ)
    [1        0        0        0;
    0         cos(θ)   0        -sin(θ);
    0         0        1        0;
    0         sin(θ)   0        cos(θ)]
end

function ZW(θ)
    [1        0        0        0;
    0         1        0        0;
    0         0        cos(θ)   -sin(θ);
    0         0        sin(θ)   cos(θ)];
end

We’ll make a handy utility function that uses one of these matrix functions to rotate an array of 4D points A in the plane defined by the matrix function:

function rotate4(A, matrixfunction)
    return map(A) do pt4
        Point4D(matrixfunction * pt4)
    end
end

So we can write rotate4(A, XW(π/2)) to rotate the array of points in A in the XW plane by π/2. We can also write rotate4(A, XW(π/2) * XY(π/2)) to rotate the points first in the XW plane and then in the XY plane.

Notice that the rotate4() function returns a new array, rather than modifying the original one. This probably has some disadvantages, but for now it means that we don’t have to keep track of mutating objects so carefully.

It moves

We’re now able to make an animation that rotates our pentachoron in one or more planes.

The frame() function generates a single frame, using the framenumber (eased_n provides a normalized version of the framenumber modulated by the easing function) to control the rotation between 0 and 2π. The scalefactor makes the unit shape big enough to fill the frame.

It’s helpful to color the faces very slightly, using poly(..., :fill)-er.

using ColorSchemes

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    sethue("black")
    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    pentachoron′  = rotate4(pentachoron, XZ(eased_n * 2π))
    pentachoron2D = flatten(pentachoron′)

    setopacity(0.2)
    for (n, face) in enumerate(pentachoronfaces)
        sethue(get(ColorSchemes.diverging_rainbow_bgymr_45_85_c67_n256,
            n/length(pentachoronfaces)))
        poly(scalefactor * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

And this function generates a series of frames and saves them to an animation:

function makemovie(w, h, fname;
        scalefactor=1000)
    movie1 = Movie(w, h, "4D movie")
    animate(movie1,
        Scene(movie1, (s, f)  -> frame(s, f, scalefactor),
            1:300,
            easingfunction=easeinoutsine),
        creategif=true,
        pathname="images/pentachoron/$(fname)")
end

makemovie(600, 600, "pentachoron-xz.gif", scalefactor=2000)

This is a rotation in the “parallel to the desk my computer’s resting on” XZ plane.

image label

It’s not a realistic rendering with hidden-surface removal—I wouldn’t know how to do that in 3D, let alone 4D. But it does make it slightly easier to follow some of the faces. Something to do one day might be to sort the “faces” to determine which are drawn first, but given the maths involved this is a challenge I’m dubious about tackling…

We can modify frame() to allow various permutations of XY, XW, XZ, YZ, YW, and ZW rotations. Here’s a version that rotates in XZ, followed by YW:

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    setopacity(0.2)

    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    pentachoron2D = flatten(
        rotate4(
            pentachoron,
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(pentachoronfaces)
        sethue(get(ColorSchemes.diverging_rainbow_bgymr_45_85_c67_n256,
            n/length(pentachoronfaces)))
        poly(scalefactor * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

makemovie(800, 800, "pentachoron-xz-yw.gif", scalefactor=2000)

image label

These are starting to make me think of kneading machines for making pizza dough.

You can draw a few different models at once, placing each one in a separate cell:

image label

This is showing 2D shadows of 3D shadows of 4D objects. Do these objects exist? I don’t know…

Other shapes

While the word “pentachoron” was unfamiliar to me until recently, the word “tesseract” is much more familiar. After all, didn’t the tesseract contain the Space Stone, one of the six Infinity Stones that apparently predate the universe and possess unlimited energy? (According to the scientists at Marvel Comics, at least…) It’s also the name of a 4-dimensional hypercube, one step up from the pentachoron.

Here are the vital numbers for the tesseract:

const tesseract = [Point4D(vertex...) for vertex in [
    [-1, -1, -1,  1],
    [ 1, -1, -1,  1],
    [ 1,  1, -1,  1],
    [-1,  1, -1,  1],
    [-1, -1,  1,  1],
    [ 1, -1,  1,  1],
    [ 1,  1,  1,  1],
    [-1,  1,  1,  1],
    [-1, -1, -1, -1],
    [ 1, -1, -1, -1],
    [ 1,  1, -1, -1],
    [-1,  1, -1, -1],
    [-1, -1,  1, -1],
    [ 1, -1,  1, -1],
    [ 1,  1,  1, -1],
    [-1,  1,  1, -1]]]

const tesseractfaces = [
    [1,   2,   3,   4],
    [1,   2,  10,   9],
    [1,   4,   8,   5],
    [1,   5,   6,   2],
    [1,   9,  12,   4],
    [2,   3,  11,  10],
    [2,   3,   7,   6],
    [3,   4,   8,   7],
    [5,   6,  14,  13],
    [5,   6,   7,   8],
    [5,   8,  16,  13],
    [6,   7,  15,  14],
    [7,   8,  16,  15],
    [9,  10,  11,  12],
    [9,  10,  14,  13],
    [9,  13,  16,  12],
    [10, 11,  15,  14],
    [13, 14,  15,  16]]

Replace pentachoron everywhere by tesseract in the above code, choosing the rotations to taste.

function frame(scene, framenumber, scalefactor=1000)
    background("black")
    setlinejoin("bevel")
    setline(10.0)
    setopacity(0.7)

    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    tesseract2D = flatten(
        rotate4(
            tesseract,
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(tesseractfaces)
        sethue([Luxor.lighter_blue, Luxor.lighter_green,
            Luxor.lighter_purple, Luxor.lighter_red][mod1(n, 4)]...)
        poly(scalefactor * tesseract2D[face], :fillpreserve, close=true)
        sethue([Luxor.darker_blue, Luxor.darker_green,
            Luxor.darker_purple, Luxor.darker_red][mod1(n, 4)]...)
        strokepath()
    end
end

makemovie(800, 800, "tesseract-xz-yw.gif", scalefactor=2000)

image label

(Just for fun I temporarily switched the color scheme over to the Julia logo colors. You can see the 3D rendering problems more clearly, as well!)

Here’s a composite of four pairs of rotations:

image label

Enough is enough

It’s tempting to go further, but these web pages do get very big with all these GIFs. But just one more… I quite like this one, the hexadecachoron:

const hexadecachoron = [Point4D(vertex...) for vertex in [
        [ 1,  0,  0,  0],
        [ 0,  1,  0,  0],
        [ 0,  0,  1,  0],
        [ 0,  0,  0,  1],
        [-1,  0,  0,  0],
        [ 0, -1,  0,  0],
        [ 0,  0, -1,  0],
        [ 0,  0,  0, -1]]]

using Combinatorics

const hexadecachoronfaces = combinations(1:8, 3)

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    setopacity(0.15)
    eased_n          =
        scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    hexadecachoron2D = flatten(
        rotate4(
            hexadecachoron,
                YZ(eased_n * 2π) *
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(hexadecachoronfaces)
        sethue(get(ColorSchemes.magma,
            n/length(hexadecachoronfaces)))
        poly(scalefactor * hexadecachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

makemovie(650, 650, "movie-hexadecachoron.gif", scalefactor=5000)

image label

There’s lots more fun to be had (combining two or more different shapes is fun), but I’m worried about your network and my brain cells, so that’s enough animated GIFs for today.

[2018-04-03]

cormullion signing off

This page was generated using Literate.jl.

juliabloggers.com

A Julia Language Blog Aggregator

Author Archives: cormullion

Pentachoron

Pentachoron

First steps

Conversion

Going up a level

Swings and roundabouts

It moves

Other shapes

Enough is enough

π day

π day