Merge pull request #1251 from LLNL/bugfix/kweiss/clip-triangle-bbox

Bugfix for primal's clip(triangle,bbox)

RELEASE-NOTES.md

@@ -50,12 +50,14 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
   returns the signed volume.
 - Primal: `intersection_volume()` operators changed from returning a signed
   volume to an unsigned volume.
+- Primal's `BoundingBox::contains(BoundingBox)`  now returns `true` when the input is empty
 
 ### Fixed
 - quest's `SamplingShaper` now properly handles material names containing underscores
 - quest's `SamplingShaper` can now be used with an mfem that is configured for (GPU) devices
 - primal's `Polygon` area computation in 3D previously only worked when the polygon was aligned with the XY-plane. It now works for arbitrary polygons.
 - Upgrades our `vcpkg` usage for automated Windows builds of our TPLs to its [2023.12.12 release](https://github.com/microsoft/vcpkg/releases/tag/2023.12.12)
+- Fixed a bug in the bounds checks for `primal::clip(Triangle, BoundingBox)`
 
 ## [Version 0.8.1] - Release date 2023-08-16
 

src/axom/mint/mesh/StructuredMesh.cpp

@@ -57,6 +57,7 @@ void StructuredMesh::setExtent(int ndims, const int64* extent)
   }
 
   SLIC_ASSERT(numNodes == getNumberOfNodes());
+  AXOM_UNUSED_VAR(numNodes);
 
 #ifdef AXOM_MINT_USE_SIDRE
   if(hasSidreGroup())

src/axom/multimat/examples/traversal.cpp

@@ -168,6 +168,8 @@ void various_traversal_methods(int nmats,
   timer.stop();
   SLIC_INFO("  Field1D: " << timer.elapsed() << " sec");
   SLIC_ASSERT(c_sum == sum);
+  AXOM_UNUSED_VAR(c_sum);
+  AXOM_UNUSED_VAR(sum);
 
   sum = 0;
   timer.reset();
@@ -203,6 +205,7 @@ void various_traversal_methods(int nmats,
   timer.stop();
   SLIC_INFO("  Field2D: " << timer.elapsed() << " sec");
   SLIC_ASSERT(x_sum == sum);
+  AXOM_UNUSED_VAR(x_sum);
 
   // ------- Dense Access ----------
   SLIC_INFO("\n -- Dense Access via map-- ");

src/axom/multimat/multimat.hpp

@@ -1055,21 +1055,21 @@ int MultiMat::addFieldArray_impl(const std::string& field_name,
   SLIC_ASSERT(m_fieldNameVec.size() == m_fieldSparsityLayoutVec.size());
   SLIC_ASSERT(m_fieldNameVec.size() == m_fieldStrideVec.size());
 
-  axom::IndexType set_size = 0;
+#ifdef AXOM_DEBUG
   if(field_mapping == FieldMapping::PER_CELL_MAT)
   {
     const BivariateSetType* s = get_mapped_biSet(data_layout, sparsity_layout);
     SLIC_ASSERT(s != nullptr);
-    set_size = s->size();
+    SLIC_ASSERT(s->size() * stride == data_arr.size());
   }
   else
   {
     SLIC_ASSERT(field_mapping == FieldMapping::PER_CELL ||
                 field_mapping == FieldMapping::PER_MAT);
     const RangeSetType& s = *getMappedRangeSet(field_mapping);
-    set_size = s.size();
+    SLIC_ASSERT(s.size() * stride == data_arr.size());
   }
-  SLIC_ASSERT(set_size * stride == data_arr.size());
+#endif
 
   m_fieldBackingVec.back() =
     std::make_unique<FieldBacking>(data_arr, owned, m_fieldAllocatorId);

src/axom/primal/CMakeLists.txt

@@ -60,6 +60,7 @@ set( primal_headers
 
     operators/detail/clip_impl.hpp
     operators/detail/compute_moments_impl.hpp
+    operators/detail/fuzzy_comparators.hpp
     operators/detail/intersect_bezier_impl.hpp
     operators/detail/intersect_bounding_box_impl.hpp
     operators/detail/intersect_impl.hpp

src/axom/primal/geometry/BoundingBox.hpp

@@ -255,6 +255,7 @@ class BoundingBox
    * \note We are allowing the other bounding box to have a different coordinate
    *  type. This should work as long as the two Ts are comparable with
    *  operator<().
+   * \note If \a otherBB is empty, we return true
    */
   template <typename OtherType>
   bool contains(const BoundingBox<OtherType, NDIMS>& otherBB) const;
@@ -438,7 +439,9 @@ template <typename T, int NDIMS>
 template <typename OtherT>
 bool BoundingBox<T, NDIMS>::contains(const BoundingBox<OtherT, NDIMS>& otherBB) const
 {
-  return this->contains(otherBB.getMin()) && this->contains(otherBB.getMax());
+  return otherBB.isValid()
+    ? this->contains(otherBB.getMin()) && this->contains(otherBB.getMax())
+    : true;
 }
 
 //------------------------------------------------------------------------------

src/axom/primal/geometry/OrientedBoundingBox.hpp

@@ -181,10 +181,9 @@ class OrientedBoundingBox
   /*!
    * \brief Expands the box so it contains the passed in box.
    * \param obb [in] OBB to add in
-   * \note If obb is invalid, makes this invalid too.
-   * \note Expands the extents so this box contains obb. This loses
-   * optimality of the box; for a better fit use merge_boxes in
-   * primal/compute_bounding_box.
+   * \note If obb is invalid, this is a no-op
+   * \note Expands the extents so this box contains obb. This loses optimality of the box
+   * for a better fit use merge_boxes in primal/compute_bounding_box.
    * \post this->contains(obb) == true
    */
   void addBox(OrientedBoxType obb);
@@ -518,19 +517,15 @@ void OrientedBoundingBox<T, NDIMS>::addPoint(Point<T, NDIMS> pt)
 template <typename T, int NDIMS>
 void OrientedBoundingBox<T, NDIMS>::addBox(OrientedBoundingBox<T, NDIMS> obb)
 {
-  SLIC_CHECK_MSG(obb.isValid(), "Passed in OBB is invalid.");
   if(!obb.isValid())
   {
     // don't do anything
     return;
   }
 
-  std::vector<Point<T, NDIMS>> res = obb.vertices();
-  int size = static_cast<int>(res.size());
-
-  for(int i = 0; i < size; i++)
+  for(const auto& vert : obb.vertices())
   {
-    this->addPoint(res[i]);
+    this->addPoint(vert);
   }
 }
 

src/axom/primal/geometry/Ray.hpp

@@ -162,4 +162,9 @@ std::ostream& operator<<(std::ostream& os, const Ray<T, NDIMS>& ray)
 }  // namespace primal
 }  // namespace axom
 
+/// Overload to format a primal::Ray using fmt
+template <typename T, int NDIMS>
+struct axom::fmt::formatter<axom::primal::Ray<T, NDIMS>> : ostream_formatter
+{ };
+
 #endif  // AXOM_PRIMAL_RAY_HPP_

src/axom/primal/operators/clip.hpp

@@ -50,16 +50,13 @@ Polygon<T, 3> clip(const Triangle<T, 3>& tri, const BoundingBox<T, 3>& bbox)
   using PolygonType = Polygon<T, 3>;
 
   // Use two polygons with pointers for 'back-buffer'-like swapping
-  const int MAX_VERTS = 6;
+  constexpr int MAX_VERTS = 6;
   PolygonType poly[2] = {PolygonType(MAX_VERTS), PolygonType(MAX_VERTS)};
   PolygonType* currentPoly = &poly[0];
   PolygonType* prevPoly = &poly[1];
 
   // First check if the triangle is contained in the bbox, if not we are empty
-  BoundingBoxType triBox;
-  triBox.addPoint(tri[0]);
-  triBox.addPoint(tri[1]);
-  triBox.addPoint(tri[2]);
+  BoundingBoxType triBox {tri[0], tri[1], tri[2]};
 
   if(!bbox.intersectsWith(triBox))
   {
@@ -82,14 +79,13 @@ Polygon<T, 3> clip(const Triangle<T, 3>& tri, const BoundingBox<T, 3>& bbox)
   {
     // Optimization note: we should be able to save some work based on
     // the clipping plane and the triangle's bounding box
-
-    if(triBox.getMax()[dim] > bbox.getMin()[dim])
+    if(triBox.getMin()[dim] < bbox.getMin()[dim])
     {
       axom::utilities::swap(prevPoly, currentPoly);
       detail::clipAxisPlane(prevPoly, currentPoly, 2 * dim + 0, bbox.getMin()[dim]);
     }
 
-    if(triBox.getMin()[dim] < bbox.getMax()[dim])
+    if(triBox.getMax()[dim] > bbox.getMax()[dim])
     {
       axom::utilities::swap(prevPoly, currentPoly);
       detail::clipAxisPlane(prevPoly, currentPoly, 2 * dim + 1, bbox.getMax()[dim]);

src/axom/primal/operators/compute_bounding_box.hpp

@@ -22,6 +22,7 @@
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/primal/geometry/Triangle.hpp"
 #include "axom/primal/geometry/Quadrilateral.hpp"
+#include "axom/primal/geometry/Polygon.hpp"
 #include "axom/primal/geometry/Vector.hpp"
 #include "axom/primal/geometry/OrientedBoundingBox.hpp"
 
@@ -169,8 +170,8 @@ template <typename T, int NDIMS>
 AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
   const Polyhedron<T, NDIMS> &poly)
 {
-  BoundingBox<T, NDIMS> res(poly[0]);
-  for(int i = 1; i < poly.numVertices(); i++)
+  BoundingBox<T, NDIMS> res;
+  for(int i = 0; i < poly.numVertices(); i++)
   {
     res.addPoint(poly[i]);
   }
@@ -180,7 +181,7 @@ AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
 /*!
  * \brief Creates a bounding box around a Tetrahedron
  *
- * \param [in] tet The Tetrahedron
+ * \param [in] tet The tetrahedron
  */
 template <typename T, int NDIMS>
 AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
@@ -189,6 +190,23 @@ AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
   return BoundingBox<T, NDIMS> {tet[0], tet[1], tet[2], tet[3]};
 }
 
+/*!
+ * \brief Creates a bounding box around a Polygon
+ *
+ * \param [in] poly The polygon
+ */
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE BoundingBox<T, NDIMS> compute_bounding_box(
+  const Polygon<T, NDIMS> &poly)
+{
+  BoundingBox<T, NDIMS> res;
+  for(int i = 0; i < poly.numVertices(); ++i)
+  {
+    res.addPoint(poly[i]);
+  }
+  return res;
+}
+
 }  // namespace primal
 }  // namespace axom
 

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

@@ -68,7 +68,7 @@ int classifyPointAxisPlane(const Point<T, NDIMS>& pt,
   // Note: we are exploiting the fact that the planes are axis aligned
   // So the dot product is +/- the given coordinate.
   // In general, we would need to call distance(pt, plane) here
-  T dist = isEven(index) ? val - pt[index / 2] : pt[index / 2] - val;
+  const T dist = isEven(index) ? val - pt[index / 2] : pt[index / 2] - val;
 
   if(dist > eps)
   {
@@ -106,10 +106,11 @@ Point<T, NDIMS> findIntersectionPoint(const Point<T, NDIMS>& a,
   // * pt = a + t * (b-a)
   // * pt[ index/2]  == val
 
-  T t = (val - a[index / 2]) / (b[index / 2] - a[index / 2]);
+  const int coord = index / 2;
+  const T t = (val - a[coord]) / (b[coord] - a[coord]);
   SLIC_ASSERT(0. <= t && t <= 1.);
 
-  PointType ret = PointType(a.array() + t * (b.array() - a.array()));
+  const auto ret = PointType::lerp(a, b, t);
   SLIC_ASSERT(classifyPointAxisPlane(ret, index, val) == ON_BOUNDARY);
 
   return ret;
@@ -142,7 +143,7 @@ void clipAxisPlane(const Polygon<T, NDIMS>* prevPoly,
   using PointType = Point<T, NDIMS>;
 
   currentPoly->clear();
-  int numVerts = prevPoly->numVertices();
+  const int numVerts = prevPoly->numVertices();
 
   if(numVerts == 0)
   {
@@ -156,7 +157,7 @@ void clipAxisPlane(const Polygon<T, NDIMS>* prevPoly,
   for(int i = 0; i < numVerts; ++i)
   {
     const PointType* b = &(*prevPoly)[i];
-    int bSide = classifyPointAxisPlane(*b, index, val);
+    const int bSide = classifyPointAxisPlane(*b, index, val);
 
     switch(bSide)
     {

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

@@ -0,0 +1,96 @@
+// Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/*!
+ * \file fuzzy_comparators.hpp
+ *
+ * This file provides helper functions for fuzzy comparisons
+ */
+
+#ifndef AXOM_PRIMAL_FUZZY_COMPARATORS_HPP_
+#define AXOM_PRIMAL_FUZZY_COMPARATORS_HPP_
+
+#include "axom/core/Macros.hpp"
+#include "axom/core/utilities/Utilities.hpp"
+
+namespace axom
+{
+namespace primal
+{
+namespace detail
+{
+/// \brief Checks if x > y, within a specified tolerance.
+AXOM_HOST_DEVICE
+inline bool isGt(double x, double y, double EPS = 1e-12)
+{
+  return ((x > y) && !(axom::utilities::isNearlyEqual(x, y, EPS)));
+}
+
+/// \brief Checks if x < y, within a specified tolerance.
+AXOM_HOST_DEVICE
+inline bool isLt(double x, double y, double EPS = 1e-12)
+{
+  return ((x < y) && !(axom::utilities::isNearlyEqual(x, y, EPS)));
+}
+
+/// \brief Checks if x <= y, within a specified tolerance.
+AXOM_HOST_DEVICE
+inline bool isLeq(double x, double y, double EPS = 1e-12)
+{
+  return !(isGt(x, y, EPS));
+}
+
+/// \brief Checks if x >= y, within a specified tolerance.
+AXOM_HOST_DEVICE
+inline bool isGeq(double x, double y, double EPS = 1e-12)
+{
+  return !(isLt(x, y, EPS));
+}
+
+/*!
+ * \brief Checks if x < y, or possibly x == y, within a specified tolerance.
+ *
+ * The check for equality is controlled by parameter includeEqual.  This
+ * lets users specify at compile time whether triangles intersecting only on
+ * border points are reported as intersecting or not.
+ *
+ * Supports checkEdge and checkVertex
+ */
+AXOM_HOST_DEVICE
+inline bool isLpeq(double x, double y, bool includeEqual = false, double EPS = 1e-12)
+{
+  if(includeEqual && axom::utilities::isNearlyEqual(x, y, EPS))
+  {
+    return true;
+  }
+
+  return isLt(x, y, EPS);
+}
+
+/*!
+ * \brief Checks if x > y, or possibly x == y, within a specified tolerance.
+ *
+ * The check for equality is controlled by parameter includeEqual.  This
+ * lets users specify at compile time whether triangles intersecting only on
+ * border points are reported as intersecting or not.
+ *
+ * Supports checkEdge and checkVertex
+ */
+AXOM_HOST_DEVICE
+inline bool isGpeq(double x, double y, bool includeEqual = false, double EPS = 1e-12)
+{
+  if(includeEqual && axom::utilities::isNearlyEqual(x, y, EPS))
+  {
+    return true;
+  }
+
+  return isGt(x, y, EPS);
+}
+
+}  // namespace detail
+}  // namespace primal
+}  // namespace axom
+
+#endif  // AXOM_PRIMAL_FUZZY_COMPARATORS_HPP_