3D Coding TAT - Projections And Corrections
The problem with the monitor screen is that it is flat, just a 2d bitmap on which 3d objects have to be drawn. It is at this stage where one dimension is lost, discarded forever. All the preceeding transformations do not lose any information (apart from precision errors) so to a great extent they can be reliably undone and their calculations reversed, but this is not the case with the perspective (projection) 3d-to-2d stage because we would need to create a dimension from nothing.
The perspective stage is normally done with a projection of some kind. The most simple is to divide the X and Y coordinates by the Z, but this is incorrect. The basic idea behind the simulated perspective is that the further away things are the smaller they should be. Using just the Z coordinate as the distance from the view point is only ever correct when X=0 and Y=0, otherwise the square-root of (Xý + Yý + Zý) should be used. This calculation is very slow because 3x multiplys and 1x square-root needs to be done PER VERTEX! Using 2 divisions is a reasonable cheat.
There are some subtle characteristics in the normal projection formulas which can cause some viewing errors and prevent any optimizations which 'might' be possible from combining the rotation and projection stages. If a 3d object model is projected with its centre origin at (0,0) then the distance errors are hardly noticable but move the object to the edge of the view-port/screen and it begins to sheer in one or two directions. Another error can be seen if an object is rotated around the view point while keeping a constant radius. Again near the edges of the view-port/screen distortion occurs, the object appears to grow in size because its Z distance has decreased and its X and/or Y has increased. This is why the square-root method gives the correct result whereas the divide-by-Z does not. The greater the resolution of the view-screen in terms of pixels, the greater the amount of this distortion. You may want to adjust the Z coordinate using a small amount of the X and Y coordinates to give a closer approximation to the true square-root or you could simply push the object slightly further away depending on it's angle to the (0,0) origin on the view-screen. So as an object gets closer to one of the screen edges you move it slightly into the distance BEFORE projecting its vertices.
Okay you might say, why not simply project the object on the (0,0) origin and THEN translate it to the correct view-port/screen position? Well this would help minimize sheering/distance errors but the parallax would be incorrect. Try holding a box directly in front of your eyes and then move it to the left or right. The sides of the object appear to shrink or grow because the front corner vertices seem to move faster than the back ones (i.e. parallax). You can think of this another way (which may or may not be correct), imagine that object has rotated slightly. Hold an object directly in front of your eyes and twist it to the right a little, now move it right a foot or so. You should hopefully see that this is almost (if not entirely) indentical to the proper square-root distancing. Because we are rotating and projecting around the object's (0,0,0) origin when it is close to the view-screen edges the sheering and size distortion should be minimize a great deal (if not completely removed).
There is some extra calculation involved to work out the amount of rotation which is needed. The amount of rotation depends on the angle between the view-port/eye (0,0) origin and the object's central origin. One way to do this it take the correct distance of JUST the object's origin (relative to the view-point) together with it's angle from the view-screen (0,0) origin and use these to rotate and project all the object's vertices.
X = 0
e.g. ³ . .
³ . object
³ . o .
³ / .
Z ---------------Å----/----------- view-plane
³ / (the screen)
³ ³a / r
³ ³ /
Å---- X ³/ PLAN VIEW
eye
In the above diagram the object's origin is marked by 'o' the object vertices are marked by '.' and 'a' is the angle between the object's-origin and the X=0 plane. The 'r' is the true length of the object's-origin to the eye-point (from 'o' to the 'eye'). This diagram only shows the X-Z plane, it should be done for the Y-Z plane as well so the true 3d distance is found and the 2 angluar adjustments.
1: The Angle Grinder
Here is a little algorithm which I use to calculate the angle between two points using the X and Y axis lengths:
IF x <> 0 THEN angle = ATN(ABS(y)/x)*180/PI
ELSE angle = 0
IF x < 0 THEN angle = angle + 180
IF y < 0 THEN angle = 360 - angle
It works for all 360° in a fairly compact way. If you use some other resolution for a complete 360° circle then just replace 360 with your resolution and 180 with 1/2 your angle resolution. For example if you use 1024° for one complete circle then use:
degrees = 1024
half = degrees / 2
IF x <> 0 THEN angle = ATN(ABS(y)/x) * half/PI
ELSE angle = 0
IF x < 0 THEN angle = angle + half
IF y < 0 THEN angle = degrees - angle
Of course the "half/PI" value can be replaced with a constant or precalculated value. This would remove a divide instruction.
NOTE: The angles are converted from DEGREES into RADIANS for the ATN function. There are 2*PI radians in a circle or 360 degrees.
So we end up with a divide and a multiply, but the real problem is the ATN (Arc TaNgent) function which is difficult to do in 80x86 (I not even sure that there are any floating-point instructions which perform it, I could be wrong). After pausing to take a quick look at the old Spectrum ROM book the correct polynomial method looks way too long and far too slow, so I suggest using a LUT (look-up table). After a little bit of experimentation I found that basing it around 45° arcs works reasonably well and using the shortest length divided by the longest length can be used to index into an correction table, something like this:
;* First create an adjustment table *
DIM adjust[45]
FOR a = 0 to 45
n = ( SIN(a*PI/180) * 45 / COS(a*PI/180) )
adjust[n] = a - n
NEXT a
;* Now here is the angle-finding code *
angle = 0
turn90:
IF x < 0 OR y < 0 THEN
SWAP x,y
y = 0 - y
angle = angle + 90
GOTO turn90
IF x = 0 THEN angle = angle + 90 : RETURN
IF y = 0 THEN RETURN
IF y > x THEN
n = 45 * x / y
angle = angle + 90 - n - adjust[n]
ELSE
n = 45 * y / x
angle = angle + n + adjust[n]
RETURN
It's far from perfect, but it works (sort of). It can be improved slightly by building two look-up tables (one for y>x and one for x>y). You could also detect 45° angles simply by comparing the x and y lengths.
2: Round the Bend
The difficulty in writting a ATN (Arc Tangent) or ARCSIN or ARCCOS function (the last two are slower) is that we must find an angle based upon an X,Y position on a circle. We need to find the reverse of a SIN (sine) or COS (cosine) function, that is the angle. There is no easy, quick way to do this. I haven't found any descriptions or algorithms which perform these functions which is suprising considering how useful they are for direction finding, anti-rotations and so on...
3: Hint, Hint, Intel
What I don't understand is why chip makers don't create floating-point instructions to performs some of the most used formulas in maths. Some functions like these would be a god send.
a = FANGLE ( x, y )
l = FMAG2D ( x, y )
l = FMAG3d ( x, y, z )
Where FANGLE returns the angle from the two operand lengths and FMAG2D returns the magnitude (SQRT of X*X + Y*Y) and FMAG3d returns the magnitude of SQRT (X*X + Y*Y + Z*Z). You can easily find the minimum values for the FMAG2D and FMAG3d functions just by finding the largest ABSolute length.
e.g.
major = ABS(x)
minor = ABS(y)
IF major < minor THEN SWAP major, minor
Now major is the minimum value of SQRT(X*X + Y*Y). You can also find the maximum value by:
2D length 3d length
major*SQRT(2) or major*SQRT(3)
4: Faster Perspective/Projection
As the 3d vertices must be converted into 2d ones there is usually the need to perform a division and a multiplication together. Say we had two calculations like the ones below.
X Y
X' = xratio * --- Y' = yratio * ---
Z Z
The 'xratio' and 'yratio' values are used to modify the angle of the projection and change the field-of-vision. Most people use a value of 128, 256 or 512 for these variables because they can be done using a very fast shift instruction, so by removing an multiply.
This should look familiar to most people as a form of these calculations are mainly used for the 3d to 2d conversion.
But we have 2x multiplies and 2x divides which are very slow compared to additions or subtractions. One way to help speed these up is to use a table of 1/n values, so instead of performing a division you perform a multiplication.
e.g.
a = b * table[c]
where each element in table[n] = 1/n (the reciprocal), so we are doing the same as:
1 b
a = b * --- which equals: a = ---
c c
Going back to the projection/perspective formulas, if the xratio equals the yratio then we can avoid the second division by reusing the result of the first.
NOTE: only IF xratio = yratio
e.g.
xratio yratio
scalar = --------- = ---------
Z Z
X' = scalar * X Y' = scalar * Y
So this turns out to be 1 division and 2 multiplies, but you may need to use fixed-point maths to keep the precision of 'scalar' to a modest level.
5: Combining Rotations & Projection
In the general render-pipe-line the projection/perspective calculations usually follow some rotations, either matrix or straight forward 3-stage X, Y then Z axis rotations. As you should already know the X and Y coordinates need to be magnified up in the projection part before being divided by the Z coordinate. These days everyone has a Pentium or at least a fast 486 slug, so floating-point instructions are finally worth using rather than integer ones. This means that using fractions and combining operations now doesn't mean losing tons of precision. So operations or formulas can be combined without too many problems.
Take a simple rotation around the Z axis (pointing into the screen) followed by a basic projection.
Z ROTATION:
x' = x * cos(a) - y * sin(a)
y' = x * sin(a) + y * cos(a)
z' = z
The above cos(a) and sin(a) values can be placed in a matrix. This helps when dealing with lots of vertices. In this case it could look like this:
x' = x * m[11] - y * m[21]
y' = x * m[12] + y * m[22]
z' = z
where:
m[11] = cos(a) m[21] = sin(a)
m[12] = sin(a) m[22] = cos(a)
Here is the basic projection stage.
PROJECTION:
x' * Xratio
h = ------------- + Xcentre
z
y' * Yratio
v = ------------- + Ycentre
z
where:
Xcentre, Ycentre are the centre coordinates of screen
Xratio, Yratio are the field of vision magnifiers
They usually have the values something like these:
Xcentre = 160
Ycentre = 100
Xratio = 500
Yratio = 500
Right enough waffle, here is the good bit.
When the rotation matrix is being built simply combine the Xratio and Yratio values into it.
e.g.
m[11] = cos(a)*Xratio m[21] = sin(a)
m[12] = sin(a)*Yratio m[22] = cos(a)
Now the rotation is,
x' = x * m[11] - y * m[21]
y' = x * m[12] + y * m[22]
z' = z
and the projection is,
x'
h = ----- + Xcentre
z
y'
v = ----- + Ycentre
z
So we have removed 2 multiply instructions for each vertex. Well actually we have moved them into the matrix set up code. This means if we rotation and project N vertices we can save 2 * (N-1) multiply instructions, not bad eh?