import bezier import numpy as np
- 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
The problem of intersecting two curves is a difficult one in computational geometry. The .Curve.intersect method (when using the ~.IntersectionStrategy.GEOMETRIC strategy) uses a combination of curve subdivision, bounding box intersection, and curve approximation (by lines) to find intersections.
intersect-1-8
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0, 0.375], ... [1.0, 0.375], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.25, 0.25], [0.75, 0.75]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.25 , 0.375], [0.75 , 0.375]])
intersect-1-9
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.5, 0.0 ], ... [0.5, 0.75], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5 , 0.6666...]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.5, 0.5]])
intersect-10-11
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [4.5, 9.0], ... [9.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0, 8.0], ... [6.0, 0.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.3333..., 0.5 ]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[3., 4.]])
intersect-8-9
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.375], ... [1.0, 0.375], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=1) >>> nodes2 = np.asfortranarray([ ... [0.5, 0.0 ], ... [0.5, 0.75], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5, 0.5]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.5 , 0.375]])
intersect-29-30
>>> nodes1 = np.asfortranarray([ ... [-1.0, 1.0], ... [ 0.5, 0.5], ... [ 0.0, 2.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [ 0.5 , 0.5 ], ... [-0.25, 1.25], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5 , 0.6666...]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0., 1.]])
For curves which intersect at exact floating point numbers, we can typically compute the intersection with zero error:
intersect-1-5
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0, 0.75], ... [0.5, -0.25], ... [1.0, 0.75], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.25, 0.25], [0.75, 0.75]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.25 , 0.375], [0.75 , 0.375]])
intersect-3-4
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [1.5, 3.0], ... [3.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [ 3.0 , 1.5 ], ... [ 2.625, -0.90625], ... [-0.75 , 2.4375 ], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.25 , 0.75 ], [0.875, 0.25 ]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.75 , 1.125 ], [2.625 , 0.65625]])
intersect-14-16
>>> nodes1 = np.asfortranarray([ ... [0.0 , 0.0 ], ... [0.375, 0.75 ], ... [0.75 , 0.375], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.25 , 0.5625], ... [0.625, 0.1875], ... [1.0 , 0.9375], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5 , 0.16666...], [0.83333..., 0.5 ]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.375 , 0.46875], [0.625 , 0.46875]])
Even for curves which don't intersect at exact floating point numbers, we can compute the intersection to machine precision:
intersect-1-2
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [1.125, 0.5], ... [0.625, -0.5], ... [0.125, 0.5], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> sq31 = np.sqrt(31.0) >>> expected_ints = np.asfortranarray([ ... [9 - sq31, 9 + sq31], ... [9 + sq31, 9 - sq31], ... ]) / 16.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) <= -53 True >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> expected_pts = np.asfortranarray([ ... [36 - 4 * sq31, 16 + sq31], ... [36 + 4 * sq31, 16 - sq31], ... ]) / 64.0 >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) -54
intersect-1-7
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0, 0.265625], ... [0.5, 0.234375], ... [1.0, 0.265625], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> sq33 = np.sqrt(33.0) >>> expected_ints = np.asfortranarray([ ... [33 - 4 * sq33, 33 - 4 * sq33], ... [33 + 4 * sq33, 33 + 4 * sq33], ... ]) / 66.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) -55 >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> expected_pts = np.asfortranarray([ ... [33 - 4 * sq33, 17], ... [33 + 4 * sq33, 17], ... ]) / 66.0 >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) -54
intersect-1-13
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0 , 0.0], ... [0.25, 2.0], ... [0.5 , -2.0], ... [0.75, 2.0], ... [1.0 , 0.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=4) >>> intersections = curve1.intersect(curve2) >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> sq7 = np.sqrt(7.0) >>> expected_ints = np.asfortranarray([ ... [7 - sq7, 7 - sq7], ... [7 + sq7, 7 + sq7], ... [ 0, 0 ], ... [ 14, 14 ], ... ]) / 14.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) <= -53 True >>> expected_pts = np.asfortranarray([ ... [7 - sq7, 6], ... [7 + sq7, 6], ... [ 0, 0], ... [ 14, 0], ... ]) / 14.0 >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) <= -53 True
intersect-21-22
>>> nodes1 = np.asfortranarray([ ... [-0.125, -0.28125], ... [ 0.5 , 1.28125], ... [ 1.125, -0.28125], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [ 1.5625, -0.0625], ... [-1.5625, 0.25 ], ... [ 1.5625, 0.5625], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> sq5 = np.sqrt(5.0) >>> expected_ints = np.asfortranarray([ ... [4 - sq5, 6 - sq5], ... [ 3, 7 ], ... [ 9, 1 ], ... [4 + sq5, 6 + sq5], ... ]) / 10.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) <= -51 True >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> expected_pts = np.asfortranarray([ ... [6 - 2 * sq5, 5 - sq5], ... [ 4, 6 ], ... [ 16, 0 ], ... [6 + 2 * sq5, 5 + sq5], ... ]) / 16.0 >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) -51
For higher degree intersections, the error starts to get a little larger.
intersect-15-25
>>> nodes1 = np.asfortranarray([ ... [0.25 , 0.625], ... [0.625, 0.25 ], ... [1.0 , 1.0 ], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0 , 0.5], ... [0.25, 1.0], ... [0.75, 1.5], ... [1.0 , 0.5], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=3) >>> intersections = curve1.intersect(curve2) >>> s_vals = np.roots([486, -3726, 13905, -18405, 6213, 1231]) >>> _, s_val, _ = np.sort(s_vals[s_vals.imag == 0].real) >>> t_vals = np.roots([4, -16, 13, 25, -28, 4]) >>> _, _, t_val = np.sort(t_vals[t_vals.imag == 0].real) >>> expected_ints = np.asfortranarray([ ... [s_val, t_val], ... ]) >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) -50 >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> x_val = (3 * s_val + 1) / 4 >>> y_val = (9 * s_val * s_val - 6 * s_val + 5) / 8 >>> expected_pts = np.asfortranarray([ ... [x_val, y_val], ... ]) >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) <= -50 True
intersect-11-26
>>> nodes1 = np.asfortranarray([ ... [0.0, 8.0], ... [6.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=1) >>> nodes2 = np.asfortranarray([ ... [0.375, 7.0], ... [2.125, 8.0], ... [3.875, 0.0], ... [5.625, 1.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=3) >>> intersections = curve1.intersect(curve2) >>> sq7 = np.sqrt(7.0) >>> expected_ints = np.asfortranarray([ ... [ 24, 24 ], ... [24 - 7 * sq7, 24 - 8 * sq7], ... [24 + 7 * sq7, 24 + 8 * sq7], ... ]) / 48.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) -53 >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> expected_pts = np.asfortranarray([ ... [ 72, 96 ], ... [72 - 21 * sq7, 96 + 28 * sq7], ... [72 + 21 * sq7, 96 - 28 * sq7], ... ]) / 24.0 >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) -50
intersect-8-27
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.375], ... [1.0, 0.375], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=1) >>> nodes2 = np.asfortranarray([ ... [0.125, 0.25 ], ... [0.375, 0.75 ], ... [0.625, 0.0 ], ... [0.875, 0.1875], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=3) >>> intersections = curve1.intersect(curve2) >>> points = curve1.evaluate_multi(intersections[:, 0]) >>> s_val2, s_val1, _ = np.sort(np.roots( ... [17920, -29760, 13512, -1691])) >>> t_val2, t_val1, _ = np.sort(np.roots([35, -60, 24, -2])) >>> expected_ints = np.asfortranarray([ ... [s_val1, t_val1], ... [s_val2, t_val2], ... ]) >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) -51 >>> expected_pts = np.asfortranarray([ ... [s_val1, 0.375], ... [s_val2, 0.375], ... ]) >>> max_err = np.max(np.abs(points - expected_pts)) >>> binary_exponent(max_err) -51
intersect-1-18
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [1.0, 0.0], ... [1.5, -1.0], ... [2.0, 0.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[1., 0.]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[1., 0.]])
intersect-1-19
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [2.0, 0.0], ... [1.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[1., 1.]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[1., 0.]])
intersect-10-17
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [4.5, 9.0], ... [9.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [11.0, 8.0], ... [ 7.0, 10.0], ... [ 3.0, 4.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.3333..., 1. ]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[3., 4.]])
intersect-12-self
>>> nodes = np.asfortranarray([ ... [ 0.0 , 2.0 ], ... [-1.0 , 0.0 ], ... [ 1.0 , 1.0 ], ... [-0.75, 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, 0 ], ... [3 - sq5, sq5], ... ]) / 3.0 >>> max_err = np.max(np.abs(intersections - expected_ints)) >>> binary_exponent(max_err) -54 >>> left.evaluate_multi(intersections[:, 0]) array([[-0.09375 , 0.828125], [-0.25 , 1.375 ]])
Intersections that occur at points of tangency are in general problematic. For example, consider
The first curve is the zero set of y − 2x(1 − x), so plugging in the second curve gives
0 = t2 + (1 − t)2 − 2t(1 − t) = (2t − 1)2.
This shows that a point of tangency is equivalent to a repeated root of a polynomial. For this example, the intersection process successfully terminates
intersect-1-6
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.0, 1.0], ... [0.5, 0.0], ... [1.0, 1.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5, 0.5]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[0.5, 0.5]])
However this library mostly avoids (for now) computing tangent intersections. For example, the curves
have a tangent intersection that this library fails to compute:
intersect-14-15
>>> nodes1 = np.asfortranarray([ ... [0.0 , 0.0 ], ... [0.375, 0.75 ], ... [0.75 , 0.375], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.25 , 0.625], ... [0.625, 0.25 ], ... [1.0 , 1.0 ], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> curve1.intersect(curve2) Traceback (most recent call last): ... NotImplementedError: Line segments parallel.
This failure comes from the fact that the linear approximations of the curves near the point of intersection are parallel.
As above, we can find some cases where tangent intersections are resolved:
intersect-10-23
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [4.5, 9.0], ... [9.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [3.0, 4.5], ... [8.0, 4.5], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=1) >>> intersections = curve1.intersect(curve2) >>> intersections array([[0.5, 0.3]]) >>> curve1.evaluate_multi(intersections[:, 0]) array([[4.5, 4.5]])
but even by rotating an intersection (from above) that we know works
we still see a failure
intersect-28-29
>>> nodes1 = np.asfortranarray([ ... [ 0.0, 0.0], ... [-0.5, 1.5], ... [ 1.0, 1.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [-1.0, 1.0], ... [ 0.5, 0.5], ... [ 0.0, 2.0], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> curve1.intersect(curve2) Traceback (most recent call last): ... NotImplementedError: The number of candidate intersections is too high.
In addition to points of tangency, coincident curve segments are (for now) not supported. For the curves
the library fails as well
intersect-1-24
>>> nodes1 = np.asfortranarray([ ... [0.0, 0.0], ... [0.5, 1.0], ... [1.0, 0.0], ... ]) >>> curve1 = bezier.Curve(nodes1, degree=2) >>> nodes2 = np.asfortranarray([ ... [0.25, 0.375], ... [0.75, 0.875], ... [1.25, -0.625], ... ]) >>> curve2 = bezier.Curve(nodes2, degree=2) >>> curve1.intersect(curve2) Traceback (most recent call last): ... NotImplementedError: The number of candidate intersections is too high.