home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 3 / goldfish_volume_3.bin / files / gfx / 3d / irit / geom_lib / geomat3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-24  |  51.7 KB  |  1,163 lines
  1. /******************************************************************************
  2. * GeoMat3d.c - Trans. Matrices , Vector computation, and Comp.geom.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. #include <math.h>
  8. #include <stdio.h>
  9. #include "irit_sm.h"
  10. #include "iritprsr.h"
  11. #include "allocate.h"
  12. #include "convex.h"
  13. #include "geomat3d.h"
  14.  
  15. /* #define DEBUG1       Print information on entry and exit of routines. */
  16.  
  17. static int CGPointRayRelation(PointType Pt, PointType PtRay, int Axes);
  18.  
  19. /*****************************************************************************
  20. * DESCRIPTION:                                                               M
  21. *  Routine to copy one vector to another:                     M
  22. *                                                                            *
  23. * PARAMETERS:                                                                M
  24. *   Vdst:      Destination vector.                                           M
  25. *   Vsrc:      Source vector.                                                M
  26. *                                                                            *
  27. * RETURN VALUE:                                                              M
  28. *   void                                                                     M
  29. *                                                                            *
  30. * KEYWORDS:                                                                  M
  31. *   GMVecCopy, copy                                                          M
  32. *****************************************************************************/
  33. void GMVecCopy(VectorType Vdst, VectorType Vsrc)
  34. {
  35.     GEN_COPY(Vdst, Vsrc, sizeof(RealType) * 3);
  36. }
  37.  
  38. /*****************************************************************************
  39. * DESCRIPTION:                                                               M
  40. *  Routine to normalize the vector length to a unit size.                    M
  41. *                                                                            *
  42. * PARAMETERS:                                                                M
  43. *   V:        To normalize.                                                  M
  44. *                                                                            *
  45. * RETURN VALUE:                                                              M
  46. *   void                                                                     M
  47. *                                                                            *
  48. * KEYWORDS:                                                                  M
  49. *   GMVecNormalize, normalize                                                M
  50. *****************************************************************************/
  51. void GMVecNormalize(VectorType V)
  52. {
  53.     int i;
  54.     RealType Len;
  55.  
  56.     Len = GMVecLength(V);
  57.     if (Len > 0.0)
  58.         for (i = 0; i < 3; i++) {
  59.         V[i] /= Len;
  60.         if (ABS(V[i]) < IRIT_EPSILON)
  61.         V[i] = 0.0;
  62.     }
  63. }
  64.  
  65. /*****************************************************************************
  66. * DESCRIPTION:                                                               M
  67. *  Routine to compute the magnitude (length) of a given 3D vector:           M
  68. *                                                                            *
  69. * PARAMETERS:                                                                M
  70. *   V:        To compute its magnitude.                                      M
  71. *                                                                            *
  72. * RETURN VALUE:                                                              M
  73. *   RealType:   Magnitude of V.                                              M
  74. *                                                                            *
  75. * KEYWORDS:                                                                  M
  76. *   GMVecLength, magnitude                                                   M
  77. *****************************************************************************/
  78. RealType GMVecLength(VectorType V)
  79. {
  80.     return sqrt(SQR(V[0]) + SQR(V[1]) + SQR(V[2]));
  81. }
  82.  
  83. /*****************************************************************************
  84. * DESCRIPTION:                                                               M
  85. *  Routine to compute the cross product of two vectors.                      M
  86. * Note Vres may be the same as V1 or V2.                                     M
  87. *                                                                            *
  88. * PARAMETERS:                                                                M
  89. *   Vres:    Result of cross product                                         M
  90. *   V1, V2:  Two vectors of the cross product.                               M
  91. *                                                                            *
  92. * RETURN VALUE:                                                              M
  93. *   void                                                                     M
  94. *                                                                            *
  95. * KEYWORDS:                                                                  M
  96. *   GMVecCrossProd, cross prod                                               M
  97. *****************************************************************************/
  98. void GMVecCrossProd(VectorType Vres, VectorType V1, VectorType V2)
  99. {
  100.     VectorType Vtemp;
  101.  
  102.     Vtemp[0] = V1[1] * V2[2] - V2[1] * V1[2];
  103.     Vtemp[1] = V1[2] * V2[0] - V2[2] * V1[0];
  104.     Vtemp[2] = V1[0] * V2[1] - V2[0] * V1[1];
  105.  
  106.     GMVecCopy(Vres, Vtemp);
  107. }
  108.  
  109. /*****************************************************************************
  110. * DESCRIPTION:                                                               M
  111. *  Verifys the colinearity of three points.                                  M
  112. *                                                                            *
  113. * PARAMETERS:                                                                M
  114. *   Pt1, Pt2, Pt3: Three points to verify for colinearity.                   M
  115. *                                                                            *
  116. * RETURN VALUE:                                                              M
  117. *   int:      TRUE if colinear, FALSE otherwise.                             M
  118. *                                                                            *
  119. * KEYWORDS:                                                                  M
  120. *   GMColinear3Pts, colinearity                                              M
  121. *****************************************************************************/
  122. int GMColinear3Pts(PointType Pt1, PointType Pt2, PointType Pt3)
  123. {
  124.     VectorType V1, V2, V3;
  125.  
  126.     PT_SUB(V1, Pt1, Pt2);
  127.     PT_SUB(V2, Pt2, Pt3);
  128.  
  129.     if (PT_SQR_LENGTH(V1) < SQR(IRIT_EPSILON) ||
  130.     PT_SQR_LENGTH(V2) < SQR(IRIT_EPSILON))
  131.     return TRUE;
  132.  
  133.     GMVecCrossProd(V3, V1, V2);
  134.     return PT_SQR_LENGTH(V3) < SQR(IRIT_EPSILON);
  135. }
  136.  
  137. /*****************************************************************************
  138. * DESCRIPTION:                                                               M
  139. *  Routine to compute the dot product of two vectors.                        M
  140. *                                                                            M
  141. * PARAMETERS:                                                                M
  142. *   V1, V2:   Two vector to compute dot product of.                          M
  143. *                                                                            *
  144. * RETURN VALUE:                                                              M
  145. *   RealType:   Resulting dot product.                                       M
  146. *                                                                            *
  147. * KEYWORDS:                                                                  M
  148. *   GMVecDotProd, dot product                                                M
  149. *****************************************************************************/
  150. RealType GMVecDotProd(VectorType V1, VectorType V2)
  151. {
  152.     return  V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2];
  153. }
  154.  
  155. /*****************************************************************************
  156. * DESCRIPTION:                                                               M
  157. * Routine to generate rotation object around the X axis in Degree degrees:   M
  158. *                                                                            *
  159. * PARAMETERS:                                                                M
  160. *   Degree:     Amount of rotation, in degrees.                              M
  161. *                                                                            *
  162. * RETURN VALUE:                                                              M
  163. *   IPObjectStruct *:   A matrix object.                                     M
  164. *                                                                            *
  165. * KEYWORDS:                                                                  M
  166. *   GMGenMatObjectRotX, rotation, transformations                            M
  167. *****************************************************************************/
  168. IPObjectStruct *GMGenMatObjectRotX(RealType *Degree)
  169. {
  170.     MatrixType Mat;
  171.  
  172.     MatGenMatRotX1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  173.  
  174.     return GenMATObject(Mat);
  175. }
  176.  
  177. /*****************************************************************************
  178. * DESCRIPTION:                                                               M
  179. * Routine to generate rotation object around the Y axis in Degree degrees:   M
  180. *                                                                            *
  181. * PARAMETERS:                                                                M
  182. *   Degree:     Amount of rotation, in degrees.                              M
  183. *                                                                            *
  184. * RETURN VALUE:                                                              M
  185. *   IPObjectStruct *:   A matrix object.                                     M
  186. *                                                                            *
  187. * KEYWORDS:                                                                  M
  188. *   GMGenMatObjectRotY, rotation, transformations                            M
  189. *****************************************************************************/
  190. IPObjectStruct *GMGenMatObjectRotY(RealType *Degree)
  191. {
  192.     MatrixType Mat;
  193.  
  194.     MatGenMatRotY1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  195.  
  196.     return GenMATObject(Mat);
  197. }
  198.  
  199. /*****************************************************************************
  200. * DESCRIPTION:                                                               M
  201. * Routine to generate rotation object around the Z axis in Degree degrees:   M
  202. *                                                                            *
  203. * PARAMETERS:                                                                M
  204. *   Degree:     Amount of rotation, in degrees.                              M
  205. *                                                                            *
  206. * RETURN VALUE:                                                              M
  207. *   IPObjectStruct *:   A matrix object.                                     M
  208. *                                                                            *
  209. * KEYWORDS:                                                                  M
  210. *   GMGenMatObjectRotZ, rotation, transformations                            M
  211. *****************************************************************************/
  212. IPObjectStruct *GMGenMatObjectRotZ(RealType *Degree)
  213. {
  214.     MatrixType Mat;
  215.  
  216.     MatGenMatRotZ1(DEG2RAD(*Degree), Mat);    /* Generate the trans. matrix. */
  217.  
  218.     return GenMATObject(Mat);
  219. }
  220.  
  221. /*****************************************************************************
  222. * DESCRIPTION:                                                               M
  223. * Routine to generate a translation object.                                  M
  224. *                                                                            *
  225. * PARAMETERS:                                                                M
  226. *   Vec:     Amount of translation, in X, Y, and Z.                          M
  227. *                                                                            *
  228. * RETURN VALUE:                                                              M
  229. *   IPObjectStruct *:   A matrix object.                                     M
  230. *                                                                            *
  231. * KEYWORDS:                                                                  M
  232. *   GMGenMatObjectTrans, translation, transformations                        M
  233. *****************************************************************************/
  234. IPObjectStruct *GMGenMatObjectTrans(VectorType Vec)
  235. {
  236.     MatrixType Mat;
  237.  
  238.     /* Generate the transformation matrix */
  239.     MatGenMatTrans(Vec[0], Vec[1], Vec[2], Mat);
  240.  
  241.     return GenMATObject(Mat);
  242. }
  243.  
  244. /*****************************************************************************
  245. * DESCRIPTION:                                                               M
  246. * Routine to generate a scaling object.                                      M
  247. *                                                                            *
  248. * PARAMETERS:                                                                M
  249. *   Vec:     Amount of scaling, in X, Y, and Z.                              M
  250. *                                                                            *
  251. * RETURN VALUE:                                                              M
  252. *   IPObjectStruct *:   A matrix object.                                     M
  253. *                                                                            *
  254. * KEYWORDS:                                                                  M
  255. *   GMGenMatObjectScale, scaling, transformations                            M
  256. *****************************************************************************/
  257. IPObjectStruct *GMGenMatObjectScale(VectorType Vec)
  258. {
  259.     MatrixType Mat;
  260.  
  261.     /* Generate the transformation matrix */
  262.     MatGenMatScale(Vec[0], Vec[1], Vec[2], Mat);
  263.  
  264.     return GenMATObject(Mat);
  265. }
  266.  
  267. /*****************************************************************************
  268. * DESCRIPTION:                                                               M
  269. * Routine to transform an object according to the transformation matrix.     M
  270. *                                                                            *
  271. * PARAMETERS:                                                                M
  272. *   PObj:      Object to be transformed.                                     M
  273. *   Mat:       Transformation matrix.                                        M
  274. *                                                                            *
  275. * RETURN VALUE:                                                              M
  276. *   IPObjectStruct *:    Transformed object.                                 M
  277. *                                                                            *
  278. * KEYWORDS:                                                                  M
  279. *   GMTransformObject, transformations                                       M
  280. *****************************************************************************/
  281. IPObjectStruct *GMTransformObject(IPObjectStruct *PObj, MatrixType Mat)
  282. {
  283.     int i;
  284.     IPObjectStruct *NewPObj;
  285.  
  286.     if (IP_IS_POLY_OBJ(PObj)) {
  287.         int IsPolygon = IP_IS_POLYGON_OBJ(PObj);
  288.         VectorType Pin, PTemp;
  289.         IPPolygonStruct *Pl;
  290.         IPVertexStruct *V, *VFirst;
  291.  
  292.         NewPObj = CopyObject(NULL, PObj, FALSE);
  293.     Pl = NewPObj -> U.Pl;
  294.  
  295.     while (Pl != NULL) {
  296.         V = VFirst = Pl -> PVertex;
  297.         PT_ADD(Pin, V -> Coord, Pl -> Plane);  /* Prepare point IN side. */
  298.         MatMultVecby4by4(Pin, Pin, Mat); /* Transformed relative to new. */
  299.  
  300.         do {
  301.         if (IsPolygon)
  302.             PT_ADD(PTemp, V -> Coord, V -> Normal);
  303.  
  304.         MatMultVecby4by4(V -> Coord, V -> Coord, Mat);/* Update pos. */
  305.  
  306.         if (IsPolygon) {
  307.             MatMultVecby4by4(PTemp, PTemp, Mat);   /* Update normal. */
  308.             PT_SUB(V -> Normal, PTemp, V -> Coord);
  309.             PT_NORMALIZE(V -> Normal);
  310.         }
  311.  
  312.         V = V -> Pnext;
  313.         }
  314.         while (V != VFirst && V != NULL);           /* VList is circular! */
  315.  
  316.         if (IsPolygon)
  317.         IritPrsrUpdatePolyPlane2(Pl, Pin);    /* Update plane eqn. */
  318.  
  319.         Pl = Pl -> Pnext;
  320.     }
  321.     }
  322.     else if (IP_IS_CRV_OBJ(PObj)) {
  323.     CagdCrvStruct *Crv;
  324.  
  325.         NewPObj = CopyObject(NULL, PObj, FALSE);
  326.  
  327.     for (Crv = NewPObj -> U.Crvs; Crv != NULL; Crv = Crv -> Pnext)
  328.         CagdCrvMatTransform(Crv, Mat);
  329.     }
  330.     else if (IP_IS_SRF_OBJ(PObj)) {
  331.     CagdSrfStruct *Srf;
  332.  
  333.         NewPObj = CopyObject(NULL, PObj, FALSE);
  334.  
  335.     for (Srf = NewPObj -> U.Srfs; Srf != NULL; Srf = Srf -> Pnext)
  336.         CagdSrfMatTransform(Srf, Mat);
  337.     }
  338.     else if (IP_IS_TRIMSRF_OBJ(PObj)) {
  339.     TrimSrfStruct *TrimSrf;
  340.  
  341.         NewPObj = CopyObject(NULL, PObj, FALSE);
  342.  
  343.     for (TrimSrf = NewPObj -> U.TrimSrfs;
  344.          TrimSrf != NULL;
  345.          TrimSrf = TrimSrf -> Pnext)
  346.         TrimSrfMatTransform(TrimSrf, Mat);
  347.     }
  348.     else if (IP_IS_TRIVAR_OBJ(PObj)) {
  349.     TrivTVStruct *Trivar;
  350.  
  351.         NewPObj = CopyObject(NULL, PObj, FALSE);
  352.  
  353.     for (Trivar = NewPObj -> U.Trivars;
  354.          Trivar != NULL;
  355.          Trivar = Trivar -> Pnext)
  356.         TrivTVMatTransform(Trivar, Mat);
  357.     }
  358.     else if (IP_IS_POINT_OBJ(PObj)) {
  359.         NewPObj = CopyObject(NULL, PObj, FALSE);
  360.  
  361.     MatMultVecby4by4(NewPObj -> U.Pt, NewPObj -> U.Pt, Mat);
  362.     }
  363.     else if (IP_IS_CTLPT_OBJ(PObj)) {
  364.     IPObjectStruct
  365.         *TmpObj = IritPrsrCoerceObjectTo(PObj, CAGD_PT_P3_TYPE);
  366.     RealType R[4];
  367.  
  368.     PT_COPY(R, &TmpObj -> U.CtlPt.Coords[1]);
  369.     R[3] = TmpObj -> U.CtlPt.Coords[0];
  370.     MatMultVecby4by4(R, R, Mat);
  371.     TmpObj -> U.CtlPt.Coords[0] = R[3];
  372.     PT_COPY(&TmpObj -> U.CtlPt.Coords[1], R);
  373.  
  374.         NewPObj = IritPrsrCoerceObjectTo(TmpObj, PObj -> U.CtlPt.PtType);
  375.     IPFreeObject(TmpObj);
  376.     }
  377.     else if (IP_IS_VEC_OBJ(PObj)) {
  378.         NewPObj = CopyObject(NULL, PObj, FALSE);
  379.  
  380.     MatMultVecby4by4(NewPObj -> U.Vec, NewPObj -> U.Vec, Mat);
  381.     }
  382.     else if (IP_IS_OLST_OBJ(PObj)) {
  383.     IPObjectStruct *PObjTmp;
  384.  
  385.     NewPObj = IPAllocObject("", IP_OBJ_LIST_OBJ, NULL);
  386.  
  387.         for (i = 0; (PObjTmp = ListObjectGet(PObj, i)) != NULL; i++)
  388.             ListObjectInsert(NewPObj, i, GMTransformObject(PObjTmp, Mat));
  389.     ListObjectInsert(NewPObj, i, NULL);
  390.     }
  391.     else {
  392.         NewPObj = CopyObject(NULL, PObj, FALSE);
  393.     }
  394.  
  395.     /* Make the name with prefix "_". */
  396.     if (PObj -> Name[0] == '_')
  397.         strcpy(NewPObj -> Name, PObj -> Name);
  398.     else {
  399.     NewPObj -> Name[0] = '_';
  400.     strcpy(&NewPObj -> Name[1], PObj -> Name);
  401.     }
  402.  
  403.     return NewPObj;
  404. }
  405.  
  406. /*****************************************************************************
  407. * DESCRIPTION:                                                               M
  408. * Routine to transform an list of objects according to a transformation      M
  409. * matrix.                                                                 M
  410. *                                                                            *
  411. * PARAMETERS:                                                                M
  412. *   PObj:       Object list to transform.                                    M
  413. *   Mat:        Transformation matrix.                                       M
  414. *                                                                            *
  415. * RETURN VALUE:                                                              M
  416. *   IPObjectStruct *:   Transformed object list.                             M
  417. *                                                                            *
  418. * KEYWORDS:                                                                  M
  419. *   GMTransformObjectList, transformations                                   M
  420. *****************************************************************************/
  421. IPObjectStruct *GMTransformObjectList(IPObjectStruct *PObj, MatrixType Mat)
  422. {
  423.     IPObjectStruct
  424.     *PTailObj = NULL,
  425.     *PTransObj = NULL;
  426.  
  427.     for ( ; PObj != NULL; PObj = PObj -> Pnext) {
  428.     if (PTailObj == NULL)
  429.         PTailObj = PTransObj = GMTransformObject(PObj, Mat);
  430.     else {
  431.         PTailObj -> Pnext = GMTransformObject(PObj, Mat);
  432.         PTailObj = PTailObj -> Pnext;
  433.     }
  434.     }
  435.  
  436.     return PTransObj;
  437. }
  438.  
  439. /*****************************************************************************
  440. * DESCRIPTION:                                                               M
  441. *   Routine to compute the distance between two 3d points.                   M
  442. *                                                                            *
  443. * PARAMETERS:                                                                M
  444. *   P1, P2:   Two points to compute the distance between.                    M
  445. *                                                                            *
  446. * RETURN VALUE:                                                              M
  447. *   RealType:    Computed distance.                                          M
  448. *                                                                            *
  449. * KEYWORDS:                                                                  M
  450. *   CGDistPointPoint, point point distance                                   M
  451. *****************************************************************************/
  452. RealType CGDistPointPoint(PointType P1, PointType P2)
  453. {
  454.     RealType t;
  455.  
  456. #ifdef DEBUG1
  457.     printf("CGDistPointPoint entered\n");
  458. #endif /* DEBUG1 */
  459.  
  460.     t = sqrt(SQR(P2[0]-P1[0]) + SQR(P2[1]-P1[1]) + SQR(P2[2]-P1[2]));
  461.  
  462. #ifdef DEBUG1
  463.     printf("CGDistPointPoint exit\n");
  464. #endif /* DEBUG1 */
  465.  
  466.      return t;
  467. }
  468.  
  469. /*****************************************************************************
  470. * DESCRIPTION:                                                               M
  471. *   Routine to construct the plane from given 3 points. If two of the points M
  472. * are the same it returns FALSE, otherwise (successful) returns TRUE.         M
  473. *                                                                            *
  474. * PARAMETERS:                                                                M
  475. *   Plane:          To compute.                                              M
  476. *   Pt1, Pt2, Pt3:  Three points to fit a plane through.                     M
  477. *                                                                            *
  478. * RETURN VALUE:                                                              M
  479. *   int:     TRUE if successful, FALSE otherwise.                            M
  480. *                                                                            *
  481. * KEYWORDS:                                                                  M
  482. *   CGPlaneFrom3Points, plane                                                M
  483. *****************************************************************************/
  484. int CGPlaneFrom3Points(PlaneType Plane,
  485.                PointType Pt1,
  486.                PointType Pt2,
  487.                PointType Pt3)
  488. {
  489.     VectorType V1, V2;
  490.  
  491. #ifdef DEBUG1
  492.     printf("CGPlaneFrom3Points entered\n");
  493. #endif /* DEBUG1 */
  494.  
  495.     if (IRIT_PT_APX_EQ(Pt1, Pt2) ||
  496.     IRIT_PT_APX_EQ(Pt2, Pt3) ||
  497.     IRIT_PT_APX_EQ(Pt1, Pt3))
  498.     return FALSE;
  499.  
  500.     PT_SUB(V1, Pt2, Pt1);
  501.     PT_SUB(V2, Pt3, Pt2);
  502.     GMVecCrossProd(Plane, V1, V2);
  503.     PT_NORMALIZE(Plane);
  504.  
  505.     Plane[3] = -DOT_PROD(Plane, Pt1);
  506.  
  507. #ifdef DEBUG1
  508.     printf("CGPlaneFrom3Points exit\n");
  509. #endif /* DEBUG1 */
  510.  
  511.     return TRUE;
  512. }
  513.  
  514. /*****************************************************************************
  515. * DESCRIPTION:                                                               M
  516. *   Routine to compute the closest point on a given 3d line to a given 3d    M
  517. * point. the line is prescribed using a point on it (Pl) and vector (Vl).    M
  518. *                                                                            *
  519. * PARAMETERS:                                                                M
  520. *   Point:         To find the closest to on the line.                       M
  521. *   Pl, Vl:        Position and direction that defines the line.             M
  522. *   ClosestPoint:  Where closest point found on the line is to be saved.     M
  523. *                                                                            *
  524. * RETURN VALUE:                                                              M
  525. *   void                                                                     M
  526. *                                                                            *
  527. * KEYWORDS:                                                                  M
  528. *   CGPointFromPointLine, point line distance                                M
  529. *****************************************************************************/
  530. void CGPointFromPointLine(PointType Point,
  531.               PointType Pl,
  532.               PointType Vl,
  533.               PointType ClosestPoint)
  534. {
  535.     int i;
  536.     PointType V1, V2;
  537.     RealType CosAlfa, VecMag;
  538.  
  539. #ifdef DEBUG1
  540.     printf("CGPointFromLinePlane entered\n");
  541. #endif /* DEBUG1 */
  542.  
  543.     for (i = 0; i < 3; i++) {
  544.         V1[i] = Point[i] - Pl[i];
  545.         V2[i] = Vl[i];
  546.     }
  547.     VecMag = GMVecLength(V1);
  548.     GMVecNormalize(V1);/* Normalized vector from Point to a point on line Pl.*/
  549.     GMVecNormalize(V2);            /* Normalized line direction vector. */
  550.  
  551.     CosAlfa = GMVecDotProd(V1, V2);/* Find the angle between the two vectors.*/
  552.  
  553.     /* Find P1 - the closest point to Point on the line: */
  554.     for (i = 0; i < 3; i++)
  555.     ClosestPoint[i] = Pl[i] + V2[i] * CosAlfa * VecMag;
  556.  
  557. #ifdef DEBUG1
  558.     printf("CGPointFromLinePlane exit\n");
  559. #endif /* DEBUG1 */
  560. }
  561.  
  562. /*****************************************************************************
  563. * DESCRIPTION:                                                               M
  564. * Routine to compute the disstance between a 3d point and a 3d line.         M
  565. *   The line is prescribed using a point on it (Pl) and vector (Vl).         M
  566. *                                                                            *
  567. * PARAMETERS:                                                                M
  568. *   Point:         To find the distance to on the line.                      M
  569. *   Pl, Vl:        Position and direction that defines the line.             M
  570. *                                                                            *
  571. * RETURN VALUE:                                                              M
  572. *   RealType:      The computed distance.                                    M
  573. *                                                                            *
  574. * KEYWORDS:                                                                  M
  575. *   CGDistPointLine, point line distance                                     M
  576. *****************************************************************************/
  577. RealType CGDistPointLine(PointType Point, PointType Pl, PointType Vl)
  578. {
  579.     RealType t;
  580.     PointType Ptemp;
  581.  
  582. #ifdef DEBUG1
  583.     printf("CGDistPointLine entered\n");
  584. #endif /* DEBUG1 */
  585.  
  586.     CGPointFromPointLine(Point, Pl, Vl, Ptemp);/* Find closest point on line.*/
  587.     t = CGDistPointPoint(Point, Ptemp);
  588.  
  589. #ifdef DEBUG1
  590.     printf("CGDistPointLine exit\n");
  591. #endif /* DEBUG1 */
  592.  
  593.     return t;
  594. }
  595.  
  596. /*****************************************************************************
  597. * DESCRIPTION:                                                               M
  598. *   Routine to compute the distance between a Point and a Plane. The Plane   M
  599. * is prescribed using its four coefficients : Ax + By + Cz + D = 0 given as  M
  600. * four elements vector.                                 M
  601. *                                                                            *
  602. * PARAMETERS:                                                                M
  603. *   Point:         To find the distance to on the plane.                     M
  604. *   Plane:         To find the distance to on the point.                     M
  605. *                                                                            *
  606. * RETURN VALUE:                                                              M
  607. *   RealType:      The computed distance.                                    M
  608. *                                                                            *
  609. * KEYWORDS:                                                                  M
  610. *   CGDistPointPlane, point plane distance                                   M
  611. *****************************************************************************/
  612. RealType CGDistPointPlane(PointType Point, PlaneType Plane)
  613. {
  614.     RealType t;
  615.  
  616. #ifdef DEBUG1
  617.     printf("CGDistPointPlane entered\n");
  618. #endif /* DEBUG1 */
  619.  
  620.     t = ((Plane[0] * Point [0] +
  621.       Plane[1] * Point [1] +
  622.       Plane[2] * Point [2] +
  623.       Plane[3]) /
  624.      sqrt(SQR(Plane[0]) + SQR(Plane[1]) + SQR(Plane[2])));
  625.  
  626. #ifdef DEBUG1
  627.     printf("CGDistPointPlane exit\n");
  628. #endif /* DEBUG1 */
  629.  
  630.     return t;
  631. }
  632.  
  633. /*****************************************************************************
  634. * DESCRIPTION:                                                               M
  635. * Routine to find the intersection point of a line and a plane (if any).     M
  636. *   The Plane is prescribed using four coefficients : Ax + By + Cz + D = 0   M
  637. * given as four elements vector. The line is define via a point on it Pl and M
  638. * a direction vector Vl. Returns TRUE only if such point exists.             M
  639. *                                                                            *
  640. * PARAMETERS:                                                                M
  641. *   Pl, Vl:        Position and direction that defines the line.             M
  642. *   Plane:         To find the intersection with the line.                   M
  643. *   InterPoint:    Where the intersection occured.                           M
  644. *   t:             Parameter along the line of the intersection location     M
  645. *                  (as Pl + Vl * t).                                         M
  646. *                                                                            *
  647. * RETURN VALUE:                                                              M
  648. *   int:            TRUE, if successful.                                     M
  649. *                                                                            *
  650. * KEYWORDS:                                                                  M
  651. *   CGPointFromLinePlane, line plane intersection                            M
  652. *****************************************************************************/
  653. int CGPointFromLinePlane(PointType Pl,
  654.              PointType Vl,
  655.              PlaneType Plane,
  656.              PointType InterPoint,
  657.              RealType *t)
  658. {
  659.     int i;
  660.     RealType DotProd;
  661.  
  662. #ifdef DEBUG1
  663.     printf("CGPointFromLinePlane entered\n");
  664. #endif /* DEBUG1 */
  665.  
  666.     /* Check to see if they are vertical - no intersection at all! */
  667.     DotProd = GMVecDotProd(Vl, Plane);
  668.     if (ABS(DotProd) < IRIT_EPSILON)
  669.     return FALSE;
  670.  
  671.     /* Else find t in line such that the plane equation plane is satisfied: */
  672.     *t = (-Plane[3] - Plane[0] * Pl[0] - Plane[1] * Pl[1] - Plane[2] * Pl[2])
  673.                                 / DotProd;
  674.  
  675.     /* And so find the intersection point which is at that t: */
  676.     for (i = 0; i < 3; i++)
  677.     InterPoint[i] = Pl[i] + *t * Vl[i];
  678.  
  679. #ifdef DEBUG1
  680.     printf("CGPointFromLinePlane exit\n");
  681. #endif /* DEBUG1 */
  682.  
  683.     return TRUE;
  684. }
  685.  
  686. /*****************************************************************************
  687. * DESCRIPTION:                                                               M
  688. * Routine to find the intersection point of a line and a plane (if any).     M
  689. *   The Plane is prescribed using four coefficients : Ax + By + Cz + D = 0   M
  690. * given as four elements vector. The line is define via a point on it Pl and M
  691. * a direction vector Vl. Returns TRUE only if such point exists.             M
  692. *   this routine accepts solutions only for t between zero and one.          M
  693. *                                                                            *
  694. * PARAMETERS:                                                                M
  695. *   Pl, Vl:        Position and direction that defines the line.             M
  696. *   Plane:         To find the intersection with the line.                   M
  697. *   InterPoint:    Where the intersection occured.                           M
  698. *   t:             Parameter along the line of the intersection location     M
  699. *                  (as Pl + Vl * t).                                         M
  700. *                                                                            *
  701. * RETURN VALUE:                                                              M
  702. *   int:            TRUE, if successful and t between zero and one.          M
  703. *                                                                            *
  704. * KEYWORDS:                                                                  M
  705. *   CGPointFromLinePlane01, line plane intersection                          M
  706. *****************************************************************************/
  707. int CGPointFromLinePlane01(PointType Pl,
  708.                PointType Vl,
  709.                PlaneType Plane,
  710.                PointType InterPoint,
  711.                RealType *t)
  712. {
  713.     int i;
  714.     RealType DotProd;
  715.  
  716. #ifdef DEBUG1
  717.     printf("CGPointFromLinePlane01 entered\n");
  718. #endif /* DEBUG1 */
  719.  
  720.     /* Check to see if they are vertical - no intersection at all! */
  721.     DotProd = GMVecDotProd(Vl, Plane);
  722.     if (ABS(DotProd) < IRIT_EPSILON)
  723.     return FALSE;
  724.  
  725.     /* Else find t in line such that the plane equation plane is satisfied: */
  726.     *t = (-Plane[3] - Plane[0] * Pl[0] - Plane[1] * Pl[1] - Plane[2] * Pl[2])
  727.                                 / DotProd;
  728.  
  729.     if ((*t < 0.0 && !IRIT_APX_EQ(*t, 0.0)) ||      /* Not in parameter range. */
  730.     (*t > 1.0 && !IRIT_APX_EQ(*t, 1.0)))
  731.     return FALSE;
  732.  
  733.     /* And so find the intersection point which is at that t: */
  734.     for (i = 0; i < 3; i++)
  735.     InterPoint[i] = Pl[i] + *t * Vl[i];
  736.  
  737. #ifdef DEBUG1
  738.     printf("CGPointFromLinePlane01 exit\n");
  739. #endif /* DEBUG1 */
  740.  
  741.     return TRUE;
  742. }
  743.  
  744. /*****************************************************************************
  745. * DESCRIPTION:                                                               M
  746. *   Routine to find the two points Pti on the lines (Pli, Vli) ,   i = 1, 2  M
  747. * with the minimal Euclidian distance between them. In other words, the      M
  748. * distance between Pt1 and Pt2 is defined as distance between the two lines. M
  749. *   The two points are calculated using the fact that if V = (Vl1 cross Vl2) M
  750. * then these two points are the intersection point between the following:    M
  751. * Point 1 - a plane (defined by V and line1) and the line line2.             M
  752. * Point 2 - a plane (defined by V and line2) and the line line1.             M
  753. *   This function returns TRUE iff the two lines are not parallel!           M
  754. *                                                                            M
  755. * PARAMETERS:                                                                M
  756. *   Pl1, Vl1:  Position and direction defining the first line.               M
  757. *   Pl2, Vl2:  Position and direction defining the second line.              M
  758. *   Pt1:       Point on Pt1 that is closest to line 2.                       M
  759. *   t1:        Parameter value of Pt1 as (Pl1 + Vl1 * t1).                   M
  760. *   Pt2:       Point on Pt2 that is closest to line 1.                       M
  761. *   t2:        Parameter value of Pt2 as (Pl2 + Vl2 * t2).                   M
  762. *                                                                            *
  763. * RETURN VALUE:                                                              M
  764. *   int:       TRUE, if successful.                                          M
  765. *                                                                            *
  766. * KEYWORDS:                                                                  M
  767. *   CG2PointsFromLineLine, line line distance                                M
  768. *****************************************************************************/
  769. int CG2PointsFromLineLine(PointType Pl1,
  770.               PointType Vl1,
  771.               PointType Pl2,
  772.               PointType Vl2,
  773.               PointType Pt1,
  774.               RealType *t1,
  775.               PointType Pt2,
  776.               RealType *t2)
  777. {
  778.     int i;
  779.     PointType Vtemp;
  780.     PlaneType Plane1, Plane2;
  781.  
  782. #ifdef DEBUG1
  783.     printf("CG2PointsFromLineLine entered\n");
  784. #endif /* DEBUG1 */
  785.  
  786.     GMVecCrossProd(Vtemp, Vl1, Vl2);   /* Check to see if they are parallel. */
  787.     if (GMVecLength(Vtemp) < IRIT_EPSILON) {
  788.     for (i = 0; i < 3; i++)
  789.         Pt1[i] = Pl1[i];             /* Pick point on line1 and find */
  790.     CGPointFromPointLine(Pl1, Pl2, Vl2, Pt2); /* closest point on line2. */
  791.         return FALSE;
  792.     }
  793.  
  794.     /* Define the two planes: 1) Vl1, Pl1, Vtemp and 2) Vl2, Pl2, Vtemp         */
  795.     /* Note this sets the first 3 elements A, B, C out of the 4...         */
  796.     GMVecCrossProd(Plane1, Vl1, Vtemp);            /* Find the A, B, C coef.'s. */
  797.     GMVecNormalize(Plane1);
  798.     GMVecCrossProd(Plane2, Vl2, Vtemp);            /* Find the A, B, C coef.'s. */
  799.     GMVecNormalize(Plane2);
  800.  
  801.     /* and now use a point on the plane to find the 4th coef. D: */
  802.     Plane1[3] = (-GMVecDotProd(Plane1, Pl1)); /* VecDotProd uses only first  */
  803.     Plane2[3] = (-GMVecDotProd(Plane2, Pl2)); /* three elements in vec.      */
  804.  
  805.     /* Thats it! now we should solve for the intersection point between a    */
  806.     /* line and a plane but we already familiar with this problem...         */
  807.     i = CGPointFromLinePlane(Pl1, Vl1, Plane2, Pt1, t1) &&
  808.     CGPointFromLinePlane(Pl2, Vl2, Plane1, Pt2, t2);
  809.  
  810. #ifdef DEBUG1
  811.     printf("CG2PointsFromLineLine exit\n");
  812. #endif /* DEBUG1 */
  813.  
  814.     return i;
  815. }
  816.  
  817. /*****************************************************************************
  818. * DESCRIPTION:                                                               M
  819. *   Routine to find the distance between two lines (Pli, Vli) ,  i = 1, 2.   M
  820. *                                                                            *
  821. * PARAMETERS:                                                                M
  822. *   Pl1, Vl1:  Position and direction defining the first line.               M
  823. *   Pl2, Vl2:  Position and direction defining the second line.              M
  824. *                                                                            *
  825. * RETURN VALUE:                                                              M
  826. *   RealType:   Distance between the two lines.                              M
  827. *                                                                            *
  828. * KEYWORDS:                                                                  M
  829. *   CGDistLineLine, line line distance                                       M
  830. *****************************************************************************/
  831. RealType CGDistLineLine(PointType Pl1,
  832.             PointType Vl1,
  833.             PointType Pl2,
  834.             PointType Vl2)
  835. {
  836.     RealType t1, t2;
  837.     PointType Ptemp1, Ptemp2;
  838.  
  839. #ifdef DEBUG1
  840.     printf("CGDistLineLine entered\n");
  841. #endif /* DEBUG1 */
  842.  
  843.     CG2PointsFromLineLine(Pl1, Vl1, Pl2, Vl2, Ptemp1, &t1, Ptemp2, &t2);
  844.     t1 = CGDistPointPoint(Ptemp1, Ptemp2);
  845.  
  846. #ifdef DEBUG1
  847.     printf("CGDistLineLine exit\n");
  848. #endif /* DEBUG1 */
  849.  
  850.     return t1;
  851. }
  852.  
  853. /*****************************************************************************
  854. * DESCRIPTION:                                                               M
  855. * Routine that implements "Jordan Theorem":                                  M
  856. *   Fire a ray from a given point and find the number of intersections of a  M
  857. * ray with the polygon, excluding the given point Pt (start of ray) itself,  M
  858. * if on polygon boundary. The ray is fired in +X (Axes == 0) or +Y if        M
  859. * (Axes == 1).                                                               M
  860. *   Only the X/Y coordinates of the polygon are taken into account, i.e. the M
  861. * orthogonal projection of the polygon on an X/Y parallel plane (equal to    M
  862. * polygon itself if on X/Y parallel plane...).                     M
  863. *   Note that if the point is on polygon boundary, the ray should not be in  M
  864. * its edge direction.                                 M
  865. *                                         M
  866. * Algorithm:                                     M
  867. * 1. 1.1. Set NumOfIntersection = 0;                         M
  868. *    1.2. Find vertex V not on Ray level and set AlgState to its level       M
  869. *         (below or above the ray level). If none goto 3;             M
  870. *    1.3. Mark VStart = V;                             M
  871. * 2. Do                                         M
  872. *    2.1. While State(V) == AlgState do                         M
  873. *        2.1.1. V = V -> Pnext;                         M
  874. *        2.1.2. If V == VStart goto 3;                     M
  875. *    2.2. IntersectionMinX = INFINITY;                         M
  876. *    2.3. While State(V) == ON_RAY do                         M
  877. *        2.3.1. IntersectionMin = MIN(IntersectionMin, V -> Coord[Axes]); M
  878. *        2.3.2. V = V -> Pnext;                         M
  879. *    2.4. If State(V) != AlgState do                         M
  880. *        2.4.1. Find the intersection point between polygon edge         M
  881. *           Vlast, V and the Ray and update IntersectionMin if         M
  882. *           lower than it.                         M
  883. *        2.4.2. If IntersectionMin is greater than Pt[Axes] increase         M
  884. *           the NumOfIntersection counter by 1.                   M
  885. *    2.5. AlgState = State(V);                             M
  886. *    2.6. goto 2.2.                                 M
  887. * 3. Return NumOfIntersection;                             M
  888. *                                                                            *
  889. * PARAMETERS:                                                                M
  890. *   Pl:        To compute "Jordan Theorem" for the given ray.                M
  891. *   PtRay:     Origin of ray.                                                M
  892. *   RayAxes:   Direction of ray. 0 for X, 1 for Y, etc.                      M
  893. *                                                                            *
  894. * RETURN VALUE:                                                              M
  895. *   int:       Number of intersections of ray with the polygon.              M
  896. *                                                                            *
  897. * KEYWORDS:                                                                  M
  898. *   CGPolygonRayInter, ray polygon intersection, Jordan theorem              M
  899. *****************************************************************************/
  900. int CGPolygonRayInter(IPPolygonStruct *Pl, PointType PtRay, int RayAxes)
  901. {
  902.     int NewState, AlgState, RayOtherAxes,
  903.     NumOfInter = 0,
  904.     Quit = FALSE;
  905.     RealType InterMin, Inter, t;
  906.     IPVertexStruct *V, *VStart,
  907.     *VLast = NULL;
  908.  
  909.     RayOtherAxes = (RayAxes == 1 ? 0 : 1);     /* Other dir: X -> Y, Y -> X. */
  910.  
  911.     /* Stage 1 - find a vertex below the ray level: */
  912.     V = VStart = Pl -> PVertex;
  913.     do {
  914.     if ((AlgState = CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes))
  915.                             != ON_RAY)
  916.         break;
  917.     V = V -> Pnext;
  918.     }
  919.     while (V != VStart && V != NULL);
  920.     if (AlgState == ON_RAY)
  921.     return 0;
  922.     VStart = V; /* Vertex Below Ray level */
  923.  
  924.     /* Stage 2 - scan the vertices and count number of intersections. */
  925.     while (!Quit) {
  926.     /* Stage 2.1. : */
  927.     while (CGPointRayRelation(V -> Coord, PtRay,
  928.                   RayOtherAxes) == AlgState) {
  929.         VLast = V;
  930.         V = V -> Pnext;
  931.         if (V == VStart) {
  932.         Quit = TRUE;
  933.         break;
  934.         }
  935.     }
  936.     InterMin = INFINITY;
  937.  
  938.     /* Stage 2.2. : */
  939.     while (CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes) == ON_RAY) {
  940.         InterMin = MIN(InterMin, V -> Coord[RayAxes]);
  941.         VLast = V;
  942.         V = V -> Pnext;
  943.         if (V == VStart)
  944.         Quit = TRUE;
  945.     }
  946.  
  947.     /* Stage 2.3. : */
  948.     if ((NewState = CGPointRayRelation(V -> Coord, PtRay, RayOtherAxes))
  949.                                 != AlgState) {
  950.         /* Stage 2.3.1 Intersection with ray is in middle of edge: */
  951.         t = (PtRay[RayOtherAxes] - V -> Coord[RayOtherAxes]) /
  952.         (VLast -> Coord[RayOtherAxes] - V -> Coord[RayOtherAxes]);
  953.         Inter = VLast -> Coord[RayAxes] * t +
  954.             V -> Coord[RayAxes] * (1.0 - t);
  955.         InterMin = MIN(InterMin, Inter);
  956.  
  957.         /* Stage 2.3.2. comp. with ray base and inc. # of inter if above.*/
  958.         if (InterMin > PtRay[RayAxes] &&
  959.         !IRIT_APX_EQ(InterMin, PtRay[RayAxes]))
  960.         NumOfInter++;
  961.     }
  962.  
  963.     AlgState = NewState;
  964.     }
  965.  
  966.     return NumOfInter;
  967. }
  968.  
  969. /*****************************************************************************
  970. * DESCRIPTION:                                                               *
  971. *   Routine to returns the relation between the ray level and a given point, *
  972. * to be used in the CGPolygonRayInter routine above.                 *
  973. *                                                                            *
  974. * PARAMETERS:                                                                *
  975. *   Pt:       Given point.                                                   *
  976. *   PtRay:    Given ray.                                                     *
  977. *   Axes:     Given axes.                                                    *
  978. *                                                                            *
  979. * RETURN VALUE:                                                              *
  980. *   int:      Pt is either above below or on the ray.                        *
  981. *****************************************************************************/
  982. static int CGPointRayRelation(PointType Pt, PointType PtRay, int Axes)
  983. {
  984.     if (IRIT_APX_EQ(PtRay[Axes], Pt[Axes]))
  985.         return ON_RAY;
  986.     else if (PtRay[Axes] < Pt[Axes])
  987.         return ABOVE_RAY;
  988.     else
  989.     return BELOW_RAY;
  990. }
  991.  
  992. /*****************************************************************************
  993. * DESCRIPTION:                                                               M
  994. * Same as CGPolygonRayInter but for arbitrary oriented polygon.             M
  995. *   The polygon is transformed into the XY plane and then CGPolygonRayInter  M
  996. * is invoked on it.                                                          M
  997. *                                                                            *
  998. * PARAMETERS:                                                                M
  999. *   Pl:        To compute "Jordan Theorem" for the given ray.                M
  1000. *   PtRay:     Origin of ray.                                                M
  1001. *   RayAxes:   Direction of ray. 0 for X, 1 for Y, etc.                      M
  1002. *                                                                            *
  1003. * RETURN VALUE:                                                              M
  1004. *   int:       Number of intersections of ray with the polygon.              M
  1005. *                                                                            *
  1006. * KEYWORDS:                                                                  M
  1007. *   CGPolygonRayInter3D, ray polygon intersection, Jordan theorem            M
  1008. *****************************************************************************/
  1009. int CGPolygonRayInter3D(IPPolygonStruct *Pl, PointType PtRay, int RayAxes)
  1010. {
  1011.     int i;
  1012.     MatrixType RotMat;
  1013.     IPVertexStruct *V, *VHead;
  1014.     IPPolygonStruct *RotPl;
  1015.     PointType RotPt;
  1016.  
  1017.     /* Make a copy of original to work on. */
  1018.     RotPl = IPAllocPolygon(1, Pl ->Tags, CopyVertexList(Pl -> PVertex), NULL);
  1019.  
  1020.     /* Make sure list is circular. */
  1021.     V = IritPrsrGetLastVrtx(RotPl -> PVertex);
  1022.     if (V -> Pnext == NULL);
  1023.     V -> Pnext = RotPl -> PVertex;
  1024.  
  1025.     /* Bring the polygon to a XY parallel plane by rotation. */
  1026.     GenRotateMatrix(RotMat, Pl -> Plane);
  1027.     V = VHead = RotPl -> PVertex;
  1028.     do {                    /* Transform the polygon itself. */
  1029.     MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);
  1030.     V = V -> Pnext;
  1031.     }
  1032.     while (V != VHead);
  1033.     MatMultVecby4by4(RotPt, PtRay, RotMat);
  1034.  
  1035.     i = CGPolygonRayInter(RotPl, RotPt, RayAxes);
  1036.  
  1037.     IPFreePolygonList(RotPl);
  1038.  
  1039.     return i;
  1040. }
  1041.  
  1042. /*****************************************************************************
  1043. * DESCRIPTION:                                                               M
  1044. *   Routine to prepar a transformation martix to do the following (in this   M
  1045. * order): scale by Scale, rotate such that the Z axis is in direction Dir    M
  1046. * and then translate by Trans.                             M
  1047. *    Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is M
  1048. * used to form the second line (the first 3 lines set the rotation), and     M
  1049. * finally Scale is used to scale first 3 lines/columns to the needed scale:  M
  1050. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord V
  1051. *                |  Bx  By  Bz  0 |  system into T, N & B as required and    V
  1052. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |  then translate it to C. T, N, B are     V
  1053. *                |  Cx  Cy  Cz  1 |  scaled by Scale.                 V
  1054. * N is exactly Dir (unit vec) but we got freedom on T & B which must be on   M
  1055. * a plane perpendicular to N and perpendicular between them but thats all!   M
  1056. * T is therefore selected using this (heuristic ?) algorithm:             M
  1057. * Let P be the axis of which the absolute N coefficient is the smallest.     M
  1058. * Let B be (N cross P) and T be (B cross N).                     M
  1059. *                                                                            *
  1060. * PARAMETERS:                                                                M
  1061. *   Mat:       To place the computed transformation.                         M
  1062. *   Trans:     Translation factor.                                           M
  1063. *   Dir:       Direction to take Z axis to.                                  M
  1064. *   Scale:     Scaling factor.                                               M
  1065. *                                                                            *
  1066. * RETURN VALUE:                                                              M
  1067. *   void                                                                     M
  1068. *                                                                            *
  1069. * KEYWORDS:                                                                  M
  1070. *    CGGenTransMatrixZ2Dir, transformations                                  M
  1071. *****************************************************************************/
  1072. void CGGenTransMatrixZ2Dir(MatrixType Mat,
  1073.                VectorType Trans,
  1074.                VectorType Dir,
  1075.                RealType Scale)
  1076. {
  1077.     int i, j;
  1078.     RealType R;
  1079.     VectorType DirN, T, B, P;
  1080.     MatrixType TempMat;
  1081.  
  1082.     PT_COPY(DirN, Dir);
  1083.     PT_NORMALIZE(DirN);
  1084.     PT_CLEAR(P);
  1085.     for (i = 1, j = 0, R = ABS(DirN[0]); i < 3; i++)
  1086.     if (R > ABS(DirN[i])) {
  1087.         R = DirN[i];
  1088.         j = i;
  1089.     }
  1090.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  1091.  
  1092.     GMVecCrossProd(B, DirN, P);                  /* Calc the bi-normal. */
  1093.     GMVecCrossProd(T, B, DirN);                /* Calc the tangent. */
  1094.  
  1095.     MatGenUnitMat(Mat);
  1096.     for (i = 0; i < 3; i++) {
  1097.     Mat[0][i] = T[i];
  1098.     Mat[1][i] = B[i];
  1099.     Mat[2][i] = DirN[i];
  1100.     }
  1101.     MatGenMatUnifScale(Scale, TempMat);
  1102.     MatMultTwo4by4(Mat, TempMat, Mat);
  1103.  
  1104.     MatGenMatTrans(Trans[0], Trans[1], Trans[2], TempMat);
  1105.     MatMultTwo4by4(Mat, Mat, TempMat);
  1106. }
  1107. /*****************************************************************************
  1108. * DESCRIPTION:                                                               M
  1109. *   Routine to prepar a transformation martix to do the following (in this   M
  1110. * order): scale by Scale, rotate such that the Z axis is in direction Dir    M
  1111. * and X axis is direction T and then translate by Trans.             M
  1112. *    Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is M
  1113. * used to form the second line (the first 3 lines set the rotation), and     M
  1114. * finally Scale is used to scale first 3 lines/columns to the needed scale:  M
  1115. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord V
  1116. *                |  Bx  By  Bz  0 |  system into T, N & B as required and    V
  1117. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |  then translate it to C. T, N, B are     V
  1118. *                |  Cx  Cy  Cz  1 |  scaled by Scale.                 V
  1119. * N is exactly Dir (unit vec) and T is exactly Dir2.                 M
  1120. *                                                                            *
  1121. * PARAMETERS:                                                                M
  1122. *   Mat:       To place the computed transformation.                         M
  1123. *   Trans:     Translation factor.                                           M
  1124. *   Dir:       Direction to take Z axis to.                                  M
  1125. *   Dir2:      Direction to take X axis to.                                  M
  1126. *   Scale:     Scaling factor.                                               M
  1127. *                                                                            *
  1128. * RETURN VALUE:                                                              M
  1129. *   void                                                                     M
  1130. *                                                                            *
  1131. * KEYWORDS:                                                                  M
  1132. *    CGGenTransMatrixZ2Dir2, transformations                                 M
  1133. *****************************************************************************/
  1134. void CGGenTransMatrixZ2Dir2(MatrixType Mat,
  1135.                 VectorType Trans,
  1136.                 VectorType Dir,
  1137.                 VectorType Dir2,
  1138.                 RealType Scale)
  1139. {
  1140.     int i;
  1141.     VectorType DirN, Dir2N, B;
  1142.     MatrixType TempMat;
  1143.  
  1144.     PT_COPY(DirN, Dir);
  1145.     PT_NORMALIZE(DirN);
  1146.     PT_COPY(Dir2N, Dir2);
  1147.     PT_NORMALIZE(Dir2N);
  1148.  
  1149.     GMVecCrossProd(B, DirN, Dir2N);              /* Calc the bi-normal. */
  1150.  
  1151.     MatGenUnitMat(Mat);
  1152.     for (i = 0; i < 3; i++) {
  1153.     Mat[0][i] = Dir2N[i];
  1154.     Mat[1][i] = B[i];
  1155.     Mat[2][i] = DirN[i];
  1156.     }
  1157.     MatGenMatUnifScale(Scale, TempMat);
  1158.     MatMultTwo4by4(Mat, TempMat, Mat);
  1159.  
  1160.     MatGenMatTrans(Trans[0], Trans[1], Trans[2], TempMat);
  1161.     MatMultTwo4by4(Mat, Mat, TempMat);
  1162. }
  1163.