home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 3 / goldfish_volume_3.bin / files / gfx / 3d / irit / geom_lib / geomvals.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-03  |  13.0 KB  |  319 lines
  1. /******************************************************************************
  2. * GeomVals.c - Area, Volume, and counts on polygonal objects.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. #include <stdio.h>
  8. #include <ctype.h>
  9. #include <math.h>
  10. #include <string.h>
  11. #include "allocate.h"
  12. #include "convex.h"
  13. #include "geomat3d.h"
  14. #include "geomvals.h"
  15.  
  16. static RealType PolygonXYArea(IPVertexStruct *VHead);
  17. static RealType Polygon3VrtxXYArea(PointType Pt1,
  18.                    PointType Pt2,
  19.                    PointType Pt3);
  20.  
  21. /*****************************************************************************
  22. * DESCRIPTION:                                                               M
  23. *  Routine to evaluate the Area of the given geom. object, in object unit.   M
  24. * Algorithm (for each polygon):                    V3         V
  25. * 1. Set Polygon Area to be zero.                   /\         V
  26. *    Make a copy of the original polygon             /      \          V
  27. *    and transform it to a XY parallel plane.           /        \V2         V
  28. *    Find the minimum Y value of the polygon           V4/         |         V
  29. *    in the XY plane.                     \         |         V
  30. * 2. Let V(0) be the first vertex, V(n) the last one.      \         |         V
  31. *    For i goes from 0 to n-1 add to Area the area            \_______|       V
  32. *    below edge V(i), V(i+1):                     V0      V1      V
  33. *    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2             V
  34. *    where V(i)y' is V(i)y - MinY, where MinY is polygon minimum Y value.    V
  35. * 3. The result of step 2 is the area of the polygon itself.              V
  36. *    However, it might be negative, so take the absolute result of step 2    V
  37. *    and add it to the global ObjectArea.                     V
  38. * Note step 2 is performed by another auxiliary routine: PolygonXYArea.      M
  39. *                                                                            *
  40. * PARAMETERS:                                                                M
  41. *   PObj:     A polyhedra object to compute its surface area.                M
  42. *                                                                            *
  43. * RETURN VALUE:                                                              M
  44. *   double:   The area of object PObj.                                       M
  45. *                                                                            *
  46. * KEYWORDS:                                                                  M
  47. *   PolyObjectArea, area                                                     M
  48. *****************************************************************************/
  49. double PolyObjectArea(IPObjectStruct *PObj)
  50. {
  51.     RealType
  52.     ObjectArea = 0.0;
  53.     MatrixType RotMat;
  54.     IPPolygonStruct *Pl;
  55.     IPVertexStruct *V, *VHead;
  56.  
  57.     if (!IP_IS_POLY_OBJ(PObj))
  58.     IritFatalError("Geometric property requested on non polygonal object");
  59.  
  60.     if (IP_IS_POLYLINE_OBJ(PObj)) {
  61.     return 0.0;
  62.     }
  63.  
  64.     Pl = PObj -> U.Pl;
  65.     while (Pl != NULL) {
  66.     /* Dont trans. original object. */
  67.     V = VHead = CopyVertexList(Pl -> PVertex);
  68.     /* Create the trans. matrix to transform the polygon to XY parallel  */
  69.     GenRotateMatrix(RotMat, Pl -> Plane);               /* plane. */
  70.     do {
  71.         MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);
  72.  
  73.         V = V -> Pnext;
  74.     }
  75.     while (V != NULL && V != VHead);
  76.  
  77.     ObjectArea += PolygonXYArea(VHead);
  78.  
  79.     IPFreeVertexList(VHead);          /* Free the vertices list. */
  80.  
  81.     Pl = Pl -> Pnext;
  82.     }
  83.  
  84.     return ObjectArea;
  85. }
  86.  
  87. /*****************************************************************************
  88. * DESCRIPTION:                                                               *
  89. *   Routine to evaluate the area of the given polygon projected on the XY    *
  90. * plane.                                                                 *
  91. *   The polygon does not have to be on a XY parallel plane, as only its XY   *
  92. * projection is considered (Z is ignored).                     *
  93. *   Returned the area of its XY parallel projection.                 *
  94. *   See GeomObjectArea above for algorithm.                     *
  95. *                                                                            *
  96. * PARAMETERS:                                                                *
  97. *   VHead:     Of vertex list to compute area for.                           *
  98. *                                                                            *
  99. * RETURN VALUE:                                                              *
  100. *   RealType:  Computed area.                                                *
  101. *****************************************************************************/
  102. static RealType PolygonXYArea(IPVertexStruct *VHead)
  103. {
  104.     RealType PolygonArea = 0.0, MinY;
  105.     IPVertexStruct *V = VHead, *Vnext;
  106.  
  107.     MinY = V -> Coord[1];
  108.     V = V -> Pnext;
  109.     while (V != VHead && V != NULL) {
  110.     if (MinY > V -> Coord[1])
  111.         MinY = V -> Coord[1];
  112.  
  113.     V = V -> Pnext;
  114.     }
  115.  
  116.     if (V == NULL)
  117.     V = VHead;
  118.     Vnext = V -> Pnext;
  119.     MinY *= 2.0;          /* Instead of subtracting twice each time. */
  120.     do {
  121.     /* Evaluate area below edge V->Vnext relative to Y level MinY. Note  */
  122.     /* it can come out negative, but thats o.k. as the sum of all these  */
  123.     /* quadraliterals should be exactly the area (up to correct sign).   */
  124.     PolygonArea += (Vnext -> Coord[0] - V -> Coord[0]) *
  125.             (Vnext -> Coord[1] + V -> Coord[1] - MinY) / 2.0;
  126.     V = Vnext;
  127.     Vnext = V -> Pnext == NULL ? VHead : V -> Pnext;
  128.     }
  129.     while (V != VHead && V != NULL);
  130.  
  131.     return ABS(PolygonArea);
  132. }
  133.  
  134. /*****************************************************************************
  135. * DESCRIPTION:                                                               M
  136. *   Routine to evaluate the Volume of the given geom object, in object unit. M
  137. *   This routine has a side effect that all non-convex polygons will be      M
  138. * splitted to convex ones.                             M
  139. * Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     M
  140. *                                V3         V
  141. * 1. Set Polygon Area to be zero.                   /\         V
  142. *    Let V(0) be the first vertex, V(n) the last.         /      \          V
  143. *    For i goes from 1 to n-1 form triangles           /        \V2         V
  144. *                   by V(0), V(i), V(i+1).             V4/         |         V
  145. *    For each such triangle di:                 \         |         V
  146. *    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |         V
  147. *         with the minimum Z - TriMinY.                 \_______|       V
  148. *    1.2. The volume below V(0), V(i), V(i+1) triangle,         V0      V1      V
  149. *      relative to ObjMinZ level, is the sum of:                 V
  150. *      1.2.1. volume of V'(0), V'(i), V'(i+1) - the                 V
  151. *         area of projection of V(0), V(i), V(i+1) on XY parallel     V
  152. *         plane, times (TriMinZ - ObjMinZ).                 V
  153. *      1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     V
  154. *         V"(i+1) be the projections of V(i) and V(i+1) on the plane  V
  155. *         Z = PolyZMin. The volume above 1.2.1. and below the polygon V
  156. *         (triangle!) will be: the area of quadraliteral V(i), V(i+1),V
  157. *         V"(i+1), V"(i), times distance of V(0) for quadraliteral    V
  158. *         plane divided by 3.                         V
  159. *    1.3. If Z component of polygon normal is negative add 1.2. result to    V
  160. *      ObjectVolume, else subtract it.                     V
  161. *                                                                            *
  162. * PARAMETERS:                                                                M
  163. *   PObj:      To compute volume for.                                        M
  164. *                                                                            *
  165. * RETURN VALUE:                                                              M
  166. *   double:    Computed volume.                                              M
  167. *                                                                            *
  168. * KEYWORDS:                                                                  M
  169. *   PolyObjectVolume, volume                                                 M
  170. *****************************************************************************/
  171. double PolyObjectVolume(IPObjectStruct *PObj)
  172. {
  173.     int PlaneExists;
  174.     RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist;
  175.     PointType Pt1;
  176.     PlaneType Plane;
  177.     IPPolygonStruct *Pl;
  178.     IPVertexStruct *V, *VHead, *Vnext, *Vtemp;
  179.  
  180.     if (!IP_IS_POLY_OBJ(PObj))
  181.     IritFatalError("Geometric property requested on non polygonal object");
  182.  
  183.     if (IP_IS_POLYLINE_OBJ(PObj)) {
  184.     return 0.0;
  185.     }
  186.  
  187.     ObjMinZ = INFINITY;     /* Find Object minimum Z value (used as min level). */
  188.     Pl = PObj -> U.Pl;
  189.     while (Pl != NULL) {
  190.     V = VHead = Pl -> PVertex;
  191.     do {
  192.         if (V -> Coord[2] < ObjMinZ)
  193.         ObjMinZ = V -> Coord[2];
  194.         V = V -> Pnext;
  195.     }
  196.     while (V != VHead && V != NULL);
  197.     if (V == NULL)
  198.         IritFatalError("For VOLUME polygons must be closed loops. Open loop detected.\n");
  199.  
  200.     Pl = Pl -> Pnext;
  201.     }
  202.  
  203.     ConvexPolyObject(PObj);           /* Make sure all polygons are convex. */
  204.     Pl = PObj -> U.Pl;
  205.     while (Pl != NULL) {
  206.     PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */
  207.     /* We set VHead to be vertex with min Z: */
  208.     V = Vtemp = VHead = Pl -> PVertex;
  209.     do {
  210.         if (V -> Coord[2] < Vtemp -> Coord[2])
  211.         Vtemp = V;
  212.         V = V -> Pnext;
  213.     }
  214.     while (V != VHead && V != NULL);
  215.     VHead = Vtemp;       /* Now VHead is the one with lowest Z in polygon! */
  216.     TriMinZ = VHead -> Coord[2];         /* Save this Z for fast access. */
  217.  
  218.     V = VHead -> Pnext;
  219.     Vnext = V -> Pnext;
  220.     do {
  221.         /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */
  222.         Area = Polygon3VrtxXYArea(VHead -> Coord, V -> Coord,
  223.                                    Vnext -> Coord);
  224.         PolygonVolume += Area * (TriMinZ - ObjMinZ);
  225.  
  226.         /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */
  227.         Area = sqrt(SQR(V -> Coord[0] - Vnext -> Coord[0]) + /* XY dist. */
  228.             SQR(V -> Coord[1] - Vnext -> Coord[1])) *
  229.            ((V -> Coord[2] + Vnext -> Coord[2]) / 2.0 - TriMinZ);
  230.         PT_COPY(Pt1, V -> Coord);
  231.         Pt1[2] = TriMinZ;
  232.         if ((PlaneExists =
  233.          CGPlaneFrom3Points(Plane, V -> Coord,
  234.                     Vnext -> Coord, Pt1)) == 0) {
  235.         /* Try second pt projected to Z = TriMinZ plane as third pt. */
  236.         PT_COPY(Pt1, Vnext -> Coord);
  237.         Pt1[2] = TriMinZ;
  238.         PlaneExists =
  239.             CGPlaneFrom3Points(Plane, V -> Coord, Vnext -> Coord,
  240.                                     Pt1);
  241.         }
  242.         if (PlaneExists) {
  243.         Dist = CGDistPointPlane(VHead -> Coord, Plane);
  244.         PolygonVolume += Area * ABS(Dist) / 3.0;
  245.         }
  246.  
  247.         V = Vnext;
  248.         Vnext = V -> Pnext;
  249.     }
  250.     while (Vnext != VHead);
  251.  
  252.     if (Pl -> Plane[2] < 0.0)
  253.         ObjVolume += PolygonVolume;
  254.     else
  255.         ObjVolume -= PolygonVolume;
  256.  
  257.     Pl = Pl -> Pnext;
  258.     }
  259.  
  260.     return ObjVolume;
  261. }
  262.  
  263. /*****************************************************************************
  264. * DESCRIPTION:                                                               *
  265. *  Routine to evaluate the area of the given triangle projected to the XY    *
  266. * plane, given as 3 Points. Only the X & Y components are considered.         *
  267. *  See PolyObjectArea above for algorithm.                     *
  268. *                                                                            *
  269. * PARAMETERS:                                                                *
  270. *   Pt1, Pt2, Pt3: Vertices of triangle to compute area for.                 *
  271. *                                                                            *
  272. * RETURN VALUE:                                                              *
  273. *   RealType:      Area of triangle in the XY plane.                         *
  274. *****************************************************************************/
  275. static RealType Polygon3VrtxXYArea(PointType Pt1,
  276.                    PointType Pt2,
  277.                    PointType Pt3)
  278. {
  279.     RealType PolygonArea = 0.0, MinY;
  280.  
  281.     MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0;
  282.  
  283.     PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0;
  284.     PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0;
  285.     PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0;
  286.  
  287.     return ABS(PolygonArea);
  288. }
  289.  
  290. /*****************************************************************************
  291. * DESCRIPTION:                                                               M
  292. *  Routine to count the number of polygons in the given geometric object.    M
  293. *                                                                            *
  294. * PARAMETERS:                                                                M
  295. *   PObj:       To count number of polygons in.                              M
  296. *                                                                            *
  297. * RETURN VALUE:                                                              M
  298. *   double:     Number of polygons, an integer value returned as double.     M
  299. *                                                                            *
  300. * KEYWORDS:                                                                  M
  301. *   PolyCountPolys                                                           M
  302. *****************************************************************************/
  303. double PolyCountPolys(IPObjectStruct *PObj)
  304. {
  305.     int Count = 0;
  306.     IPPolygonStruct *Pl;
  307.  
  308.     if (!IP_IS_POLY_OBJ(PObj))
  309.     IritFatalError("Geometric property requested on non polygonal object");
  310.  
  311.     Pl = PObj -> U.Pl;
  312.     while (Pl != NULL) {
  313.     Count++;
  314.     Pl = Pl -> Pnext;
  315.     }
  316.  
  317.     return ((double) Count);
  318. }
  319.