Real-Time 3D with Python Part III: Visibility and Shading

In the introduction I explained the idea: Using Python to program old skool demo style real-time vector graphics. In Part I I set up a simple 3D wireframe object, rotated it, and drew it on the screen, and in Part II added surfaces and perspective. This part finds out whether surfaces are visible or not, thus saving time by not drawing the invisible ones, and adds shading using a light source.

The program code for all parts can be downloaded from github. Note that the full code is not found below - to run this, use github code.

Having an Angle

To find out if a surface is pointing towards the viewer (and should be drawn) or not, we need to figure out the angle between the vector from the viewer to the plane and the vector orthogonal to the plane. In the picture below that would be vectors r and n:

Vector r from viewer and its direction vector u, crossing the plane, and vector n, the vector orthogonal to the plane. Source: superprof resources.

Now, how do we get vector n? It helps a bit that the vector is, obviously, the same at all locations on the plane, because, well, it's a plane and not some curvy surface. We can use one of the corners of the plane and the two vectors from that corner to the closest two other corners to calculate n by taking a cross product of these two vectors:

    def updateSurfaceCrossProductVector(self):
        # calculate cross product vector for each surface using rotatedNodes
        # always use vectors (1, 0) and (1, 2) (numbers representing nodes)
        # numpy "cross" was terribly slow, calculating directly as below takes about 10 % of the time.
        for surface in self.surfaces:
            vec_A = self.rotatedNodes[surface.nodes[2], 0:3] - self.rotatedNodes[surface.nodes[1], 0:3]
            vec_B = self.rotatedNodes[surface.nodes[0], 0:3] - self.rotatedNodes[surface.nodes[1], 0:3]
            vec_Cross =  ([
                vec_B[1] * vec_A[2] - vec_B[2] * vec_A[1],
                vec_B[2] * vec_A[0] - vec_B[0] * vec_A[2],
                vec_B[0] * vec_A[1] - vec_B[1] * vec_A[0]
                ])
            surface.setCrossProductVector(vec_Cross)

Numpy does have a cross function, but it uses way too much time figuring out the shapes of the vectors etc., while we already know they will always be ordinary 3D vectors. Hence, it is much faster to calculate the cross product vector as above. This is done for all surfaces of the object.

Now, you may notice that there are actually two (unit length) orthogonal vectors for the surface; one pointing "up" and one pointing "down". How do we know which one to calculate? In the previous part I explained that each surface must be defined by adding the nodes clockwise (from viewing direction) - this ensures that the vec_A and vec_B above always produce a vector pointing "out" of the object.

Next, we need to calculate the angle. Actually, we do not need the angle itself; we just need to know whether the angle is such that the surface is pointing towards the viewer. This is the math we need:

Determining the angle between two vectors. Source: superprof resources.

And this is how it is done in the code:

    def updateSurfaceAngleToViewer(self):
        # calculate acute angle between surface plane and Viewer
        # surface plane cross product vector and Viewer vector both from node 1.
        for surface in (surf for surf in self.surfaces if surf.visible == 1):
            vec_Viewer = self.rotatedNodes[surface.nodes[1], 0:3] 
            surface.setAngleToViewer(vec_Viewer)
            if surface.angleToViewer > 0 or surface.showBack == 1:
                surface.setVisible(1)
            else:
                surface.setVisible(0)

    def setAngleToViewer(self, vec_Viewer):
        if self.crossProductLen > 0 and vec_Viewer.any() != 0:
            # instead of true angle calculation using asin and vector lengths, a simple np.vdot is sufficient to find the sign (which defines if surface is visible) 
            # self.angleToViewer = math.asin(np.dot(self.crossProductVector, vec_Viewer) / (self.crossProductLen * np.linalg.norm(vec_Viewer)))
            self.angleToViewer = np.dot(self.crossProductVector, vec_Viewer) 

So, for each surface, I am defining the vector from the viewer (who sits at [0,0,0]) to the same corner of the surface I defined the orthogonal cross product vector for, and then taking a dot product. That's basically the dividend part of the lowest formula above. Since we know that the divisor, which is the product of the lengths of the vectors, is always positive, we know from the sign of the dividend whether we are taking an arcsin of a positive or a negative number. And, we also know that the number will be between -1 and 1. All this means the angle will have the same sign as the dividend. If the angle is positive, the surface is pointing towards the viewer, and it will have its Visible property set to one. We don't really know the actual size of the angle, but it does not matter for this purpose.

Note the property showBack above. If a surface is defined to "show its back" then it will always be visible, no matter what the sigh of the angle is.

Shading

We can use the same technique to shade a surface using a light source. Instead of calculating the angle between the vector from the viewer and the surface, we need to calculate the angle between a vector from some light source and the surface. The process is identical, but in this case, we are interested in the actual angle. With the angle, we can define the shading so that if the angle shows the surface is pointing directly towards the light source, we set the "shading factor" to one; if it is pointing away, we set it to something small (minShade); and in between, in proportion to the angle.

    def updateSurfaceAngleToLightSource(self, lightPosition):
        # calculate acute angle between surface plane and light source, similar to above for Viewer.
        # this is used to define shading and shadows; needed for visible surfaces using shading AND all surfaces, if shadow to be drawn.
        for surface in (surf for surf in self.surfaces if surf.visible == 1 and self.minShade < 1.0):
            surface.setLightSourceVector(self.rotatedNodes[surface.nodes[1], 0:3] - lightPosition)
            surface.setAngleToLightSource()

    def updateSurfaceColorShade(self):
        # calculate shade for surface. 
        for surface in (surf for surf in self.surfaces if surf.visible == 1):
            surface.setColorShade(self.minShade)
            if surface.showBack == 1: surface.setBackColorShade(self.minShade)

    def setAngleToLightSource(self):
        if self.crossProductLen > 0 and self.lightSourceVector.any() != 0:
            # instead of true angle calculation using asin and vector lengths, a simple np.vdot is sufficient to find the sign (which defines if surface is shadowed) 
            # self.angleToLightSource = math.asin(np.dot(self.crossProductVector, self.lightSourceVector) / (self.crossProductLen * np.linalg.norm(self.lightSourceVector))) / (np.pi / 2)
            self.angleToLightSource = np.dot(self.crossProductVector, self.lightSourceVector)
        
    def setColorShade(self, minShade):
        if self.angleToLightSource <= 0:
            self.colorShade = minShade
        else:
            self.colorShade = minShade + (1.0 - minShade) * math.asin(self.angleToLightSource / (self.crossProductLen * self.vectorLen(self.lightSourceVector))) / (np.pi / 2)

    def colorRGB(self):
        use_color = ([round(self.colorShade * x, 0) for x in self.color]) # apply shading
        return use_color

    def vectorLen(self, vector):
        # equivalent to numpy.linalg.norm for a 3D real vector, but much faster. math.sqrt is faster than numpy.sqrt.
        return math.sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2])

Now, I do the same as above with updateSurfaceAngleToLightSource and setAngleToLightSource; then when setting the actual colorShade (a factor always between 0 and 1) I use the rest of the lowest formula above: I divide the cross product vector with the lengths of the vectors, and take an arcsin. Note that I use vectorLen to get the length of the light source vector (again numpy alternative was very slow), but the length of the cross product vector is stored in crossProductLen, as it is constant (as long as the object is constant). These calculations are rather costly, so I try to skip as much of them as possible; first of all if the surface is not visible at all, and secondly, if minShade is one (no shading), and thirdly, if the angle is negative (and minShade is used in any case).

Lastly, in colorRGB, I adjust the surface color by multiplying it with colorShade.

Now I have the same cube as in the previous part, but with shading. What can not be seen is that now the surfaces in the back side not visible to the viewer are not drawn at all.

Need for Speed 2: Don't Do It

In the previous part I explained how precalculating some of the CPU intensive operations was used on the Amiga to make code run faster. Of course, as above, often the fastest code is the code which is not run at all. With Python, I am using some of these methods, but also taking care not to sacrifice how the program can be used. For example, in the original Amiga demo I dropped the ability to actually rotate the objects around all three axes and settled for the Y (vertical) axis, which has several advantages:

• calculating the rotation matrix becomes much simpler, when X and Z angles are zero - less calculations needed for rotations.
• calculating the shadows becomes much easier when the ground is always at Y coordinate zero.
• the ground needs no drawing at all, but it can be shown by setting the background color using the Copper co processor (almost no CPU time needed).

In the Python project I kept the ability to rotate around all axes, with quite a bit of added complexity and computing cost for the above. 

Next: Part IV: Cityscape.

Real-Time 3D with Python Part II: Surfaces and Perspective

In the introduction I explained the idea: Using Python to program old skool demo style real-time vector graphics. In Part I I set up a simple 3D wireframe object, rotated it, and drew it on the screen. This part is about adding surfaces to that wireframe, and how to use perspective so that more distant objects appear smaller on screen, and objects in general look more natural.

The program code for all parts can be downloaded from github. Note that the full code is not found below - to run this, use github code.

Defining surfaces

A cube object has eight nodes and six surfaces. For the program to reasonably handle a surface, it has to be flat - all its nodes must lay on a single plane. Then it is possible to define if the surface is visible or not and to draw it using a single polygon. To tell the program what those surfaces are, I define a new class:

class VectorObjectSurface:

    """
    Surfaces for a VectorObject.
    
    @author: kalle
    """
    def __init__(self):
        self.idnum = 0
        # properties set when defining the object
        self.nodes = []
        self.color = (0,0,0)
        self.edgeWidth = 0
        # the following are calculated during program execution
        self.zpos = 0.0

    def setZPos(self, zpos):
        self.zpos = zpos

Then, for the VectorObject class I add a list to hold these, and code to set one surface up and add it to the list:

    def addSurfaces(self, idnum, color, edgeWidth, node_list):
        # add a Surface, defining its properties
        surface = VectorObjectSurface()
        surface.idnum = idnum
        surface.color = color
        surface.edgeWidth = edgeWidth
        surface.nodes = node_list
        self.surfaces.append(surface)

The surfaces can then be set up after the nodes have been defined (see previous part). The colour of the surface is here defined as white for the first two, light grey for the next two, and dark grey for the last two:

    node_list = [0, 3, 2, 1] # node_list defines the four nodes forming a cube surface, in clockwise order
    vobj.addSurfaces(0, (255,255,255), 0, node_list)
    node_list = [4, 5, 6, 7]
    vobj.addSurfaces(1, (255,255,255), 0, node_list)
    node_list = [0, 1, 5, 4]
    vobj.addSurfaces(2, (150,150,150), 0, node_list)
    node_list = [3, 7, 6, 2]
    vobj.addSurfaces(3, (150,150,150), 0, node_list)
    node_list = [0, 4, 7, 3]
    vobj.addSurfaces(4, (80,80,80), 0, node_list)
    node_list = [1, 2, 6, 5]
    vobj.addSurfaces(5, (80,80,80), 0, node_list)

Note that the order of the nodes is always clockwise (looking at the object from outside). Why? We will need to figure out if the surface is visible to the viewer, and unless (in later parts) specifically defined to show also its back side, only the front side will be drawn. When defining the nodes in clockwise order, we also which side is front and which one is back.

Obviously, the nodes have to be in correct order going around the edges of the surface, not crossing it. The last node will automatically be connected back to the first one, so for example node_list = [0, 3, 2, 1] defines a polygon with edges [0, 3], [3, 2], [2, 1], and [1, 0].

A couple of other things we need to add to the VectorObject class.

    def updateSurfaceZPos(self):
        # calculate average Z position for each surface using rotatedNodes
        for surface in self.surfaces:
            zpos = sum([self.rotatedNodes[node, 2] for node in surface.nodes]) / len(surface.nodes) 
            surface.setZPos(zpos)
        
    def sortSurfacesByZPos(self):
        # sorts surfaces by Z position so that the most distant comes first in list
        self.surfaces.sort(key=lambda VectorObjectSurface: VectorObjectSurface.zpos, reverse=True)

When the object has been rotated, I want to draw the most distant surfaces first and the ones closer after them, so that if they overlap the closer surfaces will be shown over the more distant ones and not the other way around. (For a simple cube this is not a problem, but the object could be more complex.)

First I update the average distance from viewer (Z coordinate) for every surface in updateSurfaceZPos. It simply sums up the rotated Z coordinate of each of the surface's nodes and the divides the sum by the number of nodes. Then, having that stored for all surfaces, I can sort the list of surfaces by it.

Then, before drawing the object in display (see previous part), I can use these to make sure I am drawing the surface in correct order.

        # draw the actual objects
        for VectorObj in self.VectorObjs:
            # first sort object surfaces so that the most distant is first.
            VectorObj.updateSurfaceZPos()
            VectorObj.sortSurfacesByZPos()
            # then draw surface by surface. 
            for surface in VectorObj.surfaces:
                # build a list of transNodes for this surface
                node_list = ([VectorObj.transNodes[node][:2] for node in surface.nodes])
                pygame.draw.aalines(self.screen, surface.color, True, node_list)       
                pygame.draw.polygon(self.screen, surface.color, node_list, surface.edgeWidth)       

Note that unlike in the previous part, I am not drawing any circles to mark the nodes, and the aalines here is actually not required, it simply makes the filled polygon have anti-aliased edges. All surfaces are drawn, even if we know at least three of them (this being a cube) will be pointing away from the viewer on the back side, but at this stage we are still missing the code to determine which ones. (A simple solution would be to not draw the three most distant ones - but that only works for a cube-like object.)

Have some perspective, please

How about adding that perspective to the transformation from 3D to 2D? It is really simple. To make more distant objects appear smaller, I simply use the distance as a divisor. For scaling I have defined (in VectorViewer) self.zScale = width * 0.7. This can be tweaked to easily show objects bigger or smaller. Then, redefining the transform method:

    def transform(self, zScale, midScreen):
        """ 
         Flatten from 3D to 2D and add screen center.
        """
        # apply perspective using Z coordinates and add midScreen to center on screen to get to transNodes.
        # for normal objects, some of the transNodes will not be required, but possibly figuring out which are and processing them
        #   individually could take more time than this.
        self.transNodes = (self.rotatedNodes[:, 0:2] * zScale) / (self.rotatedNodes[:, 2:3]) + midScreen

As you can see (comparing to Part I), the only change here is to use the zScale as a scaling multiplier and the Z coordinate of rotatedNodes as a divisor. Hey, Presto! We have perspective.

This is where we are at this stage:

Need for Speed

On the Amiga (or more generally the MC68000 chip), division was a terribly expensive operation, taking about 13 milliseconds each. Avoiding division was thus very important, but with e.g. scaling above it is also very necessary. A common trick used here was to replace a DIVS instruction with a MULS and a SWAP,  using a precalculated table - in quite a similar manner as for the sine and cosine operators described in the previous part. In effect, in the beginning of the program, we can store in memory a table of $10000/x for all x we need as divisor (remember, integers only), and then replace all divide operations with 1) fetching the correct precalculated number from the table, 2) multiplying the dividend with it, and 3) swapping the halves of the result (taking the higher 16 bits to the lower 16 bits, in effect dividing it with $10000).

Despite the added complexity, this could save perhaps 40+ % of the CPU time needed - a MULS takes less than half the cycles a DIVS requires, while accessing the table and swapping the are relatively quick. Using tricks like this was based on detailed diagrams of instruction timings, and of course the full control of the CPU that Assembly language provides.

Next: Part III: Visibility and Shading.

Real-Time 3D with Python Part I: Rotating an Object

In the introduction I explained the idea: Using Python to program old skool demo style real-time vector graphics. This part is about setting up a simple object, rotating it, and drawing it on the screen.

The program code for all parts can be downloaded from github. Note that the full code is not found below - to run this, use github code.

Setting up 

I am going to use numpy for matrices and pygame for setting up the screen and drawing onto it, so the first thing to do is

import pygame
import numpy as np

To set up a simple object, I have defined a cube in XYZ coordinates that has the size (width, height, depth) of 200. It will be centered at (0,0,0). The nodes (the eight corners of the cube, in 3D space) are thus defined as follows:

    node_array = np.array([
            [ 100.0, 100.0, 100.0],
            [ 100.0, 100.0,-100.0],
            [ 100.0,-100.0,-100.0],
            [ 100.0,-100.0, 100.0],
            [-100.0, 100.0, 100.0],
            [-100.0, 100.0,-100.0],
            [-100.0,-100.0,-100.0],
            [-100.0,-100.0, 100.0]
            ])

To be clear, I have defined the axis so that X is the horizontal axis, Y is the vertical axis, and Z is the distance frow viewer axis. The standard Cartesian coordinates are usually presented differently:

Cartesian coordinate system. Source: Wikipedia.

Basically, I have rotated the coordinate system above so that X=y, Y=z, and Z=-x. There are other coordinates systems, like polar coordinates, but I will stick to this one.

I want to rotate objects around any of the three axes. To do that, I will define three angles, X, Y, and Z (no surprises here). It is worth noting that, in this system, the order of rotation matters; I am using the order as above. Now it is time to define a VectorObject class.

class VectorObject:

    """
    Position is the object's coordinates.
    Nodes are the predefined, static definition of object "corner points", around object position anchor point (0,0,0).
    RotatedNodes are the Nodes rotated by the given Angles and moved to Position.
    TransNodes are the RotatedNodes transformed from 3D to 2D (X.Y) screen coordinates.
    
    @author: kalle
    """
    def __init__(self):
        self.position = np.array([0.0, 0.0, 0.0, 1.0])      # position
        self.angles = np.array([0.0, 0.0, 0.0])
        self.angleScale = (2.0 * np.pi) / 360.0             # to scale degrees. 
        self.rotationMatrix = np.zeros((3,3))
        self.rotateSpeed = np.array([0.0, 0.0, 0.0])
        self.nodes = np.zeros((0, 4))                       # nodes will have unrotated X,Y,Z coordinates plus a column of ones for position handling
        self.rotatedNodes = np.zeros((0, 3))                # rotatedNodes will have X,Y,Z coordinates after rotation ("final 3D coordinates")
        self.transNodes = np.zeros((0, 2))                  # transNodes will have X,Y coordinates

    def setPosition(self, position):
        # move object by giving it a rotated position.
        self.position = position 

    def setRotateSpeed(self, angles):
        # set object rotation speed.
        self.rotateSpeed = angles 

    def addNodes(self, node_array):
        # add nodes (all at once); add a column of ones for using position in transform
        self.nodes = np.hstack((node_array, np.ones((len(node_array), 1))))
        self.rotatedNodes = node_array # initialize rotatedNodes with nodes (no added ones required)

    def increaseAngles(self):
        self.angles += self.rotateSpeed
        for i in range(3):
            if self.angles[i] >= 360: self.angles[i] -= 360
            if self.angles[i] < 0: self.angles[i] += 360

Note that there are three sets of object nodes in the class:

1. nodes: The original nodes for the object. The fourth element (in addition to X, Y, Z coordinates) will be set to 1.0 to add the position when rotating.
2. rotatedNodes: nodes as rotated to the current angles. 3D i.e. X, Y, Z.
3. transNodes: rotatedNodes transformed (or flattened) from 3D to 2D screen coordinates. 2D i.e. X, Y only.

Also, we have angles, and the respective rotationMatrix. The angles are represented in degrees so that a full circle is 360 degrees, so we need a constant angleScale to convert them to radians.

Setting up a rotation matrix involves taking a sine and cosine of all the angles. Then the successive rotations are combined in the matrix (see Wikipedia for Euler angles). Having the matrix, it is simple to convert the static nodes to their rotated form: simply perform a dot product of the nodes and the matrix (in rotate).

The matrix and rotation are defined as follows:

    def setRotationMatrix(self):
        """ Set matrix for rotation using angles. """
        
        (sx, sy, sz) = np.sin((self.angles) * self.angleScale)
        (cx, cy, cz) = np.cos((self.angles) * self.angleScale)
 
        # build a matrix for X, Y, Z rotation (in that order, see Wikipedia: Euler angles) including position shift. 
        # add a column of zeros for later position use
        self.rotationMatrix = 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 ]])

    def rotate(self):
        """ 
        Apply a rotation defined by a given rotation matrix.
        """
        matrix = np.vstack((self.rotationMatrix, self.position[0:3]))   # add position to rotation matrix to move object at the same time
        self.rotatedNodes = np.dot(self.nodes, matrix)

    def transform(self, midScreen):
        """ 
        Add screen center.
        """
        # add midScreen to center on screen to get to transNodes.
        self.transNodes = self.rotatedNodes[:, 0:2] + midScreen

To transform these rotated 3D nodes to a screen-friendly 2D plane, I (in this first phase) contend to just forgetting about the distance axis Z and simply taking the X and Y coordinates and centering the object on screen. In a way this is similar to looking at the object from very far away with a telescope, when the differences in distance from the viewer between the nodes are so small they are meaningless.

Making Something Happen

Now we need to do something with the class defined above. First, I define a couple of keys so that we can easily exit and pause the program:

key_to_function = {
    pygame.K_ESCAPE: (lambda x: x.terminate()),         # ESC key to quit
    pygame.K_SPACE:  (lambda x: x.pause())              # SPACE to pause
    }

Then let's set up the main class.

class VectorViewer:
    """
    Displays 3D vector objects on a Pygame screen.
    
    @author: kalle
    """

    def __init__(self, width, height):
        self.width = width
        self.height = height
        self.screen = pygame.display.set_mode((width,height))
        self.fullScreen = False
        pygame.display.set_caption('VectorViewer')
        self.backgroundColor = (0,0,0)
        self.VectorObjs = []
        self.midScreen = np.array([width / 2, height / 2], dtype=float)
        self.target_fps = 60                            # affects movement speeds
        self.running = True
        self.paused = False
        self.clock = pygame.time.Clock()
            
    def addVectorObj(self, VectorObj):
        self.VectorObjs.append(VectorObj)

This class basically does two things: It rotates and displays the objects. The main loop run just checks if any of the defined keys were pressed, and runs first rotate and then display:

    def run(self):
        """ Main loop. """
                  
        while self.running:
            for event in pygame.event.get():
                if event.type == pygame.QUIT:
                    self.running = False
                elif event.type == pygame.KEYDOWN:
                    if event.key in key_to_function:
                        key_to_function[event.key](self)
            
            if self.paused == True:
                pygame.time.wait(100)
            
            else:
                # main components executed here
                self.rotate()
                self.display()
                
                # release any locks on screen
                while self.screen.get_locked():
                    self.screen.unlock()
                    
                # switch between currently showed and the next screen (prepared in "buffer")
                pygame.display.flip()
                self.clock.tick(self.target_fps) # this keeps code running at max target_fps

        # exit; close display, stop music
        pygame.display.quit()

Rotation actually happens in the VectorObject class as described earlier. The display method then does the following:

1.  clears the screen by filling it with the background colour.
2.  for each node, draws a small circle, in white (RGB 255,255,255).
3.  draws a line between specified nodes to show the cube's edges.

    def rotate(self):
        """ 
        Rotate all objects. First calculate rotation matrix.
        Then apply the relevant rotation matrix with object position to each VectorObject.
        """
                
        # rotate and flatten (transform) objects
        for VectorObj in self.VectorObjs:
            VectorObj.increaseAngles()
            VectorObj.setRotationMatrix()
            VectorObj.rotate()
            VectorObj.transform(self.midScreen)
            
        
    def display(self):
        """ 
        Draw the VectorObjs on the screen. 
        """

        # lock screen for pixel operations
        self.screen.lock()

        # clear screen.
        self.screen.fill(self.backgroundColor)
                   
        # draw the actual objects
        for VectorObj in self.VectorObjs:       
            for node in VectorObj.transNodes:
                # for each node, draw a small circle (radius = 5) in white (255,255,255)
                node_int = (int(node[0]), int(node[1]))
                pygame.draw.circle(self.screen, (255,255,255), node_int, 5)
            # just for testing purposes - assuming a cube with specific node order, draw edges
            if np.shape(VectorObj.transNodes)[0] == 8:
                pygame.draw.aalines(self.screen, (255,255,255), 1, VectorObj.transNodes[0:4,:])
                pygame.draw.aalines(self.screen, (255,255,255), 1, VectorObj.transNodes[4:8,:])
                for i in range(4):
                    pygame.draw.aaline(self.screen, (255,255,255), VectorObj.transNodes[i,:], VectorObj.transNodes[i+4,:])
        
        # unlock screen
        self.screen.unlock()
 
    def terminate(self):

        self.running = False   

    def pause(self):

        if self.paused == True:
            self.paused = False
        else:
            self.paused = True

Note that for all drawing operations, transNodes are used, as they are the 2D screen coordinates. Also note that step 3, drawing the edges, is here just for testing purposes - our cube object only defines the nodes and if they would be in a different order, the edges as defined here might look very strange! This will be fixed in the following part.

Finally, this is the code to set up the cube object and initialize the VectorViewer and then run it.

if __name__ == '__main__':
    """ 
    Prepare screen, objects etc.
    """

    # set screen size
    # first check available full screen modes
    pygame.display.init()
    # disp_modes = pygame.display.list_modes(0, pygame.FULLSCREEN | pygame.DOUBLEBUF | pygame.HWSURFACE)
    # disp_size = disp_modes[4] # selecting display size from available list. Assuming the 5th element is nice...
    disp_size = (1280, 800)

    vv = VectorViewer(disp_size[0], disp_size[1])

    # set up a simple cube
    vobj = VectorObject()
    node_array = np.array([
            [ 100.0, 100.0, 100.0],
            [ 100.0, 100.0,-100.0],
            [ 100.0,-100.0,-100.0],
            [ 100.0,-100.0, 100.0],
            [-100.0, 100.0, 100.0],
            [-100.0, 100.0,-100.0],
            [-100.0,-100.0,-100.0],
            [-100.0,-100.0, 100.0]
            ])
    vobj.addNodes(node_array)
    speed_angles = np.array([1.0, -.3, 0.55])
    vobj.setRotateSpeed(speed_angles)
    position = np.array([0.0, 0.0, 1500.0, 1.0])
    vobj.setPosition(position)
 
    # add the object
    vv.addVectorObj(vobj)
     
    # run the main program
    vv.run()

This is what you should see:

When Points Do Not Float

How to do this in the Amiga? Using Assembly, we have access directly to the CPU's instruction set. The eight internal data registers are 32-bit, i.e. they can handle numbers up to 2^32-1 (4,294,967,295), but there is no such concept as decimals. Every number is and has to be an integer. There is no FPU (Floating Point Unit) either. In Python, if we need a sine or cosine of any number, we just call numpy.sin or math.sin, and there we have it. numpy.sin(numpy.pi / 4) = 0.7071067811865476.  How is it calculated? I don't know. 

In the MC68000 instruction set, there are instructions for adding, subtracting, multiplying, and dividing two integers together (bear in mind that e.g. the Commodore 64 CPU did not have a multiply nor divide instruction, so even this was great). There's certainly no instruction to get a sine. And the sine of any number is, by definition, between -1 and +1. If you only have integers, having a precision of -1, 0, or +1 for your sine and cosine certainly is is not enough for a smooth rotation. 

Solving Sine in the Amiga

First of all, a decimal number is nothing more than a big integer divided by another (e.g. 0.2 = 20,000 / 100,000), so scaling is the solution for the missing decimal point. A convenient scaling factor in the MC68000 is $10000 (hexadecimal for 65536 or 2^16), which can be de-scaled by a very quick instruction called swap. Swap simply swaps the two 16-bit parts of a 32-bit register, in effect dividing the 32-bit result with $10000 (when using only the lower 16 bits and discarding the higher 16 bits).

To get a sine for x, there are a number of methods. I used a Taylor series as follows (see link):

   

And as accuracy is not that critical a factor here, the series can basically be terminated at the fourth, or even third, term. However, care must be taken as the 32-bit data register soon overflows if you raise almost anything to the 7th power. Anyway, this would be way too slow for real-time calculation - multiplications are quite slow. A pre-calculated table is needed. Let's say we divide the full circle to 360 degrees and calculate a sine for all of them, and store the results in memory. Then, when a sine for angle 124 degrees is needed, we just pick the number which is 124 * 2 bytes (as the numbers take 2 bytes each) away from the start of our sine table! Easy and, even more importantly, very fast!

Of course, we need the cosine as well. Luckily - well, actually, it has nothing to do with luck, but still: cos(x) equals sin(x + 90 degrees), so we can use the same table just adjusting the offset.

Usually, it was more convenient to precalculate the sine table and add it to the program code as data than to actually include the precalculation routine to the code.


Next: Part II: Surfaces and Perspective.

Acknowledgement: I took inspiration from Peter Collingridge's blog.

Real-Time 3D Vector Graphics with Python

I used to program real-time graphics demos on the Commodore Amiga at the end of the 1980's and the beginning of the 1990's. The programming language of choice at the time was assembly, a.k.a. machine language, to be able to squeeze everything possible (and then some) out of the hardware.

Now, some 30 years later, I wanted to learn Python, perhaps to read data from files and perform some math operations on it. In need of a fun project to help me get a decent start, I decided to try and replicate what I did back then with assembly. Using a high-level, run-time interpreted language like Python adds a lot of overhead compared to executing assembly, and so does Windows in the background. On the other hand, my ~4 year old laptop has a 2.6Ghz 64-bit Intel i7 with four cores, 16GB RAM and a graphics card, while the Amiga 1000 had a 7.09Mhz 16/32-bit Motorola 68000, 512kB RAM and some co-processors. The i7 has more L2 cache memory than the Amiga (no cache whatsoever) has total RAM. According to the famous Moore's law (loosely interpreted) I might have to the tune of 30 000 times more computing power!

Sample of  (my) MC68000 assembly source code (1991).

So, this is where I took my inspiration: The 3D vector world I programmed for the demo Sound Vision, released by our group Reflect at the very first Assembly Demo Party in Kauniainen, Finland, in July 1992. (We finished 1st in the demo competition. Assembly is still very active, but nowadays more about gaming.) Some vector world routines had of course been presented earlier, but this was the first one with a moving light source and real time calculated shadows.

This blog is partly about how to do this with Python, and partly about how this was done 30 years ago ("old school" a.k.a. "oldskool") with assembly on the Amiga.

The whole demo, as a video, can be found at e.g. Assembly Archive. The part above I actually captured running an Amiga emulator on my laptop (WinUAE) with an original disk image of the demo, which is less than 880 kB including all the code, music, and graphics!  The Amiga (original chip set) was capable of 352 x 288 pixel resolution in 32 colors and basically there was so much happening in this part (considering the resources) that the frame rate (pictures per second) is at times as low as 12-15 - hence, it does look somewhat jerky.  

So, how to go about it with Python? First, I installed Anaconda on my laptop, and started Spyder. Quickly I found out that I needed a real-time graphics package, so I installed Pygame. Pygame can be used to set up a screen or a window, and to draw on it - something that was handled by the Copper and the Blitter co-processors in the Amiga, respectively. Also, a number of other - more standard - packages like math and numpy came in handy.

Setting up an Amiga demo

The Amiga had 512kB of "chip" RAM, ie. memory available for program code and usable to the graphics and sound co-processors. The separate ROM (read-only memory) contained the BIOS; in my Amiga 1000 even this had to be loaded from a "kickstart" disk at startup. A track loaded demo ("trackmo") would boot from a disk and then run the first 1kB stored on that disk (the "bootblock"). Basically, you could have your own code there then using the disk drive hardware directly to proceed - asking it to load "tracks" from the disk directly to memory. The Amiga OS was not needed for anything. No libraries, drivers etc., but accessing the hardware registers directly.

So once you had loaded some music, perhaps graphics like a font for a "scroller", and your actual program code, then what? The first thing to do with the demo code was to take hold of system interruptions (these "interrupt" the currently running code and give the CPU temporarily for some other code to use) and create a "copper list". The copper list would tell the Copper co-processor where the bitplanes to be shown on screen are located, how the colours (up to 32) are defined, and some other graphics related settings. Then, kick off with the music and start drawing on the screen! For the drawing operations the Blitter would be used. According to the manual, the Blitter is "the high speed line drawing and block movement component of the system" and able to clear, copy and fill areas, and to draw lines. So, to create a simple text scroller, first find the location of the first letter of your scroller text, and copy the picture of that letter to the rightmost edge of your screen. Then, in the next frame, move the screen image some bits to the left. Once the whole letter has been moved enough times, there's room for the next letter, etc. See this tutorial for an example!

More on demo programming:
    Crash Course to Amiga Assembly Programming
    Coder Aid
    Amiga Hardware Programming (video)

Next: Part I: Rotating a 3D object.