Real-Time 3D with Python Part IV: Cityscape

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. Part III was about calculating object surface visibility and shading. This part sets up an XML reader to define the objects and other properties using a common file format (XML), and adds the code needed to work with multiple objects.

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

It's a World Out There

Previously, I only had a single object shown in the middle of the screen. Now I am going to add many objects, so that they all have their own coordinates. In addition, they can have their own rotation defined, so I may need multiple "angle sets". What I will end up with all the changes is this:

To get there, let's add two new classes, VectorAngles and VectorPosition:

class VectorAngles:

    """
    Angles for rotating vector objects. For efficiency, one set of angles can be used for many objects.
    Angles are defined for axes X (horizontal), Y (vertical), Z ("distance") in degrees (360).
    
    @author: kalle
    """
    def __init__(self):
        self.angles = np.array([0.0, 0.0, 0.0])
        self.angleScale = (2.0 * np.pi) / 360.0 # to scale degrees.
        self.angName = ""
        self.rotationMatrix = np.zeros((3,3))
        self.rotateAngles = np.array([0.0, 0.0, 0.0])
        self.rotate = np.array([0.0, 1.0, 0.0])
    
    def setAngles(self, angles):
        # Set rotation angles to fixed values. 
        self.angles = angles
     
    def setRotateAngles(self):
        self.rotateAngles += self.rotate
        for i in range(3):
            if self.rotateAngles[i] >= 360: self.rotateAngles[i] -= 360
            if self.rotateAngles[i] < 0: self.rotateAngles[i] += 360

    def setRotationMatrix(self):
        # Set matrix for rotation using angles.
        
        (sx, sy, sz) = np.sin((self.angles + self.rotateAngles) * self.angleScale)
        (cx, cy, cz) = np.cos((self.angles + self.rotateAngles) * 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 ]])
    
class VectorPosition:

    """
    A vector object defining the positions of other objects in its nodes (see VectorObject).
    
    @author: kalle
    """
    def __init__(self):
        self.position = np.array([0.0, 0.0, 0.0, 1.0])
        self.angles = VectorAngles()
        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
        self.objects = []                               # connects each node to a respective VectorObject
        self.objName = ""
        
    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 with nodes
        
    def addObjects(self, object_list):
        self.objects = object_list

    def rotate(self):
        # apply a rotation defined by a given rotation matrix.
        matrix = np.vstack((self.angles.rotationMatrix, np.zeros((1, 3))))
        # apply rotation and position matrix to nodes
        self.rotatedNodes = np.dot(self.nodes, matrix) + self.position[0:3]

As you may note, I have just moved the code in VectorAngles from VectorObject, where it still was in Part III. Now that we have multiple objects, most of them share their rotation (angles), so it is better to process that part just once per rotation set.

VectorPosition is new but very simple. The nodes hold the original coordinates of each object, and rotatedNodes hold the coordinates after rotation according to angles. In addition, there's a list objects to maintain the connection between the coordinates and the respective object. In a way, VectorPosition is very similar to a VectorObject, but it does not have surfaces. Instead, each of its nodes just represents the position of a VectorObject.

Now, when rotating the objects in VectorViewer, rotation is done in two steps: first the positions above are rotated and copied into the actual objects, which are then rotated in the second step:

    def rotate(self):
        """ 
        Rotate all objects. First calculate rotation matrix.
        Then apply the relevant rotation matrix with object position to each VectorObject.
        """
                
        # calculate rotation matrices for all angle sets
        for VectorAngles in self.VectorAnglesList:
            VectorAngles.setRotateAngles()
            VectorAngles.setRotationMatrix()
        
        # rotate object positions, copy those to objects.
        self.VectorPos.rotate() 
        for (node_num, VectorObj) in self.VectorPos.objects:
            VectorObj.setPosition(self.VectorPos.rotatedNodes[node_num, :])

        # rotate and flatten (transform) objects
        for VectorObj in self.VectorObjs:
            VectorObj.updateVisiblePos(self.objMinZ) # test for object position Z
            if VectorObj.visible == 1:
                VectorObj.rotate() # rotates objects in 3D
                VectorObj.updateVisibleNodes(self.objMinZ) # test for object minimum Z
                if VectorObj.visible == 1:
                    VectorObj.transform(self.midScreen, self.zScale) # flattens to 2D, crops X,Y
                    VectorObj.updateVisibleTrans(self.midScreen) # test for outside of screen

Note how first all possible "angle sets" are processed so that the rotation matrices are ready for use. Also, when rotating the objects, I am updating a visible property; now that we have many objects moving around, it may be possible to stop processing some of them, if we know they will not be visible anyway. First, I am testing for object position Z coordinate. In my world, the viewer is sitting at the origo, so anything with a negative Z coordinate is behind her. I have actually defined an objMinZ property to not draw anything which is too close. So, if the object is indeed too close, judging by its position, it will not be shown and thus also rotating it will be unnecessary.

A second test, for the objects surviving the first, is to check is the object is on screen at all. This is done by checking the minimum and maximum transNodes i.e. X and Y coordinates against the screen boundaries. (For the actual code for these simple checks, see the github files.) 

City Planning

To set up the objects, I created an XML file. In the code, I imported some new modules to do that:

import os
from operator import itemgetter
import copy
import xml.etree.ElementTree as et

And then, before launching the VectorViewer, the part below reads the XML and adds the objects. (It's a lengthy part so this is just the beginning.) If you run the github code, you will need to change the data directory to point to the right place. Handling XML is very easy; the ElementTree.parse builds a structure which can then be iterated or searched for the data as required. Note that the code I built contains almost no error checking whatsoever for missing or incorrect data.

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

    # set data directory
    os.chdir("D:\kalle\Documents\Python")

    # 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])

    # read data file defining angles, movements and objects
    vecdata = et.parse("vectordata cityscape.xml")

    root = vecdata.getroot()
    
    for angles in root.iter('vectorangles'):
        ang = VectorAngles()
        ang.angName = angles.get('name')
        ang.angles[0] = float(angles.findtext("angleX", default="0"))
        ang.angles[1] = float(angles.findtext("angleY", default="0"))
        ang.angles[2] = float(angles.findtext("angleZ", default="0"))
        vv.addVectorAnglesList(ang)
   
    for vecobjs in root.iter('vectorobject'):
        vobj = VectorObject()
        vobj.objName = vecobjs.get('name')
            
    ...

The XML code is also rather lengthy and below is the beginning, defining a lightsource, one set of angles, and one vector object.

<?xml version="1.0"?>
<vectordata>
 <lightsource>
  <!-- set lightsource position for shading. If any object has the "lightsource" property, these will be overridden. -->
  <lightposition>
   <lightpositionX>600</lightpositionX>
   <lightpositionY>200</lightpositionY>
   <lightpositionZ>800</lightpositionZ>
  </lightposition>
 </lightsource>
 <vectorangleslist>
  <!-- set up angle sets for rotations. Values given are starting values. -->
  <vectorangles name="viewer">
   <angleX>0</angleX>
   <angleY>90</angleY>
   <angleZ>0</angleZ>
  </vectorangles>
 </vectorangleslist>
 <vectorobjectslist>
  <!-- set up vectorobjects. Each object must have a position, a reference to angles used, a list of nodes, and a list of surfaces. -->
  <!-- minshade (must be between 0 and 1) determines "shade multiplier" when surface is at 90 degree angle to light source. Use 1.0 to keep full color value ie. no shading. -->
  <!-- use "copyfrom" property (e.g. <vectorobject name="cube2" copyfrom="cube1">) to copy objects; properties like position and colors can then be re-defined, and initangles specified. -->
  <!-- object can have default values for surface color, edgewidth, showback (e.g. "defcolor"). These can be overridden for each surface. --> 
  <vectorobject name="house1_1">
   <!-- house1: a small house with gabled roof -->
   <position>
    <positionX>-120</positionX>
    <positionY>0</positionY>
    <positionZ>200</positionZ>
   </position>
   <anglesref>viewer</anglesref>
   <minshade>0.2</minshade>
   <defcolor>
    <defcolorR>255</defcolorR>
    <defcolorG>255</defcolorG>
    <defcolorB>255</defcolorB>
   </defcolor>
   <defedgewidth>0</defedgewidth>
   <defshowback>0</defshowback>
   <nodelist numnodes="10">
    <node ID="0">
     <nodeX>50</nodeX>
      <nodeY>0</nodeY>
      <nodeZ>-100</nodeZ>
    </node>
    <node ID="1">
     <nodeX>50</nodeX>
      <nodeY>0</nodeY>
      <nodeZ>100</nodeZ>
    </node>
    <node ID="2">
     <nodeX>-50</nodeX>
      <nodeY>0</nodeY>
      <nodeZ>100</nodeZ>
    </node>
    <node ID="3">
     <nodeX>-50</nodeX>
      <nodeY>0</nodeY>
      <nodeZ>-100</nodeZ>
    </node>
    <node ID="4">
     <nodeX>50</nodeX>
      <nodeY>100</nodeY>
      <nodeZ>-100</nodeZ>
    </node>
    <node ID="5">
     <nodeX>50</nodeX>
      <nodeY>100</nodeY>
      <nodeZ>100</nodeZ>
    </node>
    <node ID="6">
     <nodeX>-50</nodeX>
      <nodeY>100</nodeY>
      <nodeZ>100</nodeZ>
    </node>
    <node ID="7">
     <nodeX>-50</nodeX>
      <nodeY>100</nodeY>
      <nodeZ>-100</nodeZ>
    </node>
    <node ID="8">
     <nodeX>0</nodeX>
      <nodeY>125</nodeY>
      <nodeZ>-75</nodeZ>
    </node>
    <node ID="9">
     <nodeX>0</nodeX>
      <nodeY>125</nodeY>
      <nodeZ>75</nodeZ>
    </node>
   </nodelist>
   <surfacelist>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="0"/>
      <node order="1" refID="1"/>
      <node order="2" refID="5"/>
      <node order="3" refID="4"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="1"/>
      <node order="1" refID="2"/>
      <node order="2" refID="6"/>
      <node order="3" refID="5"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="2"/>
      <node order="1" refID="3"/>
      <node order="2" refID="7"/>
      <node order="3" refID="6"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="3"/>
      <node order="1" refID="0"/>
      <node order="2" refID="4"/>
      <node order="3" refID="7"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="4"/>
      <node order="1" refID="5"/>
      <node order="2" refID="9"/>
      <node order="3" refID="8"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="5"/>
      <node order="1" refID="6"/>
      <node order="2" refID="9"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="6"/>
      <node order="1" refID="7"/>
      <node order="2" refID="8"/>
      <node order="3" refID="9"/>
     </nodelist>
    </surface>
    <surface ID="0">
     <nodelist>
      <node order="0" refID="7"/>
      <node order="1" refID="4"/>
      <node order="2" refID="8"/>
     </nodelist>
    </surface>
   </surfacelist>
  </vectorobject>

The object above is a simple house with a gabled roof. We need ten nodes; four for the corners on the ground, four for the corners at the top of the walls, and two for the gabled roof. Furthermore, there are eight surfaces; four for the walls, and four for the roof. The color of each surface can be defined individually, but I have only defined a "defcolor", which is then the default.

The position of the object is its position in the world, and its nodes are defined relative to that being the center of the object. The "anglesref" links the object to an angle set of the same name for rotation.

There is also a copying functionality implemented. By setting e.g.

  <vectorobject name="house2_2" copyfrom="house2_1">
   <position>
    <positionX>-700</positionX>
    <positionY>0</positionY>
    <positionZ>-100</positionZ>
   </position>
  </vectorobject>

the object house2_2 will be a copy of house2_1, but its position will be set separately. It is also possible to override some other properties like surface colors, or to rotate the copied object in place before applying it. This makes re-using object data easier. 

Efficient Programming

There must be tons of books, university courses etc. on how to achieve efficiency in programming. I have not read or attended any of them. Back in the 1980's I taught myself programming first in BASIC (I think Beginners' All-purpose Symbolic Instruction Code says it all) by reading instruction manuals and magazines, and by trial and error; then later Assembly by the same method and the Amiga Hardware Reference Manual and Amiga Machine Language (still available at Amazon!).

To me, there must be at least three definitions of efficiency in programming:

1. how efficiently the program code runs - how quickly it gets the job done.
2. how efficient it is to program - how long a time it takes for the programmer to create the program for the specified purpose.
3. how efficient it is to document and maintain the program - how to keep the program running and add changes when requirements change etc.

Arguably, Assembly language excels in (1) but is miserable in (2) and (3), especially if the programmer is a 1980's self-taught demo coder. (Note that while we demo coders seldom did so, there were libraries etc. we could have used even then and in Assembly.) Python, for example, is probably the opposite - not very efficient code in comparison, but way ahead in programmability and maintainability.

Next: Part V: Ground and Roads.

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.