Implement "real" self-intersection checks for curves

[dhermes]

See: https://doi.org/10.1016/0166-3615(89)90072-9 / https://drive.google.com/file/d/1cjSSmuFfJZgujlTMaerPJzIu_5i5NERt/view

I started looking into this partly in a gist exploring angle of intersection for two planar vectors.

[dhermes]

import scipy.integrate

import bezier


def hodograph(curve):
    nodes = curve.nodes
    dp = curve.degree * (nodes[:, 1:] - nodes[:, :-1])
    return bezier.Curve(dp, degree=curve.degree - 1, _copy=False)


def turning_angle(curve):
    """Compute the exact turning angle via integration.

    This is done by considering how the angle of B'(s) = [x'(s), y'(s)] changes
    in small intervals in the parameter space; where the angle is
    a(s) = arctan(y'(s) / x'(s)). Computing a(s + ds) - a(s) as ds --> 0 leaves
    us with the **signed** angle change
    dtheta = (y'' x' - x'' y') / (x'**2 + y'**2).

    The turning angle is then INT_{s = 0}^1 |dtheta| ds.
    """
    dp = hodograph(curve)
    dpp = hodograph(dp)

    def d_theta(s):
        (xp, yp), = dp.evaluate(s).T
        (xpp, ypp), = dpp.evaluate(s).T
        return abs(ypp * xp - xpp * yp) / (xp * xp + yp * yp)

    computed, _ = scipy.integrate.quad(d_theta, 0.0, 1.0)
    return computed

and to compute an angle for 3 consecutive nodes vA, vB, vC

import numpy as np


def to_polar(vec):
    r = np.linalg.norm(vec, ord=2)
    x, y = vec
    theta = np.arctan2(y, x)
    return r, theta


def discrete_angle(xA, yA, xB, yB, xC, yC):
    # Compute actual vectors and convert to polar coordinates relative to
    # (xB, yB)
    vec0 = np.array([xB - xA, yB - yA])
    _, theta0 = to_polar(vec0)
    vec1 = np.array([xC - xB, yC - yB])
    _, theta1 = to_polar(vec1)
    # Make sure |theta1 - theta0| <= pi
    delta = theta1 - theta0
    if delta > np.pi:
        theta1 -= 2 * np.pi
    elif delta < -np.pi:
        theta1 += 2 * np.pi

    return theta1 - theta0

which for an entire curve amounts to

def abs_angle_sum(curve):
    result = 0
    for j in range(curve.degree - 1):
        xA, xB, xC = curve._nodes[0, j : j + 3]
        yA, yB, yC = curve._nodes[1, j : j + 3]
        d_theta = discrete_angle(xA, yA, xB, yB, xC, yC)
        result += abs(d_theta)

    return result

[dhermes]

Criterion: if abs_angle_sum(curve) < np.pi (to some threshold), impossible

Algorithm:

1. Check angle sum for curve C0_0
2. If angle sum >= np.pi, split in half C1_0 and C1_1
3. Now, look for regular intersections of C1_0 and C1_1 (excluding the shared subdivision point) and recursively apply self-intersection check to both C1_0 and C1_1