2. 3D Geometry

2.1 The Context Between 2D and 3D Worlds

A point in 3-space is defined by three coordinates: x, y, and z. We can't draw a point like this straight to the computer screen. Hence, we need a formula with which we can make the x and y coordinates somehow dependent on the z coordinate. Let's think about the case. A point can be multiplied by a constant the result point staying on a line that goes through the origin and the original point:

k*(X, Y, Z) = (kX, kY, kZ)

If we define the origin as the location of the camera, we can use this formula and take 1/z as k. Now the z coordinate is a constant:

1/Z * (X, Y, Z) = (X/Z, Y/Z, 1).

This kind of points lie on the plane z=1. If we want our plane to be closer / farther, we can change the formula a bit:

a/Z * (X, Y, Z) = (X*a/Z, Y*a/Z, a)

(Using a actually affects the perspective.) Now the equation of the projection plane is z=a,

 
so all points in the 3-space are transformed onto this plane in a way that they stay on a line that goes through the origin and the original point. The formula is thus the following:
X_SCREEN = X0 * SCALE / Z0
Y_SCREEN = Y0 * SCALE / Z0
SCALE tells the perspective value of the world. Usually values like 256 are used but try yourself and find the best value suitable for your purposes! When the z value is changed, we notice a 3D effect.

Note! If z reaches somehow the value zero, the program naturally crashes (division by zero). Be careful!

Example: Visualizing three dimensions with a pseudo code (plots a pixel to the center of the screen. When the z value changes, it vanishes from the screen):

x = 1
y = 1
z = 256
while (gx<320) and (gx>-1) and (gy<200) and (gy>-1)
gx = x * 256 / z + 160 ; center = the center of the screen
gy = y * 256 / z + 100 ; (mode 320x200)
putpixel(gx,gy,15)
z = z - 1
if (z = 0)
z = -1 ; just to ensure...
endif
endwhile
It's very easy to improve the code above to make it a 3D starfield, a basic and widely used effect in the demoscene half a decade ago. A good excercise  btw!

2.2 The Matrix Technique -- Foreword

NOTE! Even if you didn't understand a word about matrices, see the source code, and use some cut'n'paste to create a base of your own 3D engine. Later when you feel ready, come back and learn more about 3D geometry; using matrix technique doesn't demand too much knowledge.
In the matrix technique, the orientation of an object is presented with an object matrix, which consists of three vectors (the directions down, left, and forward). Every operation is performed to the matrix, and not until in a certain point the vertices are taken into account by multiplying them by the object matrix.
This technique has many advantages: it's faster, easier to understand, and more exact than the 'traditional' way we usually see in 3D tutorials. In what way more exact, you ask.

Let's take an airplane as an example (yeah, straight from otmmatx.doc :) The original orientation of the plane is as follows: the nose points to the direction of the z-axis, the right wing to the direction of the x-axis, and the y-axis pointing to the 'up' direction of the plane. Rotate the plane about the y-axis so that the nose points to the direction of the negative x-axis. Now when you rotate it about z-axis, does it turn around its own z-axis or the space's z-axis? It should turn around its own axis (or how could you implement a flight simulator otherwise?), but with basic rotations this is not what happens. So matrices are used.

In 3D coding, we use 3x3 matrices (or 4x4, more later) in which the data is located as follows:
 


 
In English, the first row tells the coefficients of the object's x-axis vector, second and third rows y- and z-axis, respectively. Note! The vectors are unit vectors; hence their length should always be one!

It would maybe be better to format the object matrix always so that the object lies (as in the plane example) on the xz plane, so it is naturally a unit matrix (see 1.2.1). Now the unit vectors of the object's axes are i, j, and k, the same as the vectors of the world axes.

(Nearly) all 3D operations in matrix technique are done using matrix math.

2.3 Rotating the Object Matrix

When rotating an object, we just multiply the object matrix by these matrices (or the general rotation matrix derived below). Rotations in 3D can also be thought of as changes in the coordinate systems.

About the x-axis:
 


 
About the y-axis:
 

 
About the z-axis:
 

 
cx, sx, cy, sy, cz, and sz mean the cosines and sines of the rotation angles.

We can of course precalculate a bit (just basic matrix multiplying):

So we get the rotated object matrix O from the formula
 
O = O*XYZ.
NOTE! In the versions <2.2 of 3dica I recommended using precalculated trigonometric functions to speed up the engine. The Pentium performs trigonometric functions so fast there's no sense in consuming memory to get a large enough table. For example, when using big hierarchical objects, we'd need enormous trigonometric tables to get a good enough result. This is due to cumulativity and the possibly huge size of the object. Thanks to the #coders gurus who enlightened me in this topic.

2.3.1 Rotating about an arbitrary vector

Where does one need this formula? Well, you need it only if you want the airplane in the example or any object to rotate right in your engine :) Also you need this if you want to code keyframing. So the basic formula is this:

U = n*n(t) + cos(a)*(I-n*n(t)) + sin(a)*N(x).

n = the rotation axis vector
n(t) = the transpose of n
a = rotation angle
I = the unit matrix
N(x) = so-called cross product matrix:
,
             where n* are the elements of the rotation axis vector.

I myself derived the matrix from the formula above and got the following matrix (works, tested well):
 

 
If you try this matrix using the vectors (1,0,0), (0,1,0), and (0,0,1) as rotation axes one after another, you'll notice you get the right results (the rotation matrices around the x-, y-, and z-axis). And how to use this? Easy:
 
for a=0 -> 2
axis.x = objectmatrix[a][0] ; in (y,x) syntax
axis.y = objectmatrix[a][1]
axis.z = objectmatrix[a][2]
tm[a] = rotate_with_the_matrix_above(
axis,angle_x,angle_y,angle_z)
endfor
matrix_mul(tm[0],tm[1])
matrix_mul(tm[0],tm[2])
matrix_mul(objectmatrix,tm[0])
This is certainly not a too fast way. It works perfectly, though, but if you have any better suggestions please e-mail me.

Note! Remember to normalize the object matrices after some frames! The formula requires so much float math the object matrices tend to 'bend' very quickly. Using doubles, I myself need to normalize them about every tenth frame.

2.4 The Camera

Wow, it's so easy to attach a real camera into a matrix based engine: nothing else is needed than multiplying the rotated object matrix by the camera matrix! In practice:
 
O = O * XYZ
O' = O * C,
where O' is the rotated object matrix, and C camera matrix.

Note! Remember to transform vertices with the O' matrix and normals with the O matrix (otherwise the light sources are rotated also when the camera is rotated).

There are perfect :) camera routines (including camera moving) in the source code, go check 'em out!

2.4.1 Deriving the camera matrix from a vector

Chem asked about the use of this chapter. Just to let everyone to know, cameras in 3D Studio are defined only by a direction vector, so it's kinda hard to operate with them without certain mathematical operations and the information of this chapter.

Ok. This has been The Big Question for many 3D coders. I was also interested in the topic so much, I derived a 'nice' complicated and hard-to-use way of doing the thing on a maths class in about three hours. Right after I had finished, the guy sitting next to me, Jussi Vainionpaa, asked innocently, "wouldn't this be a bit more useful way?" I went through the system he had derived in five minutes, and came to the decision that not only his way was much more elegant, but it also worked better X) So here's the technique developed by Jussi Vainionpaa. Based on straightforward vector algebra.

Deriving a rotation matrix for the camera from a single vector to whose direction the camera is pointing (z = camera forward, y = camera down, j = world down, (j dot z)*z = the projection of j on z):
 


 
The world's down vector j is known: (0,1,0). We calculate the component of j parallel to z, (j dot z)*z, and substract it from j, the result being y:

y = j - (j dot z)*z.

(vector projection, see 1.1.6). The formula is very straightforward and can be derived straight from the situation in the picture. It even shrinks due to the speciality of j:
 
(yx,yy,yz)
= (0,1,0) - (0*zx+1*zy+0*zz)*(zx,zy,zz)
= (0,1,0) - (zy*zx,zy*zy,zy*zz)
= (-zy*zx,1-zy*zy,-zy*zz).
After normalization, what we have is the final y vector of the camera (vector normalization, see 1.1.3). The z vector is known, but it must also be normalized. Finally we get the x vector as the cross product of y and z (notice the order):

x = yxz.

The x vector needen't to be normalized if y and z already were unit vectors. Now we place the vectors to the camera matrix as in the chapter 2.2 and the camera matrix is ready.

Note that if the camera z-axis is equal to the y-axis of the world (or its opposite), you get (0,0,0) as the camera y-axis vector. This isn't of course too nice. Try performing a compare and if needed, using the world z- or x-axis instead. Thanks to Arho Huttunen who noticed this first.

2.4.2 B-Splines

They sound cool, they look cool, they're cool to use, they simply are cool. What in the world are they? They are B-splines!

B-splines are actually not lines but curves (ok, clever). They are used to find a nice way through or near a specified list of key points for example to smoothen camera movement. A B-spline looks like this:
 


 
As you probably noticed, this particular B-spline was created using a formula that doesn't require the spline to go through all of the key points. I'm describing here only this one.

A B-spline is constructed using arcs. If we want to form a B-spline using n arcs, we need n+3 key points, let them be c(0),c(1),...,c(n+2). The k:th arc (k = 1,2,...,n) of the B-spline is calculated from the formula
 

rk(t) =
1/6*(1-3t+3t^2-t^3)*c(k-1) +
1/6*(4-6t^2+3t^3)*c(k) +
1/6*(1+3t+3t^2-3t^3)*c(k+1) +
1/6*t^3*c(k+2),
where t ranges linearly from 0 to 1. The more different t values, the more accurate curve. Note that the ending point of the k:th arc and the starting point of the k+1:th arc are the same point. Even better, the curve is all smooth in these most significant points, too, mainly due to the truth that the shape of an arc is determined by as many as four key points summarized (see the formula above).

A piece of pseudo that I used to create the picture above:
 

- point is the key point table
- STEPS is the number of steps
- k is the current arc index
- bl_i = bspline index (zero at first)
for k=1 -> POINTS-2 ; note! start from 1, not 0
t=0
while t<1.0
t2=t*t; ; some precalculation
t3=t*t*t; ; more precalculation :)
k1=1-3*t+3*t2-t3;
k2=4-6*t2+3*t3;
k3=1+3*t+3*t2-3*t3;
bspline[bl_i].x =
1/6.0*(k1*point[k-1].x+
k2*point[k].x+
k3*point[k+1].x+
t3*point[k+2].x);
bspline[bl_i].y =
1/6.0*(k1*point[k-1].y+
k2*point[k].y+
k3*point[k+1].y+
t3*point[k+2].y);
bspline[bl_i].z =
1/6.0*(k1*point[k-1].z+
k2*point[k].z+
k3*point[k+1].z+
t3*point[k+2].z);
bl_i++;
t=t+1.0/STEPS
endwhile
endfor

2.5 Transforming a Vertex by the Object Matrix

Now we need only to transform all the vertices by the rotation matrix, and we're finished! Transforming can easily be performed by using the information in the chapter 1.2.2.3.1, regarding a vertex as a vector:
 
X = X0*a + Y0*d + Z0*g + center_x
Y = X0*b + Y0*e + Z0*h + center_y
Z = X0*c + Y0*f + Z0*i + center_z.
We need thus nine adds and nine multiplies.

Pseudo example (rotating a point around all the three axes with the matrix technique):
 

• tm (transformation matrix), om (object matrix), t_om (transformed object matrix) (3,3)-sized float tables
• (xp,yp,zp) original point, (x,y,z) rotated point
• ox, oy, oz object origo coordinates
• ORIGO_X, ORIGO_Y screen center coordinates
reset_matrix(om) ; object matrix is reset only at the beginning
xa,ya,za = 1*PI/180 ; one degree in radians (rotation angles)
ox = 0 ; object SCALE units from the origo along the z-axis
oy = 0
oz = SCALE
xp = 50
yp = 0
zp = 30
loop for each frame:
reset_matrix(tm)
reset_matrix(t_om)
tm[0,0] = cos(xa)*cos(za)
tm[0,1] = cos(ya)*sin(za)
tm[0,2] = -sin(ya)
tm[1,0] = sin(xa)*sin(ya)*cos(za)-cos(xa)*sin(za)
tm[1,1] = sin(xa)*sin(ya)*sin(za)+cos(xa)*cos(za)
tm[1,2] = sin(xa)*cos(ya)
tm[2,0] = cos(xa)*sin(ya)*cos(za)+sin(xa)*sin(za)
tm[2,1] = cos(xa)*sin(ya)*sin(za)-sin(xa)*cos(za)
tm[2,2] = cos(xa)*sin(ya)
matrix_mul(om,tm) ; transform the object matrix
matrix_mul(t_om,om) ; transform the world matrix
; rotate the point
x = xp*t_om[0,0] + yp*t_om[1,0] + zp*t_om[2,0] + ox
y = xp*t_om[0,1] + yp*t_om[1,1] + zp*t_om[2,1] + oy
z = xp*t_om[0,2] + yp*t_om[1,2] + zp*t_om[2,2] + oz
x = x * SCALE / z + ORIGO_X ; 2D transforming
y = y * SCALE / z + ORIGO_Y
putpixel(x,y)
stop_looping
NOTE! The object rotates steadily as long as the angles are not touched, in other words it stops spinning not until the angles are set zero. Why? You see, the object matrix is not reset anywhere, so it rotates always one degree / frame if the angle is kept constant. This kinda technique is called cumulative rotating, and I strongly suggest you use it too in your own engine.

2.6 Hierarchical Transformations

Apologies to the OTM people, this part is almost straight from otmmatx.doc too. X)

Have you ever thought how in the world are implemented games where vector-based human-like objects move nicely (Tomb Raiders, Quakes, ...)? The thing is done with hierarchical transformations.

Let's take a human arm as an example. When you move your arm from the shoulder, the whole arm moves: wrist, palm, fingers, everything. If you move just your wrist, the arm is not moved above it but palm and fingers are. The different parts of an arm are thus hierarchically dependent on each other: the palm can command the fingers but not the wrist; it's located between them in the hierarchical order. With the same method we can create a linked list of all parts of the body. In hierarchical systems, an object is thus divided to parts, and different parts are dependent on other parts. But how to implement this?

Let's continue with the arm example. The wrist matrix is called R (the Finnish acronym for a wrist), the palm matrix K and finger matrices S1, S2, S3, S4, and S5. They're located in the linked list as follows:
 


 
We can rotate the wrist like this:

R = R*XYZr.

The palm must be rotated not only with its own rotation matrix but also with the formula above,
K = K*XYZk*R, and fingers with the formula

S? = S?*XYZs?*K.

Note! To use these formulas, the preceding matrix should be already transformed, so the calculating order is in the linked list from top to bottom.

2.7 Inverse Transformations

(I wrote this chapter just to annoy Altair a bit :) Using inverse transformations, we can get rid of vertex normal rotating and perform backface culling (object culling / portals if you wish, too) before rotating the vertices. This means of course a remarkable speedup when about 80% of the scene is suddenly dropped out the calculations even before rotations (that is, before doing practically anything). Even better, this all is very simple to implement! I'll take the vertex normal thing as an example.

We use vertex normals to calculate the amount of light hitting a surface. So, why bother rotating these all using the object matrix when we can achieve the same result by rotating the light(s) to the opposite direction? Think about it: you've got a cube (located in the origin) rotated 90 degrees about the y-axis, and a light source in the negative z-axis infinity pointing towards the positive z-axis. The most brightness is now got by the side that was originally facing the positive x-axis (but currently facing the negative z-axis). If we rotate the light -90 degrees, we get the light pointing towards the negative x-axis and the highest amount of light hits the same side as when rotating the vertex normals 90 degrees. So it seems to work. :)

Nice, but how to implement it? Pas de problem:  find the inverse of the object matrix, transform the light by it and voila! Since object matrices are always orthogonal (see the matrix stuff in the 1 part of the doc), we can use the transpose as the inverse.

Backface culling, object culling, and portals use the same idea. I'm describing only backface culling in detail but it can be applied to the others as well.
So, there are three steps:
1. Find the inverse of the object matrix (transponate it)
2. Transform the camera vector and position by it
3. Perform the following for each face in the object:
a) find the vector between the camera position and any vertex in the face
b) calculate the cosine of the angle between this vector and the camera vector (that is, perform a simple dot product)
c) examine its sign: if it's positive, the face can be seen, otherwise not.
Quite straightforward, isn't it? What you've maybe wondering is, why any of the face vertices suffices. Think: the idea of backface culling is that if a polygon is not facing the camera, it certainly can't be seen. So we need only one point in the plane on which the polygon is lying.

Hmm... Maybe a little picture clarifies the case a bit :)
 


 
C = inverse-rotated camera vector, alpha = the angle the cosine of which is wanted, N = the normal of the face (no actual need for it, put it there just to make the picture a bit clearlier :)

Ok, that's it. Just GO USE 'EM! :)