What is QB64?

What is 3D Rasterizing?

Triangles 101

Barycentric Coordinates

Linear Algebra for the Faint of Heart

Depth is Easier than it Seems

Correcting Barycentric Coordinates with Depth

Putting the Code All Together

Rotation with Optimization

Going Forward with Graphics

Further Resources

2: A single point with two lines jetting out from it, being filled until a terminating line.

3: A set of two triangles. With each of their terminating lines being horizontally flat.

This also makes calculating the x range much simpler. As we can set an initial x value to the x coordinate of the origin, and every time we increment the y value, we add 1/m to our x value. We use 1/m as it is the derivative of "x = (y - origin.x)/m + origin.x" (more specifically the dx/dy). Essentially, the difference between "x = (y - origin.x)/m + origin.x" and "x = (y + 1 - origin.x)/m + origin.x" is just 1/m.

Given that, we simply loop through y values, finding the left and right x values along the way, and then loop between the left and right x in order to get or set every pixel inside of the flat-terminating triangle.

RIGHT = ORIGIN_X

' Left and Right inverse slopes

LINVSLOPE = ( TERMINATOR_LEFT - ORIGIN_X ) / ( TERMINATOR_Y - ORIGIN_Y )

RINVSLOPE = ( TERMINATOR_RIGHT - ORIGIN_X ) / ( TERMINATOR_Y - ORIGIN_Y )

For Y = ORIGIN_Y To TERMINATOR_Y

LEFT += LINVSLOPE

RIGHT += RINVSLOPE

For X = LEFT To RIGHT

Set Pixel ( X , Y )

Next

Next

In order to define the two triangles, we need to find a missing point on the line opposite of our middle point. And we have two things in order to do that, we have the line that the point lies on, and the y position of the point. Using our x = (y - b) / m, we can solve for x. Thinking of our slope in terms of changes in y related to changes in x, we can find our new x with "x = (y - y1) * m + x1". With (y1, x1) being either end point of our target line, and m being the inverse slope; (x2 - x1)/(y2 - y1). This shows us our mystery point to be ((middlePoint.y-y1)*(x2-x1)/(y2-y1) + x1, middlePoint.y). While this is rather hefty in terms of a small computation, it is only needed at maximum once for every time a triangle is drawn.

'Declare inverse slops from a to b , and a to c

a_m = ( x2 - x1 ) / ( term - y1 )

b_m = ( x3 - x1 ) / ( term - y1 )

'Declare the y "scanline" , set it to the terminator and work up to the origin

y = term

'Declare the two x values that walk the lines from the terminator to the origin

lx = x2

rx = x3

'Sign dictates whether rx is in the positive-x or negative-x direction from lx

xsign = Sgn ( rx-lx )

'Same with ysign and y1-terminator

ysign = Sgn ( y1-term )

Do

'Declare the current x position that we will plot a pixel at

x = lx

Do

PSet ( x , y ) 'Set the pixel! ( QB64 method , requires 32 bit screen mode )

'Increment x in the direction of rx

x = x + xsign

'Do this until the current x surpasses rx.

'Sign is checked in order to tell if x surpasses rx when it's greater , or lower

Loop While ( ( xsign = 1 And x < rx ) Or ( xsign = -1 And x > rx ) )

'Increment y towards the origin

y = y + ysign

'Travel lx down line a , ( x2 , bottom ) -> ( x1 , y1 ) , and rx down line b

lx = lx + a_m * ysign

rx = rx + b_m * ysign

Loop While ( ( ysign = 1 And y < y1 ) Or ( ysign = -1 And y > y1 ) )

End Sub

Dim m As Double

'If any two y coordinates are equal , then draw a flat triangle

If y1 = y2 Then

Call flat_triangle ( x3 , y3 , x1 , x2 , y1 ) : GoTo ft2dskip

ElseIf y2 = y3 Then

Call flat_triangle ( x1 , y1 , x2 , x3 , y2 ) : GoTo ft2dskip

ElseIf y3 = y1 Then

Call flat_triangle ( x2 , y2 , x1 , x3 , y1 ) : GoTo ft2dskip

End If

'Find the middle y-coordinate and split the triangle into 2 flat triangles

'And for each point , find the missing x value opposite of the middle

'During the explanation , I used the top point to find the missing midpoint.

'But it can be found even if we interchange the bottom and top point ,

'so we can avoid figuring out which is the top and find the midpoint directly from

'one of the non-middle points

If ( y1 < y2 And y1 > y3 ) Or ( y1 > y2 And y1 < y3 ) Then

m = x2 + ( y1 - y2 ) * ( x3 - x2 ) / ( y3 - y2 )

Call flat_triangle ( x3 , y3 , x1 , m , y1 )

Call flat_triangle ( x2 , y2 , x1 , m , y1 )

GoTo ft2dskip

ElseIf ( y2 < y1 And y2 > y3 ) Or ( y2 > y1 And y2 < y3 ) Then

m = x1 + ( y2 - y1 ) * ( x3 - x1 ) / ( y3 - y1 )

Call flat_triangle ( x3 , y3 , x2 , m , y2 )

Call flat_triangle ( x1 , y1 , x2 , m , y2 )

GoTo ft2dskip

ElseIf ( y3 < y1 And y3 > y2 ) Or ( y3 > y1 And y3 < y2 ) Then

m = x1 + ( y3 - y1 ) * ( x2 - x1 ) / ( y2 - y1 )

Call flat_triangle ( x2 , y2 , x3 , m , y3 )

Call flat_triangle ( x1 , y1 , x3 , m , y3 )

End If

ft2dskip:

End Sub

y1 is middle

y2 is bottom (y-positive is down)

y3 is top

Etc.

The calculation of barycentric coordinates, notated as UV and W.

$$\left(U, V, W\right) = \left(\frac{\mathrm{area of} ΔPBC}{\mathrm{area of} ΔABC}, \frac{\mathrm{area of} ΔAPC}{\mathrm{area of} ΔABC}, \frac{\mathrm{area of} ΔABP}{\mathrm{area of} ΔABC}\right).$$

And for our color mixing example,

$$COLOR = RED * U + GREEN * V + BLUE * W$$

When taking the barycentric coordinates (u, v, w), we procedurally check new points of the triangle with the point we want to sample inside the triangle. In our process, that would be the x and y of the pixel we are setting; \(p_x\) and \(p_y\). In order to preserve vertex order, we must write the calculations in the way specified below. Think of it as, for every next barycentric coordinate, we increment every point we're subtracting p from; e.i. \(x_1 -> x_2; x_2 -> x_3; x_3 -> x_1\). The following gives us \(a\) as the full area of the triangle triangle, as well as \(u\), \(v\), and \(w\) as the barycentric coordinates, in the range (0-1). $$a = \frac{\left(x_3-x_1\right) * \left(y_2-y_1\right) - \left(y_3-y_1\right) * \left(x_2-x_1\right)}{2}$$ $$u = \frac{\left(x_3-p_x\right) * \left(y_2-p_y\right) - \left(y_3-p_y\right) * \left(x_2-p_x\right)}{2a}$$ $$v = \frac{\left(x_1-p_x\right) * \left(y_3-p_y\right) - \left(y_1-p_y\right) * \left(x_3-p_x\right)}{2a}$$ $$w = \frac{\left(x_2-p_x\right) * \left(y_1-p_y\right) - \left(y_2-p_y\right) * \left(x_1-p_x\right)}{2a}$$

Essentially, the attributes of any given point or pixel inside of a triangle is equal to the sum the attributes of each vertice, multiplied by their associated barycentric component. I.e; $$pixel_{attribute} = a_{attribute} * u + b_{attribute} * v + c_{attribute} * w$$ Here is a visual example, where vertex \(a\) has a color attribute of red, \(b\) with green, and \(c\) with blue.

area = ( ( cx - ax ) * ( by - ay ) - ( cy - ay ) * ( bx - ax ) ) * 0.5 ' Compute cross product divided by two

u = ( cx - px ) * ( by - px ) - ( cy - pc ) * ( bx - px )

getU = u / 2*area

End Function

Function getV ( ax , ay , bx , by , cx , cy , px , py )

area = ( ( cx - ax ) * ( by - ay ) - ( cy - ay ) * ( bx - ax ) ) * 0.5

v = ( ax - px ) * ( cy - px ) - ( ay - pc ) * ( cx - px )

getV = v / 2*area

End Function

Function getW ( ax , ay , bx , by , cx , cy , px , py )

area = ( ( cx - ax ) * ( by - ay ) - ( cy - ay ) * ( bx - ax ) ) * 0.5

w = ( bx - px ) * ( ay - px ) - ( by - pc ) * ( ax - px )

getW = w / 2*area

End Function

'For every pixel in triangle rasterizing...

'x1-2-3 and y1-2-3 are the triangle's vertices , with x and y being the current pixel's coords

Color _RGB ( 255 * getU ( x1 , y1 , x2 , y2 , x3 , y3 , x , y ) , 255 * getV ( x1 , y1 , x2 , y2 , x3 , y3 , x , y ) , 255 * getW ( x1 , y1 , x2 , y2 , x3 , y3 , x , y ) )

PSet ( x , y )

'Paints the pixels of a triangle Red , green , and blue based on the barycentric coords.

These transformations happen through matrix multiplication. You can plug any vector or coordinate in and get a resulting coordinate that follows the rules defined by the matrix. This is extremely useful in 3d graphics as it allows for all sorts of transformations and even traslations. We could then create a matrix that rotates all vectors, or one that projects all input 3d vectors into 2 dimension. This idea of projection is extremely important, and will be discussed again shortly.

Our projector will transform these coordinates into "screenspace coordinates", coordinates ranging from (-1, -1, -1) to (1, 1, 1), in a way that preserves the window's display ratio (i.e. 16:9, 4:3, 2:1, etc.) After doing so, it can bring the coordinates from 3d into 2d in a way that shows foreshortening. From there it must convert this range of numbers to actual pixel coordinates in order to show the image we desire.

The first transforming Matrix we will use, is; $$\begin{bmatrix}\color{red}F_x & 0 & 0 & 0\\0 & \color{red}F_y & 0 & 0\\0 & 0 & \color{yellow}Z_m & \color{yellow}Z_a \\ 0 & 0 & \color{blue}-1 & 0 \\\end{bmatrix}$$ For foreshortening, we will use: $$\begin{bmatrix}\color{red}1/z&0\\0&\color{red}1/z\end{bmatrix}$$ Or, simply: $$x' = x / z$$ $$y' = y / z$$ And finally, transforming the coordinates into pixel coordinates; (The input z is required to be 1; \(S_w\) and \(S_h\) are the width and height of whatever is being displayed to.) $$\begin{bmatrix}\color{red}\frac{1}{2}S_w & 0 & \color{yellow}\frac{1}{2}S_w \\ 0 & \color{red}\frac{1}{2}S_h & \color{yellow}\frac{1}{2}S_h \\ 0 & 0 & 1 \end{bmatrix}$$ Or even more simply; $$x' = \frac{S_w}{2}\left(x + 1\right)$$ $$y' = \frac{S_h}{2}\left(y + 1\right)$$ If we expand each of these transformations, we may simplify things down into an easy-to-digest calculation.

Our first matrix sets the output x to \(x * F_x\), and output y to \(y * F_y\). \(F_x\) and \(F_y\) scale the x and y values based on two factors.

Firstly, the width-height ratio of the display.

And second, the prefered field of view, or FOV.

Since we are transforming a physical space into a screen space and then stretching it onto the screen, we need to be prepared for any squishing or stretching that happens as our (-1,-1)-(1, 1) space transforms into a (0, 0)-(1920, 1080) or whatnot. If say, our display has a width-height ratio of 16:9, just like in one that's 1920 - 1080, the x values are being stretched out 16/9 times more than the y is being stretched by. In order to combat this, we can either stretch the y by width/height, or stretch the x by height/width. E.i; \(F_x = \frac{h}{w}\) or \(F_y = \frac{w}{h}\) (But not both).

But there is more to add for \(F_x\) and \(F_y\). The FOV is a scalar that scales everything equally, and this scale is chosen rather than found. But that doesn't mean it can be any number, it technically could but in order to create the effect of a "Field of View" you can calculate the FOV scalar from and angle through \(FOV = \cot(\frac{1}{2}\theta)\). Now we could combine this with the w/h ratio by directly multiplying, but instead we will multiply the \(\theta\) angle instead. Thus we can put this together to create our final \(F_x\) and \(F_y\) values. $$F_x = \cot(\frac{1}{2}\theta_{FOV})$$ $$F_y = \cot(\frac{1}{2}\theta_{FOV} \times \frac{W}{H})$$ Now we've transformed our x and y values, what more can we do? Well, \(Z_m\) and \(Z_a\) in the matrix are for setting a maximum and minimum z-value, otherwise known as the Far and Near Clipping Planes. When transforming our z-coordinates, e.i. distance to the viewpoint, we want to map Z from a range of (near-clip, far-clip) to a range of (0, far-clip). Doing this just takes multiplication to stretch or squish z-values to the right range of values. And an addition(/subtraction) to set the proper minimum value 0. $$ Z_m = \mp\frac{F}{F-N} $$ $$ Z_a = -\frac{F \times N}{F-N} $$ Whether or not \(Z_m\) is positive or negative depends on which z-direction is forward. It's typically negative, as, following the right hand coordinate system, Z-negative is "forward". \(Z_a\) is ultimately added to the input z (after multiplying it by \(Z_m\)) as, writing out the multiplication of our matrix, we get, \(z' = x \times 0 + y \times 0 + z \times Z_m + w \times Z_a\) where we've already established that w = 1. And such, we have a final Matrix of $$\begin{bmatrix}cot(\frac{\theta}{2}) & 0 & 0 & 0 \\ 0 & cot(\frac{\theta}{2} \times \frac{W}{H}) & 0 & 0 \\ 0 & 0 & -\frac{F}{F-N} & -\frac{F \times N}{F-N} \\ 0 & 0 & -1 & 0 \end{bmatrix}$$ Where \(\theta\) is the angle of FOV, \(W\) and \(H\) are the dimensions of the display, and \(F\) and \(N\) being the Near and Far clipping planes respectively.

You'll notice that there's also a -1 in the \(w' = x \times 0 + y \times 0 + z * -1 + w \times 0\) spot. This is used to store distance, as z can be thought of the distance between a point and the viewpoint. But the 'forward' direction, is z-negative, hence why we're multiplying z by -1 instead of just 1. In the end, w is the distance of a point to the viewpoint. If it's negative, it's behind the viewing area (z-positive direction), and thus shouldn't be drawn.

We have now created a selection of geometry to draw. Multiplying a vector with this matrix, it tells us whether a point will (or should) be drawn. That being, outputs with x and y values between -1 and 1, z-values between 0 and the far-clipping plane, and positive w values.

That is all the harder to grasp concepts out of the way. It also prepares us for foreshortening, as, in order to create this effect of foreshortening, we simply divide our x and y values by z. This means that if a point has a small z value below 1, dividing by z will stretch it out, making an object seem to get bigger when really close. But as z increases, dividing by z will make x and y much smaller, making objects appear smaller as they are further. With this in mind, we can write out our final computations. Check back earlier to matrix multiplication if you need to remember how to do this step.

$$ z' = z \times Z_m + Z_a $$ Done first as we need to divide x and y by it. $$ x' = \frac{x}{z'} \times cot(\frac{\theta}{2}) $$ $$ y' = \frac{y}{z'} \times cot(\frac{\theta}{2} \times \frac{W}{H}) $$ $$ w' = z * -1$$ Notice that this is z, not z'.

These gives us 2D foreshortened x and y values \(x'\) and \(y'\). We can translate these to pixel coordinates simply, but through this method any x or y outside of the (-1, 1) range will be off the screen, therefore we should check for these cases to make sure we do not attempt to wrongfully draw these cases. $$ x'' = (x' \times 0.5 + 0.5) \times W $$ $$ y'' = (y' \times 0.5 + 0.5) \times H $$ And thus (x'', y'') are our final pixel coordinates, with w' being the point's distance from the viewpoint.

Putting these concepts together, we can write code for a Function (Sub) that takes in the 3D coordinates of a point, and fills in the pixel of said point.

W = 720

H = 480

bscreen = _NewImage ( W , H , 32 ) '32 Bit color mode

Screen bscreen

Dim Shared Zm As Double , Za As Double , FOV As Double

theta = 90 * 3.14159265 / 180 ' 90 is the FOV angle in degrees , theta is in radians

FOV = 1 / Tan ( theta / 2 )

Far = 100

Near = 0.01

Zm = -Far / ( Far - Near )

Za = - ( Far * Near ) / ( Far - Near )

currentTime# = 0

Do

'Cls

Color _RGB ( Cos ( currentTime# ) * 127 + 127 , Sin ( currentTime# ) * 127 + 127 , 127 + 127 * Cos ( currentTime# * 2.718281828459 + 0.5 ) )

Call plot3DPoint ( Cos ( currentTime# ) * 50 , Sin ( currentTime# ) * 50 , -100 + 50 * Cos ( currentTime# * 2.718281828459 + 0.5 ) )

currentTime# = currentTime# + 0.001

_Display

Loop

Function projectZ ( z As Double )

projectZ = z * Zm + Za

End Function

Function projectX ( x As Double , z As Double )

projectX = ( ( x / z ) * FOV * 0.5 + 0.5 ) * W

End Function

Function projectY ( y As Double , z As Double )

projectY = ( y / z ) * FOV * W * 0.5 + 0.5 * H

'Originally , ( y/z * FOV * W/H * 0.5 + 0.5 ) *H , but simplified a little to save a division

End Function

Sub plot3DPoint ( x As Double , y As Double , z As Double )

If -z < 0 Then GoTo skip ' ( -z ) is w , checking if it is behind the camera

zPrime = projectZ ( z )

xPrime = projectX ( x , zPrime )

If xPrime < 0 Or xPrime > W Then GoTo skip 'If x position is outside screen , don't draw

yPrime = projectY ( y , zPrime )

If yPrime < 0 Or yPrime > H Then GoTo skip 'If y position is outside screen , don't draw

PSet ( xPrime , yPrime ) 'Requires QB64 Screen , like at the beginning of this program

skip:

End Sub

In order to prevent triangles from being placed in front of trangles they're supposed to behind, we do something called depth testing. Depth Testing, as we will use it, will be the storage of two arrays, with an element for every pixel. The first array will store the distance of a pixel from the viewpoint. The second, will store a number representing the last frame a pixel was updated.

I.e, our code for this should execute something along these lines; (W and H are display dimensions)

timeArray = integer array with size [W * H]

current frame = 0

For every frame/every update

increment current frame

Draw triangles and whatnot

For every pixel in a triangle

find index of current pixel

If depthArray[pixel index] < current depth or timeArray[pixel index] < current frame

depthArray[pixel index] = current depth

timeArray[pixel index] = current frame

Draw pixel

The index of a pixel is rather simple. We have arrays with room for every pixel. Width times height, area of a rectangle. Think about labeling every pixel from left to right and top to bottom. When you start, you are labeling every pixel 0, 1, 2, 3, etc. and if you think about it, you are essentially labeling every pixel with its x position. When you have hit the end of the first line, you wrap back around. The first element of the second line will have its index be the same as the Width of the display. That is because we started with index 0, so when we hit the pixel at the end of the first row, its index is Width - 1. If we go another row, our pixel at (0, 2) will be 2 * Width. After another row it's 3 * width, and etc. With each pixel, we have a unique index equal to (x + y * Width). And this is how we find the index of a pixel.

So pixel index = pixel x + pixel y * display width. What about pixel depth?

Well here's a hint, through our projection, we have the depth of the vertices of the triangle we want to draw, and we need to interpolate to find the depth per pixel. Think barycentric coordinates.

Unfortunately, unlike colors and other linear vertex attributes, we can't just do $$depth = depth_a * u + depth_b * v + depth_c * w$$ We instead need to do $$\frac{1}{depth} = u\frac{1}{depth_a} + v\frac{1}{depth_b} + w\frac{1}{depth_c}$$ You'll find that by attempting to use the first equation over the latter to calculate depth will result in strange curves when geometry intersects eachother. The reasoning for this lies in how depth actually effects coordinates. Each x and y value are divided by their depth, this ends in every x and y value not being transformed equally.

Shouldn't this depth distortion affect other attributes? Well it does, and you can see the effects and solutions in the next section

We have the ability to take the depth of every pixel, it's time we should compare them. You may find the merging of pixel indexing and depth interpolation in this following code example.

timeArray = integer array; size [W*H]

Do every frame

increment current frame

Draw triangles and whatnot

Loop

For every pixel ( x , y ) in a triangle ( x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 )

Where z1 , z2 , z3 are the depths of their respective vertex

u = getU ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

v = getV ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

w = getW ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

depth = 1 / ( u / z1 + v / z2 + w / z3 )

index = y * W + x

If depthArray[index] < depth or timeArray[index] < current frame

depthArray[index] = depth

timeArray[index] = current frame

Draw pixel ( x , y )

This problem is more apparent when using textures because images will hold straight lines, whereas simple gradients might hide misteps.

Equations that show what's going on here would look like this. Where d is the depth of each vertice and of the current pixel, and a is an attribute of a vertex like color or whatnot, and our barycentric coordinates, u, v, and w. $$d_{pixel} = \frac{1}{\frac{u}{d_1} + \frac{v}{d_2} + \frac{w}{d_3}}$$ $$a_{pixel} = \left(u\frac{a_1}{d_1} + v\frac{a_2}{d_2} + w\frac{a_3}{d_3}\right)d_{pixel}$$ This helps as the attribute not only relies on barycentric coordinates, but also depth. A pixel will not just appear like a vertice as it's barycentric coords approach it, but its depth needs to approach it as well.

The following hypothetical code will

Where z1 , z2 , z3 are the depth of each vertex

u = getU ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

v = getV ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

w = getW ( x1 , y1 , x2 , y2 , x3 , y3 , x , y )

depth = 1 / ( u/z1 + v/z2 + w/z3 )

color = depth * ( RED * u/z1 + GREEN * v/z2 + BLUE * w/z3 )

W = 720

H = 480

bscreen = _NewImage ( W , H , 32 )

Screen bscreen

Dim Shared pixelDepth ( W * H ) As Double

Dim Shared pixelTime ( W * H ) As Double

Dim Shared currentFrame As Integer

currentFrame = 0

Dim Shared FOV As Double , Zm As Double , Za As Double

theta# = 90 * 3.14159265 / 180.0 ' 90 is the FOV angle in degrees , theta is in radians

FOV = 1 / Tan ( theta# * 0.5 )

Far = 10000

Near = 0.01

Zm = -Far / ( Far - Near )

Za = - ( Far * Near ) / ( Far - Near )

'Declare variables for numbers used in multiple Functions and Subs

Dim Shared triX1 , triY1 , triZ1 , triX2 , triY2 , triZ2 , triX3 , triY3 , triZ3

Do

Cls

currentFrame = currentFrame + 1

Call triangle ( -50 , -50 , -100 , 50 , 50 , -100 , 100 , -75 , -200 )

_Display

Loop

Function projectZ ( z As Double )

projectZ = z * Zm + Za

End Function

Function projectX ( x As Double , z As Double )

projectX = ( ( x / z ) * FOV * 0.5 + 0.5 ) * W

End Function

Function projectY ( y As Double , z As Double )

projectY = ( y / z ) * FOV * W * 0.5 + 0.5 * H

End Function

Function getU ( px , py )

' Compute cross product divided by two

area = ( ( triX3 - triX1 ) * ( triY2 - triY1 ) - ( triY3 - triY1 ) * ( triX2 - triX1 ) ) * 0.5

getU = ( ( triX3 - px ) * ( triY2 - py ) - ( triY3 - py ) * ( triX2 - px ) ) / 2 * area

End Function

Function getV ( px , py )

area = ( ( triX3 - triX1 ) * ( triY2 - triY1 ) - ( triY3 - triY1 ) * ( triX2 - triX1 ) ) * 0.5

getV = ( ( triX1 - px ) * ( triY3 - py ) - ( triY1 - py ) * ( triX3 - px ) ) / 2 * area

End Function

Function getW ( px , py )

area = ( ( triX3 - triX1 ) * ( triY2 - triY1 ) - ( triY3 - triY1 ) * ( triX2 - triX1 ) ) * 0.5

getW = ( ( triX2 - px ) * ( triY1 - py ) - ( triY2 - py ) * ( triX1 - px ) ) / 2 * area

End Function

Sub rasterFlatTriangle ( x1 , y1 , x2 , x3 , term ) 'z1 , z2 , and z3 are the depths of each point

'Declare inverse slops from a to b , and a to c

Dim y As Integer , x As Integer

Dim a_m , b_m , lx , rx , xsign , ysign , bu , bv , bw , depth , index

a_m = ( x2 - x1 ) / ( term - y1 )

b_m = ( x3 - x1 ) / ( term - y1 )

'Declare the y "scanline" , set it to the terminator and work up to the origin

y = term

'Declare the two x values that walk the lines from the terminator to the origin

lx = x2

rx = x3

'Signs dictate which directions to increment in , +1 or -1

xsign = Sgn ( rx - lx )

ysign = Sgn ( y1 - term )

Do

x = lx

If y < 0 Or y >= H Then

GoTo skipy

End If

Do

bu = getU ( x , y )

bv = getV ( x , y )

bw = getW ( x , y )

depth = 1 / ( bu / triZ1 + bv / triZ2 + bw / triZ3 )

index = ( y * W ) + x

If x < 0 Or x >= W Then

GoTo skipx

End If

If pixelDepth ( index ) < depth Or pixelTime ( index ) < currentFrame Then

pixelDepth ( index ) = depth

pixelTime ( index ) = currentFrame

Color _RGB ( 255 * bu / triZ1 * depth , 255 * bv / triZ2 * depth , 255 * bw / triZ3 * depth )

'Color _RGB ( 255 , 255 , 0 )

PSet ( x , y ) 'Set the pixel! ( QB64 method , requires 32 bit screen mode )

End If

skipx:

x = x + xsign

Loop While ( ( xsign = 1 And x < rx ) Or ( xsign = -1 And x > rx ) )

skipy:

y = y + ysign

lx = lx + a_m * ysign

rx = rx + b_m * ysign

Loop While ( ( ysign = 1 And y < y1 ) Or ( ysign = -1 And y > y1 ) )

End Sub

Sub triangle ( x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 )

'Find 2d coordinates + depth from original coordinates

Dim m

triZ1 = projectZ ( z1 )

triX1 = projectX ( x1 , triZ1 )

triY1 = projectY ( y1 , triZ1 )

triZ2 = projectZ ( z2 )

triX2 = projectX ( x2 , triZ2 )

triY2 = projectY ( y2 , triZ2 )

triZ3 = projectZ ( z3 )

triX3 = projectX ( x3 , triZ3 )

triY3 = projectY ( y3 , triZ3 )

'If any two y coordinates are equal , then draw a flat triangle

If triY1 = triY2 Then

Call rasterFlatTriangle ( triX3 , triY3 , triX1 , triX2 , triY1 ) : GoTo triangleskip

ElseIf triY2 = triY3 Then

Call rasterFlatTriangle ( triX1 , triY1 , triX2 , triX3 , triY2 ) : GoTo triangleskip

ElseIf triY3 = triY1 Then

Call rasterFlatTriangle ( triX2 , triY2 , triX1 , triX3 , triY1 ) : GoTo triangleskip

End If

'Find the middle y-coordinate and split the triangle into 2 flat triangles

If ( triY1 < triY2 And triY1 > triY3 ) Or ( triY1 > triY2 And triY1 < triY3 ) Then

m = triX2 + ( triY1 - triY2 ) * ( triX3 - triX2 ) / ( triY3 - triY2 )

Call rasterFlatTriangle ( triX3 , triY3 , triX1 , m , triY1 )

Call rasterFlatTriangle ( triX2 , triY2 , triX1 , m , triY1 )

GoTo triangleskip

ElseIf ( triY2 < triY1 And triY2 > triY3 ) Or ( triY2 > triY1 And triY2 < triY3 ) Then

m = triX1 + ( triY2 - triY1 ) * ( triX3 - triX1 ) / ( triY3 - triY1 )

Call rasterFlatTriangle ( triX3 , triY3 , triX2 , m , triY2 )

Call rasterFlatTriangle ( triX1 , triY1 , triX2 , m , triY2 )

GoTo triangleskip

ElseIf ( triY3 < triY1 And triY3 > triY2 ) Or ( triY3 > triY1 And triY3 < triY2 ) Then

m = triX1 + ( triY3 - triY1 ) * ( triX2 - triX1 ) / ( triY2 - triY1 )

Call rasterFlatTriangle ( triX2 , triY2 , triX3 , m , triY3 )

Call rasterFlatTriangle ( triX1 , triY1 , triX3 , m , triY3 )

End If

triangleskip:

End Sub

Firstly, one of the biggest changes you can make to improve the look of a 3d program is to introduce rotation. Without it you would be stuck moving a still shape across the screen. Rotation is mathematically tedius, but it's much easier to understand with Matrix Transformation. If you take a vector as a point in 3d space, you can rotate it around the origin (0, 0, 0) with the following (column vector) matrices. $$R_x(\theta) = \begin{bmatrix}1 & 0 & 0\\0 & cos\theta & sin\theta\\0 & -sin\theta & cos\theta\end{bmatrix}$$ $$R_y(\theta) = \begin{bmatrix}cos\theta & 0 & -sin\theta\\0 & 1 & 0\\sin\theta & 0 & cos\theta\end{bmatrix}$$ $$R_z(\theta) = \begin{bmatrix}cos\theta & sin\theta & 0\\-sin\theta & cos\theta & 0\\0 & 0 & 1\end{bmatrix}$$ What each of these mean is that, say for \(R_x(\theta)\), given angle \(\theta\), the transformation will rotate a vector around the X axis by this angle. Thus there are 3 axes of rotation, each around their respective 3D axis.

(In order left-to-right: X-rotation; Y-rotation; Z-rotation)

Combining these matrices, you end up with one matrix you can use to rotate any vector in all 3 axes at the same time; $$R_{xyz}(\alpha, \beta, \gamma) = \begin{bmatrix}\cos\alpha \cos\beta & \sin\alpha \cos\beta & -\sin\beta \\ \cos\alpha \sin\beta \sin\gamma - \sin\alpha \cos\gamma & \sin\alpha \sin\beta \sin\gamma + \cos\alpha \cos\gamma & \cos\beta \sin\gamma \\ \cos\alpha \sin\beta \cos\gamma + \sin\alpha \sin\gamma & \sin\alpha \sin\beta \cos\gamma - \cos\alpha \sin\gamma & \cos\beta \cos\gamma \end{bmatrix}$$ Where \(\alpha\) is the angle around the Z axis, \(\beta\) around the Y axis, and \(\gamma\) around the X axis. Implementing this into BASIC is rather easy, but it is also rather tedious as there is no easy way to simplify the computations. However, we can make our code a little neater and clean using a User Defined Type (UDT). Like a struct or class for QBASIC. We will create a UDT called "vector". It will store 4 numbers, for an x, y, z and w coordinate.

And as we cannot return more than one number per function, (i.e, we cannot return a UDT in QB64, we can only return a built-in type, like a number or string), we will set a universal vector the be the set "output" of some functions. Given all that, let me show you what our rotation function will look like, with a vector UDT defined.

x As Single

y As Single

z As Single

w As Single

End Type

Dim Shared vector As vector

'Set the position you want to rotate

vector.x = 1

vector.y = 2

vector.z = -1

vector.w = 0

'Call the rotation function , the given parameters will rotate it 90 degrees around the y axis

Call rotateVector ( 0 , 3.1459265 / 2 , 0 )

Print vector.x , vector.y , vector.z , vector.w

' Prints 1 2 1 0

Sub rotateVector ( rx As Single , ry As Single , rz As Single )

Dim nvector As vector

nvector.x = ( vector.x * Cos ( rz ) * Cos ( ry ) ) + ( vector.y * Sin ( rz ) * Cos ( ry ) ) - vector.z * Sin ( ry )

nvector.y = ( vector.x * ( Cos ( rz ) * Sin ( ry ) * Sin ( rx ) - Sin ( rz ) * Cos ( rx ) ) ) + ( vector.y * ( Sin ( rz ) * Sin ( ry ) * Sin ( rx ) + Cos ( rz ) * Cos ( rx ) ) ) + ( vector.z * Cos ( ry ) * Sin ( rx ) )

nvector.z = ( vector.x * ( Cos ( rz ) * Sin ( ry ) * Cos ( rx ) + Sin ( rz ) * Sin ( rx ) ) ) + ( vector.y * ( Sin ( rz ) * Sin ( ry ) * Cos ( rx ) - Cos ( rz ) * Sin ( rx ) ) ) + ( vector.z * Cos ( ry ) * Cos ( rx ) )

vector.x = nvector.x: vector.y = nvector.y: vector.z = nvector.z: vector.w = 1

End Sub

For every triangle:

Call projection for every coordinate of all 3 points

Find flat triangles

For each flat triangle:

Calculate the slopes of the sides of the flat triangle

For every y value in the triangle's y value range:

Increment ends of x range by their slopes

For every pixel in the triangle:

Calculate u, v, and w barycentric coordinates

Calculate depth (and index) of current pixel

Compare to depth and time array

Draw pixel

Through this we can see that the heaviest load is during the pixel drawing process. As those calculations are repeated the most often: For every pixel in every triangle. But it seems as though we can't clean it up any further without sacrificing depth testing and barycentric coordinates. However, we can still simplify it using derivatives.

Take our barycentric coordinate calculations. They seem pretty complicated, however note that the only thing changing in them is the x and the y position of the pixel. And each time the values change, they change by exactly 1. If we use the derivative of these equations in respect to x and y, and depending on what we're incrementing at that current moment, add them to running barycentric coordinate variables.

Given our UVW calculations:

$$u = \frac{\left(x_3-p_x\right) * \left(y_2-p_y\right) - \left(y_3-p_y\right) * \left(x_2-p_x\right)}{2a}$$ $$v = \frac{\left(x_1-p_x\right) * \left(y_3-p_y\right) - \left(y_1-p_y\right) * \left(x_3-p_x\right)}{2a}$$ $$w = \frac{\left(x_2-p_x\right) * \left(y_1-p_y\right) - \left(y_2-p_y\right) * \left(x_1-p_x\right)}{2a}$$

We find that the derivatives of such in respect to x and y are as follows:

$$\frac{du}{dx} = \frac{y_3-y_2}{2a}$$ $$\frac{dv}{dx} = \frac{y_1-y_3}{2a}$$ $$\frac{dw}{dx} = \frac{y_2-y_1}{2a}$$ $$\frac{du}{dy} = \frac{x_2-x_3}{2a}$$ $$\frac{dv}{dy} = \frac{x_3-x_1}{2a}$$ $$\frac{dw}{dy} = \frac{x_1-x_2}{2a}$$ Given that these computations only rely on the positions of the three vertices of the triangle, they can be calculated once per triangle. Whenever we increment x, we must increment u by 'du/dx' times what our change in x is, typically 1, and the same goes for increments in y with 'du/dy'.

It is possible to take and use the derivative of the depth calculation, however it would be more trouble than it's worth.

Given rotation and our optimization, below you'll find yet another fully working version of this project. Using arrow keys to control rotation.

x As Double

y As Double

z As Double

w As Double

End Type

Dim Shared As Integer Width , Height

Dim Shared As Single triX1 , triX2 , triX3 , triY1 , triY2 , triY3 , tri1OverZ1 , tri1OverZ2 , tri1OverZ3 , triArea , tri1OverArea , triU , triV , triW

Dim Shared vector as vector

Dim Shared projection_data ( 4 ) As Double ' Store FOV * 0.5 , FOV * Width * 0.5 , Zm , and Za

Dim Shared As Single lslope , rslope , lx , rx , sy , ey , dudx , dvdx , dwdx , dudy , dvdy , dwdy , pixel_index , depth

Dim Shared pixel_x As Integer

Dim Shared pixel_y As Integer

Dim Shared current_frame As Integer

current_frame = 1

Width = 800

Height = 800

Dim Shared pixel_depth ( Width * Height ) As Double

Dim Shared pixel_time ( Width * Height ) As Double

display = _NewImage ( Width , Height , 32 )

Screen display

Call setProjection ( 90 , 10000 , 0.01 )

positionx = 0

positiony = 0

positionz = -200

rotationx = 0.0

rotationy = 0.0

Do

Cls

If _KeyDown ( 19200 ) Then

rotationy = rotationy - 0.03

End If

If _KeyDown ( 19712 ) Then

rotationy = rotationy + 0.03

End If

If _KeyDown ( 18432 ) Then

rotationx = rotationx + 0.03

End If

If _KeyDown ( 20480 ) Then

rotationx = rotationx - 0.03

End If

Call cube ( positionx , positiony , positionz , 50 , rotationx , rotationy )

current_frame = current_frame + 1

_Display

Loop

Color _RGB ( 255 , 255 , 255 )

End

Sub cube ( positionx As Single , positiony As Single , positionz As Single , size As Single , rotationx As Single , rotationy As Single )

Dim As vector lft , rft , lbt , rbt , lfb , rfb , lbb , rbb

setVector -size , size , -size , 0

rotateVector rotationx , rotationy , 0

lft.x = vector.x + positionx: lft.y = vector.y + positiony: lft.z = vector.z + positionz

setVector size , size , -size , 0

rotateVector rotationx , rotationy , 0

rft.x = vector.x + positionx: rft.y = vector.y + positiony: rft.z = vector.z + positionz

setVector -size , size , size , 0

rotateVector rotationx , rotationy , 0

lbt.x = vector.x + positionx: lbt.y = vector.y + positiony: lbt.z = vector.z + positionz

setVector size , size , size , 0

rotateVector rotationx , rotationy , 0

rbt.x = vector.x + positionx: rbt.y = vector.y + positiony: rbt.z = vector.z + positionz

setVector -size , -size , -size , 0

rotateVector rotationx , rotationy , 0

lfb.x = vector.x + positionx: lfb.y = vector.y + positiony: lfb.z = vector.z + positionz

setVector size , -size , -size , 0

rotateVector rotationx , rotationy , 0

rfb.x = vector.x + positionx: rfb.y = vector.y + positiony: rfb.z = vector.z + positionz

setVector -size , -size , size , 0

rotateVector rotationx , rotationy , 0

lbb.x = vector.x + positionx: lbb.y = vector.y + positiony: lbb.z = vector.z + positionz

setVector size , -size , size , 0

rotateVector rotationx , rotationy , 0

rbb.x = vector.x + positionx: rbb.y = vector.y + positiony: rbb.z = vector.z + positionz

tri_vector lft , rft , rfb

tri_vector rfb , lfb , lft

tri_vector rbt , lbt , lbb

tri_vector lbb , rbb , rbt

tri_vector lbt , lft , lfb

tri_vector lfb , lbb , lbt

tri_vector rft , rbt , rbb

tri_vector rbb , rfb , rft

tri_vector rft , lbt , rbt

tri_vector lbt , rft , lft

tri_vector lfb , rbb , lbb

tri_vector rbb , lfb , rfb

End Sub

' Returns the max between two numbers

Function max ( ma_x , ma_y )

max = ma_x

If ma_y > ma_x Then

max = ma_y

End If

End Function

' Returns the min between two numbers

Function min ( mi_x , mi_y )

min = mi_x

If mi_y < mi_x Then

min = mi_y

End If

End Function

Sub setVector ( x As Double , y As Double , z As Double , w As Double )

vector.x = x

vector.y = y

vector.z = z

vector.w = w

End Sub

Sub setProjection ( FOVAngle As Double , Far As Double , Near As Double )

projection_data ( 0 ) = 0.5 / Tan ( FOVAngle * 3.14159267 / 360.0 )

projection_data ( 1 ) = Width * projection_data ( 0 )

projection_data ( 2 ) = -Far / ( Far - Near )

projection_data ( 3 ) = - ( Far * Near ) / ( Far - Near )

End Sub

Sub projectVector ' Leaves z as 1/z

vector.z = 1 / ( vector.z * projection_data ( 2 ) + vector.w * projection_data ( 3 ) )

vector.x = vector.x * vector.z * projection_data ( 1 ) + vector.w * 0.5 * Width

vector.y = vector.y * vector.z * projection_data ( 1 ) + vector.w * 0.5 * Height

End Sub

Sub rotateVector ( rx As Double , ry As Double , rz As Double )

Dim nvector As vector

nvector.x = ( vector.x * Cos ( rz ) * Cos ( ry ) ) + ( vector.y * Sin ( rz ) * Cos ( ry ) ) - vector.z * Sin ( ry )

nvector.y = ( vector.x * ( Cos ( rz ) * Sin ( ry ) * Sin ( rx ) - Sin ( rz ) * Cos ( rx ) ) ) + ( vector.y * ( Sin ( rz ) * Sin ( ry ) * Sin ( rx ) + Cos ( rz ) * Cos ( rx ) ) ) + ( vector.z * Cos ( ry ) * Sin ( rx ) )

nvector.z = ( vector.x * ( Cos ( rz ) * Sin ( ry ) * Cos ( rx ) + Sin ( rz ) * Sin ( rx ) ) ) + ( vector.y * ( Sin ( rz ) * Sin ( ry ) * Cos ( rx ) - Cos ( rz ) * Sin ( rx ) ) ) + ( vector.z * Cos ( ry ) * Cos ( rx ) )

nvector.w = 1

vector.x = nvector.x: vector.y = nvector.y: vector.z = nvector.z

End Sub

Sub tri_flat_term ( x1 , y1 , x2 , x3 , term )

lx = x2 ' Set left x that will be walking line ( x2 , term ) - ( x1 , y1 )

rx = x3 ' Set right x that will be walking line ( x3 , term ) - ( x1 , y1 )

If ( x2 > x3 ) Then

rx = x2: lx = x3

End If

lslope = ( lx - x1 ) / ( term - y1 ) ' Negative left inverse slope depending on if y will be increasing or decreasing

rslope = ( rx - x1 ) / ( term - y1 ) ' Negative right inverse slope for same reason

' Set starting y and ending y

sy = term

ey = y1

If ( y1 < term ) Then

sy = y1: ey = term: lx = x1: rx = x1

End If

' Set uvw

triU = ( ( triX3 - lx ) * ( triY2 - sy ) - ( triY3 - sy ) * ( triX2 - lx ) ) * 0.5 * tri1OverArea

triV = ( ( triX1 - lx ) * ( triY3 - sy ) - ( triY1 - sy ) * ( triX3 - lx ) ) * 0.5 * tri1OverArea

triW = ( ( triX2 - lx ) * ( triY1 - sy ) - ( triY2 - sy ) * ( triX1 - lx ) ) * 0.5 * tri1OverArea

' Set duvw/dx

dudx = ( triY3 - triY2 ) * 0.5 * tri1OverArea

dvdx = ( triY1 - triY3 ) * 0.5 * tri1OverArea

dwdx = ( triY2 - triY1 ) * 0.5 * tri1OverArea

' Set duvw/dy

dudy = ( triX2 - triX3 ) * 0.5 * tri1OverArea

dvdy = ( triX3 - triX1 ) * 0.5 * tri1OverArea

dwdy = ( triX1 - triX2 ) * 0.5 * tri1OverArea

For pixel_y = sy To ey - 1

If pixel_y < 0 Or pixel_y >= Height Then

GoTo tri_flat_term_yskip

End If

For pixel_x = lx To rx

triU = triU + dudx

triV = triV + dvdx

triW = triW + dwdx

If pixel_x < 0 Or pixel_x >= Width Then

GoTo tri_flat_term_xskip

End If

pixel_index = pixel_y * Width + pixel_x

depth = 1 / ( triU * tri1OverZ1 + triV * tri1OverZ2 + triW * tri1OverZ3 )

If pixel_depth ( pixel_index ) > depth Or pixel_time ( pixel_index ) < current_frame Then

pixel_depth ( pixel_index ) = depth

pixel_time ( pixel_index ) = current_frame

Color _RGB ( triU * 256 , triV * 256 , triW * 256 )

PSet ( pixel_x , pixel_y )

End If

tri_flat_term_xskip:

Next

tri_flat_term_yskip:

triU = triU + dudy + lslope * dudx - ( CInt ( rx ) - CInt ( lx ) + 1 ) * dudx

triV = triV + dvdy + lslope * dvdx - ( CInt ( rx ) - CInt ( lx ) + 1 ) * dvdx

triW = triW + dwdy + lslope * dwdx - ( CInt ( rx ) - CInt ( lx ) + 1 ) * dwdx

lx = lx + lslope

rx = rx + rslope

Next

End Sub

Sub tri ( x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 )

'Project all 3 vertices to pixel-coordinates

setVector x1 , y1 , z1 , 1

projectVector

triX1 = vector.x: triY1 = vector.y: tri1OverZ1 = vector.z

setVector x2 , y2 , z2 , 1

projectVector

triX2 = vector.x: triY2 = vector.y: tri1OverZ2 = vector.z

setVector x3 , y3 , z3 , 1

projectVector

triX3 = vector.x: triY3 = vector.y: tri1OverZ3 = vector.z

' Compute extra triangle data

triArea = ( ( triX3 - triX1 ) * ( triY2 - triY1 ) - ( triY3 - triY1 ) * ( triX2 - triX1 ) ) * 0.5

tri1OverArea = 1 / triArea

' if any of the y coordinates are the same , draw a flat terminating triangle instead of splitting it

If triY1 = triY2 Then

tri_flat_term triX3 , triY3 , triX1 , triX2 , triY1: GoTo tri_skip

ElseIf triY2 = triY3 Then

tri_flat_term triX1 , triY1 , triX2 , triX3 , triY2: GoTo tri_skip

ElseIf triY3 = triY1 Then

tri_flat_term triX2 , triY2 , triX3 , triX1 , triY3: GoTo tri_skip

End If

Dim m

' Check which vertice is in the middle of the triangle in terms of y

If ( triY1 < triY2 And triY1 > triY3 ) Or ( triY1 > triY2 And triY1 < triY3 ) Then

m = triX2 + ( triY1 - triY2 ) * ( triX3 - triX2 ) / ( triY3 - triY2 )

tri_flat_term triX3 , triY3 , triX1 , m , triY1: tri_flat_term triX2 , triY2 , triX1 , m , triY1: GoTo tri_skip

ElseIf ( triY2 < triY1 And triY2 > triY3 ) Or ( triY2 > triY1 And triY2 < triY3 ) Then

m = triX1 + ( triY2 - triY1 ) * ( triX3 - triX1 ) / ( triY3 - triY1 )

tri_flat_term triX3 , triY3 , triX2 , m , triY2: tri_flat_term triX1 , triY1 , triX2 , m , triY2: GoTo tri_skip

Else

m = triX1 + ( triY3 - triY1 ) * ( triX2 - triX1 ) / ( triY2 - triY1 )

tri_flat_term triX2 , triY2 , triX3 , m , triY3: tri_flat_term triX1 , triY1 , triX3 , m , triY3

End If

tri_skip:

End Sub

Sub tri_vector ( a As vector , b As vector , c As vector )

tri a.x , a.y , a.z , b.x , b.y , b.z , c.x , c.y , c.z

End Sub

You might recall earlier in this paper when I provided images of cubes displaying textures. How exactly do we make this possible, and what more does that allow us to do? More often than not, barycentric coordinates are used to interpolate between texture coordinates. Unhelpfully denoted with U and V. (Sometimes even W.) This coordinate system is used to represent points on an image, typically spanning from 0 to 1.

If we give our barycentric coordinates texture coordinates to interpolate, we can take the resulting texture coordinate and use it to grab a color from an image. We can then draw this color directly to the screen, or we can do more processing with it.

A lot of modern graphics rely on clever usage of textures in order to create different lighting conditions and produce a pixel based on a combination of factors. And this is most commonly done through shaders.

A shader is a routine that is used to transform a set of input values to output values for the sake of 3d rendering. There are different types of shaders: vertex shaders which take in vertex attributes and tells the computer where it should lie on the screen or raster, geometry shaders which take in as much as whole 3d shapes and is able to transform existing and create new geometry, and finally fragment shaders, which take in anything given to them and output a color that is pushed to the raster. We have already done the work of the vertex shader, don't really need geometry shaders in our example, and pushing the barycentric coordinates to the screen is enough to constitute as a fragment shader, but lets push it further.

Fragment shaders can do as little as giving the same color to every pixel it's assigned to, or as complicated as simulating virtual light to immitate realistic lighting as best it can. All the abilities are covered quite in depth elsewhere, so I won't be digging deeper into lighting. However, take note that QB64 allows you to load textures and access their colors, and these colors can be used to manipulate normal vectors and control how strong lighting effects are. So take this last piece of BASIC code with you, and go learn more about the concepts of 3D rendering.

x As Double

y As Double

z As Double

w As Double

End Type

Type light

position As vector

intensity As Double

End Type

Dim Shared As Integer Width , Height

Dim Shared tri_data ( 29 ) As Double ' Store vertice position data , 1/z1-z3 , area , 1/area , uvw , normal , vertice phyiscal position

Dim Shared vector ( 4 ) As Double ' Store xyzw for a vector used for transformations

Dim Shared matrix ( 16 ) As Double ' Store a matrix that can be used to transform vector

Dim Shared projection_data ( 4 ) As Double ' Store FOV * 0.5 , FOV * Width * 0.5 , Zm , and Za

Dim Shared render_data ( 14 ) As Double 'left inverse slope , right inverse slope , left line x , right line x , starting y , ending y , uvw derivatives , index , depth

Dim Shared pixel_x As Integer

Dim Shared pixel_y As Integer

Dim Shared As vector shader_color , physical_position

Dim Shared current_frame As Integer

current_frame = 1

Dim Shared As Integer NUM_OF_ATTRIBUTES , NUM_OF_IMAGES , NUM_OF_LIGHTS

NUM_OF_ATTRIBUTES = 2

NUM_OF_IMAGES = 16

NUM_OF_LIGHTS = 16

Dim Shared verticeAttributes ( 3 * NUM_OF_ATTRIBUTES ) As Double ' Values that will be interpolated in a triangle

Dim Shared pixelAttributes ( NUM_OF_ATTRIBUTES ) As Double ' The resulting interpolated attributes

Dim Shared images& ( NUM_OF_IMAGES ) ' Shared images that can be used as textures and so forth , images& ( 0 ) is always screen

Dim Shared lights ( NUM_OF_LIGHTS ) As light

Width = 1200

Height = 800

Dim Shared pixel_depth ( Width * Height ) As Double

Dim Shared pixel_time ( Width * Height ) As Double

GoSub init_set_up_display

images& ( 1 ) = _LoadImage ( "metal\color.jpg" )

images& ( 2 ) = _LoadImage ( "metal\normal.jpg" )

images& ( 3 ) = _LoadImage ( "metal\metal.jpg" )

images& ( 4 ) = _LoadImage ( "metal\ao.jpg" )

Call setProjection ( 90 , 10000 , 0.01 )

'Do

'Loop

'verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 0: verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1: verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

'verticeAttributes ( 6 ) = 0: verticeAttributes ( 7 ) = 1: verticeAttributes ( 8 ) = 0: verticeAttributes ( 9 ) = 1: verticeAttributes ( 10 ) = 1: verticeAttributes ( 11 ) = 1

'verticeAttributes ( 12 ) = 0: verticeAttributes ( 13 ) = 0: verticeAttributes ( 14 ) = 1: verticeAttributes ( 15 ) = 1: verticeAttributes ( 16 ) = 1: verticeAttributes ( 17 ) = 0

lights ( 0 ) .intensity = 1

lights ( 0 ) .position.x = -100

lights ( 0 ) .position.y = -200

lights ( 0 ) .position.z = 50

lights ( 0 ) .position.w = 1

' Call tri ( -50 , -50 , -100 , 50 , 50 , -100 , 50 , -85 , -200 )

Dim As vector position , rotation , size , lightrot , lightsize

position.x = 0

position.y = 0

position.z = -200

position.w = 1

size.x = 50

size.y = 50

size.z = 50

size.w = 1

lightrot.x = 0

lightrot.y = 0

lightrot.z = 0

lightrot.w = 0

lightsize.x = 10

lightsize.y = 10

lightsize.z = 10

rotation.x = 0.7

rotation.y = 0.4

Do

Cls

'rotation.z = current_frame * 0.05

rotation.w = 1

If _KeyDown ( 119 ) Then

lights ( 0 ) .position.z = lights ( 0 ) .position.z + 10

End If

If _KeyDown ( 115 ) Then

lights ( 0 ) .position.z = lights ( 0 ) .position.z - 10

End If

If _KeyDown ( 97 ) Then

lights ( 0 ) .position.x = lights ( 0 ) .position.x + 10

End If

If _KeyDown ( 100 ) Then

lights ( 0 ) .position.x = lights ( 0 ) .position.x - 10

End If

If _KeyDown ( 101 ) Then

lights ( 0 ) .position.y = lights ( 0 ) .position.y + 10

End If

If _KeyDown ( 113 ) Then

lights ( 0 ) .position.y = lights ( 0 ) .position.y - 10

End If

If _KeyDown ( 19200 ) Then

rotation.y = rotation.y - 0.03

End If

If _KeyDown ( 19712 ) Then

rotation.y = rotation.y + 0.03

End If

If _KeyDown ( 18432 ) Then

rotation.x = rotation.x - 0.03

End If

If _KeyDown ( 20480 ) Then

rotation.x = rotation.x + 0.03

End If

Color _RGB ( 255 , 255 , 255 )

Print lights ( 0 ) .position.x; lights ( 0 ) .position.y; lights ( 0 ) .position.z

Call cube ( position , size , rotation )

'Color _RGB ( 255 , 255 , 255 )

'Print matrix ( 0 ) ; matrix ( 1 ) ; matrix ( 2 ) ; matrix ( 3 )

'Print matrix ( 4 ) ; matrix ( 5 ) ; matrix ( 6 ) ; matrix ( 7 )

'Print matrix ( 8 ) ; matrix ( 9 ) ; matrix ( 10 ) ; matrix ( 11 )

'Print matrix ( 12 ) ; matrix ( 13 ) ; matrix ( 14 ) ; matrix ( 15 )

'Call cube ( lights ( 0 ) .position , lightsize , lightrot )

current_frame = current_frame + 1

_Display

Loop

Color _RGB ( 255 , 255 , 255 )

End

init_set_up_display:

display = _NewImage ( Width , Height , 32 )

images& ( 0 ) = display

Screen display

Return

Sub cube ( position As vector , size As vector , rotation As vector )

Dim As vector lft , rft , lbt , rbt , lfb , rfb , lbb , rbb

setVector -size.x , size.y , -size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

lft.x = vector ( 0 ) + position.x: lft.y = vector ( 1 ) + position.y: lft.z = vector ( 2 ) + position.z

setVector size.x , size.y , -size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

rft.x = vector ( 0 ) + position.x: rft.y = vector ( 1 ) + position.y: rft.z = vector ( 2 ) + position.z

setVector -size.x , size.y , size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

lbt.x = vector ( 0 ) + position.x: lbt.y = vector ( 1 ) + position.y: lbt.z = vector ( 2 ) + position.z

setVector size.x , size.y , size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

rbt.x = vector ( 0 ) + position.x: rbt.y = vector ( 1 ) + position.y: rbt.z = vector ( 2 ) + position.z

setVector -size.x , -size.y , -size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

lfb.x = vector ( 0 ) + position.x: lfb.y = vector ( 1 ) + position.y: lfb.z = vector ( 2 ) + position.z

setVector size.x , -size.y , -size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

rfb.x = vector ( 0 ) + position.x: rfb.y = vector ( 1 ) + position.y: rfb.z = vector ( 2 ) + position.z

setVector -size.x , -size.y , size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

lbb.x = vector ( 0 ) + position.x: lbb.y = vector ( 1 ) + position.y: lbb.z = vector ( 2 ) + position.z

setVector size.x , -size.y , size.z , 0

rotateVector rotation.x , rotation.y , rotation.z

rbb.x = vector ( 0 ) + position.x: rbb.y = vector ( 1 ) + position.y: rbb.z = vector ( 2 ) + position.z

'Print lbt.x; lbt.y; lbt.z; lft.x; lft.y; lft.z; lfb.x; lfb.y; lfb.z , lbb.x; lbb.y; lbb.z; lfb.x; lfb.y; lfb.z; lbt.x; lbt.y; lbt.z

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector lft , rft , rfb

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector rfb , lfb , lft

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector rbt , lbt , lbb

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector lbb , rbb , rbt

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector lbt , lft , lfb

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector lfb , lbb , lbt

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector rft , rbt , rbb

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector rbb , rfb , rft

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector rft , lbt , rbt

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector lbt , rft , lft

verticeAttributes ( 0 ) = 0: verticeAttributes ( 1 ) = 0

verticeAttributes ( 2 ) = 1: verticeAttributes ( 3 ) = 0

verticeAttributes ( 4 ) = 1: verticeAttributes ( 5 ) = 1

tri_vector lfb , rbb , lbb

verticeAttributes ( 0 ) = 1: verticeAttributes ( 1 ) = 1

verticeAttributes ( 2 ) = 0: verticeAttributes ( 3 ) = 1

verticeAttributes ( 4 ) = 0: verticeAttributes ( 5 ) = 0

tri_vector rbb , lfb , rfb

End Sub

Sub shade

Dim As vector diffuseMap , normalMap , tolight , reflect , normal_pos

Dim As Single reflectMag , roughness , ao

Call texture2D ( images& ( 1 ) , pixelAttributes ( 0 ) , pixelAttributes ( 1 ) )

copyVectorTo diffuseMap

Call texture2D ( images& ( 2 ) , pixelAttributes ( 0 ) , pixelAttributes ( 1 ) )

vector ( 0 ) = vector ( 0 ) * 2 - 1: vector ( 1 ) = vector ( 1 ) * 2 - 1: vector ( 2 ) = vector ( 2 ) * 2 - 1

vector ( 3 ) = 0

applyMatrix

copyVectorTo normalMap

Call texture2D ( images& ( 3 ) , pixelAttributes ( 0 ) , pixelAttributes ( 1 ) )

roughness = vector ( 0 )

Call texture2D ( images& ( 4 ) , pixelAttributes ( 0 ) , pixelAttributes ( 1 ) )

ao = vector ( 0 )

' Apply normal map to normals

'initVector shader_color , normalMap.x , normalMap.y , normalMap.z , 1

tolight.x = lights ( 0 ) .position.x - physical_position.x

tolight.y = lights ( 0 ) .position.y - physical_position.y

tolight.z = lights ( 0 ) .position.z - physical_position.z

normalizeVector tolight

diffuse = max ( tolight.x * normalMap.x + tolight.y * normalMap.y + tolight.z * normalMap.z , 0 ) * 0.5

reflectMag = 2.0 * ( -tolight.x * normalMap.x - tolight.y * normalMap.y - tolight.z * normalMap.z )

reflect.x = -tolight.x - reflectMag * normalMap.x: reflect.y = -tolight.y - reflectMag * normalMap.y: reflect.z = -tolight.z - reflectMag * normalMap.z

normalizeVector reflect

normal_pos.x = physical_position.x: normal_pos.y = physical_position.y: normal_pos.z = physical_position.z

normalizeVector normal_pos

specular = max ( -normal_pos.x * reflect.x - normal_pos.y * reflect.y - normal_pos.z * reflect.z , 0 ) ^ 32

diffuseMap.x = diffuseMap.x * ( diffuse + specular * roughness + 0.2 ) * ao

diffuseMap.y = diffuseMap.y * ( diffuse + specular * roughness + 0.2 ) * ao

diffuseMap.z = diffuseMap.z * ( diffuse + specular * roughness + 0.2 ) * ao

'initVector shader_color , 0.1 + diffuse + specular , 0.1 + diffuse + specular , 0.1 + diffuse + specular , 1

initVector shader_color , diffuseMap.x , diffuseMap.y , diffuseMap.z , 1

'initVector shader_color , specular , specular , specular , 1

'initVector shader_color , reflect.x , reflect.y , reflect.z , 1

' Print vector ( 2 )

End Sub

' Returns the max between two numbers

Function max ( ma_x , ma_y )

max = ma_x

If ma_y > ma_x Then

max = ma_y

End If

End Function

' Returns the min between two numbers

Function min ( mi_x , mi_y )

min = mi_x

If mi_y < mi_x Then

min = mi_y

End If

End Function

Function cap ( in As Double , minimum As Double , maximum As Double )

cap = min ( max ( in , minimum ) , maximum )

End Function

Sub applyMatrix

vector ( 0 ) = vector ( 0 ) * matrix ( 0 ) + vector ( 1 ) * matrix ( 1 ) + vector ( 2 ) * matrix ( 2 ) + vector ( 3 ) * matrix ( 3 )

vector ( 1 ) = vector ( 0 ) * matrix ( 4 ) + vector ( 1 ) * matrix ( 5 ) + vector ( 2 ) * matrix ( 6 ) + vector ( 3 ) * matrix ( 7 )

vector ( 2 ) = vector ( 0 ) * matrix ( 8 ) + vector ( 1 ) * matrix ( 9 ) + vector ( 2 ) * matrix ( 10 ) + vector ( 3 ) * matrix ( 11 )

vector ( 3 ) = vector ( 0 ) * matrix ( 12 ) + vector ( 1 ) * matrix ( 13 ) + vector ( 2 ) * matrix ( 14 ) + vector ( 3 ) * matrix ( 15 )

End Sub

Sub setVector ( x As Double , y As Double , z As Double , w As Double )

vector ( 0 ) = x

vector ( 1 ) = y

vector ( 2 ) = z

vector ( 3 ) = w

End Sub

Sub initVector ( invector As vector , x As Double , y As Double , z As Double , w As Double )

invector.x = x

invector.y = y

invector.z = z

invector.w = w

End Sub

Sub setProjection ( FOVAngle As Double , Far As Double , Near As Double )

projection_data ( 0 ) = 0.5 / Tan ( FOVAngle * 3.14159267 / 360.0 )

projection_data ( 1 ) = Width * projection_data ( 0 )

projection_data ( 2 ) = -Far / ( Far - Near )

projection_data ( 3 ) = - ( Far * Near ) / ( Far - Near )

End Sub

Sub projectVector ' Leaves z as 1/z

vector ( 2 ) = 1 / ( vector ( 2 ) * projection_data ( 2 ) + vector ( 3 ) * projection_data ( 3 ) )

vector ( 0 ) = vector ( 0 ) * vector ( 2 ) * projection_data ( 1 ) + vector ( 3 ) * 0.5 * Width

vector ( 1 ) = vector ( 1 ) * vector ( 2 ) * projection_data ( 1 ) + vector ( 3 ) * 0.5 * Height

End Sub

Sub normalizeVector ( invector As vector )

mag = Sqr ( invector.x ^ 2 + invector.y ^ 2 + invector.z ^ 2 )

invector.x = invector.x / mag: invector.y = invector.y / mag: invector.z = invector.z / mag

End Sub

Sub rotateVector ( rx As Double , ry As Double , rz As Double )

Dim nvector As vector

Call initVector ( nvector , 0 , 0 , 0 , 0 )

nvector.x = ( vector ( 0 ) * Cos ( rz ) * Cos ( ry ) ) + ( vector ( 1 ) * Sin ( rz ) * Cos ( ry ) ) - vector ( 2 ) * Sin ( ry )

nvector.y = ( vector ( 0 ) * ( Cos ( rz ) * Sin ( ry ) * Sin ( rx ) - Sin ( rz ) * Cos ( rx ) ) ) + ( vector ( 1 ) * ( Sin ( rz ) * Sin ( ry ) * Sin ( rx ) + Cos ( rz ) * Cos ( rx ) ) ) + ( vector ( 2 ) * Cos ( ry ) * Sin ( rx ) )

nvector.z = ( vector ( 0 ) * ( Cos ( rz ) * Sin ( ry ) * Cos ( rx ) + Sin ( rz ) * Sin ( rx ) ) ) + ( vector ( 1 ) * ( Sin ( rz ) * Sin ( ry ) * Cos ( rx ) - Cos ( rz ) * Sin ( rx ) ) ) + ( vector ( 2 ) * Cos ( ry ) * Cos ( rx ) )

nvector.w = 1

copyToVector nvector

End Sub

Sub copyVectorTo ( invector As vector )

invector.x = vector ( 0 )

invector.y = vector ( 1 )

invector.z = vector ( 2 )

invector.w = vector ( 3 )

End Sub

Sub copyToVector ( outvector As vector )

vector ( 0 ) = outvector.x

vector ( 1 ) = outvector.y

vector ( 2 ) = outvector.z

vector ( 3 ) = outvector.w

End Sub

' Gravs color

Sub texture2D ( image& , u As Double , v As Double )

_Source image&

GRABBED_COLOR~& = Point ( Fix ( cap ( u , 0 , 1 ) * ( _Width ( image& ) - 1 ) ) , Fix ( cap ( v , 0 , 1 ) * ( _Height ( image& ) - 1 ) ) )

vector ( 0 ) = ( Fix ( GRABBED_COLOR~& / 65536 ) Mod 256 ) / 256.0

vector ( 1 ) = ( Fix ( GRABBED_COLOR~& / 256 ) Mod 256 ) / 256.0

vector ( 2 ) = ( Fix ( GRABBED_COLOR~& ) Mod 256 ) / 256.0

vector ( 3 ) = ( Fix ( GRABBED_COLOR~& / 16777216 ) Mod 256 ) / 256.0

End Sub

Sub tri_flat_term ( x1 , y1 , x2 , x3 , term )

'render_data ( 5 ) = Sgn ( y1 - term )

render_data ( 2 ) = x2 ' Set left x that will be walking line ( x2 , term ) - ( x1 , y1 )

render_data ( 3 ) = x3 ' Set right x that will be walking line ( x3 , term ) - ( x1 , y1 )

If ( x2 > x3 ) Then

render_data ( 3 ) = x2: render_data ( 2 ) = x3

End If

render_data ( 0 ) = ( render_data ( 2 ) - x1 ) / ( term - y1 ) ' Negative left inverse slope depending on if y will be increasing or decreasing

render_data ( 1 ) = ( render_data ( 3 ) - x1 ) / ( term - y1 ) ' Negative right inverse slope for same reason

render_data ( 4 ) = term

render_data ( 5 ) = y1

If ( y1 < term ) Then

render_data ( 4 ) = y1: render_data ( 5 ) = term: render_data ( 2 ) = x1: render_data ( 3 ) = x1

End If

'render_data ( 4 ) = Sgn ( render_data ( 1 ) - render_data ( 0 ) ) ' Which direction x will need to increment in , + or -

' Set uvw

tri_data ( 14 ) = ( ( tri_data ( 6 ) - render_data ( 2 ) ) * ( tri_data ( 4 ) - render_data ( 4 ) ) - ( tri_data ( 7 ) - render_data ( 4 ) ) * ( tri_data ( 3 ) - render_data ( 2 ) ) ) * 0.5 * tri_data ( 13 ) 'u = ( ( x3-starting x ) * ( y2-starting y ) - ( y3-starting y ) * ( x2-starting x ) ) * 0.5 * 1/tri area

tri_data ( 15 ) = ( ( tri_data ( 0 ) - render_data ( 2 ) ) * ( tri_data ( 7 ) - render_data ( 4 ) ) - ( tri_data ( 1 ) - render_data ( 4 ) ) * ( tri_data ( 6 ) - render_data ( 2 ) ) ) * 0.5 * tri_data ( 13 ) 'v = ( ( x1-starting x ) * ( y3-starting y ) - ( y1-starting y ) * ( x3-starting x ) ) * 0.5 * 1/tri area

tri_data ( 16 ) = ( ( tri_data ( 3 ) - render_data ( 2 ) ) * ( tri_data ( 1 ) - render_data ( 4 ) ) - ( tri_data ( 4 ) - render_data ( 4 ) ) * ( tri_data ( 0 ) - render_data ( 2 ) ) ) * 0.5 * tri_data ( 13 ) 'w = ( ( x2-starting x ) * ( y1-starting y ) - ( y2-starting y ) * ( x1-starting x ) ) * 0.5 * 1/tri area

' Set duvw/dx

render_data ( 6 ) = ( tri_data ( 7 ) - tri_data ( 4 ) ) * 0.5 * tri_data ( 13 ) ' du/dx = y3 - y2

render_data ( 7 ) = ( tri_data ( 1 ) - tri_data ( 7 ) ) * 0.5 * tri_data ( 13 ) ' du/dx = y1 - y3

render_data ( 8 ) = ( tri_data ( 4 ) - tri_data ( 1 ) ) * 0.5 * tri_data ( 13 ) ' du/dx = y2 - y1

' Set duvw/dy

render_data ( 9 ) = ( tri_data ( 3 ) - tri_data ( 6 ) ) * 0.5 * tri_data ( 13 ) ' du/dy = x2 - x3

render_data ( 10 ) = ( tri_data ( 6 ) - tri_data ( 0 ) ) * 0.5 * tri_data ( 13 ) ' du/dy = x3 - x1

render_data ( 11 ) = ( tri_data ( 0 ) - tri_data ( 3 ) ) * 0.5 * tri_data ( 13 ) ' du/dy = x1 - x2

'Print render_data ( 9 ) ; render_data ( 6 ) ; tri_data ( 14 ) , render_data ( 10 ) ; render_data ( 7 ) ; tri_data ( 15 ) , render_data ( 11 ) ; render_data ( 8 ) ; tri_data ( 16 )

For pixel_y = render_data ( 4 ) To render_data ( 5 ) - 1

uvwxinctimes = 0

If pixel_y < 0 Or pixel_y >= Height Then

GoTo tri_flat_term_yskip

End If

For pixel_x = render_data ( 2 ) To render_data ( 3 )

'If uvwxinctimes = 0 Then Print pixel_x; render_data ( 2 )

tri_data ( 14 ) = tri_data ( 14 ) + render_data ( 6 )

tri_data ( 15 ) = tri_data ( 15 ) + render_data ( 7 )

tri_data ( 16 ) = tri_data ( 16 ) + render_data ( 8 )

If pixel_x < 0 Or pixel_x >= Width Then

GoTo tri_flat_term_xskip

End If

render_data ( 12 ) = pixel_y * Width + pixel_x

render_data ( 13 ) = 1 / ( tri_data ( 14 ) * tri_data ( 9 ) + tri_data ( 15 ) * tri_data ( 10 ) + tri_data ( 16 ) * tri_data ( 11 ) )

If pixel_depth ( render_data ( 12 ) ) > render_data ( 13 ) Or pixel_time ( render_data ( 12 ) ) < current_frame Then

pixel_depth ( render_data ( 12 ) ) = render_data ( 13 )

pixel_time ( render_data ( 12 ) ) = current_frame

For i = 0 To ( NUM_OF_ATTRIBUTES - 1 )

pixelAttributes ( i ) = verticeAttributes ( 0 * NUM_OF_ATTRIBUTES + i ) * tri_data ( 14 ) * tri_data ( 9 ) * render_data ( 13 ) + verticeAttributes ( 1 * NUM_OF_ATTRIBUTES + i ) * tri_data ( 15 ) * tri_data ( 10 ) * render_data ( 13 ) + verticeAttributes ( 2 * NUM_OF_ATTRIBUTES + i ) * tri_data ( 16 ) * tri_data ( 11 ) * render_data ( 13 )

Next

physical_position.x = tri_data ( 20 ) * tri_data ( 14 ) * tri_data ( 9 ) * render_data ( 13 ) + tri_data ( 23 ) * tri_data ( 15 ) * tri_data ( 10 ) * render_data ( 13 ) + tri_data ( 26 ) * tri_data ( 16 ) * tri_data ( 11 ) * render_data ( 13 )

physical_position.y = tri_data ( 21 ) * tri_data ( 14 ) * tri_data ( 9 ) * render_data ( 13 ) + tri_data ( 24 ) * tri_data ( 15 ) * tri_data ( 10 ) * render_data ( 13 ) + tri_data ( 27 ) * tri_data ( 16 ) * tri_data ( 11 ) * render_data ( 13 )

physical_position.z = -render_data ( 13 )

physical_position.w = 1

shade

Color _RGBA ( shader_color.x * 256 , shader_color.y * 256 , shader_color.z * 256 , shader_color.w * 256 )

PSet ( pixel_x , pixel_y )

Else

'Print "Depth bypass"; pixel_x; pixel_y; render_data ( 13 )

End If

tri_flat_term_xskip:

Next

tri_flat_term_yskip:

tri_data ( 14 ) = tri_data ( 14 ) + render_data ( 9 ) + render_data ( 0 ) * render_data ( 6 ) - ( CInt ( render_data ( 3 ) ) - CInt ( render_data ( 2 ) ) + 1 ) * render_data ( 6 )

tri_data ( 15 ) = tri_data ( 15 ) + render_data ( 10 ) + render_data ( 0 ) * render_data ( 7 ) - ( CInt ( render_data ( 3 ) ) - CInt ( render_data ( 2 ) ) + 1 ) * render_data ( 7 )

tri_data ( 16 ) = tri_data ( 16 ) + render_data ( 11 ) + render_data ( 0 ) * render_data ( 8 ) - ( CInt ( render_data ( 3 ) ) - CInt ( render_data ( 2 ) ) + 1 ) * render_data ( 8 )

render_data ( 2 ) = render_data ( 2 ) + render_data ( 0 ) 'left x = left x + left inverse slope

render_data ( 3 ) = render_data ( 3 ) + render_data ( 1 ) 'right x = right x + right inverse slope

Next

End Sub

Sub tri ( x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 )

'Project all 3 vertices to pixel-coordinates

tri_data ( 20 ) = x1: tri_data ( 21 ) = y1: tri_data ( 22 ) = z1

tri_data ( 23 ) = x2: tri_data ( 24 ) = y2: tri_data ( 25 ) = z2

tri_data ( 26 ) = x3: tri_data ( 27 ) = y3: tri_data ( 28 ) = z3

setVector x1 , y1 , z1 , 1

projectVector

tri_data ( 0 ) = vector ( 0 ) : tri_data ( 1 ) = vector ( 1 ) : tri_data ( 9 ) = vector ( 2 )

setVector x2 , y2 , z2 , 1

projectVector

tri_data ( 3 ) = vector ( 0 ) : tri_data ( 4 ) = vector ( 1 ) : tri_data ( 10 ) = vector ( 2 )

setVector x3 , y3 , z3 , 1

projectVector

tri_data ( 6 ) = vector ( 0 ) : tri_data ( 7 ) = vector ( 1 ) : tri_data ( 11 ) = vector ( 2 )

' Compute extra triangle data

tri_data ( 12 ) = ( ( tri_data ( 6 ) - tri_data ( 0 ) ) * ( tri_data ( 4 ) - tri_data ( 1 ) ) - ( tri_data ( 7 ) - tri_data ( 1 ) ) * ( tri_data ( 3 ) - tri_data ( 0 ) ) ) * 0.5 'Area = ( ( x3 - x1 ) * ( y2 - y1 ) - ( y3 - y1 ) * ( x2 - x1 ) ) * 0.5

tri_data ( 13 ) = 1 / tri_data ( 12 ) '1/Area

' if any of the y coordinates are the same , draw a flat terminating triangle instead of splitting it

' if y1=y2 then tri_flat_term x3 , y3 , x1 , x2 , terminating at y1

' if y2=y3 then tri_flat_term x1 , y1 , x2 , x3 , terminating at y2

' if y3=y1 then tri_flat_term x2 , y2 , x3 , x1 , terminating at y3

If tri_data ( 1 ) = tri_data ( 4 ) Then

tri_flat_term tri_data ( 6 ) , tri_data ( 7 ) , tri_data ( 0 ) , tri_data ( 3 ) , tri_data ( 1 ) : GoTo tri_skip

ElseIf tri_data ( 4 ) = tri_data ( 7 ) Then

tri_flat_term tri_data ( 0 ) , tri_data ( 1 ) , tri_data ( 3 ) , tri_data ( 6 ) , tri_data ( 4 ) : GoTo tri_skip

ElseIf tri_data ( 7 ) = tri_data ( 1 ) Then

tri_flat_term tri_data ( 3 ) , tri_data ( 4 ) , tri_data ( 6 ) , tri_data ( 0 ) , tri_data ( 7 ) : GoTo tri_skip

End If

Dim m

' Check which vertice is in the middle of the triangle in terms of y

If ( tri_data ( 1 ) < tri_data ( 4 ) And tri_data ( 1 ) > tri_data ( 7 ) ) Or ( tri_data ( 1 ) > tri_data ( 4 ) And tri_data ( 1 ) < tri_data ( 7 ) ) Then

m = tri_data ( 3 ) + ( tri_data ( 1 ) - tri_data ( 4 ) ) * ( tri_data ( 6 ) - tri_data ( 3 ) ) / ( tri_data ( 7 ) - tri_data ( 4 ) ) 'x2+ ( y1-y2 ) * ( x3-x2 ) / ( y3-y2 )

' triangle ( x3 , y3 , x1 , m , y1 ) + triangle ( x2 , y2 , x1 , m , y1 )

tri_flat_term tri_data ( 6 ) , tri_data ( 7 ) , tri_data ( 0 ) , m , tri_data ( 1 ) : tri_flat_term tri_data ( 3 ) , tri_data ( 4 ) , tri_data ( 0 ) , m , tri_data ( 1 ) : GoTo tri_skip

ElseIf ( tri_data ( 4 ) < tri_data ( 1 ) And tri_data ( 4 ) > tri_data ( 7 ) ) Or ( tri_data ( 4 ) > tri_data ( 1 ) And tri_data ( 4 ) < tri_data ( 7 ) ) Then

m = tri_data ( 0 ) + ( tri_data ( 4 ) - tri_data ( 1 ) ) * ( tri_data ( 6 ) - tri_data ( 0 ) ) / ( tri_data ( 7 ) - tri_data ( 1 ) ) 'x1+ ( y2-y1 ) * ( x3-x1 ) / ( y3-y1 )

' triangle ( x3 , y3 , x2 , m , y2 ) + triangle ( x1 , y1 , x2 , m , y2 )

tri_flat_term tri_data ( 6 ) , tri_data ( 7 ) , tri_data ( 3 ) , m , tri_data ( 4 ) : tri_flat_term tri_data ( 0 ) , tri_data ( 1 ) , tri_data ( 3 ) , m , tri_data ( 4 ) : GoTo tri_skip

Else

m = tri_data ( 0 ) + ( tri_data ( 7 ) - tri_data ( 1 ) ) * ( tri_data ( 3 ) - tri_data ( 0 ) ) / ( tri_data ( 4 ) - tri_data ( 1 ) ) 'x1+ ( y3-y1 ) * ( x2-x1 ) / ( y2-y1 )

' triangle ( x2 , y2 , x3 , m , y3 ) + triangle ( x1 , y1 , x3 , m , y3 )

tri_flat_term tri_data ( 3 ) , tri_data ( 4 ) , tri_data ( 6 ) , m , tri_data ( 7 ) : tri_flat_term tri_data ( 0 ) , tri_data ( 1 ) , tri_data ( 6 ) , m , tri_data ( 7 )

End If

tri_skip:

End Sub

Sub tri_vector ( a As vector , b As vector , c As vector )

tri_data ( 17 ) = ( ( b.y - a.y ) * ( c.z - a.z ) - ( b.z - a.z ) * ( c.y - a.y ) )

tri_data ( 18 ) = ( ( b.z - a.z ) * ( c.x - a.x ) - ( b.x - a.x ) * ( c.z - a.z ) )

tri_data ( 19 ) = ( ( b.x - a.x ) * ( c.y - a.y ) - ( b.y - a.y ) * ( c.x - a.x ) )

mag = Sqr ( tri_data ( 17 ) ^ 2 + tri_data ( 18 ) ^ 2 + tri_data ( 19 ) ^ 2 )

tri_data ( 17 ) = tri_data ( 17 ) / mag

tri_data ( 18 ) = tri_data ( 18 ) / mag

tri_data ( 19 ) = tri_data ( 19 ) / mag

deltaUV1x = verticeAttributes ( 2 ) - verticeAttributes ( 0 ) : deltaUV1y = verticeAttributes ( 3 ) - verticeAttributes ( 1 )

deltaUV2x = verticeAttributes ( 4 ) - verticeAttributes ( 0 ) : deltaUV2y = verticeAttributes ( 5 ) - verticeAttributes ( 1 )

f = 1.0 / ( deltaUV1x * deltaUV2y - deltaUV2x * deltaUV1y )

Dim As vector tang , bitan

tang.x = f * ( deltaUV2y * ( b.x - a.x ) - deltaUV1y * ( c.x - a.x ) )

tang.y = f * ( deltaUV2y * ( b.y - a.y ) - deltaUV1y * ( c.y - a.y ) )

tang.z = f * ( deltaUV2y * ( b.z - a.z ) - deltaUV1y * ( c.z - a.z ) )

normalizeVector tang

bitan.x = f * ( -deltaUV2x * ( b.x - a.x ) + deltaUV1x * ( c.x - a.x ) )

bitan.y = f * ( -deltaUV2x * ( b.y - a.y ) + deltaUV1x * ( c.y - a.y ) )

bitan.z = f * ( -deltaUV2x * ( b.z - a.z ) + deltaUV1x * ( c.z - a.z ) )

normalizeVector bitan

matrix ( 0 ) = tang.x: matrix ( 1 ) = tang.y: matrix ( 2 ) = tang.z

matrix ( 4 ) = bitan.x: matrix ( 5 ) = bitan.y: matrix ( 6 ) = bitan.z

matrix ( 8 ) = tri_data ( 17 ) : matrix ( 9 ) = tri_data ( 18 ) : matrix ( 10 ) = tri_data ( 19 )

tri a.x , a.y , a.z , b.x , b.y , b.z , c.x , c.y , c.z

End Sub

3blue1brown

A very visual approach to linear algebra

For more on the math behind 3D Computer graphics:

Scratchapixel

An explanation of perspective matrices

An explanation of correct depth interpolation

An explanation of Barycentric coordinate correction based on depth

For more on shading, from simple to complex 3D lighting:

LearnOpenGL

While centered on OpenGL, LearnOpenGL provides a lot more details behind plenty of shading techniques

Basic lighting

Lighting maps

Normal Mapping