Texture Mapping on a Sphere

I covered simple texture mapping - affine mapping on a 3D rotated plane - in my previous post. Mapping a texture on a sphere is of course more difficult. In this post I present two Python solutions, which both are less than perfect. Also, you may want to take a look at the RGB Sphere I wrote about earlier, with more spherical 3D graphics, using polygons - it maps a world map on a sphere as well, with light sourcing, although the resolution is then quite low.

As usual, the code can be found in my github repository.

Twins - but not identical

The two mapping methods are shown below.

The first one is based on mapping each and every pixel on the sphere to a flat Mercator map of the Earth, while the second is actually rotating a lot of (ca. 165 000) 3D dots, which have been assigned a color using the same map. In the first, then, the dots are constant but their colors change; in the second, the colors are constant and the dots move...

The main drawback of the first option is that as the dots are constant, the size of the sphere is constant as well. And the problem with the second is that even that big amount of independent dots is not enough to avoid some parts of the sphere not being covered, leaving black dots between the colored ones. See the video to understand this better:

Going Forwards

The second option, rotating a lot of 3D dots, has rather little to do with texture mapping actually - it is a relatively simple exercise of 3D rotation once the dots have been set up with their respective colors.

        # generate 3D nodes for a sphere.

        a = np.pi / 180.0  # for scaling degrees to radiuses
        R = self.image_mercator_R

        # first generate "traditional" 3D nodes with ~even spacing on sphere surface
        c = int(self.point_nr / 4)  # nr of circles of points needed for half a sphere.
        for i in range(c):
            lat = 90.0 / (2 * c) + 90.0 * (i / c)  # latitude north
            rad = np.cos(lat * a)  # radius scale at this latitude
            p = int(self.point_nr * rad)  # nr of points needed at this latitude
            j = np.arange(p)[:, None]
            long = 360.0 / (2 * p) + 360.0 * (j / p)  # longitudes at this latitude; most at the equator, least close to pole(s)
            x = np.cos(long * a) * rad * self.radius
            y = np.ones((p, 1)) * np.sin(lat * a) * self.radius  # y = latitude is constant
            z = np.sin(long * a) * rad * self.radius
            # add nodes both north and south
            self.nodes = np.vstack((self.nodes, np.hstack((x, y, z)), np.hstack((x, -y, z))))

            # pick colors from image
            if self.image_type == 'Mercator':
                # Mercator mode (https://en.wikipedia.org/wiki/Mercator_projection). Crop at 0.9999 - R defines "maximum" latitude (85 deg)
                lat_coord = np.minimum(0.9999, np.ones((p, 1)) * ((1.0 + R * np.log(np.tan(np.pi / 4.0 + 0.4999 * (lat / 90.0) * np.pi / 2.0)))
                                                                  / 2.0)) * self.image_size[1]
            else:
                # Normal mode: picture is simply wrapped
                lat_coord = np.ones((p, 1)) * (lat / (90.0 * 2.0) + 0.5) * self.image_size[1]

            image_coord = (np.vstack((np.hstack((long / 360.0 * self.image_size[0], lat_coord)),
                                      np.hstack((long / 360.0 * self.image_size[0], self.image_size[1] - lat_coord - 1))))).astype(np.int16)
            self.node_colors = np.hstack((self.node_colors, self.image_array[image_coord[:, 0], image_coord[:, 1]]))

Setting up those 165 000 (self.point_nr actually depends on the window height and represents the number of dots on the sphere equator) dots as evenly as possible on a sphere involves going through latitudes (lat, based on i in the loop for i in range(c)). Number of points p is small for high latitudes and big when close to the equator, to get even spacing. lat and long (which are in degrees) are then simply converted to 3D coordinates x, y, and z; and these are used to add both the north and south (where y is negative) hemisphere nodes (dots).

Then, for each node/dot, a color must be assigned. When using a Mercator map, the conversion of latitude to image y coordinate is slightly complex; the "normal" mode is simpler and works better for any normal image. (Like a photo of your cat, for example. If using such a photo I recommend making an image with two copies of the photo side by side, like the eastern and western hemisphere.) Once the image coordinates have been calculated, the node colors are stored so that there is a color defined for each 3D node. 

The program then simply rotates these nodes and uses the colors mapped above to show the ones which are on the visible side (the "front") of the sphere. (See my earlier postings on 3D rotation if interested.)

And Backwards

The actual texture mapping is based on a fixed circular area, where each pixel is mapped to an image so that the imaginary sphere appears to be rotating in 3D space.

        for y_val in range(-self.radius + 1, 0):
            rad = int(np.sqrt(self.radius ** 2 - y_val ** 2))  # radius (pixels) at this latitude
            x = np.arange(-rad, rad + 1)[:, None]  # x simply covers all pixels on this line (y = constant)
            y = np.ones((rad * 2 + 1, 1)) * y_val  # y = constant on this line
            z = -np.sqrt(self.radius ** 2 - y_val ** 2 - x ** 2)  # z from sphere formula x**2 + y**2 + z**2 = r**2; negative as facing the viewer
            self.nodes_flat = np.vstack((self.nodes_flat, np.hstack((x, y, z)), np.hstack((x, -y - 1, z))))

        # store as integers, as anyway plotting pixels
        self.nodes_flat_x = (self.nodes_flat[:, 0]).astype(np.int16)
        self.nodes_flat_y = (self.nodes_flat[:, 1]).astype(np.int16)
        # precalculate [Mercator] conversion for each (integer) y and the distance between x's. Add precision by using 4 times the number of radius lines
        flat_y_range = np.arange(-self.radius * 4.0, self.radius * 4.0 + 1)
        self.image_flat_x = (self.image_size[0] / 4.0) / np.maximum(2.0, np.sqrt(self.radius ** 2 - (flat_y_range / 4.0) ** 2))
        if self.image_type == 'Mercator':
            self.image_flat_y = (np.minimum(0.9999, np.maximum(0.0, ((1.0 + R * np.log(np.tan(np.pi / 4.0 + 0.49999 * (np.arccos(flat_y_range /
                        (4.0 * self.radius))- np.pi / 2.0)))) / 2.0))) * self.image_size[1]).astype(np.int)
        else:
            self.image_flat_y = (np.minimum(0.9999, np.arctan2(flat_y_range, np.sqrt((self.radius * 4.0) ** 2 - flat_y_range ** 2)) / np.pi + 0.5)
                                 * self.image_size[1]).astype(np.int)

        self.rotated_nodes_flat = np.zeros((np.shape(self.nodes_flat)[0], 4), dtype=np.float)

This time, no longitudes and latitudes, but pixels to start with! That means going through half a circle vertically (y_val) and each x on that y_val horizontal line. Then, a z is calculated from these so that they (x, y and z) again represent 3D coordinates on a sphere, but in such a way that when plotted in 2D they will fill a circle with each dot representing one and exactly one pixel. These are stored as nodes_flat. 

Then, pre-calculated image mapping data are calculated and stored, as log, arccos, arctan etc. are rather cumbersome to calculate. This will help mapping later.

Now, normal 3D rotation - as in Going Forwards above - simply takes a stationary object and rotates its nodes around all three axes, resulting in a new set of nodes in 3D space. The Backwards part here means that we already know the end result i.e. the rotated nodes in 3D space - the nodes_flat above - and we instead need to know where they came from. As the nodes_flat are so defined that they neatly fill the circle we are using to plot the sphere, we can nicely map the texture this way.

    def rotate(self):

        # rotate object
        self.angles += self.rotate_speed
        matrix = self.rotate_matrix(self.angles)

        if self.mode == 1:
            # normal 3D rotation
            self.rotated_nodes[:, 0:3] = np.matmul(self.nodes, matrix) * self.size
        else:
            # invert rotation matrix - "backwards rotation" of the required 2D nodes (flat nodes) to "original" to get colors
            matrix_inv = np.linalg.inv(matrix)
            self.rotated_nodes_flat[:, 0:3] = np.matmul(self.nodes_flat, matrix_inv)

        self.measure_time("rotate")

    def rotate_matrix(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 ]])

The difference in rotation is deceptively small. Instead of using the usual rotation matrix directly, I invert it first using NumPy's linalg.inv. (If your matrix algebra needs refreshing, you can think of multiplying something with an inverse of a matrix in a same way as multiplying something with 1/x instead of x - it's a lot like dividing instead of multiplying.) 

Then, when we know where the end result nodes came from, we still need to map them to the image.

        rgb_array = pygame.surfarray.pixels2d(self.screen)
        rgb_array[self.nodes_flat_x + self.position[0], self.nodes_flat_y + self.position[1]] = self.image_array[
            ((np.arctan2(self.rotated_nodes_flat[:, 0], self.rotated_nodes_flat[:, 2]) / (2.0 * np.pi) + 0.5)
             * self.image_size[0]).astype(np.int16),  # x calculated
            # alternative X using a precalculated coordinate - but too complex to help much
            # (self.image_flat_x[(4.0 * self.rotated_nodes_flat[:, 1]).astype(np.int16) + 4 * self.radius] * self.rotated_nodes_flat[:, 0]
            #                    + (np.sign(self.rotated_nodes_flat[:, 2]) + 2) * self.image_size[0] / 4).astype(np.int16),  # x using precalcs
            # y calculated from Mercator - rather complex and hence slow
            # (np.minimum(0.9999, np.maximum(0.0, ((1.0 + R * np.log(np.tan(np.pi / 4 + 0.4999 * (np.arccos(self.rotated_nodes_flat[:, 1] / self.radius)
            #                 - np.pi / 2)))) / 2.0))) * self.image_size[1]).astype(np.int)  # y calculated
            self.image_flat_y[(4.0 * self.rotated_nodes_flat[:, 1]).astype(np.int16) + 4 * self.radius]  # y using precalcs
           ]

That is actually rather simple, or at least short, above. I use pygame's surfarray.pixels2d and set the nodes_flat pixels centered around self.position to a color found from self.image_array (another surfarray.pixels2d based on the map image). The version used above calculates the image x coordinate using self.rotated_nodes_flat, and picks the y coordinate using the pre-calculated self.image_flat_y data described above. 

The commented (#) alternative way uses pre-calculation for x as well, but it is too complex to really help; and it calculates y (Mercator only) but this is quite a bit slower than the pre-calculated alternative.

And that's all I have to say about that.

Texture Mapping (Tutorial)

While this blog has featured mostly old skool vector graphics, which were very popular some 30 years back as there was basically no hardware support for texture mapping and CPU power was what it was, I decided to try that using the usual suspects Python, NumPy and pygame. Of course, handling graphics pixel by pixel is much harder (and slower) than drawing one-color surfaces. Even so, the first real time texture mapping I am aware of (in computer graphics that is) dates as far as September 1991: the Cube-o-Matic Amiga demo by Spreadpoint. Quite a feat, I must say, even if the resolution is only about 100 x 100 pixels and updated perhaps 10-15 frames per second.

My basic PC handles about 5 million pixels per second in this implementation, perhaps 50 times more than Cube-o-Matic - a very modest improvement considering the resources. I used images from the Sound Vision demo I pythonized some time ago: credits for some of the images go to Red Baron and Overlander, and Jellybean for the music. 

As you can see from the video above, images can be mapped on any (convex) flat surface in 3D space. Besides the cube (squares), they are mapped on the faces of a dodecahedron (pentagons) and an icosahedron (triangles). 

The Python script and the images can be found in my github repository. If you run the program and wish to control the object's rotation, use cursor keys and a and z keys to rotate around the three axes.

Reading a Map, Affinely

To map a texture (image) to a 3D rotated flat surface one needs to form a new image, of different size and shape, where each pixel's color corresponds to the color of the source image's pixel at the same position. There are many advanced methods, and for advanced, read: slow and hard to implement; I chose perhaps the simplest way, which is called Affine Texture Mapping. This, according to Wikipedia, is the fastest form of texture mapping and was implemented already in e.g. the original PlayStation (1994) hardware.

Affine mapping basically means linearly interpolating the pixels on the rotated surface to the original image, and picking a color for each pixel. There are ways to correct for perspective distortions or sample the original image in a more intelligent way, but they, in this implementation and in my opinion, would bring very little benefit compared to the extra processing cost.

Linear interpolation is easy between two coordinates. For a square, even if rotated to any angle, it would be easy to linearly interpolate between the pixels on opposite sides, although due to perspective the other side might be shorter than the other. However, interpolating on lines not strictly horizontal or vertical would result in some pixels being processed many times and/or some pixels not at all. Hence, it is a much better approach to process the rotated surface line by line (horizontally or vertically), as then each pixel can be processed exactly once. I chose starting from the highest corner of the square and processing lines going down. But then one ends up with corners on the way - the process must be broken down to pieces (or segments) where there are none. 

How to do Texture Mapping

The image and video below show the tutorial mode of the program (if self.tutorial_mode is set to false, it will run like the video above). First of all, for each surface (face of the cube) I have defined the picture to be used, and also, for each node (corner) of the surface, the respective coordinate within the picture. 

See the blue triangle, red trapezoid (a trapezium in Britain), and green triangle in the image? They represent the segments which can each be processed by linear interpolation. As said I start processing the rotated surface from the top; starting point is the top of the blue triangle. Then, I can go down as far as the first corner on the left. Note that the triangle is strictly horizontal at the bottom, representing a line of pixels. To be sure to process each pixel, it is necessary to go down the blue triangle horizontal line by line, in effect increasing the y coordinate one by one, and processing all the x coordinates between the two triangle sides.

To understand the mapping, it is perhaps helpful to look at how the cube object is defined: If you wish to look into how I define and rotate 3D vector objects in more detail, see my earlier series on Real-Time 3D Vector Graphics with Python.

            # define a cube.

            self.nodes = (100.0 * np.array([
                [ 1,  1,  1],
                [ 1,  1, -1],
                [-1,  1, -1],
                [-1,  1,  1],
                [ 1, -1,  1],
                [ 1, -1, -1],
                [-1, -1, -1],
                [-1, -1,  1]
                ])).astype(np.float)
            self.surfaces = np.array([
                [0, 1, 2, 3],
                [4, 7, 6, 5],
                [0, 4, 5, 1],
                [2, 6, 7, 3],
                [1, 5, 6, 2],
                [3, 7, 4, 0]
                ])

The cube's nodes are simply 100 units from the origo in each axis; and the six surfaces are defined by the four nodes at their corners. For example, the first surface is defined by the first four [0, 1, 2, 3] nodes, each node having +100 as its y coordinate (coordinates are x, y, z). The second surface is the opposite surface: nodes [4, 7, 6, 5] have y coordinate -100. (The reason for not using [4, 5, 6, 7] is that by defining surface nodes clockwise enables later calculating if they are "facing" towards the viewer (and should be drawn) or not.)

    def map_texture(self, surface):

        # build a node array where the nodes appear twice - enabling going "over" the right edge
        nodes_x2 = np.hstack((self.surfaces[surface, :], self.surfaces[surface, :]))
        # find the topmost node (minimum y cooridnate) and maximum y coordinate
        min_y_node = np.argmin(self.trans_nodes[self.surfaces[surface, :], 1])
        max_y = np.max(self.trans_nodes[self.surfaces[surface, :], 1])
        y_beg = self.trans_nodes[nodes_x2[min_y_node], 1]
        # when going "left" and "right" through the nodes, start with the top node in both
        (left, right) = (min_y_node, min_y_node)

The map_texture function takes surface as its only argument, so for the cube object this is between 0 and 5. The rotated surface's 2D screen coordinates are in self.trans_nodes. First, I look up the highest corner min_y_node from the surface's trans_nodes (which are a pair of (x, y) coordinates, hence subscript 1). Then I get the lowest and highest corner y levels into max_y and y_beg and initialize the left and right node pointers with min_y_node (ie. the highest corner).

Even before the above I have set up a nodes_x2 variable, which is simply the surface's nodes twice in one NumPy array. As an example, for surface[0], its nodes self.surfaces[0, :] would be [0, 1, 2, 3], and nodes_x2 would then be [0, 1, 2, 3, 0, 1, 2, 3]. The reason? Let's say the min_y_node is node 3. Usually the lowest node with max_y is then node 1, the opposite corner on that face. As left and right start from node 3, I can now get to node 1 both by decreasing left and increasing right without checking whether I am already at the maximum number of items in the array.

Next, I loop through each section (in the image, the blue triangle, the red trapezoid, and the green triangle):

        # loop through each section from this y coordinate to the next y coordinate until all sections processed
        while y_beg < max_y:
            # image node depends on the node order
            img_node_beg = image_nodes[np.array([left % image_node_cnt, right % image_node_cnt]), :]
            img_node_end = image_nodes[np.array([(left - 1) % image_node_cnt, (right + 1) % image_node_cnt]), :]
            img_node_diff = img_node_end - img_node_beg
            # cube node comes from surface node list
            node_beg = self.trans_nodes[nodes_x2[np.array([left, right])], :]
            node_end = self.trans_nodes[nodes_x2[np.array([left - 1, right + 1])], :]
            node_diff = node_end - node_beg

Using the left and right node pointers, I set up img_node_beg, which contains the image coordinates to the section's top nodes on the left and right sides (in this case the same node). In the image below it is node 3. Note that for other surfaces the surface nodes could be e.g. [4, 7, 6, 5], while there are always just four image nodes for cube surfaces - that is why I am using the left and right pointers with modulo and image_node_cnt (which is 3 for triangles, 5 for pentagons etc.). Similarly, the section's bottom nodes' coordinates on the left and right sides are set up in img_node_end: in this example, they are 2 and 0, respectively - going left and right from the top node 3 in[0, 1, 2, 3, 0, 1, 2, 3].

The same is done using the coordinates of the rotated surface i.e. self.trans_nodes. Here, node_beg and node_end contain the screen coordinates of nodes 3 and 3 & 2 and 0. The img_node_diff and node_diff simply represent the vectors on the left and right side for both. 

Our section begins at y_beg.Then we need to know where it ends. That is simple:

            # find section end = y_end
            if node_end[1, 1] < node_end[0, 1]:
                # right node comes first (i.e. its Y coordinate is before left's)
                right += 1
                y_end = node_end[1, 1]
            else:
                # left node comes first (i.e. its Y coordinate is before or equal to right's)
                left -= 1
                y_end = node_end[0, 1]

Once we know which one comes first, the next node on the right or left, we can also increase the respective pointer, because that is the side we need to move forward for the next section.

Now we can get working on this section on a line by line level.

Section, Part, Portion, Piece, Segment, Fragment

Our blue triangle section has a left side beginning at node 3 and ending at node 2, and a right side beginning at node 3 but ending at point m somewhere between nodes 3 and 0. Of course, we know where m is, as it must have the same y coordinate as node 2 to keep our lines horizontal. Let's move forward.

            if y_end > y_beg:
                y = np.arange(y_beg, y_end, dtype=np.int16)
                # node multipliers for each y for left and right side. Since y_end is the first node down, node_diff[:, 1] is here always > 0
                m = (y[:, None] - node_beg[:, 1]) / node_diff[:, 1]
                # respective screen x coordinates for left and right side
                x = (np.round(node_beg[:, 0] + m * node_diff[:, 0])).astype(np.int16)
                x_cnt = np.abs(x[:, 1] - x[:, 0])  # + 1   - use +1 when using linspace method below
                # count cumulative pixel count to use as the offset when storing data
                x_cnt_cum = np.hstack((np.array([0]), np.cumsum(x_cnt))) + self.scr_cnt
                # respective image coordinates, interpolating between image nodes (usually corners)
                img_l = img_node_beg[0, :] + m[:, 0:1] * img_node_diff[0, :]
                img_r = img_node_beg[1, :] + m[:, 1:2] * img_node_diff[1, :]

First, there needs to be a check if the section has zero lines. This would happen if e.g. the top two nodes of the rotated surface would have the same y coordinate. If not, then I build the section's y coordinates into y - for example, if y_beg = 100 and y_end = 120 then y would have the 20 integers starting at 100 and ending at 119. The multipliers m are then calculated for each y and for both left and right side. In the blue triangle, the left multiplier would start from zero and end at 0.95 (19 / 20), but the right multiplier would (while starting from zero as well) end at somewhere around 0.4, as that is where point m lies on the way from node 3 to node 0.

Now we know the y coordinates for each horizontal line, and can calculate the respective left and right x coordinates using the multipliers. From those we can get the number of pixels on each horizontal line (x_cnt) and the cumulative count x_cnt_cum to use as a data pointer when storing the results. 

The image coordinates, obviously, do not lie on a horizontal line. We need to calculate the left and right end of the image data vector our horizontal line above represents, and that can be calculated using the image nodes and the same multipliers. 

Going Pixels

Now we have the screen coordinates of each horizontal line in y and x, and the respective image vector beginning and end coordinates at img_l and img_r. It's time to go through each line and figure out the surface color pixel by pixel.

                for y_line in range(y_end - y_beg):
                    # process each horizontal line, these are the x coordinates from x left to x right
                    if x_cnt[y_line] > 1:
                        # if "left" not on the left side, use negative step.
                        scr_x = np.arange(x[y_line, 0], x[y_line, 1], np.sign(x[y_line, 1] - x[y_line, 0]), dtype=np.int16)
                        # add x coordinates to self.scr_x array
                        self.scr_x[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1]] = scr_x
                        # add y coordinates similarly - y is constant
                        self.scr_y[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1]] = y_line + y_beg
                        # interpolate between line begin and end coordinates in image
                        self.img_xy[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1], :] = (img_l[y_line, :] + ((scr_x - scr_x[0]) / (scr_x[-1] - scr_x[0]))[:, None] * (img_r[y_line, :] - img_l[y_line, :])).astype(np.int16)
                        # store the color found in each interpolated pixel in self.scr_col
                        self.scr_col[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1]] = image_array[self.img_xy[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1], 0], self.img_xy[x_cnt_cum[y_line]:x_cnt_cum[y_line + 1], 1]]

Note that our "left" coordinate may not actually be on the left; it is more the direction of going through the nodes, and they are constantly rotated. Anyway, scr_x represents all the x coordinates on the horizontal line. These and the y coordinates are stored in self.scr_x and self.scr_y.

As we know the image data vector, i.e. the line within the image we want to project on this horizontal line, it is simply another interpolation required to get the same number of coordinates (as in scr_x) on that vector, evenly spaced. Those coordinates are stored in self.img_xy. 

And finally, we need the colors of those image pixels. They are stored in self.scr_col and picked from image_array, which is a pygame.surfarray.pixels2d object representing the image. 

That's it - the section has been processed and the screen x, y, and colors stored. Then it's time to set y_beg = y_end and process the next section. For the red trapezoid, the node_beg and node_end will now contain the screen coordinates of nodes 2 and 3 & 1 and 0 as the left side has moved forward while the right side has not, and y will be between the y coordinates of nodes 2 and 0.

When all surfaces have been thus processed, it is time to actually update the screen. That is done simply by updating the pygame.surfarray.pixels2d(self.screen) by setting each pixel defined by coordinates self.scr_x, self.scr_y to color self.scr_col as stored previously. It would have been possible to do this in smaller pieces when collecting the data, but it seems to be faster to store all data first and then update in one go, although it does consume more memory.

        # update screen in one go by setting the colors for the mapped pixels.
        if self.scr_cnt > 0:
            screen_array = pygame.surfarray.pixels2d(self.screen)
            screen_array[self.scr_x[0:self.scr_cnt], self.scr_y[0:self.scr_cnt]] = self.scr_col[0:self.scr_cnt]

And in case you are wondering how the animated image with the green balls and red and blue cogwheels is done, it is simply the case of having ten still images and switching between those images when getting the colors for the pixels.

Singularity

The process is very CPU intensive and more than 90 % of time is spent processing the segments line by line. It would be great to use more than one processor core to share that load, but when I tried that, the end results was slower execution. Perhaps I did not do it the right way or perhaps there are some locks I missed; or maybe the overhead, when running real time at 30 to 50 fps, was too much. 

Of course, using hardware to do the same would be way faster. Also much faster, even if not that much, would be using some optimized routines instead of high level Python code. Using e.g. OpenGL would boost performance probably a lot, although it seems rather complex when looking at PyOpenGL. Perhaps some day - although using that would mean missing all the fun of making it yourself!

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.