Oct 17, 2021

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 [0123] 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 [4567] 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, 30, 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. [4765], 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[01230123].

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

Oct 3, 2021

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.