Adds ray Bezier/NURBS surface intersection routines

src/axom/primal/CMakeLists.txt

@@ -26,6 +26,7 @@ set( primal_headers
     geometry/CurvedPolygon.hpp
     geometry/Hexahedron.hpp
     geometry/KnotVector.hpp
+    geometry/Line.hpp
     geometry/OrientedBoundingBox.hpp
     geometry/OrientationResult.hpp
     geometry/NumericArray.hpp
@@ -66,6 +67,7 @@ set( primal_headers
     operators/detail/fuzzy_comparators.hpp
     operators/detail/intersect_bezier_impl.hpp
     operators/detail/intersect_bounding_box_impl.hpp
+    operators/detail/intersect_patch_impl.hpp
     operators/detail/intersect_impl.hpp
     operators/detail/intersect_ray_impl.hpp
     operators/detail/winding_number_impl.hpp

src/axom/primal/geometry/BezierPatch.hpp

@@ -1867,10 +1867,11 @@ class BezierPatch
    * This function checks if all control points of the BezierPatch
    * are approximately on the plane defined by its four corners
    *
-   * \param [in] tol Threshold for sum of squared distances
+   * \param [in] sq_tol Threshold for sum of squared distances
+   * \param [in] EPS Threshold for nearness to zero
    * \return True if c1 is near-planar
    */
-  bool isPlanar(double tol = 1E-8) const
+  bool isPlanar(double sq_tol = 1E-8, double EPS = 1e-8) const
   {
     const int ord_u = getOrder_u();
     const int ord_v = getOrder_v();
@@ -1895,18 +1896,18 @@ class BezierPatch
     if(!axom::utilities::isNearlyEqual(
          VectorType::scalar_triple_product(v1, v2, v3),
          0.0,
-         tol))
+         EPS))
     {
       return false;
     }
 
     // Find three points that produce a nonzero normal
     Vector3D plane_normal = VectorType::cross_product(v1, v2);
-    if(axom::utilities::isNearlyEqual(plane_normal.norm(), 0.0, tol))
+    if(axom::utilities::isNearlyEqual(plane_normal.norm(), 0.0, EPS))
     {
       plane_normal = VectorType::cross_product(v1, v3);
     }
-    if(axom::utilities::isNearlyEqual(plane_normal.norm(), 0.0, tol))
+    if(axom::utilities::isNearlyEqual(plane_normal.norm(), 0.0, EPS))
     {
       plane_normal = VectorType::cross_product(v2, v3);
     }
@@ -1915,29 +1916,30 @@ class BezierPatch
     double sqDist = 0.0;
 
     // Check all control points for simplicity
-    for(int p = 0; p <= ord_u && sqDist <= tol; ++p)
+    for(int p = 0; p <= ord_u && sqDist <= sq_tol; ++p)
     {
-      for(int q = ((p == 0) ? 1 : 0); q <= ord_v && sqDist <= tol; ++q)
+      for(int q = ((p == 0) ? 1 : 0); q <= ord_v && sqDist <= sq_tol; ++q)
       {
         const double signedDist =
           plane_normal.dot(m_controlPoints(p, q) - m_controlPoints(0, 0));
         sqDist += signedDist * signedDist;
       }
     }
 
-    return (sqDist <= tol);
+    return (sqDist <= sq_tol);
   }
 
   /*!
    * \brief Predicate to check if the patch can be approximated by a polygon
    *
    * This function checks if a BezierPatch lies in a plane
-   *  and that the edged are linear up to tolerance `tol`
+   *  and that the edged are linear up to tolerance `sq_tol`
    *
    * \param [in] tol Threshold for sum of squared distances
+   * \param [in] EPS Threshold for nearness to zero
    * \return True if c1 is near-planar-polygonal
    */
-  bool isPolygonal(double tol = 1E-8) const
+  bool isPolygonal(double sq_tol = 1E-8, double EPS = 1e-8) const
   {
     const int ord_u = getOrder_u();
     const int ord_v = getOrder_v();
@@ -1956,32 +1958,95 @@ class BezierPatch
     }
 
     // Check if the patch is planar
-    if(!isPlanar(tol))
+    if(!isPlanar(sq_tol, EPS))
     {
       return false;
     }
 
     // Check if each bounding curve is linear
-    if(!isocurve_u(0).isLinear(tol))
+    if(!isocurve_u(0).isLinear(sq_tol))
     {
       return false;
     }
-    if(!isocurve_v(0).isLinear(tol))
+    if(!isocurve_v(0).isLinear(sq_tol))
     {
       return false;
     }
-    if(!isocurve_u(1).isLinear(tol))
+    if(!isocurve_u(1).isLinear(sq_tol))
     {
       return false;
     }
-    if(!isocurve_v(1).isLinear(tol))
+    if(!isocurve_v(1).isLinear(sq_tol))
     {
       return false;
     }
 
     return true;
   }
 
+  /*!
+   * \brief Predicate to check if the Bezier patch is approximately bilinear
+   *
+   * This function checks if the patch is (nearly) bilinear in terms of its shape
+   *  and it's parameterization. A necessary (but possibly not sufficient) condition
+   *  for this is that the control points are coincident with the surface of the
+   *  bilinear patch defined by its corners evaluated at uniform parameter values.   *  
+   *
+   * \param [in] sq_tol Threshold for absolute squared distances
+   * \return True if patch is bilinear
+   */
+  bool isBilinear(double sq_tol = 1e-8) const
+  {
+    const int ord_u = getOrder_u();
+    const int ord_v = getOrder_v();
+
+    if(ord_u <= 1 && ord_v <= 1)
+    {
+      return true;
+    }
+
+    // Anonymous function to evaluate the bilinear patch defined by the corners
+    auto bilinear_patch = [&](T u, T v) -> PointType {
+      PointType val;
+      for(int N = 0; N < NDIMS; ++N)
+      {
+        val[N] = axom::utilities::lerp(
+          axom::utilities::lerp(m_controlPoints(0, 0)[N],
+                                m_controlPoints(0, ord_v)[N],
+                                v),
+          axom::utilities::lerp(m_controlPoints(ord_u, 0)[N],
+                                m_controlPoints(ord_u, ord_v)[N],
+                                v),
+          u);
+      }
+      return val;
+    };
+
+    for(int u = 0; u <= ord_u; ++u)
+    {
+      for(int v = 0; v <= ord_v; ++v)
+      {
+        // Don't need to check the corners
+        if((u == 0 && v == 0) || (u == 0 && v == ord_v) ||
+           (u == ord_u && v == 0) || (u == ord_u && v == ord_v))
+        {
+          continue;
+        }
+
+        // Evaluate where the control point would be if the patch *was* bilinear
+        PointType bilinear_point =
+          bilinear_patch(u / static_cast<T>(ord_u), v / static_cast<T>(ord_v));
+
+        if(squared_distance(m_controlPoints(u, v), bilinear_point) > sq_tol)
+        {
+          return false;
+        }
+      }
+    }
+
+    return true;
+  }
+
   /*!
    * \brief Simple formatted print of a Bezier Patch instance
    *

src/axom/primal/geometry/KnotVector.hpp

@@ -191,6 +191,21 @@ class KnotVector
   /// \brief Return the number of knots in the knot vector
   axom::IndexType getNumKnots() const { return m_knots.size(); }
 
+  /// \brief Return an array of the unique knot values
+  axom::Array<T> getUniqueKnots() const
+  {
+    axom::Array<T> unique_knots;
+    for(int i = 0; i < m_knots.size(); ++i)
+    {
+      if(i == 0 || m_knots[i] != m_knots[i - 1])
+      {
+        unique_knots.push_back(m_knots[i]);
+      }
+    }
+
+    return unique_knots;
+  }
+
   /// \brief Return the number of valid knot spans
   axom::IndexType getNumKnotSpans() const
   {

src/axom/primal/geometry/Line.hpp

@@ -0,0 +1,169 @@
+// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#ifndef AXOM_PRIMAL_LINE_HPP_
+#define AXOM_PRIMAL_LINE_HPP_
+
+#include "axom/primal/geometry/Point.hpp"
+#include "axom/primal/geometry/Segment.hpp"
+#include "axom/primal/geometry/Vector.hpp"
+
+#include "axom/slic/interface/slic.hpp"
+
+#include <ostream>
+
+namespace axom
+{
+namespace primal
+{
+// Forward declare the templated classes and operator functions
+template <typename T, int NDIMS>
+class Line;
+
+/*!
+ * \brief Overloaded output operator for lines
+ */
+template <typename T, int NDIMS>
+std::ostream& operator<<(std::ostream& os, const Line<T, NDIMS>& line);
+
+/*!
+ * \class Line
+ *
+ * \brief Represents a line, \f$ L(t) \in \mathcal{R}^d \f$ , defined by an
+ *  origin point, \f$ P \f$ and a normalized direction vector, \f$ \vec{d} \f$,
+ *  \f$ \ni L(t)= P + t\vec{d} \forall t \in \mathcal{R} \f$
+ * 
+ * \tparam T the coordinate type, e.g., double, float, etc.
+ * \tparam NDIMS the number of dimensions
+ */
+template <typename T, int NDIMS>
+class Line
+{
+public:
+  using CoordType = T;
+  using PointType = Point<T, NDIMS>;
+  using SegmentType = Segment<T, NDIMS>;
+  using VectorType = Vector<T, NDIMS>;
+
+public:
+  // Disable the default constructor
+  Line() = delete;
+
+  /*!
+     * \brief Constructs a line object with the given origin and direction.
+     * \param [in] origin the origin of the line.
+     * \param [in] direction the direction of the line.
+     * \pre direction.squared_norm()!= 0.0
+     */
+  AXOM_HOST_DEVICE
+  Line(const PointType& origin, const VectorType& direction);
+
+  /*!
+    * \brief Constructs a line object from a directed segment.
+    * \params [in] S user-supplied segment 
+    */
+  explicit Line(const SegmentType& S);
+
+  /*!
+    * \brief Returns the origin of this Line instance.
+    * \return origin a point instance corresponding to the origin of the line.
+    */
+  AXOM_HOST_DEVICE
+  const PointType& origin() const { return m_origin; };
+
+  /*!
+    * \brief Returns a point along the line by evaluating \f$ L(t) \f$
+    * \param [in] t user-supplied value for L(t).
+    * \return p a point along the line.
+    */
+  AXOM_HOST_DEVICE
+  PointType at(const T& t) const;
+
+  /*!
+    * \brief Returns the direction vector of this Line instance.
+    * \return direction the direction vector of the line.
+    * \post direction.norm()==1
+    */
+  AXOM_HOST_DEVICE
+  const VectorType& direction() const { return m_direction; };
+
+  /*!
+    * \brief Simple formatted print of a line instance
+    * \param os The output stream to write to
+    * \return A reference to the modified ostream
+    */
+  std::ostream& print(std::ostream& os) const
+  {
+    os << "{origin:" << m_origin << "; direction:" << m_direction << "}";
+
+    return os;
+  }
+
+private:
+  PointType m_origin;
+  VectorType m_direction;
+};
+
+} /* namespace primal */
+} /* namespace axom */
+
+//------------------------------------------------------------------------------
+// Line Implementation
+//------------------------------------------------------------------------------
+
+namespace axom
+{
+namespace primal
+{
+//------------------------------------------------------------------------------
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE Line<T, NDIMS>::Line(const PointType& origin,
+                                      const VectorType& direction)
+  : m_origin(origin)
+  , m_direction(direction.unitVector())
+{
+  SLIC_ASSERT(m_direction.squared_norm() != 0.0);
+}
+
+//------------------------------------------------------------------------------
+template <typename T, int NDIMS>
+Line<T, NDIMS>::Line(const SegmentType& S)
+  : m_origin(S.source())
+  , m_direction(VectorType(S.source(), S.target()).unitVector())
+{
+  SLIC_ASSERT(m_direction.squared_norm() != 0.0);
+}
+
+//------------------------------------------------------------------------------
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE inline Point<T, NDIMS> Line<T, NDIMS>::at(const T& t) const
+{
+  PointType p;
+  for(int i = 0; i < NDIMS; ++i)
+  {
+    p[i] = m_origin[i] + t * m_direction[i];
+  }
+  return (p);
+}
+
+//------------------------------------------------------------------------------
+/// Free functions implementing Line's operators
+//------------------------------------------------------------------------------
+template <typename T, int NDIMS>
+std::ostream& operator<<(std::ostream& os, const Line<T, NDIMS>& line)
+{
+  line.print(os);
+  return os;
+}
+
+}  // namespace primal
+}  // namespace axom
+
+/// Overload to format a primal::Line using fmt
+template <typename T, int NDIMS>
+struct axom::fmt::formatter<axom::primal::Line<T, NDIMS>> : ostream_formatter
+{ };
+
+#endif  // AXOM_PRIMAL_LINE_HPP_