A tutorial for 2d and 3d vector and texture mapped graphics
Version 0.32ß

© 1994 by Sébastien Loisel

Note
This document is provided "as is", without guaranty of ANY kind. This document is copyrighted 1994 by Sébastien Loisel. It may be distributed freely as long as noo fee except shipping and handling is charged for it. Including this document in a commercial package requires the written permission of the author, Sébastien Loisel.

If you want to reach me, see the end of this file.

	____________________

1	An introduction to 3d transformations

First, let's introduce the idea of projections. We have to realize that we are trying to view a 3d world through a 2d «window». Thus, one full dimension is lost and it would be absurd ot think that we can reconstitute it wholly and recognizably with one less dimension. However, there are many forms of depth-cueing that we introduce in a 2d image to help us recognize that 3d dimension. Examples include but are not limited to reduced lighting or visibility at a distance, reduced size at a distance, out-of-focus blur, etc...

Our media is a screen, ideally flat, and the 3d world being represented is supposed to look as real as possible. Let's look at how the eyes see things:

Figure 1

The projection rays are the beams of light that reflect off the objects and head to your eye. Macroscopically, they are straight lines. If you're projecting these beams on a plane, you have to find the intersection of these beams with the plane.

Figure 2

The projection points are thus the intersection of the projection rays with the projection plane. Of course, projection points coming from projection lines originating from behind an object are obviously hidden, and thus comes the loss of detail from 3d to 2d. (i.e. closer objects can hide objects that lie further away from the eye).

The objects can be defined as a set of (x,y,z) points in 3d space (or 3space). Thus, a filled sphere can be defined as x**2+y**2+z**2<=1. But we need to simplify the problem right away because of machine limitations. We can safely assume the eye is always outside any object, thus we do not need do define the interior. [Note: a**b stands for «a raised to the bth power», thus, x**2 is x squared]

Now let us define a projection onto a plane. Many different projections could be defined, but the one that interests us is the following for it's simplicity. Projecting takes any given point (x,y,z) and localizes it on an (u,v) coordinate system. I.e. space (x,y,z) is mapped onto a plane (u,v). Let's transform (x,y,z) first. Multiplying (x,y,z) by a constant keeps it on a line that passes through (0,0,0) (the origin). We will define that line as the projection ray for our projection. Since all of these lines pass through the origin, we thus define the origin as our eye. Let the constant be k. Thus our transformation as of now is k(x,y,z). Now we want to set everything into a plane. If we select k as being 1/z, assuming z is nonzero, we get a projected point of (1/z)(x,y,z). This becomes (x/z, y/z, 1). Thus all points of any (x,y,z) value except z=0 will be projected with a value z'=z/z of one. Thus all points lie in the plane z=1. This projection thus fits our need for a projection with linear projector through the eye on a plane surface.

If you didn't get all that right away, you can always look it up later. The important thing to remember now is this: to project from (x,y,z) onto (u,v) we decide to use

	u=x/z
	v=y/z

But, I hear you say, this means the eye is facing parallel to the z axis, centered on the origin and is not banked. What if that is not what I desire?

First, it is obvious that if the eye is at location (a,b,c) it will not affect what one would see to move EVERYTING including the eye by (-a,-b,-c). Then the eye will be centered.

Second, and almost as obvious, if the eye is rotated around the x,y and z axis by agles of rx, ry and rz respectively, you can rotate EVERYTHING by the exact opposite angles, which will orient the eye correctly without changing what the eye would see.

Translation (moving) is a simple addition. I.e. (x,y,z)+(a,b,c) becomes (x+a,y+b,z+c). Additions are fairly quick. But there is the problem of the rotations.

1.1 Trigonometry

We will do it in two dimensions first.

Figure 3

We have the point defined by (x,y) and we want to rotate it with a counterclockwise angle of b to (x',y'). See figure 3. The radius r stays the same throughout the rotation, as per definition, rotation changes the angle but not distance to origin. Let r be the radius (not shown on figure 3). We need to express x' and y' as functions of what we know, that is x, y and b. The following can be said:

y'=sin(a+b)r
x'=cos(a+b)r

With the identities sin(a+b)=sin(a)cos(b)+sin(b)cos(a) and cos(a+b)=cos(a)cos(b)-sin(a)sin(b), we substitute.

y'=rsin(a)cos(b)+rcos(a)sin(b)
x'=rcos(a)cos(b)-rsin(a)sin(b)

But from figure 3 we know that

rsin(a)=y and
rcos(a)=x

We now substitute:

y'=ycos(b)+xsin(b)
x'=xcos(b)-ysin(b)

A special note: you do not normally calculate sinuses and cosinuses in real time. Normally, you build a lookup table, that is, an array. This array contains the precalculated value of sin(x) for all values of x within a certain range. That is, if the smallest possible angle in your model is on half of a degree, then you need a sinus table containing 720 entries, of which each entry is the value of the sinus of the corresponding angle. You can even note that a sinus is highly symmetric, that is, you need the values of sin(x) only for x between 0 and 90 degrees. Furthermore, you do not need a cosine table, as cos(x)=sin(x+90 degrees). This and fixed point mathematics may become obsolete as floating point processors become faster and faster. A good algorithm for calculating a sinus can be very quick today on a reasonably high performance platform, and it shouldn't use powers.

With that said, we will approach briefly 3d transformations.

In 3d, all rotations can be performed as a combination of 3 rotations, one about the x axis, one about the y axis and one about the z axis. If you are rotating about the z axis, use the exact same equation as above, and leave the z coordinate untouched. You can extrapolate the two other transformations from that.

Note that rotating first about the x axis and then about the y axis is not the same as rotating about the y axis and then the x axis. Consider this: point (0,0,1) rotated 90° about the x axis becomes (0,1,0). Then, rotated 90° about the y axis it remains there. Now, rotating first about the y axis yields (1,0,0), which rotated then about the x axis stays at (1,0,0), thus the order of rotations is important [(0,1,0) and (1,0,0) are not the same].

Now, if you are facing west (x axis, towards negative infinity). The user pulls the joystick up. You want the ship to dive down as a result of this. You will have to rotate the ship in the xy plane, which is the same as rotating about the z axis. Things get more complicated when the ship is facing an arbitrary angle.

	____________________

2	Polygon drawing on a raster display in two dimensions

The second main problem is to display a filled polygon on a screen. We will discuss that problem here. First, let us approach the 2d polygon.

We will define a polygon as a set of points linked by segments, with a defined inner region and outer region. The inner region need not be in a single piece.

Our rule for determining if a point is in a given region is to draw a line from infinity to that point. If that line crosses an odd number of segments (edges), the point is within the polygon. Otherwise it is not.

Figure 4

Drawing a line from infinity to any of the gray areas of the star will cross an odd number of edges. In the white areas, it will cross an even number of edges. Try it.

Figure 5

The polygon shown in a, b and c are all the same. In a, the inner region has been shaded, but this has been omitted in b and c for clarity.

The polygon shown in a can be listed as vectors, as in b or c. Let us agree that both versions are equally satisfying. In b we enumerated the edges in a counterclockwise continuous fashion, meaning that the head of each vector points to the tail of another. In c, we made each vector point generally downward, or more strictly, each vector's topmost point is its tail. These vector are noted (a,b)-(c,d) where (a,b) is the location of the tail of the vector, and (c,d) is the location of the head (arrow) of the vector. In diagram c, for all vectors d<b (or b>d, either way). (Note, we will see what to do with b=d later. That being a special case, we will not discuss it now).

As we are drawing on a screen (supposedly), it would be interesting to be able to color each pixel once and only once, going in horizontal or vertical lines (as these are easy to compute). We will choose horizontal lines for hardware reasons. Thus, we start scanning from infinitely upwards to infinitely downwards using horizontal scanlines (e.g., the «first» scanline is way up, and the last is way down). Now, our screen being finite, we need only to draw those lines that are on screen. At each pixel, we will draw an imaginary line from the infinite to the left (horizontal line) to the current pixel. If it crosses an odd number of edges, we color it. Otherwise, we don't.

This means that from infinite left to the first edge, we do not color. Then, from first edge to second edge, we color. Then, from second edge to third edge, we do not color. From third edge to fourth edge, we color. And so on. If you want, you can always draw a background bitmap left from the first edge, from second edge to third edge, etc.... Thus, you could put the picture of a beautiful blue sky behind at the same time you're drawing a polygon.

Now, let's start that over again. We are coming from infinitely upward and drawing. We do not need to consider every edge, only those that are on the current scanline (horizontal line). We will thus build an inactive edge list (those edges that have not yet been used in the drawing process) and an active edge list (those edges that are currently in use to draw the polygon, i.e. that intercept the current scanline). Thus, at first, all edges (named m,n,o,p in figure 5) would be in the inactive edge.

As the algorithm reaches the topmost line of the polygon, any relevant edge will be transfered from the inactive edge list to the active edge list. When an edge becomes unused (we're lower than the lowermost scanline), we remove it from the active list and throw it to the dogs. (# grin #)

Now you will realize that at each scanline we have to check every edge in the inactive edge list to see if we need to transfer them to the active edge list. To accelerate that, we can sort the inactive edge list in decreasing order. Thus, the ordered list would become {m,p,n,o}. Then, we need only to check the first vector at every scanline. Since they point downward, we need only to check the first point of it.

If you get any horizontal lines in the edge list, remove them from the list. Horizontal lines can be identified easily as b=d. Proving that it works is left as an exercise (# muhaha! #).

The polygon drawing algorithm (scan-line algorithm).

I=inactive edges list, initially all edges
A=active edges list, initally empty

sort I with highest first endpoint first

for each scanline i starting from the top of the screen
	does the first edge in I start at scan line i?
		Yes, transfer all edges that should start now from I to A
	find x value of all edges in A for current scanline (see the line algorithm below)
	sort all edges in A with increasing x-intercept value (the thing we jest calculated)

	previous_x_intercept<-leftmost pixel on screen
	for all edges n in A, counting by two, do
		draw background between previous_x_intercept and edge n in A
		fill polygon between edge n and edge n+1 in A
		previous_x_intercept<-x intercept of edge n+1
	end for

	draw background for the rest of the scan line (from previous_x_intercept to end of line)
	for all edges for which we have reached the bottom
		remove the edge from A
	end for
end for

2.1 A line algorithm for the polygon drawing algorithm

This section discusses a somewhat shortened version of Bresenham's line drawing algorithm. It is optimized for drawing filled polygons.

You assuredly realize you need to calculate the x coordinate of the intersection of a line with an horizontal line. Well consider the line equation that follows:

	x=dy/dx x y + k

If the line is defined with the vector (a,b)-(c,d), then we define dy as d-b, dx as c-a. Say x0 is the x coordinate of the topmost point of the segment.

x0=a

Say xn is the x intercept value of the segment for the nth line from the top. Then, xn can be defined as follows:

xn=dy/dx x n + a

Say we know xn-1 and want to know xn based on that knowledge. Well it is obvious that:

xn-x(n-1)=dy/dx x n + a - [dy/dx x (n-1) + a]
xn-x(n-1)=dy/dx
xn=dy/dx+x(n-1)

Thus, knowing xn-1 and adding dy/dx yields xn. We can thus do an addition instead of a multiply and an add to calculate the new x intercept value of a line at each scanline. An addition been at least 5 times as fast as an add, usually more in the 15 times range, this will speed things up considerably. However, it might annoy us to use floating point or fixed point, even more so with roundoff error.

To this end, we have found a solution. dy/dx can be expressed as a+b/d, b<d. a is the integer part of dy/dx. Substituting, we find that

xn=a+b/d+xn-1

Thus, we always add a at each scanline. Furthermore, when we add a sufficient amount of b/d, we will add one more to xn. We will keep a running total of the fraction part. We will keep the denominator in mind, and set the initial numerator (when n=0) to 0. Then, at each line we increase the numerator by b. If it becomes greater than d, we add one to xn and decrease the numerator by d. In short, in pseudocode, we draw a segment from (x0,y0) to (x1,y1) this way:

dx=x1-x0
dy=y1-y0
denominator=dy
numerator=dx modulo dy
add=dx divided (integer) dy
running=0
x=x0
for each line (0 to n) do
	x=x+add
	running=running+numerator
	if running>=denominator then
		running=numerator-denominator
		x=x+1
	end if
	do whatever you want with line here, such as drawing a pixel at location (x,line)
next line

	__________________________

3	3d polygon drawing

The first step is to realize that a projected 3d line with our given projection becomes a 2d line, that is it stays straight and does not bend. That can be easily proved. The simplest proof is the geometric one: the intersection of two planes in 3d is a line in the general case. The projector all lie in the plane formed by the origin ant the line to be projected, and the projection plane is, well, a plane. Thus, the projection is contained in a line.

This means that all the boundaries of a polygon can be projected as simple lines. Only their endpoints need be projected. Then, you can simply draw a line between the two endpoints of each projected edge (or more precisely, draw a 2d polygon with the given endpoints).

However, there comes the problem of overlapping polygons. Obviously, the color of a pixel must be that of the closest polygon in that point (well for our objective, that suffices). But this means that we must determine the distance of each polygon that has the possibility to be drawn in the current pixel, and that for the whole image. That can be time consuming.

We will make a few assumptions. First, we will assume no polygons are interpenetrating; this means that to come in front of another polygon, a polygon must go around it. This means that you only need to determine what polygon is topmost when you cross an edge.

Less obviously, but equally useful, if from one scanline to the next, you encounter the same edges in the same order when considering left to right, you do not need to recalculate the priorities as it has not changed.

Inevitably, though, you will need to determine the depth of the polygon (z value) in a given pixel. That can be done rather easily if you carry with the polygon its plane equation in the form Ax+By+Cz+D=0. Remember that the projection ray equation for pixel (u,v) is x=zu and y=zv. Feed that x and y in the plane equation and you find that

Azu+Bzv+Cz+D=0
z(Au+Bv+C)=-D
z=-D/(Au+Bv+C)

Thus you can easily calculate the z coordinate for any plane in a given (u,v) pixel.

The cautious reader will have noticed that it is useful to do some of these multiplications and additions incrementally from scanline to scanline. Also note that if Au+Bv+C equals zero, you cannot find z unless D is also zero. These (u,v) coordinates correspond to the escape line of the plane of the polygon (you can picture it as the line of «horizon» for the given plane). The projection ray is parallel to the plane, thus you never cross it.

This brings me to another very important thing about planes. Another way to define a plane is to give a normal vector. A vector and a point of reference generate one and only one plane perpendicular to the vector and passing through the point of reference. Think of a plane. Think of a vector sticking up perpendicular to it. That vector is called normal to the plane. It can be easily proved that this vector is (A,B,C) with A,B,C being the coefficients of the plane equation Ax+By+Cz=D. We might want to normalize the vector (A,B,C) (make it length 1), or divide it by SQRT(A2+B2+C2). But to keep the equation the same, we must divide all member by that value. The equation thus becomes A'x+B'y+C'z=D'. If we want to measure the distance perpendicularly to the plane of a point (a,b,c), it can be shown that this (signed) distance is A'a+B'b+C'c+D' (you can always take the absolute value of it if you want an unsigned distance). That might be useful for different reasons.

It is also notable that since the A, B and C coefficients are also the coords of the normal vector, and that rotating a vector can be done as we saw previously, it is easy to rotate the plane. If the vector is normalized in the first place, it will stay normalized after rotation. After rotation, you will ignore the D component of the plane, but you do know the coordinates of a point in the rotated plane. Any point in the polygon will fit the bill. Thus, you can deduce D. In fact, when performing back-face culling (see chapter 4), you will calculate D even as you are performing some optimization.

The algorithm for drawing the polygons is basically the same as for 2d polygons, with a few twists. Have fun finding them.

Another very nice algorithm is the z-buffer algorithm, discussed later. It allows interpenetrating polygon and permits mixing rendering techniques in a rather efficient way (for real-time rendering purposes).

	_________________________

4	Data Structures

A very important aspect of 3d graphics is data structures. It can speed up calculations by a factor of 6 and up, only if we put some restriction on what can be represented as 3d entities.

Suppose that everything we have in the 3d world are closed polyhedras (no self-crossing(s)!). Each polyhedra is constituted of several polygons. We will need each polygon's normal and vertexes. But how we list these is where we can improve speed.

It can be proved that each edge is shared by exactly two polygons in a closed polyhedra. Thus, if we list the edges twice, once in each polygon, we will do the transformations twice, and it will take twice as long.

To that, add that each vertex is part of at least 3 points, without maximum. If you list each point in every vertex it is part of, you will do the transformations at least thrice, which means it will take thrice the time. Three times two is six. You're doing the same transformations at least six times over.

Thus, you need to list the vertexes in every polygon as pointers to vertexes. This way, two polygons can share the same data for a shared vertex, and will allow you to do the transformations only once.

The same goes for edges, list the endpoints as pointers to endpoints, thus you will spare yourself much runtime work.

Another interesting thing. All polygons have an interior and exterior face in that model. Make it so that the normal to all polygons point towards the exterior face. It is clear that, if you are standing outside all objects, and all objects (polyhedra) are closed, you cannot see an inside face. If, when transforming polygons, you realize a polygon is facing away from the eye, you do not to transform it, just skip it. To know if it's pointing away from you do the following:

This is called back-face culling. Take the (a,b,c) absolute coordinate of any one endpoint of the polygon. Take also (u,v,w) the normal vector to the polygon. No matter what direction you are looking, do the following operation:

au+bv+cw

If the result is positive, then it is facing away from you. If it is negative, it is facing towards you. If it is null, it is parallel to you. The operation is called scalar multiply of two vectors. It is very notable that the above scalar multiplication of two vectors is interestingly enough the D part of the plane equation Ax+By+Cz=D.

You can also illuminate objects with a very distant light source (e.g. the sun) very easily. Take the vector pointing away from the light source (a,b,c) and the normal of the plane (u,v,w). Do a scalar multiplication of both as above. The result can be interpreted as the intensity of the lighting or shadow for the polygon. The approximation is good enough to give interesting results. You might also want to affect the shading of the polygon according to its distance to origin, to yield a fog result.

	_______________________

5	Texture mapping

What we are about to do would probably be more aptly named pixmap mapping onto a plane. We want to take a pixmap (or bitmap, although less appropriate a word, it is very popular) and «stick» it to the surface of a polygon, cutting away any excess. The result can be quite impressive. This, however, can be a painfully slow operation. We will discuss here of the mathematics involved, and after that, we will try to optimize them.

Let's give a polygon its own coordinate system, (i,j). Thus, all points on the said polygon can be localized in (i,j), or (x,y,z), or, once projected, in (u,v).

Figure 6

Figure 6 illustrates the problem. The blue plane is the projection plane. The green triangle is a polygon we are projecting on that plane and then scan-converting in (u,v) coordinates. What we want to do is find, what (i,j) point corresponds to any given (u,v) point when projected this way. Then, the point (i,j) can be indexed in a pixmap to see what color this pixel in (u,v) space should be.

Let the i and j vectors be expressed in (x,y,z) coordinates, which they should normally be. We want a general solution.

i=(a,b,c)
j=(d,e,f)

Let also the origin of (i,j) space be pointed to in (x,y,z) by the vector k=(m,n,o). Thus, a point at location (p,q) in (i,j) space is the same as a point at location pi+qj+k in (x,y,z) space. Furthermore, we have the projection ray defined as follow: it shoots through point (r,s) in (u,v) space. This corresponds to point (r,s,1) in (x,y,z) space. The equation of that ray is thus t(r,s,1) where t can take any real value. Simplifying, we find that the ray is (tr,ts,t), or

	x=tr
R:	y=ts
	z=t

or
	x=zr
R:	y=zs
	z=z

The point in the plane is
	x=pa+qd+m	(1)
P:	y=pb+qe+n	(2)
	z=pc+qf+o	(3)

We want to find (p,q) as a function of (r,s). Here is what i did in stupid MathCAD. [MathCAD flame deleted] I'm sorry I didn't use the same variable names in MathCAD...

#######Begin picture of mathcad screen#########
########End picture of mathcad screen#########

These equations for p and q can be optimized slightly for computers. Since it is usually faster on a computer to do a(b+c) than ab+ac, and that both are equivalent, we will try to factorize a little.

p=(Mv+Nu+O)/(Pv+Qu+R)
where M,N,O,P,Q,R are as follows:
M=CXq-AZq
N=BZq-CYq
O=AYq-BXq
P=ZqXp-ZpXq
Q=YqZp-YpZq
R=YpXq-YqXp

They are all constants, so they need only to be calculated once to draw the whole image.

Similarly for q, we have:
q=(Sv+Tu+U)/(Vv+Wu+X)
and
S=AZp-CXp
T=CYp-BZp
U=BXp-AYp
V=ZqXp-ZpXq
W=YqZp-YpZq
X=YqXp-YpXq

All these constants should be calculated only once. Now we can simplify a little with, knowing that P=V, Q=W, R=X.

p=(Mv+Nu+O)/(Pv+Qu+R)
q=(Sv+Tu+U)/(Pv+Qu+R)

This is noteworthy: (Q,P,R)(u,v,1)=uQ+vP+R, the denominator. But, (Q,P,R)=(Xq,Yq,Zq)x(Xp,Yp,Zp). This cross-product of two vectors will yield a vector that is perpendicular to the first two (i.e. a normal to the plane). If these two vectors are of length 1, then the cross product will also be of length 1. But we might already have the normal from previous transformations (e.g. if we're doing back-face culling). Then, we won't need to calculate (Q,P,R) as it is the normal to the plane. Or we can deduce the normal from these equations if we don't have it.

The denominator needs to be calculated only once per pixel, not for both p and q. The numerator for p and q is different, though. Furthermore, It is obvious that Mv+O are constant throughout a constant-v line (horizontal line). Thus, it needs to be calculated only once per line. We can spare ourselves the multiply per pixel doing it incrementally. (We're still stuck with a divide, though.) Example, let V=Mv+O, for a given v. Then, when calculating Mv+Nu+O, we can use V+Nu. If we know V+Nua, then we can calculate V+N(ua+1) as follow: V+Nua+N

Let W=V+Nua, then we get

W+N

That's a simple add instead of a multiply. But we can not get rid of the divide, so we're stuck with a divide per pixel (and at least two adds) per coordinate.

There is a way to solve this, however, and here's how to do it.

First, it is interesting to note the following. The plane equation is Ax+By+Cz=D, and x=zu, y=zv, where (u,v) is the screen coordinates of a point, and (x,y,z) is the corresponding world coordinate of the unprojected point. Thus,

Auz+Bvz+Cz=D

What happens when z is a constant? Let us express v as a function of u.

Bvz=-Auz-Cz+D
v=(-A/B)u+(-C/B)+D/(Bz)

Let M=-A/B, and N=d/(Bz)-C/B then we have

u=Mv+N

a linear equation. Furthermore, M is independant of z. Thus, a slice of constant z in the plane is projected as a line on (u,v). All of these slices project in lines that are all parallel to each other (e.g. the slope, M, is independant of the slice taken).

Now, a slice of constant z in the original Ax+By+Cz=D plane is a line. Since z is constant throughout that line, all points of that line are divided by the same number, independantly of (x,y). Thus, the projected line is merely the original line scaled up or down.

Up to now, we have been rendering polygons with a scanline algorithm where the scanlines are horizontal. First, we have to conceive that the scanlines need not be horizontal, merely parallel. If we use a slope of M for the scanlines (M=-A/B, as above), then we will be performing scans of constant z value. Two advantages are very obvious. First, it is easy to calculate the z value for a large number of pixel (it is constant throughout a scanline) and second, texture mapping the scanline becomes a mere scaling operation, which can be done incrementally. Furthermore, M need be calculated only once for the whole polygon, while N requires at most a divide and an add per scanline.

Note that if abs(A)>abs(B), you might want to express v as a function of u instead of the above, so that the slope ends up as being -B/A instead. (E.g., make certain that the slope is less than or equal to 1 in absolute value). There is a degenerate case: A=B=0. This is the case where the whole plane is of constant z value. In that case, you can use horizontal scan lines and just scale the texture by a factor of 1/z.

Thus, we're going to be rasterizing the polygon with nonhorizontal scanlines of slope M. The equation of a scanline can be expressed as:

v=Mu+Vi

Where Vi is the v coordinate when u is 0. All scanlines being parallel, M never changes, and Vi identifies uniquely each and every scanline. Say i=0 corresponds to the topmost scanline. Say we want to find the intersection of the scanline with a polygon edge. What we're really interested in is the u value of the intersection of the two lines. The edge can be represented as being:

v=Nu+K

Where K is the v value for the edge when u=0. Now, intersecting the two yields

Nu+K=Mu+Vi
(N-M)u=(Vi-K)
u=(Vi-K)/(N-M)

That is, assuming N-M is not zero, which shouldn't be because we're not including edges parallel to the scanlines in our scan line algorithm, thus N-M is never zero.

Now, V(i+1) is (Vi)+1. When we go from scanline i to scanline i+1, V is incremented by one. We thus get u', the intersection of the two lines for the next scanline, as being:

u'=(Vi+1-K)/(N-M)
u'=(Vi-K)/(N-M)+1/(N-M)
u'=u+1/(N-M)

It is thus possible to calculate u from scanline to scanline incrementally, bresenham stlye. To clarify things up a bit: if P1 and P2 delimit the edge, then let (p,q) be P2 minus P1 (that is, p is the delta u and q is the delta v). Therefore, N is q/p (by definition of slope). Now, we know M to be of the form -B/A (from the plane equation, above). We need to calculate N-M for our incremental calculations, which is relatively easy. We have

N-M=q/p-(-A/B)
N-M=q/p+A/B
N-M=(qB+Ap)/(Bp)

Thus, u' is really:
u'=u+Bp/(qB+Ap)

At one point in the algorithm, we will need to calculate which scanline contains a particular point. As an example, when sorting from topmost edge downwards, we need to use a definition of top that is perpendicular to the scanline used. Then, we have the slope of the scanline, which is, say, m. If the equation of the line that interests us is:

y=mx+b

Then,  we can sort the edges according to increasing b values in the above equation. If we want to calculate b for point (i,j), simply isolate b in the above equation.

	________________________

6	Of logarithms

If we have many multiplies or divides to make, little time and lots of memory, and that absolute precision is not utterly necessary, we may use the following tricks.

Assuming x and y are positive or null real numbers, and b a base for a logarithm, we can do the following: (note: x**y means x raised to the yth power, logby denotes base b logarithm of y).

logb(xy)=logbx+logby
b**logb(xy)=b**(logbx+logby)
xy=b**(logbx+logby)

Now, you might say that taking the log or the power of a value is very time consuming, and it usually is. However, we will reuse the trick that we used to do not so long ago, when pocket computers were not available, we will make a logarithm table. It is easy to store the values of all possible logaritms within a finite, discrete range. As an example, if we are using 16 bit fixed points or integer, we can store the result of a log of that number as a table with 65536 entries. If each entry is 32 bits wide, which should be sufficient, the resulting table will use 256kb of memory. As for the power table, you can build another lookup table using only the 16 most significant bits of the logarithm, yielding a 16 or 32 bit wide integer or fixed point. This will yield a 128 to 256kb table, for a grand total of under 512kb. With today's machines, that should not be a problem. The only problem is that you must then limit your scope to 16 bits numbers. You realize of course that you can not multiply if one or both of the numbers are negative. Well, you can do it this way. Say x<0, then let u=-x, thus u>0. Then, xy=-uy, and uy you can calculate with this method. Similarly,

x/y=b**(logbx-logby)

Thus divides are also a simple matter of looking it up in the table and subtracting.

Powers are made much more simple also. I won't say anything long winded about powers, but here is, in short:

logb(x**y)=y logbx
x**y=b**(y logbx)	(*)
And from there, even better yet:
x**y=b** ( b** [ logby+logb(logbx) ] )

For y=constant, some shortcuts are usually possible. Since ax can sometimes be very quickly calculated with left shifts, you can use (*) instead of the last equation. Very common example:
320y=(256+64)y
320y=256y+64y
320y=(y<<8)+(y<<6)

Where x<<y denotes "x, left shifted y times". A left shift of y is equivalent of multiplying by 2**y. The ancient egyptians only knew how to multiply this way (at least, as far as I know, for a while, like for the pyramids and stuff - I think they even calculated an approximation of pi during that time, and they couldn't even multiply efficiently! :-)

	______________________

7	More on data structures and clipping

Usually, the space viewed by the eye can be defined by a very high square-based pyramid. Normally, it has no bottom, but in this case it might have one for practical purposes. Let us say that the eye has a width of vision (horizontally and vertically) of 90 degrees. (Since it is square, that angle is slightly larger when taken from one corner to the opposite.) Thus, we will see 45° to both sides, up and down. This means a slope of one. The planes bounding that volume can be defined as x<z, x>-z, y<z and y>-z. To a programmer with an eye for this, it means that it easy to «clip», i.e. remove obviously hidden polygons, with these bounds. All polygons that do not satisfy either of these inequalities can be discarded. Note that these polygons must fall totally outside the bounding volume. This is called trivial rejection.

For polygons that are partially inside, partially outside that region, you may or may not want to clip (I am not sure yet of which is more efficient - I suspect it depends on the situation, i.e. number of polygons, number of edges in each polygon...)

There is also a slight problem with polygons that are partially behind, partially in front of the eye. Let us give a little proof. All of our polygon-drawing schemes are based on the assumption that a projected line is also a line. Here is the proof of that.

Let us consider the arbitrary line L, not parallel to the plane z=0, defined parametrically as:
	x= at+b
L:	y= ct+d
	z= t

Now, let us define the projection P(x,y,z) as follows:

P(x,y,z):	u=x/z
	v=y/z

Then, P(L) becomes

P(L):	u=(at+b)/t
	v=(ct+d)/t

We are using a segment here, not an infinite line. If no point with t=0 needs to be projected, we can simplify:

P(L):	u=a+b/t
	v=c+d/t

Let us measure the slope of that projection:

du/dv=(du/dt) / (dv/dt) = (-b/t2)/(-d/t2)

Since t is not 0

du/dv=b/d, a constant, thus the parametric curve defined by P(L) is a straight line if t is not zero.

However, if t is null, the result is indeterminate. We can try to find the limits, which are pretty obvious. When t tends towards zero+, then u will tend towards positive (or negative if b is negative) infinity. If t tends towards zero-, then u will tend towards negative (or positive, if b is negative) infinity. Thus, the line escapes to one infinity and comes back from the other if t can get a zero value. This means that if the segment crosses the plane z=0, the projection is not a line per se.

However, if the line is parallel to z=0, L can be defined as follows:
	x=at+b
L:	y=ct+d
	z=k

Where k is the plane that contains L. P(L) becomes

P(L):	u=(at+b)/k
	v=(ct+d)/k

The equations do not admit a null k. For a non-null k, we can find the slope:

du/dv = (du/dt)/(dv/dt) = (a/k)/(c/k)
du/dv=a/c

Still a constant, thus P(L) is a line. However, if k is null, we do not get a projection. As k gets close to zero, we find that the projection lies entirely at infinity, except for the special case b=d=0, which is an indetermination.

This will have a practical impact on our 3d models. If a cube is centered on the eye, it will be distorted unless we somehow take that into consideration. The best way to do it, with little compromise and efficient coding, is to remove all parts of polygons that lie in the volume defined by z<=0. Of course, if an edge spans both sides of the plane z=0, you will have to cut it at z=0, excluding everything defined by z<=0. Thus, you will have to make sure the point of the edge with the smaller z coordinate is strictly greater than 0. Normally, you would have to take the next higher real number (whatever that might be!) but due to hardware contraints, we will use any arbitrarily small number, such as z=.001 or even z=.1 if you are building a real-time something in which infinite precision is not required.

Lastly, you might not want to render objects that are very far away, as they may have very little impact on the global scene, or because you are using a lookup table for logarithms and you want to make sure that all values are within range. Thus, you will set an arbitrary maximum z for all objects. It might even be a fuzzy maximum z, such as «I don't wand any object to be centered past z0». Thus, some objects may have some part farther than z0, but that may be more efficient than the actual clipping of the objects.

Moreover, you might want to make several polygon representations of the objects, with different levels of details as measured by the number of polygons, edges and vertices, and show the less detailed ones at greater distances. Obviously, you would not normally need to rasterize 10 polygons for an object that will occupy a single pixel, unless you are doing precision rendering. (But then, in that case, you'll have to change the whole algorithm...)

You may want to perform your trivial rejection clipping quickly by first rotating the handle coordinates of each polyehdra (the handle coordinate is the coordinate of its local coordinates' origin), and performing trivial rejection on them. As an example, if you determine that all of a polyhedra is contained within a sphere of radius r, and that this polyhedra is «centered» at z=-3r, then it is obvious that it is entirely outside the view volume.

Normally, the view volume is exactly one sixth of the total viewable volume (if you do not take into account the bottom of the pyramid). Thus, you will statistically eliminate five sixths of all displayable polyhedras with such a clipping algorithm.

Let's make a brief review of our speed improvements over the brute-force-transform-then-see-no-think algorithm. We multiplied the speed by six with our structure of pointers to edges that point to vertexes. Then, we multiply again by two (on average) with back-face culling, and now by six again. The overall improvement is a factor of 72! This shows you how a little thinking pays. Note however that that last number is a bit skewed because some operations do not take as long as some others, and also by the fact that most polyhedra do not contain a nearly infinite amount of vertices.

We have not taken into account several optimizations we mentioned, such as the lessening of detail, a bottom to the view volume and we did not compare with other kinds of algorithm. But another appeal of this algorithm is that it draws once, and only once to each pixel, thus allowing us to superimpose several polygons with minimal performance loss, or make complicated calculations for each pixel knowing that we will not do it twice (or more!) for some pixels. Thus, this algorithm lends itself very well to texture mapping (or more properly, pixmap mapping), shading or others.

However, we have to recognize that this approach is at best an approximation. Normally, we would have to fire projection rays over the entire area covered by a pixel (an infinite number of rays), which will usually amount to an integral, and we did not perform any kind of anti-aliasing (more on that later), nor did we filter our signal, nor did we take curved or interpenetrating surfaces into account, etc... But it lends itself very well to real-time applications.

8	A few more goodies...

Anti-aliasing is something we invented to let us decrease the contrast between very different pixels. More properly, it should be considered an approximation of high frequencies in a continuous image through intensity variations in the discrete view space. Thus, to reduce the stairway effect in a diagonal line, we might want to put some gray pixels.

Figure 7

Figure 7 demonstrates this. This diagonal line a was «anti-aliased» manually (I hear the scientists scream) in b, but it should let you understand the intention behind anti-aliasing.

Scientists use filters, and so will we, but I will not discuss fully the theory of filters. Let us see it that way. We first generate a normal image, without any kind of anti-aliasing. Then we want to anti-aliase it. What we do is take each pixel, modify it so it resembles its non-anti-aliased neighbors a little more and show that on screen.

Let us further imagine that how we want to do that specifically is to take one eight of each of the four surrounding pixel (north, south, east and west) and one half of the «real» pixel, and use the resulting color for the output pixel.

If our image is generated using a 256-color palette, there comes a handy little trick. Let us imagine that we want to take one half of color a and one half of color b and add them together. Well, there are 256 possible values for a. And for each single value of a, there are 256 values of b, for a total number of 256x256 pairs, or 64k pairs. We can precalculate the result of mixing one half a with one half b and store it in a matrix 64kb large.

Now let us call the four surrounding pixels n,s,e,w and the central point o. If we use the above matrix to mix n and s, we get a color A that is one half n plus one half s, that is A=n/2+s/2. Then, mix e and w as B=e/2+w/2. Then, mix A and B. C=A/2+B/2=n/4+s/4+e/4+w/4. Lastly, mix C and o, and we get D=o/2+n/8+s/8+e/8+w/8, the exact desired result. Now, this was all done through a lookup table, thus was very fast. It can actually be done in real time.

Thus, if you want to do this, you will first create the original, un-anti-aliased image in a memory buffer, and then copy that buffer as you anti-aliase it on screen. Obviously, it will be slower a little, but you might want to consider the increased picture quality. However, you will probably want to do this in assembler, because it can be highly optimized that way (most compilers get completely confused...)

Another very nice technique is this. For pixel (u,v) on screen, mix equal parts of pixels (u,v), (u+1,v), (u,v+1), (u+1,v+1) in your rendering buffer. It's as if the screen was half a pixel of center in relation to the rendering buffer. You'll need to render an extra row and column though.

	__________________

It is also interesting to note that the silhouette of a projected sphere is a circle. The centre of the circle is the projection of the centre of the sphere. However, the radius is NOT the projection of the radius of the sphere. However, when the distance is large enough when compared to the radius of the sphere, this approximation is good enough.

	__________________

To determine if a point is contained inside a polyhedra, we extend our test from the polygons like this: if a line drawn from the point to infinity crosses an odd number of faces, the it is contained, otherwise it is not.

	__________________

If you want to approximate a continuous spectrum with 256 colors (it will not look perfect though), use base 6. That is, using the RGB model (we could use HSV or any other), you will say that red, green and blue can all take 6 different values, for a total number of combinations of 216. If you add an extra value to green or red, you get 252 colors, close to the maximum. However, 216 is not bad and it leaves you 40 colors for system purposes, e.g. palette animations, specific color requirements for, let's say, fonts and stuff. Under Microsoft Windows, 20 colors in the palette are already reserved by the system, but you still have enough room for your 216 colors. By the way, I think it's stupid that MS Windows does not allow you to do that as a standard palette allocation. You have to define your own video driver to enable programs that limit themselves to dithering to use those colors. And that's living hell.

In addition, if the colors are arranged in a logical way, it will be easier to find the closest match of the mix of two colors as per anti-aliasing, above. (A bit like cheaper 24 bit color.)

Some people give 8 values to red and green, and 4 values to blue, for a total of 256 colors, arguing that shades of blue are perceived more difficultly. Personnally, I think the result is not as good as the base six model I previously described.

Some other nice techniques: divide your palette in spans. Example, allocate palette registers 16-48 to shades of blue. Use only registers 20-44 in your drawings, but arrange so that it goes from white-blue at color #16 through black-blue at color #48. Then, shading is a simple matter of adding something to make it darker, or subtracting to make it lighter. Make a few spans for your essential colors, and maybe keep 30 for anti-aliasing colors. These should not be used in your pixmaps. They could be used only when performing final anti-aliasing of the rendered picture. Pick these colors so to minimize visual artifacts. (E.g. you have a blue span at 16-48, and a red span at 64-100. But you don't have any purple, so if blue happens to sit next to red, anti-aliasing will be hard to achieve. Then, just a bunch of your special 30 colors for shades of purple.) Another particular trick to increase the apparent number of colors (or rather, decrease the artifacts created by too few colors), you could restrict your intensity to lower values. E.g. make your darkest black/gray at 25% luminosity, and your lightest white at 75%. This will make your graphics look washed out, but your 256 colors will be spread over a narrower spectrum.

	_______________________

9	Sorting

Sorting becomes an important issue because we do it a lot in our 3d rendering process. We will discuss here briefly two general sorting algoritmthms.

The first sort we will discuss is the bubble sort. It is called this way because the elements seem to «float» into position as the sort progresses, with «lighter» elements floating up and «heavier» ones sinking down. Given the set S, composed of elements s0, s1, ..., sn-1. This set has n elements. Let us also suppose that we can compare two elements to know which of the two is «lighter», i.e. which one goes first. Let us also suppose that we can exchange two consecutive elements in the set S. To sort the set S, we will proceed as follow. Starting at element 0 up to element n-2, we will compare each element to its successor and switch if necessary. If no switches are made, the set is sorted. Otherwise, repeat the process. In pseudocode, we have:

repeat the following:
	for i varying from 0 to n-2, repeat the following
		if si should go after si+1, then exchange si and si+1
	end of the loop
until no elements are switched in the for loop (i.e the set is sorted)

This algorithm, while very inefficient, is very easy to implement, and may be sufficient for relatively small sets (n of about 15 or less). But we need a better algoritm for greater values of n.

However, we have to keep in mind that we are sorting values that are bounded; that is, these values have a finite range, and it is indeed small, probably 16 bits or 32 bits. This will be of use later. But let us introduce first the sort. It is called the radix sort.

Assume we are give a set S of numbers, all ranging from 0 to 9999 inclusively (finite range). If we can sort all si according to the number in the ones' place, then in the tens' place, then in the hundreds' place, then thousands' space, we will get a sorted list. To do this, we will make a stack for each of the 10 possible values a digit can take, then take each number in the S set and insert them in the appropriate stack, first for the rightmost digit. Then, we will regroup the 10 stacks and start over again for the second rightmost digit, etc... This can be done in pseudocode as follows:

divisor=1
while divisor<10000, do the following:
	for i varying from 0 to n-1, do the following:
		digit=(si divided (integer version) by divisor) modulo 10
		put si on stack # digit
	end of for loop
	empty set S
	for i varying from 0 to 9, do the following:
		take stack # i and add it to set S
	end of for loop
	divisor=divisor x 10
end of while loop

Now, let us assume that we want to work with hexadecimal numbers. Thus, we will need 16 stacks. In real life, what we want to sort is edges or something like that. If we limit ourselves to 16384 edges (16k, 12 bits), then our stacks need to be able to fit 16384 objects each. We could implement them as arrays and use lots of memory. A better way is to use lists.

One could also use the quicksort algorithm, which has a complexity of O(nxln(n)), which seems worse than n for this particular problem. It is also a bit more tedious to code (in my humble opinion). In our problem, we may not even need to sort for 16 bits. The only place that we need a powerful sorting algorithm is when sorting all edges from top to bottom and from left to right. If we are certain that all coordinates are bound between, say, 0 and 200 (e.g. for y coordinate), then we can sort on 8 bits. For the x coordinate, though, depending if we already clipped or not, we may have to sort on the full 16 bits. When using the z-buffer algorithm (decribed scantily elsewhere), you want to do a z-sort, then drawing the polygons that are on top FIRST so that you can avoid these nasty pixel writes.

	______________________

10	Depth-field rendering and the Z-buffer algorithm

You can approximate any 3d function z=f(x,y) with a depth field, that is, a 2d matrix containing the z- values of f(x,y). Thus, z=M[x][y]. Of course, if you want any precision at all, it has somewhat large memory requirements. However, rendering this depth-field is fairly efficient and can give very interesting results for smooth, round surfaces (like a hill range or something).

To render a depth field, find the vector V from our eye to the pixel in the projection plane in 3space. Shoot a ray, as before, using that vector. Start at the eye and move in direction V a little. If you find that you are still above the surface (depth field), then keep going (add another V and so on) until you hit or fall under the surface. When that happens, color that pixel with the color of the surface. The color of the surface can be stored in another matrix, C[x][y]. In pseudocode, we get:

M[x][y]=depth field
C[x][y]=color field
for each pixel to be rendered, do:
	V=projection ray vector
	P=Eye coordinates
	while (M[Px][Py]<Pz)
		P=P+V
	end while
	color the current pixel with C[Px][Py]
	if we are using z-buffer algorithm, then set Z[Px][Py] to Pz (see below)
end for

Of course, you will want to break the loop when P is too far from origin (i.e. when you see that nothing is close through that projection ray, stop checking).

For vector graphics, we were forced to place our eye centered on origin and looking towards positive z values. That forced us to rotate everything else so the eye would be positioned correctly. Now, obviously, we cannot rotate the whole depth field as it would be very time- and memory-consuming. Instead, the eye now can be centered anywhere (the initial value for P need not be (0,0,0)) and the vector V can be rotated into place. Even better, if we can find the vector W from the eye to the center of the projection plane, and the two base vectors u and v in 3d for the projection plane, all of the other vectors Vu,v wil be a linear combination of the form: V=W+au+bv where (a,b) is the coordinate of the pixel in the (u,v) coodinates system. If we are doing this rendering with scan lines (as we most probably are doing), then the above equation can be done incrementally, e.g. say Va,y=A, then Va+1,y=A+u. No multiplications whatsoever.

The z coordinate of each pixel is readily available in the above algorithm and it has a use. Let us store all the z coordinate values in a matrix Z[a][b]. If we then want to add a flat image with constant z value (say z=k) to the already rendered scene, we can do it very easily by writing only to pixels for which Z[a][b]>k. Thus, we can use another rendering technique along with depth-field rendering. This is called the z-buffer algorithm. Accelerating calculations for the z-buffer with polygons, you might want to remove the divide-per-pixel. That is fully discussed towards the end of the texture-mapping chapter.

A few optimizations are easily implemented. Let us assume we are drawing from the bottom of the screen to the top, going in columns instead of rows. First, if you are not banked (or banked less than the steepest slope, e.g. no «overhangs» are seen in the picture) and you have rendered pixel a. Then you know how far pixel a is from you. Then, then next pixel up, pixel b, will be at least as far as pixel a. Thus, you do not need to fire a ray starting from your eye to pixel b, you need only to start as far as pixel a was. In fact, if you are not banked, you need only to move your Pz coordinate up a little and keep going.

Therefore, if pixel a is determined to be a background pixel, so will all the pixels above it.

You might want to render all of you polygons and other entities using a z-buffer algorithm. As an example, say you have 3 polygons, A, B and C, and a sphere S. You could tell your renderer to render poly A, which it would do right away, filling the z-buffer at the same time, then render sphere S, which it would do right waway again, and lastly polygons B and C. If you are thinking about using the z-buffer algorithm in this fashion, you might want to take a look at the last part of the texture mapping chapter, where a technique for avoiding the one divide per pixel is discussed. What you gain by using this algorithm instead of a scan line algorithm is obvious. Interpenetrating polygons behave correctly, texture mapping is a mere step away, mixing rendering techniques is straightforward. There's also a double edged optimization that you can make, which follows.

Sort your polygons in generally increasing z order (of course, some polygons will have an overlap in z, but that's everybody's problem with this algorithm - your solution is just as good as anybody else's.) Now, you're scan-converting a polygon. If you can find a span (e.g. a part of the current scanline) for which both ends are behind the same convex object, then that span is wholly hidden. As an example of this, if both enpoints of the current span (on the scanline) are behind the same convex object, then the current span is entirely hidden by that convex object. If all of your objects are convex polygons, then you don't have to check for convexity of the object. Another interesting example is this:

If current span's left endpoint is hidden
	let A be the object that hides the left endpoint (of the span)
	if A hides the right endpoint
		stop drawing the scanline
	else
		start drawing from the right end of the span, leftwards, until object A is in front of the span
	end if
else
	if current span's right endpoint is hidden
		let B be the object that hides the right endpoint
		scan rightwards until object B is in front of the span
	else
		draw span normally
	end if
end if

This can be applied also to polyhedras if they are convex. This last optimization should be relatively benign (e.g. I expect a factor of two or some other constant from it).

	_____________________

11	Bitmap scaling and mixing rendering techniques

Another popular technique is called bitmap (or more properly pixmap) scaling. Say you are looking in a parallel direction to the z axis, and you are looking at a picture with constant z. Through the projection, you will realize that all that happens to it is that it is scaled by a constant factor of 1/z. That is, if you look at a photograph lying in the plane z=2, it will look exactly the same (once projected) as the same photograph at z=1, except that it will be twice as small. You could also do that with the picture of a tree. From a great distance, the approximation will be fairly good. When you get closer, you might realize that the tree looks flat as you move. But nonetheless, it can yield impressive results. Bitmap scaling often suffices for ball-like objects. You could use it for an explosion, smoke, spherical projectile, etc... The advantage of bitmap scaling is that it has only one point, probably it's center, to transform (i.e, translate and rotate into place in 3d), and scaling it is very fast.

The advantages of such a technique are obvious. First, it is very easy to scale a pixmap by a given factor. Second, the z value is a constant, thus it can easily be incorporated in the above z-buffer algorithm with depth-field rendering for very nice results. However, it will be obvious to the careful observer that you can not smoothly go around it, as it is always facing you the same way (i.e. you cannot see the back door to a house if it is rendered this way, because, obviously, the bitmap being scaled does not have a back door).

Moreover, the depth-field rendering technique blends quite smoothly with vector graphics. First, render the depth field in a buffer B. Do a scan-line algorithm for the vector graphics, but with an added twist: the «background» bitmap is B, and use the z-buffer algorithm to determine if a pixel in B lies in front of the current polygon if any.

If you do decide to use depth-field rendering along with vector graphics, and you decide not to bank your eye, then you can eliminate one full rotation from all calculations. E.g. if your eye is not banked, you do not need to rotate around the z axis.

	_________________________

Scaling a bitmap can be done fairly easily this way. The transformation is like this: we want to take every (u,v) point and send it into (u/s,v/s), where s is the scaling factor (s=2 means that we are halving the size of the bitmap). The process is linear.

(i,j)=(u/s,v/s)
(is,js)=(u,v)

Thus, if we are at pixel (i,j) on screen, it corresponds to pixel (is,js) in the bitmap. This can also be done incrementally. If un=A, then un+1=A+s. No multiplications are involved. Thus, if we are viewing a bitmap at z=2 according to the above said 3d projections, then we can render it using a scan-line algorithm with these equations and s=z=2.

If a bitmap has holes in it, we can precalculate continuous spans in it; e.g., in an O, drawing a scanline through the center shows you that you are drawing two series of pixels. If you want the center to be transparent (allow a background to show through), then you can divide each scanline of the O in two spans, one left and one right, and render them as two separate bitmaps; e.g., break the O in ( and ) parts, both of which have no holes. Note however that we have a multiplication at each beggining of a span, so this might not be efficient. Or you could select a color number as being the «transparent» color, and do a compare per pixel and not blit if it is the transparent color. Or you can allocate an extra bitplane (a field of bit, one bit per pixel) where each bit is either 1 for opaque (draw this pixel) or 0 for transparent (do not draw this pixel).

	______________________

Partially translucid primitives can be rendered pretty easily. This is an old trick, you see. You can even scavenge the lookup table you used for anti-aliasing. If you want to have a red glass on a part of the screen, the lookup table will tell you how to redden a color. That is, the item #c in the lookup table tells you what color is color c with a little red added. The lookup table need not be calculated on the fly, of course. It can be calculated as part of your initialization or whatever.

	______________________

12	About automatically generating correctly oriented normals to planes.

This bit is not meant to be a blindingly fast algorithm for region detection. It is assumed that these calculations are made before we start generating pictures. Our objective is to generate a normal for each plane in a polyhedra oriented so that it points outwards respectively to the surface. We assume all surfaces are defined using polygons. But before determining normals for planes, we will need to examine polygons in 3d.

Each polygon being constituted of a set of edges, we will later require a vector for each edge that points in the polygon for any edge. See below.

Figure 8

In figure 8, we can see in red the said vectors for all edges. As we can see, the vectors do not necessarely point towards the general center of mass. We could use the algorithm we used for rendering polygons in two dimensions to determine the direction of the arrows, but that is a bit complicated and involves trying to find a direction that is not parallel to any of the above lines. Thus, we will try another approach.

First, pick any point in the polygon. Then, from that point, take the farthest point from it. See below:

Figure 9

(Please, do not measure it with a ruler :-) )
Now say we started with the point marked by the big red square. Now, both point a and b are the farthest points from the red square. You can pick any of the two. The objective here is to find three consecutive points around the polygon for which the inside angle is less than 180 degrees. If all three are on the circle (worst case), their angle is strictly smaller than 180 degrees (unless they are all the same point, silly).

You could also pick the topmost point. If several are at the same height, use the leftmost (or rightmost, whatever) one in the topmost points.

Now we are certain that the small angle is the inside angle. Let us draw the general situation:

Figure 10

Now, in a we see the above mentionned farthest point and the two adjacent edges. If we take these edges as vectors as in b and then sum them as in c, we find a point (pointed by the red vector in c) which is necessarely on the in side for both edges. Now it is easy to find a perpendicular vector for both edges and make sure that they are in the good general direction. A way would be to take the perpendicular vector and do a scalar multiplication with the above red vector, the result of which must be strictly positive. If it is not, the you must take the opposite vector for perpendicular vector.

Now we know the situation for at least one edge. We need some way to deduce the other edges' status from that edge's.

If we take Figure 10.c as our general case, the point at the end of the red vector is not necessarely on the in side for both edges. However, if it is on the in side for any of the two edges, it is on the in side for both, and if it is on the out side for any edge, it is on the out side for both.

Normally, we will never find two colinear edges in a polygon, but if we do, it is fairly easy to solve it anyway (left as an exercise).

Now, briefly, here is how to determine if a normal to a plane points in the right direction. Pick the point with the greatest z value (or x, or y). If many points have the same z value, break ties with greatest x value, then greatest y value. Then, take all connected edges. Make them into vectors, pointing from the above point to the other endpoints (that is, the z component is less than or equal to zero). Make them unit vectors, that is, length of one. Find the vector with the most positive (or less negative) z value. Break ties with the greates x value, and then greatest y value. Take its associated edge. Now that edge is said to be outermost.

Now, from that outermost edge, take the two associated polygons. The edge has an associated vector for each connected polygon to indicate the "in" direction. Starting from any of the two points of the edge, add both these vectors. The resulting point is inside. Adjust the normal so that it is on the right side of the plane for both polygons.

Now we have the normal for at least one plane straight. We will deduce correct side for the normal to the other polygons in a similar manner than we did for the edges of a polygon, above.

Take two adjacent polygons, one for which you do know the normal, the other for which you want to calculate the correct normal. Now take a shared edge (any edge will do). That edge has two "in" vectors for the two polygons. Starting from any of the two points of the edge, add the two «in» vectors. If the resulting point is on the «in» side for a plane, it is also on the «in» side of the other plane. If it is on the «out» side of a plane, it is on the «out» side for both planes.

To get a general picture of what I'm trying to say, here's a drawing. Basically, where two lines intersect, you get four distinct regions. One region lies entirely on the «in» side of both lines and another region lies entirely on the «out» side of both lines (if the two lines define a polygon). The two other regions are mixed region (they lie on the «in» side of a line and on the «out» side of the other line). The «mixed» regions may or may not be part of the actual polygon, depending on whether the angle is lesser or greater than 180 degrees.

Figure 11

So we are trying to find a point that is in the «in» or «out» region, but not in the «mixed» region. If, as above, we take the edges as vectors and add the together, we end up in the either the «in» region or the «out» region, but not in the «mixed» region. That can be proved very easily, geometrically (yet strictly), but I don't feel like drawing it. :) Anyway, it should be simple enough.

	__________________

13	Reducing a polygon to a mesh of triangles

Now that's fairly simple. The simplest is when you have a convex polygon. Pick any vertex. Now, take its two adjacent vertexes. That's a triangle. In fact, that's your first triangle in your mesh of triangles. Remove it from your polygon. Now you have a smaller polygon. Repeat the same operation again until you are left with nothing. Example:

Figure 12

For the above polygon, we pick the green edges in 1. The triangle is shown in 2. When we remove it from the polygon, we're left with what we see in 3.

This will take a n sided polygon and transform into n-2 triangles.

Now, if the polygon is not convex, it's a bit harder, and there's still a proof that's missing (if you know it, please mail it to me). First, we will assume that for any finite polygon that does not intersect itself, there is a vertex for which the two connected edges are convex, AND the triangle thus formed does NOT contain any other vertex besides the three that form the triangle. (That is what I haven't proved yet...) The proof exists, such a vertex is called the ear of a polygon.

To reduce a polygon to a mesh of triangles, find that vertex, form the triangle and remove it from the polygon. Repeat untill you're finished. It should be fairly easy to find the vertex. First, determine the «in» sides to all edges, as per Chapter 12. Now, find all edges for which the polygon is locally convex:

Figure 13

[The arrows point towards the «in» side]

The two possible cases are shown in figure 13, either it is or it is not convex. In 1, it is convex, as in 2 it is concave. To determine if it is or not convex, take the line that passes through a and u, above. Now, if v stands on the «in» side, it is convex. Otherwise, it is concave.

Now, of all vertex that are locally convex, pick any one for which the triangle (triangle auv in figure 13) does not contain any of the other vertexes of the polygon. A point is in a triangle if and only if it is on the «in» side of the three edges of the triangle.

	__________________

14	Gouraud shading

First, let us discuss a simple way for shading polygons. Imagine we have a light source that's infinitely far. Imagine it's at z=positive infinity, x=0, y=0. If it's not, you can rotate everything so that it is (in practice, you'll only need to rotate normals). Now, the z component of the polygons normals give us a clue as to the angle between the polygon and the light rays. The z component goes from -1, facing away from the light source (should be darkest) to 0, facing perpendicular to the light source (should be in half light/penumbra), to +1, facing the light source (should be the brightest). From that, you can shade it linearly or any way you want.

With this model, a polygon is uniformly shaded. Approximating round objects through polygons will always look rough.

Gouraud shading is what we call an interpolated shading. It is not physically correct, but it looks good. What it does is this: it makes the shade of a polygon change smoothly to match the shade of the adjacent polygon. This makes the polyhedra look smooth and curved. However, it does not change the actual edges to curves, thus the silhouette will remain unsmoothed. You can help that by allowing curved edges if you feel like it.

First, we will interpolate a vertex normal for all vertexes. That is, we will give each vertex a «normal» (note that a point does not have a normal, so you could call it pseudo-normal). To interpolate that vertex «normal», averaging all connected polygons normals would be a way. Now, find the «shade» for each vertex.

Figure 14

Now, say point a is at (ax,ay) and point b is at (bx,by). We're scan-converting the polygon, and the gray-and-red line represents the current sanline, y0. Point P is the pixel we are currently drawing, located at (x0,y0). Points m and n define the span we're currently drawing. Now, we will interpolate the color for point m using a and b, and the color for n using a and c, then interpolate the color for P using the color for m and n. Our interpolation will be linear. Say color for point a is Ca, color for point b is Cb, for point c it is Cc, Cm for m, Cn for n and CP for P.

We will say that color at point m, Cm, is as follows: 

Cm=(y0-ay)/(by-ay) x (Cb-Ca)+Ca

That is, is point m is halfway between a and b, we'll use half ca and half cb. If it's two-third the way towards b, we'll use 1/3Ca and 2/3Cb. You get the picture. Same for n:

Cn=(y0-ay)/(cy-ay) x (Cc-Ca)+Ca.

Then we do the same for CP.

CP=(x0-mx)/(nx-mx) x (Cn-Cm)+Cn.

If you think a little, these calculations are exactly the same as Bresenham's line drawing algorithm, seen previously in chapter 2.1. It is thus possible to do them incrementally. Example:

Say we start with the topmost scanline. Color for point m is at first Cb. Then, it will change linearly. When point m reaches point a, it will be color Ca. Now, say dy=by-ay, dC=Cb-Ca, and u=y0-ay.

Cm=dC/dy x u+Ca.

Then, when y0 increases by 1, u increases by one and we get

Cm'=dC/dy x (u+1)+Ca=dC/dy x u+Ca + dC/dy
Cm'=Cm+dC/dy

Same as Bresenham. Thus, when going from one scanline to the next, we simply need to add dC/dy to Cm, no multiplications or divisions are involved. The same goes for Cn. CP is done incrementally from pixel to pixel.

	____________________

Version history

Original version (no version number): original draft, chapters 1 through 8.

Version 0.10 beta: added chapters 9 through 11, added version history. Still no proofreading, spellchecking or whatever.

Version 0.11 beta: modified some chapters, noticed than (Q,P,R) is the normal vector (chapter 5). A little proofreading was done.

Version 0.12 beta: I just noticed that this document needs a revision badly. Still, I don't have time. I just added a few things about the scalar multiplication of two vectors being the D coefficient in the plane equation. Added chapter 12. Also corrected a part about using arrays instead of lists. The overhead is just as bad. Better to use lists. Ah, and I have to remember to add a part on calculating the plane equation from more than three points to reduce roundoff error.

Version 0.20 beta: Still need to add the part about calculating A,B and C in plane eq. from many points. However, corrected several of the mistakes I talked about in the preceding paragraph (i.e. revision). Ameliorated the demonstration for 2d rotations. Started work on 3d rotations (this had me thinking...) 

Version 0.21 beta: converted to another word processor. Will now save in RTF instead of WRI.

Version 0.22 beta: corrected a few goofs. 3d rotations still on the workbench. Might add a part about subdividing a polygon in triangles (a task which seems to fascinate people though I'm not certain why). Will also add the very simple Gouraud shading algorithm in the next version (and drop a line about bump mapping). This thing needs at least a table of contents, geez... I ran the grammar checker on this. Tell me if it left anything out.

Version 0.23 beta: well, I made the version number at the start of this doc correct! :-) Did chapter 13 (subdividing into triangles) and chapter 14 (Gouraud shading).

Version 0.24 beta: removed all GIFs from file because of recent CompuServe-Unisys thingy. Don't ask me for them, call CompuServe or Unisys [sigh].

Version 0.30 beta: added the very small yet powerful technique for texture mapping that eliminates the division per pixel. It can also help for z-buffering. Corrected the copyright notice.

Version 0.31 beta: I, err, added yet another bit to texture mapping [did this right after posting it to x2ftp.oulu, silly me]. Corrected a typo in chapter 13: a n-sided polygon turns into n-2 triangles, not n triangles. Was bored tonight, added something to the logs chapter (I know, it's not that useful, but hey, why would I remove it?) Added «Canada» to my address below (it seems I almost forgot internet is worldwide tee hee hee). Addendum: I'm not certain that the GIF thing (see the paragraph about v0.24beta above) applies to me, but I'm not taking any chances.

Version 0.32 beta: added bits and pieces here and there. The WWW version supposedly just became available today. Special thanks to Lasse Rasinen for converting it. Special thanks also to Antti Rasinen. I can't seem to be able to connect to the WWW server though, I just hope it really works. Corrected a goof in the free directional texture mapping section.

About the author

Sébastien Loisel is a student in university. He is studying computer sciences and mathematics and has been doing all sorts of programs, both down-to-earth and hypothetical, theoritical models. He's been thinking a lot about computer graphics lately and so decided to write this because he wanted to get his ideas straight. After a while, he decided to distribute it on Internet (why the hell am I using third person?)

The author would be happy to hear your comments, suggestions, critics, ideas or receive a postcard, money, new computer and/or hardware or whatever (grin). If you do want to get in touch with me, well here's my snail mail address:

Sébastien Loisel
1 J.K. Laflamme
Lévis, Québec
Canada
G6V 3R1

At the time of this writing, I can also be reached through e-mail at
rloisel@vm1.ulaval.ca

Addendum: see enclosed file, README.3D.

I would like to thank Lasse Rasinen making the WWW version available by converting this doc to the proper format. I would also like to thank Antti Rasinen for the WWW version. Thanks, guys.