home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / IRIT / IRITS.ZIP / GEOMVALS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-05  |  9.9 KB  |  271 lines
  1. /*****************************************************************************
  2. *   "Irit" - the 3d polygonal solid modeller.                     *
  3. *                                         *
  4. * Written by:  Gershon Elber                Ver 0.2, Mar. 1990   *
  5. ******************************************************************************
  6. *  Module to evaluate geometric properties of the objects such as area,      *
  7. * volume, number of polygons etc.                         *
  8. *****************************************************************************/
  9.  
  10. #ifdef __MSDOS__
  11. #include <graphics.h>
  12. #endif /* __MSDOS__ */
  13.  
  14. #include <stdio.h>
  15. #include <ctype.h>
  16. #include <math.h>
  17. #include <string.h>
  18. #include "program.h"
  19. #include "objects.h"
  20. #include "windowsg.h"
  21. #include "geomat3d.h"
  22. #include "geomvall.h"
  23. #include "geomvalg.h"
  24. #include "allocatg.h"
  25. #include "convexg.h"
  26.  
  27. #ifndef __MSDOS__
  28. #include "xgraphic.h"
  29. #endif /* __MSDOS__ */
  30.  
  31. /*****************************************************************************
  32. *  Routine to evaluate the Area of the given geom. object, in object unit.   *
  33. * Algorithm (for each polygon):                    V3         *
  34. * 1. Set Polygon Area to be zero.                   /\         *
  35. *    Make a copy of the original polygon             /      \          *
  36. *    and transform it to a XY parallel plane.           /        \V2         *
  37. *    Find the minimum Y value of the polygon           V4/         |         *
  38. *    in the XY plane                     \         |         *
  39. * 2. Let V(0) be the first vertex, V(n) the last one       \         |         *
  40. *    for i goes from 0 to n-1 add to Area the area of         \_______|       *
  41. *    below edge V(i), V(i+1):                     V0      V1      *
  42. *    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2             *
  43. *    where V(i)y' is V(i)y - MinY, where MinY is the polygon minimum Y value *
  44. * 3. Note that the result of step 2 is the area of the polygon itself.       *
  45. *    However it might be negative, so take the absolute result of step 2 and *
  46. *    add it to the global ObjectArea.                         *
  47. * Note step 2 is performed by another aux. routine below: PolygonXYArea.     *
  48. *****************************************************************************/
  49. double GeomObjectArea(ObjectStruct *PObj)
  50. {
  51.     RealType ObjectArea = 0.0;
  52.     MatrixType RotMat;
  53.     PolygonStruct *Pl;
  54.     VertexStruct *V, *VHead;
  55.  
  56.     if (!IS_GEOM_OBJ(PObj))
  57.     FatalError("Geometric property requested on non geometric object");
  58.  
  59.     if (IS_POLYLINE_GEOM_OBJ(PObj)) {
  60.     WndwInputWindowPutStr("Warning: geometric object is polyline - zero area",
  61.                                     RED);
  62.     return 0.0;
  63.     }
  64.  
  65.     Pl = PObj -> U.Pl;
  66.     while (Pl != NULL) {
  67.     V = VHead = CopyVList(Pl -> V);      /* Dont transform original object. */
  68.     /* Create the trans. matrix to transform the polygon to XY parallel  */
  69.     GenRotateMatrix(RotMat, Pl -> Plane);               /* plane. */
  70.     do {
  71.         MatMultVecby4by4(V -> Pt, V -> Pt, RotMat);
  72.  
  73.         V = V -> Pnext;
  74.     }
  75.     while (V != NULL && V != VHead);
  76.  
  77.     ObjectArea += PolygonXYArea(VHead);
  78.  
  79.     MyFree((char *) VHead, VERTEX_TYPE);      /* Free the vertices list. */
  80.  
  81.     Pl = Pl -> Pnext;
  82.     }
  83.  
  84.     return ObjectArea;
  85. }
  86.  
  87. /*****************************************************************************
  88. *  Routine to evaluate the area of the given polygon projection on the XY    *
  89. * plane. Note the polygon does not have to be on a XY parallel plane, as     *
  90. * only its XY projection is considered (Z is ignored). Returns the area of   *
  91. * its XY parallel projection.                             *
  92. *  See GeomObjectArea above for algorithm:                     *
  93. *****************************************************************************/
  94. static RealType PolygonXYArea(VertexStruct *VHead)
  95. {
  96.     RealType PolygonArea = 0.0, MinY;
  97.     VertexStruct *V = VHead, *Vnext;
  98.  
  99.     MinY = V -> Pt[1];
  100.     V = V -> Pnext;
  101.     while (V != VHead && V != NULL /* Should not happen! */) {
  102.     if (MinY > V -> Pt[1]) MinY = V -> Pt[1];
  103.  
  104.     V = V -> Pnext;
  105.     }
  106.  
  107.     Vnext = V -> Pnext;
  108.     MinY *= 2.0;          /* Instead of subtracting twice each time. */
  109.     do {
  110.     /* Evaluate area below edge V-Vnext relative to Y level MinY. Note   */
  111.     /* it can come out negative, but thats o.k. as the sum of all these  */
  112.     /* quadraliterals should be exactly the area (up to correct sign).   */
  113.     PolygonArea += (Vnext -> Pt[0] - V -> Pt[0]) *
  114.                 (Vnext -> Pt[1] + V -> Pt[1] - MinY) / 2.0;
  115.     V = Vnext;
  116.     Vnext = V -> Pnext;
  117.     }
  118.     while (V != VHead && V != NULL /* Should not happen! */);
  119.  
  120.     return ABS(PolygonArea);
  121. }
  122.  
  123. /*****************************************************************************
  124. *   Routine to evaluate the Volume of the given geom object, in object unit. *
  125. *   This routine has a side effect that all non-convex polygons will be      *
  126. * splitted to convex ones.                             *
  127. * Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     *
  128. *                                V3         *
  129. * 1. Set Polygon Area to be zero.                   /\         *
  130. *    Let V(0) be the first vertex, V(n) the last.         /      \          *
  131. *    For i goes from 1 to n-1 form triangles           /        \V2         *
  132. *    by V(0), V(i), V(i+1). For each such           V4/         |         *
  133. *    triangle do:                     \         |         *
  134. *    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |         *
  135. *         with the minimum Z - TriMinY.                 \_______|       *
  136. *    1.2. The volume below V(0), V(i), V(i+1) triangle,         V0      V1      *
  137. *      relative to ObjMinZ level, is the sum of:                 *
  138. *      1.2.1. volume of V'(0), V'(i), V'(i+1) - the                 *
  139. *         area of projection of V(0), V(i), V(i+1) on XY parallel     *
  140. *         plane, times (TriMinZ - ObjMinZ).                 *
  141. *      1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     *
  142. *         V"(i+1) be the projections of V(i) and V(i+1) on the plane  *
  143. *         Z = PolyZMin. The volume above 1.2.1. and below the polygon *
  144. *         (triangle!) will be: the area of quadraliteral V(i), V(i+1),*
  145. *         V"(i+1), V"(i), times distance of V(0) for quadraliteral    *
  146. *         plane divided by 3.                         *
  147. *    1.3. If Z component of polygon normal is negative add 1.2. result to    *
  148. *      ObjectVolume, else subtract it.                     *
  149. *****************************************************************************/
  150. double GeomObjectVolume(ObjectStruct *PObj)
  151. {
  152.     int PlaneExists;
  153.     RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist;
  154.     PointType Pt1;
  155.     PlaneType Plane;
  156.     PolygonStruct *Pl;
  157.     VertexStruct *V, *VHead, *Vnext, *Vtemp;
  158.  
  159.     if (!IS_GEOM_OBJ(PObj))
  160.     FatalError("Geometric property requested on non geometric object");
  161.  
  162.     if (IS_POLYLINE_GEOM_OBJ(PObj)) {
  163.     WndwInputWindowPutStr("Warning: geometric object is polyline - zero area",
  164.                                     RED);
  165.     return 0.0;
  166.     }
  167.  
  168.     ObjMinZ = INFINITY;     /* Find Object minimum Z value (used as min level). */
  169.     Pl = PObj -> U.Pl;
  170.     while (Pl != NULL) {
  171.     V = VHead = Pl -> V;
  172.     do {
  173.         if (V -> Pt[2] < ObjMinZ) ObjMinZ = V -> Pt[2];
  174.         V = V -> Pnext;
  175.     }
  176.     while (V != VHead && V != NULL);
  177.     Pl = Pl -> Pnext;
  178.     }
  179.  
  180.     ConvexPolyObject(PObj);           /* Make sure all polygons are convex. */
  181.     Pl = PObj -> U.Pl;
  182.     while (Pl != NULL) {
  183.     PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */
  184.     V = Vtemp = VHead = Pl -> V;/* We set VHead to be vertex with min Z: */
  185.     do {
  186.         if (V -> Pt[2] < Vtemp -> Pt[2]) Vtemp = V;
  187.         V = V -> Pnext;
  188.     }
  189.     while (V != VHead && V != NULL);
  190.     VHead = Vtemp;       /* Now VHead is the one with lowest Z in polygon! */
  191.     TriMinZ = VHead -> Pt[2];         /* Save this Z for fast access. */
  192.  
  193.     V = VHead -> Pnext;
  194.     Vnext = V -> Pnext;
  195.     do {
  196.         /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */
  197.         Area = Polygon3VrtxXYArea(VHead -> Pt, V -> Pt, Vnext -> Pt);
  198.         PolygonVolume += Area * (TriMinZ - ObjMinZ);
  199.  
  200.         /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */
  201.         Area = sqrt(SQR(V -> Pt[0] - Vnext -> Pt[0]) +   /* XY distance. */
  202.             SQR(V -> Pt[1] - Vnext -> Pt[1])) *
  203.            ((V -> Pt[2] + Vnext -> Pt[2]) / 2.0 - TriMinZ);
  204.         PT_COPY(Pt1, V -> Pt);
  205.         Pt1[2] = TriMinZ;
  206.         if ((PlaneExists =
  207.          CGPlaneFrom3Points(Plane, V -> Pt, Vnext -> Pt, Pt1)) == 0) {
  208.         /* Try second pt projected to Z = TriMinZ plane as third pt. */
  209.         PT_COPY(Pt1, Vnext -> Pt);
  210.         Pt1[2] = TriMinZ;
  211.         PlaneExists =
  212.             CGPlaneFrom3Points(Plane, V -> Pt, Vnext -> Pt, Pt1);
  213.         }
  214.         if (PlaneExists) {
  215.         Dist = CGDistPointPlane(VHead -> Pt, Plane);
  216.         PolygonVolume += Area * ABS(Dist) / 3.0;
  217.         }
  218.  
  219.         V = Vnext;
  220.         Vnext = V -> Pnext;
  221.     }
  222.     while (Vnext != VHead);
  223.  
  224.     if (Pl -> Plane[2] < 0.0)
  225.          ObjVolume += PolygonVolume;
  226.     else ObjVolume -= PolygonVolume;
  227.  
  228.     Pl = Pl -> Pnext;
  229.     }
  230.  
  231.     return ObjVolume;
  232. }
  233.  
  234. /*****************************************************************************
  235. *  Routine to evaluate the area of the given triangle projected to the XY    *
  236. * plane, given as 3 Points. Only the X & Y components are considered.         *
  237. *  See GeomObjectArea above for algorithm:                     *
  238. *****************************************************************************/
  239. static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3)
  240. {
  241.     RealType PolygonArea = 0.0, MinY;
  242.  
  243.     MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0;
  244.  
  245.     PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0;
  246.     PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0;
  247.     PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0;
  248.  
  249.     return ABS(PolygonArea);
  250. }
  251.  
  252. /*****************************************************************************
  253. *  Routine to count number of polygons in given geometric object.         *
  254. *****************************************************************************/
  255. double GeomCountPolys(ObjectStruct *PObj)
  256. {
  257.     int Count = 0;
  258.     PolygonStruct *Pl;
  259.  
  260.     if (!IS_GEOM_OBJ(PObj))
  261.     FatalError("Geometric property requested on non geometric object");
  262.  
  263.     Pl = PObj -> U.Pl;
  264.     while (Pl != NULL) {
  265.     Count++;
  266.     Pl = Pl -> Pnext;
  267.     }
  268.  
  269.     return ((double) Count);
  270. }
  271.