home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / pc / UGPRG.ZIP / 3D / 3D.DOC next >
Encoding:
Text File  |  1996-07-27  |  121.4 KB  |  3,053 lines
  1.  
  2.                               3DGPL  Version 1.0
  3.                              --------------------
  4.                                GRAPHICS TUTORIAL
  5.                      ------------------------------------
  6.                              FOR GAME PROGRAMMERS
  7.                             ----------------------
  8.  
  9.  
  10.  
  11.            Copyright (c) 1995 Sergei Savchenko (savs@cs.mcgill.ca)
  12.  
  13. ------------------------------------------------------------------------------                
  14.         Original Archive is available @x2ftp.oulu.fi
  15.     Please note that this document does not contain the source!
  16. ------------------------------------------------------------------------------                
  17.  
  18.                                    * * *
  19.  
  20.                              3-D transformations.
  21.                             ----------------------
  22.  
  23.  
  24.  
  25.                                     "Several large, artificial constructions
  26.                                     are approaching us," ZORAC
  27.                                     announced after a short pause.
  28.                                     "The designs are not familiar, but they
  29.                                     are obviously the products of
  30.                                     intelligence. Implications: we have been
  31.                                     intercepted deliberately by a means
  32.                                     unknown, for a purpose unknown, and
  33.                                     transferred to a place unknown by
  34.                                     a form of intelligence unknown.
  35.                                     Apart from the unknowns, everything
  36.                                     is obvious."
  37.                                     -- James P. Hogan, "Giants Star"
  38.  
  39.  
  40.  
  41. 3D transformations would assume the central position in the engine.
  42. After all, this is something which would allow us to see objects from
  43. different sides, manage the three dimensional space etc. The
  44. simplest but by no means least important transformations are
  45. translation and scaling. The rotation transformation would be little
  46. bit more complicated to understand and implement but fortunately
  47. not much. Talking of all three transformations we can either consider
  48. them from the point of an object (collection of points in 3D)
  49. being transformed, or a space itself being transformed. In
  50. some cases it is more productive to think one way in some another way.
  51.  
  52.  
  53.  
  54. Translation and scaling.
  55. ------------------------
  56. Translation. This transformation moves points of an object to a new
  57. specified position pic 2.1, On the other hand we may think of it as a
  58. system of references being moved in the opposite direction pic 2.2:
  59.  
  60.  
  61.  
  62.           Y ^                                                 ^
  63.             |          * A'
  64.             |          ^ * B'                                 |
  65.             |          | ^                           Y'^
  66.             |A*--------+ | y component                 |      |A*
  67.             |  B*--------+                             |         B*
  68.             |    x component                           |      |
  69.    ---------+-------------------> X        - - - - - - | - - - - - - -> X
  70.            0|                                          |     0|
  71.             |                               -----------+------------->X'
  72.             |                                        0'|      |
  73.             |                                          |
  74.             |                                          |
  75.             |                                          |
  76.  
  77.          pic 2.1                                    pic 2.2
  78.  
  79. This way if it is a collection of points we are moving, Distances
  80. between them would not change. Any move across 2D space can separate
  81. into 2 components - move along X and move along Y. It would take
  82. 3 components in 3D space - along X, Y and Z. Important place we
  83. may need translation is to , for example , move (0,0) of the screen
  84. from upper left corner of the screen to it's physical centre. So
  85. the translation function might look like:
  86.  
  87.  
  88. void T_translation(register int *from,register int *into,int number,
  89.                    int addx,int addy,int addz
  90.                   )
  91. {
  92.  register int i;
  93.  
  94.  for(i=0;i<number;i++)
  95.  {
  96.   (*into++)=(*from++)+addx;
  97.   (*into++)=(*from++)+addy;
  98.   (*into++)=(*from++)+addz;                 /* translation */
  99.  }
  100. }
  101.  
  102.  
  103. Another transformation is scaling: blowing up of distances between
  104. points. Or we can think of it as a contraction or blowing of the space
  105. itself pic2.3:
  106.  
  107.  
  108.                  A'*      Y ^
  109.                     \       |
  110.                      \      |
  111.                       *   B'*    *C'
  112.                      A      ^   /
  113.                            B*  *C
  114.                             |
  115.               --------------+-----------------> X
  116.                            0|
  117.                             |
  118.                             |
  119.  
  120.                          pic 2.3
  121.  
  122. The obvious usage of this transformation: trying to fit 2000x2000 object
  123. into 200x200 window - evidently every point of this object would have to
  124. have coordinates scaled to fit in the window range. In this case scaling
  125. is just multiplication by 1/10. The object is contracted that way. On the
  126. other hand we can think that it was actually the window space expanding to
  127. enclose the object (we would have to abstract far enough from computer
  128. hardware today on the market to do that, however). In general case
  129. we can introduce separate factors for every axis, if this factor is
  130. greater then 0 and less then 1 we would contract the object if it would
  131. be greater then 1 the object would be expanded. The code is once again
  132. trivial:
  133.  
  134.  
  135. void T_scaling(register int *from,register int *to,int number,
  136.                float sx,float sy,float sz
  137.               )
  138. {
  139.  register int i;
  140.  
  141.  for(i=0;i<number;i++)
  142.  {
  143.   (*to++)=(int)(*from++)*sx;
  144.   (*to++)=(int)(*from++)*sy
  145.   (*to++)=(int)(*from++)*sz;
  146.  }                                          /* scaling */
  147. }
  148.  
  149.  
  150.  
  151. Planar case rotation.
  152. ---------------------
  153. Of those three transformations rotation is by far the most complicated
  154. one, let's first consider planar - 2-D case of rotation. there is some
  155. simple trigonometry involved, sin and cos functions. just for those
  156. who don't deal with this subject every day sin and cos for a triangle
  157. on pic 2.4 are defined as:
  158.  
  159.  
  160.                                   +
  161.                                  /|
  162.                                l/ |
  163.                                /  |y
  164.                               /alp|
  165.                              +----+
  166.                                 x
  167.  
  168.                              pic 2.4
  169.  
  170.  
  171.                             x                  y
  172.                 sin(alp) = ---     cos(alp) = ---
  173.                             l                  l
  174.  
  175. And so if we do know sin(alp) or cos(alp) we can calculate segments x and y
  176. knowing segment l:
  177.  
  178.                  x = l*sin(alp)    y = l*cos(alp)
  179.  
  180. This can be considered as seeking projections of a point for which we
  181. know the distance from the beginning of coordinates onto orthogonal
  182. axes pic 2.5:
  183.  
  184.                                ^
  185.                              Y |
  186.                                +----*
  187.                                |   /|
  188.                               y|  /l|
  189.                                | /  |
  190.                                |/alp|
  191.                          ------+----+---->X
  192.                                |  x
  193.                                |
  194.  
  195.                              pic 2.5
  196.  
  197. Now, let's consider the situation when we rotate the system of references
  198. counterclockwise by the angle alp. What would be coordinates of point
  199. A in the turned system of references Y'X' if it's coordinates in the
  200. original system of references were (x,y)?
  201.  
  202.  
  203.                                ^ Y            X'
  204.                                |             /
  205.                     Y'       Y/*\--*A      /
  206.                       \     /  |  \|     /
  207.                         \ /    |   |\  /
  208.                        Yy'\    |   | /Yx'
  209.                             \  |   |
  210.                               \| /\ Xx'
  211.                   -------------+---+-------------------> X
  212.                               /| \/ X
  213.                             / Xy'  \
  214.                           /    |     \
  215.                         /      |       \
  216.                       /        |         \
  217.                                |
  218.  
  219.                              pic 2.6
  220.  
  221.  
  222. Using the sin/cos formulas from above we can find projections of
  223. x and y (that is of the projections of the point onto original axis's)
  224. onto new axis's X'and Y'. Let's sum projection of x onto X' and projection
  225. of y onto same axis's and do the same thing with projections of x and y
  226. onto Y'. Those sums would be just what we are looking for, coordinates 
  227. of the point A in the new system of references X'Y'
  228.  
  229. i.e.
  230.  
  231.                X' = Yx'+Xx' = y*sin(alp) +    x*cos(alp)
  232.                Y' = Yy'+Xy' = y*cos(alp) +  (-x*sin(alp))
  233.  
  234. Note the sign of Xy', just from looking on the pic2.6 one can see that
  235. in this case x is being projected onto the negative side of Y'. That's
  236. why the negative sign. So those are the transformation formulas for
  237. counterclockwise rotating system of reference by an angle alp:
  238.  
  239.                      x'= y*sin(alp) + x*cos(alp)
  240.                      y'= y*cos(alp) - x*sin(alp)
  241.  
  242. On the other hand we can consider it as point A itself rotating
  243. around zero clockwise.
  244.  
  245.  
  246.  
  247. 3-D rotation.
  248. -------------
  249. This subject is complicated only because it is a bit hard to visualize
  250. what is happening. The idea on the other hand is simple: applying three
  251. 2-D rotations one after another so that each next works with the results
  252. obtained on the previous stage. This way every time we are applying new
  253. transformation we are working not with the original object but the
  254. object already transformed by all previous rotations.
  255.  
  256. Each of 2D rotations this way is bound to some axis around which two
  257. other axes would rotate.
  258.  
  259. There are few important considerations influencing the way we would
  260. find the formulas: 1) What kind of system of references we have. What
  261. is the positive direction of the rotations we would be applying; and
  262. 2) In what order we want to apply 2D rotations.
  263.  
  264.  
  265.  
  266. System of references:
  267. ---------------------
  268. Well, we do have freedom of choosing directions of axes, for one
  269. thing, we are working in orthogonal system of references where each
  270. axes is orthogonal to the plane composed of remaining two. As
  271. to where the positive direction of every axis point we have freedom
  272. to decide. It is customary for example to have positive side of Y
  273. directed upside. We don't actually have to follow customs if we really
  274. don't want to, and remembering that the colourmap, we would be doing all
  275. drawings into, have Y directed down we might want to choose
  276. corresponding direction here. However having it pointing up is more
  277. natural for us humans, so we might save time on debugging a while
  278. afterwards. Let's point it up. As to the Z let's direct it away from
  279. us and X to the right.
  280.  
  281.  
  282.                                  Y^    Z
  283.                                   |   /
  284.                                   |  / alp
  285.                                  /|<---+
  286.                                 | |/   |
  287.                          -------|-+-----------> X
  288.                             bet | |   __
  289.                                 V/|   /| gam
  290.                                 /----+
  291.                                /  |
  292.                               /   |
  293.                                   |
  294.  
  295.                                pic 2.7
  296.  
  297.  
  298. As to the rotations, let's call the angle to turn XY plane around Z axis
  299. as gam, ZY around X as bet and ZX around Y as gam, pic 2.7
  300.  
  301.  
  302.  
  303. Transformation order.
  304. ---------------------
  305. The problem with transformation order is that the point consequently
  306. turned by gam-bet-alp won't be necessarily at the same place with
  307. where it might go if turned alp-bet-gam. The reason is our original
  308. assumption: each next transformation works on already transformed
  309. point by previous 2D rotations. That is, we are rotating coordinates
  310. around moving axes. On the picture 2.7 we can see that rotation
  311. bet is performed around axis x but the next rotation alp is performed
  312. around axis z which after application of previous transformation
  313. would move.
  314.  
  315.  
  316.                              +              +----+
  317.                     z /     /|             /    / alp
  318.                     /      + V bet            \
  319.        ------------+-------|---- x     ---------+------------ x
  320.                  /         |                      \
  321.            +----+                                  z
  322.            | /  |
  323.            |    V alp
  324.  
  325.                                 pic 2.8
  326.  
  327.  
  328. The implication is: we have to know with respect to what each
  329. next rotation is performed, that is what is the order of applications
  330. of 2D rotations. Let's think how we are coordinating the surrounding
  331. world? First we think of the direction. Then we can tilt our head up
  332. and down, and finally from there we can shake it from left to right.
  333. Now, when we are tilting head we have already chosen the direction.
  334. If we first would tilt head then the directional rotation to be
  335. performed after would not be parallel to the ground but rather
  336. perpendicular to the imaginary line being continuation of our
  337. tilted neck. ( No responsibility hereby assumed for all direct
  338. or consequential damage or injury resulted from doing that ) So
  339. it all comes to directional rotation first - gam in our case.
  340. pitch second - bet;  roll last - alp. Let's do just that using
  341. obtained formula for 2D rotation:
  342.  
  343.  
  344.   ^Z                     ^Y'                    ^Y"
  345.   |                      |                      |
  346.   |<---+                 |<---+                 |<---+
  347.   | gam|                 | bet|                 | alp|
  348.   |    |                 |    |                 |    |
  349.  -+----+--->X           -+----+--->Z'          -+----+--->X"
  350.   |                      |                      |
  351.  
  352.  
  353.  x'=z*sin(gam)+x*cos(gam)
  354.  y'=y
  355.  z'=z*cos(gam)-x*sin(gam)
  356.  
  357.                       x"=x
  358.                       y"=y*cos(bet)-z*sin(bet)
  359.                       z"=y*sin(bet)+z*cos(bet)
  360.  
  361.                                               x"'=y*sin(alp)+x*cos(alp)
  362.                                               y"'=y*cos(alp)-x*sin(alp)
  363.                                               z"'=z
  364.  
  365.                            pic 2.9
  366.  
  367. That's basically it, using those formulas we can apply 3D rotation
  368. to coordinates x,y,z and get x"',y"',z"'. We can see here 12
  369. multiplication's, The question is can we reduce this number? Let's
  370. try getting rid of temporary variables x',y',z',x",y",z". We can do it
  371. by just substituting each occurrence of them by their real meaning:
  372.  
  373. First let's try obtaining x",y" and z" expressed via x,y,z:
  374.  
  375.  
  376. y"=y'*cos(bet) - z'*sin(bet)
  377. y"=y*cos(bet) - (z*cos(gam)-x*sin(gam)) *sin(bet)
  378. y"=y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet)
  379. - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  380. x"=x'
  381. x"= z*sin(gam) + x*cos(gam)
  382. - - - - - - - - - - - - - - -
  383. z"=y'*sin(bet) + z'*cos(bet)
  384. z"=y*sin(bet) + (z*cos(gam)-x*sin(gam)) *cos(bet)
  385. z"=y*sin(bet) + z*cos(gam)*cos(bet) - x*sin(gam)*cos(bet)
  386. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  387.  
  388.  
  389. Now using that we can express x"',y"' and z"' via x,y,z:
  390.  
  391.  
  392. y"'=y"*cos(alp) - x"*sin(alp)
  393.  
  394. y"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*cos(alp) -
  395.     (z*sin(gam) + x*cos(gam)) *sin(alp))
  396.  
  397. y"'=y*cos(bet)*cos(alp)-z*cos(gam)*sin(bet)*cos(alp)+x*sin(gam)*sin(bat)*cos(alp)
  398.     -z*sin(gam)*sin(alp) - x*cos(gam)*sin(alp)
  399.  
  400. y"'=x* [ sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)] +
  401.     y* [ cos(bet)*cos(alp) ] +
  402.     z* [-cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)]
  403. -----------------------------------------------------------
  404.  
  405. x"'= y"*sin(alp) + x"*cos(alp)
  406.  
  407. x"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*sin(alp) +
  408.     (z*sin(gam) + x*cos(gam)) *cos(alp)
  409.  
  410. x"'=y*cos(bet)*sin(alp)-z*cos(gam)*sin(bet)*sin(alp)+x*sin(gam)*sin(bet)*sin(alp)
  411.     z*sin(gam)*cos(alp) + x*cos(gam)*cos(alp)
  412.  
  413. x"'=x* [sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)] +
  414.     y* [cos(bet)*sin(alp) ] +
  415.     z* [sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp) ]
  416. ----------------------------------------------------------
  417.  
  418. z"'=z"
  419. z"'=x* [- sin(gam)*cos(bet)] +
  420.     y* [sin(bet)] +
  421.     z* [cos(gam)*cos(bet)]
  422. ----------------------------------------------------------------
  423.  
  424.  
  425. This really appear to be more complicated that what we had before.
  426. It appear to have much more multiplications than we had. But if we
  427. look on those resulting formulas we can see that all coefficients
  428. in square brackets can be calculated just once so that transformation
  429. of a point would look like:
  430.  
  431.  
  432.                    x"'=x*mx1 + y*my1 + z*mz1
  433.                    y"'=x*mx2 + y*my2 + z*mz2
  434.                    z"'=x*mx3 + y*my3 + z*mz3
  435.  
  436. It can be seen that it takes just 9 multiplication. Of course calculating
  437. all the coefficients takes 16 multiplications:
  438.  
  439. mx1= sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)
  440. my1= cos(bet)*sin(alp)
  441. mz1= sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp)
  442.  
  443. mx2= sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)
  444. my2= cos(bet)*cos(alp)
  445. mz2= -cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)
  446.  
  447. mx3= -sin(gam)*cos(bet)
  448. my3= sin(bet)
  449. mz3= cos(gam)*cos(bet)
  450.  
  451. But if we have 100 points to turn. Original method would take
  452. 12*100=1200 multiplications. This one would take 9*100+16=916 since
  453. coefficients are calculated just once for all points to transform. One
  454. more consideration having to do with optimization, in most cases it makes
  455. sense to measure angles as integers. We don't usually care of fractions
  456. smaller then 1 degree. So we don't really need sin and cos for all the
  457. real number range. That's why we can calculate sin/cos just once for
  458. all 360 degrees before any rotations are performed and put the obtained
  459. values into arrays so that next time we would need any of them we
  460. wouldn't have to call a functions which probably uses a power series
  461. with few multiplications to calculate each. (Processors nowadays come
  462. with floating point units which can calculate sin/cos pretty fast
  463. but unlikely faster then single array access (this might become false
  464. though). By the way we don't really need 360 degrees. This number is
  465. not very convenient. We can go with full angle divided into 256 pseudo
  466. degrees. The reason we would need just one unsigned char to store the
  467. angle and whenever we go beyond 255 we simple wrap around to zero,
  468. saving this way a conditional statement we would need otherwise in
  469. the case of 360 degrees.
  470.  
  471.  
  472. #define T_RADS                  40.7436642  /* pseudo grads into rads */
  473.  
  474. float T_mx1,T_mx2,T_mx3,T_my1,T_my2,T_my3,T_mz1,T_mz2,T_mz3;
  475. float T_sin[256],T_cos[256];                /* precalculated */
  476.  
  477. void T_init_math(void)
  478. {
  479.  int i;
  480.  
  481.  for(i=0;i<256;i++)
  482.  {
  483.   T_sin[i]=sin(i/T_RADS);
  484.   T_cos[i]=cos(i/T_RADS);
  485.  }
  486. }
  487.  
  488.  
  489. Now finally to the code to perform 3D rotation. First function is building
  490. set of rotation coefficients. Note that T_mx1, T_mx2 etc. are global 
  491. variables. The other function: T_rotation uses those coefficients, so 
  492. they better be ready when later is called.
  493.  
  494.  
  495. void T_set_rotation_gam_bet_alp(int alp,int bet,int gam)
  496. {
  497.  T_cosalp=T_cos[alp];
  498.  T_sinalp=T_sin[alp];
  499.  T_cosbet=T_cos[bet];
  500.  T_sinbet=T_sin[bet];
  501.  T_cosgam=T_cos[gam];
  502.  T_singam=T_sin[gam];                       /* initializing */
  503.  
  504.  T_mx1=T_singam*T_sinbet*T_sinalp + T_cosgam*T_cosalp;
  505.  T_my1=T_cosbet*T_sinalp;
  506.  T_mz1=T_singam*T_cosalp - T_cosgam*T_sinbet*T_sinalp;
  507.  
  508.  T_mx2=T_singam*T_sinbet*T_cosalp - T_cosgam*T_sinalp;
  509.  T_my2=T_cosbet*T_cosalp;
  510.  T_mz2=-T_cosgam*T_sinbet*T_cosalp - T_singam*T_sinalp;
  511.  
  512.  T_mx3=-T_singam*T_cosbet;
  513.  T_my3=T_sinbet;
  514.  T_mz3=T_cosgam*T_cosbet;                   /* calculating the matrix */
  515. }
  516.  
  517.  
  518. void T_rotation(int *from,register int *to,int number)
  519. {
  520.  register int i;
  521.  register int xt,yt,zt;
  522.  
  523.  for(i=0;i<number;i++)
  524.  {
  525.   xt=*from++;
  526.   yt=*from++;
  527.   zt=*from++;
  528.  
  529.   *to++=(int)(T_mx1*xt+T_my1*yt+T_mz1*zt);
  530.   *to++=(int)(T_mx2*xt+T_my2*yt+T_mz2*zt);
  531.   *to++=(int)(T_mx3*xt+T_my3*yt+T_mz3*zt);  /* performing the operation */
  532.  }
  533. }
  534.  
  535.  
  536.  
  537. Representing transformations in matrix form.
  538. --------------------------------------------
  539. Matrix is a table of values:
  540.  
  541.       | a b c |
  542.     A=| d e f |
  543.       | g i j |
  544.  
  545. Special cases are column and row vectors:
  546.  
  547.                     | x |
  548.     V=| a b c |   C=| y |
  549.                     | z |
  550.  
  551.  
  552. Row can by multiplied by a column producing a value, a scalar:
  553.  
  554.               | x |
  555.     | a b c |*| y | = a*x + b*y + c*z
  556.               | z |
  557.  
  558. Matrix can be multiplied by another matrix: each value of the 
  559. result was obtained by multiplying respective rows and columns 
  560. of the arguments (since columns and vectors are to have the same 
  561. dimension in order for multiplication on them to be defined it 
  562. automatically places restrictions on dimensions of matrices that 
  563. can be multiplied):
  564.  
  565.                                  | x |            | k |
  566.                      | | a b c |*| y |  | a b c |*| l | |
  567.                      |           | z |            | m | |
  568.                      |                                  |
  569.  | a b c | | x k |   |           | x |            | k | |
  570.  | d e f |*| y l | = | | d e f |*| y |  | d e f |*| l | |
  571.  | g i j | | z m |   |           | z |            | m | |
  572.                      |                                  |
  573.                      |           | x |            | k | |
  574.                      | | g i j |*| y |  | g i j |*| l | |
  575.                                  | z |            | m |
  576.  
  577. We may be interested in multiplying a row vector by a 3x3 matrix:
  578.  
  579.            | a b c |
  580.  | x y z |*| d e f | = | x*a+y*d+z*g  x*b+y*e+z*i  x*c+y*f+z*j |
  581.            | g i j |
  582.  
  583. Looks familiar? Yes, the rotation coefficients we've calculated
  584. lie nicely in the 3x3 matrix and so we can think of 3D rotation
  585. as of matrix multiplication of "rotation" matrix onto original
  586. vector producing result vector. We can also try fitting translation 
  587. and scaling into a matrix, to do that for translation, however, we 
  588. would have to increase dimensions of the coefficient matrices and 
  589. the argument vector:
  590.  
  591. Translation:
  592.  
  593.              | 1  0  0  0 |
  594.  | x y z 1 |*| 0  1  0  0 | = | x+tx  y+ty  z+tz  1 |
  595.              | 0  0  1  0 |
  596.              | tx ty tz 1 |
  597.  
  598. Scaling:
  599.  
  600.              | sx  0  0 0 |
  601.  | x y z 1 |*|  0 sy  0 0 | = | x*sx  y*sy  z*sz  1 |
  602.              |  0  0 sz 0 |
  603.              |  0  0  0 1 |
  604.  
  605. The advantage of the matrix approach is that if we have several
  606. transformations each represented in matrix form:
  607.  
  608.             | X'| = (| X |*| A |) *| B |
  609.  
  610. Where A, and B are transformation matrixes and | X | the argument
  611. vector to be transformed we can:
  612.  
  613.              | X'| = | X |* ( | A |*| B | )
  614.  
  615. calculate the "concatenated" matrix | K | = | A |*| B | and each
  616. following coordinate just transform as:
  617.  
  618.                     | X'| = | X |*| K |
  619.  
  620. By the way by representing each of 2-D rotations in 3x3 matrix form we
  621. can get the formulas we have found for 3-D rotation by calculating a
  622. concatenated matrix.
  623.  
  624. So, matrix approaches might be very effective in instances where
  625. there are a lot of consecutive transformation. For example if we
  626. have consecutive rotations it makes sense to calculate their
  627. concatenated matrix and use it from then on. It would save
  628. us quite a few multiplication, one verses several consecutive matrix
  629. multiplication for each point to transform.
  630.  
  631. And in the run-time it is preferably to use pre-calculated
  632. formulas for concatenated matrix coefficients then trying to get it
  633. every time by multiplying respective transformational matrices.
  634. The reason 3 consecutive matrix multiplication of, for instance, 3x3
  635. matrices take 9+9+9=27 multiplication's. But since lot of values in
  636. those matrices are 0, it is taken into account in pre-calculated
  637. coefficient formulas which allow to produce the very same results
  638. with just 16 multiplications.
  639.  
  640. As to the translation transformation matrix and especially increasing
  641. the dimension from 3 to 4 it looks much less natural. 4x4 matrix means
  642. more multiplications too, of course some of the value in 4x4 would
  643. always be zero so we can optimize (cripple) the routine to original 9
  644. multiplications and 3 additions. In one word, there's a freedom
  645. to decide what approach to take.
  646.  
  647.  
  648.  
  649. perspective.
  650. ------------
  651. I believe, it is not hard to visualize what perspective is, perhaps the
  652. first association lot of us would have, is this empty straight street
  653. with identical buildings on both sides disappearing in infinity. ( At
  654. least people residing in newer neighborhoods of Moscow or Toronto know
  655. what I am talking about ) Why do we need this transformation? Just to
  656. achieve landscape realism, after all this same street would look pretty
  657. weird if the road would not collapse into a dot and buildings further
  658. away from us won't look smaller. So weird it might be even hard to
  659. picture.
  660.  
  661. Perspective might not be needed for CAD style applications, but even
  662. there it would, most likely, be implemented in one viewing mode or
  663. another.
  664.  
  665. Before we can implement this transformation it is a good idea to understand
  666. physics of what is happening. A ray from some point in space is being
  667. caught by the viewer's eye. This ray had intersected a plane located
  668. in front of the viewer. If we can find an intersection point and plot the
  669. dot there, for the viewer it would be the same as if the ray from the
  670. plotted dot was actually coming from the initial point in space.
  671.  
  672.  
  673.                              * A
  674.                        |   /
  675.                      A'| /
  676.                        *
  677.                      / | B'
  678.                    *---*-----------------------* B
  679.                      \ |
  680.                        *\
  681.                      C'|   \
  682.                        |      \
  683.                                 * C
  684.  
  685.                               pic 2.10
  686.  
  687.  
  688. Rays from different points in space , in this model, all converge
  689. in one point - viewer's eye. The math here is quite straightforward 
  690. pic 2.11:
  691.  
  692.  
  693.                   screen plane->|   *
  694.                                 |  /| A (X,Z)
  695.                                 | / |
  696.                                 |/  |
  697.                             ^   *   |
  698.                             |  /|X' |
  699.                               / |   |
  700.                             |/  |   |
  701.                             *---+---+- - - - ->
  702.                             0   |   Z
  703.                             |   |
  704.            focus distance ->|---|<-
  705.  
  706.                              pic 2.11
  707.  
  708.  
  709. Let's consider 2-D case first: the viewer's eye is located at point 0
  710. the distance between viewer's eye and the screen would be called -
  711. "focus". The problem is to determine at what point the ray going from
  712. point A into the viewer's eye would intersect the screen. We would have
  713. to plot the dot just there. This is very simple, just yet another
  714. case of solving similar triangles. Just from the fact that the angle
  715. at 0 is the same for both triangles, and if two angles are the same their
  716. tangents would be too. We have:
  717.  
  718.            X'      X                          X*focus
  719.          ----- = -----         =>        X'= -----------
  720.          focus     Z                             Z
  721.  
  722. Same considerations apply for Y dimension, together describing 3-D case:
  723.  
  724.           Y'=Y*focus/Z                   X'=X*focus/Z
  725.  
  726. It is a bit of a mystery what is happening with Z coordinate. Basically
  727. after the perspective transformation we would be performing rendering
  728. onto the screen which is 2-D surface we would need X and Y but hardly
  729. Z, so we might as well have Z'=0. However when trying to render multiface
  730. objects we would need to know the depth (Z) of all polygons so that
  731. we can decide which are visible and which are obscured. So we might
  732. have Z'=Z that is just leave old value. Then again by doing transformation
  733. on X and Y and leaving Z unchanged we actually may allow for depth
  734. relation to change. That is a polygon obscuring another one in the original
  735. space, before projection, can actually come out behind or partly behind
  736. with this kind of transformation pic 2.21.
  737.  
  738.             ^ Z
  739.   same      |                   *                     *
  740.  deapth     |               * /                 * /
  741.  before     |           A / /                //
  742.   and       |           * /               /* A
  743.  after      |           /B            / B
  744.             |          *            *
  745.             |        before                 after
  746.  
  747.                               pic 2.12
  748.  
  749. To avoid that we can apply some non linear transformation on Z,
  750. like Z'=C1+C2/Z, where C1, C2 are some constants.
  751.  
  752.  
  753. In the enclosed code only X and Y are being transformed. The implementation
  754. of all those things is that easy I won't even put a code example here.
  755. However, there are few possible questions: first: what should be the value
  756. for "focus distance" ?
  757.  
  758.  
  759.          \                    /          \              /
  760.            \                /             \            /
  761.              \            /                \          /
  762.            ---I\--------/I---           ---I\--------/I---
  763.                  \    /                      \      /
  764.                    \/                         \    /
  765.                                                \  /
  766.                                                 \/
  767.  
  768.                                pic 2.13
  769.  
  770.  
  771. pic 2.13 shows it's physical sense, focus basically determines the view
  772. angle. when it is smaller - angle is wider, when it is bigger angle
  773. is smaller. It is helpful to think of this distance as being measured
  774. in screen pixels, this way for 320 pixel display if focus is chosen
  775. to be 180 the view angle we are getting is 90'. Perspective
  776. transformation might produce at first, weird looking results, with
  777. distortions sometimes similar to that of wide-range photographic
  778. pictures. One has to experiment a bit to find what's looking best,
  779. values giving view angle somewhere between 75-85' would probably
  780. look not that terrible, but that depends of scene and screen geometries.
  781.  
  782.  
  783.  
  784. clipping problem...
  785. -------------------
  786. Where transforms point with Z equal to 0? since we are dividing by
  787. Z that would be infinity, or at least something producing divide error,
  788. in this computer-down-to-earth-world. The only way to deal with this
  789. is to guarantee that there would be no coordinates with such Zs. another
  790. thing, calculations for points with negative Z would produce negative
  791. coordinates, what we would see is objects (or even parts thereof) flipped
  792. over, but then again objects with negative coordinates are behind the
  793. viewing plane and effectively behind the viewer, so this is something
  794. we can't see and better make sure this is not being rendered. The way
  795. to do it is by applying 3-D clipping.
  796.  
  797. The perspective transformation  can also be represented as a matrix,
  798. however, let's think of all transformations discussed up to this moment
  799. a point would have to go through in order to be rendered. First we would
  800. rotate the world with respect to viewer orientation, then before we can
  801. apply perspective transformation we would have to do clipping and
  802. at least get rid of the vertices behind the viewing plane, and only then
  803. can we apply perspective transformation. Representing clipping in
  804. matrix form is not the most pleasant thing there is. And we've already 
  805. discussed that matrices would be of use when we have several consecutive 
  806. transformations so that all of them would be represented by one matrix. 
  807. In this case however rotation and perspective are separated by the clipping. 
  808. If we are sure that none of the coordinates would ever get behind the
  809. viewing plane, we don't actually need clipping, that's the instance to
  810. represent everything in the matrix form, otherwise it might not be
  811. a good idea neither in terms of code complicity nor performance.
  812.  
  813. fixed point math.
  814. -----------------
  815. In previous chapters in several places we didn't have another choice but to
  816. use floating point multiplications. (sin and cos are after all real-valued
  817. functions) However this is still fairly expensive operation, integer
  818. multiplications and divisions would usually be faster. However since we
  819. have to represent numbers with fractions it appeared we didn't have another
  820. choice but to use floating point operations. The question is, can we
  821. represent fractions using just integers? One might think there's a
  822. contradiction in the question itself, however we would see that, as in
  823. many other cases, we can use one mechanism to represent something else.
  824. This is the case here, We may use fixed point numbers representing them
  825. as integers, and with small amendments we can use regular integer additions
  826. multiplications etc, to perform fixed point operations. First of all how do 
  827. we represent integers?
  828.  
  829.  
  830.                    +------+------+------+------+------+
  831.                    |   1  |   0  |  0   |  0   |   1  |
  832.                    +------+------+------+------+------+
  833.                      4       3       2      1      0
  834.                     2 =16   2 =8    2 =4   2 =2   2 =1
  835.  
  836.                                 pic 2.14
  837.  
  838. With any counting system, digits in multi digit numbers would have
  839. different weight, so to speak. With decimal numbers this weight
  840. would be some power of ten, that is 102 is actually:
  841.  
  842.        10**2*1 + 10**1*0 + 10**0*2 == 100+0+2 ==102
  843.  
  844. Same with binary numbers, the only difference, weight would be
  845. some power of 2. this way number on pic 2.13 would be:
  846.  
  847.  2**4*1 + 2**3*0 + 2**2*0 + 2**1*0 + 2**0*1 == 16+0+0+0+1 == 17
  848.  
  849. Now, it is quite natural for us to place a decimal point into a decimal
  850. number, placing a binary point into a binary number should not
  851. seem more difficult. How do we interpret the decimal point? It is
  852. just that for the numbers to the right of the decimal point we
  853. pick negative power of ten. That is:
  854.  
  855.   102.15 == 10**2*1 + 10**1*0 + 10**0*2 + 10**-1*1 + 10**-2*5 ==
  856.                     1      5             15
  857.   == 100 + 0 + 2 + ---- + ---- == 102 + -----
  858.                     10    100            100
  859.  
  860. In this representation the point separates negative powers from
  861. zero power and unlike numbers represented in exponential form - the
  862. floating point numbers, in this case the point never changes it's
  863. position.
  864.  
  865. Nothing should stop us from extending the same scheme onto binary
  866. numbers:
  867.  
  868.                    +------+------+------+------+------+
  869.                    |  1   |   0  |   0  |   1  |   1  |
  870.                    +------+------+------+------+------+
  871.                      2       1     0     -1     -2
  872.                     2 =4    2 =2  2 =1   2 =1/2 2 =1/4
  873.  
  874.                                 pic 2.15
  875.  
  876. Using the very same technique we can see that number represented
  877. on the pic 2.15 can be considered also as:
  878.  
  879.   2**2*1 + 2**1*0 + 2**0*0 + 2**-1*1 + 2**-2*1 == 4+0+0+1/2+1/4
  880.  
  881.   == 4 + 3/4
  882.  
  883. Numbers in what range can we represent? Resorting once again to
  884. analogy with decimal numbers it can be seen that two digits
  885. to the right from the decimal point can cover the range of 0.99
  886. in 0.01 steps. Numbers smaller then minimal 0.01 precision step
  887. just can't be represented without increasing number of digits
  888. after the decimal point. Very same thing is happening with binaries
  889. in the example at pic 2.15 since we have two binary digits- minimal
  890. number we can represent is 1/4 ( 0.01 (bin) == 1/4 (dec) ) and
  891. again numbers with higher precision can be represented only by
  892. increasing number of binary digits after the binary point.
  893.  
  894. Now to the representation of negative numbers, there actually
  895. several ways to do that but effectively today one has prevailed,
  896. this is so called 2's compliment form. The idea is to negate
  897. the weight of the leftmost digit pic 2.16:
  898.  
  899.                    +------+------+------+------+------+
  900.                    |   1  |   1  |   1  |   1  |   1  |
  901.                    +------+------+------+------+------+
  902.                      4       3       2      1      0
  903.                    -2 =-16  2 =8    2 =4   2 =2   2 =1
  904.  
  905.                                 pic 2.16
  906.  
  907. the number represented above is considered to be:
  908.  
  909.   -2**4*1 + 2**3*1 + 2**2*1 + 2**1*1 + 2**0*1== -16+8+4+2+1 = -1
  910.  
  911. The advantage of this approach, we don't have to change addition
  912. and subtraction algorithm working for positive integers in order
  913. to accommodate 2's compliment numbers. Why is that? just by the
  914. virtue of natural overwrap of integer numbers the way they are
  915. represented in the computers. If we add 1 to the maximum number
  916. we can represent we would obtain a carry from the most significant
  917. (leftmost) bit, which is ignored, and the value of exactly 0.
  918. -1 in signed representation is exactly the maximum possible value
  919. to represent in unsigned representation, -1+1==0 indeed. (One can
  920. check this to be true for all numbers and by more formal means, but
  921. since this is not the point I would ignore it) Again, addition and
  922. subtraction don't have to be changed to accommodate 2's compliment
  923. negative numbers (multiplication however should, that's why there
  924. are instructions for both signed and unsigned multiplications
  925. in most computers)
  926.  
  927. Since the leftmost digit in 2's compliment representation carries
  928. negative weight and that's exactly the one with highest power the
  929. minimum negative number possible to represent in the example
  930. above would be 10000 (bin) = -16 all the other digits have
  931. positive weights so maximum possible positive number would be
  932. 01111 (bin) = 15 but this slight asymmetry is not a problem in
  933. the majority of cases.
  934.  
  935. However, values of sin and cosin functions are between 1 and -1 so
  936. to represent them we might choose format with just one integer
  937. field:
  938.  
  939.                    +------+------+------+------+------+
  940.                    |   0  |   1  |   1  |   1  |   1  |
  941.                    +------+------+------+------+------+
  942.                      0      -1      -2     -3     -4
  943.                    -2 =-1   2 =1/2  2 =1/4 2 =1/8 2 =1/16
  944.  
  945.                                  pic 2.17
  946.  
  947. but, because of the asymmetry described above there would be a
  948. value for -1 (1000) but there would be no pure 1, just it's
  949. approximation (01111): (pic 2.17) 1/2+1/4+1/8+1/16, Then again, for
  950. most graphics applications when we have, say, 15 fields representing
  951. fractions this would be close enough and won't cause a lot of problems.
  952. As you can see, the very same rules work for 2's compliment whether
  953. with fixed point or without.
  954.  
  955.  
  956.  
  957. multiplication of fixed point numbers.
  958. --------------------------------------
  959. Let's consider what is happening when we are multiplying two
  960. decimal numbers, for example we have to multiply an integer
  961. by a number with a decimal point and we would need just
  962. an integer as a result:
  963.  
  964.                             1.52
  965.                         *     11
  966.                            -----
  967.                             1 52
  968.                         +  15 2
  969.                            -----
  970.                            16.72
  971.  
  972. Actual result of multiplications has as many digits after
  973. the decimal point as there were in both arguments. Since
  974. we need just an integer as a result we would discard numbers
  975. after the point effectively shift digits to the right by
  976. 2 in this case. We can do the same with fixed point binary
  977. numbers, it means that we have to readjust precision after
  978. the multiplication. Say, if we are multiplying an integer
  979. by a number having 8bit considered to be after the binary
  980. point the result would also 8bit after the point. If it
  981. is just the integer part we are interested in, we have to
  982. shift the result right 8 bit destroying all fractional
  983. information. Similar things are happening with division
  984. except that if we want to divide two integers and obtain
  985. a fixed point result the argument to be divided has to
  986. be shifted left- added fixed point fields, so that effectively
  987. we are dividing a fixed point number by an integer.
  988.  
  989. How to implement all that? Well, first of all we have to
  990. decide what kind of values and what operations on them
  991. we would need. Sometimes it is nice just to have general
  992. purpose fixed point math library, on the other hand we
  993. may choose several different formats of fixed point numbers
  994. one for every special place. Another thing we have to
  995. consider is: would we be able to represent all the operations
  996. using just high level C constructs or it would be necessary
  997. to descend into assembly instructions. The reason we might
  998. consider later is: in most assemblies the result of
  999. multiplication has twice as many bits as each of arguments,
  1000. and since we would occupy some bits with fractions we really
  1001. need all bits we can get. Another thing, the result of
  1002. multiplication would often be located in two registers, So
  1003. we can organize operations the way so that instead of adjusting
  1004. result by shifting we would just take it from the register
  1005. carrying higher bits, -effectively doing zero-cost shifting
  1006. right by the bit length of the register. On the other hand,
  1007. do we really want to compromise code portability by coding
  1008. in assembly? Obviously, there are valid reasons to answer that
  1009. both ways. In particular case of this text and associated
  1010. code we are interested in good portability so let's try to
  1011. stay within straight C boundaries.
  1012.  
  1013. Where would we be using fixed point? first of all - in 3-D
  1014. transformations, to keep rotation matrix coefficients. We
  1015. would need two kinds of multiplication: In calculation of
  1016. the coefficients: fixed * fixed producing same precision
  1017. fixed and in actual coordinate transformation fixed * int
  1018. producing int. lets define fixed to be:
  1019.  
  1020.  
  1021.  
  1022. typedef signed_32_bit fixed;                /* better be 32bit machine */
  1023. #define T_P                             10  /* fixed math precision */
  1024.  
  1025.  
  1026.  
  1027. the T_P would be the number of bits we are allocating to
  1028. keep fractions. Let's make a little function to transform the
  1029. floating point values in the range [-1,1] into fixed point
  1030. values:
  1031.  
  1032.  
  1033.  
  1034. fixed TI_float2fixed(float a)
  1035. {
  1036.  fixed res=0;
  1037.  int i,dv,sign;
  1038.  
  1039.  if((sign=a<0.0)==1) a=-a;
  1040.  
  1041.  for(i=1,dv=2;i<=T_P;i++,dv*=2)
  1042.   if(a>=1.0/dv) { a-=1.0/dv; res|=(0x1<<(T_P-i)); };
  1043.  
  1044.  if(sign==1) res=-res;
  1045.  
  1046.  return(res);
  1047. }
  1048.  
  1049.  
  1050.  
  1051. In the function above we are iteratively subtracting powers
  1052. of two from the argument and using the result to set the bits
  1053. in the fixed point number. We would use this function to
  1054. precalculate the array of sin/cos values. 
  1055.  
  1056. When calculating the rotational coefficients we would multiply
  1057. one fixed by another and obtain same precision fixed result.
  1058. again, actual result would have twice as many fractional bits
  1059. as each of the argument so same precision result would be
  1060. (a*b)>>T_P where 'a' and 'b' fixed point values and T_P number
  1061. of bits after the binary point in each. This way code would
  1062. look something like:
  1063.  
  1064.                                 ...
  1065.                                 ...
  1066.  
  1067.  T_wx1=((singam*((sinbet*sinalp)>>T_P))>>T_P) + ((cosgam*cosalp)>>T_P);
  1068.  T_wy1=((cosbet*sinalp)>>T_P);
  1069.  T_wz1=((singam*cosalp)>>T_P) - ((cosgam*((sinbet*sinalp)>>T_P))>>T_P);
  1070.  
  1071.                                 ...
  1072.                                 ...
  1073.  
  1074.  
  1075. In rotation function ints would have to be multiplied by a fixed
  1076. coefficient, the result would have T_P fractional bits, since we
  1077. are interested in just integer part same rule as above applies:
  1078.  
  1079.  
  1080.                                 ...
  1081.                                 ...
  1082.  
  1083.  *to++=((T_wx1*xt)>>T_P) + ((T_wy1*yt)>>T_P) + ((T_wz1*zt)>>T_P);
  1084.  *to++=((T_wx2*xt)>>T_P) + ((T_wy2*yt)>>T_P) + ((T_wz2*zt)>>T_P);
  1085.  *to++=((T_wx3*xt)>>T_P) + ((T_wy3*yt)>>T_P) + ((T_wz3*zt)>>T_P);
  1086.  
  1087.                                 ...
  1088.                                 ...
  1089.  
  1090. That's about it. The only thing, we really have to be careful with
  1091. the range of numbers the operations would work for. If we are working
  1092. with 32bit numbers and would choose 16 of them for fractional part
  1093. and one for sign bit then only 15bit would carry the integer part.
  1094. after multiplication of two fixed numbers since result has as many
  1095. fractional bits as in both arguments all 32 bits would actually carry
  1096. fractions. The integer part would be gone in the dark realm of left-carry,
  1097. and there's no way in standard straight C to get it out of there.
  1098. (when using assembly we don't have this problem since result
  1099. of multiplication has twice bits of each argument, so we really can
  1100. get trough to the integer part, however, most C compilers define the
  1101. result of multiplication to be of the same bit size as arguments
  1102. int*int==int, some compilers try to be smarter then that but we really
  1103. have to consider worst-case scenario)
  1104.  
  1105. And finally is it really true that integer multiplication plus shifts
  1106. is faster then floating multiplication? What follows are results of
  1107. some very approximate tests:
  1108.  
  1109.  
  1110.  
  1111. sparc:   floating faster fixed   about  1.3 times
  1112. 68040:   floating faster fixed   about  1.5 times
  1113. rs4000:  floating faster fixed   about  1.1 times
  1114.  
  1115. rs6000:  fixed faster floating   about  1.1 times
  1116. i386sx:  fixed faster floating   about  5.1 times
  1117.  
  1118. floating in the test was:    float a,b;    a*b;
  1119. fixed in the test was:       int a,b;      (a*b)>>15;
  1120.  
  1121.  
  1122.  
  1123. As one can see processors with fast floating point units nowadays
  1124. do floating operations faster then integer once. Then again in terms
  1125. of portability integer multiplication would always be fast
  1126. enough whereas floating multiplication especially on processors
  1127. without FPUs (i386sx) might get really slow. The decision whether
  1128. to implement fixed point, this way, might get a bit tough to make.
  1129.  
  1130. In the provided library's source transformations are implemented 
  1131. both ways but with the same function prototypes, so that one can
  1132. decide which implementation is more suitable to link in.
  1133.  
  1134.                                   * * *
  1135.  
  1136.  
  1137.  
  1138.  
  1139.                              2-D rasterization.
  1140.                             --------------------
  1141.  
  1142.  
  1143.  
  1144.  
  1145.  
  1146.  
  1147.                                 Re graphics: A picture is worth 10K words -
  1148.                                 but only those to describe the picture.
  1149.                                 Hardly any sets of 10K words can be adequately
  1150.                                 described with pictures.
  1151.  
  1152.  
  1153.  
  1154.  
  1155.  
  1156.  
  1157. Rasterization, the term itself refers to the technology we are using to
  1158. represent images: turning it into mosaic of small square pixel cells set
  1159. to different colours. The algorithms performing just that are covered in
  1160. this section.
  1161.  
  1162. A dot.
  1163. ------
  1164. Rasterizing dot is very straightforward, I am not even sure it is
  1165. warranted to call it rasterization at all. All it involves is just setting
  1166. some location in the screen memory to particular colour. Let's recall
  1167. that the very first location of the colourmap describes top rightmost
  1168. pixel, that determines the system of references which is convenient in
  1169. this case with X axis directed to the right and Y axis directed downward.
  1170.  
  1171.  
  1172.                                HW_X_SIZE
  1173.  
  1174.                    (0,0)+-+-+-+- ... -+------->  X
  1175.                         | | |         |
  1176.                         +-+-+         +
  1177.                         |
  1178.             HW_Y_SIZE    ...
  1179.  
  1180.                         |
  1181.                         +-+-      +-+
  1182.                         | |       |*|<------y
  1183.                         +-+-+     +-+
  1184.                         |          ^
  1185.                         |          |
  1186.                         V          |
  1187.                       Y            x
  1188.  
  1189.  
  1190.                                 pic 3.1
  1191.  
  1192.  
  1193. Assuming that each pixel is represented by one byte, and that the dimensions
  1194. of the output surface is HW_SCREEN_X_SIZE by HW_SCREEN_Y_SIZE it is the
  1195. y*HW_SCREEN_X_SIZE+x location we have to access in order to plot a dot
  1196. at (x,y). Here's some code to do just that:
  1197.  
  1198.  
  1199. unsigned char *G_buffer;                    /* colourmap- offscreen buffer */
  1200.  
  1201. void G_dot(int x,int y,unsigned char colour)
  1202. {
  1203.  G_buffer[y*HW_SCREEN_X_SIZE+x]=colour;
  1204. }
  1205.  
  1206. A line.
  1207. -------
  1208. Line rasterization is a very important peace of code for the library. It
  1209. would also be used for polygon rendering, and few other routines. What we
  1210. have to do is set to some specified colour all dots on the line's path.
  1211. We can easily calculate coordinates of all dots belonging to this line:
  1212.  
  1213.                                          * (xend,yend) ---
  1214.                                        /                 |
  1215.                                       * (x,y)            dy
  1216.                                     / |
  1217.                                   /   |y                 |
  1218.                 (xstart,ystart) *-----+          ---------
  1219.                                    x
  1220.                                          |
  1221.                                 |        |
  1222.                                 |- -dx- -|
  1223.  
  1224.                                   pic 3.2
  1225.  
  1226. Remembering a little bit of high (high?) school geometry, it can be
  1227. seen that if:
  1228.                               dx = xend-xstart
  1229.                               dy = yend-ystart
  1230. then:
  1231.                       x     dx               x * dy
  1232.                      --- = ----    and  y = --------
  1233.                       y     dy                 dx
  1234.  
  1235. Since we want to set all pixels along the path we just have to take
  1236. all Xs in the range of [xstart, xend] and calculate respective Y. There
  1237. is a bit of a catch, we would have only as much as xend-xstart
  1238. calculated coordinates. But if the Y range was actually longer the path
  1239. we have just found won't be continuous pic 3.3. In this case we should have
  1240. been calculating coordinates taking values from the longer range
  1241. [ystart, yend] and getting respective X, pic 3.4 as:
  1242.  
  1243.                                   y * dx
  1244.                              x = --------
  1245.                                     dy
  1246.  
  1247.               ....*                             ....*
  1248.               .....                             ....*
  1249.               ...*.                             ...*.
  1250.               .....                             ...*.
  1251.               ..*..  yrange                     ..*..  yrange
  1252.               .....                             ..*..
  1253.               .*...                             .*...
  1254.               .....                             .*...
  1255.               *....                             *....
  1256.  
  1257.               xrange                            xrange
  1258.  
  1259.              pic 3.3                           pic 3.4
  1260.  
  1261.  
  1262. What we just did is quite logical and easily programmable algorithm,
  1263. but... it is not that commonly used, the reason has to do with the way
  1264. we did the calculation, it involved division, and that is by far the
  1265. slowest operation any computer does. What we would do instead doesn't
  1266. involve a single division. It is called a Bresenham's integer based
  1267. algorithm (Bresenham, 1965) and here we find all points using few 
  1268. logical operations and integer addition. The idea is to find a bigger 
  1269. range among [xstart,xend] and [ystart,yend] go along points in it and 
  1270. have a variable saying when it is time to advance in the other range.
  1271. Let's consider what is happening at some stage in line rasterization,
  1272. let's assume X is a longer range (dx>dy) and that with the line we are
  1273. rendering xstart<xend  and  ystart<yend:
  1274.  
  1275.                                   |
  1276.                               |       | H(x+1,y+1)
  1277.                              -+---|---*-
  1278.                               |       |h/
  1279.                         |     |   |   * I(x+1,real_y)
  1280.                               |      /|
  1281.                         + - - - - + - - - - + -
  1282.                               |   /   |
  1283.                         |     | / |   |l    |
  1284.                               /       |
  1285.                              -*---|---*-L(x+1,y)
  1286.                               |P(x,y) |
  1287.  
  1288.  
  1289.                                pic 3.5
  1290.  
  1291.  
  1292. Pic 3.5 shows the situation where we just rendered a pixel P (Previous) at
  1293. coordinates (x,y) and now have to make a decision where to go to next, to
  1294. pixel L(x+1,y) (Lower) or H(x+1,y+1) (Higher). points P,H,L actually represent
  1295. middles of respective pixels. Since actual line would pass somewhere in the
  1296. middle between those two points we must plot the one whose centre is closer
  1297. to the intersection point I(x+1,real_y). this can be measured by comparing
  1298. segments h and l, resulted from the intersection: remembering dependency
  1299. used in the previous method we can find that:
  1300.  
  1301.                                      dy * (x+1)
  1302.                            real_y = ------------
  1303.                                          dx
  1304. and:
  1305.                                                       dy * (x+1)
  1306.                  h = y+1-real_y     =>     h = y+1 - ------------
  1307.                                                           dx
  1308.  
  1309.                                                 dy * (x+1)
  1310.                  l = real_y-y       =>     l = ------------ - y
  1311.                                                     dx
  1312. Now, we are interested in comparing l and h:
  1313.  
  1314.                                      dy * (x+1)
  1315.                           l - h = 2 ------------  - 2y-1
  1316.                                         dx
  1317.  
  1318.  
  1319. So, if l-h>0 it means that l>h the intersection point L is closer to point
  1320. H and the later should be rendered, otherwise L should be rendered. let's
  1321. multiply both sides by dx:
  1322.  
  1323.                        dx(l - h)= 2dyx + 2dy - 2dxy - dx
  1324.  
  1325. Now, since dx assumed bigger then 0, signs of (l-h) and dx(l-h) would be
  1326. the same. Let's call dx(l-h) as d and find it's sign and value at some step
  1327. i and next step i+1:
  1328.  
  1329.                       d    = 2dyx    -2dxy    + 2dy - dx
  1330.                        i         i-1      i-1
  1331.  
  1332.                       d    = 2dyx    -2dxy    + 2dy - dx
  1333.                        i+1       i        i
  1334.  
  1335. We can also find the initial value d, in the very first point since before
  1336. that x==0 and y==0:
  1337.  
  1338.                                d  = 2dy-dx
  1339.                                 1
  1340.                               -------------
  1341.  
  1342. Since sign of d determines which point to move next (H or L) lets find
  1343. value of d at some step i assuming we know it's value at i-1 from
  1344. previous calculations:
  1345.  
  1346.               d   - d  = 2dyx - 2dxy - 2dyx    + 2dxy
  1347.                i+1   i       i      i      i-1       i-1
  1348.  
  1349.               d   - d  = 2dy(x - x   ) - 2dx(y - y   )
  1350.                i+1   i        i   i-1         i   i-1
  1351.  
  1352. Depending on the decision taken in the previous point value of d in the
  1353. next one would be either:
  1354.  
  1355.                when H was chosen since x - x    =1 and y - y   =1
  1356.                                          i   i-1         i   i-1
  1357.  
  1358.                                 d   - d = 2dy - 2dx
  1359.                                  i+1   i
  1360.  
  1361.                                 d   = d + 2dy - 2dx
  1362.                                  i+1   i
  1363.                                ---------------------
  1364.  
  1365. or:
  1366.  
  1367.                when L was chosen since x - x   =1 and y - y   =0
  1368.                                          i   i-1        i   i-1
  1369.  
  1370.                                 d   -d = 2dy
  1371.                                  i+1  i
  1372.                                 d   =d + 2dy
  1373.                                  i+1  i
  1374.                                --------------
  1375.  
  1376. This may seem a little complicated but the implementation code certainly
  1377. doesn't look this way. We just have to find the initial value of d, and
  1378. then in the loop depending on it's sign add to it either 2dy-2dx, or
  1379. just 2dy. since those are constants, they can be calculated before the
  1380. actual loop. In the proof above we have assumed xstart<yend ystrat<yend,
  1381. however, in real life we can't guarantee that. The way we can take care
  1382. of it, is just changing increments to decrements when above assumptions
  1383. don't hold. We also have to remember that in the loop when L point is
  1384. picked it is only value belonging to the bigger range which is incremented
  1385. (decremented) whereas when picking H point both X and Y have to be changed
  1386. since that is when we are advancing along both axes.
  1387.  
  1388.  
  1389. void G_line(int x,int y,int x2,int y2,unsigned char colour)
  1390. {
  1391.  int dx,dy,long_d,short_d;
  1392.  int d,add_dh,add_dl;
  1393.  register int inc_xh,inc_yh,inc_xl,inc_yl;  
  1394.  register int i;                            
  1395.  
  1396.  dx=x2-x; dy=y2-y;                          /* ranges */
  1397.  
  1398.  if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}    /* making sure dx and dy >0 */
  1399.  else    {        inc_xh=1;  inc_xl=1; }    /* adjusting increments */
  1400.  if(dy<0){dy=-dy; inc_yh=-1; inc_yl=-1;}
  1401.  else    {        inc_yh=1;  inc_yl=1; }
  1402.  
  1403.  if(dx>dy){long_d=dx; short_d=dy; inc_yl=0;}/* long range,&making sure either */
  1404.  else     {long_d=dy; short_d=dx; inc_xl=0;}/* x or y is changed in L case */
  1405.  
  1406.  d=2*short_d-long_d;                        /* initial value of d */
  1407.  add_dl=2*short_d;                          /* d adjustment for H case */
  1408.  add_dh=2*short_d-2*long_d;                 /* d adjustment for L case */
  1409.  
  1410.  for(i=0;i<=long_d;i++)                     /* for all points in longer range */
  1411.  {
  1412.   G_dot(x,y,colour);                        /* rendering */
  1413.  
  1414.   if(d>=0){x+=inc_xh; y+=inc_yh; d+=add_dh;}/* previous point was H type */
  1415.   else    {x+=inc_xl; y+=inc_yl; d+=add_dl;}/* previous point was L type */
  1416.  }
  1417. }
  1418.  
  1419.  
  1420. As you can see this algorithm doesn't involve divisions at all. It does
  1421. have multiplication by two, but almost any modern compiler would optimize
  1422. it to 1 bit left shift - a very cheap operation. Or if we don't have faith
  1423. in the compiler we can do it ourselves. Lets optimize this code a little
  1424. bit. Whenever trying to optimize something we first have to go after
  1425. performance of the code in the loop since it is something being done
  1426. over and over. In this source above we can see the function call to G_dot,
  1427. let's remember what is inside there,.... oh, it is y*HW_SCREEN_X_SIZE+x a
  1428. multiplication... an operation similar in length to division which
  1429. we just spent so much effort to avoid. Lets think back to the
  1430. representation of the colourmap, if we just rendered a point and
  1431. now want to render another one next to it how their addresses in the
  1432. colourmap would be different? pic 3.6
  1433.  
  1434.                             +-+-+-
  1435.                            A|*->|B
  1436.                             +|+-+-
  1437.                            C|V| |
  1438.                             +-+-
  1439.                             | |
  1440.  
  1441.                             pic 3.6
  1442.  
  1443. Well, if it is along X axis that we want to advance we just have to add 1,
  1444. to the address of pixel A to get to pixel B, if we want to advance along
  1445. Y we have to add to the address of A number of bytes in the colourmap
  1446. separating A and C. this number is exactly HW_SCREEN_X_SIZE that is length
  1447. of one horizontal line of pixels in memory. Now, using that, let's try
  1448. instead of coordinates (x,y) to calculate it's address in the colourmap:
  1449.  
  1450.  
  1451. void G_line(int x,int y,int x2,int y2,unsigned char colour)
  1452. {
  1453.  int dx,dy,long_d,short_d;
  1454.  int d,add_dh,add_dl;
  1455.  int inc_xh,inc_yh,inc_xl,inc_yl;
  1456.  register int inc_ah,inc_al;          
  1457.  register int i;                      
  1458.  register unsigned char *adr=G_buffer;
  1459.  
  1460.  dx=x2-x; dy=y2-y;                          /* ranges */
  1461.  
  1462.  if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}    /* making sure dx and dy >0 */
  1463.  else    {        inc_xh=1;  inc_xl=1; }    /* adjusting increments */
  1464.  if(dy<0){dy=-dy; inc_yh=-HW_SCREEN_X_SIZE;
  1465.                   inc_yl=-HW_SCREEN_X_SIZE;}/* to get to the neighboring */
  1466.  else    {        inc_yh= HW_SCREEN_X_SIZE; /* point along Y have to add */
  1467.                   inc_yl= HW_SCREEN_X_SIZE;}/* or subtract this */
  1468.  
  1469.  if(dx>dy){long_d=dx; short_d=dy; inc_yl=0;}/* long range,&making sure either */
  1470.  else     {long_d=dy; short_d=dx; inc_xl=0;}/* x or y is changed in L case */
  1471.  
  1472.  inc_ah=inc_xh+inc_yh;
  1473.  inc_al=inc_xl+inc_yl;                      /* increments for point adress */
  1474.  adr+=y*HW_SCREEN_X_SIZE+x;                 /* address of first point */
  1475.  
  1476.  d=2*short_d-long_d;                        /* initial value of d */
  1477.  add_dl=2*short_d;                          /* d adjustment for H case */
  1478.  add_dh=2*short_d-2*long_d;                 /* d adjustment for L case */
  1479.  
  1480.  for(i=0;i<=long_d;i++)                     /* for all points in longer range */
  1481.  {
  1482.   *adr=colour;                              /* rendering */
  1483.  
  1484.   if(d>=0){adr+=inc_ah; d+=add_dh;}         /* previous point was H type */
  1485.   else    {adr+=inc_al; d+=add_dl;}         /* previous point was L type */
  1486.  }
  1487. }
  1488.  
  1489.  
  1490. We can see from this code that although the complexity of the preliminary
  1491. part of the algorithm went up a bit, the code inside the loop is much shorter
  1492. and simpler. There is just one thing left to add in order to make it all
  1493. completely usable - clipping that is making sure our drawing doesn't go
  1494. outside the physical boundaries of the colourmap, but this is covered in
  1495. the next chapter.
  1496.  
  1497.  
  1498.  
  1499. Ambient polygons
  1500. ----------------
  1501. What we have to do is to fill with some colour inside the polygon's edges.
  1502. The easiest and perhaps fastest way to do it is by using "scan line"
  1503. algorithm. The idea behind it is to convert a polygon into a collection
  1504. of usually horizontal lines and then render them one at a time. The lines
  1505. have to be horizontal only because the most common colourmap would have
  1506. horizontal lines contingent in memory which automatically allows to render
  1507. them faster then vertical lines since increment by 1 is usually faster
  1508. then addition. There is one inherited limitation to the algorithm, because
  1509. at each Y heights there would be just one horizontal line only concave
  1510. polygons can be rendered correctly. It can be seen that non concave polygons
  1511. might need two or more lines sometimes pic 3.7 (that's by the way, in case
  1512. you wondered, determines the difference in definitions between concave and
  1513. non-concave polygons):
  1514.  
  1515.                              *\
  1516.                              |  \    /*
  1517.                             y|---\  /-|
  1518.                              |    \/  |
  1519.                              *--------*
  1520.  
  1521.                                pic 3.7
  1522.  
  1523. How can we represent the collection of horizontal lines? obviously we
  1524. just need to know at which Y heights it is, X coordinate of the point
  1525. where it starts and another X coordinate of the point where it ends.
  1526. So lets have one "start" and one "end" location for every Y possible
  1527. on the screen. Since every line is limited by two points each belonging
  1528. to certain edge, we can take one edge at a time find all points belonging
  1529. to it and for each of them change the "start" and "end" at particular Y
  1530. heights. So that at the end we would have coordinates of all scan lines
  1531. in the polygon pic 3.8:
  1532.  
  1533.                  start end
  1534.                 +-----+---+              1 2 3 4 5 6 7 8
  1535.                 |  1  | 5 |- - - - - ->  *------\
  1536.                 |  2  | 8 |               \------------*
  1537.                 |  2  | 7 |- - - - - - ->  \----------/
  1538.                 |  3  | 7 |                  \------/
  1539.                 |  4  | 6 |- - - - - - - ->   *----*
  1540.                 +-----+---+
  1541.                                 pic 3.8
  1542.  
  1543. Initially we would put the biggest possible value in all locations of
  1544. "start" array and the smallest possible in the locations of "end" array.
  1545. (Those are by the way are defined in <limits.h> and it is not a bad
  1546. practice to use them). Each time we've calculated the point on the edge
  1547. we have to go to "start" and "end" locations at Y heights and if X of
  1548. the point less then what we have in "start" write it there. Similar to
  1549. that if this same X is bigger then value in "end" location we have to
  1550. put it there too. Because of the way the arrays were initialized first
  1551. point calculated for every Y height would change both "start" and "end"
  1552. location. Let's look at the code to find Xs for each edge, so called
  1553. scanning function:
  1554.  
  1555. void GI_scan(int *edges,int dimension,int skip)
  1556. {
  1557.  signed_32_bit cur_v[C_MAX_DIMENSIONS];     /* initial values for Z+ dims */
  1558.  signed_32_bit inc_v[C_MAX_DIMENSIONS];     /* increment for Z+ dimensions */
  1559.  signed_32_bit tmp;
  1560.  int dx,dy,long_d,short_d;
  1561.  register int d,add_dh,add_dl;              
  1562.  register int inc_xh,inc_yh,inc_xl,inc_yl;  
  1563.  int x,y,i,j;
  1564.  int *v1,*v2;                               /* first and second verteces */
  1565.  
  1566.  v1=edges; edges+=skip; v2=edges;           /* length ints in each */
  1567.  
  1568.  if(C_line_y_clipping(&v1,&v2,dimension))   /* vertical clipping */
  1569.  {
  1570.   dx=*v2++; dy=*v2++;                       /* extracting 2-D coordinates */
  1571.   x=*v1++; y=*v1++;                         /* v2/v1 point remaining dim-2 */
  1572.   dimension-=2;
  1573.  
  1574.   if(y<G_miny) G_miny=y;
  1575.   if(y>G_maxy) G_maxy=y;
  1576.   if(dy<G_miny) G_miny=dy;
  1577.   if(dy>G_maxy) G_maxy=dy;                  /* updating vertical size */
  1578.  
  1579.   dx-=x; dy-=y;                             /* ranges */
  1580.  
  1581.   if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}   /* making sure dx and dy >0 */
  1582.   else    {        inc_xh=1;  inc_xl=1; }   /* adjusting increments */
  1583.   if(dy<0){dy=-dy; inc_yh=-1; inc_yl=-1;}
  1584.   else    {        inc_yh=1;  inc_yl=1; }
  1585.  
  1586.   if(dx>dy){long_d=dx;short_d=dy;inc_yl=0;} /* long range,&making sure either */
  1587.   else     {long_d=dy;short_d=dx;inc_xl=0;} /* x or y is changed in L case */
  1588.  
  1589.   d=2*short_d-long_d;                       /* initial value of d */
  1590.   add_dl=2*short_d;                         /* d adjustment for H case */
  1591.   add_dh=2*short_d-2*long_d;                /* d adjustment for L case */
  1592.  
  1593.   for(i=0;i<dimension;i++)                  /* for all remaining dimensions */
  1594.   {
  1595.    tmp=v1[i]; tmp<<=G_P; cur_v[i]=tmp;      /* f.p. start value */
  1596.    tmp=v2[i]-v1[i]; tmp<<=G_P;              /* f.p. increment */
  1597.    if(long_d>0) inc_v[i]=tmp/long_d;        /* if needed (non 0 lines) */
  1598.   }
  1599.  
  1600.   for(i=0;i<=long_d;i++)                    /* for all points in longer range */
  1601.   {
  1602.    if(x<G_start[y])                         /* further then rightmost */
  1603.    {
  1604.     G_start[y]=x;                           /* the begining of scan line */
  1605.     for(j=0;j<dimension;j++)
  1606.      G_rest_start[j][y]=cur_v[j];            /* all other dimensions */
  1607.    }
  1608.  
  1609.    if(G_end[y]<x)                           /* further the leftmost */
  1610.    {
  1611.     G_end[y]=x;                             /* the end of scan line */
  1612.     for(j=0;j<dimension;j++)
  1613.      G_rest_end[j][y]=cur_v[j];              /* and for other dimension */
  1614.    }
  1615.  
  1616.    if(d>=0){x+=inc_xh;y+=inc_yh;d+=add_dh;} /* previous dot was H type */
  1617.    else    {x+=inc_xl;y+=inc_yl;d+=add_dl;} /* previous dot was L type */
  1618.    for(j=0;j<dimension;j++)
  1619.     cur_v[j]+=inc_v[j];                     /* for all other dimensions */
  1620.   }
  1621.  }
  1622. }
  1623.  
  1624. As you can see this edge scanning function doesn't have much differences
  1625. with our most basic line rendering code. The only thing we are updating the
  1626. vertical size of polygon so that after all edges would go through this
  1627. function we would have the correct vertical dimensions in G_miny and
  1628. G_maxy. There is another difference: scan function is designed to process,
  1629. other dimensions, that is it would interpolate X and similar to that
  1630. other values specified in vertices, finding them for all Ys in range.
  1631. We would need that when adding interpolative shading feature. As to the
  1632. ambient polygon rendering, After scanning is over we would just have
  1633. to render all horizontal scan lines in the range [G_miny,G_maxy]:
  1634.  
  1635. void G_ambient_polygon(int *edges,int length,unsigned char colour)
  1636. {
  1637.  int new_edges[G_MAX_TUPLES];               /* verticaly clipped polygon */
  1638.  int new_length;                            /* although no edges there yet */
  1639.  register unsigned char *l_adr,*adr;
  1640.  register int beg,end,i;
  1641.  
  1642.  GI_boarder_array_init();                   /* initializing the arrays */
  1643.  
  1644.  new_length=C_polygon_x_clipping(edges,new_edges,2,length);
  1645.  
  1646.  for(i=0;i<new_length;i++)
  1647.   GI_scan(&new_edges[i*2],2,2);             /* Searching polygon boarders */
  1648.  
  1649.  if(G_miny<=G_maxy)                         /* nothing to do? */
  1650.  {
  1651.   l_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  1652.  
  1653.   for(; G_miny<=G_maxy; G_miny++,
  1654.       l_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  1655.   {
  1656.    adr=l_adr+(beg=G_start[G_miny]);         /* addr of the current line */
  1657.    end=G_end[G_miny];                       /* ends here */
  1658.  
  1659.    for(;beg<=end;beg++) *adr++=colour;      /* rendering single line */
  1660.   }
  1661.  }
  1662. }
  1663.  
  1664. As to the optimization it doesn't appear we can do a lot without
  1665. leaving comfortable boundaries of C. On the other hand when there's
  1666. really nothing else we can do to speed things up, when all manuals
  1667. are read and all friends are out of ideas let's do it, let's go
  1668. for assembly. I don't think rewriting all functions this way is
  1669. appealing. Moreover only rewriting the real execution bottleneck
  1670. would give considerable speed increase. That's, perhaps, one of
  1671. the reason why modern compilers have such a nice feature as inline
  1672. assembly, allowing to add low-level code directly into C source.
  1673.  
  1674. Looking at the code above we can see this line rendering loop. This is
  1675. an obvious candidate, since it executes fairly often and does quite a
  1676. lot of iterations, especially compared with surrounding code.
  1677.  
  1678. Different C compilers have different provision for inline assembly.
  1679. GNU C compiler has quite powerful syntax allowing to specify assembly
  1680. code. It is:
  1681.  
  1682.      asm("assembly_command %1,%0":"=constraint" (result)
  1683.                                  :" constraint" (argument),
  1684.                                   " constraint" (argument),
  1685.                                            ...
  1686.                                  :" clobbered_reg", "clobbered_reg" ...
  1687.         );
  1688.  
  1689. The assembly instruction is specified with aliases allowing
  1690. to add names of C variables into it. it also allows to specify
  1691. the CPU registers the command, we are inserting, might clobber.
  1692. It is particularly important if we are expecting to compile with
  1693. optimizing option or have register variables in the code. We
  1694. don't want the program trying to use contents of a register we
  1695. just destroyed by inline assembly instruction. If we are specifying
  1696. the clobbered registers, on the other hand, compiler would make
  1697. sure that no useful information is in this registers at the
  1698. moment inline assembly code starts to execute.
  1699.  
  1700. WATCOM C compiler has different provision for inline assembly:
  1701.  
  1702. return_type c_function(parameter1,parametr2, ...);
  1703. #pragma aux c_fucntion=        \
  1704.            "assembly_command"  \
  1705.            "assembly_command"  \
  1706.            parm [reg] [reg]    \
  1707.            value [reg]         \
  1708.            modify [reg reg];
  1709.  
  1710. It is implemented by providing sort of a pseudo c-function which
  1711. is treated by the compiler pretty much like a macro. "parm" specifies
  1712. which registers take in parameters from function's argument list,
  1713. "value"- from what register, if any, the return value is taken and
  1714. "modify", specifies clobbered registers.
  1715.  
  1716. for DJGPP compiler code substituting line rendering loop:
  1717.  
  1718.    for(;beg<=end;beg++) *adr++=colour;      /* rendering single line */
  1719.  
  1720. may look like:
  1721.  
  1722.    asm("movl  %0,%%ecx" :: "g" (end):"%ecx");
  1723.    asm("subl  %0,%%ecx" :: "g" (beg):"%ecx");
  1724.    asm("movl  %0,%%edi" :: "g" (adr):"%edi");
  1725.    asm("movl  %0,%%eax" :: "g" (colour):"%eax");
  1726.    asm("cld");                              /* increment when stosb */
  1727.    asm("rep");
  1728.    asm("stosb %al,(%edi)");                 /* rendering line here */
  1729.  
  1730. for WATCOM:
  1731.  
  1732. void a_loop(unsigned char *adr,int end,int beg,int colour);
  1733. #pragma aux a_loop=                      \
  1734.         "sub ecx,esi"                    \
  1735.         "cld"                            \
  1736.         "rep stosb"                      \
  1737.         parm [edi] [ecx] [esi] [al];
  1738.  
  1739. Not a terribly complicated code to write but something giving
  1740. considerable speed increase. Then again we might have used a C
  1741. library function memset for filling chunks of memory:
  1742.  
  1743.    void *memset(void *s,int c,size_t n);
  1744.  
  1745.    memset(l_adr,colour,end-beg);            /* filling scan line */
  1746.  
  1747. This function might have had just that code, we wrote, inside of
  1748. it. let's not believe assumptions neither pro nor con and go check
  1749. that out. with DJGPP or, in this matter, most UNIX compilers we can
  1750. just extract and disassemble memset code from libc.a (standard
  1751. c library) and see for ourselves. let's go:
  1752.  
  1753.      ar -x libc.a memset.o
  1754.      objdump -d memset.o
  1755.  
  1756. We would see something like:
  1757.  
  1758. memset.o:     file format coff-go32
  1759.  
  1760. Disassembly of section .text:
  1761. 00000000 <_memset> pushl  %ebp
  1762. 00000001 <_memset+1> movl   %esp,%ebp
  1763. 00000003 <_memset+3> pushl  %edi
  1764. 00000004 <_memset+4> movl   0x8(%ebp),%edi
  1765. 00000007 <_memset+7> movl   0xc(%ebp),%eax
  1766. 0000000a <_memset+a> movl   0x10(%ebp),%ecx
  1767. 0000000d <_memset+d> jcxz   00000011 <_memset+11>
  1768. 0000000f <_memset+f> repz stosb %al,%es:(%edi)
  1769. 00000011 <_memset+11> popl   %edi
  1770. 00000012 <_memset+12> movl   0x8(%ebp),%eax
  1771. 00000015 <_memset+15> leave
  1772. 00000016 <_memset+16> ret
  1773. 00000017 <_memset+17> Address 0x18 is out of bounds.
  1774. Disassembly of section .data:
  1775.  
  1776. Very similar to what we wrote? Other compilers might have
  1777. worse library code? it might be so, let's try WATCOM C,
  1778. let's go:
  1779.  
  1780.      wlib clibs.lib *memset
  1781.      wdisasm memset.obj
  1782.  
  1783. that's what we would see:
  1784.  
  1785. Module: memset
  1786. Group: 'DGROUP' CONST,CONST2,_DATA,_BSS
  1787.  
  1788. Segment: '_TEXT' BYTE USE32  00000023 bytes
  1789.  0000  57                memset          push    edi
  1790.  0001  55                                push    ebp
  1791.  0002  89 e5                             mov     ebp,esp
  1792.  0004  8b 45 10                          mov     eax,+10H[ebp]
  1793.  0007  8b 4d 14                          mov     ecx,+14H[ebp]
  1794.  000a  8b 7d 0c                          mov     edi,+0cH[ebp]
  1795.  000d  06                                push    es
  1796.  000e  1e                                push    ds
  1797.  000f  07                                pop     es
  1798.  0010  57                                push    edi
  1799.  0011  88 c4                             mov     ah,al
  1800.  0013  d1 e9                             shr     ecx,1
  1801.  0015  f2 66 ab                          repne   stosw
  1802.  0018  11 c9                             adc     ecx,ecx
  1803.  001a  f2 aa                             repne   stosb
  1804.  001c  5f                                pop     edi
  1805.  001d  07                                pop     es
  1806.  001e  89 f8                             mov     eax,edi
  1807.  0020  c9                                leave
  1808.  0021  5f                                pop     edi
  1809.  0022  c3                                ret
  1810.  
  1811. No disassembly errors
  1812.  
  1813. ------------------------------------------------------------
  1814.  
  1815. Is it not very similar to what we just did using inline assembly?
  1816. WATCOM C code is even trying to be smarter then us by storing as
  1817. much as it can by word store instruction, since every word store
  1818. in this case is equivalent to two char stores, there should be
  1819. twice less iterations. Of course there would be some time spent
  1820. on the function call itself but most likely performance of the
  1821. inline assembly code we wrote and memset call would be very very
  1822. similar. It is just yet another example of how easy it might be
  1823. in some cases to achieve results by very simple means instead of
  1824. spending sleepless hours in front of the glowing box (unfortunately
  1825. only sleepless hours teach us how to do things by simple means but
  1826. that's the topic for another story only marginally dedicated to
  1827. 3-D graphics).
  1828.  
  1829. By the way, in the provided source code there is an independent
  1830. approach to fast memory moves and fills. That is functions to perform
  1831. those are brought out into hardware interface as: HW_set_char ;
  1832. HW_set_int; HW_copy_char; HW_copy_int. And their implementation
  1833. may actually be either of the ways described above depending on
  1834. the particular platform.
  1835.  
  1836.  
  1837. Shaded polygons.
  1838. ----------------
  1839.  
  1840. Shaded polygons are used to smooth appearance of faceted models,
  1841. that is those where we approximate real, curved surfaces by a
  1842. set of plane polygons. Interpolative or Gouraud shading (Gouraud, 1971), 
  1843. we would discuss, is the easiest to implement, it is also not very 
  1844. expensive in terms of speed. The idea is that we carry a colour 
  1845. intensity value in every vertex of the polygon and linearly interpolate 
  1846. those values to find colour for each pixel inside the polygon.
  1847.  
  1848. Finally this is the place where implemented N-dimensional scanning
  1849. procedure comes in handy. Indeed colour intensity can be pretty
  1850. much handled as just an extra dimension as far as edge scanning
  1851. and clipping are concerned.
  1852.  
  1853.                              *I3
  1854.                             /  \
  1855.                            /      \
  1856.                         I0*          \
  1857.                            \           *I2
  1858.                             \    I    /
  1859.                           IT1*---*---*IT2
  1860.                               \     /
  1861.                                \   /
  1862.                                 \ /
  1863.                                  *I1
  1864.  
  1865.                               pic 3.9
  1866.  
  1867. We can obtain values on the left and right boarders of a polygon (IT1,
  1868. IT2 on pic 3.9) by interpolating colour intensity value in every edge
  1869. during edge scanning, just as "start/end' X values were found for every
  1870. horizontal line. Afterwards when rendering certain horizontal scan line we
  1871. can further interpolate "start" and "end" intensities for this line finding
  1872. colour for pixels belonging to it (I on pic 3.9). Fixed point math is
  1873. probably a very convenient way to implement this interpolation:
  1874.  
  1875. void G_shaded_polygon(int *edges,int length)
  1876. {
  1877.  int new_edges[G_MAX_SHADED_TUPLES];        /* verticaly clipped polygon */
  1878.  int new_length;
  1879.  register unsigned char *l_adr,*adr;
  1880.  register int beg,end,i;
  1881.  register signed_32_bit cur_c,inc_c;        /* current colour and it's inc */
  1882.  
  1883.  GI_boarder_array_init();                   /* initializing the array */
  1884.  
  1885.  new_length=C_polygon_x_clipping(edges,new_edges,3,length);
  1886.  
  1887.  for(i=0;i<new_length;i++)
  1888.   GI_scan(&new_edges[i*3],3,3);             /* Searching polygon boarders */
  1889.  
  1890.  if(G_miny<=G_maxy)                         /* nothing to do? */
  1891.  {
  1892.   l_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  1893.  
  1894.   for(; G_miny<=G_maxy; G_miny++,
  1895.       l_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  1896.   {
  1897.    adr=l_adr+(beg=G_start[G_miny]);         /* addr of the current line */
  1898.    end=G_end[G_miny];                       /* ends here */
  1899.  
  1900.    cur_c=G_0_start[G_miny];                 /* colour starts with this value */
  1901.    inc_c=G_0_end[G_miny]-cur_c;
  1902.    if(end>beg) inc_c/=end-beg+1;            /* f.p. increment */
  1903.  
  1904.    for(;beg<=end;beg++)                     /* render this line */
  1905.    {
  1906.     *adr++=cur_c>>G_P;                      /* rendering single point */
  1907.     cur_c+=inc_c;                           /* incrementing colour */
  1908.    }
  1909.   }
  1910.  }
  1911. }
  1912.  
  1913. As you see code above looks not that much different from ambient polygon
  1914. rendering source.
  1915.  
  1916. There is one small catch. I said that intensities can be handled pretty
  1917. much as an extra dimension. In fact we can consider shaded polygon on
  1918. the screen to take 2 space dimensions, and to have one colour dimension.
  1919. But would all vertexes of this polygon belong to one plane in such 3-D space?
  1920. If they would not, shading would look rather weird especially when this
  1921. polygon would start to rotate. The common solution is limit polygons to just
  1922. triangles. Why? well three points always belong to the same plane in
  1923. 3-D space.
  1924.  
  1925.  
  1926.  
  1927. Textured polygons.
  1928. ------------------
  1929.  
  1930. It might seem that we can extend interpolation technique used to
  1931. calculate colour intensity for shading, onto texture rendering. Just
  1932. keep texture X and Y in every vertex of the polygon, interpolate
  1933. those two along edges and then along horizontal lines obtaining
  1934. this way texture (X,Y) for every pixel to be rendered. This however 
  1935. does not work... Or to be more exact it stops working when perspective 
  1936. projection is used, the reason is simple: perspective transformation 
  1937. is not linear.
  1938.  
  1939. One may ask, why did it work for interpolative shading, after all
  1940. we were using perspective transformation all along? The answer is
  1941. it did not work, but visual aspect of neglecting perspective effect
  1942. for shading is quite suitable, neglecting it for texture rendering
  1943. however is not. I believe it has to do with different image frequencies.
  1944.  
  1945. Just consider the following situation: a rectangular polygon is
  1946. displayed on the screen so that it's one side is much closer
  1947. to us then the other one, which is almost disappearing in the
  1948. infinity: where do you think in the texture map lines A,B,C and D
  1949. would be mapped from by linear interpolating method? how do you think
  1950. the texture on the polygon would look like?
  1951.  
  1952.  
  1953.             +\                       +--/---/+1
  1954.            A|--\                    A|/  /   |
  1955.            B|----\                  B|/      |
  1956.            C|----/                  C|\    <------ Missed area
  1957.            D|--/                    D|\  \   |
  1958.             +/                       +--\---\+2
  1959.  
  1960.        Polygon on the                 Texture
  1961.           Screen                        Map
  1962.  
  1963.  
  1964.                          pic 3.10
  1965.  
  1966. Well it would look like if someone has carved triangular area
  1967. marked "missed" on the picture above, and connected points 1 and
  1968. 2 together, real ugly... In less dramatic cases texture rendered this
  1969. way might look nicer, being perfect when effect of perspective
  1970. transformation is negligible: all points are at almost the same
  1971. distance from the viewer, or very far away from the viewing plane.
  1972.  
  1973. To get always nice looking results we would have to consider
  1974. what is really happening when texture is being rendered:
  1975.  
  1976.  
  1977.                                  u^                  Texture Space
  1978.                                   | *T(u,v)
  1979.                                   +---->
  1980.                                  o     v
  1981.  
  1982.                          y ^
  1983.                            |       z
  1984.                            |  U \ /*W(x,y,z)/ V      World Space
  1985.                            |    / \       /
  1986.                            |  /     \   /
  1987.                            |/         * O
  1988.                            +--------------> x
  1989.  
  1990.                                 j^
  1991.                                  |
  1992.                                  |    *S(i,j)        Screen Space
  1993.                                  |
  1994.                                  +------->i
  1995.  
  1996.  
  1997.                                    pic 3.11
  1998.  
  1999. Consider three spaces pictured above: the texture space, world space (the
  2000. one before projection) and finally contents of the screen - screen or image
  2001. space. We may think of them just as original polygon/texture map in case
  2002. of texture space, polygon/texture map turned somehow in case of world
  2003. space and finally same projected onto the screen.
  2004.  
  2005. If we know mapping of origin and two orthogonal vectors from texture
  2006. space into the world space, we can express mapping of any point through
  2007. them:
  2008.  
  2009.                            x = Ox + v*Vx + u*Ux
  2010.                            y = Oy + v*Vy + u*Uy
  2011.                            z = Oz + v*Vz + u*Uz
  2012.  
  2013. Where Vx...Uz are corresponding components of the respective vectors.
  2014. Further point in the world space W(x,y,z) is being perspectively
  2015. transformed into the screen space S(i,j):
  2016.  
  2017.                                  i = x / z
  2018.                                  j = y / z
  2019.  
  2020. (the actual transformation used, would most likely also contain a
  2021. multiplier - focus, but let's worry about particularities on some later
  2022. stage). In order to perform mapping, on the other hand, we need u,v
  2023. texture coordinates as function of screen coordinates i,j
  2024.  
  2025.                                  x = i * z
  2026.                                  y = i * z
  2027.  
  2028. or:
  2029.                   Ox + v*Vx + u*Ux = i * [ Oz + v*Vz + u*Uz ]
  2030.                   Oy + v*Vy + u*Uy = j * [ Oz + v*Vz + u*Uz ]
  2031.  
  2032. trying to express u,v through i,j:
  2033.  
  2034.                     v*(Vx-i*Vz) + u*(Ux-i*Uz) = i*Oz - Ox
  2035.                     v*(Vy-j*Vz) + u*(Uy-j*Uz) = j*Oz - Oy
  2036.  
  2037. further:
  2038.                           (i*Oz - Ox) - u*(Ux - i*Uz)
  2039.                      v = -----------------------------
  2040.                                  (Vx - i*Vz)
  2041.  
  2042.                           (i*Oz - Ox) - v*(Vx - i*Vz)
  2043.                      u = -----------------------------
  2044.                                  (Ux - i*Uz)
  2045. and:
  2046.  
  2047.      (i*Oz - Ox) - u*(Ux - i*Uz)
  2048.     ----------------------------- * (Vy - j*Vz) + u*(Uy - j*Uz) = j*Oz - Oy
  2049.               (Vx - i*Vz)
  2050.  
  2051.                      (i*Oz - Ox) - v*(Vx - i*Vz)
  2052.     v*(Vy - j*Vz) + ----------------------------- * (Uy - j*Uz) = j*Oz - Oy
  2053.                               (Ux - i*Uz)
  2054.  
  2055. or in nicer form:
  2056.  
  2057.           i*(Vz*Oy-Vy*Oz) + j*(Vx*Oz-Vz*Ox) + (Vy*Ox-Vx*Oy)
  2058.       u = --------------------------------------------------
  2059.           i*(Vy*Uz-Vz*Uy) + j*(Vz*Ux-Vx*Uz) + (Vx*Uy-Vy*Ux)
  2060.  
  2061.           i*(Uy*Oz-Uz*Oy) + j*(Uz*Ox-Ux*Oz) + (Ux*Oy-Uy*Ox)
  2062.       v = --------------------------------------------------
  2063.           i*(Vy*Uz-Vz*Uy) + j*(Vz*Ux-Vx*Uz) + (Vx*Uy-Vy*Ux)
  2064.  
  2065. At this point we have formulas of from where in the texture space the
  2066. point in the screen space originates. There are several implementational
  2067. problems, first few are simple, the actual perspective transformation
  2068. is  i=x*focus/z rather then i=x/z, plus, we are performing screen centre
  2069. translation so that screen space (0,0) is in the middle of the screen
  2070. and not in the upper left corner. Hence original dependency should have
  2071. been:
  2072.  
  2073.                i = (x * focus) / z + HW_SCREEN_X_CENTRE
  2074.                j = (y * focus) / z + HW_SCREEN_Y_CENTRE
  2075.  
  2076.  
  2077. And so finally, so called, reverse mapping formulas have to be amended 
  2078. respectively:
  2079.  
  2080.                  ((i-HW_SCREEN_X_CENTRE)/focus) * (Vz* ...
  2081.             u = -------------------------------------------
  2082.                                 ...
  2083.  
  2084.                                 ...
  2085.  
  2086.  
  2087. Another, a bit difficult part has to do with texture scaling.
  2088. If the size of texture vectors was picked to be 1 that would
  2089. assume no scaling, that is unit size along the polygon in the
  2090. world space would correspond to unit size in the texture space.
  2091. What we however want to do in a lot of instances is to map
  2092. smaller texture onto a bigger polygon or the other way around,
  2093. bigger texture onto a smaller polygon. Let's try scaling the
  2094. texture space:
  2095.  
  2096.        T
  2097.     *---------->
  2098.    v| *        |         S
  2099.     | | \          *--------------------->
  2100.     | |    \   | v'| *                   |       v       T
  2101.     | *-- ---*     | | \                       ----- = -----
  2102.     V- - - - - |   | |   \               |       v'      S
  2103.                    | |      \
  2104.                    | |mapping \          |           v'*T
  2105.                    | |of the      \              v = -----
  2106.                    | |polygon      \     |             S
  2107.                    | *---------------*
  2108.                    V- - - - - - - - - - -|
  2109.  
  2110.                              pic 3.12
  2111.  
  2112.  
  2113.  
  2114. now we can put expression for v into direct mapping equations:
  2115.  
  2116.  
  2117.                           Vx*T         Ux*T
  2118.             x = Ox + v'* ------ + u'* ------
  2119.                            S            S
  2120.  
  2121.                                ...
  2122.  
  2123. Vx and Ux multiplied by T in fact are just projections of vectors 
  2124. enclosing the hole texture space, let's call them V'x == Vx*T, 
  2125. U'x == Ux*T etc.
  2126.  
  2127.  
  2128.                           V'x         U'x
  2129.             x = Ox + v'* ----- + u'* -----
  2130.                            S           S
  2131.                                ...
  2132.  
  2133.  
  2134. considering this reverse mapping equations would look like:
  2135.  
  2136.  
  2137.                 ((i-HW_SCREEN_X_CENTRE)/focus) * (V'z* ...
  2138.        u'= S * --------------------------------------------
  2139.                                 ...
  2140.  
  2141.                                 ...
  2142.  
  2143. This would allow us to scale texture on the polygon by changing S
  2144. multiplier.
  2145.  
  2146. Another problem has to do with point O - mapping of the origin of the
  2147. texture space into the world space. The thing is, if the mapping of this
  2148. point gets onto Z==0 plane mapping equations above no longer hold. Just 
  2149. due to the fact that perspective transformation is having singularity 
  2150. there. In general we deal with this problem of perspective transformation
  2151. by clipping the polygon against some plane in front of Z==0. And do we
  2152. really need mapping of exactly O of the texture space? Not necessarily,
  2153. it may be any point in texture space assuming we change a bit direct
  2154. mapping equations:
  2155.  
  2156.  
  2157.                       x' = Ox + (v'-v'0)*V'x + (u-u0)*U'x
  2158.                                    ...
  2159.  
  2160.  
  2161. where (v'0,u'0) are coordinates in the texture space of the point we
  2162. have a mapping for in the world space that is (Ox,Oy,Oz). And the
  2163. reverse mapping equations would be then:
  2164.  
  2165.                    ((i-HW_SCREEN_X_CENTRE)/focus) * (V'z* ...
  2166.            u'= S * ------------------------------------------- + u'0
  2167.                                    ...
  2168.  
  2169.                                    ...
  2170.  
  2171. How can we obtain a mapping of, now, any point from the texture
  2172. space into the world space? If we would associate with every
  2173. polygon's vertex besides just x,y,z also u,v of where this vertex
  2174. is in the texture, we can treat that as 5-D polygon and perform
  2175. clipping on it, obtaining vertexes with coordinates suitable
  2176. for perspective transformation, hence any of them would be fine 
  2177. for the mapping equations. (How to do clipping on polygons is 
  2178. covered in the next chapter).
  2179.  
  2180. Last thing, in order to render regular patterns, those repeating
  2181. themselves, we may want to introduce texture size parameter,
  2182. and make sure that if we want to access texture outside this
  2183. size, this access would wrap around to pick a colour from somewhere
  2184. inside the texture. This can be easily done if texture size
  2185. is chosen to be some power of 2, for in this case we can
  2186. perform range checking with just logical "and" operation.
  2187.  
  2188. However, using above equations for texture mapping appear to be
  2189. quite expensive, there are few divisions involved per each rendered
  2190. pixel. How to optimize that? The, perhaps, easiest way is: horizontal
  2191. line subdivision, that is calculating real texture mapping every N pixels
  2192. and linearly interpolating in between, this method is used in the
  2193. implementation below:
  2194.  
  2195.  
  2196. void G_textured_polygon(int *edges,int length,int *O,int *u,int *v,
  2197.             unsigned char *texture,int log_texture_size,
  2198.                            int log_texture_scale
  2199.                )
  2200. {
  2201.  int new_edges[G_MAX_SHADED_TUPLES];        /* verticaly clipped polygon */
  2202.  int new_length;
  2203.  signed_32_bit Vx,Vy,Vz;
  2204.  signed_32_bit Ux,Uy,Uz;                    /* extracting vectors */
  2205.  signed_32_bit x0,y0,z0;
  2206.  signed_32_bit ui,uj,uc;
  2207.  signed_32_bit vi,vj,vc;
  2208.  signed_32_bit zi,zj,zc;                    /* back to texture coeficients */
  2209.  signed_32_bit v0,u0;
  2210.  signed_32_bit xx,yy,zz,zzz;
  2211.  int xstart,xend;
  2212.  signed_32_bit txstart,tystart;
  2213.  signed_32_bit txend,tyend;
  2214.  signed_32_bit txinc,tyinc;
  2215.  register unsigned char *g_adr,*adr;
  2216.  register int i,x,y;
  2217.  signed_32_bit txtrmasc=(0x1<<(log_texture_size+G_P))-0x1;
  2218.  log_texture_scale+=G_P;
  2219.  
  2220.  GI_boarder_array_init();                   /* initializing the array */
  2221.  
  2222.  new_length=C_polygon_x_clipping(edges,new_edges,4,length);
  2223.  
  2224.  for(i=0;i<new_length;i++)
  2225.   GI_scan(&new_edges[i*4],2,4);             /* Searching polygon boarders */
  2226.  
  2227.  if(G_miny<=G_maxy)                         /* nothing to do? */
  2228.  {
  2229.   x0=O[0]; y0=O[1]; z0=O[2];
  2230.   u0=O[3]<<G_P; v0=O[4]<<G_P;               /* world point <-> texture point */
  2231.  
  2232.   Vx=v[0]; Vy=v[1]; Vz=v[2];
  2233.   Ux=u[0]; Uy=u[1]; Uz=u[2];                /* extracting vectors */
  2234.  
  2235.   ui=(Vz*y0)-(Vy*z0);
  2236.   uj=(Vx*z0)-(Vz*x0);
  2237.   uc=(Vy*x0)-(Vx*y0);
  2238.   vi=(Uy*z0)-(Uz*y0);
  2239.   vj=(Uz*x0)-(Ux*z0);
  2240.   vc=(Ux*y0)-(Uy*x0);
  2241.   zi=(Vy*Uz)-(Vz*Uy);
  2242.   zj=(Vz*Ux)-(Vx*Uz);
  2243.   zc=(Vx*Uy)-(Vy*Ux);                       /* back to texture */
  2244.  
  2245.   g_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  2246.  
  2247.   for(; G_miny<=G_maxy; G_miny++,
  2248.       g_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  2249.   {
  2250.    xstart=G_start[G_miny];
  2251.    adr=g_adr+xstart;
  2252.    xstart-=HW_SCREEN_X_CENTRE;
  2253.    x=xstart;
  2254.    xend=G_end[G_miny]-HW_SCREEN_X_CENTRE;
  2255.    y=G_miny-HW_SCREEN_Y_CENTRE;
  2256.  
  2257.    xx=((y*uj)>>T_LOG_FOCUS)+uc;
  2258.    yy=((y*vj)>>T_LOG_FOCUS)+vc;
  2259.    zz=((y*zj)>>T_LOG_FOCUS)+zc;             /* valid for the hole run */
  2260.  
  2261.    if(( zzz=zz+((x*zi)>>T_LOG_FOCUS) )!=0)
  2262.    {
  2263.     txend=( (xx+ ((x*ui)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +u0;
  2264.     tyend=( (yy+ ((x*vi)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +v0;
  2265.    }
  2266.  
  2267.    for(;xstart<xend;)
  2268.    {
  2269.     x+=G_LINEAR; if(x>xend) x=xend;         /* size of linear run */
  2270.     txstart=txend;
  2271.     tystart=tyend;
  2272.  
  2273.     if(( zzz=zz+((x*zi)>>T_LOG_FOCUS) )!=0)
  2274.     {
  2275.      txend=( (xx+ ((x*ui)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +u0;
  2276.      tyend=( (yy+ ((x*vi)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +v0;
  2277.     }
  2278.  
  2279.     length=x-xstart;                        /* ends here */
  2280.     if(length!=0) { txinc=(txend-txstart)/length;
  2281.             tyinc=(tyend-tystart)/length;
  2282.           }
  2283.  
  2284.     for(;xstart<x;xstart++)                 /* linear run */
  2285.     {
  2286.      txstart&=txtrmasc; tystart&=txtrmasc;  /* wrap around */
  2287.      *adr++=*(texture+((tystart>>G_P)<<log_texture_size)+(txstart>>G_P));
  2288.      txstart+=txinc; tystart+=tyinc;        /* new position in the texture */
  2289.     }
  2290.    }
  2291.   }
  2292.  }
  2293. }
  2294.  
  2295.  
  2296.  
  2297.  
  2298. Other possible speed-up techniques are: area subdivision involving
  2299. 2-D interpolation, faking texture mapping with polynomials (later
  2300. has not very much to do with the true mapping equations described
  2301. here, however) and rendering texture along constant z lines (and
  2302. not along horizontal line) the advantage gained by the former is
  2303. that along every constant Z line perspective transformation is
  2304. linear ( perspective transformation is i=x/z if z==const we
  2305. can write it as i=coef*x where coef=1/z which is linear of course)
  2306.  
  2307. The function above might be extended with interpolative shading
  2308. as well. To do that special consideration has to be given to
  2309. palette's layout, that is where and how to keep different intensities
  2310. of the same colour. But once decided coding is quite straightforward.
  2311.  
  2312.                                    * * *
  2313.  
  2314.  
  2315.  
  2316.  
  2317.  
  2318.                                2-D clipping.
  2319.                               ---------------
  2320.  
  2321.  
  2322.  
  2323.  
  2324.  
  2325.  
  2326.                                          "Don't discount flying pigs before
  2327.                                          you have good air defense."
  2328.                                          -- jvh@clinet.fi
  2329.  
  2330.  
  2331.  
  2332.  
  2333.  
  2334.  
  2335. Screen boundaries clipping.
  2336. ---------------------------
  2337. Throughout the last chapter we've assumed that primitives we were rendering
  2338. are completely within the screen boundaries. This is something we can't
  2339. really guarantee. The worst thing- if we would try using those functions
  2340. to render an element outside the screen, they won't much complain and
  2341. would most likely try accessing memory outside the colourmap's allocated
  2342. space. And I guess: segmentation fault, core dumped is hardly most
  2343. graceful way to exit a program. So we don't have another choice but
  2344. to guarantee rendering algorithms that the element passed is indeed
  2345. inside the screen.
  2346.  
  2347.  
  2348.  
  2349. A dot.
  2350. ------
  2351. Clipping looks trivial in case of a dot, we just have to add few
  2352. comparisons checking if the coordinates are in the screen range and
  2353. proceed only if that is the case:
  2354.  
  2355. void G_dot(int x,int y,unsigned char colour)
  2356. {
  2357.  if( (x>=0)&&(x<HW_SCREEN_X_SIZE) && (y>=0)&&(y<HW_SCREEN_Y_SIZE) )
  2358.  {
  2359.   G_buffer[y*HW_SCREEN_X_SIZE+x]=colour;
  2360.  }
  2361. }
  2362.  
  2363. And of course there if we are always clipping to a rectangular with a
  2364. diagonal (0,0-something_x,something_y) we can do it with just 2
  2365. comparisons instead of 4, just by making sure x and y passed are
  2366. considered by unsigned comparisons, (negative numbers after all 
  2367. would be thought of as just very big positive).
  2368.  
  2369.  
  2370. A line.
  2371. -------
  2372. We can extend above method for line and check every time before
  2373. plotting a dot whether it is within the boundaries. Unfortunately by
  2374. doing that we push up the complexity of the code within the loop.
  2375. And moreover our optimized line routine work with addresses of
  2376. points within the colourmap rather then with the coordinates. What
  2377. we would have to construct is a function that would take in
  2378. arbitrary coordinates of some line and return back the coordinates
  2379. of a line's part which is within the screen boundaries. This
  2380. process is called clipping. The clipping routine must basically
  2381. find coordinates of an intersection point between the line and
  2382. some of the screen edges.
  2383.  
  2384.                        * A
  2385.                         \\
  2386.                          \ \
  2387.                           \  \I1
  2388.                            \ +-o---------------- Y=0
  2389.                             \|   \
  2390.                            I2o     \
  2391.                              |\      * B
  2392.                              | *C
  2393.                              |
  2394.                             X=0
  2395.  
  2396.                           pic 4.1
  2397.  
  2398. The problem is that we are not sure on what edge the intersection
  2399. point lies, for example in the picture above line A-B intersects
  2400. screen at the edge Y==0, whereas line A-C intersects it at the
  2401. edge X==0. Generally we can avoid thinking where the intersection
  2402. is located by consecutive clipping a line against first, say,
  2403. vertical boundaries limits and then horizontal:
  2404.  
  2405.                         * A  |
  2406.                           \  |
  2407.                             \oI1
  2408.                              |\
  2409.                        ------+--o--------------- Y=0
  2410.           * G                | I2 \
  2411.            \                 |      \
  2412.              \               |        * B     *E
  2413.               *H        C*---o--*D             \
  2414.                              |                  *F
  2415.                              |
  2416.                             X=0
  2417.  
  2418.                           pic 4.2
  2419.  
  2420. In the example above after first call to clipping routine line
  2421. A-B would become I1-B, after the second it would assume the
  2422. final acceptable form I2-B. It can be seen also that some lines
  2423. going through this clipping method won't actually be clipped at
  2424. all (E-F) some would be clipped just once (C-D), some twice
  2425. (A-B) and finally some would be completely clipped, that is,
  2426. would be totally outside the screen area (G-H). Now, how to find
  2427. an intersection point? obviously it is not a very complicated
  2428. case of solving similar triangles:
  2429.  
  2430.  
  2431.                           |
  2432.               (x1,y1) *   |
  2433.                       | \ |
  2434.                       |   * (xm,ym)
  2435.                       |   | \
  2436.                       +---|---* (x2,y2)
  2437.                           |
  2438.  
  2439.                         pic 4.3
  2440.  
  2441.  
  2442.         y2-ym     y2-y1                    (y2-y1)*(x2-xm)
  2443.        ------- = -------   =>   ym = y2 - -----------------
  2444.         x2-xm     x2-x1                         x2-x1
  2445.  
  2446.  
  2447. but this method involves multiplications and divisions, so
  2448. beholding to our tradition let's try to avoid it. The alternative
  2449. method is using binary search:
  2450.  
  2451.  
  2452.                       A(-3,1)
  2453.                         *    |
  2454.                          \   |
  2455.                      (-1,3)* |
  2456.                              o I(0,?)
  2457.                              | \
  2458.                              |   * B(1,5)
  2459.                              |
  2460.                             X=0
  2461.  
  2462.                           pic 4.4
  2463.  
  2464.  
  2465. Let's see how it works on a simple example: suppose we have
  2466. to clip a line A(-3,1)-B(1,5) by an edge X=0 , we have to
  2467. find Y of the intersection point. let's find midpoint of
  2468. the line A-B:
  2469.  
  2470.                    Ax+Bx              Ay+By
  2471.              midx=-------       midy=-------
  2472.                      2                  2
  2473.  
  2474.              midx=(-3+1)/2=-1   midy=(1+5)/2=3
  2475.  
  2476. Now, let's see where the boundary lies with respect to the midpoint?
  2477. It is still to the right of it, so of two lines A-MidPoint MidPoint-B,
  2478. edge intersects MidPoint-B. Let's rename MidPoint into A and start
  2479. the midpoint search over again:
  2480.  
  2481.              midx=(-1+1)/2=0    midy=(3+5)/2=4
  2482.  
  2483. Here, the intersection came precisely onto the edge. So midy is the
  2484. Y coordinate of Intersection point we were looking for. It appears
  2485. that this method involve a lot of calculations, but they are all
  2486. very cheap, just an addition and division by two (which is
  2487. actually a 1 bit right shift). Practice shows binary search works
  2488. indeed faster then calculation resulting from solving similar
  2489. triangles.
  2490.  
  2491. When dealing with divisions performed by shits however, one has to be 
  2492. aware of truncation problems that might arise. Since -1>>1=-1 that means 
  2493. that if we would try to use algorithm described above to clipp (-1,y1)-(0,y2) 
  2494. line by X==0 boundary chances are we would loop forever, at each iteration 
  2495. finding that -1>>1 is still -1. (and O(for ever) is hardly, the time 
  2496. complexity contributing towards high frame-rate). As in the code below 
  2497. this situation should be thought of.
  2498.  
  2499. Let's create a function which would perform clipping against
  2500. two vertical screen edges: that is C_X_CLIPPING_MIN and C_X_CLIPPING_MAX.
  2501.  
  2502.  
  2503.  
  2504. int C_line_x_clipping(int **vertex1,int **vertex2,int dimension)
  2505. {
  2506.  register int i;
  2507.  register int whereto;
  2508.  register int *l,*r,*m,*t;                  /* left right and midle and tmp */
  2509.  static int g_store0[C_MAX_DIMENSIONS];     /* static stores for clipped vxes */
  2510.  static int g_store1[C_MAX_DIMENSIONS];
  2511.  static int g_store2[C_MAX_DIMENSIONS];
  2512.  static int g_store3[C_MAX_DIMENSIONS];
  2513.  static int g_store4[C_MAX_DIMENSIONS];
  2514.  static int g_store5[C_MAX_DIMENSIONS];
  2515.  int **vmn,**vmx;                           /* so that *vmn[0] < *vmx[0] */
  2516.  int swap;                                  /* we coordinates swaped? */
  2517.  
  2518.  C_2D_clipping=0;                           /* default no clipping yet */
  2519.  
  2520.  if((*vertex1)[0]<(*vertex2)[0])
  2521.  { swap=0; vmn=vertex1; vmx=vertex2; }      /* so that *vmn[0] < *vmx[0] */
  2522.  else
  2523.  { swap=1; vmn=vertex2; vmx=vertex1; }
  2524.  
  2525.  if(((*vmn)[0]>C_X_CLIPPING_MAX)||((*vmx)[0]<C_X_CLIPPING_MIN)) return(0);
  2526.  else
  2527.  {
  2528.   if((*vmn)[0]<=C_X_CLIPPING_MIN)           /* clipping */
  2529.   {
  2530.    HW_copy_int(*vmn,l=g_store0,dimension);  /* copying old vertixes */
  2531.    HW_copy_int(*vmx,m=g_store1,dimension);
  2532.    r=g_store2;                              /* let middle be there */
  2533.  
  2534.    whereto=0;
  2535.    while(m[0]!=C_X_CLIPPING_MIN)
  2536.    {
  2537.     if(whereto==1) { t=l; l=m; m=t; }
  2538.     else           { t=r; r=m; m=t; }
  2539.     for(i=0;i<dimension;i++) m[i]=(l[i]+r[i])>>1;
  2540.     whereto=m[0]<C_X_CLIPPING_MIN;
  2541.    }
  2542.    *vmn=m;                                  /* that is why m[] is static */
  2543.    C_2D_clipping=swap^1;
  2544.   }
  2545.  
  2546.   if((*vmx)[0]>=C_X_CLIPPING_MAX)           /* clipping */
  2547.   {
  2548.    HW_copy_int(*vmn,l=g_store3,dimension);  /* copying old vertixes */
  2549.    HW_copy_int(*vmx,m=g_store4,dimension);
  2550.    r=g_store5;                              /* let middle be here */
  2551.  
  2552.    whereto=0;
  2553.    while(m[0]!=C_X_CLIPPING_MAX)
  2554.    {
  2555.     if(whereto==1) { t=l; l=m; m=t; }
  2556.     else           { t=r; r=m; m=t; }
  2557.     for(i=0;i<dimension;i++) m[i]=(l[i]+r[i])>>1;
  2558.     whereto=m[0]<C_X_CLIPPING_MAX;
  2559.    }
  2560.    *vmx=m;                                  /* that is why m[] is static */
  2561.    C_2D_clipping|=swap&1;
  2562.   }
  2563.  }
  2564.  return(1);                                 /* partialy or not clipped */
  2565. }
  2566.  
  2567.  
  2568.  
  2569. The return value of this function is zero if line is completely
  2570. clipped otherwise it is 1. The purpose of global C_2D_clipping variable
  2571. would become clear when we would be talking about polygon clipping.
  2572. As to the Y clipping we can either modify the above code to
  2573. do both functions or just duplicate most of the code changing few
  2574. things, like constants names ( C_Y_CLIP... instead of C_X_CLIPP...)
  2575. and indexes of clipped variable ( 1 instead of 0).
  2576.  
  2577. As you can see this code is generic enough to allow clipping of
  2578. N-dimensional lines ( so that to allow for X Y Z Colour Intensity etc. to
  2579. be processed at the same time). To do it effectively static arrays are
  2580. being set to keep temporary left (l) middle (m) and right (r) N-tuples.
  2581. Whenever making decision which segment to take for the next iteration
  2582. we can just swap pointers, so that array keeping coordinate of segment
  2583. being discarded could be set to accept middle point at the next iteration.
  2584.  
  2585.             <-----------  m
  2586.          l  ----------->                   r
  2587.          |                |                |
  2588.          v                v                v
  2589.          +---------+      +---------+      +---------+
  2590.          |x,y,z... |      |x,y,z... |      |x,y,z... |
  2591.          +---------+      +---------+      +---------+
  2592.  
  2593.                             pic 4.5
  2594.  
  2595. This is being done so that at every iteration only pointers to
  2596. actual data are moved not the data itself. (The point I am trying to
  2597. make: (and I am sure everybody knows that, just a bit of reinforcement)
  2598. it's ok to physically copy small amounts of data, but when we have a lot
  2599. of it, it is better to move pointers instead)
  2600.  
  2601.  
  2602. Let's insert calls to clipping function into the G_line routine:
  2603.  
  2604.  
  2605.                                   ...
  2606.                                   ...
  2607.  
  2608.  v1=vertex1;
  2609.  v2=vertex2;
  2610.  
  2611.  if(C_line_x_clipping(&v1,&v2,2))           /* horizontal clipping */
  2612.  {
  2613.   if(C_line_y_clipping(&v1,&v2,2))          /* vertical clipping */
  2614.   {
  2615.    dx=v2[0]-v1[0]; dy=v2[1]-v1[1];          /* ranges */
  2616.  
  2617.                                   ...
  2618.                                   ...
  2619.  
  2620.  
  2621. So, whenever line is completely clipped we wouldn't go any further
  2622. in the rasterization function.
  2623.  
  2624.  
  2625.  
  2626. A polygon.
  2627. ----------
  2628. Now, remembering that we render polygon by scanning it's edges
  2629. using rasterization function pretty much like the one above, we
  2630. may think that adding clipping calls into GI_scan would solve our
  2631. problem with polygon clipping, unfortunately, it is only half true
  2632. (literary half true). Lets consider pictures pic 4.6 and 4.7:
  2633.  
  2634.  
  2635.                    B                          * B
  2636.                  | *====                    I/
  2637.                  |/====                -----*------
  2638.                 I*======                   /======
  2639.                 /|??????                  /========
  2640.                / |?????                  /========
  2641.              A*  |??????               A*==========
  2642.  
  2643.  
  2644.                 pic 4.6                   pic 4.7
  2645.  
  2646. In the cases above edge A-B contribute to the left side of the
  2647. polygon, clipping present no problem for case on pic 4.7. Clipped
  2648. part I-B of the edge can be discarded since there's nothing to
  2649. form outside the screen. On the other hand in the case on
  2650. pic 4.6 we can't simply discard clipped part A-I since we would
  2651. loose left boundary of all horizontal lines marked "???". The
  2652. observation we can make is that polygon, being scanned edge at
  2653. a time, may not be clipped against horizontal boundaries but
  2654. must be clipped against vertical, this way the code within GI_scan
  2655. would look like:
  2656.  
  2657.                                   ...
  2658.                                   ...
  2659.  
  2660.  v1=edges; edges+=skip; v2=edges;           /* length ints in each */
  2661.  
  2662.  if(C_line_y_clipping(&v1,&v2,dimension))   /* vertical clipping */
  2663.  {
  2664.   dx=*v2++; dy=*v2++;                       /* extracting 2-D coordinates */
  2665.   x=*v1++; y=*v1++;                         /* v2/v1 point remaining dim-2 */
  2666.  
  2667.                                   ...
  2668.                                   ...
  2669.  
  2670.  
  2671. Creating a vertically clipped polygon on the other hand is a bit
  2672. more complicated. The problem is that the clipped polygon may
  2673. have different from the original one number of edges. let's
  2674. consider the situation on the pic 4.8:
  2675.  
  2676.                                                /* 3
  2677.                      |                      | / |
  2678.                      |5'                    |/  |
  2679.                   5*-*-----*6               *   |
  2680.                   /  |      \              /|2''|
  2681.                 4*   |       *1         2 * |   |
  2682.                   \  |      /              \|   |
  2683.                   3*-*-----*2               *\  |
  2684.                      |3'                    |2'\* 1
  2685.  
  2686.  
  2687.                    pic 4.8                pic 4.9
  2688.  
  2689. It can be seen that original polygon has 6 edges which becomes 5
  2690. after clipping. Some other polygon pic 4.9 can have 1 more edge
  2691. after clipping. It is obvious that we would have to create a
  2692. temporary structure to hold the clipped polygon. Now, we know how
  2693. to clip a single edge, let's try to find pattern how clipped
  2694. edges compose clipped polygon. Let's be considering clipping
  2695. against a single boundary, say X==0. it is obvious that if an edge
  2696. is completely outside it's vertexes won't be present in the new
  2697. polygon. On the other hand if an edge is within the boundary
  2698. both points would be present. Any point on the intersection
  2699. line would be present too. One more consideration is that each
  2700. point in the polygon belong to two edges, so when the edge is
  2701. not clipped we may put into the new structure just the second
  2702. point assuming that the first one is already taken care of when
  2703. we were processing previous edge.
  2704.  
  2705. Let's list the patterns:
  2706.  
  2707. (1) If edge is not clipped or the second point is clipped put into
  2708.     new structure just the second point;
  2709.  
  2710. (2) When first point is changed or both are changed put both;
  2711.  
  2712. (3) When both are outside put none.
  2713.  
  2714. Surprisingly enough this algorithm doesn't have to be changed
  2715. when we consider simultaneous clipping against two parallel
  2716. boundaries:
  2717.  
  2718.  
  2719.                              | *-----*1|
  2720.                              |/5      \|
  2721.                            4'*         *2'
  2722.                             /|         |\
  2723.                           4* |         | \
  2724.                            | |         |  \
  2725.                            *-*---------*---*2
  2726.                            3 |3'       |2''
  2727.                              |         |
  2728.  
  2729.                                pic 4.10
  2730.  
  2731. edge 1-2   Second point clipped - pattern (1)  putting in 2'
  2732. edge 2-3   Both points changed  - pattern (2)  putting in 2'' and 3'
  2733. edge 3-4   Both points outside  - pattern (3)  ignoring
  2734. edge 4-5   First point clipped  - pattern (2)  putting in 4'  and 5
  2735. edge 5-1   No clipping          - pattern (1)  putting in 1
  2736.  
  2737. The result is: 2'-2''-3'-4'-5-1 just looking at the picture assures
  2738. us that what we got is indeed right. Now finally the purpose of
  2739. C_2D_clipping variable being set in C_line_x_clipping must be clear.
  2740. It would be set to 1 whenever first or both points were changed,
  2741. and would be 0 if none or just second one were changed. Knowing
  2742. this all let's write some code:
  2743.  
  2744.  
  2745.  
  2746. int C_polygon_x_clipping(register int *from,register int *to,
  2747.                          int dimension,int length
  2748.                         )
  2749. {
  2750.  register int i;
  2751.  int *v1,*v2,new_lng=0;
  2752.  int *first_vrtx=to;                        /* begining of the source */
  2753.  
  2754.  for(i=0;i<length;i++)                      /* for all edges */
  2755.  {
  2756.   v1=from; from+=dimension; v2=from;        /* taking two vertexes */
  2757.  
  2758.   if(C_line_x_clipping(&v1,&v2,dimension))  /* clipping */
  2759.   {
  2760.    if(C_2D_clipping)                        /* depends which one was clipped */
  2761.    {
  2762.     HW_copy_int(v1,to,dimension); to+=dimension;
  2763.     HW_copy_int(v2,to,dimension); to+=dimension;
  2764.     new_lng+=2;                             /* first point clipped */
  2765.    }
  2766.    else
  2767.    {
  2768.     HW_copy_int(v2,to,dimension); to+=dimension;
  2769.     new_lng++;                              /* second point clipped */
  2770.    }
  2771.   }
  2772.  }
  2773.  HW_copy_int(first_vrtx,to,dimension);      /* looping the polygon vertexes */
  2774.  
  2775.  return(new_lng);
  2776. }
  2777.  
  2778.  
  2779.  
  2780.  
  2781. Again the way we represent the polygon the very first point has to
  2782. be the last too since every point belong to two edges, and we can
  2783. have a pointer to a new edge by just advancing it by "dimension" of
  2784. a single vertex in the list of polygon vertices.
  2785.  
  2786.  
  2787. View plane clipping.
  2788. --------------------
  2789. As discussed before if we want to do perspective transformation we
  2790. have to make sure we are actually applying it to something that we know
  2791. produces a valid result, hence there should be no points with negative
  2792. or zero Z coordinates. Zero would produce divide errors, Negative once
  2793. would appear flipped over. There is another good reason to clip out 
  2794. vertices with negative Z, we just can't see things which are behind 
  2795. us (at least majority of us can't).
  2796.  
  2797. The technique we would be using is just an expansion of 2-D clipping
  2798. that was used to make sure 2-D primitives fit physical screen before
  2799. rendering. As before, we can find the intersection point from solving
  2800. similar triangles, or using binary-search techniques. The single clipping
  2801. plane would be specified to be somewhere in front of Z==0. (And no, not
  2802. quite as far as "focus" distance used in perspective transformation,
  2803. viewing plane is not necessarily a clipping plane).
  2804.  
  2805. The function to clip polygons would be almost identical too. (By
  2806. the way, we can indeed try writing generic line and polygon clipping
  2807. functions performing it for both 2D screen boundaries and view plane).
  2808.  
  2809.  
  2810.  
  2811. Volume clipping.
  2812. ----------------
  2813. First of all what view volume is? This is the set of all points in
  2814. the "world" space which can appear on screen under certain projection 
  2815. method. For example if we are using just the parallel projection 
  2816. (let's forget for a second about perspective) we may see that only 
  2817. points from depicted below rectangular area would appear on the screen.
  2818.  
  2819.  
  2820.                         |                |
  2821.                         |                |
  2822.                         |                |    * B
  2823.                         |        * C     |
  2824.              Screen     |        |       |
  2825.              boundaries |        v       |
  2826.                    -----I--------*-------I----- Viewing plane.
  2827.  
  2828.                             * A
  2829.  
  2830.  
  2831.                               pic 4.11
  2832.  
  2833. Point A can't appear on screen as it is behind the viewing plane - it would
  2834. be clipped by the algorithm described above in this chapter. point B on the
  2835. other hand is in front of the viewing plane, so it would pass the first test
  2836. but since it's projection onto the viewing plane won't fit the screen it would
  2837. be clipped by 2-D screen boundaries clipping routines. Point C is lucky to be
  2838. inside the view volume. The conclusion here is that by doing first viewing
  2839. plane clipping and then screen boundaries clipping we actually performed
  2840. volume clipping, that is we selected from space the set of all points that
  2841. can be projected and then did actual projection (Yes, by passing to the
  2842. rendering routines just X and Y of points and discarding Z we basically did
  2843. parallel projection). Should this be changed somehow to accommodate
  2844. perspective projection? after all by clipping against view plane we guarantee
  2845. that there would be no points with Z==0. However the problem is that even
  2846. having points close to the viewing plane is not very good. pic 4.12
  2847.  
  2848.  
  2849.  
  2850.  
  2851.                                         *--->
  2852.  
  2853.                                         *------------->
  2854.                                         *-------------------------->
  2855.                                                                      A*- - ->
  2856.                           -----I-----I-----
  2857.  
  2858.  
  2859.  
  2860.                                pic 4.12
  2861.  
  2862. By applying perspective transformation we increase the absolute values
  2863. of coordinate components depending on inverse of it's distance to the viewer,
  2864. so if Z==0 transforms into infinity numbers close to that transform into
  2865. something bitsize of numbers stored in computers might not be able to handle,
  2866. and no, we can't quite solve it by moving the clipping plane slightly forward,
  2867. since there are still those nasty points which are slightly in front of the
  2868. viewing plane but already have big absolute value of X or Y. (pic 4.12,
  2869. point A). The overflow problem that may result from perspectively transforming
  2870. this point is this: positive values may become negative, and if it would
  2871. happen to just one point of two off-screen points representing a line we
  2872. may actually see this line all across the screen instead of not seeing it
  2873. at all.
  2874.  
  2875. The conclusion when using perspective transformation, it is better to apply
  2876. it to the points we know would produce a valid result. And since what we
  2877. ultimately want is to project to the screen we are coming back to the
  2878. view volumes since those are exactly sets of points that would be
  2879. projected inside the screen. Hence we do need to consider actual viewing
  2880. volume for perspective projection.
  2881.  
  2882. What the view volume for perspective projection would look like?
  2883. The way we modeled this transformation - rays of all visible points
  2884. converge in viewer's eye. If we would cast back into space lines connecting
  2885. the eye and the screen boundaries, we would obtain the pyramid marking points
  2886. from space that would project onto screen edges. what's inside the
  2887. pyramid would project inside the screen what's outside - outside the screen.
  2888. adding the clipping plane somewhere in front of the viewer we are obtaining the
  2889. view-volume for perspective projection, pic 4.13.
  2890.  
  2891.  
  2892.                            \             /
  2893.                             \           /
  2894.                              \         /
  2895.                               \       /
  2896.                                \     /
  2897.                           -----I\---/I-----
  2898.                                  \ /
  2899.                                   *
  2900.  
  2901.                                pic 4.13
  2902.  
  2903. The only problem - we now have to clip against planes which are
  2904. directed almost arbitrary in space (the exact geometry of this view volume
  2905. depends on the "focus" parameter - the distance between the viewer and the
  2906. viewing plane (again, not quite the same as clipping plane). Although 
  2907. conceptually easy to achieve clipping against arbitrary directed plane 
  2908. would have higher complexity.
  2909.  
  2910. There is number of solutions one can consider: obvious one: indeed implement
  2911. real volume clipping, although expansive we would be able to completely
  2912. get rid of 2-D clipping routines which overall might be quite descent.
  2913. Second: implement volume clipping with special kind of perspective view
  2914. volume, the one having 90' angle. The reason can be seen from the pic 4.14:
  2915.  
  2916.  
  2917.                        \                     /
  2918.                          \                 /
  2919.                     x=-z   \             /   x=z
  2920.                              \         /
  2921.                           ----I\-----/I-----
  2922.                                  \ /
  2923.                                   *
  2924.  
  2925.                                pic 4.14
  2926.  
  2927. Planes composing this perspective view volume have pretty simple equations:
  2928. x==z , y==z, x==-z and y==-z clipping against those is way easier then
  2929. against arbitrary directed planes. The method however we would discuss is a
  2930. bit different, that is, we would not do clipping at all, to be more exact we
  2931. would only throw away polygons which are for sure outside the view volume and
  2932. leave those even partially inside to be further clipped against screen
  2933. boundaries after being projected. ( Clipping against viewing plane, is a
  2934. sacred matter that one can hardly get rid of ) One should understand that this
  2935. method works on assumption that there is no terribly big polygons around,
  2936. since a one part of a really big polygon can be within the view volume whereas
  2937. another can be both close to viewing plane and have really big absolute value
  2938. of X or Y, just something we are trying to exclude from being perspectively
  2939. projected.
  2940.  
  2941. How can one figure out whether a point is outside the pyramid with 90'
  2942. angle? From the pic 4.14 we can see that parts of the space separated
  2943. by Z==X and Z==-X planes have obvious relationships among X and Z of
  2944. every point, if Z<X or Z<-X point is for sure outside this area. it
  2945. should be realized we can't quite extend this scheme onto polygons by 
  2946. saying if all points of a polygon are outside it is outside also, 
  2947. example is a polygon AB on the pic 4.15, it is clear that although 
  2948. both points composing it are outside there supposed to be some part
  2949. of this polygon inside the view volume.
  2950.  
  2951.  
  2952.                                   ^ Z
  2953.                                   |
  2954.                        \          |          /
  2955.                          \   Z>-X | Z>X    /
  2956.                            \      |      /
  2957.                      A *---------------------* B
  2958.                                \  |  /
  2959.                           Z<-X   \ /   Z<X
  2960.                     --------------*-----------------> X
  2961.  
  2962.                                pic 4.15
  2963.  
  2964.  
  2965. We would make decisions, on the other hand, using notion of polygon's 
  2966. extends. pic 4.16 Extend is a cube completely enclosing within itself 
  2967. certain polygon. So coordinates of extend planes are that of maximum and 
  2968. minimum among the polygon vertices along all axes.
  2969.  
  2970.                             xmin       xmax
  2971.                             --------------  ymin
  2972.                              |\\        |
  2973.                              | \  \     |
  2974.                              |   \    \ |
  2975.                              |     \   /|
  2976.                              |       \/ |
  2977.                              |------------ ymax
  2978.  
  2979.                                  pic 4.16
  2980.  
  2981. This way by considering for example (xmin,zmax) point of the extend box we
  2982. can make a decision whether polygon's extend is outside the x=z plane.
  2983.  
  2984.  
  2985.  
  2986.                                   ^ Z
  2987.                                   |
  2988.                        \          |          /
  2989.                          \        |        /
  2990.               (xmax,zmax)  \      |      /  (xmin,zmax)
  2991.                      ----+   \    |    /  +----+
  2992.                          |     \  |  /    |
  2993.                        --+       \ /      +---
  2994.                     --------------*-----------------> X
  2995.                                      +-------+ (zmax)
  2996.                                      |       |
  2997.  
  2998.  
  2999.                                pic 4.17
  3000.  
  3001.  
  3002. Let's list all the other cases:
  3003.  
  3004.                xmin > zmax
  3005.                ymin > zmax
  3006.               -xmax > zmax
  3007.               -ymax > zmax
  3008.  
  3009. On the same stage we can check if the polygon is completely behind
  3010. the view plane as well.
  3011.  
  3012.  
  3013. int C_volume_clipping(register int *vxes,int number)
  3014. {
  3015.  int xmin,ymin,zmin,xmax,ymax,zmax;
  3016.  int i;
  3017.  
  3018.  ymin=xmin=zmin=INT_MAX;
  3019.  ymax=xmax=zmax=INT_MIN;                    /* initializing searching */
  3020.  
  3021.  for(i=0;i<number;i++)
  3022.  {
  3023.   if(*vxes>xmax) xmax=*vxes;
  3024.   if(*vxes<xmin) xmin=*vxes;
  3025.   vxes++;
  3026.   if(*vxes>ymax) ymax=*vxes;
  3027.   if(*vxes<ymin) ymin=*vxes;
  3028.   vxes++;
  3029.   if(*vxes>zmax) zmax=*vxes;
  3030.   if(*vxes<zmin) zmin=*vxes;                /* searching max/min */
  3031.   vxes++;
  3032.  }
  3033.  
  3034.  if((zmax<xmin)||(zmax<ymin)||(zmax<-xmax)||
  3035.    (zmax<-ymax)||(zmax<C_plane))
  3036.   return(0);
  3037.  else
  3038.   if(zmin<C_plane) return(-1);              /* behind clipping plane */
  3039.   else return(1);
  3040. }
  3041.  
  3042.  
  3043. It should be realized that, extend testing is approximate, there
  3044. might be polygons outside the view volume whose extends are partly 
  3045. inside, but since we are doing clipping to screen boundaries 
  3046. afterwards they would eventually be taken care of. Besides, clipping 
  3047. view volume is a pyramid with 90' angle, real view volume on the other 
  3048. hand depends on the perspective transformation formula, and would 
  3049. likely have angle less then 90'. And again, there should be no very 
  3050. big polygons, since we are doing "full element" volume clipping.
  3051.  
  3052.                                      * * *
  3053.