Update 2D GWN algorithm to match paper

[jcs15c]

Summary

This PR is an upgrade of the existing 2D generalized winding number (GWN) algorithm so that the implementation more closely matches the pseudocode in the paper "Robust Containment Queries over Collections of Rational Parametric Curves via Generalized Winding Numbers". The updated algorithm is more analytically tractable and computationally efficient, significantly reducing the number of trigonometric operations needed to compute the GWN for near-field points.

This upgrade does not change the external API for existing GWN methods, but does alter function calls in the detail namespace, i.e., replacing curve_winding_number_recursive with construct_approximating_polygon.

Additionally, two additional functions are added to the interface provided by winding_number.hpp:

• winding_number<axom::Array<BezierCurve>>
• A version of winding_number<Polygon> that has an additional reference parameter to report points located on the polygon boundary.

Algorithm Details

Specifically, the current implementation recursively discretizes a Bezier curve into (potentially many) linear segments and computes the GWN for each segment, evaluating an arccosine for each. In contrast, the updated algorithm only evaluates a single GWN value for each point. In this way, the new algorithm more directly decomposes GWN evaluation of Bezier curves into the following sequence of actions:

1. Check the point against an axis-aligned bounding-box for the curve. If outside, we can treat curve as a line between its endpoints, and exit the algorithm.
2. Close the curve with a straight line connecting the endpoints, and compute the generalized winding number of the line at the point.
3. Adaptively construct an axom::Polygon object that provably has the same integer winding number as the closed curve at the point, and evaluate the integer with a point-in-polygon algorithm.
4. Subtract the generalized winding number from the integer winding number from the resulting integer the generalized winding number of the closing line.

Extra bugfix via refactor of winding_number<polygon>

During review, an additional bug was found relating to the evaluation of the winding number for points that are near to the closing line. Specifically, it was possible for a point to be considered as interior/exterior to the approximating polygon object, but coincident with the closing curve during the computation of it's generalized winding number, causing an incorrect value to be subtracted from the GWN.

This "off-by-0.5" error, unacceptably high for points near the closing curve, is due to a different physical interpretation of numerical tolerances between winding_number<polygon> and winding_number<line>, where the former checks for colinearity of an edge's endpoints and the point using a determinant formula, and the latter checks for physical distance between the edge and the point.

To fix this, the various numerical tolerances in winding_number<polygon> (which also includes an Linf distances to vertices) have been replaced with checks for physical distances between the query point and edges. This approach is nominally more costly, but makes the method considerably more usable in conjunction with other winding_number methods. Beyond changing variable names, this change does not alter the interface of the existing method.

An additional test has been added to primal_winding_number to verify the corrected behavior.

[kennyweiss]

Thanks @jcs15c !

Before merging, I'd like to check that our winding number example application still works on all the test shapes from our data repo -- https://github.com/LLNL/axom_data/tree/main/contours/svg

(I haven't gotten around to adding automated tests for these yet, but it's on my todo list).

[rhornung67]

Thank you @jcs15c Very nice summary of the changes in the PR.

[agcapps]

Thank you, Jacob!

[jcs15c]

Before merging, I'd like to check that our winding number example application still works on all the test shapes from our data repo -- https://github.com/LLNL/axom_data/tree/main/contours/svg

@kennyweiss I checked against these shapes, and although the results are consistent with the develop branch in terms of performance and (mostly) classification accuracy, I did find through this check a handful of discrepancies. This led to the discovery of a much larger class of errors owing to some inconsistent numerical tolerances used by different winding number methods. I've fixed the issue by modifying the tolerances used by winding_number<polygon>, which I've described in more detail in the main PR text.

Either way, it'd definitely be worth it to automate those tests somehow, these issues would have stuck out immediately.

[kennyweiss]

Awesome -- thanks @jcs15c. Now that the issue's resolved, feel free to merge this PR. We'll automate those tests in a follow-up PR.

[kennyweiss]

This is a great test!
Would it makes sense to also test Point2D({0, -diff}) ?

[jcs15c]

Probably! The internal tolerance checks don't care which side of the closure the point is on, but might as well check it here too.

[jcs15c]

It actually made a lot of sense! There was the exact same issue with the BoundingBox check, since the closure is coincident with an edge. This means it would be counted as "outside" the bounding box but "on" the edge, and therefore give the wrong answer. I'm actually surprised that this didn't come up as an issue in the previous versions of winding_number<Bezier>, since nothing new to this PR has changed this behavior. Neat!

This one's a lot easier to fix, since it just involves adding the same edge_tol parameter to the bounding box checks via their .expand() method.