home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / Libraries / 3DGPL 1.0 / TEXT / 2.TXT < prev    next >
Encoding:
Text File  |  1995-06-22  |  43.5 KB  |  1,117 lines  |  [TEXT/ttxt]
  1.  
  2.  
  3.                              3-D transformations.
  4.                             ----------------------
  5.  
  6.  
  7.  
  8.                                     "Several large, artificial constructions
  9.                                     are approaching us," ZORAC
  10.                                     announced after a short pause.
  11.                                     "The designs are not familiar, but they
  12.                                     are obviously the products of
  13.                                     intelligence. Implications: we have been
  14.                                     intercepted deliberately by a means
  15.                                     unknown, for a purpose unknown, and
  16.                                     transferred to a place unknown by
  17.                                     a form of intelligence unknown.
  18.                                     Apart from the unknowns, everything
  19.                                     is obvious."
  20.                                     -- James P. Hogan, "Giants Star"
  21.  
  22.  
  23.  
  24. 3D transformations would assume the central position in the engine.
  25. After all, this is something which would allow us to see objects from
  26. different sides, manage the three dimensional space etc. The
  27. simplest but by no means least important transformations are
  28. translation and scaling. The rotation transformation would be little
  29. bit more complicated to understand and implement but fortunately
  30. not much. Talking of all three transformations we can either consider
  31. them from the point of an object (collection of points in 3D)
  32. being transformed, or a space itself being transformed. In
  33. some cases it is more productive to think one way in some another way.
  34.  
  35.  
  36.  
  37. Translation and scaling.
  38. ------------------------
  39. Translation. This transformation moves points of an object to a new
  40. specified position pic 2.1, On the other hand we may think of it as a
  41. system of references being moved in the opposite direction pic 2.2:
  42.  
  43.  
  44.  
  45.           Y ^                                                 ^
  46.             |          * A'
  47.             |          ^ * B'                                 |
  48.             |          | ^                           Y'^
  49.             |A*--------+ | y component                 |      |A*
  50.             |  B*--------+                             |         B*
  51.             |    x component                           |      |
  52.    ---------+-------------------> X        - - - - - - | - - - - - - -> X
  53.            0|                                          |     0|
  54.             |                               -----------+------------->X'
  55.             |                                        0'|      |
  56.             |                                          |
  57.             |                                          |
  58.             |                                          |
  59.  
  60.          pic 2.1                                    pic 2.2
  61.  
  62. This way if it is a collection of points we are moving, Distances
  63. between them would not change. Any move across 2D space can separate
  64. into 2 components - move along X and move along Y. It would take
  65. 3 components in 3D space - along X, Y and Z. Important place we
  66. may need translation is to , for example , move (0,0) of the screen
  67. from upper left corner of the screen to it's physical centre. So
  68. the translation function might look like:
  69.  
  70.  
  71. void T_translation(register int *from,register int *into,int number,
  72.                    int addx,int addy,int addz
  73.                   )
  74. {
  75.  register int i;
  76.  
  77.  for(i=0;i<number;i++)
  78.  {
  79.   (*into++)=(*from++)+addx;
  80.   (*into++)=(*from++)+addy;
  81.   (*into++)=(*from++)+addz;                 /* translation */
  82.  }
  83. }
  84.  
  85.  
  86. Another transformation is scaling: blowing up of distances between
  87. points. Or we can think of it as a contraction or blowing of the space
  88. itself pic2.3:
  89.  
  90.  
  91.                  A'*      Y ^
  92.                     \       |
  93.                      \      |
  94.                       *   B'*    *C'
  95.                      A      ^   /
  96.                            B*  *C
  97.                             |
  98.               --------------+-----------------> X
  99.                            0|
  100.                             |
  101.                             |
  102.  
  103.                          pic 2.3
  104.  
  105. The obvious usage of this transformation: trying to fit 2000x2000 object
  106. into 200x200 window - evidently every point of this object would have to
  107. have coordinates scaled to fit in the window range. In this case scaling
  108. is just multiplication by 1/10. The object is contracted that way. On the
  109. other hand we can think that it was actually the window space expanding to
  110. enclose the object (we would have to abstract far enough from computer
  111. hardware today on the market to do that, however). In general case
  112. we can introduce separate factors for every axis, if this factor is
  113. greater then 0 and less then 1 we would contract the object if it would
  114. be greater then 1 the object would be expanded. The code is once again
  115. trivial:
  116.  
  117.  
  118. void T_scaling(register int *from,register int *to,int number,
  119.                float sx,float sy,float sz
  120.               )
  121. {
  122.  register int i;
  123.  
  124.  for(i=0;i<number;i++)
  125.  {
  126.   (*to++)=(int)(*from++)*sx;
  127.   (*to++)=(int)(*from++)*sy
  128.   (*to++)=(int)(*from++)*sz;
  129.  }                                          /* scaling */
  130. }
  131.  
  132.  
  133.  
  134. Planar case rotation.
  135. ---------------------
  136. Of those three transformations rotation is by far the most complicated
  137. one, let's first consider planar - 2-D case of rotation. there is some
  138. simple trigonometry involved, sin and cos functions. just for those
  139. who don't deal with this subject every day sin and cos for a triangle
  140. on pic 2.4 are defined as:
  141.  
  142.  
  143.                                   +
  144.                                  /|
  145.                                l/ |
  146.                                /  |y
  147.                               /alp|
  148.                              +----+
  149.                                 x
  150.  
  151.                              pic 2.4
  152.  
  153.  
  154.                             x                  y
  155.                 sin(alp) = ---     cos(alp) = ---
  156.                             l                  l
  157.  
  158. And so if we do know sin(alp) or cos(alp) we can calculate segments x and y
  159. knowing segment l:
  160.  
  161.                  x = l*sin(alp)    y = l*cos(alp)
  162.  
  163. This can be considered as seeking projections of a point for which we
  164. know the distance from the beginning of coordinates onto orthogonal
  165. axes pic 2.5:
  166.  
  167.                                ^
  168.                              Y |
  169.                                +----*
  170.                                |   /|
  171.                               y|  /l|
  172.                                | /  |
  173.                                |/alp|
  174.                          ------+----+---->X
  175.                                |  x
  176.                                |
  177.  
  178.                              pic 2.5
  179.  
  180. Now, let's consider the situation when we rotate the system of references
  181. counterclockwise by the angle alp. What would be coordinates of point
  182. A in the turned system of references Y'X' if it's coordinates in the
  183. original system of references were (x,y)?
  184.  
  185.  
  186.                                ^ Y            X'
  187.                                |             /
  188.                     Y'       Y/*\--*A      /
  189.                       \     /  |  \|     /
  190.                         \ /    |   |\  /
  191.                        Yy'\    |   | /Yx'
  192.                             \  |   |
  193.                               \| /\ Xx'
  194.                   -------------+---+-------------------> X
  195.                               /| \/ X
  196.                             / Xy'  \
  197.                           /    |     \
  198.                         /      |       \
  199.                       /        |         \
  200.                                |
  201.  
  202.                              pic 2.6
  203.  
  204.  
  205. Using the sin/cos formulas from above we can find projections of
  206. x and y (that is of the projections of the point onto original axis's)
  207. onto new axis's X'and Y'. Let's sum projection of x onto X' and projection
  208. of y onto same axis's and do the same thing with projections of x and y
  209. onto Y'. Those sums would be just what we are looking for, coordinates 
  210. of the point A in the new system of references X'Y'
  211.  
  212. i.e.
  213.  
  214.                X' = Yx'+Xx' = y*sin(alp) +    x*cos(alp)
  215.                Y' = Yy'+Xy' = y*cos(alp) +  (-x*sin(alp))
  216.  
  217. Note the sign of Xy', just from looking on the pic2.6 one can see that
  218. in this case x is being projected onto the negative side of Y'. That's
  219. why the negative sign. So those are the transformation formulas for
  220. counterclockwise rotating system of reference by an angle alp:
  221.  
  222.                      x'= y*sin(alp) + x*cos(alp)
  223.                      y'= y*cos(alp) - x*sin(alp)
  224.  
  225. On the other hand we can consider it as point A itself rotating
  226. around zero clockwise.
  227.  
  228.  
  229.  
  230. 3-D rotation.
  231. -------------
  232. This subject is complicated only because it is a bit hard to visualize
  233. what is happening. The idea on the other hand is simple: applying three
  234. 2-D rotations one after another so that each next works with the results
  235. obtained on the previous stage. This way every time we are applying new
  236. transformation we are working not with the original object but the
  237. object already transformed by all previous rotations.
  238.  
  239. Each of 2D rotations this way is bound to some axis around which two
  240. other axes would rotate.
  241.  
  242. There are few important considerations influencing the way we would
  243. find the formulas: 1) What kind of system of references we have. What
  244. is the positive direction of the rotations we would be applying; and
  245. 2) In what order we want to apply 2D rotations.
  246.  
  247.  
  248.  
  249. System of references:
  250. ---------------------
  251. Well, we do have freedom of choosing directions of axes, for one
  252. thing, we are working in orthogonal system of references where each
  253. axes is orthogonal to the plane composed of remaining two. As
  254. to where the positive direction of every axis point we have freedom
  255. to decide. It is customary for example to have positive side of Y
  256. directed upside. We don't actually have to follow customs if we really
  257. don't want to, and remembering that the colourmap, we would be doing all
  258. drawings into, have Y directed down we might want to choose
  259. corresponding direction here. However having it pointing up is more
  260. natural for us humans, so we might save time on debugging a while
  261. afterwards. Let's point it up. As to the Z let's direct it away from
  262. us and X to the right.
  263.  
  264.  
  265.                                  Y^    Z
  266.                                   |   /
  267.                                   |  / alp
  268.                                  /|<---+
  269.                                 | |/   |
  270.                          -------|-+-----------> X
  271.                             bet | |   __
  272.                                 V/|   /| gam
  273.                                 /----+
  274.                                /  |
  275.                               /   |
  276.                                   |
  277.  
  278.                                pic 2.7
  279.  
  280.  
  281. As to the rotations, let's call the angle to turn XY plane around Z axis
  282. as gam, ZY around X as bet and ZX around Y as gam, pic 2.7
  283.  
  284.  
  285.  
  286. Transformation order.
  287. ---------------------
  288. The problem with transformation order is that the point consequently
  289. turned by gam-bet-alp won't be necessarily at the same place with
  290. where it might go if turned alp-bet-gam. The reason is our original
  291. assumption: each next transformation works on already transformed
  292. point by previous 2D rotations. That is, we are rotating coordinates
  293. around moving axes. On the picture 2.7 we can see that rotation
  294. bet is performed around axis x but the next rotation alp is performed
  295. around axis z which after application of previous transformation
  296. would move.
  297.  
  298.  
  299.                              +              +----+
  300.                     z /     /|             /    / alp
  301.                     /      + V bet            \
  302.        ------------+-------|---- x     ---------+------------ x
  303.                  /         |                      \
  304.            +----+                                  z
  305.            | /  |
  306.            |    V alp
  307.  
  308.                                 pic 2.8
  309.  
  310.  
  311. The implication is: we have to know with respect to what each
  312. next rotation is performed, that is what is the order of applications
  313. of 2D rotations. Let's think how we are coordinating the surrounding
  314. world? First we think of the direction. Then we can tilt our head up
  315. and down, and finally from there we can shake it from left to right.
  316. Now, when we are tilting head we have already chosen the direction.
  317. If we first would tilt head then the directional rotation to be
  318. performed after would not be parallel to the ground but rather
  319. perpendicular to the imaginary line being continuation of our
  320. tilted neck. ( No responsibility hereby assumed for all direct
  321. or consequential damage or injury resulted from doing that ) So
  322. it all comes to directional rotation first - gam in our case.
  323. pitch second - bet;  roll last - alp. Let's do just that using
  324. obtained formula for 2D rotation:
  325.  
  326.  
  327.   ^Z                     ^Y'                    ^Y"
  328.   |                      |                      |
  329.   |<---+                 |<---+                 |<---+
  330.   | gam|                 | bet|                 | alp|
  331.   |    |                 |    |                 |    |
  332.  -+----+--->X           -+----+--->Z'          -+----+--->X"
  333.   |                      |                      |
  334.  
  335.  
  336.  x'=z*sin(gam)+x*cos(gam)
  337.  y'=y
  338.  z'=z*cos(gam)-x*sin(gam)
  339.  
  340.                       x"=x
  341.                       y"=y*cos(bet)-z*sin(bet)
  342.                       z"=y*sin(bet)+z*cos(bet)
  343.  
  344.                                               x"'=y*sin(alp)+x*cos(alp)
  345.                                               y"'=y*cos(alp)-x*sin(alp)
  346.                                               z"'=z
  347.  
  348.                            pic 2.9
  349.  
  350. That's basically it, using those formulas we can apply 3D rotation
  351. to coordinates x,y,z and get x"',y"',z"'. We can see here 12
  352. multiplication's, The question is can we reduce this number? Let's
  353. try getting rid of temporary variables x',y',z',x",y",z". We can do it
  354. by just substituting each occurrence of them by their real meaning:
  355.  
  356. First let's try obtaining x",y" and z" expressed via x,y,z:
  357.  
  358.  
  359. y"=y'*cos(bet) - z'*sin(bet)
  360. y"=y*cos(bet) - (z*cos(gam)-x*sin(gam)) *sin(bet)
  361. y"=y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet)
  362. - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  363. x"=x'
  364. x"= z*sin(gam) + x*cos(gam)
  365. - - - - - - - - - - - - - - -
  366. z"=y'*sin(bet) + z'*cos(bet)
  367. z"=y*sin(bet) + (z*cos(gam)-x*sin(gam)) *cos(bet)
  368. z"=y*sin(bet) + z*cos(gam)*cos(bet) - x*sin(gam)*cos(bet)
  369. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  370.  
  371.  
  372. Now using that we can express x"',y"' and z"' via x,y,z:
  373.  
  374.  
  375. y"'=y"*cos(alp) - x"*sin(alp)
  376.  
  377. y"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*cos(alp) -
  378.     (z*sin(gam) + x*cos(gam)) *sin(alp))
  379.  
  380. y"'=y*cos(bet)*cos(alp)-z*cos(gam)*sin(bet)*cos(alp)+x*sin(gam)*sin(bat)*cos(alp)
  381.     -z*sin(gam)*sin(alp) - x*cos(gam)*sin(alp)
  382.  
  383. y"'=x* [ sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)] +
  384.     y* [ cos(bet)*cos(alp) ] +
  385.     z* [-cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)]
  386. -----------------------------------------------------------
  387.  
  388. x"'= y"*sin(alp) + x"*cos(alp)
  389.  
  390. x"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*sin(alp) +
  391.     (z*sin(gam) + x*cos(gam)) *cos(alp)
  392.  
  393. x"'=y*cos(bet)*sin(alp)-z*cos(gam)*sin(bet)*sin(alp)+x*sin(gam)*sin(bet)*sin(alp)
  394.     z*sin(gam)*cos(alp) + x*cos(gam)*cos(alp)
  395.  
  396. x"'=x* [sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)] +
  397.     y* [cos(bet)*sin(alp) ] +
  398.     z* [sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp) ]
  399. ----------------------------------------------------------
  400.  
  401. z"'=z"
  402. z"'=x* [- sin(gam)*cos(bet)] +
  403.     y* [sin(bet)] +
  404.     z* [cos(gam)*cos(bet)]
  405. ----------------------------------------------------------------
  406.  
  407.  
  408. This really appear to be more complicated that what we had before.
  409. It appear to have much more multiplications than we had. But if we
  410. look on those resulting formulas we can see that all coefficients
  411. in square brackets can be calculated just once so that transformation
  412. of a point would look like:
  413.  
  414.  
  415.                    x"'=x*mx1 + y*my1 + z*mz1
  416.                    y"'=x*mx2 + y*my2 + z*mz2
  417.                    z"'=x*mx3 + y*my3 + z*mz3
  418.  
  419. It can be seen that it takes just 9 multiplication. Of course calculating
  420. all the coefficients takes 16 multiplications:
  421.  
  422. mx1= sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)
  423. my1= cos(bet)*sin(alp)
  424. mz1= sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp)
  425.  
  426. mx2= sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)
  427. my2= cos(bet)*cos(alp)
  428. mz2= -cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)
  429.  
  430. mx3= -sin(gam)*cos(bet)
  431. my3= sin(bet)
  432. mz3= cos(gam)*cos(bet)
  433.  
  434. But if we have 100 points to turn. Original method would take
  435. 12*100=1200 multiplications. This one would take 9*100+16=916 since
  436. coefficients are calculated just once for all points to transform. One
  437. more consideration having to do with optimization, in most cases it makes
  438. sense to measure angles as integers. We don't usually care of fractions
  439. smaller then 1 degree. So we don't really need sin and cos for all the
  440. real number range. That's why we can calculate sin/cos just once for
  441. all 360 degrees before any rotations are performed and put the obtained
  442. values into arrays so that next time we would need any of them we
  443. wouldn't have to call a functions which probably uses a power series
  444. with few multiplications to calculate each. (Processors nowadays come
  445. with floating point units which can calculate sin/cos pretty fast
  446. but unlikely faster then single array access (this might become false
  447. though). By the way we don't really need 360 degrees. This number is
  448. not very convenient. We can go with full angle divided into 256 pseudo
  449. degrees. The reason we would need just one unsigned char to store the
  450. angle and whenever we go beyond 255 we simple wrap around to zero,
  451. saving this way a conditional statement we would need otherwise in
  452. the case of 360 degrees.
  453.  
  454.  
  455. #define T_RADS                  40.7436642  /* pseudo grads into rads */
  456.  
  457. float T_mx1,T_mx2,T_mx3,T_my1,T_my2,T_my3,T_mz1,T_mz2,T_mz3;
  458. float T_sin[256],T_cos[256];                /* precalculated */
  459.  
  460. void T_init_math(void)
  461. {
  462.  int i;
  463.  
  464.  for(i=0;i<256;i++)
  465.  {
  466.   T_sin[i]=sin(i/T_RADS);
  467.   T_cos[i]=cos(i/T_RADS);
  468.  }
  469. }
  470.  
  471.  
  472. Now finally to the code to perform 3D rotation. First function is building
  473. set of rotation coefficients. Note that T_mx1, T_mx2 etc. are global 
  474. variables. The other function: T_rotation uses those coefficients, so 
  475. they better be ready when later is called.
  476.  
  477.  
  478. void T_set_rotation_gam_bet_alp(int alp,int bet,int gam)
  479. {
  480.  T_cosalp=T_cos[alp];
  481.  T_sinalp=T_sin[alp];
  482.  T_cosbet=T_cos[bet];
  483.  T_sinbet=T_sin[bet];
  484.  T_cosgam=T_cos[gam];
  485.  T_singam=T_sin[gam];                       /* initializing */
  486.  
  487.  T_mx1=T_singam*T_sinbet*T_sinalp + T_cosgam*T_cosalp;
  488.  T_my1=T_cosbet*T_sinalp;
  489.  T_mz1=T_singam*T_cosalp - T_cosgam*T_sinbet*T_sinalp;
  490.  
  491.  T_mx2=T_singam*T_sinbet*T_cosalp - T_cosgam*T_sinalp;
  492.  T_my2=T_cosbet*T_cosalp;
  493.  T_mz2=-T_cosgam*T_sinbet*T_cosalp - T_singam*T_sinalp;
  494.  
  495.  T_mx3=-T_singam*T_cosbet;
  496.  T_my3=T_sinbet;
  497.  T_mz3=T_cosgam*T_cosbet;                   /* calculating the matrix */
  498. }
  499.  
  500.  
  501. void T_rotation(int *from,register int *to,int number)
  502. {
  503.  register int i;
  504.  register int xt,yt,zt;
  505.  
  506.  for(i=0;i<number;i++)
  507.  {
  508.   xt=*from++;
  509.   yt=*from++;
  510.   zt=*from++;
  511.  
  512.   *to++=(int)(T_mx1*xt+T_my1*yt+T_mz1*zt);
  513.   *to++=(int)(T_mx2*xt+T_my2*yt+T_mz2*zt);
  514.   *to++=(int)(T_mx3*xt+T_my3*yt+T_mz3*zt);  /* performing the operation */
  515.  }
  516. }
  517.  
  518.  
  519.  
  520. Representing transformations in matrix form.
  521. --------------------------------------------
  522. Matrix is a table of values:
  523.  
  524.       | a b c |
  525.     A=| d e f |
  526.       | g i j |
  527.  
  528. Special cases are column and row vectors:
  529.  
  530.                     | x |
  531.     V=| a b c |   C=| y |
  532.                     | z |
  533.  
  534.  
  535. Row can by multiplied by a column producing a value, a scalar:
  536.  
  537.               | x |
  538.     | a b c |*| y | = a*x + b*y + c*z
  539.               | z |
  540.  
  541. Matrix can be multiplied by another matrix: each value of the 
  542. result was obtained by multiplying respective rows and columns 
  543. of the arguments (since columns and vectors are to have the same 
  544. dimension in order for multiplication on them to be defined it 
  545. automatically places restrictions on dimensions of matrices that 
  546. can be multiplied):
  547.  
  548.                                  | x |            | k |
  549.                      | | a b c |*| y |  | a b c |*| l | |
  550.                      |           | z |            | m | |
  551.                      |                                  |
  552.  | a b c | | x k |   |           | x |            | k | |
  553.  | d e f |*| y l | = | | d e f |*| y |  | d e f |*| l | |
  554.  | g i j | | z m |   |           | z |            | m | |
  555.                      |                                  |
  556.                      |           | x |            | k | |
  557.                      | | g i j |*| y |  | g i j |*| l | |
  558.                                  | z |            | m |
  559.  
  560. We may be interested in multiplying a row vector by a 3x3 matrix:
  561.  
  562.            | a b c |
  563.  | x y z |*| d e f | = | x*a+y*d+z*g  x*b+y*e+z*i  x*c+y*f+z*j |
  564.            | g i j |
  565.  
  566. Looks familiar? Yes, the rotation coefficients we've calculated
  567. lie nicely in the 3x3 matrix and so we can think of 3D rotation
  568. as of matrix multiplication of "rotation" matrix onto original
  569. vector producing result vector. We can also try fitting translation 
  570. and scaling into a matrix, to do that for translation, however, we 
  571. would have to increase dimensions of the coefficient matrices and 
  572. the argument vector:
  573.  
  574. Translation:
  575.  
  576.              | 1  0  0  0 |
  577.  | x y z 1 |*| 0  1  0  0 | = | x+tx  y+ty  z+tz  1 |
  578.              | 0  0  1  0 |
  579.              | tx ty tz 1 |
  580.  
  581. Scaling:
  582.  
  583.              | sx  0  0 0 |
  584.  | x y z 1 |*|  0 sy  0 0 | = | x*sx  y*sy  z*sz  1 |
  585.              |  0  0 sz 0 |
  586.              |  0  0  0 1 |
  587.  
  588. The advantage of the matrix approach is that if we have several
  589. transformations each represented in matrix form:
  590.  
  591.             | X'| = (| X |*| A |) *| B |
  592.  
  593. Where A, and B are transformation matrixes and | X | the argument
  594. vector to be transformed we can:
  595.  
  596.              | X'| = | X |* ( | A |*| B | )
  597.  
  598. calculate the "concatenated" matrix | K | = | A |*| B | and each
  599. following coordinate just transform as:
  600.  
  601.                     | X'| = | X |*| K |
  602.  
  603. By the way by representing each of 2-D rotations in 3x3 matrix form we
  604. can get the formulas we have found for 3-D rotation by calculating a
  605. concatenated matrix.
  606.  
  607. So, matrix approaches might be very effective in instances where
  608. there are a lot of consecutive transformation. For example if we
  609. have consecutive rotations it makes sense to calculate their
  610. concatenated matrix and use it from then on. It would save
  611. us quite a few multiplication, one verses several consecutive matrix
  612. multiplication for each point to transform.
  613.  
  614. And in the run-time it is preferably to use pre-calculated
  615. formulas for concatenated matrix coefficients then trying to get it
  616. every time by multiplying respective transformational matrices.
  617. The reason 3 consecutive matrix multiplication of, for instance, 3x3
  618. matrices take 9+9+9=27 multiplication's. But since lot of values in
  619. those matrices are 0, it is taken into account in pre-calculated
  620. coefficient formulas which allow to produce the very same results
  621. with just 16 multiplications.
  622.  
  623. As to the translation transformation matrix and especially increasing
  624. the dimension from 3 to 4 it looks much less natural. 4x4 matrix means
  625. more multiplications too, of course some of the value in 4x4 would
  626. always be zero so we can optimize (cripple) the routine to original 9
  627. multiplications and 3 additions. In one word, there's a freedom
  628. to decide what approach to take.
  629.  
  630.  
  631.  
  632. perspective.
  633. ------------
  634. I believe, it is not hard to visualize what perspective is, perhaps the
  635. first association lot of us would have, is this empty straight street
  636. with identical buildings on both sides disappearing in infinity. ( At
  637. least people residing in newer neighborhoods of Moscow or Toronto know
  638. what I am talking about ) Why do we need this transformation? Just to
  639. achieve landscape realism, after all this same street would look pretty
  640. weird if the road would not collapse into a dot and buildings further
  641. away from us won't look smaller. So weird it might be even hard to
  642. picture.
  643.  
  644. Perspective might not be needed for CAD style applications, but even
  645. there it would, most likely, be implemented in one viewing mode or
  646. another.
  647.  
  648. Before we can implement this transformation it is a good idea to understand
  649. physics of what is happening. A ray from some point in space is being
  650. caught by the viewer's eye. This ray had intersected a plane located
  651. in front of the viewer. If we can find an intersection point and plot the
  652. dot there, for the viewer it would be the same as if the ray from the
  653. plotted dot was actually coming from the initial point in space.
  654.  
  655.  
  656.                              * A
  657.                        |   /
  658.                      A'| /
  659.                        *
  660.                      / | B'
  661.                    *---*-----------------------* B
  662.                      \ |
  663.                        *\
  664.                      C'|   \
  665.                        |      \
  666.                                 * C
  667.  
  668.                               pic 2.10
  669.  
  670.  
  671. Rays from different points in space , in this model, all converge
  672. in one point - viewer's eye. The math here is quite straightforward 
  673. pic 2.11:
  674.  
  675.  
  676.                   screen plane->|   *
  677.                                 |  /| A (X,Z)
  678.                                 | / |
  679.                                 |/  |
  680.                             ^   *   |
  681.                             |  /|X' |
  682.                               / |   |
  683.                             |/  |   |
  684.                             *---+---+- - - - ->
  685.                             0   |   Z
  686.                             |   |
  687.            focus distance ->|---|<-
  688.  
  689.                              pic 2.11
  690.  
  691.  
  692. Let's consider 2-D case first: the viewer's eye is located at point 0
  693. the distance between viewer's eye and the screen would be called -
  694. "focus". The problem is to determine at what point the ray going from
  695. point A into the viewer's eye would intersect the screen. We would have
  696. to plot the dot just there. This is very simple, just yet another
  697. case of solving similar triangles. Just from the fact that the angle
  698. at 0 is the same for both triangles, and if two angles are the same their
  699. tangents would be too. We have:
  700.  
  701.            X'      X                          X*focus
  702.          ----- = -----         =>        X'= -----------
  703.          focus     Z                             Z
  704.  
  705. Same considerations apply for Y dimension, together describing 3-D case:
  706.  
  707.           Y'=Y*focus/Z                   X'=X*focus/Z
  708.  
  709. It is a bit of a mystery what is happening with Z coordinate. Basically
  710. after the perspective transformation we would be performing rendering
  711. onto the screen which is 2-D surface we would need X and Y but hardly
  712. Z, so we might as well have Z'=0. However when trying to render multiface
  713. objects we would need to know the depth (Z) of all polygons so that
  714. we can decide which are visible and which are obscured. So we might
  715. have Z'=Z that is just leave old value. Then again by doing transformation
  716. on X and Y and leaving Z unchanged we actually may allow for depth
  717. relation to change. That is a polygon obscuring another one in the original
  718. space, before projection, can actually come out behind or partly behind
  719. with this kind of transformation pic 2.21.
  720.  
  721.             ^ Z
  722.   same      |                   *                     *
  723.  deapth     |               * /                 * /
  724.  before     |           A / /                //
  725.   and       |           * /               /* A
  726.  after      |           /B            / B
  727.             |          *            *
  728.             |        before                 after
  729.  
  730.                               pic 2.12
  731.  
  732. To avoid that we can apply some non linear transformation on Z,
  733. like Z'=C1+C2/Z, where C1, C2 are some constants.
  734.  
  735.  
  736. In the enclosed code only X and Y are being transformed. The implementation
  737. of all those things is that easy I won't even put a code example here.
  738. However, there are few possible questions: first: what should be the value
  739. for "focus distance" ?
  740.  
  741.  
  742.          \                    /          \              /
  743.            \                /             \            /
  744.              \            /                \          /
  745.            ---I\--------/I---           ---I\--------/I---
  746.                  \    /                      \      /
  747.                    \/                         \    /
  748.                                                \  /
  749.                                                 \/
  750.  
  751.                                pic 2.13
  752.  
  753.  
  754. pic 2.13 shows it's physical sense, focus basically determines the view
  755. angle. when it is smaller - angle is wider, when it is bigger angle
  756. is smaller. It is helpful to think of this distance as being measured
  757. in screen pixels, this way for 320 pixel display if focus is chosen
  758. to be 180 the view angle we are getting is 90'. Perspective
  759. transformation might produce at first, weird looking results, with
  760. distortions sometimes similar to that of wide-range photographic
  761. pictures. One has to experiment a bit to find what's looking best,
  762. values giving view angle somewhere between 75-85' would probably
  763. look not that terrible, but that depends of scene and screen geometries.
  764.  
  765.  
  766.  
  767. clipping problem...
  768. -------------------
  769. Where transforms point with Z equal to 0? since we are dividing by
  770. Z that would be infinity, or at least something producing divide error,
  771. in this computer-down-to-earth-world. The only way to deal with this
  772. is to guarantee that there would be no coordinates with such Zs. another
  773. thing, calculations for points with negative Z would produce negative
  774. coordinates, what we would see is objects (or even parts thereof) flipped
  775. over, but then again objects with negative coordinates are behind the
  776. viewing plane and effectively behind the viewer, so this is something
  777. we can't see and better make sure this is not being rendered. The way
  778. to do it is by applying 3-D clipping.
  779.  
  780. The perspective transformation  can also be represented as a matrix,
  781. however, let's think of all transformations discussed up to this moment
  782. a point would have to go through in order to be rendered. First we would
  783. rotate the world with respect to viewer orientation, then before we can
  784. apply perspective transformation we would have to do clipping and
  785. at least get rid of the vertices behind the viewing plane, and only then
  786. can we apply perspective transformation. Representing clipping in
  787. matrix form is not the most pleasant thing there is. And we've already 
  788. discussed that matrices would be of use when we have several consecutive 
  789. transformations so that all of them would be represented by one matrix. 
  790. In this case however rotation and perspective are separated by the clipping. 
  791. If we are sure that none of the coordinates would ever get behind the
  792. viewing plane, we don't actually need clipping, that's the instance to
  793. represent everything in the matrix form, otherwise it might not be
  794. a good idea neither in terms of code complicity nor performance.
  795.  
  796. fixed point math.
  797. -----------------
  798. In previous chapters in several places we didn't have another choice but to
  799. use floating point multiplications. (sin and cos are after all real-valued
  800. functions) However this is still fairly expensive operation, integer
  801. multiplications and divisions would usually be faster. However since we
  802. have to represent numbers with fractions it appeared we didn't have another
  803. choice but to use floating point operations. The question is, can we
  804. represent fractions using just integers? One might think there's a
  805. contradiction in the question itself, however we would see that, as in
  806. many other cases, we can use one mechanism to represent something else.
  807. This is the case here, We may use fixed point numbers representing them
  808. as integers, and with small amendments we can use regular integer additions
  809. multiplications etc, to perform fixed point operations. First of all how do 
  810. we represent integers?
  811.  
  812.  
  813.                    +------+------+------+------+------+
  814.                    |   1  |   0  |  0   |  0   |   1  |
  815.                    +------+------+------+------+------+
  816.                      4       3       2      1      0
  817.                     2 =16   2 =8    2 =4   2 =2   2 =1
  818.  
  819.                                 pic 2.14
  820.  
  821. With any counting system, digits in multi digit numbers would have
  822. different weight, so to speak. With decimal numbers this weight
  823. would be some power of ten, that is 102 is actually:
  824.  
  825.        10**2*1 + 10**1*0 + 10**0*2 == 100+0+2 ==102
  826.  
  827. Same with binary numbers, the only difference, weight would be
  828. some power of 2. this way number on pic 2.13 would be:
  829.  
  830.  2**4*1 + 2**3*0 + 2**2*0 + 2**1*0 + 2**0*1 == 16+0+0+0+1 == 17
  831.  
  832. Now, it is quite natural for us to place a decimal point into a decimal
  833. number, placing a binary point into a binary number should not
  834. seem more difficult. How do we interpret the decimal point? It is
  835. just that for the numbers to the right of the decimal point we
  836. pick negative power of ten. That is:
  837.  
  838.   102.15 == 10**2*1 + 10**1*0 + 10**0*2 + 10**-1*1 + 10**-2*5 ==
  839.                     1      5             15
  840.   == 100 + 0 + 2 + ---- + ---- == 102 + -----
  841.                     10    100            100
  842.  
  843. In this representation the point separates negative powers from
  844. zero power and unlike numbers represented in exponential form - the
  845. floating point numbers, in this case the point never changes it's
  846. position.
  847.  
  848. Nothing should stop us from extending the same scheme onto binary
  849. numbers:
  850.  
  851.                    +------+------+------+------+------+
  852.                    |  1   |   0  |   0  |   1  |   1  |
  853.                    +------+------+------+------+------+
  854.                      2       1     0     -1     -2
  855.                     2 =4    2 =2  2 =1   2 =1/2 2 =1/4
  856.  
  857.                                 pic 2.15
  858.  
  859. Using the very same technique we can see that number represented
  860. on the pic 2.15 can be considered also as:
  861.  
  862.   2**2*1 + 2**1*0 + 2**0*0 + 2**-1*1 + 2**-2*1 == 4+0+0+1/2+1/4
  863.  
  864.   == 4 + 3/4
  865.  
  866. Numbers in what range can we represent? Resorting once again to
  867. analogy with decimal numbers it can be seen that two digits
  868. to the right from the decimal point can cover the range of 0.99
  869. in 0.01 steps. Numbers smaller then minimal 0.01 precision step
  870. just can't be represented without increasing number of digits
  871. after the decimal point. Very same thing is happening with binaries
  872. in the example at pic 2.15 since we have two binary digits- minimal
  873. number we can represent is 1/4 ( 0.01 (bin) == 1/4 (dec) ) and
  874. again numbers with higher precision can be represented only by
  875. increasing number of binary digits after the binary point.
  876.  
  877. Now to the representation of negative numbers, there actually
  878. several ways to do that but effectively today one has prevailed,
  879. this is so called 2's compliment form. The idea is to negate
  880. the weight of the leftmost digit pic 2.16:
  881.  
  882.                    +------+------+------+------+------+
  883.                    |   1  |   1  |   1  |   1  |   1  |
  884.                    +------+------+------+------+------+
  885.                      4       3       2      1      0
  886.                    -2 =-16  2 =8    2 =4   2 =2   2 =1
  887.  
  888.                                 pic 2.16
  889.  
  890. the number represented above is considered to be:
  891.  
  892.   -2**4*1 + 2**3*1 + 2**2*1 + 2**1*1 + 2**0*1== -16+8+4+2+1 = -1
  893.  
  894. The advantage of this approach, we don't have to change addition
  895. and subtraction algorithm working for positive integers in order
  896. to accommodate 2's compliment numbers. Why is that? just by the
  897. virtue of natural overwrap of integer numbers the way they are
  898. represented in the computers. If we add 1 to the maximum number
  899. we can represent we would obtain a carry from the most significant
  900. (leftmost) bit, which is ignored, and the value of exactly 0.
  901. -1 in signed representation is exactly the maximum possible value
  902. to represent in unsigned representation, -1+1==0 indeed. (One can
  903. check this to be true for all numbers and by more formal means, but
  904. since this is not the point I would ignore it) Again, addition and
  905. subtraction don't have to be changed to accommodate 2's compliment
  906. negative numbers (multiplication however should, that's why there
  907. are instructions for both signed and unsigned multiplications
  908. in most computers)
  909.  
  910. Since the leftmost digit in 2's compliment representation carries
  911. negative weight and that's exactly the one with highest power the
  912. minimum negative number possible to represent in the example
  913. above would be 10000 (bin) = -16 all the other digits have
  914. positive weights so maximum possible positive number would be
  915. 01111 (bin) = 15 but this slight asymmetry is not a problem in
  916. the majority of cases.
  917.  
  918. However, values of sin and cosin functions are between 1 and -1 so
  919. to represent them we might choose format with just one integer
  920. field:
  921.  
  922.                    +------+------+------+------+------+
  923.                    |   0  |   1  |   1  |   1  |   1  |
  924.                    +------+------+------+------+------+
  925.                      0      -1      -2     -3     -4
  926.                    -2 =-1   2 =1/2  2 =1/4 2 =1/8 2 =1/16
  927.  
  928.                                  pic 2.17
  929.  
  930. but, because of the asymmetry described above there would be a
  931. value for -1 (1000) but there would be no pure 1, just it's
  932. approximation (01111): (pic 2.17) 1/2+1/4+1/8+1/16, Then again, for
  933. most graphics applications when we have, say, 15 fields representing
  934. fractions this would be close enough and won't cause a lot of problems.
  935. As you can see, the very same rules work for 2's compliment whether
  936. with fixed point or without.
  937.  
  938.  
  939.  
  940. multiplication of fixed point numbers.
  941. --------------------------------------
  942. Let's consider what is happening when we are multiplying two
  943. decimal numbers, for example we have to multiply an integer
  944. by a number with a decimal point and we would need just
  945. an integer as a result:
  946.  
  947.                             1.52
  948.                         *     11
  949.                            -----
  950.                             1 52
  951.                         +  15 2
  952.                            -----
  953.                            16.72
  954.  
  955. Actual result of multiplications has as many digits after
  956. the decimal point as there were in both arguments. Since
  957. we need just an integer as a result we would discard numbers
  958. after the point effectively shift digits to the right by
  959. 2 in this case. We can do the same with fixed point binary
  960. numbers, it means that we have to readjust precision after
  961. the multiplication. Say, if we are multiplying an integer
  962. by a number having 8bit considered to be after the binary
  963. point the result would also 8bit after the point. If it
  964. is just the integer part we are interested in, we have to
  965. shift the result right 8 bit destroying all fractional
  966. information. Similar things are happening with division
  967. except that if we want to divide two integers and obtain
  968. a fixed point result the argument to be divided has to
  969. be shifted left- added fixed point fields, so that effectively
  970. we are dividing a fixed point number by an integer.
  971.  
  972. How to implement all that? Well, first of all we have to
  973. decide what kind of values and what operations on them
  974. we would need. Sometimes it is nice just to have general
  975. purpose fixed point math library, on the other hand we
  976. may choose several different formats of fixed point numbers
  977. one for every special place. Another thing we have to
  978. consider is: would we be able to represent all the operations
  979. using just high level C constructs or it would be necessary
  980. to descend into assembly instructions. The reason we might
  981. consider later is: in most assemblies the result of
  982. multiplication has twice as many bits as each of arguments,
  983. and since we would occupy some bits with fractions we really
  984. need all bits we can get. Another thing, the result of
  985. multiplication would often be located in two registers, So
  986. we can organize operations the way so that instead of adjusting
  987. result by shifting we would just take it from the register
  988. carrying higher bits, -effectively doing zero-cost shifting
  989. right by the bit length of the register. On the other hand,
  990. do we really want to compromise code portability by coding
  991. in assembly? Obviously, there are valid reasons to answer that
  992. both ways. In particular case of this text and associated
  993. code we are interested in good portability so let's try to
  994. stay within straight C boundaries.
  995.  
  996. Where would we be using fixed point? first of all - in 3-D
  997. transformations, to keep rotation matrix coefficients. We
  998. would need two kinds of multiplication: In calculation of
  999. the coefficients: fixed * fixed producing same precision
  1000. fixed and in actual coordinate transformation fixed * int
  1001. producing int. lets define fixed to be:
  1002.  
  1003.  
  1004.  
  1005. typedef signed_32_bit fixed;                /* better be 32bit machine */
  1006. #define T_P                             10  /* fixed math precision */
  1007.  
  1008.  
  1009.  
  1010. the T_P would be the number of bits we are allocating to
  1011. keep fractions. Let's make a little function to transform the
  1012. floating point values in the range [-1,1] into fixed point
  1013. values:
  1014.  
  1015.  
  1016.  
  1017. fixed TI_float2fixed(float a)
  1018. {
  1019.  fixed res=0;
  1020.  int i,dv,sign;
  1021.  
  1022.  if((sign=a<0.0)==1) a=-a;
  1023.  
  1024.  for(i=1,dv=2;i<=T_P;i++,dv*=2)
  1025.   if(a>=1.0/dv) { a-=1.0/dv; res|=(0x1<<(T_P-i)); };
  1026.  
  1027.  if(sign==1) res=-res;
  1028.  
  1029.  return(res);
  1030. }
  1031.  
  1032.  
  1033.  
  1034. In the function above we are iteratively subtracting powers
  1035. of two from the argument and using the result to set the bits
  1036. in the fixed point number. We would use this function to
  1037. precalculate the array of sin/cos values. 
  1038.  
  1039. When calculating the rotational coefficients we would multiply
  1040. one fixed by another and obtain same precision fixed result.
  1041. again, actual result would have twice as many fractional bits
  1042. as each of the argument so same precision result would be
  1043. (a*b)>>T_P where 'a' and 'b' fixed point values and T_P number
  1044. of bits after the binary point in each. This way code would
  1045. look something like:
  1046.  
  1047.                                 ...
  1048.                                 ...
  1049.  
  1050.  T_wx1=((singam*((sinbet*sinalp)>>T_P))>>T_P) + ((cosgam*cosalp)>>T_P);
  1051.  T_wy1=((cosbet*sinalp)>>T_P);
  1052.  T_wz1=((singam*cosalp)>>T_P) - ((cosgam*((sinbet*sinalp)>>T_P))>>T_P);
  1053.  
  1054.                                 ...
  1055.                                 ...
  1056.  
  1057.  
  1058. In rotation function ints would have to be multiplied by a fixed
  1059. coefficient, the result would have T_P fractional bits, since we
  1060. are interested in just integer part same rule as above applies:
  1061.  
  1062.  
  1063.                                 ...
  1064.                                 ...
  1065.  
  1066.  *to++=((T_wx1*xt)>>T_P) + ((T_wy1*yt)>>T_P) + ((T_wz1*zt)>>T_P);
  1067.  *to++=((T_wx2*xt)>>T_P) + ((T_wy2*yt)>>T_P) + ((T_wz2*zt)>>T_P);
  1068.  *to++=((T_wx3*xt)>>T_P) + ((T_wy3*yt)>>T_P) + ((T_wz3*zt)>>T_P);
  1069.  
  1070.                                 ...
  1071.                                 ...
  1072.  
  1073. That's about it. The only thing, we really have to be careful with
  1074. the range of numbers the operations would work for. If we are working
  1075. with 32bit numbers and would choose 16 of them for fractional part
  1076. and one for sign bit then only 15bit would carry the integer part.
  1077. after multiplication of two fixed numbers since result has as many
  1078. fractional bits as in both arguments all 32 bits would actually carry
  1079. fractions. The integer part would be gone in the dark realm of left-carry,
  1080. and there's no way in standard straight C to get it out of there.
  1081. (when using assembly we don't have this problem since result
  1082. of multiplication has twice bits of each argument, so we really can
  1083. get trough to the integer part, however, most C compilers define the
  1084. result of multiplication to be of the same bit size as arguments
  1085. int*int==int, some compilers try to be smarter then that but we really
  1086. have to consider worst-case scenario)
  1087.  
  1088. And finally is it really true that integer multiplication plus shifts
  1089. is faster then floating multiplication? What follows are results of
  1090. some very approximate tests:
  1091.  
  1092.  
  1093.  
  1094. sparc:   floating faster fixed   about  1.3 times
  1095. 68040:   floating faster fixed   about  1.5 times
  1096. rs4000:  floating faster fixed   about  1.1 times
  1097.  
  1098. rs6000:  fixed faster floating   about  1.1 times
  1099. i386sx:  fixed faster floating   about  5.1 times
  1100.  
  1101. floating in the test was:    float a,b;    a*b;
  1102. fixed in the test was:       int a,b;      (a*b)>>15;
  1103.  
  1104.  
  1105.  
  1106. As one can see processors with fast floating point units nowadays
  1107. do floating operations faster then integer once. Then again in terms
  1108. of portability integer multiplication would always be fast
  1109. enough whereas floating multiplication especially on processors
  1110. without FPUs (i386sx) might get really slow. The decision whether
  1111. to implement fixed point, this way, might get a bit tough to make.
  1112.  
  1113. In the provided library's source transformations are implemented 
  1114. both ways but with the same function prototypes, so that one can
  1115. decide which implementation is more suitable to link in.
  1116.  
  1117.                                   * * *