Merge pull request #1215 from LLNL/feature/dayton8/tet_plane_clip_2

Adds a clip operator for a tet and the half-space defined by a plane

RELEASE-NOTES.md

@@ -34,6 +34,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Primal: Adds a `Quadrilateral` primitive
 - Primal: Adds a `compute_bounding_box()` operator for computing the bounding
   box of a `Quadrilateral`
+- Primal: Adds a `clip()` operator for clipping a tetrahedron against the
+  half-space defined by a plane
 
 ### Fixed
 - quest's `SamplingShaper` now properly handles material names containing underscores

src/axom/primal/operators/clip.hpp

@@ -304,6 +304,90 @@ AXOM_HOST_DEVICE Polyhedron<T, 3> clip(const Tetrahedron<T, 3>& tet1,
   return detail::clipTetrahedron(tet1, tet2, eps, checkSign);
 }
 
+/*!
+ * \brief Clips a 3D tetrahedron against the half-space defined by a plane
+ *        and returns the resulting polyhedron
+ *
+ *  This function clips a tetrahedron against the half-space defined by a
+ *  plane. This involves finding new vertices at the intersection of the
+ *  polyhedron edges and the plane, removing vertices from the polyhedron
+ *  that are below the plane, and redefining the neighbors for each vertex
+ *  (a vertex is a neighbor of another vertex if there is an edge between
+ *  them).
+ *
+ * \param [in] tet The tetrahedron to clip
+ * \param [in] plane The plane defining the half-space used to clip the tetrahedron
+ * \param [in] eps The tolerance for plane point orientation
+ * \param [in] tryFixOrientation If true and the tetrahedron has a negative
+ *             signed volume, swaps the order of some vertices in the
+ *             tetrathedron to try to obtain a nonnegative signed volume.
+ *             Defaults to false.
+ *
+ * \return The polyhedron obtained from clipping a tetrahedron against
+ *         the half-space defined by a plane.
+ *
+ * \note Function is based off clipPolyhedron() in Mike Owen's PolyClipper.
+ *
+ * \warning tryFixOrientation flag does not guarantee the tetrahedron's vertex
+ *          order will be valid. It is the responsiblity of the caller to pass
+ *          a tetrahedron with a valid vertex order. Otherwise, the returned
+ *          polyhedron will have a non-positive and/or unexpected volume.
+ *
+ * \warning If the tryFixOrientation flag is false and the tetrahedron has
+ *          a negative signed volume, the returned polyhedron will have a
+ *          non-positive and/or unexpected volume.
+ */
+template <typename T>
+AXOM_HOST_DEVICE Polyhedron<T, 3> clip(const Tetrahedron<T, 3>& tet,
+                                       const Plane<T, 3>& plane,
+                                       double eps = 1.e-10,
+                                       bool tryFixOrientation = false)
+{
+  return detail::clipTetrahedron(tet, plane, eps, tryFixOrientation);
+}
+
+/*!
+ * \brief Clips a 3D tetrahedron against the half-space defined by a plane
+ *        and returns the resulting polyhedron
+ *
+ *  This function clips a tetrahedron against the half-space defined by a
+ *  plane. This involves finding new vertices at the intersection of the
+ *  polyhedron edges and the plane, removing vertices from the polyhedron
+ *  that are below the plane, and redefining the neighbors for each vertex
+ *  (a vertex is a neighbor of another vertex if there is an edge between
+ *  them).
+ *
+ * \param [in] plane The plane defining the half-space used to clip the tetrahedron
+ * \param [in] tet The tetrahedron to clip
+ * \param [in] eps The tolerance for plane point orientation
+ * \param [in] tryFixOrientation If true and the tetrahedron has a negative
+ *             signed volume, swaps the order of some vertices in the
+ *             tetrathedron to try to obtain a nonnegative signed volume.
+ *             Defaults to false.
+ *
+ * \return The polyhedron obtained from clipping a tetrahedron against
+ *         the half-space defined by a plane.
+ *
+ * \note Function is based off clipPolyhedron() in Mike Owen's PolyClipper.
+ *
+ * \warning tryFixOrientation flag does not guarantee the tetrahedron's vertex
+ *          order will be valid. It is the responsiblity of the caller to pass
+ *          a tetrahedron with a valid vertex order. Otherwise, the returned
+ *          polyhedron will have a non-positive and/or unexpected volume.
+ *
+ * \warning If the tryFixOrientation flag is false and the tetrahedron has
+ *          a negative signed volume, the returned polyhedron will have a
+ *          non-positive and/or unexpected volume.
+ */
+template <typename T>
+AXOM_HOST_DEVICE Polyhedron<T, 3> clip(const Plane<T, 3>& plane,
+                                       const Tetrahedron<T, 3>& tet,
+                                       double eps = 1.e-10,
+                                       bool tryFixOrientation = false)
+{
+  return clip(tet, plane, eps, tryFixOrientation);
+}
+
 }  // namespace primal
 }  // namespace axom
 

src/axom/primal/operators/detail/clip_impl.hpp

@@ -400,83 +400,90 @@ AXOM_HOST_DEVICE void poly_clip_reindex(Polyhedron<T, NDIMS>& poly,
 }
 
 /*!
- * \brief Finds the clipped intersection Polyhedron between Polyhedron
- *        poly and an array of Planes.
- *
- * \param [in] poly The polyhedron
- * \param [in] planes The array of planes
- * \param [in] eps The tolerance for plane point orientation.
- *
- * \return The Polyhedron formed from clipping the polyhedron with a set of planes.
+ * \brief Clips a polyhedron against a half-space defined by a plane
  *
+ * \param [inout] poly The polyhedron to clip
+ * \param [in] plane The plane defining the half-space used to clip the polyhedron
+ * \param [in] eps The tolerance for plane point orientation
  */
 template <typename T, int NDIMS>
-AXOM_HOST_DEVICE Polyhedron<T, NDIMS> clipPolyhedron(
-  Polyhedron<T, NDIMS>& poly,
-  axom::ArrayView<Plane<T, NDIMS>> planes,
-  double eps)
+AXOM_HOST_DEVICE void clipPolyhedron(Polyhedron<T, NDIMS>& poly,
+                                     const Plane<T, NDIMS>& plane,
+                                     double eps)
 {
-  using PointType = Point<T, NDIMS>;
   using BoxType = BoundingBox<T, NDIMS>;
-  using PlaneType = Plane<T, NDIMS>;
-
-  //Bounding Box of Polyhedron
-  BoxType polyBox(&poly[0], poly.numVertices());
 
-  //Clip Polyhedron by each plane
-  for(PlaneType plane : planes)
+  // Check that plane intersects Polyhedron
+  if(intersect(plane, BoxType(&poly[0], poly.numVertices()), true, eps))
   {
-    // Check that plane intersects Polyhedron
-    if(intersect(plane, polyBox, true, eps))
-    {
-      int numVerts = poly.numVertices();
+    int numVerts = poly.numVertices();
 
-      // Each bit value indicates if that Polyhedron vertex is formed from
-      // Polyhedron clipping with a plane.
-      unsigned int clipped = 0;
+    // Each bit value indicates if that Polyhedron vertex is formed from
+    // Polyhedron clipping with a plane.
+    unsigned int clipped = 0;
 
-      // Clip polyhedron against current plane, generating extra vertices
-      // where edges meet the plane.
-      poly_clip_vertices(poly, plane, eps, clipped);
+    // Clip polyhedron against current plane, generating extra vertices
+    // where edges meet the plane.
+    poly_clip_vertices(poly, plane, eps, clipped);
 
-      // Adjust connectivity to link up newly-generated vertices.
-      poly_clip_fix_nbrs(poly, plane, numVerts, eps, clipped);
+    // Adjust connectivity to link up newly-generated vertices.
+    poly_clip_fix_nbrs(poly, plane, numVerts, eps, clipped);
 
-      // Reindex polyhedron connectivity by removing vertices on the negative
-      // side of the plane.
-      poly_clip_reindex(poly, clipped);
+    // Reindex polyhedron connectivity by removing vertices on the negative
+    // side of the plane.
+    poly_clip_reindex(poly, clipped);
+  }
 
-      // Generate new bounding box for polyhedron
-      polyBox = BoxType();
+  // If entire polyhedron is below a plane (points can be on the plane),
+  // it is completely removed.
+  else
+  {
+    bool completeClip = true;
 
-      for(int i = 0; i < poly.numVertices(); i++)
+    for(int i = 0; i < poly.numVertices(); i++)
+    {
+      if(plane.getOrientation(poly[i], eps) == ON_POSITIVE_SIDE)
       {
-        polyBox.addPoint(PointType(poly[i]));
+        completeClip = false;
+        break;
       }
     }
 
-    // If entire polyhedron is below a plane (points can be on the plane), it is completely removed.
-    else
+    if(completeClip)
     {
-      bool completeClip = true;
-      for(int i = 0; i < poly.numVertices(); i++)
-      {
-        if(plane.getOrientation(poly[i], eps) == ON_POSITIVE_SIDE)
-        {
-          completeClip = false;
-          break;
-        }
-      }
+      poly.clear();
+    }
+  }
 
-      if(completeClip)
-      {
-        poly.clear();
-        return poly;
-      }
+  return;
+}
+
+/*!
+ * \brief Clips a polyhedron against an array of planes
+ *
+ * \param [inout] poly The polyhedron to clip
+ * \param [in] planes The array of planes
+ * \param [in] eps The tolerance for plane point orientation
+ */
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE void clipPolyhedron(Polyhedron<T, NDIMS>& poly,
+                                     axom::ArrayView<Plane<T, NDIMS>> planes,
+                                     double eps)
+{
+  using PlaneType = Plane<T, NDIMS>;
+
+  // Clip Polyhedron by each plane
+  for(const PlaneType& plane : planes)
+  {
+    clipPolyhedron(poly, plane, eps);
+
+    if(poly.numVertices() == 0)
+    {
+      return;
     }
   }
 
-  return poly;
+  return;
 }
 
 /*!
@@ -523,7 +530,8 @@ AXOM_HOST_DEVICE Polyhedron<T, NDIMS> clipHexahedron(
   axom::StackArray<IndexType, 1> planeSize = {4};
   axom::ArrayView<PlaneType> planesView(planes, planeSize);
 
-  return clipPolyhedron(poly, planesView, eps);
+  clipPolyhedron(poly, planesView, eps);
+  return poly;
 }
 
 /*!
@@ -570,7 +578,8 @@ AXOM_HOST_DEVICE Polyhedron<T, NDIMS> clipOctahedron(
   axom::StackArray<IndexType, 1> planeSize = {4};
   axom::ArrayView<PlaneType> planesView(planes, planeSize);
 
-  return clipPolyhedron(poly, planesView, eps);
+  clipPolyhedron(poly, planesView, eps);
+  return poly;
 }
 
 /*!
@@ -617,7 +626,45 @@ AXOM_HOST_DEVICE Polyhedron<T, NDIMS> clipTetrahedron(
   axom::StackArray<IndexType, 1> planeSize = {4};
   axom::ArrayView<PlaneType> planesView(planes, planeSize);
 
-  return clipPolyhedron(poly, planesView, eps);
+  clipPolyhedron(poly, planesView, eps);
+  return poly;
+}
+
+/*!
+ * \brief Clips a tetrahedron against the half-space defined by a plane
+ *        and returns the resulting polyhedron
+ *
+ * \param [in] tet The tetrahedron to clip
+ * \param [in] plane The plane defining the half-space used to clip the tetrahedron
+ * \param [in] eps The tolerance for plane point orientation
+ * \param [in] tryFixOrientation If true and the tetrahedron has a negative
+ *             signed volume, swaps the order of some vertices in the
+ *             tetrathedron to try to obtain a nonnegative signed volume.
+ *             Defaults to false.
+ *
+ * \return The polyhedron obtained from clipping a tetrahedron against
+ *         the half-space defined a plane
+ *
+ * \warning tryFixOrientation flag does not guarantee the tetrahedron's vertex
+ *          order will be valid. It is the responsiblity of the caller to pass
+ *          a tetrahedron with a valid vertex order. Otherwise, the returned
+ *          polyhedron will have a non-positive and/or unexpected volume.
+ *
+ * \warning If the tryFixOrientation flag is false and the tetrahedron has
+ *          a negative signed volume, the returned polyhedron will have a
+ *          non-positive and/or unexpected volume.
+ */
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE Polyhedron<T, NDIMS> clipTetrahedron(
+  const Tetrahedron<T, NDIMS>& tet,
+  const Plane<T, NDIMS>& plane,
+  double eps,
+  bool tryFixOrientation)
+{
+  using PolyhedronType = Polyhedron<T, NDIMS>;
+  PolyhedronType poly = PolyhedronType::from_primitive(tet, tryFixOrientation);
+  clipPolyhedron(poly, plane, eps);
+  return poly;
 }
 
 }  // namespace detail