home *** CD-ROM | disk | FTP | other *** search
/ Collection of Hack-Phreak Scene Programs / cleanhpvac.zip / cleanhpvac / FAQSYS18.ZIP / FAQS.DAT / 3D.ASC < prev    next >
Text File  |  1996-06-02  |  137KB  |  3,357 lines

  1.                                 
  2.    A tutorial for 2d and 3d vector and texture mapped graphics
  3.                           Version 0.60ß
  4.  
  5.  
  6. c 1994, 1995 by Sébastien Loisel
  7.  
  8. Acknowledgements:
  9.  
  10.      I am very thankful to the following people for providing
  11.      support and helping in the development of this document.
  12.      
  13.      Kevin Hunter for doing his best to explain a bit of linear
  14.      algebra to the rest of the world.
  15.      
  16.      Simos Hadjiyiannis for providing everyone with WWW access to
  17.      this document.
  18.      
  19.      For more information about these people, see the end of this
  20.      file.
  21.  
  22.  
  23. Note:
  24.                    This package is Shareware!
  25.  
  26. This means that you can enjoy this package freely, pass it on to
  27. your friends, and learn from it. However, if you find that this
  28. package (3d document and C source files) has been useful to you,
  29. you have to register yourself to me. Registering costs US$10.00
  30. and will allow you to use the source code in any way you like,
  31. even commercial purpose if you want to. You will be registered
  32. for life, which means that you can get any future versions of
  33. Zed3D freely and use it without feeling bad in the morning. You
  34. will be allowed to use the source code in any way you want.
  35. (Note, however, that you will have to get the Zed3D package
  36. yourself.) Sending me money will encourage me in giving you more
  37. quality products like this. Please support the shareware concept
  38. and support the shareware authors.
  39.  
  40. If you want to register, simply fill in the following
  41. registration form and print it, then mail it to my snail-mail
  42. address along with a cheque or money order for US$10.00 (no other
  43. currency accepted).
  44.  
  45.      Yes, I want to register my copy of the Zed3D document
  46.      and all future versions. As a bonus, I will be able the
  47.      Zed3D source code in any way I please. Enclosed, you
  48.      will find my US$10 contribution to the shareware
  49.      concept.
  50.      
  51.      My name is: __________________________________________
  52.      
  53.      My (snail mail) address is:
  54.      
  55.           _____________________________________________
  56.      
  57.           _____________________________________________
  58.      
  59.           _____________________________________________
  60.      
  61.      
  62.      My e-mail address is: _________________________________
  63.      
  64.      
  65.      I got Zed3D from: ____________________________________
  66.      
  67.      
  68.      I would like to see this in future versions of Zed3D:
  69.      
  70.      __________________________________________________
  71.      
  72.      __________________________________________________
  73.      
  74.      __________________________________________________
  75.      
  76.      __________________________________________________
  77.      
  78.      
  79.      Signature: __________________________________________
  80.      
  81.      Mail to:       Sebastien Loisel / Zed3D registration
  82.                     1 J.K. Laflamme
  83.                     Levis, Quebec, Canada
  84.                     Postal Code G6V 3R1
  85.      Cheque should go to Sebastien Loisel
  86. Disclaimer & copyright notice:
  87.  
  88. This document is provided "as is", without guaranty of ANY kind.
  89. This document is copyrighted 1994, 1995 by Sébastien Loisel. It
  90. may be distributed freely as long as no fee except shipping and
  91. handling is charged for it. Including this document in a
  92. commercial package requires the written permission of the author,
  93. Sébastien Loisel. The term "document" refers both to this
  94. document and the accompanying source code.
  95.  
  96. If you want to reach me, see the end of this file.
  97.  
  98.  
  99.      ____________________
  100.  
  101.  
  102.  
  103. 1    An introduction to 3d transformations
  104.  
  105.  
  106.  
  107. First, let's introduce the idea of projections. We have to
  108. realize that we are trying to view a 3d world through a 2d
  109. «window». Thus, one full dimension is lost and it would be absurd
  110. ot think that we can reconstitute it wholly and recognizably with
  111. one less dimension. However, there are many forms of depth-cueing
  112. that we introduce in a 2d image to help us recognize that 3d
  113. dimension. Examples include but are not limited to reduced
  114. lighting or visibility at a distance, reduced size at a distance,
  115. out-of-focus blur, etc...
  116.  
  117. Our media is a screen, ideally flat, and the 3d world being
  118. represented is supposed to look as real as possible. Let's look
  119. at how the eyes see things:
  120.  
  121.  
  122.  
  123. Figure 1
  124.  
  125.  
  126. The projection rays are the beams of light that reflect off the
  127. objects and head to your eye. Macroscopically, they are straight
  128. lines. If you're projecting these beams on a plane, you have to
  129. find the intersection of these beams with the plane.
  130.  
  131.  
  132.  
  133. Figure 2
  134.  
  135. The projection points are thus the intersection of the projection
  136. rays with the projection plane. Of course, projection points
  137. coming from projection lines originating from behind an object
  138. are obviously hidden, and thus comes the loss of detail from 3d
  139. to 2d. (i.e. closer objects can hide objects that lie further
  140. away from the eye).
  141.  
  142. The objects can be defined as a set of (x,y,z) points in 3d space
  143. (or 3space). Thus, a filled sphere can be defined as
  144. x**2+y**2+z**2<=1. But we need to simplify the problem right away
  145. because of machine limitations. We can safely assume the eye is
  146. always outside any object, thus we do not need do define the
  147. interior. [Note: a**b stands for «a raised to the bth power»,
  148. thus, x**2 is x squared]
  149.  
  150. Now let us define a projection onto a plane. Many different
  151. projections could be defined, but the one that interests us is
  152. the following for it's simplicity. Projecting takes any given
  153. point (x,y,z) and localizes it on an (u,v) coordinate system.
  154. I.e. space (x,y,z) is mapped onto a plane (u,v). Let's transform
  155. (x,y,z) first. Multiplying (x,y,z) by a constant keeps it on a
  156. line that passes through (0,0,0) (the origin). We will define
  157. that line as the projection ray for our projection. Since all of
  158. these lines pass through the origin, we thus define the origin as
  159. our eye. Let the constant be k. Thus our transformation as of now
  160. is k(x,y,z). Now we want to set everything into a plane. If we
  161. select k as being 1/z, assuming z is nonzero, we get a projected
  162. point of (1/z)(x,y,z). This becomes (x/z, y/z, 1). Thus all
  163. points of any (x,y,z) value except z=0 will be projected with a
  164. value z'=z/z of one. Thus all points lie in the plane z=1. This
  165. projection thus fits our need for a projection with linear
  166. projector through the eye on a plane surface.
  167.  
  168. If you didn't get all that right away, you can always look it up
  169. later. The important thing to remember now is this: to project
  170. from (x,y,z) onto (u,v) we decide to use
  171.  
  172.      u=x/z
  173.      v=y/z
  174.  
  175.  
  176. But, I hear you say, this means the eye is facing parallel to the
  177. z axis, centered on the origin and is not banked. What if that is
  178. not what I desire?
  179.  
  180. First, it is obvious that if the eye is at location (a,b,c) it
  181. will not affect what one would see to move EVERYTING including
  182. the eye by (-a,-b,-c). Then the eye will be centered.
  183.  
  184. Second, and almost as obvious, if the eye is rotated around the
  185. x,y and z axis by agles of rx, ry and rz respectively, you can
  186. rotate EVERYTHING by the exact opposite angles, which will orient
  187. the eye correctly without changing what the eye would see.
  188.  
  189. Translation (moving) is a simple addition. I.e. (x,y,z)+(a,b,c)
  190. becomes (x+a,y+b,z+c). Additions are fairly quick. But there is
  191. the problem of the rotations.
  192.  
  193.  
  194.  
  195.  
  196. 1.1 Trigonometry
  197.  
  198.  
  199. We will do it in two dimensions first.
  200.  
  201.  
  202.  
  203. Figure 3
  204.  
  205.  
  206. We have the point defined by (x,y) and we want to rotate it with
  207. a counterclockwise angle of b to (x',y'). See figure 3. The
  208. radius r stays the same throughout the rotation, as per
  209. definition, rotation changes the angle but not distance to
  210. origin. Let r be the radius (not shown on figure 3). We need to
  211. express x' and y' as functions of what we know, that is x, y and
  212. b. The following can be said:
  213.  
  214. y'=sin(a+b)r
  215. x'=cos(a+b)r
  216.  
  217. With the identities sin(a+b)=sin(a)cos(b)+sin(b)cos(a) and
  218. cos(a+b)=cos(a)cos(b)-sin(a)sin(b), we substitute.
  219.  
  220. y'=rsin(a)cos(b)+rcos(a)sin(b)
  221. x'=rcos(a)cos(b)-rsin(a)sin(b)
  222.  
  223. But from figure 3 we know that
  224.  
  225. rsin(a)=y and
  226. rcos(a)=x
  227.  
  228. We now substitute:
  229.  
  230. y'=ycos(b)+xsin(b)
  231. x'=xcos(b)-ysin(b)
  232.  
  233. A special note: you do not normally calculate sinuses and
  234. cosinuses in real time. Normally, you build a lookup table, that
  235. is, an array. This array contains the precalculated value of
  236. sin(x) for all values of x within a certain range. That is, if
  237. the smallest possible angle in your model is on half of a degree,
  238. then you need a sinus table containing 720 entries, of which each
  239. entry is the value of the sinus of the corresponding angle. You
  240. can even note that a sinus is highly symmetric, that is, you need
  241. the values of sin(x) only for x between 0 and 90 degrees.
  242. Furthermore, you do not need a cosine table, as cos(x)=sin(x+90
  243. degrees). This and fixed point mathematics may become obsolete as
  244. floating point processors become faster and faster. A good
  245. algorithm for calculating a sinus can be very quick today on a
  246. reasonably high performance platform, and it shouldn't use
  247. powers.
  248.  
  249. With that said, we will approach briefly 3d transformations.
  250.  
  251. In 3d, all rotations can be performed as a combination of 3
  252. rotations, one about the x axis, one about the y axis and one
  253. about the z axis. If you are rotating about the z axis, use the
  254. exact same equation as above, and leave the z coordinate
  255. untouched. You can extrapolate the two other transformations from
  256. that.
  257.  
  258. Note that rotating first about the x axis and then about the y
  259. axis is not the same as rotating about the y axis and then the x
  260. axis. Consider this: point (0,0,1) rotated 90° about the x axis
  261. becomes (0,1,0). Then, rotated 90° about the y axis it remains
  262. there. Now, rotating first about the y axis yields (1,0,0), which
  263. rotated then about the x axis stays at (1,0,0), thus the order of
  264. rotations is important [(0,1,0) and (1,0,0) are not the same].
  265.  
  266. Now, if you are facing west (x axis, towards negative infinity).
  267. The user pulls the joystick up. You want the ship to dive down as
  268. a result of this. You will have to rotate the ship in the xy
  269. plane, which is the same as rotating about the z axis. Things get
  270. more complicated when the ship is facing an arbitrary angle.
  271.  
  272. Here I include a document authored by Kevin Hunter. Since he did
  273. a very nice job of it, I didn't see a reason to re-do it,
  274. especially when considering the fact that it would probably not
  275. have been as tidy. :-) So, from here to chapter two has been
  276. entirely written by Kevin Hunter.
  277.  
  278.  
  279.  
  280.  
  281.      _____________________
  282.  
  283.  
  284.  
  285.  
  286. Unit Vector Coordinate Translations
  287. By Kevin Hunter
  288.  
  289.  
  290. One of the more common means of doing "world to eye" coordinate
  291. transformations involves maintaining a system of "eye" unit
  292. vectors, and using a property of the vector dot product.  The
  293. advantage of this approach, if your "eye" frame of reference is
  294. rotated with respect to the "world" frame of reference, is that
  295. you can do all the rotational calculation as simple projections
  296. to the unit vectors.
  297.  
  298. The dot product of two vectors yields a scalar quantity:
  299.  
  300.      R dot S = |R| |S| cosa
  301.  
  302. where |R| and |S| are the lengths of the vectors R and S,
  303. respectively, and a is the angle between the two vectors.  If S
  304. is a unit vector (i.e. a vector whose length is 1.0), this
  305. reduces to
  306.  
  307.      R dot S = |R| cosa
  308.  
  309. However, a bit of geometry shows that |R| cosa is the portion of
  310. the R vector that is parallel to the S vector.  (|R| sina is the
  311. part perpendicular).  This means that if R is a vector with its
  312. origin at the eye position, then R can be converted from "world"
  313. space into "eye" space by calculating the dot product of R with
  314. the eye X, Y, and Z unit vectors.  Getting a "eye origin" vector
  315. is trivial - all you need to do is subtract the "world"
  316. coordinates for the eye point from the "world" coordinates for
  317. the point to be transformed.  If all this sounds a bit foggy,
  318. here's a more concrete example, along with a picture.  To make
  319. the diagram easier, we'll consider the two-dimensional case
  320. first.
  321.  
  322.                                 
  323. Figure 3a
  324.  
  325. In this figure, the "world" coordinates have been shown in their
  326. typical orientation.  The "eye" is displaced up and left, and
  327. rotated to the right.  Thus, the EyeX vector would have world-
  328. space coordinates something like (0.866, -0.500) [I've picked a
  329. 30 degree rotation just for the numbers, although my less-than-
  330. perfect drawing talents don't make it look like it).  The EyeY
  331. vector, which is perpendicular to the EyeX vector, has world-
  332. space coordinates (0.500, 0.866).  The eye origin point is at
  333. something like (1, 2).
  334. Based on these numbers, any point in world space can be turned
  335. into eye-relative coordinates by the following calculations:
  336.  
  337.      EyeSpace.x   = (Px - 1.0, Py - 2.0) dot (0.866, -0.500)
  338.              = (0.866Px - 0.866) + (-0.500Py + 1.00)
  339.              = 0.866Px - 0.500Py + 0.134
  340.      EyeSpace.y   = (Px - 1.0, Py - 2.0) dot (0.500, 0.866)
  341.              = (0.500Px - 0.500) + (0.866Py - 1.732)
  342.              = 0.500Px + 0.866Py - 2.232
  343.  
  344. The simple way of implementing this transformation is a two-step
  345. approach.  First, translate the point by subtracting the eye
  346. origin.  This gives you the eye-relative vector.  Then, perform
  347. each of the three dot products to calculate the resulting X, Y,
  348. and Z coordinates in the eye space.  Doing this takes the
  349. following operations for each point:
  350.  
  351.      1.Three subtractions to translate the point
  352.      2.Three multiplications and two additions for the X
  353.        coordinate
  354.      3.Three multiplications and two additions for the Y
  355.        coordinate
  356.      4.Three multiplications and two additions for the Z
  357.        coordinate
  358.  
  359. Since additions and subtractions cost the same in virtually any
  360. system, this reduces to a total of nine multiplications and nine
  361. additions to translate a point.
  362.  
  363. If you're striving for performance, however, you can make this
  364. run a little better.  Note, in the two-dimensional example above,
  365. that we were able to combine two constants that fell out of the
  366. second line into a single constant in the third line.  That last
  367. constant turns out to be the dot product of the particular eye
  368. coordinate vector and the negative of the eye origin point.  If
  369. you're smart, you can evaluate this constant once at the
  370. beginning of a set of point transformations and use it over and
  371. over again.
  372.  
  373. This doesn't change the total number of operations - the three
  374. subtractions you take away at the beginning you put back later.
  375. The advantage, however, is that the constant is known ahead of
  376. time, which may let you keep it in a register or cache area.
  377. Basically, you've reduced the number of constants needed to
  378. transform the coordinates and the number of intermediate values,
  379. which may let you keep things in registers or the cache area.
  380. For example, this approach lets you use a stacked-based FPU in
  381. the following manner:
  382.  
  383.      push P.x
  384.      push UnitVector.x
  385.      mult
  386.      push P.y
  387.      push UnitVector.y
  388.      mult
  389.      add
  390.      push P.z
  391.      push UnitVector.z
  392.      mult
  393.      add
  394.      push CombinedConstant
  395.      add
  396.      pop  result
  397.  
  398. This algorithm never has more than three floating point constants
  399. on the FPU stack at a time, and doesn't need any of the
  400. intermediate values to be removed from the FPU.  This type of
  401. approach can cut down on FPU-to-memory transfers, which take
  402. time.
  403.  
  404. Of course, if you don't have a stack-based FPU, this second
  405. approach may or may not help.  The best thing to do is to try
  406. both algorithms and see which works better for you.
  407.  
  408.  
  409.  
  410. Maintaining Systems of Unit Vectors
  411.  
  412. Spinning In Space
  413. To effectively use the system discussed in the previous section,
  414. you need to be able to generate and maintain a set of eye-space
  415. unit vectors.  Generating the vectors is almost never a problem.
  416. If nothing else, you can use the world X, Y and Z vectors to
  417. start.  Life gets more interesting, however, when you try to
  418. maintain these vectors over time.
  419.  
  420. Translations of the eye point are never a problem.  Translations
  421. are implemented by adding or subtracting values to or from the
  422. eye origin point.  Translations don't affect the unit vectors.
  423. Rotations, however, are another matter.  This kind of makes
  424. sense, since handling rotations is what maintaining this set of
  425. vectors is all about.
  426.  
  427. To be effective, the set of three unit vectors have to maintain
  428. two basic properties.  First, they have to stay unit vectors.
  429. This means that we can't let their lengths be anything other than
  430. 1.0.  Second, the three vectors have to stay mutually
  431. perpendicular.
  432.  
  433. Ensuring that a vector is a unit vector is very straight-forward
  434. conceptually.  All you need to do is to calculate the length of
  435. the vector and then divide each of its individual coordinates by
  436. that length.  The length, in turn, can be easily calculated as
  437. SQRT(x^2 + y^2 + z^2).  The only interesting part of this is the
  438. SQRT() function.  We'll come back to that one later.
  439.  
  440. Keeping the set of vectors perpendicular takes a bit more work.
  441. If you restrict rotations of the vector set to sequences of
  442. rotations about the individual vectors, it's not too hard to do a
  443. decent job of this.  For example, if you assume that, in eye
  444. space, Z is to the front, X is to the right and Y is down (a
  445. common convention since it lets you use X and Y for screen
  446. coordinates in a 0,0 == top left system), rotation of the eye to
  447. the right involves rotating around the Y axis, with Z moving
  448. toward X, and X toward -Z.  The calculations here are very
  449. simple:
  450.  
  451.      Znew = Zoldcosa + Xoldsina
  452.      Xnew = Xoldcosa + (-Zold)sina
  453.  
  454. Of course, while you're doing this you need to save at least the
  455. original Z vector so that you don't mess up the X calculation by
  456. overwriting the Z vector before you get a chance to use it in the
  457. second equation, but this is a detail.
  458.  
  459. In a perfect world, when you're done with this rotation, X, Y,
  460. and Z are still all unit vectors and still mutually
  461. perpendicular.  Unfortunately, there's this thing called round-
  462. off error that creeps into calculations eventually.  With good
  463. floating point implementations, this won't show up quickly.  If
  464. you're having to do fixed-point approximations (for speed, I
  465. assume) it will show up a lot faster, since you're typically
  466. storing many fewer mantissa bits than standard FP formats.
  467.  
  468. To get things back square again there are a couple of things you
  469. can do.  The most straight forward is to use the cross product.
  470. The cross product of two vectors yields a third vector which is
  471. guaranteed (within your round-off error) to be perpendicular to
  472. both of the input vectors.  The cross product is typically
  473. represented in determinant form:
  474.                    x   y   z
  475.          A cross   Ax  Ay  Az
  476.          B =
  477.                    Bx  By  Bz
  478.  
  479. In more conventional format, the vector resulting from A cross B
  480. is
  481.  
  482.      (AyBz - AzBy, AzBx - AxBz, AxBy - AyBx)
  483.  
  484. One thing about the cross product is that you have to watch out
  485. for the sign of things and the right-handedness or left-
  486. handedness of your coordinate system.  A cross B is not equal to
  487. B cross A - it's equal to a vector in exactly the opposite
  488. direction.  In a right-handed coordinate system, X cross Y is Z,
  489. Y cross Z is X, and Z cross X is Y.  In a left-handed coordinate
  490. system, the arguments of the crosses have to be exchanged.
  491.  
  492. One tactic, then, is to rotate one of the vectors, and then,
  493. instead of rotating the other, to recalculate the other using the
  494. cross product.  That guarantees that the third vector is
  495. perpendicular to the first two.  It does not, of course,
  496. guarantee that the vector that we rotated stayed perpendicular to
  497. the stationary vector.  To handle this, we could then adjust
  498. either the stationary vector or the one we originally rotated by
  499. another cross product if we wanted to.  In practice, this
  500. frequently isn't necessary, because we rarely are rotating only
  501. about one axis.  If, for example, we pitch and then roll, we can
  502. combine rotations and crosses so that no vector has undergone
  503. more than one or two transformations since it was last readjusted
  504. via a cross product.  This keeps us from building up too much
  505. error as we go along.
  506.  
  507. One other thing to remember about the cross product.  The length
  508. of the generated vector is the product of the lengths of the two
  509. input vectors.  This means that, over time, your calculated
  510. "cross" vectors will tend to drift away from "unitness" and will
  511. need to be readjusted in length.
  512.  
  513. Faster!
  514. It is unusual in interactive graphics for absolute and unwavering
  515. precision to be our goal.  More often (much more often) we're
  516. willing to give up a tiny bit of perfection for speed.  This is
  517. particularly the case with relatively low-resolution screens,
  518. where a slight error in scale will be lost in pixel rounding.  As
  519. a result, we'd like a measure of how often we need to readjust
  520. the perpendicularity of our unit vectors.  As it happens, our
  521. friend the dot product can help us out.
  522.  
  523. Remember that we said that the dot product of two vectors
  524. involves the cosine of the angle between them.  If our unit
  525. vectors are actually perpendicular, the dot product of any pair
  526. of them should then be zero, since the cosine of a right angle is
  527. zero.  This lets us take the approach of measuring just how badly
  528. things have gotten out of alignment, and correcting it only if a
  529. tolerable threshold has been exceeded.  If we take this approach,
  530. we "spend" three multiplications and two additions, hoping to
  531. save six multiplications and three additions.  If we're within
  532. our tolerance more often than not, we're in good shape.
  533.  
  534. Foley et al point out that we can make this tradeoff work even
  535. better, however.  Remember that the dot product represents the
  536. projection of one vector along another.  Thus, if we find that X
  537. dot Y is not zero, the resulting value is the projection of the X
  538. vector along the Y vector (or the Y vector along the X vector, if
  539. you want).  But if we subtract from the X vector a vector
  540. parallel to the Y vector whose length is this projection, we end
  541. up with X being perpendicular to Y again.  The figure may help to
  542. make this clearer
  543.                                 
  544. Figure 3b
  545.  
  546. This lets use the following algorithm:
  547.      A = X dot Y
  548.      if (abs(A) > threshold)
  549.      {
  550.           X -= Y times A
  551.      }
  552. This approach uses three multiples and two adds if we don't need
  553. to adjust, and six multiplies and five adds if we do.  This
  554. represents a savings of three multiplications over the "test and
  555. then cross product" approach.  Note that, like the cross product
  556. approach, this technique will tend to change the length of the
  557. vector being adjusted, so vectors will need to be "reunitized"
  558. periodically.
  559.  
  560. Perfection?
  561. If we want to use our unit vector system, but we need to model
  562. rotations that aren't around one of our unit axes, life gets more
  563. difficult.  Basically, except in a very few circumstances (such
  564. as if the rotation is around one of the world axes) the
  565. calculations are much harder.  In this case, we have at least two
  566. options.  The first is to approximate the rotation by incremental
  567. rotations that we can handle.  This has to be done carefully,
  568. however, because rotations are not interchangeable - rotating
  569. about X then Y doesn't give the same answer as rotating around Y
  570. then X unless the angles are infinitely small.
  571.  
  572. The second option is to use quaternions to do our rotations.
  573. Quaternions are a very different way of representing rotations
  574. and rotational orientations in space.  They use four coordinates
  575. rather than three, in order to avoid singularities, and can model
  576. arbitrary rotations very well.  Unfortunately, I'm not completely
  577. up on them at this point, so we'll have to defer discussion of
  578. them until some future issue when I've had a chance to do my
  579. studying.
  580.  
  581.  
  582.  
  583. Efficient Computation of Square Roots
  584.  
  585. We mentioned earlier that we'd need to look at square roots to
  586. "unitize" vectors.
  587.  
  588. Let's assume for the moment that we need to calculate the square
  589. root of a number, and we don't have the hardware to do it for us.
  590. Let us assume further than multiplications and divisions are
  591. cheap, comparitively speaking.  How can we go about calculating
  592. the square root of a number?
  593.  
  594. [Note:  Intel-heads like me think this is an unnatural situation,
  595. but it does, in fact, occur out there.  For example, the Power PC
  596. can do FP add/subtract/multiply/divide in a single cycle, but
  597. lacks sqrt() or trig function support.]
  598.  
  599. Isaac Strikes Again
  600. Perhaps the simplest approach is to use Newton's method.
  601. Newton's method is a technique for successively approximating
  602. roots to an equation based on the derivative of that equation.
  603. To calculate the square root of N, we use the polynomial equation
  604. x2 - N = 0.  Assuming that we have a guess at the square root,
  605. which we'll call xi, then next approximation can be calculated by
  606. the formula
  607.  
  608.      xi+1 = xi + (N - xi2) / (2xi)
  609.  
  610. This begs two questions.  First, how do we get started?  Second,
  611. when do we get there?  To answer the first part, Newton's method
  612. works really well for a well-behaved function like the square
  613. root function.  In fact, it will eventually converge to the
  614. correct answer no matter what initial value we pick.  To answer
  615. the second part, "Well, it all depends on where you start."
  616.  
  617. To get an idea of how this iteration performs, let's expand it a
  618. bit.  Let's assume that our current guess (xi) is off by a value
  619. "e" (in other words, the square root of N is actually xi + e.)
  620. In this case,
  621.  
  622.      xi+1    = xi + ((xi + e)2 - xi2) / (2xi)
  623.         = xi + (xi2 +2exi + e2 - xi2) / (2xi)
  624.         = xi + (2exi + e2) / (2xi)
  625.         = xi + e + (e2 / 2xi)
  626.  
  627. But since "xi + e" is the answer we were looking for, this says
  628. that the error of the next iteration is "e2/2xi".  If the error
  629. is small, the operation converges very quickly, because the
  630. number of accurate bits doubles with each iteration.  If the
  631. error is not small, however, it's not so obvious what happens.
  632. I'll spare you the math, and simply tell you that if the error is
  633. large, the error typically halves each iteration (except for the
  634. first iteration, in which the error can actually increase).  So,
  635. while a bad initial guess won't prevent us from finding an
  636. answer, it can slow us down a lot.
  637.  
  638. So how do we find a good initial guess?  Well, in many cases, we
  639. have a good idea what the value ought to be.  If we're re-
  640. unitizing a vector, we suspect that the vector we're working with
  641. is already somewhere near a length of 1.0, so an initial guess of
  642. 1.0 will typically work very well.  In this situation, four
  643. iterations or so usually gets you a value as accurate as you
  644. need.
  645.  
  646. What if we're trying to convert a random vector into a unit
  647. vector, or just taking the square root of a number out of the
  648. blue?  In this case, we can't be sure that 1.0 will be a very
  649. good guess.  Here are a few things you can do:
  650.  
  651.      1.If you have easy access to the floating point number's
  652.        exponent, halve it.  The square root of A x 2M is SQRT(M)
  653.        x 2M/2.  Using A x 2M/2 is a really good first
  654.        approximation.
  655.      2.If you are taking the length of a vector (x,y,z),
  656.        abs(x)+abs(y)+abs(z) won't be off by more than a factor
  657.        of 2 - also a good start.
  658.      3.It turns out that the operation will converge faster if
  659.        the initial guess is higher than the real answer.  If all
  660.        else fails, using N/2 will converge pretty quickly.
  661.  
  662. Don't go crazy trying to get a good first guess.  Things will
  663. converge fast enough even if you're somewhat off that a fast
  664. guess and an extra iteration may be faster than taking extra time
  665. to refine your guess.
  666.  
  667. If You Don't Like Newton...
  668.  
  669. If you're looking for a really accurate answer, Newton's method
  670. may not be for you.  First, you can't, in general, just iterate a
  671. fixed number of times.  Rather, you have to check to see if the
  672. last two iterations differ by less than a threshold, which adds
  673. calculation time to each iteration.  If you're in fixed-point
  674. land, the assumption we made earlier that multiplications were
  675. cheap may not be the case.  In this case, there is an alternate
  676. closed-form algorithm you can use that will result in the "exact"
  677. answer in a predefined amount of time.  This algorithm also uses
  678. an incremental approach to the calculation.
  679.  
  680. Since (x + a)2 = x2 + 2xa + a2, we can iterate from x2 to (x+a)2
  681. by adding 2xa + a2 to our working sum.  If a is2N, however, this
  682. translates into adding 2N+1 x + 22N.  Since we can multiply by 2N
  683. by shifting to the left by N bits, we can take a square root by
  684. use of shifts and subtractions.  To take the square root of a 32-
  685. bit number, for example, this results in an algorithm that looks
  686. something like this:
  687.  
  688.      remainder = valueSquared;
  689.      answer = 0;
  690.      for (iteration = N; iteration >= 0; iteration--)
  691.      {
  692.           if (remainder >= (answer << iteration+1) + (1 <<
  693.      iteration * 2))
  694.           {
  695.                remainder -= (answer << iteration+1) + (1 <<
  696.      iteration * 2);
  697.                answer += 1 << iteration;
  698.           }
  699.      }
  700.  
  701. A "real" algorithm won't repeatedly perform the shifts on the 5th
  702. line.  Rather, an incremental approach can be taken in which the
  703. answer and "1 bit" is kept in pre-shifted form and shifted to the
  704. right as the algorithm proceeds.  With appropriate temporary
  705. variable space, this algorithm shifts two values and does two
  706. subtractions and either one or zero add or "or" per iteration,
  707. and uses one iteration for each two bits in the original value.
  708.  
  709. This algorithm can be used very effectively on integer values,
  710. and can also be adapted to take the square root of the mantissa
  711. portion of a floating-point number.  In the latter case, care has
  712. to be taken to adjust the number into a form in which the
  713. exponent is even before the mantissa square root is taken, or the
  714. answer won't be correct.
  715.  
  716. Table It
  717. If you don't need absolute precision, there is a table-driven
  718. square root algorithm that you can use as well.  The code for the
  719. 32-bit FP format ("float") is in Graphics Gems, and the 64-bit FP
  720. format ("double") is in Graphics Gems III.  The source code is
  721. online at ftp:://Princeton.edu/pub/Graphics/GraphicsGems.  (Note
  722. the capital `G'.  There is also a /pub/graphics (small `g') which
  723. is not what you want...)
  724.  
  725. The way this algorithm works is to reach inside the IEEE floating
  726. point format.  It checks the exponent to see if it is odd or
  727. even, modifies the mantissa appropriately for odd exponents, and
  728. then halves the exponent.  The mantissa is then converted to a
  729. table lookup with a shift and mask operation, and the new
  730. mantissa found in a table.  The version in Graphics Gems III is
  731. set up to allow you to adjust the size of the lookup table, and
  732. thus the accuracy of the operation.
  733.  
  734. Whether or not this method is appropriate for you depends on your
  735. tolerance for accuracy vs. table size tradeoffs.  The algorithm
  736. runs in constant time, like the previous one, which is an
  737. advantage in many situations.  It's also possible to combine
  738. them, using a table lookup to get the first N bit guess, and then
  739. finishing it out using the bit shifting approach above.
  740.  
  741. Trig Functions
  742. Unfortunately, trig functions like sine and cosine don't lend
  743. themselves nicely to methods like the square root algorithms
  744. above.  One of the reasons for this is that sine and cosine are
  745. interrelated.  If you try to use Newton's method to find a
  746. cosine, you need to know the sine for your guess.  But to figure
  747. out the sine, you need cosine.  Department of Redundancy
  748. Department.
  749.  
  750. One alternate approach you can take is to use a Taylor series.  A
  751. Taylor series is a polynomial that will approximate a sine or
  752. cosine to any desired accuracy.  The problem with this is that
  753. the Taylor series for sine and cosine don't converge really
  754. quickly.  They don't do too bad, but you still have to do a bunch
  755. of computation.
  756.  
  757. A better approach in this case is probably to use a combination
  758. of computation and a look-up table, like the last square root
  759. algorithm.  Since sine and cosine are cyclic, and interrelated,
  760. you only really need to be able to calculate one of them over the
  761. range 0 - p/2.  Given this, the rest of the values can be
  762. determined.  Depending on how much calculation you want to do,
  763. how much table space you want to spend, and how accurate your
  764. answer needs to be, you can break up the region above into a
  765. fixed number of regions, and use some form of interpolation
  766. between points.  If you're in Floating Point land, over a
  767. relatively short region, you can fit a quadratic or cubic
  768. function to either the sine or cosine function with pretty good
  769. accuracy.  This gives you a second or third order polynomial to
  770. evaluate and a little bit of table lookup overhead.  For all but
  771. the lowest accuracy cases, this will probably net out faster than
  772. calculating the number of Taylor terms you'd need.
  773.  
  774. If you're in Fixed Point Land, you may be better off with a table
  775. lookup and linear interpolation, rather than using a higher-order
  776. polynomial.  For comparable accuracy, this approach will probably
  777. need more table space, but if you're avoiding floating point,
  778. you'll probably be avoiding integer multiplications as well.
  779. (Particularly because you have the decimal point adjustments to
  780. worry about).  The linear interpolation can be reduced to a
  781. single multiply by storing delta information in the table, so all
  782. you have to do is multiply the bits you stripped off while
  783. building the table index (the "fractional part" by this delta
  784. value and add it to the table entry.
  785.  
  786. One trick you can use in either case to make your table lookups
  787. easier is to measure your angles in 216 units rather than radians
  788. or degrees. 216 units basically consist of the following mapping:
  789.      Degrees Radians     216 units
  790.      0         0         0
  791.      45        p/4       0x2000
  792.      90        p/2       0x4000
  793.      180       p         0x8000
  794.      270       3p/2      0xC000
  795.      359.999   2p-e      0xFFFF
  796. The neat thing about this set of units is that they wrap just the
  797. way that radians do (and degrees are supposed to).  They give you
  798. resolution of about 0.005 degrees, and they make table lookups
  799. very easy, since you can shift-and-mask them down into table
  800. indices very easily.
  801.  
  802. Interestingly enough, this is something I'd devised for myself
  803. late last year (feeling very proud of myself) only to see the
  804. technique published in Embedded Systems a month or so later.
  805. Just goes to show you...
  806.  
  807.  
  808.  
  809.  
  810.  
  811.      ____________________
  812.  
  813.  
  814.  
  815.  
  816.  
  817. 2    Polygon drawing on a raster display in two dimensions
  818.  
  819.  
  820. The second main problem is to display a filled polygon on a
  821. screen. We will discuss that problem here. First, let us approach
  822. the 2d polygon.
  823.  
  824. We will define a polygon as a set of points linked by segments,
  825. with a defined inner region and outer region. The inner region
  826. need not be in a single piece.
  827.  
  828. Our rule for determining if a point is in a given region is to
  829. draw a line from infinity to that point. If that line crosses an
  830. odd number of segments (edges), the point is within the polygon.
  831. Otherwise it is not.
  832.  
  833.  
  834.  
  835. Figure 4
  836.  
  837.  
  838. Drawing a line from infinity to any of the gray areas of the star
  839. will cross an odd number of edges. In the white areas, it will
  840. cross an even number of edges. Try it.
  841.  
  842.  
  843.  
  844. Figure 5
  845.  
  846. The polygon shown in a, b and c are all the same. In a, the inner
  847. region has been shaded, but this has been omitted in b and c for
  848. clarity.
  849.  
  850. The polygon shown in a can be listed as vectors, as in b or c.
  851. Let us agree that both versions are equally satisfying. In b we
  852. enumerated the edges in a counterclockwise continuous fashion,
  853. meaning that the head of each vector points to the tail of
  854. another. In c, we made each vector point generally downward, or
  855. more strictly, each vector's topmost point is its tail. These
  856. vector are noted (a,b)-(c,d) where (a,b) is the location of the
  857. tail of the vector, and (c,d) is the location of the head (arrow)
  858. of the vector. In diagram c, for all vectors d<b (or b>d, either
  859. way). (Note, we will see what to do with b=d later. That being a
  860. special case, we will not discuss it now).
  861.  
  862. As we are drawing on a screen (supposedly), it would be
  863. interesting to be able to color each pixel once and only once,
  864. going in horizontal or vertical lines (as these are easy to
  865. compute). We will choose horizontal lines for hardware reasons.
  866. Thus, we start scanning from infinitely upwards to infinitely
  867. downwards using horizontal scanlines (e.g., the «first» scanline
  868. is way up, and the last is way down). Now, our screen being
  869. finite, we need only to draw those lines that are on screen. At
  870. each pixel, we will draw an imaginary line from the infinite to
  871. the left (horizontal line) to the current pixel. If it crosses an
  872. odd number of edges, we color it. Otherwise, we don't.
  873.  
  874. This means that from infinite left to the first edge, we do not
  875. color. Then, from first edge to second edge, we color. Then, from
  876. second edge to third edge, we do not color. From third edge to
  877. fourth edge, we color. And so on. If you want, you can always
  878. draw a background bitmap left from the first edge, from second
  879. edge to third edge, etc.... Thus, you could put the picture of a
  880. beautiful blue sky behind at the same time you're drawing a
  881. polygon.
  882.  
  883. Now, let's start that over again. We are coming from infinitely
  884. upward and drawing. We do not need to consider every edge, only
  885. those that are on the current scanline (horizontal line). We will
  886. thus build an inactive edge list (those edges that have not yet
  887. been used in the drawing process) and an active edge list (those
  888. edges that are currently in use to draw the polygon, i.e. that
  889. intercept the current scanline). Thus, at first, all edges (named
  890. m,n,o,p in figure 5) would be in the inactive edge.
  891.  
  892. As the algorithm reaches the topmost line of the polygon, any
  893. relevant edge will be transfered from the inactive edge list to
  894. the active edge list. When an edge becomes unused (we're lower
  895. than the lowermost scanline), we remove it from the active list
  896. and throw it to the dogs. (# grin #)
  897.  
  898. Now you will realize that at each scanline we have to check every
  899. edge in the inactive edge list to see if we need to transfer them
  900. to the active edge list. To accelerate that, we can sort the
  901. inactive edge list in decreasing order. Thus, the ordered list
  902. would become {m,p,n,o}. Then, we need only to check the first
  903. vector at every scanline. Since they point downward, we need only
  904. to check the first point of it.
  905.  
  906. If you get any horizontal lines in the edge list, remove them
  907. from the list. Horizontal lines can be identified easily as b=d.
  908. Proving that it works is left as an exercise (# muhaha! #).
  909.  
  910. The polygon drawing algorithm (scan-line algorithm).
  911.  
  912. I=inactive edges list, initially all edges
  913. A=active edges list, initally empty
  914.  
  915. sort I with highest first endpoint first
  916.  
  917. for each scanline i starting from the top of the screen
  918.      does the first edge in I start at scan line i?
  919.           Yes, transfer all edges that should start now from I to
  920. A
  921.      find x value of all edges in A for current scanline (see the
  922. line algorithm below)
  923.      sort all edges in A with increasing x-intercept value (the
  924. thing we jest calculated)
  925.  
  926.      previous_x_intercept<-leftmost pixel on screen
  927.      for all edges n in A, counting by two, do
  928.           draw background between previous_x_intercept and edge n
  929. in A
  930.           fill polygon between edge n and edge n+1 in A
  931.           previous_x_intercept<-x intercept of edge n+1
  932.      end for
  933.  
  934.      draw background for the rest of the scan line (from
  935. previous_x_intercept to end of line)
  936.      for all edges for which we have reached the bottom
  937.           remove the edge from A
  938.      end for
  939. end for
  940.  
  941.  
  942.  
  943.  
  944. 2.1 A line algorithm for the polygon drawing algorithm
  945.  
  946.  
  947. This section discusses a somewhat shortened version of
  948. Bresenham's line drawing algorithm. It is optimized for drawing
  949. filled polygons.
  950.  
  951. You assuredly realize you need to calculate the x coordinate of
  952. the intersection of a line with an horizontal line. Well consider
  953. the line equation that follows:
  954.  
  955.      x=dy/dx x y + k
  956.  
  957. If the line is defined with the vector (a,b)-(c,d), then we
  958. define dy as d-b, dx as c-a. Say x0 is the x coordinate of the
  959. topmost point of the segment.
  960.  
  961. x0=a
  962.  
  963. Say xn is the x intercept value of the segment for the nth line
  964. from the top. Then, xn can be defined as follows:
  965.  
  966. xn=dy/dx x n + a
  967.  
  968. Say we know xn-1 and want to know xn based on that knowledge.
  969. Well it is obvious that:
  970.  
  971. xn-x(n-1)=dy/dx x n + a - [dy/dx x (n-1) + a]
  972. xn-x(n-1)=dy/dx
  973. xn=dy/dx+x(n-1)
  974.  
  975. Thus, knowing xn-1 and adding dy/dx yields xn. We can thus do an
  976. addition instead of a multiply and an add to calculate the new x
  977. intercept value of a line at each scanline. An addition been at
  978. least 5 times as fast as an add, usually more in the 15 times
  979. range, this will speed things up considerably. However, it might
  980. annoy us to use floating point or fixed point, even more so with
  981. roundoff error.
  982.  
  983. To this end, we have found a solution. dy/dx can be expressed as
  984. a+b/d, b<d. a is the integer part of dy/dx. Substituting, we find
  985. that
  986.  
  987. xn=a+b/d+xn-1
  988.  
  989. Thus, we always add a at each scanline. Furthermore, when we add
  990. a sufficient amount of b/d, we will add one more to xn. We will
  991. keep a running total of the fraction part. We will keep the
  992. denominator in mind, and set the initial numerator (when n=0) to
  993. 0. Then, at each line we increase the numerator by b. If it
  994. becomes greater than d, we add one to xn and decrease the
  995. numerator by d. In short, in pseudocode, we draw a segment from
  996. (x0,y0) to (x1,y1) this way:
  997.  
  998. dx=x1-x0
  999. dy=y1-y0
  1000. denominator=dy
  1001. numerator=dx modulo dy
  1002. add=dx divided (integer) dy
  1003. running=0
  1004. x=x0
  1005. for each line (0 to n) do
  1006.      x=x+add
  1007.      running=running+numerator
  1008.      if running>=denominator then
  1009.           running=numerator-denominator
  1010.           x=x+1
  1011.      end if
  1012.      do whatever you want with line here, such as drawing a pixel
  1013. at location (x,line)
  1014. next line
  1015.  
  1016.  
  1017.  
  1018.      __________________________
  1019.  
  1020.  
  1021.  
  1022.  
  1023. 3    3d polygon drawing
  1024.  
  1025.  
  1026.  
  1027. The first step is to realize that a projected 3d line with our
  1028. given projection becomes a 2d line, that is it stays straight and
  1029. does not bend. That can be easily proved. The simplest proof is
  1030. the geometric one: the intersection of two planes in 3d is a line
  1031. in the general case. The projector all lie in the plane formed by
  1032. the origin ant the line to be projected, and the projection plane
  1033. is, well, a plane. Thus, the projection is contained in a line.
  1034.  
  1035. This means that all the boundaries of a polygon can be projected
  1036. as simple lines. Only their endpoints need be projected. Then,
  1037. you can simply draw a line between the two endpoints of each
  1038. projected edge (or more precisely, draw a 2d polygon with the
  1039. given endpoints).
  1040.  
  1041. However, there comes the problem of overlapping polygons.
  1042. Obviously, the color of a pixel must be that of the closest
  1043. polygon in that point (well for our objective, that suffices).
  1044. But this means that we must determine the distance of each
  1045. polygon that has the possibility to be drawn in the current
  1046. pixel, and that for the whole image. That can be time consuming.
  1047.  
  1048. We will make a few assumptions. First, we will assume no polygons
  1049. are interpenetrating; this means that to come in front of another
  1050. polygon, a polygon must go around it. This means that you only
  1051. need to determine what polygon is topmost when you cross an edge.
  1052.  
  1053. Less obviously, but equally useful, if from one scanline to the
  1054. next, you encounter the same edges in the same order when
  1055. considering left to right, you do not need to recalculate the
  1056. priorities as it has not changed.
  1057.  
  1058. Inevitably, though, you will need to determine the depth of the
  1059. polygon (z value) in a given pixel. That can be done rather
  1060. easily if you carry with the polygon its plane equation in the
  1061. form Ax+By+Cz+D=0. Remember that the projection ray equation for
  1062. pixel (u,v) is x=zu and y=zv. Feed that x and y in the plane
  1063. equation and you find that
  1064.  
  1065. Azu+Bzv+Cz+D=0
  1066. z(Au+Bv+C)=-D
  1067. z=-D/(Au+Bv+C)
  1068.  
  1069. Thus you can easily calculate the z coordinate for any plane in a
  1070. given (u,v) pixel.
  1071.  
  1072. The cautious reader will have noticed that it is useful to do
  1073. some of these multiplications and additions incrementally from
  1074. scanline to scanline. Also note that if Au+Bv+C equals zero, you
  1075. cannot find z unless D is also zero. These (u,v) coordinates
  1076. correspond to the escape line of the plane of the polygon (you
  1077. can picture it as the line of «horizon» for the given plane). The
  1078. projection ray is parallel to the plane, thus you never cross it.
  1079.  
  1080. This brings me to another very important thing about planes.
  1081. Another way to define a plane is to give a normal vector. A
  1082. vector and a point of reference generate one and only one plane
  1083. perpendicular to the vector and passing through the point of
  1084. reference. Think of a plane. Think of a vector sticking up
  1085. perpendicular to it. That vector is called normal to the plane.
  1086. It can be easily proved that this vector is (A,B,C) with A,B,C
  1087. being the coefficients of the plane equation Ax+By+Cz=D. We might
  1088. want to normalize the vector (A,B,C) (make it length 1), or
  1089. divide it by SQRT(A2+B2+C2). But to keep the equation the same,
  1090. we must divide all member by that value. The equation thus
  1091. becomes A'x+B'y+C'z=D'. If we want to measure the distance
  1092. perpendicularly to the plane of a point (a,b,c), it can be shown
  1093. that this (signed) distance is A'a+B'b+C'c+D' (you can always
  1094. take the absolute value of it if you want an unsigned distance).
  1095. That might be useful for different reasons.
  1096.  
  1097. It is also notable that since the A, B and C coefficients are
  1098. also the coords of the normal vector, and that rotating a vector
  1099. can be done as we saw previously, it is easy to rotate the plane.
  1100. If the vector is normalized in the first place, it will stay
  1101. normalized after rotation. After rotation, you will ignore the D
  1102. component of the plane, but you do know the coordinates of a
  1103. point in the rotated plane. Any point in the polygon will fit the
  1104. bill. Thus, you can deduce D. In fact, when performing back-face
  1105. culling (see chapter 4), you will calculate D even as you are
  1106. performing some optimization.
  1107.  
  1108. Normally, you will build an inactive edge table (IET) in which
  1109. you will store all edges of all polygons to be rendered. Also
  1110. initialize an empty list called the active edge table (AET). In
  1111. the IET, make sure that the edges are sorted in increasing values
  1112. of y for the topmost endpoint of any edge. Also, as per the 2d
  1113. polygon drawing algorithm, discard any horizontal edge.
  1114.  
  1115. As you scan from the top of the screen to the bottom of the
  1116. screen, move the edges from the IET to the AET when you encounter
  1117. them. You do not need to check every edge in the IET as they are
  1118. sorted in increasing y order, thus you need only to check the
  1119. first edge on the IET.
  1120.  
  1121. As you are scanning from left to right, each time you cross an
  1122. edge, you toggle the associated polys on or off. That is, if poly
  1123. X is currently "off", it means you are not inside poly X. Then,
  1124. as you cross an edge that happens to be one of poly X's edge,
  1125. turn it "on". Of all polys that are "on", draw the topmost only.
  1126. If polys are not self-intersecting, you only need to recalculate
  1127. which one's on top when you cross edges. Your edge lists must
  1128. contain pointers to thir parent polygons. Here's a pseudo-code
  1129. for the algorithm:
  1130.  
  1131. Let polylist be the list of all poly
  1132. Put all edges of all polygon in IET
  1133. Empty AET
  1134. for Y varying from the topmost scanline to the last
  1135.      while the edge on top of the IET's topmost endpoint's Y
  1136. coodrinate equals Y do
  1137.           take the edge off the IET and insert it in the AET with
  1138. an insertion sort
  1139.           so that all edges in the AET are in increasing X
  1140. values.
  1141.      end_while
  1142.  
  1143.      discard any edge in AET for which the second endpoint's Y
  1144. value is Y
  1145.  
  1146.      previous_X=minus infinity
  1147.      current_poly=none
  1148.      for i=first edge through last edge in AET
  1149.           for X varying from previous_X to i's X coordinate
  1150.                draw a pixel at position (X,Y), attribute is taken
  1151. from current_poly
  1152.           end_for
  1153.           toggle polygons that are parent to i
  1154.           current_poly=find_topmost_poly()
  1155.      end_for
  1156.  
  1157.      update intercept values of all edges in AET and sort them by
  1158. increasing X intercept
  1159. end_for
  1160.  
  1161. A recent uproar about "s-buffers" made me explain this algorithm
  1162. a little more carefully. Paul Nettle wrote a FAQ explaining
  1163. something that he calls "s-buffering". It is essentially an
  1164. algorithm where, instead of calculating the spans on the fly as
  1165. above, you buffer them, that is, you precalculate all spans and
  1166. send them to a span buffer. Then you give the span buffer to a
  1167. very simple rasterizing function. It's a bit harder to take
  1168. advantage of coherence in that case (an example of coherence is
  1169. when we note that if the order of the edges do not change from
  1170. one scanline to the next, we don't need to recalculate which
  1171. spans are topmost, assuming polygons aren't interpenetrating).
  1172.  
  1173. Another very nice algorithm is the z-buffer algorithm, discussed
  1174. later. It allows interpenetrating polygon and permits mixing
  1175. rendering techniques in a rather efficient way (for real-time
  1176. rendering purposes).
  1177.  
  1178.  
  1179.  
  1180.  
  1181.      _________________________
  1182.  
  1183.  
  1184.  
  1185.  
  1186.  
  1187.  
  1188. 4    Data Structures
  1189.  
  1190.  
  1191.  
  1192.  
  1193. A very important aspect of 3d graphics is data structures. It can
  1194. speed up calculations by a factor of 6 and up, only if we put
  1195. some restriction on what can be represented as 3d entities.
  1196.  
  1197. Suppose that everything we have in the 3d world are closed
  1198. polyhedras (no self-crossing(s)!). Each polyhedra is constituted
  1199. of several polygons. We will need each polygon's normal and
  1200. vertexes. But how we list these is where we can improve speed.
  1201.  
  1202. It can be proved that each edge is shared by exactly two polygons
  1203. in a closed polyhedra. Thus, if we list the edges twice, once in
  1204. each polygon, we will do the transformations twice, and it will
  1205. take twice as long.
  1206.  
  1207. To that, add that each vertex is part of at least 3 points,
  1208. without maximum. If you list each point in every vertex it is
  1209. part of, you will do the transformations at least thrice, which
  1210. means it will take thrice the time. Three times two is six.
  1211. You're doing the same transformations at least six times over.
  1212.  
  1213. Thus, you need to list the vertexes in every polygon as pointers
  1214. to vertexes. This way, two polygons can share the same data for a
  1215. shared vertex, and will allow you to do the transformations only
  1216. once.
  1217.  
  1218. The same goes for edges, list the endpoints as pointers to
  1219. endpoints, thus you will spare yourself much runtime work.
  1220.  
  1221. Another interesting thing. All polygons have an interior and
  1222. exterior face in that model. Make it so that the normal to all
  1223. polygons point towards the exterior face. It is clear that, if
  1224. you are standing outside all objects, and all objects (polyhedra)
  1225. are closed, you cannot see an inside face. If, when transforming
  1226. polygons, you realize a polygon is facing away from the eye, you
  1227. do not to transform it, just skip it. To know if it's pointing
  1228. away from you do the following:
  1229.  
  1230. This is called back-face culling. Take the (a,b,c) vector from
  1231. the eye to any one endpoint of the polygon. (E.g., if the eye is
  1232. at (m,n,o) and a vertex of the poly is at (i,j,k), then the said
  1233. vector would be (i-m,j-n,k-o).) Take also (u,v,w) the normal
  1234. vector to the polygon. No matter what direction you are looking,
  1235. do the following operation:
  1236.  
  1237. au+bv+cw
  1238.  
  1239. If the result is positive, then it is facing away from you. If it
  1240. is negative, it is facing towards you. If it is null, it is
  1241. parallel to you. The operation is called scalar multiply of two
  1242. vectors. It is very notable that the above scalar multiplication
  1243. of two vectors is interestingly enough the D part of the plane
  1244. equation Ax+By+Cz=D.
  1245.  
  1246. You can also illuminate objects with a very distant light source
  1247. (e.g. the sun) very easily. Take the vector pointing away from
  1248. the light source (a,b,c) and the normal of the plane (u,v,w). Do
  1249. a scalar multiplication of both as above. The result can be
  1250. interpreted as the intensity of the lighting or shadow for the
  1251. polygon. The approximation is good enough to give interesting
  1252. results. You might also want to affect the shading of the polygon
  1253. according to its distance to origin, to yield a fog result.
  1254.  
  1255.  
  1256.  
  1257.  
  1258.      _______________________
  1259.  
  1260.  
  1261.  
  1262.  
  1263.  
  1264.  
  1265. 5    Texture mapping
  1266.  
  1267.  
  1268.  
  1269. What we are about to do would probably be more aptly named pixmap
  1270. mapping onto a plane. We want to take a pixmap (or bitmap,
  1271. although less appropriate a word, it is very popular) and «stick»
  1272. it to the surface of a polygon, cutting away any excess. The
  1273. result can be quite impressive. This, however, can be a painfully
  1274. slow operation. We will discuss here of the mathematics involved,
  1275. and after that, we will try to optimize them.
  1276.  
  1277. Let's give a polygon its own coordinate system, (i,j). Thus, all
  1278. points on the said polygon can be localized in (i,j), or (x,y,z),
  1279. or, once projected, in (u,v).
  1280.  
  1281.  
  1282.  
  1283. Figure 6
  1284.  
  1285. Figure 6 illustrates the problem. The blue plane is the
  1286. projection plane. The green triangle is a polygon we are
  1287. projecting on that plane and then scan-converting in (u,v)
  1288. coordinates. What we want to do is find, what (i,j) point
  1289. corresponds to any given (u,v) point when projected this way.
  1290. Then, the point (i,j) can be indexed in a pixmap to see what
  1291. color this pixel in (u,v) space should be.
  1292.  
  1293. Let the i and j vectors be expressed in (x,y,z) coordinates,
  1294. which they should normally be. We want a general solution.
  1295.  
  1296. i=(a,b,c)
  1297. j=(d,e,f)
  1298.  
  1299. Let also the origin of (i,j) space be pointed to in (x,y,z) by
  1300. the vector k=(m,n,o). Thus, a point at location (p,q) in (i,j)
  1301. space is the same as a point at location pi+qj+k in (x,y,z)
  1302. space. Furthermore, we have the projection ray defined as follow:
  1303. it shoots through point (r,s) in (u,v) space. This corresponds to
  1304. point (r,s,1) in (x,y,z) space. The equation of that ray is thus
  1305. t(r,s,1) where t can take any real value. Simplifying, we find
  1306. that the ray is (tr,ts,t), or
  1307.  
  1308.      x=tr
  1309. R:   y=ts
  1310.      z=t
  1311.  
  1312. or
  1313.      x=zr
  1314. R:   y=zs
  1315.      z=z
  1316.  
  1317. The point in the plane is
  1318.      x=pa+qd+m (1)
  1319. P:   y=pb+qe+n (2)
  1320.      z=pc+qf+o (3)
  1321.  
  1322. We want to find (p,q) as a function of (r,s). Here is what i did
  1323. in stupid MathCAD. [MathCAD flame deleted] I'm sorry I didn't use
  1324. the same variable names in MathCAD...
  1325.  
  1326. #######Begin picture of mathcad screen#########
  1327.  
  1328. ########End picture of mathcad screen#########
  1329.  
  1330. These equations for p and q can be optimized slightly for
  1331. computers. Since it is usually faster on a computer to do a(b+c)
  1332. than ab+ac, and that both are equivalent, we will try to
  1333. factorize a little.
  1334.  
  1335. p=(Mv+Nu+O)/(Pv+Qu+R)
  1336. where M,N,O,P,Q,R are as follows:
  1337. M=CXq-AZq
  1338. N=BZq-CYq
  1339. O=AYq-BXq
  1340. P=ZqXp-ZpXq
  1341. Q=YqZp-YpZq
  1342. R=YpXq-YqXp
  1343.  
  1344. They are all constants, so they need only to be calculated once
  1345. to draw the whole image.
  1346.  
  1347. Similarly for q, we have:
  1348. q=(Sv+Tu+U)/(Vv+Wu+X)
  1349. and
  1350. S=AZp-CXp
  1351. T=CYp-BZp
  1352. U=BXp-AYp
  1353. V=ZqXp-ZpXq
  1354. W=YQZP-YPZQ
  1355. X=YqXp-YpXq
  1356.  
  1357. All these constants should be calculated only once. Now we can
  1358. simplify a little with, knowing that P=V, Q=W, R=X.
  1359.  
  1360. p=(Mv+Nu+O)/(Pv+Qu+R)
  1361. q=(Sv+Tu+U)/(Pv+Qu+R)
  1362.  
  1363. This is noteworthy: (Q,P,R)(u,v,1)=uQ+vP+R, the denominator. But,
  1364. (Q,P,R)=(Xq,Yq,Zq)x(Xp,Yp,Zp). This cross-product of two vectors
  1365. will yield a vector that is perpendicular to the first two (i.e.
  1366. a normal to the plane). If these two vectors are of length 1,
  1367. then the cross product will also be of length 1. But we might
  1368. already have the normal from previous transformations (e.g. if
  1369. we're doing back-face culling). Then, we won't need to calculate
  1370. (Q,P,R) as it is the normal to the plane. Or we can deduce the
  1371. normal from these equations if we don't have it.
  1372.  
  1373. The denominator needs to be calculated only once per pixel, not
  1374. for both p and q. The numerator for p and q is different, though.
  1375. Furthermore, It is obvious that Mv+O are constant throughout a
  1376. constant-v line (horizontal line). Thus, it needs to be
  1377. calculated only once per line. We can spare ourselves the
  1378. multiply per pixel doing it incrementally. (We're still stuck
  1379. with a divide, though.) Example, let V=Mv+O, for a given v. Then,
  1380. when calculating Mv+Nu+O, we can use V+Nu. If we know V+Nua, then
  1381. we can calculate V+N(ua+1) as follow: V+Nua+N
  1382.  
  1383. Let W=V+Nua, then we get
  1384.  
  1385. W+N
  1386.  
  1387. That's a simple add instead of a multiply. But we can not get rid
  1388. of the divide, so we're stuck with a divide per pixel (and at
  1389. least two adds) per coordinate.
  1390.  
  1391. There is a way to solve this, however, and here's how to do it.
  1392.  
  1393. First, it is interesting to note the following. The plane
  1394. equation is Ax+By+Cz=D, and x=zu, y=zv, where (u,v) is the screen
  1395. coordinates of a point, and (x,y,z) is the corresponding world
  1396. coordinate of the unprojected point. Thus,
  1397.  
  1398. Auz+Bvz+Cz=D
  1399.  
  1400. What happens when z is a constant? Let us express v as a function
  1401. of u.
  1402.  
  1403. Bvz=-Auz-Cz+D
  1404. v=(-A/B)u+(-C/B)+D/(Bz)
  1405.  
  1406. Let M=-A/B, and N=d/(Bz)-C/B then we have
  1407.  
  1408. u=Mv+N
  1409.  
  1410. a linear equation. Furthermore, M is independant of z. Thus, a
  1411. slice of constant z in the plane is projected as a line on (u,v).
  1412. All of these slices project in lines that are all parallel to
  1413. each other (e.g. the slope, M, is independant of the slice
  1414. taken).
  1415.  
  1416. Now, a slice of constant z in the original Ax+By+Cz=D plane is a
  1417. line. Since z is constant throughout that line, all points of
  1418. that line are divided by the same number, independantly of (x,y).
  1419. Thus, the projected line is merely the original line scaled up or
  1420. down.
  1421.  
  1422. Up to now, we have been rendering polygons with a scanline
  1423. algorithm where the scanlines are horizontal. First, we have to
  1424. conceive that the scanlines need not be horizontal, merely
  1425. parallel. If we use a slope of M for the scanlines (M=-A/B, as
  1426. above), then we will be performing scans of constant z value. Two
  1427. advantages are very obvious. First, it is easy to calculate the z
  1428. value for a large number of pixel (it is constant throughout a
  1429. scanline) and second, texture mapping the scanline becomes a mere
  1430. scaling operation, which can be done incrementally. Furthermore,
  1431. M need be calculated only once for the whole polygon, while N
  1432. requires at most a divide and an add per scanline.
  1433.  
  1434. Note that if abs(A)>abs(B), you might want to express v as a
  1435. function of u instead of the above, so that the slope ends up as
  1436. being -B/A instead. (E.g., make certain that the slope is less
  1437. than or equal to 1 in absolute value). There is a degenerate
  1438. case: A=B=0. This is the case where the whole plane is of
  1439. constant z value. In that case, you can use horizontal scan lines
  1440. and just scale the texture by a factor of 1/z.
  1441.  
  1442. Thus, we're going to be rasterizing the polygon with
  1443. nonhorizontal scanlines of slope M. The equation of a scanline
  1444. can be expressed as:
  1445.  
  1446. v=Mu+Vi
  1447.  
  1448. Where Vi is the v coordinate when u is 0. All scanlines being
  1449. parallel, M never changes, and Vi identifies uniquely each and
  1450. every scanline. Say i=0 corresponds to the topmost scanline. Say
  1451. we want to find the intersection of the scanline with a polygon
  1452. edge. What we're really interested in is the u value of the
  1453. intersection of the two lines. The edge can be represented as
  1454. being:
  1455.  
  1456. v=Nu+K
  1457.  
  1458. Where K is the v value for the edge when u=0. Now, intersecting
  1459. the two yields
  1460.  
  1461. Nu+K=Mu+Vi
  1462. (N-M)u=(Vi-K)
  1463. u=(Vi-K)/(N-M)
  1464.  
  1465. That is, assuming N-M is not zero, which shouldn't be because
  1466. we're not including edges parallel to the scanlines in our scan
  1467. line algorithm, thus N-M is never zero.
  1468.  
  1469. Now, V(i+1) is (Vi)+1. When we go from scanline i to scanline
  1470. i+1, V is incremented by one. We thus get u', the intersection of
  1471. the two lines for the next scanline, as being:
  1472.  
  1473. u'=(Vi+1-K)/(N-M)
  1474. u'=(Vi-K)/(N-M)+1/(N-M)
  1475. u'=u+1/(N-M)
  1476.  
  1477. It is thus possible to calculate u from scanline to scanline
  1478. incrementally, bresenham stlye. To clarify things up a bit: if P1
  1479. and P2 delimit the edge, then let (p,q) be P2 minus P1 (that is,
  1480. p is the delta u and q is the delta v). Therefore, N is q/p (by
  1481. definition of slope). Now, we know M to be of the form -B/A (from
  1482. the plane equation, above). We need to calculate N-M for our
  1483. incremental calculations, which is relatively easy. We have
  1484.  
  1485. N-M=q/p-(-A/B)
  1486. N-M=q/p+A/B
  1487. N-M=(qB+Ap)/(Bp)
  1488.  
  1489. Thus, u' is really:
  1490. u'=u+Bp/(qB+Ap)
  1491.  
  1492. At one point in the algorithm, we will need to calculate which
  1493. scanline contains a particular point. As an example, when sorting
  1494. from topmost edge downwards, we need to use a definition of top
  1495. that is perpendicular to the scanline used. Then, we have the
  1496. slope of the scanline, which is, say, m. If the equation of the
  1497. line that interests us is:
  1498.  
  1499. y=mx+b
  1500.  
  1501. Then, we can sort the edges according to increasing b values in
  1502. the above equation. If we want to calculate b for point (i,j),
  1503. simply isolate b in the above equation.
  1504.  
  1505. This last bit is called constant-z texture mapping. It is not
  1506. particularly useful in the general case, for a number of reasons,
  1507. mainly added aliasing, and span generation problems. However, it
  1508. has a very well known and very used application. If we have
  1509. vertical walls and horizontal floors exclusively, and that the
  1510. eye is never banked, then the following is always true: the
  1511. constant-z lines are always either horizontal or vertical. For
  1512. example, the walls' constant-z lines are vertical while the
  1513. floors' and ceilings' are horizontal. So basically, you can
  1514. render floors and walls using the techniques described in this
  1515. chapter in a very straightforward manner, if you use both
  1516. horizontal and vertical scanlines. Just use the last bit about
  1517. constant-z texture mapping and optimize it for horizontal and
  1518. vertical scanlines. It is to note that looking up or down will
  1519. not affect the fact that walls will have vertical scan-lines and
  1520. floors and ceilings will have horizontal scan-lines.
  1521.  
  1522.  
  1523.  
  1524.  
  1525.      ________________________
  1526.  
  1527.  
  1528.  
  1529.  
  1530.  
  1531. 6    Of logarithms
  1532.  
  1533.  
  1534.  
  1535. If we have many multiplies or divides to make, little time and
  1536. lots of memory, and that absolute precision is not utterly
  1537. necessary, we may use the following tricks.
  1538.  
  1539. Assuming x and y are positive or null real numbers, and b a base
  1540. for a logarithm, we can do the following: (note: x**y means x
  1541. raised to the yth power, logby denotes base b logarithm of y).
  1542.  
  1543. logb(xy)=logbx+logby
  1544. b**logb(xy)=b**(logbx+logby)
  1545. xy=b**(logbx+logby)
  1546.  
  1547. Now, you might say that taking the log or the power of a value is
  1548. very time consuming, and it usually is. However, we will reuse
  1549. the trick that we used to do not so long ago, when pocket
  1550. computers were not available, we will make a logarithm table. It
  1551. is easy to store the values of all possible logaritms within a
  1552. finite, discrete range. As an example, if we are using 16 bit
  1553. fixed points or integer, we can store the result of a log of that
  1554. number as a table with 65536 entries. If each entry is 32 bits
  1555. wide, which should be sufficient, the resulting table will use
  1556. 256kb of memory. As for the power table, you can build another
  1557. lookup table using only the 16 most significant bits of the
  1558. logarithm, yielding a 16 or 32 bit wide integer or fixed point.
  1559. This will yield a 128 to 256kb table, for a grand total of under
  1560. 512kb. With today's machines, that should not be a problem. The
  1561. only problem is that you must then limit your scope to 16 bits
  1562. numbers. You realize of course that you can not multiply if one
  1563. or both of the numbers are negative. Well, you can do it this
  1564. way. Say x<0, then let u=-x, thus u>0. Then, xy=-uy, and uy you
  1565. can calculate with this method. Similarly,
  1566.  
  1567. x/y=b**(logbx-logby)
  1568.  
  1569. Thus divides are also a simple matter of looking it up in the
  1570. table and subtracting.
  1571.  
  1572. Powers are made much more simple also. I won't say anything long
  1573. winded about powers, but here is, in short:
  1574.  
  1575. logb(x**y)=y logbx
  1576. x**y=b**(y logbx)   (*)
  1577. And from there, even better yet:
  1578. x**y=b** ( b** [ logby+logb(logbx) ] )
  1579.  
  1580. For y=constant, some shortcuts are usually possible. Since ax can
  1581. sometimes be very quickly calculated with left shifts, you can
  1582. use (*) instead of the last equation. Very common example:
  1583. 320y=(256+64)y
  1584. 320y=256y+64y
  1585. 320y=(y<<8)+(y<<6)
  1586.  
  1587. Where x<<y denotes "x, left shifted y times". A left shift of y
  1588. is equivalent of multiplying by 2**y. The ancient egyptians only
  1589. knew how to multiply this way (at least, as far as I know, for a
  1590. while, like for the pyramids and stuff - I think they even
  1591. calculated an approximation of pi during that time, and they
  1592. couldn't even multiply efficiently! :-)
  1593.  
  1594.  
  1595.  
  1596.      ______________________
  1597.  
  1598.  
  1599.  
  1600.  
  1601. 7    More on data structures and clipping
  1602.  
  1603.  
  1604.  
  1605. Usually, the space viewed by the eye can be defined by a very
  1606. high square-based pyramid. Normally, it has no bottom, but in
  1607. this case it might have one for practical purposes. Let us say
  1608. that the eye has a width of vision (horizontally and vertically)
  1609. of 90 degrees. (Since it is square, that angle is slightly larger
  1610. when taken from one corner to the opposite.) Thus, we will see
  1611. 45° to both sides, up and down. This means a slope of one. The
  1612. planes bounding that volume can be defined as x<z, x>-z, y<z and
  1613. y>-z. To a programmer with an eye for this, it means that it easy
  1614. to «clip», i.e. remove obviously hidden polygons, with these
  1615. bounds. All polygons that do not satisfy either of these
  1616. inequalities can be discarded. Note that these polygons must fall
  1617. totally outside the bounding volume. This is called trivial
  1618. rejection.
  1619.  
  1620. For polygons that are partially inside, partially outside that
  1621. region, you may or may not want to clip (I am not sure yet of
  1622. which is more efficient - I suspect it depends on the situation,
  1623. i.e. number of polygons, number of edges in each polygon...)
  1624.  
  1625. There is also a slight problem with polygons that are partially
  1626. behind, partially in front of the eye. Let us give a little
  1627. proof. All of our polygon-drawing schemes are based on the
  1628. assumption that a projected line is also a line. Here is the
  1629. proof of that.
  1630.  
  1631. Let us consider the arbitrary line L, not parallel to the plane
  1632. z=0, defined parametrically as:
  1633.      x= at+b
  1634. L:   y= ct+d
  1635.      z= t
  1636.  
  1637. Now, let us define the projection P(x,y,z) as follows:
  1638.  
  1639. P(x,y,z): u=x/z
  1640.      v=y/z
  1641.  
  1642. Then, P(L) becomes
  1643.  
  1644. P(L):     u=(at+b)/t
  1645.      v=(ct+d)/t
  1646.  
  1647. We are using a segment here, not an infinite line. If no point
  1648. with t=0 needs to be projected, we can simplify:
  1649.  
  1650. P(L):     u=a+b/t
  1651.      v=c+d/t
  1652.  
  1653. Let us measure the slope of that projection:
  1654.  
  1655. du/dv=(du/dt) / (dv/dt) = (-b/t2)/(-d/t2)
  1656.  
  1657. Since t is not 0
  1658.  
  1659. du/dv=b/d, a constant, thus the parametric curve defined by P(L)
  1660. is a straight line if t is not zero.
  1661.  
  1662. However, if t is null, the result is indeterminate. We can try to
  1663. find the limits, which are pretty obvious. When t tends towards
  1664. zero+, then u will tend towards positive (or negative if b is
  1665. negative) infinity. If t tends towards zero-, then u will tend
  1666. towards negative (or positive, if b is negative) infinity. Thus,
  1667. the line escapes to one infinity and comes back from the other if
  1668. t can get a zero value. This means that if the segment crosses
  1669. the plane z=0, the projection is not a line per se.
  1670.  
  1671. However, if the line is parallel to z=0, L can be defined as
  1672. follows:
  1673.      x=at+b
  1674. L:   y=ct+d
  1675.      z=k
  1676.  
  1677. Where k is the plane that contains L. P(L) becomes
  1678.  
  1679. P(L):     u=(at+b)/k
  1680.      v=(ct+d)/k
  1681.  
  1682. The equations do not admit a null k. For a non-null k, we can
  1683. find the slope:
  1684.  
  1685. du/dv = (du/dt)/(dv/dt) = (a/k)/(c/k)
  1686. du/dv=a/c
  1687.  
  1688. Still a constant, thus P(L) is a line. However, if k is null, we
  1689. do not get a projection. As k gets close to zero, we find that
  1690. the projection lies entirely at infinity, except for the special
  1691. case b=d=0, which is an indetermination.
  1692.  
  1693. This will have a practical impact on our 3d models. If a cube is
  1694. centered on the eye, it will be distorted unless we somehow take
  1695. that into consideration. The best way to do it, with little
  1696. compromise and efficient coding, is to remove all parts of
  1697. polygons that lie in the volume defined by z<=0. Of course, if an
  1698. edge spans both sides of the plane z=0, you will have to cut it
  1699. at z=0, excluding everything defined by z<=0. Thus, you will have
  1700. to make sure the point of the edge with the smaller z coordinate
  1701. is strictly greater than 0. Normally, you would have to take the
  1702. next higher real number (whatever that might be!) but due to
  1703. hardware contraints, we will use any arbitrarily small number,
  1704. such as z=.001 or even z=.1 if you are building a real-time
  1705. something in which infinite precision is not required.
  1706.  
  1707. Lastly, you might not want to render objects that are very far
  1708. away, as they may have very little impact on the global scene, or
  1709. because you are using a lookup table for logarithms and you want
  1710. to make sure that all values are within range. Thus, you will set
  1711. an arbitrary maximum z for all objects. It might even be a fuzzy
  1712. maximum z, such as «I don't wand any object to be centered past
  1713. z0». Thus, some objects may have some part farther than z0, but
  1714. that may be more efficient than the actual clipping of the
  1715. objects.
  1716.  
  1717. Moreover, you might want to make several polygon representations
  1718. of the objects, with different levels of details as measured by
  1719. the number of polygons, edges and vertices, and show the less
  1720. detailed ones at greater distances. Obviously, you would not
  1721. normally need to rasterize 10 polygons for an object that will
  1722. occupy a single pixel, unless you are doing precision rendering.
  1723. (But then, in that case, you'll have to change the whole
  1724. algorithm...)
  1725.  
  1726. You may want to perform your trivial rejection clipping quickly
  1727. by first rotating the handle coordinates of each polyehdra (the
  1728. handle coordinate is the coordinate of its local coordinates'
  1729. origin), and performing trivial rejection on them. As an example,
  1730. if you determine that all of a polyhedra is contained within a
  1731. sphere of radius r, and that this polyhedra is «centered» at z=-
  1732. 3r, then it is obvious that it is entirely outside the view
  1733. volume.
  1734.  
  1735. Normally, the view volume is exactly one sixth of the total
  1736. viewable volume (if you do not take into account the bottom of
  1737. the pyramid). Thus, you will statistically eliminate five sixths
  1738. of all displayable polyhedras with such a clipping algorithm.
  1739.  
  1740. Let's make a brief review of our speed improvements over the
  1741. brute-force-transform-then-see-no-think algorithm. We multiplied
  1742. the speed by six with our structure of pointers to edges that
  1743. point to vertexes. Then, we multiply again by two (on average)
  1744. with back-face culling, and now by six again. The overall
  1745. improvement is a factor of 72! This shows you how a little
  1746. thinking pays. Note however that that last number is a bit skewed
  1747. because some operations do not take as long as some others, and
  1748. also by the fact that most polyhedra do not contain a nearly
  1749. infinite amount of vertices.
  1750.  
  1751. We have not taken into account several optimizations we
  1752. mentioned, such as the lessening of detail, a bottom to the view
  1753. volume and we did not compare with other kinds of algorithm. But
  1754. another appeal of this algorithm is that it draws once, and only
  1755. once to each pixel, thus allowing us to superimpose several
  1756. polygons with minimal performance loss, or make complicated
  1757. calculations for each pixel knowing that we will not do it twice
  1758. (or more!) for some pixels. Thus, this algorithm lends itself
  1759. very well to texture mapping (or more properly, pixmap mapping),
  1760. shading or others.
  1761.  
  1762. However, we have to recognize that this approach is at best an
  1763. approximation. Normally, we would have to fire projection rays
  1764. over the entire area covered by a pixel (an infinite number of
  1765. rays), which will usually amount to an integral, and we did not
  1766. perform any kind of anti-aliasing (more on that later), nor did
  1767. we filter our signal, nor did we take curved or interpenetrating
  1768. surfaces into account, etc... But it lends itself very well to
  1769. real-time applications.
  1770.  
  1771.  
  1772.  
  1773.  
  1774.  
  1775. 8    A few more goodies...
  1776.  
  1777.  
  1778.  
  1779. Anti-aliasing is something we invented to let us decrease the
  1780. contrast between very different pixels. More properly, it should
  1781. be considered an approximation of high frequencies in a
  1782. continuous image through intensity variations in the discrete
  1783. view space. Thus, to reduce the stairway effect in a diagonal
  1784. line, we might want to put some gray pixels.
  1785.  
  1786.  
  1787.  
  1788. Figure 7
  1789.  
  1790. Figure 7 demonstrates this. This diagonal line a was «anti-
  1791. aliased» manually (I hear the scientists scream) in b, but it
  1792. should let you understand the intention behind anti-aliasing.
  1793.  
  1794. Scientists use filters, and so will we, but I will not discuss
  1795. fully the theory of filters. Let us see it that way. We first
  1796. generate a normal image, without any kind of anti-aliasing. Then
  1797. we want to anti-aliase it. What we do is take each pixel, modify
  1798. it so it resembles its non-anti-aliased neighbors a little more
  1799. and show that on screen.
  1800.  
  1801. Let us further imagine that how we want to do that specifically
  1802. is to take one eight of each of the four surrounding pixel
  1803. (north, south, east and west) and one half of the «real» pixel,
  1804. and use the resulting color for the output pixel.
  1805.  
  1806. If our image is generated using a 256-color palette, there comes
  1807. a handy little trick. Let us imagine that we want to take one
  1808. half of color a and one half of color b and add them together.
  1809. Well, there are 256 possible values for a. And for each single
  1810. value of a, there are 256 values of b, for a total number of
  1811. 256x256 pairs, or 64k pairs. We can precalculate the result of
  1812. mixing one half a with one half b and store it in a matrix 64kb
  1813. large.
  1814.  
  1815. Now let us call the four surrounding pixels n,s,e,w and the
  1816. central point o. If we use the above matrix to mix n and s, we
  1817. get a color A that is one half n plus one half s, that is
  1818. A=n/2+s/2. Then, mix e and w as B=e/2+w/2. Then, mix A and B.
  1819. C=A/2+B/2=n/4+s/4+e/4+w/4. Lastly, mix C and o, and we get
  1820. D=o/2+n/8+s/8+e/8+w/8, the exact desired result. Now, this was
  1821. all done through a lookup table, thus was very fast. It can
  1822. actually be done in real time.
  1823.  
  1824. Thus, if you want to do this, you will first create the original,
  1825. un-anti-aliased image in a memory buffer, and then copy that
  1826. buffer as you anti-aliase it on screen. Obviously, it will be
  1827. slower a little, but you might want to consider the increased
  1828. picture quality. However, you will probably want to do this in
  1829. assembler, because it can be highly optimized that way (most
  1830. compilers get completely confused...)
  1831.  
  1832. Another very nice technique is this. For pixel (u,v) on screen,
  1833. mix equal parts of pixels (u,v), (u+1,v), (u,v+1), (u+1,v+1) in
  1834. your rendering buffer. It's as if the screen was half a pixel of
  1835. center in relation to the rendering buffer. You'll need to render
  1836. an extra row and column though.
  1837.  
  1838.      __________________
  1839.  
  1840.  
  1841. It is also interesting to note that the silhouette of a projected
  1842. sphere is a circle. The centre of the circle is the projection of
  1843. the centre of the sphere. However, the radius is NOT the
  1844. projection of the radius of the sphere. However, when the
  1845. distance is large enough when compared to the radius of the
  1846. sphere, this approximation is good enough.
  1847.  
  1848.      __________________
  1849.  
  1850.  
  1851. To determine if a point is contained inside a polyhedra, we
  1852. extend our test from the polygons like this: if a line drawn from
  1853. the point to infinity crosses an odd number of faces, the it is
  1854. contained, otherwise it is not.
  1855.  
  1856.      __________________
  1857.  
  1858.  
  1859. If you want to approximate a continuous spectrum with 256 colors
  1860. (it will not look perfect though), use base 6. That is, using the
  1861. RGB model (we could use HSV or any other), you will say that red,
  1862. green and blue can all take 6 different values, for a total
  1863. number of combinations of 216. If you add an extra value to green
  1864. or red, you get 252 colors, close to the maximum. However, 216 is
  1865. not bad and it leaves you 40 colors for system purposes, e.g.
  1866. palette animations, specific color requirements for, let's say,
  1867. fonts and stuff. Under Microsoft Windows, 20 colors in the
  1868. palette are already reserved by the system, but you still have
  1869. enough room for your 216 colors. By the way, I think it's stupid
  1870. that MS Windows does not allow you to do that as a standard
  1871. palette allocation. You have to define your own video driver to
  1872. enable programs that limit themselves to dithering to use those
  1873. colors. And that's living hell.
  1874.  
  1875. In addition, if the colors are arranged in a logical way, it will
  1876. be easier to find the closest match of the mix of two colors as
  1877. per anti-aliasing, above. (A bit like cheaper 24 bit color.)
  1878.  
  1879. Some people give 8 values to red and green, and 4 values to blue,
  1880. for a total of 256 colors, arguing that shades of blue are
  1881. perceived more difficultly. Personnally, I think the result is
  1882. not as good as the base six model I previously described.
  1883.  
  1884. Some other nice techniques: divide your palette in spans.
  1885. Example, allocate palette registers 16-48 to shades of blue. Use
  1886. only registers 20-44 in your drawings, but arrange so that it
  1887. goes from white-blue at color #16 through black-blue at color
  1888. #48. Then, shading is a simple matter of adding something to make
  1889. it darker, or subtracting to make it lighter. Make a few spans
  1890. for your essential colors, and maybe keep 30 for anti-aliasing
  1891. colors. These should not be used in your pixmaps. They could be
  1892. used only when performing final anti-aliasing of the rendered
  1893. picture. Pick these colors so to minimize visual artifacts. (E.g.
  1894. you have a blue span at 16-48, and a red span at 64-100. But you
  1895. don't have any purple, so if blue happens to sit next to red,
  1896. anti-aliasing will be hard to achieve. Then, just a bunch of your
  1897. special 30 colors for shades of purple.) Another particular trick
  1898. to increase the apparent number of colors (or rather, decrease
  1899. the artifacts created by too few colors), you could restrict your
  1900. intensity to lower values. E.g. make your darkest black/gray at
  1901. 25% luminosity, and your lightest white at 75%. This will make
  1902. your graphics look washed out, but your 256 colors will be spread
  1903. over a narrower spectrum.
  1904.  
  1905.  
  1906.  
  1907.  
  1908.      _______________________
  1909.  
  1910.  
  1911.  
  1912.  
  1913.  
  1914. 9    Sorting
  1915.  
  1916.  
  1917.  
  1918. Sorting becomes an important issue because we do it a lot in our
  1919. 3d rendering process. We will discuss here briefly two general
  1920. sorting algoritmthms.
  1921.  
  1922. The first sort we will discuss is the bubble sort. It is called
  1923. this way because the elements seem to «float» into position as
  1924. the sort progresses, with «lighter» elements floating up and
  1925. «heavier» ones sinking down. Given the set S, composed of
  1926. elements s0, s1, ..., sn-1. This set has n elements. Let us also
  1927. suppose that we can compare two elements to know which of the two
  1928. is «lighter», i.e. which one goes first. Let us also suppose that
  1929. we can exchange two consecutive elements in the set S. To sort
  1930. the set S, we will proceed as follow. Starting at element 0 up to
  1931. element n-2, we will compare each element to its successor and
  1932. switch if necessary. If no switches are made, the set is sorted.
  1933. Otherwise, repeat the process. In pseudocode, we have:
  1934.  
  1935. repeat the following:
  1936.      for i varying from 0 to n-2, repeat the following
  1937.           if si should go after si+1, then exchange si and si+1
  1938.      end of the loop
  1939. until no elements are switched in the for loop (i.e the set is
  1940. sorted)
  1941.  
  1942. This algorithm, while very inefficient, is very easy to
  1943. implement, and may be sufficient for relatively small sets (n of
  1944. about 15 or less). But we need a better algoritm for greater
  1945. values of n.
  1946.  
  1947. However, we have to keep in mind that we are sorting values that
  1948. are bounded; that is, these values have a finite range, and it is
  1949. indeed small, probably 16 bits or 32 bits. This will be of use
  1950. later. But let us introduce first the sort. It is called the
  1951. radix sort.
  1952.  
  1953. Assume we are give a set S of numbers, all ranging from 0 to 9999
  1954. inclusively (finite range). If we can sort all si according to
  1955. the number in the ones' place, then in the tens' place, then in
  1956. the hundreds' place, then thousands' space, we will get a sorted
  1957. list. To do this, we will make a stack for each of the 10
  1958. possible values a digit can take, then take each number in the S
  1959. set and insert them in the appropriate stack, first for the
  1960. rightmost digit. Then, we will regroup the 10 stacks and start
  1961. over again for the second rightmost digit, etc... This can be
  1962. done in pseudocode as follows:
  1963.  
  1964. divisor=1
  1965. while divisor<10000, do the following:
  1966.      for i varying from 0 to n-1, do the following:
  1967.           digit=(si divided (integer version) by divisor) modulo
  1968. 10
  1969.           put si on stack # digit
  1970.      end of for loop
  1971.      empty set S
  1972.      for i varying from 0 to 9, do the following:
  1973.           take stack # i and add it to set S
  1974.      end of for loop
  1975.      divisor=divisor x 10
  1976. end of while loop
  1977.  
  1978.  
  1979. Now, let us assume that we want to work with hexadecimal numbers.
  1980. Thus, we will need 16 stacks. In real life, what we want to sort
  1981. is edges or something like that. If we limit ourselves to 16384
  1982. edges (16k, 12 bits), then our stacks need to be able to fit
  1983. 16384 objects each. We could implement them as arrays and use
  1984. lots of memory. A better way is to use lists.
  1985.  
  1986.  
  1987.  
  1988. One could also use the quicksort algorithm, which has a
  1989. complexity of O(nxln(n)), which seems worse than n for this
  1990. particular problem. It is also a bit more tedious to code (in my
  1991. humble opinion). In our problem, we may not even need to sort for
  1992. 16 bits. The only place that we need a powerful sorting algorithm
  1993. is when sorting all edges from top to bottom and from left to
  1994. right. If we are certain that all coordinates are bound between,
  1995. say, 0 and 200 (e.g. for y coordinate), then we can sort on 8
  1996. bits. For the x coordinate, though, depending if we already
  1997. clipped or not, we may have to sort on the full 16 bits. When
  1998. using the z-buffer algorithm (decribed scantily elsewhere), you
  1999. want to do a z-sort, then drawing the polygons that are on top
  2000. FIRST so that you can avoid these nasty pixel writes.
  2001.  
  2002.  
  2003.  
  2004.  
  2005.  
  2006.      ______________________
  2007.  
  2008.  
  2009.  
  2010.  
  2011.  
  2012.  
  2013.  
  2014. 10   Depth-field rendering and the Z-buffer algorithm
  2015.  
  2016.  
  2017.  
  2018.  
  2019. You can approximate any 3d function z=f(x,y) with a depth field,
  2020. that is, a 2d matrix containing the z- values of f(x,y). Thus,
  2021. z=M[x][y]. Of course, if you want any precision at all, it has
  2022. somewhat large memory requirements. However, rendering this depth-
  2023. field is fairly efficient and can give very interesting results
  2024. for smooth, round surfaces (like a hill range or something).
  2025.  
  2026. To render a depth field, find the vector V from our eye to the
  2027. pixel in the projection plane in 3space. Shoot a ray, as before,
  2028. using that vector. Start at the eye and move in direction V a
  2029. little. If you find that you are still above the surface (depth
  2030. field), then keep going (add another V and so on) until you hit
  2031. or fall under the surface. When that happens, color that pixel
  2032. with the color of the surface. The color of the surface can be
  2033. stored in another matrix, C[x][y]. In pseudocode, we get:
  2034.  
  2035.  
  2036. M[x][y]=depth field
  2037. C[x][y]=color field
  2038. for each pixel to be rendered, do:
  2039.      V=projection ray vector
  2040.      P=Eye coordinates
  2041.      while (M[Px][Py]<Pz)
  2042.           P=P+V
  2043.      end while
  2044.      color the current pixel with C[Px][Py]
  2045.      if we are using z-buffer algorithm, then set Z[Px][Py] to Pz
  2046. (see below)
  2047. end for
  2048.  
  2049.  
  2050. Of course, you will want to break the loop when P is too far from
  2051. origin (i.e. when you see that nothing is close through that
  2052. projection ray, stop checking).
  2053.  
  2054. For vector graphics, we were forced to place our eye centered on
  2055. origin and looking towards positive z values. That forced us to
  2056. rotate everything else so the eye would be positioned correctly.
  2057. Now, obviously, we cannot rotate the whole depth field as it
  2058. would be very time- and memory-consuming. Instead, the eye now
  2059. can be centered anywhere (the initial value for P need not be
  2060. (0,0,0)) and the vector V can be rotated into place. Even better,
  2061. if we can find the vector W from the eye to the center of the
  2062. projection plane, and the two base vectors u and v in 3d for the
  2063. projection plane, all of the other vectors Vu,v wil be a linear
  2064. combination of the form: V=W+au+bv where (a,b) is the coordinate
  2065. of the pixel in the (u,v) coodinates system. If we are doing this
  2066. rendering with scan lines (as we most probably are doing), then
  2067. the above equation can be done incrementally, e.g. say Va,y=A,
  2068. then Va+1,y=A+u. No multiplications whatsoever.
  2069.  
  2070. The z coordinate of each pixel is readily available in the above
  2071. algorithm and it has a use. Let us store all the z coordinate
  2072. values in a matrix Z[a][b]. If we then want to add a flat image
  2073. with constant z value (say z=k) to the already rendered scene, we
  2074. can do it very easily by writing only to pixels for which
  2075. Z[a][b]>k. Thus, we can use another rendering technique along
  2076. with depth-field rendering. This is called the z-buffer
  2077. algorithm. Accelerating calculations for the z-buffer with
  2078. polygons, you might want to remove the divide-per-pixel. That is
  2079. fully discussed towards the end of the texture-mapping chapter.
  2080.  
  2081. A few optimizations are easily implemented. Let us assume we are
  2082. drawing from the bottom of the screen to the top, going in
  2083. columns instead of rows. First, if you are not banked (or banked
  2084. less than the steepest slope, e.g. no «overhangs» are seen in the
  2085. picture) and you have rendered pixel a. Then you know how far
  2086. pixel a is from you. Then, then next pixel up, pixel b, will be
  2087. at least as far as pixel a. Thus, you do not need to fire a ray
  2088. starting from your eye to pixel b, you need only to start as far
  2089. as pixel a was. In fact, if you are not banked, you need only to
  2090. move your Pz coordinate up a little and keep going.
  2091.  
  2092. Therefore, if pixel a is determined to be a background pixel, so
  2093. will all the pixels above it.
  2094.  
  2095. You might want to render all of you polygons and other entities
  2096. using a z-buffer algorithm. As an example, say you have 3
  2097. polygons, A, B and C, and a sphere S. You could tell your
  2098. renderer to render poly A, which it would do right away, filling
  2099. the z-buffer at the same time. That is, for each screen pixel
  2100. that the buffer A covers, it would calculate the z-value of the
  2101. polygon in that pixel, then check with the z-buffer if it stands
  2102. in front of whatever is already drawn. If so, it would draw the
  2103. pixel and store the polygon's z-value for that pixel in the z-
  2104. buffer. Then render sphere S, which it would do right away again,
  2105. calculating the z-value for each pixel and displying only the
  2106. visible pixels using the z-buffer, and updating the z-buffer
  2107. accordingly. Lastly polygons B and C would be displayed.
  2108.  
  2109. What you gain by using this algorithm instead of a scan line
  2110. algorithm is obvious. Interpenetrating polygons behave correctly,
  2111. mixing rendering techniques is straightforward. There is but one
  2112. annoying thing: calculating the z value for a polygon in any
  2113. point on the screen involves a division, hence a division per
  2114. pixel if using z-buffer. But we can do better.
  2115.  
  2116. We are not really interested in the z value of a polygon in a
  2117. determined pixel, what we want really is to know what is visible.
  2118. To that purpose, we can obviously calculate z and show the object
  2119. with the smallest z value. Or we could evaluate z^2 and show the
  2120. one with the smallest z^2 value. If we want whatever has the
  2121. smallest z value, this is the same as saying that we want
  2122. whatever has the largest 1/z value. Fortunately 1/z is way easier
  2123. to calculate. Let the plane equation be Ax+By+Cz=D, and the
  2124. projection equations be u=x/z and v=y/z, or
  2125.  
  2126. Ax+By+Cz=D
  2127. x=uz
  2128. y=vz
  2129.  
  2130. or
  2131.  
  2132. A(uz)+B(vz+)Cz=D
  2133.  
  2134. or
  2135.  
  2136. z(Au+Bv+C)=D
  2137.  
  2138. (Au+Bv+C)/D=1/z
  2139.  
  2140. 1/z= (A/D)u+(B/D)v+(C/D)
  2141.  
  2142. Let M=A/D (a constant), N=B/D (another constant) and K=C/D (a
  2143. third constant), then
  2144.  
  2145. 1/z=Mu+Nv+K
  2146.  
  2147. As we can see plainly, we don't have a divide. This is linear so
  2148. can be calculated incrementally. E.g. for a given scanline (v=k),
  2149. the equation is
  2150.  
  2151. 1/z=Mu+Nk+K
  2152.  
  2153. Let P=Nk+K (a constant)
  2154.  
  2155. 1/z=Mu+P
  2156.  
  2157. A linear equation. When going from pixel u to pixel u+1, we just
  2158. need to add M to the value of 1/z. A single add per pixel.
  2159.  
  2160. There's also a double edged optimization that you can make, which
  2161. follows.
  2162.  
  2163. Sort your polygons in generally increasing z order (of course,
  2164. some polygons will have an overlap in z, but that's everybody's
  2165. problem with this algorithm - your solution is just as good as
  2166. anybody else's.) Now, you're scan-converting a polygon. If you
  2167. can find a span (e.g. a part of the current scanline) for which
  2168. both ends are behind the same convex object, then that span is
  2169. wholly hidden. As an example of this, if both enpoints of the
  2170. current span (on the scanline) are behind the same convex object,
  2171. then the current span is entirely hidden by that convex object.
  2172. If all of your objects are convex polygons, then you don't have
  2173. to check for convexity of the object. Another interesting example
  2174. is this:
  2175.  
  2176. If current span's left endpoint is hidden
  2177.      let A be the object that hides the left endpoint (of the
  2178. span)
  2179.      if A hides the right endpoint
  2180.           stop drawing the scanline
  2181.      else
  2182.           start drawing from the right end of the span,
  2183. leftwards, until object A is in front of the span
  2184.      end if
  2185. else
  2186.      if current span's right endpoint is hidden
  2187.           let B be the object that hides the right endpoint
  2188.           scan rightwards until object B is in front of the span
  2189.      else
  2190.           draw span normally
  2191.      end if
  2192. end if
  2193.  
  2194. This can be applied also to polyhedras if they are convex. This
  2195. last optimization should be relatively benign (e.g. I expect a
  2196. factor of two or some other constant from it).
  2197.  
  2198.  
  2199.  
  2200.  
  2201.  
  2202.      _____________________
  2203.  
  2204.  
  2205.  
  2206.  
  2207.  
  2208.  
  2209. 11   Bitmap scaling and mixing rendering techniques
  2210.  
  2211.  
  2212.  
  2213.  
  2214. Another popular technique is called bitmap (or more properly
  2215. pixmap) scaling. Say you are looking in a parallel direction to
  2216. the z axis, and you are looking at a picture with constant z.
  2217. Through the projection, you will realize that all that happens to
  2218. it is that it is scaled by a constant factor of 1/z. That is, if
  2219. you look at a photograph lying in the plane z=2, it will look
  2220. exactly the same (once projected) as the same photograph at z=1,
  2221. except that it will be twice as small. You could also do that
  2222. with the picture of a tree. From a great distance, the
  2223. approximation will be fairly good. When you get closer, you might
  2224. realize that the tree looks flat as you move. But nonetheless, it
  2225. can yield impressive results. Bitmap scaling often suffices for
  2226. ball-like objects. You could use it for an explosion, smoke,
  2227. spherical projectile, etc... The advantage of bitmap scaling is
  2228. that it has only one point, probably it's center, to transform
  2229. (i.e, translate and rotate into place in 3d), and scaling it is
  2230. very fast.
  2231.  
  2232. The advantages of such a technique are obvious. First, it is very
  2233. easy to scale a pixmap by a given factor. Second, the z value is
  2234. a constant, thus it can easily be incorporated in the above z-
  2235. buffer algorithm with depth-field rendering for very nice
  2236. results. However, it will be obvious to the careful observer that
  2237. you can not smoothly go around it, as it is always facing you the
  2238. same way (i.e. you cannot see the back door to a house if it is
  2239. rendered this way, because, obviously, the bitmap being scaled
  2240. does not have a back door).
  2241.  
  2242. Moreover, the depth-field rendering technique blends quite
  2243. smoothly with vector graphics. First, render the depth field in a
  2244. buffer B. Do a scan-line algorithm for the vector graphics, but
  2245. with an added twist: the «background» bitmap is B, and use the z-
  2246. buffer algorithm to determine if a pixel in B lies in front of
  2247. the current polygon if any.
  2248.  
  2249. If you do decide to use depth-field rendering along with vector
  2250. graphics, and you decide not to bank your eye, then you can
  2251. eliminate one full rotation from all calculations. E.g. if your
  2252. eye is not banked, you do not need to rotate around the z axis.
  2253.  
  2254.  
  2255.  
  2256.  
  2257.      _________________________
  2258.  
  2259.  
  2260.  
  2261.  
  2262.  
  2263. Scaling a bitmap can be done fairly easily this way. The
  2264. transformation is like this: we want to take every (u,v) point
  2265. and send it into (u/s,v/s), where s is the scaling factor (s=2
  2266. means that we are halving the size of the bitmap). The process is
  2267. linear.
  2268.  
  2269. (i,j)=(u/s,v/s)
  2270. (is,js)=(u,v)
  2271.  
  2272. Thus, if we are at pixel (i,j) on screen, it corresponds to pixel
  2273. (is,js) in the bitmap. This can also be done incrementally. If
  2274. un=A, then un+1=A+s. No multiplications are involved. Thus, if we
  2275. are viewing a bitmap at z=2 according to the above said 3d
  2276. projections, then we can render it using a scan-line algorithm
  2277. with these equations and s=z=2.
  2278.  
  2279. If a bitmap has holes in it, we can precalculate continuous spans
  2280. in it; e.g., in an O, drawing a scanline through the center shows
  2281. you that you are drawing two series of pixels. If you want the
  2282. center to be transparent (allow a background to show through),
  2283. then you can divide each scanline of the O in two spans, one left
  2284. and one right, and render them as two separate bitmaps; e.g.,
  2285. break the O in ( and ) parts, both of which have no holes. Note
  2286. however that we have a multiplication at each beggining of a
  2287. span, so this might not be efficient. Or you could select a color
  2288. number as being the «transparent» color, and do a compare per
  2289. pixel and not blit if it is the transparent color. Or you can
  2290. allocate an extra bitplane (a field of bit, one bit per pixel)
  2291. where each bit is either 1 for opaque (draw this pixel) or 0 for
  2292. transparent (do not draw this pixel).
  2293.  
  2294.  
  2295.  
  2296.  
  2297.  
  2298.      ______________________
  2299.  
  2300.  
  2301.  
  2302.  
  2303.  
  2304. Partially translucid primitives can be rendered pretty easily.
  2305. This is an old trick, you see. You can even scavenge the lookup
  2306. table you used for anti-aliasing. If you want to have a red glass
  2307. on a part of the screen, the lookup table will tell you how to
  2308. redden a color. That is, the item #c in the lookup table tells
  2309. you what color is color c with a little red added. The lookup
  2310. table need not be calculated on the fly, of course. It can be
  2311. calculated as part of your initialization or whatever.
  2312.  
  2313.  
  2314.  
  2315.  
  2316.      ______________________
  2317.  
  2318.  
  2319.  
  2320.  
  2321.  
  2322.  
  2323. 12   About automatically generating correctly oriented normals to
  2324. planes.
  2325.  
  2326.  
  2327.  
  2328. This bit is not meant to be a blindingly fast algorithm for
  2329. region detection. It is assumed that these calculations are made
  2330. before we start generating pictures. Our objective is to generate
  2331. a normal for each plane in a polyhedra oriented so that it points
  2332. outwards respectively to the surface. We assume all surfaces are
  2333. defined using polygons. But before determining normals for
  2334. planes, we will need to examine polygons in 3d.
  2335.  
  2336. Each polygon being constituted of a set of edges, we will later
  2337. require a vector for each edge that points in the polygon for any
  2338. edge. See below.
  2339.  
  2340.  
  2341.  
  2342. Figure 8
  2343.  
  2344. In figure 8, we can see in red the said vectors for all edges. As
  2345. we can see, the vectors do not necessarely point towards the
  2346. general center of mass. We could use the algorithm we used for
  2347. rendering polygons in two dimensions to determine the direction
  2348. of the arrows, but that is a bit complicated and involves trying
  2349. to find a direction that is not parallel to any of the above
  2350. lines. Thus, we will try another approach.
  2351.  
  2352. First, pick any point in the polygon. Then, from that point, take
  2353. the farthest point from it. See below:
  2354.  
  2355.  
  2356.  
  2357. Figure 9
  2358.  
  2359. (Please, do not measure it with a ruler :-) )
  2360. Now say we started with the point marked by the big red square.
  2361. Now, both point a and b are the farthest points from the red
  2362. square. You can pick any of the two. The objective here is to
  2363. find three consecutive points around the polygon for which the
  2364. inside angle is less than 180 degrees. If all three are on the
  2365. circle (worst case), their angle is strictly smaller than 180
  2366. degrees (unless they are all the same point, silly).
  2367.  
  2368. You could also pick the topmost point. If several are at the same
  2369. height, use the leftmost (or rightmost, whatever) one in the
  2370. topmost points.
  2371.  
  2372. Now we are certain that the small angle is the inside angle. Let
  2373. us draw the general situation:
  2374.  
  2375.  
  2376.  
  2377. Figure 10
  2378.  
  2379. Now, in a we see the above mentionned farthest point and the two
  2380. adjacent edges. If we take these edges as vectors as in b and
  2381. then sum them as in c, we find a point (pointed by the red vector
  2382. in c) which is necessarely on the in side for both edges. Now it
  2383. is easy to find a perpendicular vector for both edges and make
  2384. sure that they are in the good general direction. A way would be
  2385. to take the perpendicular vector and do a scalar multiplication
  2386. with the above red vector, the result of which must be strictly
  2387. positive. If it is not, the you must take the opposite vector for
  2388. perpendicular vector.
  2389.  
  2390. Now we know the situation for at least one edge. We need some way
  2391. to deduce the other edges' status from that edge's.
  2392.  
  2393. If we take Figure 10.c as our general case, the point at the end
  2394. of the red vector is not necessarely on the in side for both
  2395. edges. However, if it is on the in side for any of the two edges,
  2396. it is on the in side for both, and if it is on the out side for
  2397. any edge, it is on the out side for both.
  2398.  
  2399.  
  2400. Normally, we will never find two colinear edges in a polygon, but
  2401. if we do, it is fairly easy to solve it anyway (left as an
  2402. exercise).
  2403.  
  2404.  
  2405.  
  2406. Now, briefly, here is how to determine if a normal to a plane
  2407. points in the right direction. Pick the point with the greatest z
  2408. value (or x, or y). If many points have the same z value, break
  2409. ties with greatest x value, then greatest y value. Then, take all
  2410. connected edges. Make them into vectors, pointing from the above
  2411. point to the other endpoints (that is, the z component is less
  2412. than or equal to zero). Make them unit vectors, that is, length
  2413. of one. Find the vector with the most positive (or less negative)
  2414. z value. Break ties with the greates x value, and then greatest y
  2415. value. Take its associated edge. Now that edge is said to be
  2416. outermost.
  2417.  
  2418. Now, from that outermost edge, take the two associated polygons.
  2419. The edge has an associated vector for each connected polygon to
  2420. indicate the "in" direction. Starting from any of the two points
  2421. of the edge, add both these vectors. The resulting point is
  2422. inside. Adjust the normal so that it is on the right side of the
  2423. plane for both polygons.
  2424.  
  2425. Now we have the normal for at least one plane straight. We will
  2426. deduce correct side for the normal to the other polygons in a
  2427. similar manner than we did for the edges of a polygon, above.
  2428.  
  2429. Take two adjacent polygons, one for which you do know the normal,
  2430. the other for which you want to calculate the correct normal. Now
  2431. take a shared edge (any edge will do). That edge has two "in"
  2432. vectors for the two polygons. Starting from any of the two points
  2433. of the edge, add the two «in» vectors. If the resulting point is
  2434. on the «in» side for a plane, it is also on the «in» side of the
  2435. other plane. If it is on the «out» side of a plane, it is on the
  2436. «out» side for both planes.
  2437.  
  2438.  
  2439.  
  2440. To get a general picture of what I'm trying to say, here's a
  2441. drawing. Basically, where two lines intersect, you get four
  2442. distinct regions. One region lies entirely on the «in» side of
  2443. both lines and another region lies entirely on the «out» side of
  2444. both lines (if the two lines define a polygon). The two other
  2445. regions are mixed region (they lie on the «in» side of a line and
  2446. on the «out» side of the other line). The «mixed» regions may or
  2447. may not be part of the actual polygon, depending on whether the
  2448. angle is lesser or greater than 180 degrees.
  2449.  
  2450.  
  2451.  
  2452. Figure 11
  2453.  
  2454.  
  2455. So we are trying to find a point that is in the «in» or «out»
  2456. region, but not in the «mixed» region. If, as above, we take the
  2457. edges as vectors and add the together, we end up in the either
  2458. the «in» region or the «out» region, but not in the «mixed»
  2459. region. That can be proved very easily, geometrically (yet
  2460. strictly), but I don't feel like drawing it. :) Anyway, it should
  2461. be simple enough.
  2462.  
  2463.  
  2464.  
  2465.  
  2466.      __________________
  2467.  
  2468.  
  2469.  
  2470.  
  2471.  
  2472. 13   Reducing a polygon to a mesh of triangles
  2473.  
  2474.  
  2475.  
  2476. Now that's fairly simple. The simplest is when you have a convex
  2477. polygon. Pick any vertex. Now, take its two adjacent vertexes.
  2478. That's a triangle. In fact, that's your first triangle in your
  2479. mesh of triangles. Remove it from your polygon. Now you have a
  2480. smaller polygon. Repeat the same operation again until you are
  2481. left with nothing. Example:
  2482.  
  2483.  
  2484. Figure 12
  2485.  
  2486. For the above polygon, we pick the green edges in 1. The triangle
  2487. is shown in 2. When we remove it from the polygon, we're left
  2488. with what we see in 3.
  2489.  
  2490. This will take a n sided polygon and transform into n-2
  2491. triangles.
  2492.  
  2493. If the polygon is concave, we can subdivide it until it is convex
  2494. and here is how we do that.
  2495.  
  2496. Take any vertex where the polygon is concave (i.e. the angle at
  2497. that vertex is more than 180 degrees) and call that vertex A.
  2498. From that vertex, it can be proved geometrically that there
  2499. exists another vertex that is not connected by an edge to which
  2500. you can extend a line segment without intersecting any of the
  2501. polygon's edges. In short, from vertex A, find another vertex,
  2502. not immediately connected by an edge (there are two vertices
  2503. connected to A by an edge, don't pick any one of them). Call that
  2504. new vertex B. Make certain that AB does not intersect any of the
  2505. other edges. If it does, pick another B. Once you're settled on
  2506. your choice of B, split your polygon in two smaller polygons by
  2507. inserting edge AB. Repeat the operation with the two smaller
  2508. polygons until you either end up with triangles or end up with
  2509. convex polygons which you can split as above. It can again be
  2510. proved (though a little more difficultly) that you will end up
  2511. with n-2 triangles if you had a n sided polygon.
  2512.  
  2513. Here's a tip on how to find a locally concave vertex.
  2514.  
  2515.  
  2516. Figure 13
  2517.  
  2518. [The arrows point towards the «in» side]
  2519.  
  2520. The two possible cases are shown in figure 13, either it is or it
  2521. is not convex. In 1, it is convex, as in 2 it is concave. To
  2522. determine if it is or not convex, take the line that passes
  2523. through a and u, above. Now, if v stands on the «in» side of the
  2524. au line, it is convex. Otherwise, it is concave.
  2525.  
  2526.  
  2527.  
  2528.  
  2529.  
  2530.      __________________
  2531.  
  2532.  
  2533.  
  2534.  
  2535.  
  2536. 14   Gouraud shading
  2537.  
  2538.  
  2539.  
  2540. First, let us discuss a simple way for shading polygons. Imagine
  2541. we have a light source that's infinitely far. Imagine it's at
  2542. z=positive infinity, x=0, y=0. If it's not, you can rotate
  2543. everything so that it is (in practice, you'll only need to rotate
  2544. normals). Now, the z component of the polygons normals give us a
  2545. clue as to the angle between the polygon and the light rays. The
  2546. z component goes from -1, facing away from the light source
  2547. (should be darkest) to 0, facing perpendicular to the light
  2548. source (should be in half light/penumbra), to +1, facing the
  2549. light source (should be the brightest). From that, you can shade
  2550. it linearly or any way you want.
  2551.  
  2552. With this model, a polygon is uniformly shaded. Approximating
  2553. round objects through polygons will always look rough.
  2554.  
  2555. Gouraud shading is what we call an interpolated shading. It is
  2556. not physically correct, but it looks good. What it does is this:
  2557. it makes the shade of a polygon change smoothly to match the
  2558. shade of the adjacent polygon. This makes the polyhedra look
  2559. smooth and curved. However, it does not change the actual edges
  2560. to curves, thus the silhouette will remain unsmoothed. You can
  2561. help that by allowing curved edges if you feel like it.
  2562.  
  2563. First, we will interpolate a vertex normal for all vertexes. That
  2564. is, we will give each vertex a «normal» (note that a point does
  2565. not have a normal, so you could call it pseudo-normal). To
  2566. interpolate that vertex «normal», averaging all connected
  2567. polygons normals would be a way. Now, find the «shade» for each
  2568. vertex.
  2569.  
  2570.  
  2571. Figure 14
  2572.  
  2573. Now, say point a is at (ax,ay) and point b is at (bx,by). We're
  2574. scan-converting the polygon, and the gray-and-red line represents
  2575. the current sanline, y0. Point P is the pixel we are currently
  2576. drawing, located at (x0,y0). Points m and n define the span we're
  2577. currently drawing. Now, we will interpolate the color for point m
  2578. using a and b, and the color for n using a and c, then
  2579. interpolate the color for P using the color for m and n. Our
  2580. interpolation will be linear. Say color for point a is Ca, color
  2581. for point b is Cb, for point c it is Cc, Cm for m, Cn for n and
  2582. CP for P.
  2583.  
  2584. We will say that color at point m, Cm, is as follows:
  2585.  
  2586. Cm=(y0-ay)/(by-ay) x (Cb-Ca)+Ca
  2587.  
  2588. That is, is point m is halfway between a and b, we'll use half ca
  2589. and half cb. If it's two-third the way towards b, we'll use 1/3Ca
  2590. and 2/3Cb. You get the picture. Same for n:
  2591.  
  2592. Cn=(y0-ay)/(cy-ay) x (Cc-Ca)+Ca.
  2593.  
  2594. Then we do the same for CP.
  2595.  
  2596. CP=(x0-mx)/(nx-mx) x (Cn-Cm)+Cn.
  2597.  
  2598. If you think a little, these calculations are exactly the same as
  2599. Bresenham's line drawing algorithm, seen previously in chapter
  2600. 2.1. It is thus possible to do them incrementally. Example:
  2601.  
  2602. Say we start with the topmost scanline. Color for point m is at
  2603. first Cb. Then, it will change linearly. When point m reaches
  2604. point a, it will be color Ca. Now, say dy=by-ay, dC=Cb-Ca, and
  2605. u=y0-ay.
  2606.  
  2607. Cm=dC/dy x u+Ca.
  2608.  
  2609. Then, when y0 increases by 1, u increases by one and we get
  2610.  
  2611. Cm'=dC/dy x (u+1)+Ca=dC/dy x u+Ca + dC/dy
  2612. Cm'=Cm+dC/dy
  2613.  
  2614. Same as Bresenham. Thus, when going from one scanline to the
  2615. next, we simply need to add dC/dy to Cm, no multiplications or
  2616. divisions are involved. The same goes for Cn. CP is done
  2617. incrementally from pixel to pixel.
  2618.  
  2619.  
  2620.  
  2621.  
  2622.  
  2623.  
  2624.      ____________________
  2625.  
  2626.  
  2627.  
  2628.  
  2629.  
  2630.  
  2631.  
  2632. 15   Clipping polygons to planes
  2633.  
  2634.  
  2635.  
  2636. Eventually, you might need to clip your polygon, minimally to the
  2637. z>0 volume. Several approaches can be used. We will discuss here
  2638. the Sutherlan-Hodgman algorithm. But first, let us speak of
  2639. trivial rejection/acceptation.
  2640.  
  2641. If a polygon does not intersect the clipping plane, it can be
  2642. either trivially accepted or rejected. For example, if every
  2643. single point in a polygon are in z>0 and we are clipping with
  2644. z>0, then the whole polygon can be accepted. If every single
  2645. point is in z<=0, the whole polygon can be trivially rejected.
  2646. The real problem comes with cases where some points are in z>0
  2647. and some are in z<=0. Another nifty thing you can do is pre-
  2648. calculate the bounding sphere of a poly with its center on a
  2649. given point in the poly. Let O be the sphere's center and R it's
  2650. radius. If Oz-R>0, then the poly can be trivially accepted even
  2651. faster (no need to check every single point). If Oz+R<=0, the
  2652. poly can be trivially refused. You can extend this to polyhedras.
  2653. You could also check wether the intersection of the sphere and
  2654. the polygon plane is in z>0 (or z<=0), which might be slightly
  2655. better than checking for the whole sphere.
  2656.  
  2657. Here comes the Sutherland-Hodgman algorithm. Start with a point
  2658. that you know is going to be accepted. Now, move clockwise (or
  2659. counter-clockwise) around the poly, accepting all edges that
  2660. should be trivially accepted. Now, when an edge crosses from
  2661. acceptance to rejection (AtoR), find the intersection with the
  2662. clipping plane and replace the offending edge by a shortened one
  2663. that can be trivially accepted. Then, keep moving until you find
  2664. an edge that crosses from rejection to acceptance (RtoA), clip
  2665. the edge and keep the acceptable half. Add an edge between the
  2666. two edges AtoR and RtoA. Keep going until you are done.
  2667.  
  2668.  
  2669. Figure 15
  2670.  
  2671. Figure 15 illustrates the process of clipping the abcde polygon
  2672. to the clipping plane. Of note is this: when performing this kind
  2673. of clipping on convex polygons, the results are clear.
  2674. Furthermore, a convex polygon always produce a convex polygon
  2675. when clipped. However, concave polygons introduce the issue of
  2676. degenerate edges and whether they should be there in a correctly
  2677. clipped polygon. Figure 16 shows such a case of degenerate edge
  2678. generation.
  2679.  
  2680.  
  2681. Figure 16
  2682.  
  2683. In figure 16, the polygon to the left is clipped to the plane
  2684. shown, the result is on the right. The bold edge is double. That
  2685. is, two edges pass there. This is what we call a degenerate edge.
  2686. Degenerate edges don't matter if what interests us is the area of
  2687. the polygon, or if they happen to fall outside of the view
  2688. volume. What interests us is indeed the area of the polygon, but
  2689. due to roundoff error, we could end up with the faint trace of an
  2690. edge on screen. On could eliminate these degenerate edges pretty
  2691. easily.
  2692.  
  2693. To do so, do not add an edge between the points where you clip at
  2694. first (edge xy in figure 15). Instead, once you know all clipping
  2695. points (x and y in the above example), sort them in increasing or
  2696. decreasing order according to one of the coordinate (e.g. you
  2697. could sort the clipping point in increasing x). Then, between
  2698. first and second point, add an edge, between third and fourth,
  2699. add an edge, between fifth and sixth, add an edge and so on. This
  2700. is based on the same idea that we used for our polygon drawing
  2701. algorithm in an earlier chapter. That is, when you intersect an
  2702. edge, you either come out of the polygon or go in the polygon.
  2703. Then, between the first and second clipped point belongs an edge,
  2704. etc...
  2705.  
  2706. However, since we are clipping to z=0, the degenerate edge has to
  2707. be outside the view volume. (Division by a very small z, when
  2708. making projections, ensures that the projected edge is far from
  2709. the view window). Since we are clipping to z>0, let us pull out
  2710. our limit mathematics.
  2711.  
  2712. Let us first examine the case where the edge lies entirely in the
  2713. z=0 plane.
  2714.  
  2715. Let us assume we are displaying our polygon using a horizontal
  2716. scan-line algorithm. Now, if the both of the offending edge's
  2717. endpoints are at y>0, then the projected edge will end up far up
  2718. of the view window. Same goes if both endpoints are at y>0. If
  2719. they are on either side of y=0, then all projected points will
  2720. end up above or below the view window, except the exact point on
  2721. the edge where y=0. If that point's x<0, then x/z will end up far
  2722. to the left, and if x>0, the point will end up far to the right.
  2723. If x=0, then we have the case x=y=0 and z tends towards positive
  2724. zero (don't flame me for that choice of words). In the case x=0,
  2725. x/z will be zero (because 0/a for any nonzero a is 0) and so for
  2726. y. Thus, the projected edge should pass through the center of the
  2727. view window, and go to infinity in both directions. It's slope
  2728. should be the same as the one it has in the z=0 plane. Just to
  2729. make it simple, if we have the plane equation (as we should) in
  2730. the form Ax+By+Cz+D=0, we can find the projected slope by fixing
  2731. z=1: Ax+By+C+D=0, or y=-A/Bx-(C+D). Thus the slope of the
  2732. projected line is -A/B. If B is zero, we have a vertical line.
  2733.  
  2734. What we conclude is that if the edge's z coordinate is 0 and both
  2735. of its endpoints are on the same side of plane y=0, then the edge
  2736. will not intercept any horizontal scanline and can thus be fully
  2737. ignored. If the endpoints are on either side of y=0, then it will
  2738. have an intersection with all horizontal scanlines. The projected
  2739. edge will end up left if, when evaluated at y=0 and z=0 its x
  2740. coordinate is less than 0 (from the plane equation, x=-D/A). If
  2741. the x coordinate is positive, the edge ends up far to the right.
  2742. If x=y=0 when z=0, the edge will be projected to a line passing
  2743. through the center of the screen with a slope of -A/B.
  2744.  
  2745. In the case where one of the edge's endpoints is in the z=0
  2746. plane, that endpoint will be projected to infinity and the other
  2747. will not, we will end up with a half line in the general case,
  2748. sometimes a line segment.
  2749.  
  2750. Let us name the endpoint whose z coordinate is 0 P0. If P0's x<0,
  2751. then it will be projected to the left of the window. If P0's x>0,
  2752. it will be projected to the right. Likewise, if P0's y<0, it will
  2753. be projected upwards, if y>0 it will be projected downwards.
  2754.  
  2755. The slope will be (v0-v1)/(u0-u1) where (u0,v0) is projected P0
  2756. and (u1,v1) is projected P1. This can be written as
  2757. u0=P0x/p0z     u1=P1x/P1z     v0=P0y/P0z     v1=P1y/P1z
  2758.  
  2759. It is of note that u1 and v1 are both relatively small numbers
  2760. compared to u0 and v0. The slope is therefore:
  2761.  
  2762. m=(P0y/P0z-P1y/P1z)/(P0x/P0z-P1x/P1z)
  2763. m=(P0y/P0z)/(P0x/P0z-P1x/P1z)-(P1y/P1z)/(P0x/P0z-P1x/P1z)
  2764. m=(P0y/P0z)/([P0xP1z-P1xP0z]/P0zP1z)-0
  2765. m=(P0y/P0z)(P0zP1z/[P0xP1z-0])
  2766. m=P0yP1z/(P0xP1z)
  2767. m=P0y/P0x.
  2768.  
  2769. By simmetry we could have found the inverse slope p=P0x/P0y.
  2770. Thus, in the case where P0x=0, we are facing a vertical
  2771. projection. Else, the projection's slope is P0y/P0x. Anyway, the
  2772. line passes through (u1,v1).
  2773.  
  2774. In the case where P0=0, (u0,v0)=(0,0). Then the projected edge
  2775. remains a line segment.
  2776.  
  2777.  
  2778.  
  2779.  
  2780.  
  2781.      _________________________
  2782.  
  2783.  
  2784.  
  2785.  
  2786.  
  2787.  
  2788. 16   Interpolations and forward differencing
  2789.  
  2790.  
  2791.  
  2792.  
  2793. Oftentimes, us graphics programmer want to interpolate stuff,
  2794. especially when the calculations are very complex. For example,
  2795. the texture mapping equations involve at least two divides per
  2796. pixel, which can be just a bit too much for small machines. Most
  2797. of the time, we have a scanline v=k on screen that we are
  2798. displaying, and we want to display a span of pixels, e.g. a
  2799. scanline in a triangle or something. The calculations for these
  2800. pixels can sometimes be complicated. If, for example, we have to
  2801. draw pixels in the span from u=2 to u=16 (a 14 pixel wide span)
  2802. for a given scanline, we would appreciate reducing the per-pixel
  2803. calculations to the minimum.
  2804.  
  2805. Let's take for our example, the z value of a plane in a given
  2806. pixel. Let's say the plane equation is x+2y+z=1. We have already
  2807. seen that z can be calculated with z=D/(Au+Bv+C), or
  2808. z=1/(u+2v+1). Let us imagine that we are rendering a polygon, and
  2809. that the current scanline is v=1. Then, z=1/(u+3) for the current
  2810. scanline. If that a span goes from u=0 to u=20 (that is, on line
  2811. v=.5, the polygon covers pixels u=0 through u=20). We can see
  2812. form the equation for z that when u=0, z=1/3, and when u=20,
  2813. z=1/23.
  2814.  
  2815. Calculating the other values is a bit complex because of the
  2816. division. (Divisions on computers are typically slower than other
  2817. operations.) As of such, we could perhaps use an approximation of
  2818. z instead. The most naive one would be the following: just assume
  2819. z is really of the form z=mu+b and passes through the points u=0,
  2820. z=1/3 and u=20, z=1/23. Thus:
  2821.  
  2822. z=mu+b verifies
  2823. 1/3=m(0)+b and
  2824. 1/23=m(20)+b
  2825.  
  2826. thus
  2827.  
  2828. 1/3=b
  2829. and
  2830. 1/23=20m+b
  2831. or
  2832. 1/23=20m+1/3
  2833. 1/23-1/3=20m
  2834. m=1/460-1/60 or approximately -0.01449275362319.
  2835.  
  2836. Thus, instead of using z=1/(u+3), we could use z=-
  2837. 0.01449275362319u+1/3. In this particular case, this would be
  2838. somewhat accurate because m is small.
  2839.  
  2840. Example, we know from the real equation that z=1/(u+3), thus when
  2841. u=10, z=1/13, approximately 0.07692307692308. From the
  2842. interpolated equation, when u=10, we get z=-
  2843. 0.01449275362319*10+1/3, or approximately 0.1884057971014. The
  2844. absolute error is approximately 0.11, which may or may not be
  2845. large according to the units used for z. However, the relative
  2846. error is around 30%, which is quite high, and it is possible to
  2847. create cases where the absolute error is quite drastic.
  2848.  
  2849. To reduce the error, we could use a higher degree polynomial.
  2850. E.g. instead of using z=mx+b, we could for instance use
  2851. z=Ax^2+Bx+C, a second degree polynomial. We would expect this
  2852. polynomial to reduce the error. The only problem with this and
  2853. higher degree polynomials is to generate the coefficients (A, B
  2854. and C in this case) which will minimize the error. Even defining
  2855. what "minimizing the error" is can be troublesome, and depending
  2856. on the exact definition, we will find different coefficients.
  2857.  
  2858. A nice way of "minimizing the error" is using a Taylor series.
  2859. Since this is not a calculus document, I will not teach it in
  2860. detail, but merely remind you of the sum. It can be demonstrated
  2861. that all continuous functions can be expressed as a sum of
  2862. polynomials for a given range. We can translate that function so
  2863. as to center that range about any real, and such other things. In
  2864. brief, the Taylor series of f(x) around a is:
  2865.  
  2866. T=f(a)+f'(a)(x-a)/1!+f''(a)(x-a)^2/2!+f'''(a)(x-
  2867. a)^3/3!+f''''(a)(x-a)^4/4!+...
  2868.  
  2869. This may look a bit scary, so here is the form of the general
  2870. term:
  2871.  
  2872. Ti=fi(a)(x-a)^i/i!
  2873.  
  2874. and
  2875.  
  2876. T=T0+T1+T2+T3+...+Ti+...
  2877.  
  2878. Where fi(a) is the ith derivative of f evaluated at point a, and
  2879. i! is the factorial of i, or 1x2x3x4x....xi. As an example, the
  2880. factorial of 4 is 1x2x3x4=24. It is possible to study this serie
  2881. and determine for which values of x it converges and such things,
  2882. and it is very interesting to note that for values of x for which
  2883. it does not converge, it diverges badly usually.
  2884.  
  2885. Taylor series usually do a good job of converging relatively
  2886. quickly in general. For particular cases, one can usually find
  2887. better solutions than a Taylor series, but in the general case,
  2888. they are quite handy. Saying that they converge "relatively
  2889. quickly" means that you could evaluate, say, the first 10 terms
  2890. and be quite confident that the error is small. There are
  2891. exceptions and special cases, of course, but this is generally
  2892. true.
  2893.  
  2894. What interests us in Taylor series is that they give us a very
  2895. convenient way of generating a polynomial of any degree to
  2896. approximate any given function for a certain range. In
  2897. particular, the z=1/(mu+b) equation could be approximated with a
  2898. Taylor series, and it's radius of convergence determined, etc...
  2899. (Please look carefully at the taylor series and you will see that
  2900. it is indeed a polynomial).
  2901.  
  2902. Here are a few well known and very used Taylor series expansions:
  2903.  
  2904. exp(x)=1+x+x^2/2+x^3/3!+x^4/4!+x^5/5!+x^6/6!+....
  2905. sin(x)=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!-...
  2906. cos(x)=1-x^2/2!+x^4/4!-x^6/6!+x^8/8!-...
  2907.  
  2908. Taylor series for multivariable functions also exist, but they
  2909. will not be discussed here. For a good description of these, you
  2910. can get a good advanced calculus book. My personnal reference is
  2911. Advanced engineering mathematics, seventh edition by Erwin
  2912. Kreysig, ISBN 0-471-55380-8. As for derivatives, limits et al,
  2913. grab an introductory level calculus book.
  2914.  
  2915. Now the only problem is evalutating a polynomial at any point. As
  2916. you know, evaluating Ax^2+Bx+C at point x0 requires that we
  2917. square x0, multiply it by A, add B times x0 and then finally add
  2918. C. If we're doing that once per pixel, we are losing a lot of
  2919. time. However, the same idea as in incremental calculations can
  2920. be recuperated and used for any degree of polynomial. Generally
  2921. speaking, we're going to try to find what happens to f(x) as we
  2922. increase x by one. Or, we're trying to find f(x+1)-f(x) to see
  2923. the difference between one pixel and the next. Then, all we need
  2924. to do is add that value as we move from one pixel to the next.
  2925. Let's do an example with f(x)=mx+b.
  2926.  
  2927. f(x+1)-f(x)=[m(x+1)+b]-[mx+b]
  2928. f(x+1)-f(x)=m
  2929.  
  2930. So, as we move from a pixel to the next, we just need to add m to
  2931. the previous value of f(x).
  2932.  
  2933. If we have a second order polynomial, the calculations are
  2934. similar:
  2935.  
  2936. f(x)=Ax^2+Bx+C
  2937.  
  2938. f(x+1)-f(x)=[A(x+1)^2+B(x+1)+C]-[Ax^2+Bx+C]
  2939.      =[A(x^2+2x+1)+Bx+B+C]-[Ax^2+Bx+C]
  2940.      =[Ax^2+(2A+B)x+A+B+C]-[Ax^2+Bx+C]
  2941.      =2Ax+A+B
  2942.  
  2943. Let's name that last equation g(x). So, as we move from x to x+1,
  2944. we just need to add g(x) to the value of f(x). However,
  2945. calculating g(x) involves two multiplications (one of which which
  2946. can be optimized out by using bit shifts) and at least an add.
  2947. But g(x) is a linear equation, so we can apply forward
  2948. differencing to it again and calculate:
  2949.  
  2950. g(x+1)-g(x)=[2A(x+1)+A+B]-[2AX+A+B]
  2951.      =2A
  2952.  
  2953. So what we can do as we move from x to x+1 is first add g(x) to
  2954. f(x) and then add 2A to g(x), only two adds per pixel.
  2955.  
  2956. This can be extended to any degree of polynomial. In particular,
  2957. NURBS (not described in this document) can be optimized this way.
  2958. (I do not intend to discuss NURBS, but this optimization is
  2959. considered less efficient that subdivision, but is worth
  2960. mentionning.)
  2961.  
  2962. So what is the use for all this? Well, you could use a Taylor
  2963. series to do texture mapping instead of two divisions per pixel.
  2964. This is a generalization of linear approximation for texture
  2965. mapping. You could also use it for evaluating the z value of a
  2966. polygon in any pixel, as we were doing in the introduction to
  2967. this chapter. This however might not be very useful since you
  2968. might not need the actual z value, and the 1/z value might
  2969. suffice which, as we have already seen, can be calculated
  2970. incrementally without approximation (e.g. no error except
  2971. roundoff). This would be useful in visible surface determination.
  2972.  
  2973.  
  2974.  
  2975.  
  2976.      _________________________
  2977.  
  2978.  
  2979.  
  2980.  
  2981.  
  2982.  
  2983. 17   Specular highlights and the Phong illumination model
  2984.  
  2985.  
  2986.  
  2987.  
  2988.  
  2989. Up to now, we have been using a very straight illumination model.
  2990. We have assumed that the "quantity of light at a given point" can
  2991. be calculated with some form of dot product of the light vector
  2992. with the normal vector. (Let me remind you that (a,b,c) dot
  2993. (d,e,f) is ad+be+cf.) You should normalize the light vector for
  2994. this, and the plane normal should be of unit length too. Let us
  2995. assume the plane normal is (A,B,C) and that it is already
  2996. normalized. Let us further assume that the light vector (vector
  2997. from light source to point that is lighted in 3d) is (P,Q,R).
  2998. Then, the normalized light vector is (P,Q,R) times
  2999. 1/sqrt(P*P+Q*Q+R*R). As you can see, normalizing the light vector
  3000. is very annoying. (See chapter 2 for a more lengthy discussion of
  3001. related topics).
  3002.  
  3003. Now, we have been imagining that no matter what angle you look at
  3004. the object from, the intensity of the light is the same. In
  3005. reality, that may not be the case. If the object is somewhat
  3006. shiny and a bright light shines on it and you are looking at the
  3007. object from an appropriate angle, there should be a very intense
  3008. spot. If you move your eye, the spot moves, and if you move your
  3009. eye enough, the spot disappears entirely. (Try it with an apple
  3010. in a room with a single, intense light source, no mirrors et al).
  3011. This means that the light is reflected more in one general
  3012. direction than in the others.
  3013.  
  3014. Let us look at a little drawing. Let N be the surface normal, L
  3015. the light vector, R the reflection vector and V the eye-to-object
  3016. vector.
  3017.  
  3018.  
  3019. Figure 17
  3020.  
  3021. In the above picture, we see clearly that R is L mirrored by N.
  3022. Please not that the angle between L and R is not 90 degrees as it
  3023. may appear, but rather twice the angle between L and N. It can be
  3024. demonstrated that:
  3025.  
  3026. R=2N(N dot L)-L
  3027.  
  3028. Or, in more detail, if R=(A,B,C), N=(D,E,F) and L=(G,H,I), then N
  3029. dot L is DG+EH+FI, and thus we have
  3030.  
  3031. A=2D(DG+EH+FI)-G
  3032. B=2E(DG+EH+FI)-H
  3033. C=2F(DG+EH+FI)-I
  3034.  
  3035. It is of note that, if the light is at infinity, for a flat
  3036. polygon, N dot L is constant. If the light is not at infinity or
  3037. if the surface is curved, then N dot L varies.
  3038.  
  3039. Now, if we do not take into account the angle b in figure 17, the
  3040. amount of light perceived should be something like
  3041. I=Kcos(a)
  3042.  
  3043. where K is some constant that depends on the material and a is
  3044. the angle shown in figure 17. If to that we want to add the
  3045. specular highlight, we should add a term that depend on angle b
  3046. of figure 17. I becomes
  3047.  
  3048. I=Kcos(a)+Ccos^n(b)
  3049.  
  3050. Where C is a constant that depends on the material. The constant
  3051. n is a value that will specify how sharp is the reflection. The
  3052. higher the value for n, the more localized the reflection will
  3053. be. Lower values of n can be used for matte surfaces, and higher
  3054. values can be used for more reflective surfaces. Most light
  3055. sources don't shine as bright from a distance (but some do shine
  3056. as bright at any small distance, such as the sun, but that light
  3057. can be considered to lay "infinitely" far). In that case, we
  3058. could want to divide the above intensity by some function of the
  3059. distance. We know from physics that the quantity of energy at a
  3060. distance r of a light source is inversely proportional to the
  3061. square of the radius. However, this may not look good in computer
  3062. graphics. In short, you might want to use the square of the
  3063. length of the light vector, or the length of the light vector, or
  3064. some other formula. Let us assume we are using the length of the
  3065. light vector. Then, the intensity becomes:
  3066.  
  3067. I=1/|L| * [Kcos(a)+Ccos^n(b)]
  3068.  
  3069. Furthermore, it is known that if t is the angle between the
  3070. vectors A and B, then A dot B is |A||B|cos(t). Then,
  3071.  
  3072. cos(t)=(A dot B)/(|A||B|)
  3073.  
  3074. If |A| and |B| are both one, we can ignore them. So, let L' be
  3075. normalized L, and R' be normalized R and V' be normalized V, and
  3076. of course N is of unit length, then the equation is:
  3077.  
  3078. I=D * (K[N dot L']+C[R' dot V'])
  3079.  
  3080. Where D=1/|L|. We could have used some other function of |V|.
  3081.  
  3082. To that equation above, you may want to add some ambiant light,
  3083. e.g. light that shines from everywhere at the same time, such as
  3084. when you are standing outside (you wouldn't see a very dark
  3085. corner in daylight, normally, so the light that is still present
  3086. in a "dark corner" is the ambiant light). If you have multiple
  3087. light sources, just calculate the intensity from each light
  3088. source and add them together.
  3089.  
  3090. This of course assumes that a is comprised between 0 and 90
  3091. degrees, otherwise the light is on the wrong side of the surface
  3092. and doesn't affect it at all. (In which case the only light
  3093. present would be the light of intensity A).
  3094.  
  3095.  
  3096.  
  3097.  
  3098.  
  3099.      _________________________
  3100.  
  3101.  
  3102.  
  3103.  
  3104.  
  3105.  
  3106. 18   Phong shading
  3107.  
  3108.  
  3109.  
  3110.  
  3111.  
  3112. Polygons are flat, and as such have a constant normal across
  3113. their surface. When we used pseudo-normals for gouraud shading.
  3114. This was because we were attempting to model a curved surface
  3115. with flat polygons. We can use the Phong illumination model with
  3116. flat polygons if we want. That is, forget about the pseudo-
  3117. normals, use the plane normal everywhere in the Phong
  3118. illumination equation. This will yield very nice looking flat
  3119. surfaces, very realistic. However, most of the time, we want to
  3120. model curved surfaces and the polygons are merely an
  3121. approximation to that curved surface. In a curved surface, the
  3122. normal would not be constant for a large area, it would change
  3123. smoothly to reflect the curvature. In our polygons, of course,
  3124. the normal does not change. However, if we really want to model
  3125. curved surfaces with polygons, we find that the sharp contrasts
  3126. between facets is visually disturbing and does not represent a
  3127. curved surface well enough. To this end, we made the intensity
  3128. change gradually with Gouraud shading. There is however a problem
  3129. with Gouraud shading.
  3130.  
  3131. It is possible to demonstrate that no light inside a Gouraud
  3132. shaded polygon can be of higher intensity than any of the
  3133. vertices. As of such, if a highlight were to fall inside the
  3134. polygon (as is most probably the case) and not on a vertex, we
  3135. would miss it perhaps entierly. To this end, we might want to
  3136. interpolate the polygon normal instead of intensity from the
  3137. vertices pseudo-normals. The idea is the same as with Gouraud
  3138. shading, except that instead of calculating a pseudo-normal at
  3139. the vertices, calculating the intensities at the vertices and
  3140. then interpolating the intensity linearly, we will calculate
  3141. pseudo-normals at the vertices, interpolate linearly the pseudo-
  3142. normals and then calculate the intensity. The intensity can be
  3143. calculated using any illumination model, and in particular the
  3144. Phong illumination model can be visually pleasent, though some
  3145. surfaces don't have any specular highlights. (In which case you
  3146. can remove the specular factor from the equation entirely).
  3147.  
  3148. The calculations for Phong shading are very expensive per pixel,
  3149. in fact too expensive even for dedicated hardware oftentimes. It
  3150. would be absurd to think that one could do real-time true Phong
  3151. shading on today's platforms. But in the end, what is true Phong
  3152. shading? Nothing but an approximation to what would happen if the
  3153. polygon mesh really were meant to represent a curved surface.
  3154. I.e. Phong shading is an approximation to an approximation. The
  3155. principal objective of Phong shading is to yield nonlinear
  3156. falloff and not miss any specular highlights. Who is there to say
  3157. that no other function will achieve that?
  3158.  
  3159. One of the solutions to this has been proposed in SIGGRAPH '86,
  3160. included with this document. What they do is pretty
  3161. straightforward (but probably required a lot of hard work to
  3162. achieve). They make a Taylor series out of Phong shading and use
  3163. only the first two terms of the series for the approximation.
  3164. Thus, once the initialization phase is completed, Phong shading
  3165. requires but two adds per pixel, a little more with specular
  3166. highlights. The main problem, though, is the initialization
  3167. phase, which includes many multiplications, divisions, and even a
  3168. square root. It will doubtlessly be faster than exact Phong
  3169. shading, but I am not so confident that it will be sufficiently
  3170. fast to become very popular. Using Gouraud shading and a very
  3171. large number of polygons can yield comparable results at a very
  3172. acceptable speed.
  3173.  
  3174.  
  3175. Approximating Phong Shading
  3176.  
  3177.  
  3178. There are a few approaches to approximating phong shading. They
  3179. all boil down to a simple technique. Find the spherical
  3180. coordinates equivalent of the normal vector at each vertex and
  3181. interpolate this value linearly along edges and inside the
  3182. polygon. Or, you can interpolate the (x,y) value of the normal
  3183. vector linearly, and then, instead of normalizing the (x,y,z)
  3184. vector, use a lookup table to find what exact z value will make
  3185. the particular (x,y,z) vector normal.
  3186.  
  3187. If you are using a lookup-table technique, the lookup-table can
  3188. be viewed as a texture mapped onto the polygon. In this case, you
  3189. could envision placing something other than z value or shade
  3190. value in the lookup table. Nice effect such as environment
  3191. mapping can be achieved this way.
  3192.  
  3193.  
  3194.  
  3195.  
  3196.  
  3197.      _________________________
  3198.  
  3199.  
  3200.  
  3201.  
  3202. Version history
  3203.  
  3204.  
  3205. Original version (no version number): original draft, chapters 1
  3206. through 8.
  3207.  
  3208. Version 0.10 beta: added chapters 9 through 11, added version
  3209. history. Still no proofreading, spellchecking or whatever.
  3210.  
  3211. Version 0.11 beta: modified some chapters, noticed than (Q,P,R)
  3212. is the normal vector (chapter 5). A little proofreading was done.
  3213.  
  3214. Version 0.12 beta: I just noticed that this document needs a
  3215. revision badly. Still, I don't have time. I just added a few
  3216. things about the scalar multiplication of two vectors being the D
  3217. coefficient in the plane equation. Added chapter 12. Also
  3218. corrected a part about using arrays instead of lists. The
  3219. overhead is just as bad. Better to use lists. Ah, and I have to
  3220. remember to add a part on calculating the plane equation from
  3221. more than three points to reduce roundoff error.
  3222.  
  3223. Version 0.20 beta: Still need to add the part about calculating
  3224. A,B and C in plane eq. from many points. However, corrected
  3225. several of the mistakes I talked about in the preceding paragraph
  3226. (i.e. revision). Ameliorated the demonstration for 2d rotations.
  3227. Started work on 3d rotations (this had me thinking...)
  3228.  
  3229. Version 0.21 beta: converted to another word processor. Will now
  3230. save in RTF instead of WRI.
  3231.  
  3232. Version 0.22 beta: corrected a few goofs. 3d rotations still on
  3233. the workbench. Might add a part about subdividing a polygon in
  3234. triangles (a task which seems to fascinate people though I'm not
  3235. certain why). Will also add the very simple Gouraud shading
  3236. algorithm in the next version (and drop a line about bump
  3237. mapping). This thing needs at least a table of contents, geez...
  3238. I ran the grammar checker on this. Tell me if it left anything
  3239. out.
  3240.  
  3241. Version 0.23 beta: well, I made the version number at the start
  3242. of this doc correct! :-) Did chapter 13 (subdividing into
  3243. triangles) and chapter 14 (Gouraud shading).
  3244.  
  3245. Version 0.24 beta: removed all GIFs from file because of recent
  3246. CompuServe-Unisys thingy. Don't ask me for them, call CompuServe
  3247. or Unisys [sigh].
  3248.  
  3249. Version 0.30 beta: added the very small yet powerful technique
  3250. for texture mapping that eliminates the division per pixel. It
  3251. can also help for z-buffering. Corrected the copyright notice.
  3252.  
  3253. Version 0.31 beta: I, err, added yet another bit to texture
  3254. mapping [did this right after posting it to x2ftp.oulu, silly
  3255. me]. Corrected a typo in chapter 13: a n-sided polygon turns into
  3256. n-2 triangles, not n triangles. Was bored tonight, added
  3257. something to the logs chapter (I know, it's not that useful, but
  3258. hey, why would I remove it?) Added «Canada» to my address below
  3259. (it seems I almost forgot internet is worldwide tee hee hee).
  3260. Addendum: I'm not certain that the GIF thing (see the paragraph
  3261. about v0.24beta above) applies to me, but I'm not taking any
  3262. chances.
  3263.  
  3264. Version 0.32 beta: added bits and pieces here and there. The WWW
  3265. version supposedly just became available today. Special thanks to
  3266. Lasse Rasinen for converting it. Special thanks also to Antti
  3267. Rasinen. I can't seem to be able to connect to the WWW server
  3268. though, I just hope it really works. Corrected a goof in the free
  3269. directional texture mapping section.
  3270.  
  3271. Version 0.40 beta: Kevin Hunter joined forces with me and made
  3272. that part in chapter two about changes in coordinates system,
  3273. finding the square root et al. Modified the thing about
  3274. triangulating a polygon. Added the Sutherland-Hodgman clipping
  3275. algorithm and related material (the whole chapter 15). This
  3276. document deserves a bibliography/suggested readings/whatever, but
  3277. I don't think I have the time for that. Might add a pointer to
  3278. the 3d books list if I can find where it is. Anyway, if you're
  3279. really interested and want to know more, I'm using Computer
  3280. Graphics, Principles and Practice by Foley, van Dam, Feiner and
  3281. Hughes, from Addison-Wesley ISBN 0-201-12110-7. There. It's not
  3282. in the right format, but it's there. MY E-MAIL ADDRESS CHANGED!
  3283. PLEASE USE THIS ONE FROM NOW ON!
  3284.  
  3285. Version 0.41 beta: added to chapter 3. Gave a more detailed
  3286. discussion of a scan-line algorithm. Put in a line about Paul
  3287. Nettle's sbuffering in the same chapter. Will now include a
  3288. "where can I get this doc from" in the readme, and I'll try to
  3289. put the last few version history entries in my announcements for
  3290. future versions of the doc.
  3291.  
  3292. Version 0.42 beta: removed the buggy PS file. Added an ASCII file
  3293. containing most of the pics of the document.
  3294.  
  3295. Version 0.50 beta: added chapters 16, 17 and 18 about Phong
  3296. shading and illumination model, interpolations forward
  3297. differencing and Taylor series. I am also including an extract
  3298. from SIGGRAPH '86 for fast Phong shading. Source code is coming
  3299. out good, so it just might be in the next version.
  3300.  
  3301. Version 0.60 beta: added source code and Phong approximations in
  3302. chapter 18. Zed3D is now SHAREWARE. It means you should give me a
  3303. small contribution if you find it useful. See the beginning of
  3304. the document for more details.
  3305.  
  3306.  
  3307.  
  3308.  
  3309.      _____________________________
  3310.  
  3311.  
  3312.  
  3313.  
  3314.  
  3315.  
  3316. About the author
  3317.  
  3318.  
  3319.  
  3320. Sébastien Loisel is a student in university. He is studying
  3321. computer sciences and mathematics and has been doing all sorts of
  3322. programs, both down-to-earth and hypothetical, theoritical
  3323. models. He's been thinking a lot about computer graphics lately
  3324. and so decided to write this because he wanted to get his ideas
  3325. straight. After a while, he decided to distribute it on Internet
  3326. (why the hell am I using third person?)
  3327.  
  3328. The author would be happy to hear your comments, suggestions,
  3329. critics, ideas or receive a postcard, money, new computer and/or
  3330. hardware or whatever (grin). If you do want to get in touch with
  3331. me, well here's my snail mail address:
  3332.  
  3333. Sébastien Loisel
  3334. 1 J.K. Laflamme
  3335. Lévis, Québec
  3336. Canada
  3337. G6V 3R1
  3338.  
  3339. MY E-MAIL ADDRESS CHANGED! PLEASE USE THIS ONE FROM NOW ON!
  3340. loiselse@ift.ulaval.ca
  3341.  
  3342. See also the enclosed file, README.3D.
  3343.  
  3344. I would like to thank Lasse and Antti Rasinen making a WWW
  3345. version available in Finland.
  3346.  
  3347. I would like to thank Simos Hadjiyiannis for making the WWW
  3348. available to the rest of the world. Mr. Hadjiyiannis can be
  3349. reached at hadjiyia@cs.colostate.edu. The WWW page is
  3350. http://www.cs.colostate.edu/~hadjiyia/3d.html.
  3351.  
  3352. I would like to thank Kevin Hunter for writing a big chunk of
  3353. chapter 2. Thank you Kevin. Kevin Hunter can be reached (at the
  3354. time of this writing) as HUNTER@symbol.com.
  3355.  
  3356.  
  3357.