Jelly Cubes

This is something I planned on programming some 30 years back on the Amiga and Assembly, but never ended up doing. Now did, with Python and Pygame. The idea here is to have two semi-transparent, co-centric cubes rotating so that one will grow out of the other, showing the parts outside the bigger cube being "cut". It's not realistic in any way but a nice exercise in vector graphics. Unfortunately, the routine is not perfect, either.

The code can be found in my github repository.

Rotating Back and Forth

It is quite possible to calculate if and how a 3D vector cuts through a 3D plane of any definition. However, it is rather much easier to calculate if and how a 3D vector cuts through a 3D plane if one of the 3D plane's coordinates is constant. For example, picture a cube where all coordinates (X, Y, Z) are -1 or +1; each of the six surfaces of that cube has one coordinate (X, Y or Z) which is constant (either -1 or +1). 

Then picture the other cube, which is smaller, but only a little - so that, when rotated, one or more of its vertices (or nodes or corners) will be outside the bigger cube. It will be quite simple to see when that happens - if and only if the vertex has a coordinate > 1 (or < -1). Also calculating the point where the edge vector cuts through the other cube's surface is simple - it is the point where that coordinate equals 1 (or -1).

There are two ways to achieve this. One is to first rotate the smaller cube, keeping the bigger cube stationary, then calculating the cuts, and rotating everything again to make also the bigger cube move. However, this causes difficulties when the cubes change size and the smaller cube at one point becomes the bigger one. The other option is to rotate the cubes, and then "rotate back" the smaller cube using the inverse of the big cube's rotation matrix. This will result in a version of the smaller cube which is rotated as if the bigger cube stayed stationary. Then the cuts can be calculated using this version, but applied directly to the "fully rotated" version - avoiding rotating the cut points.

    def rotate(self):

        if self.cube_sizes[0] < self.cube_sizes[1] and self.cube_small == 1:
            self.cube_small = 0
            self.cube_big = 1
        elif self.cube_sizes[0] > self.cube_sizes[1] and self.cube_big == 1:
            self.cube_small = 1
            self.cube_big = 0

        matrix_small = self.rotate_matrix_XYZ(self.angles[self.cube_small, :])
        matrix_big = self.rotate_matrix_XYZ(self.angles[self.cube_big, :])
        matrix_big_inv = np.linalg.inv(matrix_big)

        # rotate small cube
        self.rotated_nodes[0:8, :] = np.matmul(self.nodes * self.cube_sizes[self.cube_small], matrix_small)
        # rotate big cube
        self.rotated_nodes[8:16, :] = np.matmul(self.nodes * self.cube_sizes[self.cube_big], matrix_big)
        # make a special rotated copy of the small cube by rotating it with the inverse of big cube matrix.
        # the result is small cube rotated as if big cube was stationary (or not rotated at all).
        self.rotated_nodes_small = np.matmul(self.rotated_nodes[0:8, :] * self.cube_sizes[self.cube_small], matrix_big_inv)

    def rotate_matrix_XYZ(self, angles):

        # define rotation matrix for given angles
        (sx, sy, sz) = np.sin(angles)
        (cx, cy, cz) = np.cos(angles)

        # build a matrix for X, Y, Z rotation (in that order, see Wikipedia: Euler angles).
        return np.array([[cy * cz               , -cy * sz              , sy      ],
                         [cx * sz + cz * sx * sy, cx * cz - sx * sy * sz, -cy * sx],
                         [sx * sz - cx * cz * sy, cz * sx + cx * sy * sz, cx * cy ]])

What is needed for this "rotating back" is the inverse of the rotation matrix. Inverting a matrix can be a very costly operation and requires the matrix to be of a suitable form; fortunately, a 3 x 3 rotation matrix is rather quick and easy to invert.

Cutting Edge

In the simple case, a vertex protrudes the other cube's surface so that it forms a small pyramid, consisting of the three surfaces using that vertex, but cut at the +1 (or -1) level mentioned. (See the image above, the red vertex closest to the center of the image.) So we need to calculate the three cut points for the three edges using the vertex, which are simply the weighted averages of the edges' vertices, using the relative lengths of the edge within the cube and outside the cube as weights. So, if e.g. X = 1.1 for the vertex outside the cube, and X = 0.7 at the other end of the edge, the cut point (where X = 1.0) would be (0.1 / 0.4) x IN + (0.3 / 0.4) x OUT, where IN and OUT represent the inside and outside vertices, respectively. And, as said, the calculation can be applied directly to the fully rotated smaller cube vertices, even if the weights are calculated using the one which was "rotated back".

Of course, life is never simple. If the cube sizes are close to each other, there will be cases when the whole edge, i.e. both the vertices defining it and all points on the edge, are outside the bigger cube. (See the image above, the red edge at lower left side.) Then there is no cut point at all, and the "cut" is a much more complex area than a pyramid. The code handles this by building such cut surfaces separately. Also, the code calculates all the cut points in one go, utilising NumPy operations, so it is not quite as simple as the example above - but still it has to build all the cut surfaces in a loop, one by one.