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.

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.