Update 2D GWN algorithm to match paper

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

@@ -31,7 +31,7 @@ namespace primal
 namespace detail
 {
 /*
- * \brief Compute the winding number with respect to a line segment
+ * \brief Compute the generalized winding number with respect to a line segment
  *
  * \param [in] q The query point to test
  * \param [in] c0 The initial point of the line segment
@@ -42,7 +42,7 @@ namespace detail
  * is the signed angle subtended by the query point to each endpoint.
  * Colinear points return 0 for their winding number.
  *
- * \return double The winding number
+ * \return double The GWN
  */
 template <typename T>
 double linear_winding_number(const Point<T, 2>& q,
@@ -176,54 +176,55 @@ double convex_endpoint_winding_number(const Point<T, 2>& q,
 }
 
 /*!
- * \brief Recursively compute the winding number for a query point with respect
- *        to a single Bezier curve.
+ * \brief Recursively construct a polygon with the same *integer* winding number 
+ *        as the closed original Bezier curve at a given query point.
  *
  * \param [in] q The query point at which to compute winding number
- * \param [in] c The BezierCurve object along which to compute the winding number
- * \param [in] isConvexControlPolygon Boolean flag if the input Bezier curve 
+ * \param [in] c A BezierCurve subcurve of the curve along which to compute the winding number
+ * \param [in] isConvexControlPolygon Boolean flag if the input Bezier subcurve 
                                       is already convex
  * \param [in] edge_tol The physical distance level at which objects are 
  *                      considered indistinguishable
  * \param [in] EPS Miscellaneous numerical tolerance for nonphysical distances, used in
  *   isLinear, isNearlyZero, in_polygon, is_convex
+ * \param [out] approximating_polygon The Polygon that, by termination of recursion,
+ *   has the same integer winding number as the original closed curve
+ * \param [out] endpoint_gwn A running sum for the exact GWN if the point is at the 
+ *   endpoint of a subcurve
+ *
+ * By the termination of the recursive algorithm, `approximating_polygon` contains
+ *  a polygon that has the same *integer* winding number as the original curve.
  * 
- * Use a recursive algorithm that checks if the query point is exterior to
- * some convex shape containing the Bezier curve, in which case we have a direct 
- * formula for the winding number. If not, we bisect our curve and run the algorithm on 
- * each half. Use the proximity of the query point to endpoints and approximate
- * linearity of the Bezier curve as base cases.
- * 
- * \return double The winding number.
+ * Upon entering this algorithm, the closing line of `c` is already an
+ *  edge of the approximating polygon.
+ * If q is outside a convex shape that contains the entire curve, the 
+ *  integer winding number for the *closed* curve `c` is zero, 
+ *  and the algorithm terminates.
+ * If the shape is not convex or we're inside it, instead add the midpoint 
+ *  as a vertex and repeat the algorithm. 
  */
 template <typename T>
-double curve_winding_number_recursive(const Point<T, 2>& q,
-                                      const BezierCurve<T, 2>& c,
-                                      bool isConvexControlPolygon,
-                                      double edge_tol = 1e-8,
-                                      double EPS = 1e-8)
+void construct_approximating_polygon(const Point<T, 2>& q,
+                                     const BezierCurve<T, 2>& c,
+                                     bool isConvexControlPolygon,
+                                     double edge_tol,
+                                     double EPS,
+                                     Polygon<T, 2>& approximating_polygon,
+                                     double& endpoint_gwn,
+                                     bool& isCoincident)
 {
   const int ord = c.getOrder();
-  if(ord <= 0)
-  {
-    return 0.0;  // Catch degenerate cases
-  }
-
-  // If q is outside a convex shape that contains the entire curve, the winding
-  //   number for the shape connected at the endpoints with straight lines is zero.
-  //   We then subtract the contribution of this line segment.
 
   // Simplest convex shape containing c is its bounding box
-  BoundingBox<T, 2> bBox(c.boundingBox());
-  if(!bBox.contains(q))
+  if(!c.boundingBox().contains(q))
   {
-    return 0.0 - linear_winding_number(q, c[ord], c[0], edge_tol);
+    return;
   }
 
-  // Use linearity as base case for recursion.
+  // Use linearity as base case for recursion
   if(c.isLinear(EPS))
   {
-    return linear_winding_number(q, c[0], c[ord], edge_tol);
+    return;
   }
 
   // Check if our control polygon is convex.
@@ -236,28 +237,51 @@ double curve_winding_number_recursive(const Point<T, 2>& q,
   {
     isConvexControlPolygon = is_convex(controlPolygon, EPS);
   }
-  else  // Formulas for winding number only work if shape is convex
+
+  // Formulas for winding number only work if shape is convex
+  if(isConvexControlPolygon)
   {
     // Bezier curves are always contained in their convex control polygon
     if(!in_polygon(q, controlPolygon, includeBoundary, useNonzeroRule, EPS))
     {
-      return 0.0 - linear_winding_number(q, c[ord], c[0], edge_tol);
+      return;
     }
 
-    // If the query point is at either endpoint, use direct formula
-    if((squared_distance(q, c[0]) <= edge_tol * edge_tol) ||
-       (squared_distance(q, c[ord]) <= edge_tol * edge_tol))
+    // If the query point is at either endpoint...
+    if(squared_distance(q, c[0]) <= edge_tol * edge_tol ||
+       squared_distance(q, c[ord]) <= edge_tol * edge_tol)
     {
-      return convex_endpoint_winding_number(q, c, edge_tol, EPS);
+      // ...we can use a direct formula for the GWN at the endpoint
+      endpoint_gwn += convex_endpoint_winding_number(q, c, edge_tol, EPS);
+      isCoincident = true;
+
+      return;
     }
   }
 
   // Recursively split curve until query is outside some known convex region
   BezierCurve<T, 2> c1, c2;
   c.split(0.5, c1, c2);
 
-  return curve_winding_number_recursive(q, c1, isConvexControlPolygon, edge_tol, EPS) +
-    curve_winding_number_recursive(q, c2, isConvexControlPolygon, edge_tol, EPS);
+  construct_approximating_polygon(q,
+                                  c1,
+                                  isConvexControlPolygon,
+                                  edge_tol,
+                                  EPS,
+                                  approximating_polygon,
+                                  endpoint_gwn,
+                                  isCoincident);
+  approximating_polygon.addVertex(c2[0]);
+  construct_approximating_polygon(q,
+                                  c2,
+                                  isConvexControlPolygon,
+                                  edge_tol,
+                                  EPS,
+                                  approximating_polygon,
+                                  endpoint_gwn,
+                                  isCoincident);
+
+  return;
 }
 
 /// Type to indicate orientation of singularities relative to surface

src/axom/primal/operators/winding_number.hpp

@@ -27,6 +27,7 @@
 #include "axom/primal/geometry/CurvedPolygon.hpp"
 #include "axom/primal/geometry/BoundingBox.hpp"
 #include "axom/primal/geometry/OrientedBoundingBox.hpp"
+
 #include "axom/primal/operators/detail/winding_number_impl.hpp"
 
 // C++ includes
@@ -109,16 +110,19 @@ int winding_number(const Point<T, 2>& q,
 template <typename T>
 int winding_number(const Point<T, 2>& R,
                    const Polygon<T, 2>& P,
+                   bool& isOnEdge,
                    bool includeBoundary = false,
                    double EPS = 1e-8)
 {
   const int nverts = P.numVertices();
+  isOnEdge = false;
 
   // If the query is a vertex, return a value interpreted
   //  as "inside" by evenodd or nonzero protocols
   if(axom::utilities::isNearlyEqual(P[0][0], R[0], EPS) &&
      axom::utilities::isNearlyEqual(P[0][1], R[1], EPS))
   {
+    isOnEdge = true;
     return includeBoundary;
   }
 
@@ -131,10 +135,12 @@ int winding_number(const Point<T, 2>& R,
     {
       if(axom::utilities::isNearlyEqual(P[j][0], R[0], EPS))
       {
+        isOnEdge = true;
         return includeBoundary;  // On vertex
       }
       else if(P[i][1] == R[1] && ((P[j][0] > R[0]) == (P[i][0] < R[0])))
       {
+        isOnEdge = true;
         return includeBoundary;  // On horizontal edge
       }
     }
@@ -152,12 +158,13 @@ int winding_number(const Point<T, 2>& R,
         {
           // clang-format off
           double det = axom::numerics::determinant(P[i][0] - R[0], P[j][0] - R[0],
-                                                  P[i][1] - R[1], P[j][1] - R[1]);
+                                                   P[i][1] - R[1], P[j][1] - R[1]);
           // clang-format on
 
           // On edge
           if(axom::utilities::isNearlyEqual(det, 0.0, EPS))
           {
+            isOnEdge = true;
             return includeBoundary;
           }
 
@@ -180,6 +187,7 @@ int winding_number(const Point<T, 2>& R,
           // On edge
           if(axom::utilities::isNearlyEqual(det, 0.0, EPS))
           {
+            isOnEdge = true;
             return includeBoundary;
           }
 
@@ -197,15 +205,40 @@ int winding_number(const Point<T, 2>& R,
 }
 
 /*!
- * \brief Computes the generalized winding number for a single 2D Bezier curve
+ * \brief Computes the solid angle winding number for a 2D polygon
+ *
+ * Overload function without additional returning parameter
+ */
+template <typename T>
+int winding_number(const Point<T, 2>& R,
+                   const Polygon<T, 2>& P,
+                   bool includeBoundary = false,
+                   double EPS = 1e-8)
+{
+  bool isOnEdge = false;
+  return winding_number(R, P, isOnEdge, includeBoundary, EPS);
+}
+
+/*!
+ * \brief Computes the generalized winding number (GWN) for a single 2D Bezier curve
  *
  * \param [in] query The query point to test
  * \param [in] c The Bezier curve object 
  * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable
  * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances
  *
- * Computes the winding number using a recursive, bisection algorithm,
- * using nearly-linear Bezier curves as a base case.
+ * Computes the GWN using a recursive, bisection algorithm
+ * that constructs a polygon with the same *integer* WN as 
+ * the curve closed with a linear segment. The *generalized* WN 
+ * of the closing line is then subtracted from the integer WN to
+ * return the GWN of the original curve. (See )
+ *  
+ * Nearly-linear Bezier curves are the base case for recursion.
+ * 
+ * See Algorithm 2 in
+ *  Jacob Spainhour, David Gunderman, and Kenneth Weiss. 2024. 
+ *  Robust Containment Queries over Collections of Rational Parametric Curves via Generalized Winding Numbers. 
+ *  ACM Trans. Graph. 43, 4, Article 38 (July 2024)
  * 
  * \return double the generalized winding number.
  */
@@ -215,7 +248,64 @@ double winding_number(const Point<T, 2>& q,
                       double edge_tol = 1e-8,
                       double EPS = 1e-8)
 {
-  return detail::curve_winding_number_recursive(q, c, false, edge_tol, EPS);
+  const int ord = c.getOrder();
+  if(ord <= 0) return 0.0;
+
+  // Early return is possible for must points + curves
+  if(!c.boundingBox().contains(q))
+  {
+    return detail::linear_winding_number(q, c[0], c[ord], edge_tol);
+  }
+
+  // The first vertex of the polygon is the t=0 point of the curve
+  Polygon<T, 2> approximating_polygon(1);
+  approximating_polygon.addVertex(c[0]);
+
+  // Need to keep a running total of the GWN to account for
+  //  the winding number of coincident points
+  double gwn = 0.0;
+  bool isCoincident = false;
+  detail::construct_approximating_polygon(q,
+                                          c,
+                                          false,
+                                          edge_tol,
+                                          EPS,
+                                          approximating_polygon,
+                                          gwn,
+                                          isCoincident);
+
+  // The last vertex of the polygon is the t=1 point of the curve
+  approximating_polygon.addVertex(c[ord]);
+
+  // Compute the integer winding numberof the closed curve
+  bool isOnEdge = false;
+  double closed_curve_wn =
+    winding_number(q, approximating_polygon, isOnEdge, false, edge_tol);
+
+  // Compute the fractional value of the closed curve
+  int n = approximating_polygon.numVertices();
+  double closure_wn = detail::linear_winding_number(q,
+                                                    approximating_polygon[n - 1],
+                                                    approximating_polygon[0],
+                                                    edge_tol);
+
+  // If the point is on the boundary of the approximating polygon,
+  //  or coincident with the curve (rare), then winding_number<polygon>
+  //  doesn't return the right half-integer. Have to go edge-by-edge
+  if(isCoincident || isOnEdge)
+  {
+    closed_curve_wn = closure_wn;
+    for(int i = 1; i < n; ++i)
+    {
+      closed_curve_wn +=
+        detail::linear_winding_number(q,
+                                      approximating_polygon[i - 1],
+                                      approximating_polygon[i],
+                                      edge_tol);
+    }
+  }
+
+  return gwn + closed_curve_wn - closure_wn;
 }
 
 /*!
@@ -239,13 +329,26 @@ double winding_number(const Point<T, 2>& q,
   double ret_val = 0.0;
   for(int i = 0; i < cpoly.numEdges(); i++)
   {
-    ret_val +=
-      detail::curve_winding_number_recursive(q, cpoly[i], false, edge_tol, EPS);
+    ret_val += winding_number(q, cpoly[i], edge_tol, EPS);
   }
 
   return ret_val;
 }
 
+template <typename T>
+double winding_number(const Point<T, 2>& q,
+                      const axom::Array<BezierCurve<T, 2>>& carray,
+                      double edge_tol = 1e-8,
+                      double EPS = 1e-8)
+{
+  double ret_val = 0.0;
+  for(int i = 0; i < carray.size(); i++)
+  {
+    ret_val += winding_number(q, carray[i], false, edge_tol);
+  }
+
+  return ret_val;
+}
 //@}
 
 //@{