NURBS curve intersections in 2D

[dylan-copeland]

Summary

• This PR adds a small function to find intersections of two NURBS curves in 2D, assuming finitely many intersections.
• A unit test is included for two simple curves with a single intersection.

[kennyweiss]

Thanks for the contribution @dylan-copeland !

Please add some additional doxygen/code comments and some additional tests per my PR comments.

[kennyweiss]

Presumably, the entries in p1 and p2 correspond to the same point in space. Please mention this in the doxygen.

[dylan-copeland]

Done in 9ddf55b.

[jcs15c]

It's probably also worth mentioning in the doxygen that this method also parameterizes the curves as half-open (since the underlying bezier-bezier intersection routine does too), such that intersections at either terminal endpoint will be ignored.

[dylan-copeland]

I added a comment about endpoints in abdb66d.

[kennyweiss]

Consider adding a comment that it is inside the range of knots[i] and knots[i+1] because the know array was created with getUniqueKnots()

[dylan-copeland]

Done in 9ddf55b.

[kennyweiss]

Please add some more tests, including:

• curves that don't intersect
• curves that intersect more than once
• curves that intersect more than once at the same physical points
• curves that intersect at a knot vector

[dylan-copeland]

I added tests for the first two suggestions in 9ddf55b. For the third, I am not sure of a good way to make curves in axom intersecting more than once at the same points. There is a way to make a circle up to 2*pi. Making something more interesting would need to be done manually (to find control points, weights, etc), unless someone has suggestions. For the fourth, I made a test that seems to expose issues with intersect_bezier_curves. The test was two circles of radius 1, centered at (0,0) and (2,0), which should touch at (1,0), which should be a knot for both circles. The circles need not be complete (can be arcs). For several variations of this, it failed to find any intersection. Unless I am missing something, intersect_bezier_curves needs to be made more robust. It seems to recursively split curves in half to check each half, using a linear approximation when the curve is within the prescribed tolerance of being linear. This misses some intersections, e.g. two circles that touch at one point.

[jcs15c]

I don't know if we expect intersect_bezier_curves to find intersections at points of tangencies. It says as much in the doxygen for intersect( bezier, bezier ), but not for intersect_bezier_curves, so I think that's an absence in the documentation more than anything. That's not to say that there's philosophically wrong with tangent intersections, it's just a considerably hard edge case to handle, and we've been kicking it down the road for a while now.

[dylan-copeland]

I added more comments in abdb66d. I also tried to create a test with intersection of two circular arcs of 180 degrees, which have knots 0, 0.5, 1. The arcs intersect at their midpoints, and the intersection is found correctly, but it is not at the expected parameter of 0.5, matching a knot. Evaluating the arcs at points in [0,1] shows the parameterization is unusual, and the location of knot 0.5 is some point far from the midpoint. So it wasn't easy to make a test intersecting at a knot, but it should be possible for a manually defined curve. I'm not sure it is worth the effort to make these tests.

[kennyweiss]

I was thinking of something like this:

NURBSCurve2D make_cubic_shape()
{
  // Open cubic NURBS curve (degree 3), non-rational.
  // 24 control points, 28 knots (n + p + 2: 24 + 3 + 1 = 28).
  return NURBSCurve2D(
    axom::Array<Point2D> {
      Point2D {-5.0000000000, +3.6000000000}, Point2D {-3.4000000000, +3.0000000000},
      Point2D {-1.8000000000, +2.2000000000}, Point2D {-0.4000000000, +1.0000000000},
      Point2D {+1.0000000000, -0.4000000000}, Point2D {+2.4000000000, -1.6000000000},
      Point2D {+3.6000000000, -2.4000000000}, Point2D {+4.4000000000, -1.0000000000},
      Point2D {+3.6000000000, +0.8000000000}, Point2D {+2.0000000000, +2.0000000000},
      Point2D {+0.0000000000, +2.8000000000}, Point2D {-2.0000000000, +2.0000000000},
      Point2D {-3.2000000000, +0.8000000000}, Point2D {-2.6000000000, -0.8000000000},
      Point2D {-1.0000000000, -1.6000000000}, Point2D {+0.4000000000, -2.2000000000},
      Point2D {+1.6000000000, -2.6000000000}, Point2D {+0.0000000000, -3.4000000000},
      Point2D {-2.2000000000, -3.0000000000}, Point2D {-2.8000000000, -1.4000000000},
      Point2D {-1.0000000000, +0.4000000000}, Point2D {+1.8000000000, +1.8000000000},
      Point2D {+3.4000000000, +3.0000000000}, Point2D {+5.0000000000, +3.6000000000}
    },
    axom::Array<double> {
      0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,
      0.0476190476, 0.0952380952, 0.1428571429, 0.1904761905,
      0.2380952381, 0.2857142857, 0.3333333333, 0.3809523810,
      0.4285714286, 0.4761904762, 0.5238095238, 0.5714285714,
      0.6190476190, 0.6666666667, 0.7142857143, 0.7619047619,
      0.8095238095, 0.8571428571, 0.9047619048, 0.9523809524,
      1.0000000000, 1.0000000000, 1.0000000000, 1.0000000000
    },
    3 /* degree */);
}

NURBSCurve2D make_ellipse_curve()
{
  // Quadratic rational NURBS ellipse (degree 2), 9-control-point 4-arc construction. 
  // Center, semi-axes:
  //   center      = (+0.055509, +0.246636)
  //   semi-axes   = (2.506091, 2.506234)   (a == b -> circle in this case)
  //   rotation    = +0.0000 degrees
  // Corner control-point weights are w = cos(pi/4) = sqrt(2)/2.
  return NURBSCurve2D(
    axom::Array<Point2D> {
      Point2D {+2.5615994286, +0.2466357797}, Point2D {+2.5615994286, +2.7528699841}, Point2D {+0.0555086484, +2.7528699841},
      Point2D {-2.4505821318, +2.7528699841}, Point2D {-2.4505821318, +0.2466357797}, Point2D {-2.4505821318, -2.2595984246},
      Point2D {+0.0555086484, -2.2595984246}, Point2D {+2.5615994286, -2.2595984246}, Point2D {+2.5615994286, +0.2466357797}
    },
    axom::Array<double> {
      1.0000000000, 0.7071067812, 1.0000000000,
      0.7071067812, 1.0000000000, 0.7071067812,
      1.0000000000, 0.7071067812, 1.0000000000
    },
    axom::Array<double> {
      0.0000000000, 0.0000000000, 0.0000000000, 0.2500000000,
      0.2500000000, 0.5000000000, 0.5000000000, 0.7500000000,
      0.7500000000, 1.0000000000, 1.0000000000, 1.0000000000
    },
    2 /* degree */);
}

Per a python script generated by claude, the intersections should be:

  "approximate_curve_curve_intersections": [
    [
      -1.7060766733448747,
      2.0292202966044366
    ],
    [
      2.044376271360025,
      -1.2782097206437852
    ],
    [
      1.884074237582652,
      1.9604430112255964
    ],
    [
      -2.1644784828526533,
      -0.9161966686461217
    ],
    [
      0.5039968961448462,
      -2.2191102081421588
    ]
  ]

[dylan-copeland]

Thanks @kennyweiss for this nice example. I added it in a812569.

[kennyweiss]

Do we need to handle the case where the intersection is at a knot vector, and could be listed multiple times in the output? (e.g. by deduplicating the list). We should mention how this is handled in the doxygen.

[jcs15c]

The intersect(bezier, bezier) routine assumes a half-open parameter space (since bezier-bezier intersections use recursive subdivision to get all the intersections), so there shouldn't need to be a deduplicating. If an intersection is listed twice, that's almost always going to suggest some numerical issue.

[kennyweiss]

Thanks for contribution and the updates @dylan-copeland !

[jcs15c]

Looks great, thanks for clarifying in the documentation about tangencies!