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.