bezier

bezier package

Helper for Bézier Curves, Triangles, and Higher Order Objects.

Intended to perform basic operations on Bézier objects such as intersections, length/area/etc. computations, subdivision, implicitization and other relevant information.

Plotting utilities will also be provided.

Submodules

bezier.curve module

Helper for Bézier Curves.

See Curve-Curve Intersection for examples using the Curve class to find intersections.

bezier._intersection_helpers.linearization_error(curve)

Compute the maximum error of a linear approximation.

We use the line

\[L(s) = v_0 (1 - s) + v_n s\]

and compute a bound on the maximum error

\[\max_{s \in \left[0, 1\right]} \|B(s) - L(s)\|_2.\]

Rather than computing the actual maximum (a tight bound), we use an upper bound via the remainder from Lagrange interpolation in each component. This leaves us with \(\frac{s(s - 1)}{2!}\) times the second derivative in each component.

The second derivative curve is degree \(d = n - 2\) and is given by

\[B''(s) = n(n - 1) \sum_{j = 0}^{d} \binom{d}{j} s^j (1 - s)^{d - j} \cdot \Delta^2 v_j\]

Due to this form (and the convex combination property of Bézier Curves) we know each component of the second derivative will be bounded by the maximum of that component among the \(\Delta^2 v_j\).

For example, the curve

\[\begin{split}B(s) = \left[\begin{array}{c} 0 \\ 0 \end{array}\right] (1 - s)^2 + \left[\begin{array}{c} 1 \\ 3 \end{array}\right] 2s(1 - s) + \left[\begin{array}{c} -2 \\ 9 \end{array}\right] s^2\end{split}\]

has \(B''(s) \equiv \left[\begin{array}{c} -8 \\ 6 \end{array}\right]\) which has norm \(10\) everywhere, hence the maximum error is

\[\left.\frac{s(1 - s)}{2!} \cdot 10\right|_{s = \frac{1}{2}} = \frac{5}{4}.\]
>>> nodes = np.array([
...     [ 0.0, 0.0],
...     [ 1.0, 3.0],
...     [-2.0, 9.0],
... ])
>>> curve = bezier.Curve(nodes)
>>> linearization_error(curve)
1.25
Parameters:curve (Curve) – A curve to be approximated by a line.
Returns:The maximum error between the curve and the linear approximation.
Return type:float
bezier._intersection_helpers.newton_refine(s, curve1, t, curve2)

Apply one step of 2D Newton’s method.

We want to use Newton’s method on the function

\[F(s, t) = B_1(s) - B_2(t)\]

to refine \(\left(s_{\ast}, t_{\ast}\right)\). Using this, and the Jacobian \(DF\), we “solve”

\[\begin{split}\left[\begin{array}{c} 0 \\ 0 \end{array}\right] \approx F\left(s_{\ast} + \Delta s, t_{\ast} + \Delta t\right) \approx F\left(s_{\ast}, t_{\ast}\right) + \left[\begin{array}{c c} B_1'\left(s_{\ast}\right) & - B_2'\left(t_{\ast}\right) \end{array}\right] \left[\begin{array}{c} \Delta s \\ \Delta t \end{array}\right]\end{split}\]

and refine with the component updates \(\Delta s\) and \(\Delta t\).

Note

This implementation assumes curve1 and curve2 live in \(\mathbf{R}^2\).

For example, the curves

\[\begin{split}\begin{align*} B_1(s) &= \left[\begin{array}{c} 0 \\ 0 \end{array}\right] (1 - s)^2 + \left[\begin{array}{c} 2 \\ 4 \end{array}\right] 2s(1 - s) + \left[\begin{array}{c} 4 \\ 0 \end{array}\right] s^2 \\ B_2(t) &= \left[\begin{array}{c} 2 \\ 0 \end{array}\right] (1 - t) + \left[\begin{array}{c} 0 \\ 3 \end{array}\right] t \end{align*}\end{split}\]

intersect at the point \(B_1\left(\frac{1}{4}\right) = B_2\left(\frac{1}{2}\right) = \frac{1}{2} \left[\begin{array}{c} 2 \\ 3 \end{array}\right]\). However, starting from the wrong point we have

\[\begin{split}\begin{align*} F\left(\frac{3}{8}, \frac{1}{4}\right) &= \frac{1}{8} \left[\begin{array}{c} 0 \\ 9 \end{array}\right] \\ DF\left(\frac{3}{8}, \frac{1}{4}\right) &= \left[\begin{array}{c c} 4 & 2 \\ 2 & -3 \end{array}\right] \\ \Longrightarrow \left[\begin{array}{c} \Delta s \\ \Delta t \end{array}\right] &= \frac{9}{64} \left[\begin{array}{c} -1 \\ 2 \end{array}\right]. \end{align*}\end{split}\]
>>> nodes1 = np.array([
...     [0.0, 0.0],
...     [2.0, 4.0],
...     [4.0, 0.0],
... ])
>>> curve1 = bezier.Curve(nodes1)
>>> nodes2 = np.array([
...     [2.0, 0.0],
...     [0.0, 3.0],
... ])
>>> curve2 = bezier.Curve(nodes2)
>>> s, t = 0.375, 0.25
>>> new_s, new_t = newton_refine(s, curve1, t, curve2)
>>> 64.0 * (new_s - s)
-9.0
>>> 64.0 * (new_t - t)
18.0

For “typical” curves, we converge to a solution quadratically. This means that the number of correct digits doubles every iteration (until machine precision is reached).

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.25, 2.0],
...     [0.5, -2.0],
...     [0.75, 2.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 1.0],
...     [0.25, 0.5],
...     [0.5, 0.5],
...     [0.75, 0.5],
...     [1.0, 0.0],
... ]))
>>> # The expected intersection is the only real root of
>>> # 28 s^3 - 30 s^2 + 9 s - 1.
>>> rts = np.roots([28, -30, 9, -1])
>>> expected, = rts[rts.imag == 0.0].real
>>> s, t = 0.625, 0.625
>>> np.log2(abs(expected - s))
-4.399...
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-7.901...
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-16.010...
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-32.110...
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-50.0

However, when the intersection occurs at a point of tangency, the convergence becomes linear. This means that the number of correct digits added each iteration is roughly constant.

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 0.5],
...     [1.0, 0.5],
... ]))
>>> expected = 0.5
>>> s, t = 0.375, 0.375
>>> np.log2(abs(expected - s))
-3.0
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-4.0
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-5.0
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-6.0
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-7.0
>>> s, t = newton_refine(s, curve1, t, curve2)
>>> np.log2(abs(expected - s))
-8.0
Parameters:
  • s (float) – Parameter of a near-intersection along curve1.
  • curve1 (Curve) – First curve forming intersection.
  • t (float) – Parameter of a near-intersection along curve2.
  • curve2 (Curve) – Second curve forming intersection.
Returns:

The refined parameters from a single Newton step.

Return type:

Tuple [ float, float ]

bezier._intersection_helpers.segment_intersection(start0, end0, start1, end1)

Determine the intersection of two line segments.

Assumes each line is parametric

\[\begin{split}\begin{alignat*}{2} L_0(s) &= S_0 (1 - s) + E_0 s &&= S_0 + s \Delta_0 \\ L_1(t) &= S_1 (1 - t) + E_1 t &&= S_1 + t \Delta_1. \end{alignat*}\end{split}\]

To solve \(S_0 + s \Delta_0 = S_1 + t \Delta_1\), we use the cross product:

\[\left(S_0 + s \Delta_0\right) \times \Delta_1 = \left(S_1 + t \Delta_1\right) \times \Delta_1 \Longrightarrow s \left(\Delta_0 \times \Delta_1\right) = \left(S_1 - S_0\right) \times \Delta_1.\]

Similarly

\[\Delta_0 \times \left(S_0 + s \Delta_0\right) = \Delta_0 \times \left(S_1 + t \Delta_1\right) \Longrightarrow \left(S_1 - S_0\right) \times \Delta_0 = \Delta_0 \times \left(S_0 - S_1\right) = t \left(\Delta_0 \times \Delta_1\right).\]

Note

Since our points are in \(\mathbf{R}^2\), the “traditional” cross-product in \(\mathbf{R}^3\) will always point in the \(z\) direction, so in the above we mean the \(z\) component of the cross product, rather than the entire vector.

For example, the diagonal lines

\[\begin{split}\begin{align*} L_0(s) &= \left[\begin{array}{c} 0 \\ 0 \end{array}\right] (1 - s) + \left[\begin{array}{c} 2 \\ 2 \end{array}\right] s \\ L_1(t) &= \left[\begin{array}{c} -1 \\ 2 \end{array}\right] (1 - t) + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] t \end{align*}\end{split}\]

intersect at \(L_0\left(\frac{1}{4}\right) = L_1\left(\frac{3}{4}\right) = \frac{1}{2} \left[\begin{array}{c} 1 \\ 1 \end{array}\right]\).

>>> start0 = np.array([[0.0, 0.0]])
>>> end0 = np.array([[2.0, 2.0]])
>>> start1 = np.array([[-1.0, 2.0]])
>>> end1 = np.array([[1.0, 0.0]])
>>> s, t, _ = segment_intersection(start0, end0, start1, end1)
>>> s
0.25
>>> t
0.75

Taking the parallel (but different) lines

\[\begin{split}\begin{align*} L_0(s) &= \left[\begin{array}{c} 1 \\ 0 \end{array}\right] (1 - s) + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] s \\ L_1(t) &= \left[\begin{array}{c} -1 \\ 3 \end{array}\right] (1 - t) + \left[\begin{array}{c} 3 \\ -1 \end{array}\right] t \end{align*}\end{split}\]

we should be able to determine that the lines don’t intersect, but this function is not meant for that check:

>>> start0 = np.array([[1.0, 0.0]])
>>> end0 = np.array([[0.0, 1.0]])
>>> start1 = np.array([[-1.0, 3.0]])
>>> end1 = np.array([[3.0, -1.0]])
>>> _, _, success = segment_intersection(start0, end0, start1, end1)
>>> success
False

Instead, we use parallel_different():

>>> parallel_different(start0, end0, start1, end1)
True
Parameters:
  • start0 (numpy.ndarray) – A 1x2 NumPy array that is the start vector \(S_0\) of the parametric line \(L_0(s)\).
  • end0 (numpy.ndarray) – A 1x2 NumPy array that is the end vector \(E_0\) of the parametric line \(L_0(s)\).
  • start1 (numpy.ndarray) – A 1x2 NumPy array that is the start vector \(S_1\) of the parametric line \(L_1(s)\).
  • end1 (numpy.ndarray) – A 1x2 NumPy array that is the end vector \(E_1\) of the parametric line \(L_1(s)\).
Returns:

Pair of \(s_{\ast}\) and \(t_{\ast}\) such that the lines intersect: \(L_0\left(s_{\ast}\right) = L_1\left(t_{\ast}\right)\) and then a boolean indicating if an intersection was found.

Return type:

Tuple [ float, float, bool ]

bezier._intersection_helpers.parallel_different(start0, end0, start1, end1)

Checks if two parallel lines ever meet.

Meant as a back-up when segment_intersection() fails.

Note

This function assumes but never verifies that the lines are parallel.

In the case that the segments are parallel and lie on the exact same line, finding a unique intersection is not possible. However, if they are parallel but on different lines, then there is a guarantee of no intersection.

In segment_intersection(), we utilized the normal form of the lines (via the cross product):

\[\begin{split}\begin{align*} L_0(s) \times \Delta_0 &\equiv S_0 \times \Delta_0 \\ L_1(t) \times \Delta_1 &\equiv S_1 \times \Delta_1 \end{align*}\end{split}\]

So, we can detect if \(S_1\) is on the first line by checking if

\[S_0 \times \Delta_0 \stackrel{?}{=} S_1 \times \Delta_0.\]

If it is not on the first line, then we are done, the segments don’t meet:

>>> # Line: y = 1
>>> start0 = np.array([[0.0, 1.0]])
>>> end0 = np.array([[1.0, 1.0]])
>>> # Vertical shift up: y = 2
>>> start1 = np.array([[-1.0, 2.0]])
>>> end1 = np.array([[3.0, 2.0]])
>>> parallel_different(start0, end0, start1, end1)
True

If \(S_1\) is on the first line, we want to check that \(S_1\) and \(E_1\) define parameters outside of \(\left[0, 1\right]\). To compute these parameters:

\[L_1(t) = S_0 + s_{\ast} \Delta_0 \Longrightarrow s_{\ast} = \frac{\Delta_0^T \left( L_1(t) - S_0\right)}{\Delta_0^T \Delta_0}.\]

For example, the segments \(\left[0, 1\right]\) and \(\left[\frac{3}{2}, 2\right]\) (via \(S_1 = S_0 + \frac{3}{2} \Delta_0\) and \(E_1 = S_0 + 2 \Delta_0\)) don’t meet:

>>> start0 = np.array([[0.0, 1.0]])
>>> delta0 = np.array([[-1.0, 2.0]])
>>> end0 = start0 + 1.0 * delta0
>>> start1 = start0 + 1.5 * delta0
>>> end1 = start0 + 2.0 * delta0
>>> parallel_different(start0, end0, start1, end1)
True

but if the segments, overlap, like \(\left[0, 1\right]\) and \(\left[-1, \frac{1}{2}\right]\), the segments meet:

>>> start1 = start0 - 1.0 * delta0
>>> end1 = start0 + 0.5 * delta0
>>> parallel_different(start0, end0, start1, end1)
False

Similarly, if the second interval completely contains the first, the segments meet:

>>> start1 = start0 + 3.0 * delta0
>>> end1 = start0 - 2.0 * delta0
>>> parallel_different(start0, end0, start1, end1)
False

Note

This method doesn’t currently allow wiggle room around the desired value, i.e. the two values must be bitwise identical. However, the most “correct” version of this function likely should allow for some round off.

Parameters:
  • start0 (numpy.ndarray) – A 1x2 NumPy array that is the start vector \(S_0\) of the parametric line \(L_0(s)\).
  • end0 (numpy.ndarray) – A 1x2 NumPy array that is the end vector \(E_0\) of the parametric line \(L_0(s)\).
  • start1 (numpy.ndarray) – A 1x2 NumPy array that is the start vector \(S_1\) of the parametric line \(L_1(s)\).
  • end1 (numpy.ndarray) – A 1x2 NumPy array that is the end vector \(E_1\) of the parametric line \(L_1(s)\).
Returns:

Indicating if the lines are different.

Return type:

bool

class bezier.curve.Curve(nodes, start=0.0, end=1.0, root=None, _copy=True)

Bases: bezier._base.Base

Represents a Bézier curve.

We take the traditional definition: a Bézier curve is a mapping from \(s \in \left[0, 1\right]\) to convex combinations of points \(v_0, v_1, \ldots, v_n\) in some vector space:

\[B(s) = \sum_{j = 0}^n \binom{n}{j} s^j (1 - s)^{n - j} \cdot v_j\]
>>> import bezier
>>> nodes = np.array([
...     [0.0  , 0.0],
...     [0.625, 0.5],
...     [1.0  , 0.5],
... ])
>>> curve = bezier.Curve(nodes)
>>> curve
<Curve (degree=2, dimension=2)>
Parameters:
  • nodes (numpy.ndarray) – The nodes in the curve. The rows represent each node while the columns are the dimension of the ambient space.
  • start (Optional [ float ]) – The beginning of the sub-interval that this curve represents.
  • end (Optional [ float ]) – The end of the sub-interval that this curve represents.
  • root (Optional [ Curve ]) – The root curve that contains this current curve.
  • _copy (bool) – Flag indicating if the nodes should be copied before being stored. Defaults to True since callers may freely mutate nodes after passing in.
__repr__()

Representation of current object.

Returns:Object representation.
Return type:str
length

float: The length of the current curve.

start

float: Start of sub-interval this curve represents.

This value is used to track the current curve in the re-parameterization / subdivision process. The curve is still defined on the unit interval, but this value illustrates how this curve relates to a “parent” curve. For example:

>>> nodes = np.array([
...     [0.0, 0.0],
...     [1.0, 2.0],
... ])
>>> curve = bezier.Curve(nodes)
>>> curve
<Curve (degree=1, dimension=2)>
>>> left, right = curve.subdivide()
>>> left
<Curve (degree=1, dimension=2, start=0, end=0.5)>
>>> right
<Curve (degree=1, dimension=2, start=0.5, end=1)>
>>> _, mid_right = left.subdivide()
>>> mid_right
<Curve (degree=1, dimension=2, start=0.25, end=0.5)>
>>> mid_right.nodes
array([[ 0.25, 0.5 ],
       [ 0.5 , 1.  ]])
end

float: End of sub-interval this curve represents.

See start for more information.

root

Curve: The “root” curve that contains the current curve.

This indicates that the current curve is a section of the “root” curve. For example:

>>> _, right = curve.subdivide()
>>> right
<Curve (degree=2, dimension=2, start=0.5, end=1)>
>>> right.root is curve
True
>>> right.evaluate(0.0) == curve.evaluate(0.5)
array([ True, True], dtype=bool)
>>>
>>> mid_left, _ = right.subdivide()
>>> mid_left
<Curve (degree=2, dimension=2, start=0.5, end=0.75)>
>>> mid_left.root is curve
True
>>> mid_left.evaluate(1.0) == curve.evaluate(0.75)
array([ True, True], dtype=bool)
evaluate(s)

Evaluate \(B(s)\) along the curve.

See evaluate_multi() for more details.

>>> nodes = np.array([
...     [0.0  , 0.0],
...     [0.625, 0.5],
...     [1.0  , 0.5],
... ])
>>> curve = bezier.Curve(nodes)
>>> curve.evaluate(0.75)
array([ 0.796875, 0.46875 ])
Parameters:s (float) – Parameter along the curve.
Returns:The point on the curve (as a one dimensional NumPy array).
Return type:numpy.ndarray
evaluate_multi(s_vals)

Evaluate \(B(s)\) for multiple points along the curve.

This is done by first evaluating each member of the Bernstein basis at each value in s_vals and then applying those to the control points for the current curve.

This is done instead of using de Casteljau’s algorithm. Implementing de Casteljau is problematic because it requires a choice between one of two methods:

  • vectorize operations of the form \((1 - s)v + s w\), which requires a copy of the curve’s control points for each value in s_vals
  • avoid vectorization and compute each point in serial

Instead, we can use vectorized operations to build up the Bernstein basis values.

>>> nodes = np.array([
...     [0.0, 0.0, 0.0],
...     [1.0, 2.0, 3.0],
... ])
>>> curve = bezier.Curve(nodes)
>>> curve
<Curve (degree=1, dimension=3)>
>>> s_vals = np.linspace(0.0, 1.0, 5)
>>> curve.evaluate_multi(s_vals)
array([[ 0.  , 0.  , 0.  ],
       [ 0.25, 0.5 , 0.75],
       [ 0.5 , 1.  , 1.5 ],
       [ 0.75, 1.5 , 2.25],
       [ 1.  , 2.  , 3.  ]])
Parameters:s_vals (numpy.ndarray) – Parameters along the curve (as a 1D array).
Returns:The points on the curve. As a two dimensional NumPy array, with the rows corresponding to each s value and the columns to the dimension.
Return type:numpy.ndarray
plot(num_pts, ax=None, show=False)

Plot the current curve.

Parameters:
Returns:

The axis containing the plot. This may be a newly created axis.

Return type:

matplotlib.artist.Artist

Raises:

NotImplementedError – If the curve’s dimension is not 2.

subdivide()

Split the curve \(\gamma(s)\) into a left and right half.

Takes the interval \(\left[0, 1\right]\) and splits the curve into \(\gamma_1 = \gamma\left(\left[0, \frac{1}{2}\right]\right)\) and \(\gamma_2 = \gamma\left(\left[\frac{1}{2}, 1\right]\right)\). In order to do this, also reparameterizes the curve, hence the resulting left and right halves have new nodes.

>>> nodes = np.array([
...     [0.0 , 0.0],
...     [1.25, 3.0],
...     [2.0 , 1.0],
... ])
>>> curve = bezier.Curve(nodes)
>>> left, right = curve.subdivide()
>>> left
<Curve (degree=2, dimension=2, start=0, end=0.5)>
>>> left.nodes
array([[ 0.   , 0.   ],
       [ 0.625, 1.5  ],
       [ 1.125, 1.75 ]])
>>> right
<Curve (degree=2, dimension=2, start=0.5, end=1)>
>>> right.nodes
array([[ 1.125, 1.75 ],
       [ 1.625, 2.   ],
       [ 2.   , 1.   ]])
Returns:The left and right sub-curves.
Return type:Tuple [ Curve, Curve ]
intersect(other)

Find the points of intersection with another curve.

>>> nodes1 = np.array([
...     [0.0  , 0.0 ],
...     [0.375, 0.75 ],
...     [0.75 , 0.375],
... ])
>>> curve1 = bezier.Curve(nodes1)
>>> nodes2 = np.array([
...     [0.5, 0.0],
...     [0.5, 3.0],
... ])
>>> curve2 = bezier.Curve(nodes2)
>>> curve1.intersect(curve2)
array([[ 0.5, 0.5]])
Parameters:

other (Curve) – Other curve to intersect with.

Returns:

Array of intersection points (possibly empty).

Return type:

numpy.ndarray

Raises:
__eq__(other)

Check equality against another shape.

Returns:Boolean indicating if the shapes are the same.
Return type:bool
__ne__(other)

Check inequality against another shape.

Returns:Boolean indicating if the shapes are not the same.
Return type:bool
copy()

Make a copy of the current shape.

Returns:Instance of the current shape.
degree

int: The degree of the current shape.

dimension

int: The dimension that the shape lives in.

For example, if the shape lives in \(\mathbf{R}^3\), then the dimension is 3.

elevate()

Return a degree-elevated version of the current curve.

Does this by converting the current nodes \(v_0, \ldots, v_n\) to new nodes \(w_0, \ldots, w_{n + 1}\) where

\[\begin{split}\begin{align*} w_0 &= v_0 \\ w_j &= \frac{j}{n + 1} v_{j - 1} + \frac{n + 1 - j}{n + 1} v_j \\ w_{n + 1} &= v_n \end{align*}\end{split}\]
Returns:The degree-elevated curve.
Return type:Curve
nodes

numpy.ndarray: The nodes that define the current shape.

bezier.surface module

Helper for Bézier Surfaces / Triangles.

class bezier.surface.Surface(nodes, _copy=True)

Bases: bezier._base.Base

Represents a Bézier surface.

We define a Bézier triangle as a mapping from the unit simplex in 2D (i.e. the unit triangle) onto a surface in an arbitrary dimension. We use barycentric coordinates

\[\lambda_1 = 1 - s - t, \lambda_2 = s, \lambda_3 = t\]

for points in

\[\left\{(s, t) \mid 0 \leq s, t, s + t \leq 1\right\}.\]

As with curves, using these weights we get convex combinations of points \(v_{i, j, k}\) in some vector space:

\[B\left(\lambda_1, \lambda_2, \lambda_3\right) = \sum_{i + j + k = d} \binom{d}{i \, j \, k} \lambda_1^i \lambda_2^j \lambda_3^k \cdot v_{i, j, k}\]

Note

We assume the nodes are ordered from left-to-right and from bottom-to-top. So for example, the linear triangle:

(0,0,1)

(1,0,0)  (0,1,0)

is ordered as

\[\begin{split}\left[\begin{array}{c c c} v_{1,0,0} & v_{0,1,0} & v_{0,0,1} \end{array}\right]^T\end{split}\]

the quadratic triangle:

(0,0,2)

(1,0,1)  (0,1,1)

(2,0,0)  (1,1,0)  (0,2,0)

is ordered as

\[\begin{split}\left[\begin{array}{c c c c c c} v_{2,0,0} & v_{1,1,0} & v_{0,2,0} & v_{1,0,1} & v_{0,1,1} & v_{0,0,2} \end{array}\right]^T\end{split}\]

the cubic triangle:

(0,0,3)

(1,0,2)  (0,1,2)

(2,0,1)  (1,1,1)  (0,2,1)

(3,0,0)  (2,1,0)  (1,2,0)  (0,3,0)

is ordered as

\[\begin{split}\left[\begin{array}{c c c c c c c c c c} v_{3,0,0} & v_{2,1,0} & v_{1,2,0} & v_{0,3,0} & v_{2,0,1} & v_{1,1,1} & v_{0,2,1} & v_{1,0,2} & v_{0,1,2} & v_{0,0,3} \end{array}\right]^T\end{split}\]

and so on.

>>> import bezier
>>> nodes = np.array([
...     [0.0 , 0.0 ],
...     [1.0 , 0.25],
...     [0.25, 1.0 ],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface
<Surface (degree=1, dimension=2)>
Parameters:nodes (numpy.ndarray) – The nodes in the surface. The rows represent each node while the columns are the dimension of the ambient space.
area

float: The area of the current surface.

Raises:NotImplementedError – If the area isn’t already cached.
edges

tuple: The edges of the surface.

>>> nodes = np.array([
...     [0.0   ,  0.0   ],
...     [0.5   , -0.1875],
...     [1.0   ,  0.0   ],
...     [0.1875,  0.5   ],
...     [0.625 ,  0.625 ],
...     [0.0   ,  1.0   ],
... ])
>>> surface = bezier.Surface(nodes)
>>> edge1, _, _ = surface.edges
>>> edge1
<Curve (degree=2, dimension=2)>
>>> edge1.nodes
array([[ 0.  ,  0.    ],
       [ 0.5 , -0.1875],
       [ 1.  ,  0.    ]])
Returns:The edges of the surface.
Return type:Tuple [ Curve, Curve, Curve ]
evaluate_barycentric(lambda1, lambda2, lambda3)

Compute a point on the surface.

Evaluates \(B\left(\lambda_1, \lambda_2, \lambda_3\right)\).

>>> nodes = np.array([
...     [0.0 , 0.0 ],
...     [1.0 , 0.25],
...     [0.25, 1.0 ],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface.evaluate_barycentric(0.125, 0.125, 0.75)
array([ 0.3125 , 0.78125])

However, this can’t be used for points outside the reference triangle:

>>> surface.evaluate_barycentric(-0.25, 0.75, 0.5)
Traceback (most recent call last):
  ...
ValueError: ('Parameters must be positive', -0.25, 0.75, 0.5)

or for non-Barycentric coordinates;

>>> surface.evaluate_barycentric(0.25, 0.25, 0.25)
Traceback (most recent call last):
  ...
ValueError: ('Values do not sum to 1', 0.25, 0.25, 0.25)
Parameters:
  • lambda1 (float) – Parameter along the reference triangle.
  • lambda2 (float) – Parameter along the reference triangle.
  • lambda3 (float) – Parameter along the reference triangle.
Returns:

The point on the surface (as a one dimensional NumPy array).

Return type:

numpy.ndarray

Raises:
  • ValueError – If the weights are not valid barycentric coordinates, i.e. they don’t sum to 1.
  • ValueError – If some weights are negative.
evaluate_cartesian(s, t)

Compute a point on the surface.

Evaluates \(B\left(1 - s - t, s, t\right)\) by calling evaluate_barycentric().

Parameters:
  • s (float) – Parameter along the reference triangle.
  • t (float) – Parameter along the reference triangle.
Returns:

The point on the surface (as a one dimensional NumPy array).

Return type:

numpy.ndarray

evaluate_multi(param_vals)

Compute multiple points on the surface.

If param_vals has two columns, this method treats them as Cartesian:

>>> nodes = np.array([
...     [ 0., 0. ],
...     [ 2., 1. ],
...     [-3., 2. ],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface
<Surface (degree=1, dimension=2)>
>>> param_vals = np.array([
...     [0.0, 0.0],
...     [1.0, 0.0],
...     [0.5, 0.5],
... ])
>>> surface.evaluate_multi(param_vals)
array([[ 0. , 0. ],
       [ 2. , 1. ],
       [-0.5, 1.5]])

and if param_vals has three columns, treats them as Barycentric:

>>> nodes = np.array([
...     [ 0. , 0.  ],
...     [ 1. , 0.75],
...     [ 2. , 1.  ],
...     [-1.5, 1.  ],
...     [-0.5, 1.5 ],
...     [-3. , 2.  ],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface
<Surface (degree=2, dimension=2)>
>>> param_vals = np.array([
...     [0.   , 0.25, 0.75 ],
...     [1.   , 0.  , 0.   ],
...     [0.25 , 0.5 , 0.25 ],
...     [0.375, 0.25, 0.375],
... ])
>>> surface.evaluate_multi(param_vals)
array([[-1.75  , 1.75    ],
       [ 0.    , 0.      ],
       [ 0.25  , 1.0625  ],
       [-0.625 , 1.046875]])

Note

This currently just uses evaluate_cartesian() and evaluate_barycentric() so is less performant than it could be.

Parameters:

param_vals (numpy.ndarray) – Array of parameter values (as a 2D array).

Returns:

The point on the surface.

Return type:

numpy.ndarray

Raises:
  • ValueError – If param_vals is not a 2D array.
  • ValueError – If param_vals doesn’t have 2 or 3 columns.
plot(pts_per_edge, ax=None, show=False)

Plot the current surface.

Parameters:
Returns:

The axis containing the plot. This may be a newly created axis.

Return type:

matplotlib.artist.Artist

Raises:

NotImplementedError – If the surface’s dimension is not 2.

subdivide()

Split the surface into four sub-surfaces.

Takes the reference triangle

\[T = \left\{(s, t) \mid 0 \leq s, t, s + t \leq 1\right\}\]

and splits it into four sub-triangles

\[\begin{split}\begin{align*} A &= \left\{(s, t) \mid 0 \leq s, t, s + t \leq \frac{1}{2}\right\} \\ B &= -A + \left(\frac{1}{2}, \frac{1}{2}\right) \\ C &= A + \left(\frac{1}{2}, 0\right) \\ D &= A + \left(0, \frac{1}{2}\right). \end{align*}\end{split}\]

These are the lower left (\(A\)), central (\(B\)), lower right (\(C\)) and upper left (\(D\)) sub-triangles.

>>> nodes = np.array([
...     [ 0. , 0. ],
...     [ 2. , 0. ],
...     [ 0. , 4. ],
... ])
>>> surface = bezier.Surface(nodes)
>>> _, _, _, sub_surface_d = surface.subdivide()
>>> sub_surface_d
<Surface (degree=1, dimension=2)>
>>> sub_surface_d.nodes
array([[ 0., 2.],
       [ 1., 2.],
       [ 0., 4.]])
Returns:The lower left, central, lower right and upper left sub-surfaces (in that order).
Return type:Tuple [ Surface, Surface, Surface, Surface ]
is_valid

bool: Flag indicating if the surface no singularites.

This checks if the Jacobian of the map from the reference triangle is nonzero. For example, a linear “surface” with collinear points is invalid:

>>> nodes = np.array([
...     [0.0, 0.0],
...     [1.0, 1.0],
...     [2.0, 2.0],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface.is_valid
False

while a quadratic surface with one straight side:

>>> nodes = np.array([
...     [ 0.0  , 0.0  ],
...     [ 0.5  , 0.125],
...     [ 1.0  , 0.0  ],
...     [-0.125, 0.5  ],
...     [ 0.5  , 0.5  ],
...     [ 0.0  , 1.0  ],
... ])
>>> surface = bezier.Surface(nodes)
>>> surface.is_valid
True
__eq__(other)

Check equality against another shape.

Returns:Boolean indicating if the shapes are the same.
Return type:bool
__ne__(other)

Check inequality against another shape.

Returns:Boolean indicating if the shapes are not the same.
Return type:bool
__repr__()

Representation of current object.

Returns:Object representation.
Return type:str
copy()

Make a copy of the current shape.

Returns:Instance of the current shape.
degree

int: The degree of the current shape.

dimension

int: The dimension that the shape lives in.

For example, if the shape lives in \(\mathbf{R}^3\), then the dimension is 3.

nodes

numpy.ndarray: The nodes that define the current shape.

Curve-Curve Intersection

The problem of intersecting two curves is a difficult one in computational geometry. The Curve.intersect() method uses a combination of curve subdivision, bounding box intersection, and curve approximation (by lines) to find intersections.

Curve-Line Intersection

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 0.375],
...     [1.0, 0.375],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.25 , 0.375],
       [ 0.75 , 0.375]])
_images/test_curves1_and_8.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.5, 0.0 ],
...     [0.5, 0.75],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.5, 0.5]])
_images/test_curves1_and_9.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [4.5, 9.0],
...     [9.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 8.0],
...     [6.0, 0.0],
... ]))
>>> curve1.intersect(curve2)
array([[ 3., 4.]])
_images/test_curves10_and_11.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.375],
...     [1.0, 0.375],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.5, 0.0 ],
...     [0.5, 0.75],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.5 , 0.375]])
_images/test_curves8_and_9.png
>>> curve1 = bezier.Curve(np.array([
...     [-1.0, 1.0],
...     [ 0.5, 0.5],
...     [ 0.0, 2.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [ 0.5 , 0.5 ],
...     [-0.25, 1.25],
... ]))
>>> curve1.intersect(curve2)
array([[ 0., 1.]])
_images/test_curves29_and_30.png

Curved Intersections

For curves which intersect at exact floating point numbers, we can typically compute the intersection with zero error:

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0,  0.75],
...     [0.5, -0.25],
...     [1.0,  0.75],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.25 , 0.375],
       [ 0.75 , 0.375]])
_images/test_curves1_and_5.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [1.5, 3.0],
...     [3.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [ 3.0  ,  1.5    ],
...     [ 2.625, -0.90625],
...     [-0.75 ,  2.4375 ],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.75  , 1.125  ],
       [ 2.625 , 0.65625]])
_images/test_curves3_and_4.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0  , 0.0  ],
...     [0.375, 0.75 ],
...     [0.75 , 0.375],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.25 , 0.5625],
...     [0.625, 0.1875],
...     [1.0  , 0.9375],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.375 , 0.46875],
       [ 0.625 , 0.46875]])
_images/test_curves14_and_16.png

Even for curves which don’t intersect at exact floating point numbers, we can compute the intersection to machine precision:

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [1.125,  0.5],
...     [0.625, -0.5],
...     [0.125,  0.5],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> sq31 = np.sqrt(31.0)
>>> expected = np.array([
...     [36 - 4 * sq31, 16 + sq31],
...     [36 + 4 * sq31, 16 - sq31],
... ]) / 64.0
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-54.0
_images/test_curves1_and_2.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 0.265625],
...     [0.5, 0.234375],
...     [1.0, 0.265625],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> sq33 = np.sqrt(33.0)
>>> expected = np.array([
...     [33 - 4 * sq33, 17],
...     [33 + 4 * sq33, 17],
... ]) / 66.0
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-54.0
_images/test_curves1_and_7.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0 ,  0.0],
...     [0.25,  2.0],
...     [0.5 , -2.0],
...     [0.75,  2.0],
...     [1.0 ,  0.0],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> sq7 = np.sqrt(7.0)
>>> expected = np.array([
...     [7 - sq7, 6],
...     [7 + sq7, 6],
...     [      0, 0],
...     [     14, 0],
... ]) / 14.0
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-54.0
_images/test_curves1_and_13.png
>>> curve1 = bezier.Curve(np.array([
...     [-0.125, -0.28125],
...     [ 0.5  ,  1.28125],
...     [ 1.125, -0.28125],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [ 1.5625, -0.0625],
...     [-1.5625,  0.25  ],
...     [ 1.5625,  0.5625],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> sq5 = np.sqrt(5.0)
>>> expected = np.array([
...     [6 - 2 * sq5, 5 - sq5],
...     [          4, 6      ],
...     [         16, 0      ],
...     [6 + 2 * sq5, 5 + sq5],
... ]) / 16.0
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-50.415...
_images/test_curves21_and_22.png

For higher degree intersections, the error starts to get a little larger.

>>> curve1 = bezier.Curve(np.array([
...     [0.25 , 0.625],
...     [0.625, 0.25 ],
...     [1.0  , 1.0  ],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0 , 0.5],
...     [0.25, 1.0],
...     [0.75, 1.5],
...     [1.0 , 0.5],
... ]))
>>> 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)
>>> x_val = (3 * s_val + 1) / 4
>>> y_val = (9 * s_val * s_val - 6 * s_val + 5) / 8
>>> expected = np.array([
...     [x_val, y_val],
... ])
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-49.678...
_images/test_curves15_and_25.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 8.0],
...     [6.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.375, 7.0],
...     [2.125, 8.0],
...     [3.875, 0.0],
...     [5.625, 1.0],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> sq7 = np.sqrt(7.0)
>>> expected = np.array([
...     [           72, 96           ],
...     [72 - 21 * sq7, 96 + 28 * sq7],
...     [72 + 21 * sq7, 96 - 28 * sq7],
... ]) / 24.0
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-50.0
_images/test_curves11_and_26.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.375],
...     [1.0, 0.375],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.125, 0.25  ],
...     [0.375, 0.75  ],
...     [0.625, 0.0   ],
...     [0.875, 0.1875],
... ]))
>>> intersections = curve1.intersect(curve2)
>>> s_val1, s_val2, _ = np.sort(np.roots(
...     [17920, -29760, 13512, -1691]))
>>> expected = np.array([
...     [s_val2, 0.375],
...     [s_val1, 0.375],
... ])
>>> max_err = np.max(np.abs(intersections - expected))
>>> np.log2(max_err)
-50.678...
_images/test_curves8_and_27.png

Intersections at Endpoints

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [1.0,  0.0],
...     [1.5, -1.0],
...     [2.0,  0.0],
... ]))
>>> curve1.intersect(curve2)
array([[ 1., 0.]])
_images/test_curves1_and_18.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [2.0, 0.0],
...     [1.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve1.intersect(curve2)
array([[ 1., 0.]])
_images/test_curves1_and_19.png
>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [4.5, 9.0],
...     [9.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [11.0,  8.0],
...     [ 7.0, 10.0],
...     [ 3.0,  4.0],
... ]))
>>> curve1.intersect(curve2)
array([[ 3., 4.]])
_images/test_curves10_and_17.png

Detecting Self-Intersections

>>> curve1 = bezier.Curve(np.array([
...     [ 0.0 , 0.0  ],
...     [-1.0 , 0.0  ],
...     [ 1.0 , 1.0  ],
...     [-0.75, 1.625],
... ]))
>>> left, right = curve1.subdivide()
>>> left.intersect(right)
array([[-0.09375 , 0.578125]])
_images/test_curve12_self_crossing.png

Limitations

Intersections that occur at points of tangency are in general problematic. For example, consider

\[\begin{split}B_1(s) = \left[ \begin{array}{c} s \\ 2s(1 - s)\end{array}\right], \quad B_2(t) = \left[ \begin{array}{c} t \\ t^2 + (1 - t)^2 \end{array}\right]\end{split}\]

The first curve is the zero set of \(y - 2x(1 - x)\), so plugging in the second curve gives

\[0 = t^2 + (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

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.0, 1.0],
...     [0.5, 0.0],
...     [1.0, 1.0],
... ]))
>>> curve1.intersect(curve2)
array([[ 0.5, 0.5]])
_images/test_curves1_and_6.png

However this library mostly avoids (for now) computing tangent intersections. For example, the curves

_images/test_curves14_and_15.png

have a tangent intersection that this library fails to compute:

>>> curve1 = bezier.Curve(np.array([
...     [0.0  , 0.0  ],
...     [0.375, 0.75 ],
...     [0.75 , 0.375],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.25 , 0.625],
...     [0.625, 0.25 ],
...     [1.0  , 1.0  ],
... ]))
>>> 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:

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [4.5, 9.0],
...     [9.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [3.0, 4.5],
...     [8.0, 4.5],
... ]))
>>> curve1.intersect(curve2)
array([[ 4.5, 4.5]])
_images/test_curves10_and_23.png

but even by rotating an intersection (from above) that we know works

_images/test_curves28_and_29.png

we still see a failure

>>> curve1 = bezier.Curve(np.array([
...     [ 0.0, 0.0],
...     [-0.5, 1.5],
...     [ 1.0, 1.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [-1.0, 1.0],
...     [ 0.5, 0.5],
...     [ 0.0, 2.0],
... ]))
>>> 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

_images/test_curves1_and_24.png

the library fails as well

>>> curve1 = bezier.Curve(np.array([
...     [0.0, 0.0],
...     [0.5, 1.0],
...     [1.0, 0.0],
... ]))
>>> curve2 = bezier.Curve(np.array([
...     [0.25,  0.375],
...     [0.75,  0.875],
...     [1.25, -0.625],
... ]))
>>> curve1.intersect(curve2)
Traceback (most recent call last):
  ...
NotImplementedError: The number of candidate intersections is too high.

bezier is a Python helper for Bézier Curves, Triangles, and Higher Order Objects.

This library provides:

Installing

bezier can be installed with pip:

$ pip install --upgrade bezier

bezier is open-source, so you can alternatively grab the source code from GitHub and install from source.

Usage

The Bézier Package documentation provides API-level documentation.

License

bezier is made available under the Apache License. For more details, see LICENSE