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
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
- Raises
ValueError – If the curve is a point.
UnsupportedDegree – If the degree is not 1, 2 or 3.
-
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 ofnodes1
. Then plugst
into the second parametric curve to get anx
- andy
-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
-
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 bynodes2
.- Parameters
nodes1 (numpy.ndarray) – The nodes in the first curve.
nodes2 (numpy.ndarray) – The nodes in the second curve.
- Returns
Array of coefficients.
- Return type
- Raises
NotImplementedError – If the degree pair is not
1-1
,1-2
,1-3
,1-4
,2-2
,2-3
,2-4
or3-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
-
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 (
Optional
[float
]) – 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
-
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
Tuple
[numpy.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\]>>> 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
-
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
Tuple
[numpy.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 (
Optional
[float
]) – The value to check that the polynomial evaluates to. Defaults to0.0
.
- Returns
Indicates if \(f\left(s_{\ast}\right) = r\) (where \(s_{\ast}\) is
s_val
and \(r\) isrhs_val
).- Return type
-
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
-
bezier.hazmat.algebraic_intersection.
intersect_curves
(nodes1, nodes2)¶ Intersect two parametric Bézier curves.
- Parameters
nodes1 (numpy.ndarray) – The nodes in the first curve.
nodes2 (numpy.ndarray) – The nodes in the second curve.
- 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
- 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
- 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
-
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
sinceintersect_curves()
fails if coincident curves are detected.)
- Return type
Tuple
[numpy.ndarray
,bool
]