/****************************************************************************** * GeomVals.c - Area, Volume, and counts on polygonal objects. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, March 1990. * ******************************************************************************/ #include #include #include #include #include "allocate.h" #include "convex.h" #include "geomat3d.h" #include "geomvals.h" static RealType PolygonXYArea(IPVertexStruct *VHead); static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3); /***************************************************************************** * DESCRIPTION: M * Routine to evaluate the Area of the given geom. object, in object unit. M * Algorithm (for each polygon): V3 V * 1. Set Polygon Area to be zero. /\ V * Make a copy of the original polygon / \ V * and transform it to a XY parallel plane. / \V2 V * Find the minimum Y value of the polygon V4/ | V * in the XY plane. \ | V * 2. Let V(0) be the first vertex, V(n) the last one. \ | V * For i goes from 0 to n-1 add to Area the area \_______| V * below edge V(i), V(i+1): V0 V1 V * PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2 V * where V(i)y' is V(i)y - MinY, where MinY is polygon minimum Y value. V * 3. The result of step 2 is the area of the polygon itself. V * However, it might be negative, so take the absolute result of step 2 V * and add it to the global ObjectArea. V * Note step 2 is performed by another auxiliary routine: PolygonXYArea. M * * * PARAMETERS: M * PObj: A polyhedra object to compute its surface area. M * * * RETURN VALUE: M * double: The area of object PObj. M * * * KEYWORDS: M * PolyObjectArea, area M *****************************************************************************/ double PolyObjectArea(IPObjectStruct *PObj) { RealType ObjectArea = 0.0; IPPolygonStruct *Pl; if (!IP_IS_POLY_OBJ(PObj)) IritFatalError("Geometric property requested on non polygonal object"); if (IP_IS_POLYLINE_OBJ(PObj)) { return 0.0; } for (Pl = PObj -> U.Pl; Pl != NULL; Pl = Pl -> Pnext) ObjectArea += PolyOnePolyArea(Pl); return ObjectArea; } /***************************************************************************** * DESCRIPTION: M * Computes the area of a single polygon. M * * * PARAMETERS: M * Pl: To compute its area. M * * * RETURN VALUE: M * double: area of polygon Pl M * * * SEE ALSO: M * PolyObjectArea M * * * KEYWORDS: M * PolyOnePolyArea M *****************************************************************************/ double PolyOnePolyArea(IPPolygonStruct *Pl) { MatrixType RotMat; IPVertexStruct *VHead = CopyVertexList(Pl -> PVertex),/* Dont trans original object.*/ *V = VHead; RealType ObjectArea; /* Create the trans matrix to transform the polygon to XY parallel plane */ GenRotateMatrix(RotMat, Pl -> Plane); do { MatMultVecby4by4(V -> Coord, V -> Coord, RotMat); V = V -> Pnext; } while (V != NULL && V != VHead); ObjectArea = PolygonXYArea(VHead); IPFreeVertexList(VHead); /* Free the vertices list. */ return ObjectArea; } /***************************************************************************** * DESCRIPTION: * * Routine to evaluate the area of the given polygon projected on the XY * * plane. * * The polygon does not have to be on a XY parallel plane, as only its XY * * projection is considered (Z is ignored). * * Returned the area of its XY parallel projection. * * See GeomObjectArea above for algorithm. * * * * PARAMETERS: * * VHead: Of vertex list to compute area for. * * * * RETURN VALUE: * * RealType: Computed area. * *****************************************************************************/ static RealType PolygonXYArea(IPVertexStruct *VHead) { RealType PolygonArea = 0.0, MinY; IPVertexStruct *V = VHead, *Vnext; MinY = V -> Coord[1]; V = V -> Pnext; while (V != VHead && V != NULL) { if (MinY > V -> Coord[1]) MinY = V -> Coord[1]; V = V -> Pnext; } if (V == NULL) V = VHead; Vnext = V -> Pnext; MinY *= 2.0; /* Instead of subtracting twice each time. */ do { /* Evaluate area below edge V->Vnext relative to Y level MinY. Note */ /* it can come out negative, but thats o.k. as the sum of all these */ /* quadraliterals should be exactly the area (up to correct sign). */ PolygonArea += (Vnext -> Coord[0] - V -> Coord[0]) * (Vnext -> Coord[1] + V -> Coord[1] - MinY) / 2.0; V = Vnext; Vnext = V -> Pnext == NULL ? VHead : V -> Pnext; } while (V != VHead && V != NULL); return ABS(PolygonArea); } /***************************************************************************** * DESCRIPTION: M * Routine to evaluate the Volume of the given geom object, in object unit. M * This routine has a side effect that all non-convex polygons will be M * splitted to convex ones. M * Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y): M * V3 V * 1. Set Polygon Area to be zero. /\ V * Let V(0) be the first vertex, V(n) the last. / \ V * For i goes from 1 to n-1 form triangles / \V2 V * by V(0), V(i), V(i+1). V4/ | V * For each such triangle di: \ | V * 1.1. Find the vertex (out of V(0), V(i), V(i+1)) \ | V * with the minimum Z - TriMinY. \_______| V * 1.2. The volume below V(0), V(i), V(i+1) triangle, V0 V1 V * relative to ObjMinZ level, is the sum of: V * 1.2.1. volume of V'(0), V'(i), V'(i+1) - the V * area of projection of V(0), V(i), V(i+1) on XY parallel V * plane, times (TriMinZ - ObjMinZ). V * 1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and V * V"(i+1) be the projections of V(i) and V(i+1) on the plane V * Z = PolyZMin. The volume above 1.2.1. and below the polygon V * (triangle!) will be: the area of quadraliteral V(i), V(i+1),V * V"(i+1), V"(i), times distance of V(0) for quadraliteral V * plane divided by 3. V * 1.3. If Z component of polygon normal is negative add 1.2. result to V * ObjectVolume, else subtract it. V * * * PARAMETERS: M * PObj: To compute volume for. M * * * RETURN VALUE: M * double: Computed volume. M * * * KEYWORDS: M * PolyObjectVolume, volume M *****************************************************************************/ double PolyObjectVolume(IPObjectStruct *PObj) { int PlaneExists; RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist; PointType Pt1; PlaneType Plane; IPPolygonStruct *Pl; IPVertexStruct *V, *VHead, *Vnext, *Vtemp; if (!IP_IS_POLY_OBJ(PObj)) IritFatalError("Geometric property requested on non polygonal object"); if (IP_IS_POLYLINE_OBJ(PObj)) { return 0.0; } ObjMinZ = IRIT_INFNTY; /* Find Object minimum Z value (used as min level). */ Pl = PObj -> U.Pl; while (Pl != NULL) { V = VHead = Pl -> PVertex; do { if (V -> Coord[2] < ObjMinZ) ObjMinZ = V -> Coord[2]; V = V -> Pnext; } while (V != VHead && V != NULL); if (V == NULL) IritFatalError("For VOLUME polygons must be closed loops. Open loop detected.\n"); Pl = Pl -> Pnext; } ConvexPolyObject(PObj); /* Make sure all polygons are convex. */ Pl = PObj -> U.Pl; while (Pl != NULL) { PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */ /* We set VHead to be vertex with min Z: */ V = Vtemp = VHead = Pl -> PVertex; do { if (V -> Coord[2] < Vtemp -> Coord[2]) Vtemp = V; V = V -> Pnext; } while (V != VHead && V != NULL); VHead = Vtemp; /* Now VHead is the one with lowest Z in polygon! */ TriMinZ = VHead -> Coord[2]; /* Save this Z for fast access. */ V = VHead -> Pnext; Vnext = V -> Pnext; do { /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */ Area = Polygon3VrtxXYArea(VHead -> Coord, V -> Coord, Vnext -> Coord); PolygonVolume += Area * (TriMinZ - ObjMinZ); /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */ Area = sqrt(SQR(V -> Coord[0] - Vnext -> Coord[0]) + /* XY dist. */ SQR(V -> Coord[1] - Vnext -> Coord[1])) * ((V -> Coord[2] + Vnext -> Coord[2]) / 2.0 - TriMinZ); PT_COPY(Pt1, V -> Coord); Pt1[2] = TriMinZ; if ((PlaneExists = CGPlaneFrom3Points(Plane, V -> Coord, Vnext -> Coord, Pt1)) == 0) { /* Try second pt projected to Z = TriMinZ plane as third pt. */ PT_COPY(Pt1, Vnext -> Coord); Pt1[2] = TriMinZ; PlaneExists = CGPlaneFrom3Points(Plane, V -> Coord, Vnext -> Coord, Pt1); } if (PlaneExists) { Dist = CGDistPointPlane(VHead -> Coord, Plane); PolygonVolume += Area * ABS(Dist) / 3.0; } V = Vnext; Vnext = V -> Pnext; } while (Vnext != VHead); if (Pl -> Plane[2] < 0.0) ObjVolume += PolygonVolume; else ObjVolume -= PolygonVolume; Pl = Pl -> Pnext; } return ObjVolume; } /***************************************************************************** * DESCRIPTION: * * Routine to evaluate the area of the given triangle projected to the XY * * plane, given as 3 Points. Only the X & Y components are considered. * * See PolyObjectArea above for algorithm. * * * * PARAMETERS: * * Pt1, Pt2, Pt3: Vertices of triangle to compute area for. * * * * RETURN VALUE: * * RealType: Area of triangle in the XY plane. * *****************************************************************************/ static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3) { RealType PolygonArea = 0.0, MinY; MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0; PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0; PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0; PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0; return ABS(PolygonArea); } /***************************************************************************** * DESCRIPTION: M * Routine to count the number of polygons in the given geometric object. M * * * PARAMETERS: M * PObj: To count number of polygons in. M * * * RETURN VALUE: M * double: Number of polygons, an integer value returned as double. M * * * KEYWORDS: M * PolyCountPolys M *****************************************************************************/ double PolyCountPolys(IPObjectStruct *PObj) { int Count = 0; IPPolygonStruct *Pl; if (!IP_IS_POLY_OBJ(PObj)) IritFatalError("Geometric property requested on non polygonal object"); Pl = PObj -> U.Pl; while (Pl != NULL) { Count++; Pl = Pl -> Pnext; } return ((double) Count); }