home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / IRIT / IRITS.ZIP / PRIMITIV.C < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-05  |  37.7 KB  |  972 lines
  1. /*****************************************************************************
  2. *   "Irit" - the 3d polygonal solid modeller.                     *
  3. *                                         *
  4. * Written by:  Gershon Elber                Ver 0.2, Mar. 1990   *
  5. ******************************************************************************
  6. *   Module to generate the geometric primitives defined in the system. The   *
  7. * primitives currently defined are:                         *
  8. * 1. BOX - main planes parallel box.                         *
  9. * 2. GBOX - generalized box - 6 arbitrary planes.                 *
  10. * 3. CYLIN - cylinder with any main direction.                     *
  11. * 4. CONE - cone with any main direction.                     *
  12. * 5. SPHERE                                     *
  13. * 6. TORUS - with any main direction.                         *
  14. * 7. PLANE - non closed, single polygon object: circle with resolution edges *
  15. * 8. POLY - directly define single polygon object by specifing its vertices. *
  16. *   In addition the following lower level operations are defined to create   *
  17. * objects - EXTRUDE, and SURFREV, both require a polygon and a vector to     *
  18. * extrude/rotate the polygon along.                         *
  19. *****************************************************************************/
  20.  
  21. #ifdef __MSDOS__
  22. #include <graphics.h>
  23. #endif /* __MSDOS__ */
  24.  
  25. #include <stdio.h>
  26. #include <math.h>
  27. #include "program.h"
  28. #include "allocatg.h"
  29. #include "objects.h"
  30. #include "primitvl.h"
  31. #include "primitvg.h"
  32. #include "geomat3d.h"
  33. #include "windowsg.h"
  34. #include "convexg.h"
  35.  
  36. #ifndef __MSDOS__
  37. #include "xgraphic.h"
  38. #endif /* __MSDOS__ */
  39.  
  40. /*****************************************************************************
  41. *   Routine to create a BOX geometric object defined by Pt - the minimun     *
  42. * 3d point, and Width - Dx Dy & Dz vector.        4             *
  43. * Order of vertices is as                           5       7             *
  44. * follows in the picture:                           |   6   |             *
  45. *                            |   |   |             *
  46. * (Note vertex 0 is hidden behind edge 2-6)        |    |   |             *
  47. *                            1   |   3                *
  48. *                            2             *
  49. *****************************************************************************/
  50. struct ObjectStruct * GenBOXObject(VectorType Pt, RealType *WidthX,
  51.                     RealType *WidthY, RealType *WidthZ)
  52. {
  53.     VectorType Dir1, Dir2, Dir3;
  54.  
  55.     PT_CLEAR(Dir1);    Dir1[0] = (*WidthX);   /* Prepare direction vectors. */
  56.     PT_CLEAR(Dir2);    Dir2[1] = (*WidthY);       /* Parallel to main axes. */
  57.     PT_CLEAR(Dir3);    Dir3[2] = (*WidthZ);           /* For GBOX call. */
  58.  
  59.     return GenGBOXObject(Pt, Dir1, Dir2, Dir3);
  60. }
  61.  
  62. /*****************************************************************************
  63. *   Routine to create a GBOX geometric object defined by Pt - the minimun    *
  64. * 3d point, and 3 direction Vectors Dir1, Dir2, Dir3. If two of the         *
  65. * direction vectors are parallel the GBOX converges to zero volume. A NULL   *
  66. * pointer is returned in that case!                         *
  67. *                             4             *
  68. * Order of vertices is as                           5       7             *
  69. * follows in the picture:                           |   6   |             *
  70. *                            |   |   |             *
  71. * (Note vertex 0 is hidden behind edge 2-6)        |    |   |             *
  72. *                            1   |   3                *
  73. *                            2             *
  74. *****************************************************************************/
  75. struct ObjectStruct * GenGBOXObject(VectorType Pt,
  76.             VectorType Dir1, VectorType Dir2, VectorType Dir3)
  77. {
  78.     int i;
  79.     VectorType Temp;
  80.     VectorType V[8];                  /* Hold 8 vertices of BOX. */
  81.     struct ObjectStruct *PBox;
  82.  
  83.     VecCrossProd(Temp, Dir1, Dir2);
  84.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  85.     VecCrossProd(Temp, Dir2, Dir3);
  86.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  87.     VecCrossProd(Temp, Dir3, Dir1);
  88.     if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
  89.  
  90.     /* Also the 0..7 sequence is binary decoded such that bit 0 is Dir1, */
  91.     /* bit 1 Dir2, and bit 2 is Dir3 increment:                 */
  92.     for (i=0; i<8; i++) {
  93.     PT_COPY(V[i], Pt);
  94.  
  95.     if (i & 1) { PT_ADD(V[i], V[i], Dir1); }
  96.     if (i & 2) { PT_ADD(V[i], V[i], Dir2); }
  97.     if (i & 4) { PT_ADD(V[i], V[i], Dir3); }
  98.     }
  99.  
  100.     PBox = GenGeomObject("", NULL, NULL); /* Generate the BOX object itself: */
  101.  
  102.     /* And generate the 6 polygons (Bottom, top and 4 sides in this order):  */
  103.     PBox -> U.Pl = GenPolygon4Vrtx(V[0], V[1], V[3], V[2], V[4], PBox -> U.Pl);
  104.     PBox -> U.Pl = GenPolygon4Vrtx(V[6], V[7], V[5], V[4], V[0], PBox -> U.Pl);
  105.     PBox -> U.Pl = GenPolygon4Vrtx(V[4], V[5], V[1], V[0], V[2], PBox -> U.Pl);
  106.     PBox -> U.Pl = GenPolygon4Vrtx(V[5], V[7], V[3], V[1], V[0], PBox -> U.Pl);
  107.     PBox -> U.Pl = GenPolygon4Vrtx(V[7], V[6], V[2], V[3], V[1], PBox -> U.Pl);
  108.     PBox -> U.Pl = GenPolygon4Vrtx(V[6], V[4], V[0], V[2], V[3], PBox -> U.Pl);
  109.  
  110.     SET_OBJECT_COLOR(PBox, PrimColor);           /* Set its default color. */
  111.  
  112.     return PBox;
  113. }
  114.  
  115. /*****************************************************************************
  116. *  Routine to fetch the resolution parameter from the RESOLUTION object.     *
  117. *****************************************************************************/
  118. static int GetResolution(void)
  119. {
  120.     int Resolution;
  121.     struct ObjectStruct *PObj = GetObject("RESOLUTION");
  122.  
  123.     if (PObj == NULL || !IS_NUM_OBJ(PObj)) {
  124.     WndwInputWindowPutStr("No numeric object name RESOLUTION is defined",
  125.                                     RED);
  126.     Resolution = DEFAULT_RESOLUTION;
  127.     }
  128.     else Resolution = MAX(((int) (PObj -> U.R)), MIN_RESOLUTION);
  129.  
  130.     Resolution = (Resolution / 2) * 2;        /* Make sure its an even number. */
  131.  
  132.     return Resolution;
  133. }
  134.  
  135. /*****************************************************************************
  136. *   Routine to create a CONE geometric object defined by Pt - the base       *
  137. * 3d center point, Dir - the cone direction and length, and base radius R.   *
  138. *****************************************************************************/
  139. struct ObjectStruct * GenCONEObject(VectorType Pt, VectorType Dir, RealType *R)
  140. {
  141.     int i, Resolution;
  142.     RealType Angle, AngleStep;
  143.     PointType LastCirclePt, CirclePt, ApexPt;
  144.     MatrixType Mat;
  145.     struct VertexStruct *VBase;
  146.     struct PolygonStruct *PBase;
  147.     struct ObjectStruct *PCone;
  148.  
  149.     Resolution = GetResolution();     /* Get refinement factor of object. */
  150.  
  151.     GenTransformMatrix(Mat, Pt, Dir, *R);     /* Transform from unit circle. */
  152.  
  153.     PT_COPY(ApexPt, Pt);           /* Find the apex point: Pt + Dir. */
  154.     PT_ADD(ApexPt, ApexPt, Dir);
  155.  
  156.     PCone = GenGeomObject("", NULL, NULL);   /* Gen. the CONE object itself: */
  157.     /* Also allocate the base polygon header with first vertex on it: */
  158.     PBase = AllocPolygon(0, 0, VBase = AllocVertex(0, 0, NULL, NULL), NULL);
  159.  
  160.     LastCirclePt[0] = 1.0;        /* First point is allways Angle = 0. */
  161.     LastCirclePt[1] = 0.0;
  162.     LastCirclePt[2] = 0.0;
  163.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  164.  
  165.     PT_COPY(VBase -> Pt, LastCirclePt);  /* Update first pt in base polygon. */
  166.  
  167.     AngleStep = M_PI * 2 / Resolution;
  168.  
  169.     for (i=1; i<=Resolution; i++) {          /* Pass the whole base circle. */
  170.     Angle = AngleStep * i;             /* Prevent from additive error. */
  171.  
  172.     CirclePt[0] = cos(Angle);
  173.     CirclePt[1] = sin(Angle);
  174.     CirclePt[2] = 0.0;
  175.  
  176.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  177.  
  178.     PCone -> U.Pl = GenPolygon3Vrtx(LastCirclePt, ApexPt, CirclePt,
  179.                             Pt, PCone -> U.Pl);
  180.     /* And add this vertex to base polygon: */
  181.     if (i == Resolution)           /* Its last point - make it circular. */
  182.         VBase -> Pnext = PBase -> V;
  183.     else {
  184.         VBase -> Pnext = AllocVertex(0, 0, NULL, NULL);
  185.         VBase = VBase -> Pnext;
  186.         PT_COPY(VBase -> Pt, CirclePt);
  187.     }
  188.  
  189.     PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
  190.     }
  191.  
  192.     UpdatePolyPlane(PBase, ApexPt);   /* Update base polygon plane equation. */
  193.     SET_CONVEX_POLY(PBase);               /* Mark it as convex polygon. */
  194.     PBase -> Pnext = PCone -> U.Pl;  /* And stick it into the cone polygons. */
  195.     PCone -> U.Pl = PBase;
  196.  
  197.     SET_OBJECT_COLOR(PCone, PrimColor);           /* Set its default color. */
  198.  
  199.     return PCone;
  200. }
  201.  
  202. /*****************************************************************************
  203. *   Routine to create a CYLINder geometric object defined by Pt - the base   *
  204. * 3d center point, Dir - the cylin direction and length, and base radius R.  *
  205. *   The second base is defined from first one by translating it by vector    *
  206. * Dir, and its points are prefixed with T.                     *
  207. *****************************************************************************/
  208. struct ObjectStruct * GenCYLINObject(VectorType Pt, VectorType Dir, RealType *R)
  209. {
  210.     int i, Resolution;
  211.     RealType Angle, AngleStep;
  212.     PointType LastCirclePt, CirclePt, TLastCirclePt, TCirclePt, TPt;
  213.     MatrixType Mat;
  214.     struct VertexStruct *VBase1, *VBase2 = NULL;
  215.     struct PolygonStruct *PBase1, *PBase2;
  216.     struct ObjectStruct *PCylin;
  217.  
  218.     Resolution = GetResolution();     /* Get refinement factor of object. */
  219.  
  220.     GenTransformMatrix(Mat, Pt, Dir, *R);     /* Transform from unit circle. */
  221.  
  222.     PCylin = GenGeomObject("", NULL, NULL); /* Gen. the CYLIN object itself: */
  223.     /* Also allocate the bases polygon header with first vertex on it: */
  224.     PBase1 = AllocPolygon(0, 0, VBase1 = AllocVertex(0, 0, NULL, NULL), NULL);
  225.     PBase2 = AllocPolygon(0, 0, VBase2 = AllocVertex(0, 0, NULL, NULL), NULL);
  226.  
  227.     PT_ADD(TPt, Pt, Dir);           /* Translated circle center (by Dir). */
  228.  
  229.     LastCirclePt[0] = 1.0;        /* First point is allways Angle = 0. */
  230.     LastCirclePt[1] = 0.0;
  231.     LastCirclePt[2] = 0.0;
  232.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  233.     PT_COPY(VBase1 -> Pt, LastCirclePt);/* Update first pt in base1 polygon. */
  234.     PT_ADD(TLastCirclePt, LastCirclePt, Dir); /* Translated circle (by Dir). */
  235.     PT_COPY(VBase2 -> Pt, TLastCirclePt);/* Update first pt in base2 polygon.*/
  236.  
  237.     AngleStep = M_PI * 2 / Resolution;
  238.  
  239.     for (i=1; i<=Resolution; i++) {          /* Pass the whole base circle. */
  240.     Angle = AngleStep * i;             /* Prevent from additive error. */
  241.  
  242.     CirclePt[0] = cos(Angle);
  243.     CirclePt[1] = sin(Angle);
  244.     CirclePt[2] = 0.0;
  245.  
  246.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  247.     PT_ADD(TCirclePt, CirclePt, Dir);     /* Translated circle (by Dir). */
  248.  
  249.     PCylin -> U.Pl = GenPolygon4Vrtx(TLastCirclePt, TCirclePt, CirclePt,
  250.                      LastCirclePt, Pt, PCylin -> U.Pl);
  251.     /* And add this vertices to the two cylinder bases:              */
  252.     /* Note Base1 is build forward, while Base2 is build backward so it  */
  253.     /* will be consistent - cross product of 2 consecutive edges will    */
  254.     /* point into the model.                         */
  255.     if (i == Resolution) {           /* Its last point - make it circular. */
  256.         VBase1 -> Pnext = PBase1 -> V;
  257.         VBase2 -> Pnext = PBase2 -> V;
  258.     }
  259.     else {
  260.         VBase1 -> Pnext = AllocVertex(0, 0, NULL, NULL);
  261.         VBase1 = VBase1 -> Pnext;
  262.         PT_COPY(VBase1 -> Pt, CirclePt);
  263.         PBase2 -> V = AllocVertex(0, 0, NULL, PBase2 -> V);
  264.         PT_COPY(PBase2 -> V -> Pt, TCirclePt);
  265.     }
  266.  
  267.     PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
  268.     PT_COPY(TLastCirclePt, TCirclePt);
  269.     }
  270.  
  271.     UpdatePolyPlane(PBase1, TPt);     /* Update base polygon plane equation. */
  272.     SET_CONVEX_POLY(PBase1);               /* Mark it as convex polygon. */
  273.     PBase1 -> Pnext = PCylin -> U.Pl;   /* And stick it into cylin polygons. */
  274.     PCylin -> U.Pl = PBase1;
  275.     UpdatePolyPlane(PBase2, Pt);      /* Update base polygon plane equation. */
  276.     SET_CONVEX_POLY(PBase2);               /* Mark it as convex polygon. */
  277.     PBase2 -> Pnext = PCylin -> U.Pl;   /* And stick it into cylin polygons. */
  278.     PCylin -> U.Pl = PBase2;
  279.  
  280.     SET_OBJECT_COLOR(PCylin, PrimColor);       /* Set its default color. */
  281.  
  282.     return PCylin;
  283. }
  284.  
  285. /*****************************************************************************
  286. *   Routine to create a SPHERE geometric object defined by Center - sphere   *
  287. * 3d center point, and R its radius.                         *
  288. *   Note polygons on the poles are triangles, while others are rectangles    *
  289. *   Teta is horizontal circle angle, Fee is the vertical (spherical coords.) *
  290. *   The vertical axes here is assumed to be the Z axes.                 *
  291. *****************************************************************************/
  292. struct ObjectStruct * GenSPHEREObject(VectorType Center, RealType *R)
  293. {
  294.     int i, j, Resolution;
  295.     RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
  296.     CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
  297.     PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt;
  298.     struct ObjectStruct *PSphere;
  299.  
  300.     Resolution = GetResolution();     /* Get refinement factor of object. */
  301.  
  302.     PSphere = GenGeomObject("", NULL, NULL);/* Gen the SPHERE object itself: */
  303.  
  304.     TetaAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  305.     FeeAngleStep = M_PI * 2.0 / Resolution;    /* Runs from -PI/2 yo +PI/2. */
  306.  
  307.     /* Generate the lowest (south pole) triangular polygons: */
  308.     FeeAngle = (-M_PI/2.0) + FeeAngleStep; /* First circle above south pole. */
  309.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  310.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  311.     PT_COPY(LastCirclePt, Center);        /* Calculate the south pole. */
  312.     LastCirclePt[2] -= (*R);
  313.     PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
  314.     CircleLastPt[0] += CosFeeAngle1;
  315.     CircleLastPt[2] += SinFeeAngle1;
  316.  
  317.     for (i=1; i<=Resolution; i++) {   /* Pass the whole (horizontal) circle. */
  318.     TetaAngle = TetaAngleStep * i;         /* Prevent from additive error. */
  319.  
  320.     PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
  321.     CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  322.     CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  323.     CirclePt[2] += SinFeeAngle1;
  324.  
  325.     PSphere -> U.Pl = GenPolygon3Vrtx(LastCirclePt, CircleLastPt,
  326.                     CirclePt, Center, PSphere -> U.Pl);
  327.  
  328.     PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
  329.     }
  330.  
  331.     /* Generate the middle rectangular polygons: */
  332.     for (i=1; i<Resolution/2-1; i++) { /* For all the horizontal circles do. */
  333.     FeeAngle = (-M_PI/2.0) + FeeAngleStep * i;
  334.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  335.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  336.     FeeAngle = (-M_PI/2.0) + FeeAngleStep * (i + 1);
  337.     CosFeeAngle2 = cos(FeeAngle) * (*R);
  338.     SinFeeAngle2 = sin(FeeAngle) * (*R);
  339.     PT_COPY(CircleLastPt, Center);/* Calc. last point on current circle. */
  340.     CircleLastPt[0] += CosFeeAngle2;
  341.     CircleLastPt[2] += SinFeeAngle2;
  342.     PT_COPY(LastCircleLastPt, Center);/* Calc. last point on last circle.*/
  343.     LastCircleLastPt[0] += CosFeeAngle1;
  344.     LastCircleLastPt[2] += SinFeeAngle1;
  345.  
  346.     for (j=1; j<=Resolution; j++) {/* Pass the whole (horizontal) circle.*/
  347.         TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */
  348.  
  349.         PT_COPY(CirclePt, Center);/* Calc. current pt on current circle. */
  350.         CirclePt[0] += cos(TetaAngle) * CosFeeAngle2;
  351.         CirclePt[1] += sin(TetaAngle) * CosFeeAngle2;
  352.         CirclePt[2] += SinFeeAngle2;
  353.         PT_COPY(LastCirclePt, Center);/* Calc. current pt on last circle.*/
  354.         LastCirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  355.         LastCirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  356.         LastCirclePt[2] += SinFeeAngle1;
  357.  
  358.         PSphere -> U.Pl = GenPolygon4Vrtx(LastCirclePt, LastCircleLastPt,
  359.             CircleLastPt, CirclePt, Center, PSphere -> U.Pl);
  360.  
  361.         PT_COPY(CircleLastPt, CirclePt);          /* Save pt in last pt. */
  362.         PT_COPY(LastCircleLastPt, LastCirclePt);
  363.     }
  364.     }
  365.  
  366.     /* Generate the upper most (north pole) triangular polygons: */
  367.     FeeAngle = (M_PI/2.0) - FeeAngleStep;  /* First circle below north pole. */
  368.     CosFeeAngle1 = cos(FeeAngle) * (*R);
  369.     SinFeeAngle1 = sin(FeeAngle) * (*R);
  370.     PT_COPY(LastCirclePt, Center);        /* Calculate the north pole. */
  371.     LastCirclePt[2] += (*R);
  372.     PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
  373.     CircleLastPt[0] += CosFeeAngle1;
  374.     CircleLastPt[2] += SinFeeAngle1;
  375.  
  376.     for (i=1; i<=Resolution; i++) {   /* Pass the whole (horizontal) circle. */
  377.     TetaAngle = TetaAngleStep * i;         /* Prevent from additive error. */
  378.  
  379.     PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
  380.     CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
  381.     CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
  382.     CirclePt[2] += SinFeeAngle1;
  383.  
  384.     PSphere -> U.Pl = GenPolygon3Vrtx(LastCirclePt, CirclePt, CircleLastPt,
  385.                     Center, PSphere -> U.Pl);
  386.  
  387.     PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
  388.     }
  389.  
  390.     SET_OBJECT_COLOR(PSphere, PrimColor);       /* Set its default color. */
  391.  
  392.     return PSphere;
  393. }
  394.  
  395. /*****************************************************************************
  396. *   Routine to create a TORUS geometric object defined by Center - torus 3d  *
  397. * center point, the main torus plane normal Normal, major radius Rmajor and  *
  398. * minor radius Rminor (Tube radius).                         *
  399. *   Teta runs on the major circle, Fee on the minor one (Before Mat trans.): *
  400. * X = (Rmajor + Rminor * cos(Fee)) * cos(Teta)                     *
  401. * Y = (Rmajor + Rminor * cos(Fee)) * sin(Teta)                     *
  402. * Z = Rminor * sin(Fee)                                 *
  403. *****************************************************************************/
  404. struct ObjectStruct * GenTORUSObject(VectorType Center, VectorType Normal,
  405.                     RealType *Rmajor, RealType *Rminor)
  406. {
  407.     int i, j, Resolution;
  408.     RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
  409.     CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
  410.     PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt, InPt;
  411.     MatrixType Mat;
  412.     struct ObjectStruct *PTorus;
  413.  
  414.     Resolution = GetResolution();     /* Get refinement factor of object. */
  415.  
  416.     GenTransformMatrix(Mat, Center, Normal, 1.0); /* Trans from unit circle. */
  417.  
  418.     PTorus = GenGeomObject("", NULL, NULL); /* Gen. the Torus object itself: */
  419.  
  420.     TetaAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  421.     FeeAngleStep = M_PI * 2.0 / Resolution;         /* Runs from 0 to 2*PI. */
  422.  
  423.     for (i=1; i<=Resolution; i++) {
  424.     FeeAngle = FeeAngleStep * (i - 1);
  425.     CosFeeAngle1 = cos(FeeAngle) * (*Rminor);
  426.     SinFeeAngle1 = sin(FeeAngle) * (*Rminor);
  427.     FeeAngle = FeeAngleStep * i;
  428.     CosFeeAngle2 = cos(FeeAngle) * (*Rminor);
  429.     SinFeeAngle2 = sin(FeeAngle) * (*Rminor);
  430.     LastCircleLastPt[0] = (*Rmajor) + CosFeeAngle1;
  431.     LastCircleLastPt[1] = 0.0;
  432.     LastCircleLastPt[2] = SinFeeAngle1;
  433.     MatMultVecby4by4(LastCircleLastPt, LastCircleLastPt, Mat);
  434.     LastCirclePt[0] = (*Rmajor) + CosFeeAngle2;
  435.     LastCirclePt[1] = 0.0;
  436.     LastCirclePt[2] = SinFeeAngle2;
  437.     MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
  438.  
  439.     for (j=1; j<=Resolution; j++) {
  440.         TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */
  441.  
  442.         CircleLastPt[0] = ((*Rmajor) + CosFeeAngle1) * cos(TetaAngle);
  443.         CircleLastPt[1] = ((*Rmajor) + CosFeeAngle1) * sin(TetaAngle);
  444.         CircleLastPt[2] = SinFeeAngle1;
  445.         MatMultVecby4by4(CircleLastPt, CircleLastPt, Mat);
  446.         CirclePt[0] = ((*Rmajor) + CosFeeAngle2) * cos(TetaAngle);
  447.         CirclePt[1] = ((*Rmajor) + CosFeeAngle2) * sin(TetaAngle);
  448.         CirclePt[2] = SinFeeAngle2;
  449.         MatMultVecby4by4(CirclePt, CirclePt, Mat);
  450.         /* Point inside the object relative to this polygon: */
  451.         InPt[0] = (*Rmajor) * cos(TetaAngle);
  452.         InPt[1] = (*Rmajor) * sin(TetaAngle);
  453.         InPt[2] = 0.0;
  454.         MatMultVecby4by4(InPt, InPt, Mat);
  455.  
  456.         PTorus -> U.Pl = GenPolygon4Vrtx(CirclePt, CircleLastPt,
  457.         LastCircleLastPt, LastCirclePt, InPt, PTorus -> U.Pl);
  458.  
  459.         PT_COPY(LastCirclePt, CirclePt);          /* Save pt in last pt. */
  460.         PT_COPY(LastCircleLastPt, CircleLastPt);
  461.     }
  462.     }
  463.  
  464.     SET_OBJECT_COLOR(PTorus, PrimColor);       /* Set its default color. */
  465.  
  466.     return PTorus;
  467. }
  468.  
  469. /*****************************************************************************
  470. *   Routine to create a PLANE geometric object defined by the normal N and   *
  471. * the translation vector T. The object is a FINITE plane (a circle of         *
  472. * RESOLUTION point in it...) and its radius is equal to R.             *
  473. *   The normal direction is assumed to point to the inside of the object.    *
  474. *****************************************************************************/
  475. struct ObjectStruct * GenPLANEObject(VectorType N, VectorType T, RealType *R)
  476. {
  477.     int i, Resolution;
  478.     RealType Angle, AngleStep;
  479.     PointType CirclePt;
  480.     MatrixType Mat;
  481.     struct VertexStruct *V;
  482.     struct PolygonStruct *PCirc;
  483.     struct ObjectStruct *PPlane;
  484.  
  485.     Resolution = GetResolution();     /* Get refinement factor of object. */
  486.  
  487.     GenTransformMatrix(Mat, T, N, *R);          /* Transform from unit circle. */
  488.  
  489.     PPlane = GenGeomObject("", NULL, NULL); /* Gen. the PLANE object itself: */
  490.     PCirc = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), NULL);
  491.     PPlane -> U.Pl = PCirc;
  492.  
  493.     CirclePt[0] = 1.0;            /* First point is allways Angle = 0. */
  494.     CirclePt[1] = 0.0;
  495.     CirclePt[2] = 0.0;
  496.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  497.     PT_COPY(V -> Pt, CirclePt);           /* Update first pt in circle polygon. */
  498.  
  499.     AngleStep = M_PI * 2 / Resolution;
  500.  
  501.     for (i=1; i<=Resolution; i++) {          /* Pass the whole base circle. */
  502.     Angle = AngleStep * i;             /* Prevent from additive error. */
  503.  
  504.     CirclePt[0] = cos(Angle);
  505.     CirclePt[1] = sin(Angle);
  506.     CirclePt[2] = 0.0;
  507.  
  508.     MatMultVecby4by4(CirclePt, CirclePt, Mat);
  509.  
  510.     /* And add this vertices to the two cylinder bases: */
  511.     if (i == Resolution) {           /* Its last point - make it circular. */
  512.         V -> Pnext = PCirc -> V;
  513.     }
  514.     else {
  515.         V -> Pnext = AllocVertex(0, 0, NULL, NULL);
  516.         V = V -> Pnext;
  517.         PT_COPY(V -> Pt, CirclePt);
  518.     }
  519.     }
  520.  
  521.     PT_ADD(CirclePt, CirclePt, N);   /* Make a point "IN" the circle object. */
  522.     UpdatePolyPlane(PCirc, CirclePt); /* Update base polygon plane equation. */
  523.     SET_CONVEX_POLY(PCirc);               /* Mark it as convex polygon. */
  524.  
  525.     SET_OBJECT_COLOR(PPlane, PrimColor);       /* Set its default color. */
  526.  
  527.     return PPlane;
  528. }
  529.  
  530. /*****************************************************************************
  531. *   Routine to create a POLYgon directly from its specifed vertices.         *
  532. *   The validity of list elements is tested to make sure they are vectors as *
  533. * nobody else checked it till now (lists can hold any object type!).         *
  534. *   No test is made to make sure all vertices are on one plane, and that no  *
  535. * two vertices are similar.                             *
  536. *****************************************************************************/
  537. struct ObjectStruct *GenPOLYObject(ObjectStruct *PObjList)
  538. {
  539.     int i, NumVertices = 0;
  540.     VertexStruct *V, *VHead, *VTail;
  541.     PolygonStruct *PPoly;
  542.     ObjectStruct *PObj, *PObjPoly;
  543.  
  544.     if (!IS_OLST_OBJ(PObjList))
  545.     FatalError("GenPOLYObject: Not object list object!\n");
  546.  
  547.     i = 0;
  548.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  549.     if (!IS_VEC_OBJ(PObj)) {
  550.         WndwInputWindowPutStr("POLY: None vector object found in list, empty object result.", RED);
  551.         return NULL;
  552.     }
  553.     NumVertices++;
  554.     }
  555.  
  556.     if (NumVertices < 3) {
  557.     WndwInputWindowPutStr("POLY: Less than 3 vertices, empty object result.", RED);
  558.     return NULL;
  559.     }
  560.  
  561.     PPoly = AllocPolygon(0, 0, VHead = NULL, NULL);
  562.     i = 0;
  563.     while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
  564.     V = AllocVertex(0, 0, NULL, NULL);
  565.     PT_COPY(V -> Pt, PObj -> U.Vec);
  566.     if (VHead == NULL)
  567.     {
  568.         PPoly -> V = VHead = VTail = V;
  569.     }
  570.     else
  571.     {
  572.         VTail -> Pnext = V;
  573.         VTail = V;
  574.     }
  575.     }
  576.     VTail -> Pnext = VHead;              /* Close the vertex list loop. */
  577.  
  578.     CGPlaneFrom3Points(PPoly -> Plane, VHead -> Pt, VHead -> Pnext -> Pt,
  579.                         VHead -> Pnext -> Pnext -> Pt);
  580.  
  581.     PObjPoly = GenGeomObject("", NULL, NULL);/* Gen. the POLY object itself: */
  582.     PObjPoly -> U.Pl = PPoly;
  583.     SET_OBJECT_COLOR(PObjPoly, PrimColor);       /* Set its default color. */
  584.  
  585.     return PObjPoly;
  586. }
  587.  
  588. /*****************************************************************************
  589. *   Routine to a cross section polygon, by interactively allowing the user   *
  590. * to create it. This polygon is returned as a one polygon object. This       *
  591. * polygon is mainly for Surface of Revolution (SURFREV) and extrusion        *
  592. * (EXTRUDE) operations, although it may be used as an open object by itself. *
  593. *****************************************************************************/
  594. struct ObjectStruct * GenCROSSECObject(ObjectStruct *PObj)
  595. {
  596.     if (PObj && !IS_GEOM_OBJ(PObj))
  597.     FatalError("CrossSec: operation on none geometric object\n");
  598.  
  599.     WndwInputWindowPutStr("GenCrossSecObject not implemented", RED);
  600.  
  601.     return NULL;
  602. }
  603.  
  604. /*****************************************************************************
  605. *   Routine to a create surface of revolution by rotating the given cross    *
  606. * section along Z axes. The polygon plane must not be perpendicular to Z.    *
  607. *****************************************************************************/
  608. struct ObjectStruct * GenSURFREVObject(ObjectStruct *Cross)
  609. {
  610.     int Resolution, i;
  611.     MatrixType Mat;               /* Rotation Matrix around Z axes. */
  612.     VertexStruct *V1, *V1Head, *V2, *V2Head, *VIn, *VInHead;
  613.     PolygonStruct *Pl1, *Pl2, *PlIn, *PlNew = NULL;
  614.     ObjectStruct *PSurfRev;
  615.  
  616.     if (!IS_GEOM_OBJ(Cross))
  617.     FatalError("SurfRev: cross section is not geometric object\n");
  618.  
  619.     if (APX_EQ(Cross -> U.Pl -> Plane[0], 0.0) &&
  620.     APX_EQ(Cross -> U.Pl -> Plane[1], 0.0)) {
  621.     WndwInputWindowPutStr("SurfRev: cross-section perpendicular to Z. Empty object result", RED);
  622.     return NULL;
  623.     }
  624.  
  625.     Pl1 = AllocPolygon(0, 0, V1Head = CopyVList(Cross -> U.Pl -> V), NULL);
  626.     PLANE_COPY(Pl1 -> Plane, Cross -> U.Pl -> Plane);
  627.     Pl2 = AllocPolygon(0, 0, V2Head = CopyVList(Cross -> U.Pl -> V), NULL);
  628.     PLANE_COPY(Pl2 -> Plane, Cross -> U.Pl -> Plane);
  629.     PlIn = GenInsidePoly(Pl1);
  630.     VInHead = PlIn -> V;
  631.     Resolution = GetResolution();     /* Get refinement factor of object. */
  632.     MatGenMatRotZ1(2.0 * M_PI / Resolution, Mat);
  633.  
  634.     for (i=0; i<Resolution; i++)
  635.     {
  636.     V2 = V2Head;
  637.     do {
  638.         MatMultVecby4by4(V2 -> Pt, V2 -> Pt , Mat);
  639.         V2 = V2 -> Pnext;
  640.     }
  641.     while (V2 != NULL && V2 != V2Head);
  642.  
  643.     V1 = V1Head;
  644.     if (i < Resolution - 1) /* If this is the last loop use the original */
  645.          V2 = V2Head; /* polygon as we might accumulate error during the */
  646.     else V2 = Cross -> U.Pl -> V;   /* transformations along the circle. */
  647.     VIn = VInHead;
  648.     do
  649.     {
  650.         PlNew = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
  651.                     V2 -> Pnext -> Pt, V2 -> Pt,
  652.                     VIn -> Pt, PlNew);
  653.  
  654.         VIn = VIn -> Pnext;
  655.         V1 = V1 -> Pnext;
  656.         V2 = V2 -> Pnext;
  657.     }
  658.     while (V1 -> Pnext != NULL && V1 != V1Head);
  659.  
  660.     V1 = V1Head;
  661.     do {
  662.         MatMultVecby4by4(V1 -> Pt, V1 -> Pt , Mat);
  663.         V1 = V1 -> Pnext;
  664.     }
  665.     while (V1 != NULL && V1 != V1Head);
  666.     VIn = VInHead;
  667.     do {
  668.         MatMultVecby4by4(VIn -> Pt, VIn -> Pt , Mat);
  669.         VIn = VIn -> Pnext;
  670.     }
  671.     while (VIn != NULL && VIn != VInHead);
  672.     }
  673.  
  674.     MyFree((char *) PlIn, POLYGON_TYPE);
  675.     MyFree((char *) Pl1, POLYGON_TYPE);
  676.     MyFree((char *) Pl2, POLYGON_TYPE);
  677.  
  678.     PSurfRev = GenGeomObject("", NULL, NULL);
  679.     PSurfRev -> U.Pl = PlNew;
  680.     SET_OBJECT_COLOR(PSurfRev, PrimColor);       /* Set its default color. */
  681.  
  682.     return PSurfRev;
  683. }
  684.  
  685. /*****************************************************************************
  686. *   Routine to a create surface of extrusion out of the given cross section  *
  687. * and the given direction. A NULL object will be returned, if the direction  *
  688. * vector is in the cross section plane.                         *
  689. *****************************************************************************/
  690. struct ObjectStruct * GenEXTRUDEObject(ObjectStruct *Cross, VectorType Dir)
  691. {
  692.     int i;
  693.     RealType R;
  694.     VertexStruct *V1, *V1Head, *V2, *VIn;
  695.     PolygonStruct *PBase1, *PBase2, *Pl, *PlIn;
  696.     ObjectStruct *PExtrude;
  697.  
  698.     if (!IS_GEOM_OBJ(Cross))
  699.     FatalError("Extrude: cross section is not geometric object\n");
  700.  
  701.     R = DOT_PROD(Cross -> U.Pl -> Plane, Dir);
  702.     if (APX_EQ(R, 0.0)) {
  703.     WndwInputWindowPutStr("Extrusion direction in cross-section plane. Empty object result", RED);
  704.     return NULL;
  705.     }
  706.  
  707.     /* Prepare the two bases (update their plane normal to point INDISE): */
  708.     PBase1 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl -> V), NULL);
  709.     Pl = PBase2 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl -> V), PBase1);
  710.     V1 = V1Head = PBase2 -> V;
  711.     do {
  712.     PT_ADD(V1 -> Pt, Dir, V1 -> Pt);
  713.     V1 = V1 -> Pnext;
  714.     }
  715.     while (V1 != NULL && V1 != V1Head);
  716.     if (R > 0.0) {
  717.     PLANE_COPY(PBase1 -> Plane, Cross -> U.Pl -> Plane);
  718.     for (i=0; i<3; i++) PBase2 -> Plane[i] = (-Cross -> U.Pl -> Plane[i]);
  719.     PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
  720.     }
  721.     else {
  722.     for (i=0; i<4; i++) PBase1 -> Plane[i] = (-Cross -> U.Pl -> Plane[i]);
  723.     PLANE_COPY(PBase2 -> Plane, Cross -> U.Pl -> Plane);
  724.     PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
  725.     }
  726.  
  727.     /* Now generate all the 4 corner polygon between the two bases: */
  728.     V1 = V1Head = PBase1 -> V;
  729.     V2 = PBase2 -> V;
  730.     PlIn = GenInsidePoly(PBase1);
  731.     VIn = PlIn -> V;
  732.     do {
  733.     Pl = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
  734.                  V2 -> Pnext -> Pt, V2 -> Pt, VIn -> Pt, Pl);
  735.     VIn = VIn -> Pnext;
  736.     V1 = V1 -> Pnext;
  737.     V2 = V2 -> Pnext;
  738.     }
  739.     while (V1 -> Pnext != NULL && V1 != V1Head);
  740.  
  741.     MyFree((char *) PlIn, POLYGON_TYPE);
  742.  
  743.     PExtrude = GenGeomObject("", NULL, NULL);
  744.     PExtrude -> U.Pl = Pl;
  745.     SET_OBJECT_COLOR(PExtrude, PrimColor);       /* Set its default color. */
  746.  
  747.     return PExtrude;
  748. }
  749.  
  750. /*****************************************************************************
  751. *   Routine to create a pseudo polygon out of a given polygon such that each *
  752. * vertex Vi is in the inside side of the corresponding edge ViVi+1 in the    *
  753. * given polygon. Used in polygon generation for EXTRUDE/SURFREV operations.  *
  754. *****************************************************************************/
  755. static struct PolygonStruct *GenInsidePoly(PolygonStruct *Pl)
  756. {
  757.     int Axes;
  758.     RealType Dx, Dy;
  759.     PointType Pt;
  760.     MatrixType Mat;
  761.     PolygonStruct *PlIn;
  762.     VertexStruct *VHead, *V, *Vnext, *VInHead, *VIn = NULL;
  763.  
  764.     PlIn = AllocPolygon(0, 0, VInHead = NULL, NULL);
  765.  
  766.     /* Generate transformation matrix to bring polygon to a XY parallel      */
  767.     /* plane, and transform a copy of the polygon to that plane.         */
  768.     GenRotateMatrix(Mat, Pl -> Plane);
  769.     VHead = V = CopyVList(Pl -> V);      /* We dont want to modify original! */
  770.     Pl = AllocPolygon(0, 0, VHead, NULL);
  771.     do {
  772.     MatMultVecby4by4(V -> Pt, V -> Pt, Mat);
  773.     V = V -> Pnext;
  774.     }
  775.     while (V != NULL && V != VHead);
  776.  
  777.     V = VHead;
  778.     do {
  779.     Vnext = V -> Pnext;
  780.     Dx = ABS(V -> Pt[0] - Vnext -> Pt[0]);
  781.     Dy = ABS(V -> Pt[1] - Vnext -> Pt[1]);
  782.     Pt[0] = (V -> Pt[0] + Vnext -> Pt[0]) / 2.0;/* Prepare middle point. */
  783.     Pt[1] = (V -> Pt[1] + Vnext -> Pt[1]) / 2.0;
  784.     Pt[2] = V -> Pt[2];
  785.     /* If Dx > Dy fire ray in +Y direction, otherwise in +X direction    */
  786.     /* and if number of intersections is even (excluding the given point */
  787.     /* itself) then that direction is the outside, otherwise, its inside.*/
  788.     Axes = (Dx > Dy ? 1 : 0);
  789.     if (CGPolygonRayInter(Pl, Pt, Axes) % 2 == 0)
  790.     {
  791.         /* The amount we move along Axes is not of a big meaning as long */
  792.         /* as it is not zero, so MAX(Dx, Dy) guarantee non zero value... */
  793.         Pt[Axes] -= MAX(Dx, Dy);
  794.     }
  795.     else {
  796.         Pt[Axes] += MAX(Dx, Dy);
  797.     }
  798.  
  799.     /* Now Pt holds point which is in the inside part of vertex V, Vnext.*/
  800.     /* Put it in the pseudo inside polygon PlIn:                 */
  801.     if (VInHead) {
  802.         VIn -> Pnext = AllocVertex(0, 0, NULL, NULL);
  803.         VIn = VIn -> Pnext;
  804.     }
  805.     else {
  806.         PlIn -> V = VInHead = VIn = AllocVertex(0, 0, NULL, NULL);
  807.     }
  808.     PT_COPY(VIn ->Pt, Pt);
  809.  
  810.     V = Vnext;
  811.     }
  812.     while (V != NULL && V != VHead);
  813.     VIn -> Pnext = VInHead;
  814.  
  815.     MyFree((char *) Pl, POLYGON_TYPE);/* Free copied (and trans.) vrtx list. */
  816.  
  817.     /* Transform PlIn to the plane where original Pl is... */
  818.     if (!MatInverseMatrix(Mat, Mat))         /* Find the inverse matrix. */
  819.     FatalError("GenInsidePoly: Inverse matrix does not exits");
  820.     VIn = VInHead;
  821.     do {
  822.     MatMultVecby4by4(VIn -> Pt, VIn -> Pt, Mat);
  823.     VIn = VIn -> Pnext;
  824.     }
  825.     while (VIn != NULL && VIn != VInHead);
  826.  
  827.     return PlIn;
  828. }
  829.  
  830. /*****************************************************************************
  831. *   Routine to create a polygon out of a list of 4 vertices V1/2/3/4         *
  832. * The fifth vertex is inside (actually, this is not true, as this point will *
  833. * be in the positive part of the plane, which only locally in the object...) *
  834. * the object, so the polygon normal direction can be evaluated uniquely.     *
  835. *  No test is made to make sure the 4 points are co-planar...             *
  836. *****************************************************************************/
  837. static struct PolygonStruct *GenPolygon4Vrtx(VectorType V1, VectorType V2,
  838.     VectorType V3, VectorType V4, VectorType Vin, PolygonStruct *Pnext)
  839. {
  840.     struct PolygonStruct *PPoly;
  841.     struct VertexStruct *V;
  842.  
  843.     PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
  844.     PT_COPY(V -> Pt, V1);
  845.  
  846.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  847.     PT_COPY(V -> Pt, V2);
  848.  
  849.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  850.     PT_COPY(V -> Pt, V3);
  851.  
  852.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  853.     PT_COPY(V -> Pt, V4);
  854.  
  855.     V -> Pnext = PPoly -> V;           /* Make the Vertex list circular. */
  856.  
  857.     UpdatePolyPlane(PPoly, Vin);           /* Update plane equation. */
  858.  
  859.     SET_CONVEX_POLY(PPoly);               /* Mark it as convex polygon. */
  860.  
  861.     return PPoly;
  862. }
  863.  
  864. /*****************************************************************************
  865. *   Routine to create a polygon out of a list of 3 vertices V1/2/3         *
  866. * The forth vertex is inside (actually, this is not true, as this point will *
  867. * be in the positive part of the plane, which only locally in the object...) *
  868. * the object, so the polygon normal direction can be evaluated uniquely.     *
  869. *****************************************************************************/
  870. static struct PolygonStruct *GenPolygon3Vrtx(VectorType V1, VectorType V2,
  871.             VectorType V3, VectorType Vin, PolygonStruct *Pnext)
  872. {
  873.     struct PolygonStruct *PPoly;
  874.     struct VertexStruct *V;
  875.  
  876.     PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
  877.     PT_COPY(V -> Pt, V1);
  878.  
  879.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  880.     PT_COPY(V -> Pt, V2);
  881.  
  882.     V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
  883.     PT_COPY(V -> Pt, V3);
  884.  
  885.     V -> Pnext = PPoly -> V;           /* Make the Vertex list circular. */
  886.  
  887.     UpdatePolyPlane(PPoly, Vin);           /* Update plane equation. */
  888.  
  889.     SET_CONVEX_POLY(PPoly);               /* Mark it as convex polygon. */
  890.  
  891.     return PPoly;
  892. }
  893.  
  894. /*****************************************************************************
  895. *   Routine to preper transformation martix to do the following (in this     *
  896. * order): scale by Scale, rotate such that the Z axis is in direction Dir    *
  897. * and then translate by Trans.                             *
  898. *    Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is *
  899. * used to form the second line (the first 3 lines set the rotation), and     *
  900. * finally Scale is used to scale first 3 lines/columns to the needed scale:  *
  901. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord *
  902. *                |  Bx  By  Bz  0 |  system into T, N & B as required and    *
  903. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |  then translate it to C. T, N, B are     *
  904. *                |  Cx  Cy  Cz  1 |  scaled by Scale.                 *
  905. * N is exactly Dir (unit vec) but we got freedom on T & B - T & B must be on *
  906. * a plane perpendicular to N and perpendicular between them but thats all!   *
  907. * T is therefore selected using this (heuristic ?) algorithm:             *
  908. * Let P be the axis of which the absolute N coefficient is the smallest.     *
  909. * Let B be (N cross P) and T be (B cross N).                     *
  910. *****************************************************************************/
  911. static void GenTransformMatrix(MatrixType Mat, VectorType Trans,
  912.                         VectorType Dir, RealType Scale)
  913. {
  914.     int i, j;
  915.     RealType R;
  916.     VectorType DirN, T, B, P;
  917.     MatrixType TempMat;
  918.  
  919.     PT_COPY(DirN, Dir);
  920.     PT_NORMALIZE(DirN);
  921.     PT_CLEAR(P);
  922.     for (i=1, j=0, R = ABS(DirN[0]); i<3; i++) if (R > ABS(DirN[i])) {
  923.     R = DirN[i];
  924.     j = i;
  925.     }
  926.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  927.  
  928.     VecCrossProd(B, DirN, P);                  /* Calc the bi-normal. */
  929.     VecCrossProd(T, B, DirN);                /* Calc the tangent. */
  930.  
  931.     MatGenUnitMat(Mat);
  932.     for (i=0; i<3; i++) {
  933.     Mat[0][i] = T[i];
  934.     Mat[1][i] = B[i];
  935.     Mat[2][i] = DirN[i];
  936.     }
  937.     MatGenMatScale(Scale, Scale, Scale, TempMat);
  938.     MatMultTwo4by4(Mat, TempMat, Mat);
  939.  
  940.     MatGenMatTrans(Trans[0], Trans[1], Trans[2], TempMat);
  941.     MatMultTwo4by4(Mat, Mat, TempMat);
  942. }
  943.  
  944. /*****************************************************************************
  945. *   Routine to update the Plane equation of the given polygon such that the  *
  946. * Vin vertex will be in the positive side of it.                 *
  947. *   It is assumed the polygon has at list 3 points...                 *
  948. *****************************************************************************/
  949. void UpdatePolyPlane(PolygonStruct *PPoly, VectorType Vin)
  950. {
  951.     int i;
  952.     VectorType V1, V2;
  953.     RealType Len;
  954.     struct VertexStruct *V;
  955.  
  956.     V = PPoly -> V;    PT_SUB(V1, V -> Pt, V -> Pnext -> Pt);
  957.     V = V -> Pnext;    PT_SUB(V2, V -> Pt, V -> Pnext -> Pt);
  958.  
  959.     VecCrossProd(PPoly -> Plane, V1, V2);           /* Find Plane Normal. */
  960.     /* Normalize the plane such that the normal has length of 1: */
  961.     Len = PT_LENGTH(PPoly -> Plane);
  962.     for (i=0; i<3; i++) PPoly -> Plane[i] /= Len;
  963.  
  964.     PPoly -> Plane[3] = (-DOT_PROD(PPoly -> Plane, PPoly -> V -> Pt));
  965.  
  966.     if (DOT_PROD(PPoly -> Plane, Vin) + PPoly -> Plane[3] < 0) {
  967.     /* Flip plane normal and reverse the vertex list. */
  968.     ReverseVrtxList(PPoly);
  969.     for (i=0; i<4; i++) PPoly -> Plane[i] = (-PPoly -> Plane[i]);
  970.     }
  971. }
  972.