bezier.hazmat.algebraic_intersection module

Helpers for intersecting Bézier curves via algebraic methods.

Primarily helps implicitize Bézier curves.

Uses the resultant to evaluate the implicitized algebraic curve. In order to do this on Bézier curves without translating to a power basis, we utilize the work of Farouki and Rajan to compute a modified Sylvester determinant.

Given two parametric curves \((x_1(s), y_1(s))\) and \((x_2(t), y_2(t))\), we can determine an “intersection polynomial” for both \(s\) and \(t\). For example, by implicitizing the first curve, we determine \(f_1(x, y)\) and plugging the second curve into this we find

\[g(t) = f_1\left(x_2(t), y_2(t)\right) = 0\]

is the “intersection polynomial” for \(t\).

bezier.hazmat.algebraic_intersection.evaluate(nodes, x_val, y_val)

Evaluate the implicitized bivariate polynomial containing the curve.

Assumes the algebraic curve containing \(B(s, t)\) is given by \(f(x, y) = 0\). This function evaluates \(f(x, y)\).

Note

This assumes, but doesn’t check, that nodes has 2 rows.

Note

This assumes, but doesn’t check, that nodes is not degree-elevated. If it were degree-elevated, then the Sylvester matrix will always have zero determinant.

Parameters:

  • nodes (numpy.ndarray) – 2 x N array of nodes in a curve.

  • x_val (float) – x-coordinate for evaluation.

  • y_val (float) – y-coordinate for evaluation.

Returns:

The computed value of \(f(x, y)\).

Return type:

float

Raises:
bezier.hazmat.algebraic_intersection.eval_intersection_polynomial(nodes1, nodes2, t)

Evaluates a parametric curve on an implicitized algebraic curve.

Uses evaluate() to evaluate \(f_1(x, y)\), the implicitization of nodes1. Then plugs t into the second parametric curve to get an x- and y-coordinate and evaluate the intersection polynomial:

\[g(t) = f_1\left(x_2(t), y_2(t)\right)\]
Parameters:

  • nodes1 (numpy.ndarray) – The nodes in the first curve.

  • nodes2 (numpy.ndarray) – The nodes in the second curve.

  • t (float) – The parameter along nodes2 where we evaluate the function.

Returns:

The computed value of \(f_1(x_2(t), y_2(t))\).

Return type:

float

bezier.hazmat.algebraic_intersection.to_power_basis(nodes1, nodes2)

Compute the coefficients of an intersection polynomial.

Note

This requires that the degree of the curve given by nodes1 is less than or equal to the degree of that given by nodes2.

Parameters:
Returns:

Array of coefficients.

Return type:

numpy.ndarray

Raises:

NotImplementedError – If the degree pair is not 1-1, 1-2, 1-3, 1-4, 2-2, 2-3, 2-4 or 3-3.

bezier.hazmat.algebraic_intersection.polynomial_norm(coeffs)

Computes \(L_2\) norm of polynomial on \(\left[0, 1\right]\).

We have

\[\left\langle f, f \right\rangle = \sum_{i, j} \int_0^1 c_i c_j x^{i + j} \, dx = \sum_{i, j} \frac{c_i c_j}{i + j + 1} = \sum_{i} \frac{c_i^2}{2 i + 1} + 2 \sum_{j > i} \frac{c_i c_j}{i + j + 1}.\]
Parameters:

coeffs (numpy.ndarray) – d + 1-array of coefficients in monomial / power basis.

Returns:

The \(L_2\) norm of the polynomial.

Return type:

float

bezier.hazmat.algebraic_intersection.normalize_polynomial(coeffs, threshold=9.094947017729282e-13)

Normalizes a polynomial in the \(L_2\) sense.

Does so on the interval \(\left[0, 1\right]\) via polynomial_norm().

Parameters:

  • coeffs (numpy.ndarray) – d + 1-array of coefficients in monomial / power basis.

  • threshold (Optionalfloat) – The point \(\tau\) below which a polynomial will be considered to be numerically equal to zero, applies to all \(f\) with \(\| f \|_{L_2} < \tau\).

Returns:

The normalized polynomial.

Return type:

numpy.ndarray

bezier.hazmat.algebraic_intersection.bernstein_companion(coeffs)

Compute a companion matrix for a polynomial in Bernstein basis.

Note

This assumes the caller passes in a 1D array but does not check.

This takes the polynomial

\[f(s) = \sum_{j = 0}^n b_{j, n} \cdot C_j.\]

and uses the variable \(\sigma = \frac{s}{1 - s}\) to rewrite as

\[f(s) = (1 - s)^n \sum_{j = 0}^n \binom{n}{j} C_j \sigma^j.\]

This converts the Bernstein coefficients \(C_j\) into “generalized Bernstein” coefficients \(\widetilde{C_j} = \binom{n}{j} C_j\).

Parameters:

coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis.

Returns:

A triple of

  • 2D NumPy array with the companion matrix.

  • the “full degree” based on the size of coeffs

  • the “effective degree” determined by the number of leading zeros.

Return type:

Tuplenumpy.ndarray, int, int

bezier.hazmat.algebraic_intersection.bezier_roots(coeffs)

Compute polynomial roots from a polynomial in the Bernstein basis.

Note

This assumes the caller passes in a 1D array but does not check.

This takes the polynomial

\[f(s) = \sum_{j = 0}^n b_{j, n} \cdot C_j.\]

and uses the variable \(\sigma = \frac{s}{1 - s}\) to rewrite as

\[\begin{split}\begin{align*} f(s) &= (1 - s)^n \sum_{j = 0}^n \binom{n}{j} C_j \sigma^j \\ &= (1 - s)^n \sum_{j = 0}^n \widetilde{C_j} \sigma^j. \end{align*}\end{split}\]

Then it uses an eigenvalue solver to find the roots of

\[g(\sigma) = \sum_{j = 0}^n \widetilde{C_j} \sigma^j\]

and convert them back into roots of \(f(s)\) via \(s = \frac{\sigma}{1 + \sigma}\).

For example, consider

\[\begin{split}\begin{align*} f_0(s) &= 2 (2 - s)(3 + s) \\ &= 12(1 - s)^2 + 11 \cdot 2s(1 - s) + 8 s^2 \end{align*}\end{split}\]

First, we compute the companion matrix for

\[g_0(\sigma) = 12 + 22 \sigma + 8 \sigma^2\]
>>> import numpy as np
>>> coeffs0 = np.asfortranarray([12.0, 11.0, 8.0])
>>> companion0, _, _ = bernstein_companion(coeffs0)
>>> companion0
array([[-2.75, -1.5 ],
       [ 1.  ,  0.  ]])

then take the eigenvalues of the companion matrix:

>>> sigma_values0 = np.linalg.eigvals(companion0)
>>> sigma_values0
array([-2.  , -0.75])

after transforming them, we have the roots of \(f(s)\):

>>> sigma_values0 / (1.0 + sigma_values0)
array([ 2., -3.])
>>> bezier_roots(coeffs0)
array([ 2., -3.])

In cases where \(s = 1\) is a root, the lead coefficient of \(g\) would be \(0\), so there is a reduction in the companion matrix.

\[\begin{split}\begin{align*} f_1(s) &= 6 (s - 1)^2 (s - 3) (s - 5) \\ &= 90 (1 - s)^4 + 33 \cdot 4s(1 - s)^3 + 8 \cdot 6s^2(1 - s)^2 \end{align*}\end{split}\]
>>> coeffs1 = np.asfortranarray([90.0, 33.0, 8.0, 0.0, 0.0])
>>> companion1, degree1, effective_degree1 = bernstein_companion(
...     coeffs1)
>>> companion1
array([[-2.75 , -1.875],
       [ 1.   ,  0.   ]])
>>> degree1
4
>>> effective_degree1
2

so the roots are a combination of the roots determined from \(s = \frac{\sigma}{1 + \sigma}\) and the number of factors of \((1 - s)\) (i.e. the difference between the degree and the effective degree):

>>> bezier_roots(coeffs1)
array([3., 5., 1., 1.])

In some cases, a polynomial is represented with an “elevated” degree:

\[\begin{split}\begin{align*} f_2(s) &= 3 (s^2 + 1) \\ &= 3 (1 - s)^3 + 3 \cdot 3s(1 - s)^2 + 4 \cdot 3s^2(1 - s) + 6 s^3 \end{align*}\end{split}\]

This results in a “point at infinity” \(\sigma = -1 \Longleftrightarrow s = \infty\):

>>> coeffs2 = np.asfortranarray([3.0, 3.0, 4.0, 6.0])
>>> companion2, _, _ = bernstein_companion(coeffs2)
>>> companion2
array([[-2. , -1.5, -0.5],
       [ 1. ,  0. ,  0. ],
       [ 0. ,  1. ,  0. ]])
>>> sigma_values2 = np.linalg.eigvals(companion2)
>>> sigma_values2
array([-1. +0.j , -0.5+0.5j, -0.5-0.5j])

so we drop any values \(\sigma\) that are sufficiently close to \(-1\):

>>> expected2 = np.asfortranarray([1.0j, -1.0j])
>>> roots2 = bezier_roots(coeffs2)
>>> np.allclose(expected2, roots2, rtol=8 * machine_eps, atol=0.0)
True
Parameters:

coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis.

Returns:

A 1D array containing the roots.

Return type:

numpy.ndarray

bezier.hazmat.algebraic_intersection.lu_companion(top_row, value)

Compute an LU-factored \(C - t I\) and its 1-norm.

Note

The output of this function is intended to be used with dgecon from LAPACK. dgecon expects both the 1-norm of the matrix and expects the matrix to be passed in an already LU-factored form (via dgetrf).

The companion matrix \(C\) is given by the top_row, for example, the polynomial \(t^3 + 3 t^2 - t + 2\) has a top row of -3, 1, -2 and the corresponding companion matrix is:

\[\begin{split}\left[\begin{array}{c c c} -3 & 1 & -2 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right]\end{split}\]

After doing a full cycle of the rows (shifting the first to the last and moving all other rows up), row reduction of \(C - t I\) yields

\[\begin{split}\left[\begin{array}{c c c} 1 & -t & 0 \\ 0 & 1 & -t \\ -3 - t & 1 & -2 \end{array}\right] = \left[\begin{array}{c c c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -3 - t & 1 + t(-3 - t) & 1 \end{array}\right] \left[\begin{array}{c c c} 1 & -t & 0 \\ 0 & 1 & -t \\ 0 & 0 & -2 + t(1 + t(-3 - t)) \end{array}\right]\end{split}\]

and in general, the terms in the bottom row correspond to the intermediate values involved in evaluating the polynomial via Horner’s method.

>>> top_row = np.asfortranarray([-3.0, 1.0, -2.0])
>>> t_val = 0.5
>>> lu_mat, one_norm = lu_companion(top_row, t_val)
>>> lu_mat
array([[ 1.   , -0.5  ,  0.   ],
       [ 0.   ,  1.   , -0.5  ],
       [-3.5  , -0.75 , -2.375]])
>>> one_norm
4.5
>>> l_mat = np.tril(lu_mat, k=-1) + np.eye(3)
>>> u_mat = np.triu(lu_mat)
>>> a_mat = l_mat.dot(u_mat)
>>> a_mat
array([[ 1. , -0.5,  0. ],
       [ 0. ,  1. , -0.5],
       [-3.5,  1. , -2. ]])
>>> np.linalg.norm(a_mat, ord=1)
4.5
Parameters:

  • top_row (numpy.ndarray) – 1D array, top row of companion matrix.

  • value (float) – The \(t\) value used to form \(C - t I\).

Returns:

Pair of

  • 2D array of LU-factored form of \(C - t I\), with the non-diagonal part of \(L\) stored in the strictly lower triangle and \(U\) stored in the upper triangle (we skip the permutation matrix, as it won’t impact the 1-norm)

  • the 1-norm the matrix \(C - t I\)

As mentioned above, these two values are meant to be used with dgecon.

Return type:

Tuplenumpy.ndarray, float

bezier.hazmat.algebraic_intersection.bezier_value_check(coeffs, s_val, rhs_val=0.0)

Check if a polynomial in the Bernstein basis evaluates to a value.

This is intended to be used for root checking, i.e. for a polynomial \(f(s)\) and a particular value \(s_{\ast}\):

Is it true that \(f\left(s_{\ast}\right) = 0\)?

Does so by re-stating as a matrix rank problem. As in bezier_roots(), we can rewrite

\[f(s) = (1 - s)^n g\left(\sigma\right)\]

for \(\sigma = \frac{s}{1 - s}\) and \(g\left(\sigma\right)\) written in the power / monomial basis. Now, checking if \(g\left(\sigma\right) = 0\) is a matter of checking that \(\det\left(C_g - \sigma I\right) = 0\) (where \(C_g\) is the companion matrix of \(g\)).

Due to issues of numerical stability, we’d rather ask if \(C_g - \sigma I\) is singular to numerical precision. A typical approach to this using the singular values (assuming \(C_g\) is \(m \times m\)) is that the matrix is singular if

\[\sigma_m < m \varepsilon \sigma_1 \Longleftrightarrow \frac{1}{\kappa_2} < m \varepsilon\]

(where \(\kappa_2\) is the 2-norm condition number of the matrix). Since we also know that \(\kappa_2 < \kappa_1\), a stronger requirement would be

\[\frac{1}{\kappa_1} < \frac{1}{\kappa_2} < m \varepsilon.\]

This is more useful since it is much easier to compute the 1-norm condition number / reciprocal condition number (and methods in LAPACK are provided for doing so).

Parameters:

  • coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis representing a polynomial.

  • s_val (float) – The value to check on the polynomial: \(f(s) = r\).

  • rhs_val (Optionalfloat) – The value to check that the polynomial evaluates to. Defaults to 0.0.

Returns:

Indicates if \(f\left(s_{\ast}\right) = r\) (where \(s_{\ast}\) is s_val and \(r\) is rhs_val).

Return type:

bool

bezier.hazmat.algebraic_intersection.roots_in_unit_interval(coeffs)

Compute roots of a polynomial in the unit interval.

Parameters:

coeffs (numpy.ndarray) – A 1D array (size d + 1) of coefficients in monomial / power basis.

Returns:

N-array of real values in \(\left[0, 1\right]\).

Return type:

numpy.ndarray

bezier.hazmat.algebraic_intersection.intersect_curves(nodes1, nodes2)

Intersect two parametric Bézier curves.

Parameters:
Returns:

2 x N array of intersection parameters. Each row contains a pair of values \(s\) and \(t\) (each in \(\left[0, 1\right]\)) such that the curves intersect: \(B_1(s) = B_2(t)\).

Return type:

numpy.ndarray

Raises:

NotImplementedError – If the “intersection polynomial” is all zeros – which indicates coincident curves.

bezier.hazmat.algebraic_intersection.poly_to_power_basis(bezier_coeffs)

Convert a Bézier curve to polynomial in power basis.

Note

This assumes, but does not verify, that the “Bézier degree” matches the true degree of the curve. Callers can guarantee this by calling full_reduce().

Parameters:

bezier_coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis.

Returns:

1D array of coefficients in monomial basis.

Return type:

numpy.ndarray

Raises:

UnsupportedDegree – If the degree of the curve is not among 0, 1, 2 or 3.

bezier.hazmat.algebraic_intersection.locate_point(nodes, x_val, y_val)

Find the parameter corresponding to a point on a curve.

Note

This assumes that the curve \(B(s, t)\) defined by nodes lives in \(\mathbf{R}^2\).

Parameters:

  • nodes (numpy.ndarray) – The nodes defining a Bézier curve.

  • x_val (float) – The \(x\)-coordinate of the point.

  • y_val (float) – The \(y\)-coordinate of the point.

Returns:

The parameter on the curve (if it exists).

Return type:

Optionalfloat

bezier.hazmat.algebraic_intersection.all_intersections(nodes_first, nodes_second)

Find the points of intersection among a pair of curves.

Note

This assumes both curves are in \(\mathbf{R}^2\), but does not explicitly check this. However, functions used here will fail if that assumption fails.

Parameters:

  • nodes_first (numpy.ndarray) – Control points of a curve to be intersected with nodes_second.

  • nodes_second (numpy.ndarray) – Control points of a curve to be intersected with nodes_first.

Returns:

An array and a flag:

  • A 2 x N array of intersection parameters. Each row contains a pair of values \(s\) and \(t\) (each in \(\left[0, 1\right]\)) such that the curves intersect: \(B_1(s) = B_2(t)\).

  • Flag indicating if the curves are coincident. (For now, this will always be False since intersect_curves() fails if coincident curves are detected.)

Return type:

Tuplenumpy.ndarray, bool