Adding Curve.self_intersections().

docs/algorithms/curve-curve-intersection.rst

@@ -567,35 +567,6 @@ Intersections at Endpoints
 .. image:: ../images/curves10_and_17.png
    :align: center
 
-Detecting Self-Intersections
-----------------------------
-
-.. doctest:: intersect-12-self
-   :options: +NORMALIZE_WHITESPACE
-
-   >>> nodes = np.asfortranarray([
-   ...     [0.0, -1.0, 1.0, -0.75 ],
-   ...     [2.0,  0.0, 1.0,  1.625],
-   ... ])
-   >>> curve = bezier.Curve(nodes, degree=3)
-   >>> left, right = curve.subdivide()
-   >>> intersections = left.intersect(right)
-   >>> sq5 = np.sqrt(5.0)
-   >>> expected_ints = np.asfortranarray([
-   ...     [3, 3 - sq5],
-   ...     [0, sq5    ],
-   ... ]) / 3.0
-   >>> max_err = np.max(np.abs(intersections - expected_ints))
-   >>> binary_exponent(max_err)
-   -53
-   >>> s_vals = np.asfortranarray(intersections[0, :])
-   >>> left.evaluate_multi(s_vals)
-   array([[-0.09375 , -0.25  ],
-          [ 0.828125,  1.375 ]])
-
-.. image:: ../images/curves42_and_43.png
-   :align: center
-
 Limitations
 -----------
 

docs/make_images.py

@@ -355,6 +355,54 @@ def curve_intersect(curve1, curve2, s_vals):
     save_image(ax.figure, "curve_intersect.png")
 
 
+def curve_self_intersect2(curve, self_intersections):
+    """Image for :meth`.Curve.self_intersections` docstring."""
+    if NO_IMAGES:
+        return
+
+    ax = curve.plot(256, color=BLUE)
+    if self_intersections.shape != (2, 1):
+        raise ValueError("Unexpected shape", self_intersections)
+
+    s1_val = self_intersections[0, 0]
+    intersection_xy = curve.evaluate(s1_val)
+    ax.plot(
+        intersection_xy[0, :],
+        intersection_xy[1, :],
+        color="black",
+        linestyle="None",
+        marker="o",
+    )
+    ax.axis("scaled")
+    ax.set_xlim(-0.8125, 0.0625)
+    ax.set_ylim(0.75, 2.125)
+    save_image(ax.figure, "curve_self_intersect2.png")
+
+
+def curve_self_intersect3(curve, self_intersections):
+    """Image for :meth`.Curve.self_intersections` docstring."""
+    if NO_IMAGES:
+        return
+
+    ax = curve.plot(256, color=BLUE)
+    if self_intersections.shape != (2, 2):
+        raise ValueError("Unexpected shape", self_intersections)
+
+    s1_vals = np.asfortranarray(self_intersections[0, :])
+    intersection_xy = curve.evaluate_multi(s1_vals)
+    ax.plot(
+        intersection_xy[0, :],
+        intersection_xy[1, :],
+        color="black",
+        linestyle="None",
+        marker="o",
+    )
+    ax.axis("scaled")
+    ax.set_xlim(-330.0, 330.0)
+    ax.set_ylim(0.125, 266.0)
+    save_image(ax.figure, "curve_self_intersect3.png")
+
+
 def triangle_constructor(triangle):
     """Image for :class`.Triangle` docstring."""
     if NO_IMAGES:

docs/releases/latest.rst

@@ -52,6 +52,13 @@ Additive Changes
    (`#251 <https://github.com/dhermes/bezier/pull/251>`__). For example, in
    ``Curve.intersect()``
    (`doc <https://bezier.readthedocs.io/en/2021.2.13.dev1/python/reference/bezier.curve.html#bezier.curve.Curve.intersect>`__)
+-  Adding provisional support for self-intersection checks in planar curves
+   (`#265 <https://github.com/dhermes/bezier/pull/265>`__,
+   `#267 <https://github.com/dhermes/bezier/pull/267>`__). Fixed
+   `#165 <https://github.com/dhermes/bezier/issues/165>`__ and
+   `#171 <https://github.com/dhermes/bezier/issues/171>`__.
+   Supported via ``Curve.self_intersections()``
+   `method <https://bezier.readthedocs.io/en/2021.2.13.dev1/python/reference/bezier.curve.html#bezier.curve.Curve.self_intersections>`__.
 
 Documentation
 -------------

src/python/bezier/curve.py

@@ -22,6 +22,13 @@
 
    import numpy as np
    import bezier
+
+   def binary_exponent(value):
+       if value == 0.0:
+           return -np.inf
+       _, result = np.frexp(value)
+       # Shift [1/2, 1) --> [1, 2) borrows one from exponent
+       return result - 1
 """
 
 import numpy as np
@@ -32,6 +39,7 @@
 from bezier import _plot_helpers
 from bezier import _symbolic
 from bezier.hazmat import algebraic_intersection
+from bezier.hazmat import geometric_intersection
 from bezier.hazmat import intersection_helpers
 
 
@@ -457,6 +465,123 @@ def intersect(
         st_vals, _ = all_intersections(self._nodes, other._nodes)
         return st_vals
 
+    def self_intersections(
+        self, strategy=IntersectionStrategy.GEOMETRIC, verify=True
+    ):
+        """Find the points where the curve intersects itself.
+
+        For curves in general position, there will be no self-intersections:
+
+        .. doctest:: curve-self-intersect1
+           :options: +NORMALIZE_WHITESPACE
+
+           >>> nodes = np.asfortranarray([
+           ...     [0.0, 1.0, 0.0],
+           ...     [0.0, 1.0, 2.0],
+           ... ])
+           >>> curve = bezier.Curve(nodes, degree=2)
+           >>> curve.self_intersections()
+           array([], shape=(2, 0), dtype=float64)
+
+        However, some curves do have self-intersections. Consider a cubic
+        with
+
+        .. math::
+
+           B\\left(\\frac{3 - \\sqrt{5}}{6}\\right) =
+               B\\left(\\frac{3 + \\sqrt{5}}{6}\\right)
+
+        .. image:: ../../images/curve_self_intersect2.png
+           :align: center
+
+        .. doctest:: curve-self-intersect2
+           :options: +NORMALIZE_WHITESPACE
+
+           >>> nodes = np.asfortranarray([
+           ...     [0.0, -1.0, 1.0, -0.75 ],
+           ...     [2.0,  0.0, 1.0,  1.625],
+           ... ])
+           >>> curve = bezier.Curve(nodes, degree=3)
+           >>> self_intersections = curve.self_intersections()
+           >>> sq5 = np.sqrt(5.0)
+           >>> expected = np.asfortranarray([
+           ...     [3 - sq5],
+           ...     [3 + sq5],
+           ... ]) / 6.0
+           >>> max_err = np.max(np.abs(self_intersections - expected))
+           >>> binary_exponent(max_err)
+           -53
+
+        .. testcleanup:: curve-self-intersect2
+
+           import make_images
+           make_images.curve_self_intersect2(curve, self_intersections)
+
+        Some (somewhat pathological) curves can have multiple
+        self-intersections, though the number possible is largely constrained
+        by the degree. For example, this degree six curve has two
+        self-intersections:
+
+        .. image:: ../../images/curve_self_intersect3.png
+           :align: center
+
+        .. doctest:: curve-self-intersect3
+           :options: +NORMALIZE_WHITESPACE
+
+           >>> nodes = np.asfortranarray([
+           ...     [-300.0, 227.5 ,  -730.0,    0.0 ,   730.0, -227.5 , 300.0],
+           ...     [ 150.0, 953.75, -2848.0, 4404.75, -2848.0,  953.75, 150.0],
+           ... ])
+           >>> curve = bezier.Curve(nodes, degree=6)
+           >>> self_intersections = curve.self_intersections()
+           >>> 6.0 * self_intersections
+           array([[1., 4.],
+                  [2., 5.]])
+           >>> curve.evaluate_multi(self_intersections[:, 0])
+           array([[-150., -150.],
+                  [  75.,   75.]])
+           >>> curve.evaluate_multi(self_intersections[:, 1])
+           array([[150., 150.],
+                  [ 75.,  75.]])
+
+        .. testcleanup:: curve-self-intersect3
+
+           import make_images
+           make_images.curve_self_intersect3(curve, self_intersections)
+
+        Args:
+            strategy (Optional[ \
+                ~bezier.hazmat.intersection_helpers.IntersectionStrategy]): The
+                intersection algorithm to use. Defaults to geometric.
+            verify (Optional[bool]): Indicates if extra caution should be
+                used to verify assumptions about the current curve. Can be
+                disabled to speed up execution time. Defaults to :data:`True`.
+
+        Returns:
+            numpy.ndarray: ``2 x N`` array of ``s1``- and ``s2``-parameters
+            where self-intersections occur (possibly empty). For each pair
+            we have :math:`s_1 \\neq s_2` and :math:`B(s_1) = B(s_2)`.
+
+        Raises:
+            NotImplementedError: If at the curve isn't two-dimensional
+                (and ``verify=True``).
+            NotImplementedError: If ``strategy`` is not
+                :attr:`~.IntersectionStrategy.GEOMETRIC`.
+        """
+        if strategy != IntersectionStrategy.GEOMETRIC:
+            raise NotImplementedError(
+                "Only geometric strategy for self-intersection detection"
+            )
+        if verify:
+            if self._dimension != 2:
+                raise NotImplementedError(
+                    "Self-intersection only implemented in 2D",
+                    "Current dimension",
+                    self._dimension,
+                )
+
+        return geometric_intersection.self_intersections(self._nodes)
+
     def elevate(self):
         r"""Return a degree-elevated version of the current curve.
 

src/python/bezier/hazmat/curve_helpers.py

@@ -1056,6 +1056,10 @@ def full_reduce(nodes):
 def discrete_turning_angle(nodes):
     r"""Determine the absolute sum of B |eacute| zier node angles.
 
+    .. note::
+
+       This assumes, but does not check, that ``nodes`` is ``2 x N``.
+
     For the set of vectors :math:`v_j` given by ``nodes``, the discrete
     angles :math:`\theta_j` at each internal node is given by
 
@@ -1087,19 +1091,8 @@ def discrete_turning_angle(nodes):
 
     Returns:
         float: The (discrete) turning angle.
-
-    Raises:
-        NotImplementedError: If the curve's dimension is not ``2``.
     """
-    dimension, num_nodes = nodes.shape
-
-    if dimension != 2:
-        raise NotImplementedError(
-            "2D is the only supported dimension",
-            "Current dimension",
-            dimension,
-        )
-
+    _, num_nodes = nodes.shape
     if num_nodes < 3:
         return 0.0
 

tests/unit/hazmat/test_curve_helpers.py

@@ -1069,24 +1069,6 @@ def _call_function_under_test(nodes):
 
         return curve_helpers.discrete_turning_angle(nodes)
 
-    def test_invalid_dimension(self):
-        nodes = np.asfortranarray(
-            [
-                [0.0, 1.0, 3.0],
-                [0.0, 4.0, 5.0],
-                [0.0, 2.0, 1.0],
-            ]
-        )
-        with self.assertRaises(NotImplementedError) as exc_info:
-            self._call_function_under_test(nodes)
-
-        expected = (
-            "2D is the only supported dimension",
-            "Current dimension",
-            3,
-        )
-        self.assertEqual(expected, exc_info.exception.args)
-
     def test_linear(self):
         nodes = np.asfortranarray(
             [

tests/unit/test_curve.py

@@ -23,6 +23,9 @@
 from tests.unit import utils
 
 
+SPACING = np.spacing  # pylint: disable=no-member
+
+
 class TestCurve(utils.NumPyTestCase):
 
     ZEROS = np.zeros((2, 2), order="F")
@@ -274,6 +277,58 @@ def test_intersect_unsupported_dimension(self):
         with self.assertRaises(NotImplementedError):
             curve2.intersect(curve1)
 
+    def test_self_intersections_unsupported_dimension(self):
+        nodes = np.asfortranarray([[0.0, 1.0, 2.0]])
+        curve = self._make_one(nodes, 2)
+
+        with self.assertRaises(NotImplementedError) as exc_info:
+            curve.self_intersections()
+
+        expected = (
+            "Self-intersection only implemented in 2D",
+            "Current dimension",
+            1,
+        )
+        self.assertEqual(expected, exc_info.exception.args)
+
+    def test_self_intersections_unsupported_strategy(self):
+        from bezier.hazmat import intersection_helpers
+
+        nodes = np.asfortranarray([[0.0, 1.0, 2.0], [0.0, 1.0, 2.0]])
+        curve = self._make_one(nodes, 2)
+
+        strategy = intersection_helpers.IntersectionStrategy.ALGEBRAIC
+        with self.assertRaises(NotImplementedError) as exc_info:
+            curve.self_intersections(strategy=strategy)
+
+        expected = ("Only geometric strategy for self-intersection detection",)
+        self.assertEqual(expected, exc_info.exception.args)
+
+    def test_self_intersections_valid(self):
+        nodes = np.asfortranarray(
+            [
+                [0.0, -1.0, 1.0, -0.75],
+                [2.0, 0.0, 1.0, 1.625],
+            ]
+        )
+        curve = self._make_one(nodes, 3)
+        intersections = curve.self_intersections(verify=False)
+        self.assertEqual((2, 1), intersections.shape)
+        intersections_verified = curve.self_intersections(verify=True)
+        self.assertEqual(intersections, intersections_verified)
+
+        sq5 = np.sqrt(5.0)
+        expected_s = 0.5 - sq5 / 6.0
+        local_eps = 3 * abs(SPACING(expected_s))
+        self.assertAlmostEqual(
+            expected_s, intersections[0, 0], delta=local_eps
+        )
+        expected_t = 0.5 + sq5 / 6.0
+        local_eps = abs(SPACING(expected_t))
+        self.assertAlmostEqual(
+            expected_t, intersections[1, 0], delta=local_eps
+        )
+
     def test_elevate(self):
         nodes = np.asfortranarray([[0.0, 1.0, 3.0, 3.5], [0.5, 1.0, 2.0, 4.0]])
         curve = self._make_one(nodes, 3)