On Intersections of B-Spline Curves

Abstract

Bézier and B-spline curves are foundational tools for curve representation in computer graphics and computer-aided geometric design, with their intersection computation presenting a fundamental challenge in geometric modeling. This study introduces an innovative algorithm that quickly and effectively resolves intersections between Bézier and B-spline curves. The number of intersections between the two input curves within a specified region is initially determined by applying the resultant of a polynomial system and Sturm’s theorem. Subsequently, the potential region of the intersection is established through the utilization of the pseudo-curvature-based subdivision scheme and the bounding box detection technique. The projected Gauss-Newton method is ultimately employed to efficiently converge to the intersection. The robustness and efficiency of the proposed algorithm are demonstrated through numerical experiments, demonstrating a speedup of 3 to 150 times over traditional methods.

1. Introduction

Geometric modeling is a computer-based technology for the representation, manipulation, analysis, and design of solids [1]. Curves/surfaces modeling technology, which originated from shaping the outer appearance of products such as automobiles, ship fuselages, aircraft wings, and other products, is a significant research area in computer-aided design (CAD) and computer-aided geometric design (CAGD). Its applications span a broad spectrum of industrial engineering fields, encompassing CNC tooling, alignment of roads, robot path planning, design of luggage shells, and clearance detection for sheet metal parts, among others [2,3,4,5].

Hoffmann pointed out that computing the intersections of parametric curves and surfaces represents a foundational challenge within the field of CAGD and geometric modeling [6]. The resolution of intersections between curves and surfaces is critical in applications such as Boolean operations on solids [7], surface rendering via ray tracing [8,9], and collision detection [10]. In the process of geometric modeling, intersections play a significant role in determining the shape and structure of the model, so it is important to accurately identify and represent these features to ensure the model is correctly interpreted and manufactured [11,12]. As CAD/CAGD/CAM systems evolve and scientific and technological advancements continue, the computational demands and data volumes required to tackle intersection challenges are rapidly escalating. Therefore, the development of efficient and robust methods to address these challenges has become crucial.

Bézier curves defined by Bernstein basis functions and control points have a simple and intuitive mathematical expression and are easy to calculate [1]. B-spline curves, as an extension of Bézier curves, are a curve representation method renowned for their excellent local control capacity and smoothness. The flexibility to adjust the shape of B-spline curves through manipulation of knot vectors and control points has made them popular in computer graphics for curve and surface modeling, animation, interpolation and shape editing, and font design. With their different characteristics and advantages, Bézier curves and B-spline curves become important tools for representing and manipulating curves in computer graphics, and many scholars have conducted research on the intersection problems of Bézier [13,14] and B-spline curves/surfaces [15,16,17] .

In fact, the problem of curve/curve intersection is somewhat equivalent to finding the roots of polynomials, leading some researchers to tackle these issues concurrently. The Bézier clipping algorithm was first proposed by Sederberg to robustly and quickly calculate not only the intersection points but also the tangent points of two Bézier curves [18]. The convergence rate of the Bézier clipping algorithm was proved to be second-order in [19]. Since then, various improved algorithms based on the Bézier clipping algorithm have been developed to address both polynomial root finding and curve intersection challenges. Based on degree reduction, the quadratic and cubic polynomials were generated to enclose the graph of the polynomial within the interval of interest, respectively, by Bartoň and Jüttler [20] and Liu [21] et al. for computing the all roots of a univariate polynomial equation. Additionally, the cubic clipping was proved to have at least a second-order convergence rate and used to compute the intersections of two Bézier curves [22]. A geometric interval algorithm that can tightly bind a curve/surface or contain a point on a curve/surface was proposed by North [23]. Further advancements were made by Yuan [24], who developed the cubic hybrid clipping method with a fourth-order convergence rate for root finding of univariate polynomial equations. This method was subsequently adapted by Wu and Li to address curve intersection problems [25].

This study focuses on resolving the intricate challenge of identifying all intersection points between Bézier curves and B-spline curves. We introduce a robust and efficient approach for determining the intersections for both Bézier curves and B-spline curves. First, a robust algorithm for determining the number of intersections between two Bézier curves is proposed. The core of this problem lies in the task of determining the roots of a system of bivariate polynomial equations. To this end, the resultants of bivariate polynomials and Sturm’s sequence are employed. Second, a fast computation algorithm based on pseudo-curvature and subdivision is proposed in this paper. While Newton’s iteration is a well-known technique for solving non-linear equations, it is known to encounter difficulties, settling on local solutions in scenarios where multiple solutions exist. To address this, we subdivide the input curves by considering the pseudo-curvature of the curves and the intersection of their bounding boxes. Once segments are deemed to be sufficiently straight, we ensure that any two intersect at most once. Eventually, the Gauss-Newton method is employed to quickly converge on the intersection points. Numerical experiments validate the superior performance of our approach compared to traditional algorithms, achieving a speedup of 3 to 150 times.

The remainder of this paper is structured as follows: Section 2 presents essential definitions of Bézier and B-spline curves, along with Sturm’s theorem. In Section 3, we outline a method for computing the number of intersections between polynomial parametric curves, focusing primarily on Bézier curves within the unit square. Section 4 introduces our proposed subdivision-based method tailored for efficiently computing all intersections between B-spline curves. Section 5 presents numerical experiments assessing the performance of the proposed method. Finally, Section 6 offers concluding remarks summarizing our contributions and suggesting potential avenues for future research.

2. Preliminary

In this section, we commence with a brief overview of the concepts underlying Bézier and B-spline curves. Subsequently, we introduce Sturm’s theorem, an essential tool for isolating the real roots of higher-order polynomial equations.

2.1. Bézier and B-Spline Curves

A degree-n Bézier curve is expressed as

$$\mathcal{C}\left(\xi \right)=\sum _{i=0}^n{\mathit{P}}_i{B}_i^n\left(\xi \right),\xi \in [0,1],$$

(1)

where ${\mathit{P}}_i\in {\mathbb{R}}^d$ ($d=2,3$ typically) denote the control points, and $B_i^n\left(\xi \right)=\left(\genfrac{}{}{0pt}{}{n}{i}\right){\xi }^i{(1-\xi )}^{n-i}$$(i=0,1,\cdots ,n)$ are Bernstein basis functions [1]. These curves are characterized by their high degree of smoothness, ease of computation, and the property that they are contained within the convex hull of their control points. Figure 1 illustrates quartic Bernstein basis functions and its associated Bézier curve.

Building on the concept of Bézier curves, B-spline curves offer increased flexibility and local modification capability through the introduction of a knot vector. The knot vector, a non-decreasing sequence denoted as:

$$\Xi =\{{\xi }_0,{\xi }_1,\dots ,{\xi }_{n+p+1}\},$$

(2)

where ${\xi }_i\le {\xi }_{i+1}$ for $i=1,2,\dots ,n+p$ serves as the basis for defining B-spline curves [1].

The B-spline basis functions defined over $\Xi$ can be derived recursively by the following de Boor-Cox formula

$$\left\{\begin{array}{c}{N}_{i,0}\left(\xi \right)=\left\{\begin{array}{cc}1,\hfill & \xi \in [{\xi }_i,{\xi }_{i+1}),\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \\ N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}\left(\xi \right),p\ge 1.\hfill \end{array}\right.$$

(3)

Consequently, a B-spline curve of degree p is formulated as the linear combination of its basis functions, $N_{i,p}\left(\xi \right)$, for $i=1,2,\dots ,n$, and the corresponding control points, ${\mathit{P}}_i$, where ${\mathit{P}}_i\in {\mathbb{R}}^d$ with $d=2,3$. Thus, a B-spline curve can be expressed as:

$$\mathcal{C}\left(\xi \right)=\sum _{i=1}^n{\mathit{P}}_i{N}_{i,p}\left(\xi \right),\xi \in [{\xi }_p,{\xi }_{n+1}].$$

(4)

Figure 2 shows quadratic B-spline basis functions defined over the knot vector $\Xi =\{0,0,0,0.25,0.5,0.75,0.75,1,1,1\}$ and its associated B-spline curve.

B-spline basis functions are endowed with essential desirable properties such as non-negativity and partition of unity, which contribute significantly to the stability and manipulability of the curves they define. Therefore, B-spline curves, closely related to these basis functions, exhibit the so-called convex-hull property, ensuring they remain within the convex hull defined by their control points. These features contribute to their widespread application in CAGD and CG for their stability and ease of manipulation. Furthermore, a degree-p Bézier curve can be considered a specific form of a B-spline curve, characterized by a particular knot vector $\Xi =\{\underset{\underset{p+1}{︸}}{0,0,\dots ,0},\underset{\underset{p+1}{︸}}{1,1,\dots ,1}\}$. This connection means that Bézier curves inherit many of the advantageous properties of B-spline curves, enhancing their utility and desirability in various applications.

2.2. Resultant and Sturm’s Theorem

To lay a solid theoretical foundation for the ensuing discussion, this subsection introduces the basic concept of the resultant of two polynomials, the Sturm sequence associated with a polynomial, and Sturm’s theorem.

Definition 1

([26]). Let $f_1\left(x\right)={\sum }_{i=0}^m{a}_i{x}^i$ and $f_2\left(x\right)={\sum }_{i=0}^n{b}_i{x}^i$ represent two polynomials of degree m and n, respectively. The resultant of $f_1$ and $f_2$, $Res(f_1,f_2)$, is defined by a determinant that encapsulates the coefficients of these polynomials, i.e.,

$$Res(f_1,f_2)=\left|\begin{array}{cccccccc}\hfill a_m& a_{m-1}& \hfill \cdots \hfill & a_1& a_0& 0& \cdots & 0\\ \hfill 0& a_m& \hfill a_{m-1}\hfill & \cdots & a_1& a_0& \cdots & 0\\ \hfill ⋮& \ddots & \hfill \ddots \hfill & \ddots & \ddots & \ddots & \ddots & ⋮\\ \hfill 0& \cdots & \hfill a_m\hfill & a_{m-1}& \cdots & \cdots & a_1& a_0\\ \hfill b_n& b_{n-1}& \hfill \cdots \hfill & b_1& b_0& 0& \cdots & 0\\ \hfill 0& b_n& \hfill b_{n-1}\hfill & \cdots & b_1& b_0& \cdots & 0\\ \hfill ⋮& \ddots & \hfill \ddots \hfill & \ddots & \ddots & \ddots & \ddots & ⋮\\ \hfill 0& \cdots & \hfill b_n\hfill & b_{n-1}& \cdots & \cdots & b_1& b_0\end{array}\right|\begin{array}{c}\mathit{row}1\hfill \\ \mathit{row}2\hfill \\ ⋮\hfill \\ \mathit{row}n\hfill \\ \mathit{row}n+1\hfill \\ \mathit{row}n+2\hfill \\ ⋮\hfill \\ \mathit{row}n+m\hfill \end{array}.$$

(5)

On a natural extension of Definition 1, when considering the resultant of two bivariate polynomials concerning one of the variables, the other variable is treated as a constant. A simple example is given to illustrate below.

Example 1.

Given two bivariate polynomials $f_1(x,y)=5x^2-6xy+5y^2-16$ and $f_2(x,y)=2x^2-(1+y)x+y^2-y-4$. The resultant defined by $f_1,f_2$ with respect to variable x is

$$Res(f_1,f_2;x)=\left|\begin{array}{cccc}5& -6y& 5y^2-16& 0\\ 0& 5& -6y& 5y^2-16\\ 2& -(1+y)& y^2-y-4& 0\\ 0& 2& -(1+y)& y^2-y-4\end{array}\right|,$$

which is essentially a polynomial of the variable y.

Consider a polynomial $f\left(x\right)$ with real coefficients and free of any square factors. The Sturm sequence for $f\left(x\right)$ is constructed as follows

$$\begin{array}{c}{\mathbf{Sturm}}^{\left(0\right)}=f\left(x\right)\hfill \\ {\mathbf{Sturm}}^{\left(1\right)}=f^{\prime }\left(x\right)\hfill \\ {\mathbf{Sturm}}^{\left(\mathrm{i}\right)}=-\mathbf{rem}\left({\mathbf{Sturm}}^{(\mathrm{i}-2)},{\mathbf{Sturm}}^{(\mathrm{i}-1)}\right),\mathrm{for}\mathrm{i}=2,3,\dots ,\mathrm{r},\hfill \end{array}$$

(6)

where $\mathbf{rem}$ denotes the polynomial remainder operator. In this sequence, each successive term is the negated remainder from dividing the two preceding terms. The construction of this sequence concludes when a constant term or zero is reached at ${\mathbf{Sturm}}^{\left(r\right)}$. The Sturm sequence, denoted as $\mathbf{Sturm}=\{{\mathbf{Sturm}}^{\left(0\right)},{\mathbf{Sturm}}^{\left(1\right)},\dots ,{\mathbf{Sturm}}^{\left(r\right)}\}$, serves as an effective tool for identifying the count of zeros that the polynomial $f\left(x\right)$ possesses within a given interval.

Consider the sequence $\mathbf{\Gamma }=\{{\gamma }_0,{\gamma }_1,\dots ,{\gamma }_t\}$. Denote by $pre\left(\mathbf{\Gamma }\right)$ the number of sign-preserving occurrences within $\mathbf{\Gamma }$. This is quantified by examining each pair of consecutive values, ${\gamma }_i$ and ${\gamma }_{i+1}$. If these values share the same sign, the pair $({\gamma }_i,{\gamma }_{i+1})$ is counted as contributing 1 to $pre\left(\mathbf{\Gamma }\right)$. Notably, a term ${\gamma }_i=0$ is treated as though it possesses the opposite sign of its immediate predecessor ${\gamma }_{i-1}$ [27,28].

Example 2.

Suppose ${\mathbf{\Gamma }}_1=[1,2,1,2,-1,-2]$ and ${\mathbf{\Gamma }}_2=[1,0,1,2,-1,-2]$. Then, we have the numbers of sign-preserving occurrences $pre\left({\mathbf{\Gamma }}_1\right)=4$ and $pre\left({\mathbf{\Gamma }}_2\right)=2$.

Integrating the definition of the resultant for two bivariate polynomials with the concept of sign-preserving within a Sturm sequence of a specified polynomial, Sturm’s theorem emerges as a powerful method to ascertain the count of real roots for a bivariate polynomial system within a defined interval.

Theorem 1

([29]). Consider a compact set $\mathbb{D}=[x_{min},x_{max}]\times [y_{min},y_{max}]\subset {\mathbb{R}}^2$. Let $f_1(x,y)$ and $f_2(x,y)$ be two bivariate polynomials with the degrees of x as $n_1$ and $n_2$, respectively. Define a polynomial in terms of α by

$$RES(\alpha ,{\alpha }_1,{\alpha }_2)=\frac{Res(Res(f_1,g;x),Res(f_2,g;x);y)}{{\alpha }^n},$$

(7)

where $g(x,y)=\alpha +(x-{\alpha }_1)(y-{\alpha }_2)$ and $n=n_1\times n_2$.

Let $\left\{\mathbf{Sturm}(\alpha ,{\alpha }_1,{\alpha }_2)\right\}$ stand for the Sturm sequence of the polynomial $RES(\alpha ,{\alpha }_1,{\alpha }_2)$ of α. Here, $Num\_pre(a,b)$ represents the count of sign-preserving occurrences in the Sturm sequence $\left\{\mathbf{Sturm}(\alpha ,{\alpha }_1,{\alpha }_2)\right\}$ when $\alpha =0$, ${\alpha }_1=a$, and ${\alpha }_2=b$. Consequently, the total number of real roots for the system of equations $f_1=f_2=0$ within the compact set $\mathbb{D}$ can be determined by

$$\begin{array}{cc}\hfill Num(f_1,f_2,\mathbb{D})=& \frac{1}{2}[Num\_pre(x_{min},y_{min})-Num\_pre(x_{min},y_{max})\hfill \\ & -Num\_pre(x_{max},y_{min})+Num\_pre(x_{max},y_{max})].\hfill \end{array}$$

(8)

Remark 1.

In cases where the first term in a sequence is 0, it is stated that $pre\left(\right[0,1\left]\right)=pre\left(\right[0,-1\left]\right)=1/2$, as detailed in [29]. This convention implies that roots located on the boundaries of a specified compact set $\mathbb{D}$ are represented by $1/2$ in Equation (8). Similarly, roots situated at the vertices of the compact set $\mathbb{D}$ are attributed a value of $1/4$.

Theorem 1 provides a theoretical foundation for determining the exact number of roots for a pair of bivariate polynomial equations $f_1(x,y)=f_2(x,y)=0$ within a specified compact set $\mathbb{D}$. In the subsequent discussion, we will apply Theorem 1 to ascertain the number of intersections between two polynomial parametric curves. This step can be regarded as a preparatory phase, essential for enhancing the precision and reliability of the subsequent calculation of intersections.

3. Determining the Number of Intersections

In this section, we provide a method for counting the number of intersections between two polynomial parametric curves, with an emphasis on Bézier curves within the interval $[0,1]\times [0,1]$ for clarity and brevity.

Given two Bézier curves

$${\mathcal{C}}_A\left(\xi \right)=\sum _{i=0}^n{\mathit{P}}_i{B}_i^n\left(\xi \right)$$

and

$${\mathcal{C}}_B\left(\eta \right)=\sum _{j=0}^m{\mathit{Q}}_j{B}_j^m\left(\eta \right),$$

where ${\mathit{P}}_i={(x_i,y_i,z_i)}^{\mathrm{T}}$, ${\mathit{Q}}_j={({\tilde{x}}_j,{\tilde{y}}_j,{\tilde{z}}_j)}^{\mathrm{T}}$ are their control points. Computing the intersection of ${\mathcal{C}}_A\left(\xi \right)$ and ${\mathcal{C}}_B\left(\eta \right)$ involves identifying the parameter pair ($\xi$, $\eta$) that fulfills the condition ${\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)=\mathbf{0}$.

From the partition of the unity property of Bernstein basis functions, we have

$$\begin{array}{cc}{\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)\hfill & =\sum _{i=0}^n{\mathit{P}}_i{B}_i^n\left(\xi \right)-\sum _{j=0}^m{\mathit{Q}}_j{B}_j^m\left(\eta \right)\hfill \\ & =\sum _{j=0}^m{B}_j^m\left(\eta \right)·\sum _{i=0}^n{\mathit{P}}_i{B}_i^n\left(\xi \right)-\sum _{i=0}^n{B}_i^n\left(\xi \right)·\sum _{j=0}^m{\mathit{Q}}_j{B}_j^m\left(\eta \right)\hfill \\ & =\sum _{i=0}^n\sum _{j=0}^m({\mathit{P}}_i-{\mathit{Q}}_j)B_i^n\left(\xi \right)B_j^m\left(\eta \right).\hfill \end{array}$$

(9)

To simplify the notation and facilitate computation, let us denote

$$\begin{array}{cc}\hfill {\mathbf{B}}^n\left(\xi \right)& \hfill ={\left(B_0^n\left(\xi \right),B_1^n\left(\xi \right),\dots ,B_n^n\left(\xi \right)\right)}^{\mathrm{T}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathbf{B}}^m\left(\eta \right)& \hfill ={\left(B_0^m\left(\eta \right),B_1^m\left(\eta \right),\dots ,B_m^m\left(\eta \right)\right)}^{\mathrm{T}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\Delta }_i& \hfill ={\left({\mathit{P}}_i-{\mathit{Q}}_0,{\mathit{P}}_i-{\mathit{Q}}_1,\dots ,{\mathit{P}}_i-{\mathit{Q}}_m\right)}^{\mathrm{T}},\hfill \end{array}$$

(12)

and

$$\mathbf{R}={\left({\mathbf{R}}_0,{\mathbf{R}}_1,\dots ,{\mathbf{R}}_n\right)}^{\mathrm{T}},$$

(13)

where ${\mathbf{R}}_i={\sum }_{j=0}^m({\mathit{P}}_i-{\mathit{Q}}_j)B_j^m\left(\eta \right)={\Delta }_i^{\mathrm{T}}·{\mathbf{B}}^m\left(\eta \right)$.

Therefore, Equation (9) can be expressed in a matrix form as

$${\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)={\mathbf{R}}^{\mathrm{T}}·{\mathbf{B}}^n\left(\xi \right).$$

(14)

Subsequently, according to the relationship between the Bernstein basis functions and power basis functions, we have

$${\mathbf{B}}^n\left(\xi \right)=\mathit{M}·\mathit{\xi },$$

where $\mathit{M}$ is the basis transformation matrix and $\mathit{\xi }={(1,\xi ,{\xi }^2,\dots ,{\xi }^n)}^{\mathrm{T}}$ is the collection of power basis functions. Therefore, we have

$${\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)=\left(\begin{array}{c}x(\xi ,\eta )\\ y(\xi ,\eta )\\ z(\xi ,\eta )\end{array}\right)={\mathbf{R}}^{\mathrm{T}}·\mathit{M}·\mathit{\xi }.$$

(15)

Recall that computing the intersections between ${\mathcal{C}}_A\left(\xi \right)$ and ${\mathcal{C}}_B\left(\eta \right)$ is equivalent to solving the roots of ${\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)=\mathbf{0}$. That is,

$$\left\{\begin{array}{c}x(\xi ,\eta )=0\\ y(\xi ,\eta )=0\\ z(\xi ,\eta )=0\end{array}\right..$$

(16)

Given that Theorem 1 is applicable solely for isolating the real roots of bivariate polynomials, Equation (16) can be equivalently converted into the following form

$$\left\{\begin{array}{c}{f}_1(\xi ,\eta )=x^2(\xi ,\eta )+y^2(\xi ,\eta )=0\\ f_2(\xi ,\eta )=y^2(\xi ,\eta )+z^2(\xi ,\eta )=0\end{array}\right..$$

(17)

Denote by ${\Delta }_i^1$ and ${\Delta }_i^2$ the first two lines and the last two lines of ${\Delta }_i$, respectively. Similarly, we have

$${\mathbf{R}}^{\phi }={\left({\mathbf{R}}_0^{\phi },{\mathbf{R}}_1^{\phi },\dots ,{\mathbf{R}}_n^{\phi }\right)}^{\mathrm{T}},$$

(18)

where ${\mathbf{R}}_i^{\phi }={\left({\Delta }_i^{\phi }\right)}^{\mathrm{T}}·{\mathbf{B}}^m\left(\eta \right)$ with $\phi =1,2$. Accordingly, the Equation (17) to compute the intersection of ${\mathcal{C}}_A\left(\xi \right)$ and ${\mathcal{C}}_B\left(\eta \right)$ can thus be reformulated as follows

$$\left\{\begin{array}{c}{f}_1(\xi ,\eta )={\left({\mathbf{R}}^1\mathit{M}\mathit{\xi }\right)}^{\mathrm{T}}\ast \left({\mathbf{R}}^1\mathit{M}\mathit{\xi }\right)={\mathit{\xi }}^{\mathrm{T}}\left({\mathit{M}}^{\mathrm{T}}{\left({\mathbf{R}}^1\right)}^{\mathrm{T}}{\mathbf{R}}^1\mathit{M}\right)\mathit{\xi },\\ f_2(\xi ,\eta )={\left({\mathbf{R}}^2\mathit{M}\mathit{\xi }\right)}^{\mathrm{T}}\ast \left({\mathbf{R}}^2\mathit{M}\mathit{\xi }\right)={\mathit{\xi }}^{\mathrm{T}}\left({\mathit{M}}^{\mathrm{T}}{\left({\mathbf{R}}^2\right)}^{\mathrm{T}}{\mathbf{R}}^2\mathit{M}\right)\mathit{\xi }.\end{array}\right.$$

(19)

Remark 2.

Referring to Equation (7) from Theorem 1, the computation of the polynomial $RES(\alpha ,{\alpha }_1,{\alpha }_2)$ necessitates initially deriving the resultant of the polynomials $f_i$ ($i=1,2$) and g with respect to the variable ξ, treating η as a constant during this process. Consequently, we reformulate Equation (17) into a quadratic equation in terms of ξ.

Remark 3.

In the case of plane polynomial parametric curves, the approach simplifies significantly. We can directly set $f_1(\xi ,\eta )=x(\xi ,\eta )$ and $f_2(\xi ,\eta )=y(\xi ,\eta )$, thereby streamlining the process.

For clearer understanding, we present an algorithm that determines the number of intersections between two Bézier curves based on Theorem 1, followed by an illustrative example in Example 3. Moreover, the proposed algorithm may be naturally applicable to B-spline curves by combining with the Bézier extraction technique (Algorithm 1).

Algorithm 1: Number of intersections between two Bézier curves

Example 3.

Let ${\mathit{P}}_0=(0,0)$, ${\mathit{P}}_1=(1,1)$, ${\mathit{P}}_2=(2,0)$ and ${\mathit{Q}}_0=(0.5,0)$, ${\mathit{Q}}_1=(1,1)$, ${\mathit{Q}}_2=(0,2)$ be the control points of two quadratic Bézier curves ${\mathcal{C}}_A\left(\xi \right)$ and ${\mathcal{C}}_B\left(\eta \right)$, as shown in Figure 3, respectively. Thus computing the intersections between ${\mathcal{C}}_A\left(\xi \right)$ and ${\mathcal{C}}_B\left(\eta \right)$ is equivalent to the roots of the following bivariate polynomials

$$\left\{\begin{array}{c}{f}_1(\xi ,\eta )=2\xi +\frac{3}{2}{\eta }^2-\eta -\frac{1}{2}\hfill \\ f_2(\xi ,\eta )=-2{\xi }^2+2\xi -2\eta \hfill \end{array}\right..$$

Now, we propose the detailed procedure for determining the roots of $f_1(\xi ,\eta )=f_2(\xi ,\eta )=0$ in $[0,1]\times [0,1]$ by Theorem 1.

First, compute $Num\_pre(0,0)$, $Num\_pre(0,1)$, $Num\_pre(1,0)$, and $Num\_pre(1,1)$ in Theorem 1, respectively. A comprehensive description of the process for calculating $Num\_pre(0,0)$ is shown below, which can also be applied to appraising $Num\_pre(0,1)$, $Num\_pre(1,0)$, and $Num\_pre(1,1)$.

1. 

Let ${\alpha }_1=0,{\alpha }_2=0$, then $g(\xi ,\eta )=\alpha +\xi \eta =\eta \xi +\alpha$. Similar to Example 1, we can state that the resultant $Res(f_1,g;\xi )$ and $Res(f_2,g;\xi )$ in Equation (7) as

$$Res(f_1,g;\xi )=\left|\begin{array}{cc}2& \frac{3}{2}{\eta }^2-\eta -\frac{1}{2}\\ \eta & \alpha \end{array}\right|=-\frac{3}{2}{\eta }^3+{\eta }^2+\frac{1}{2}\eta +2\alpha ,$$

(20)

and

$$Res(f_2,g;\xi )=\left|\begin{array}{ccc}-2& 2& -2\eta \\ \eta & \alpha & 0\\ 0& \eta & \alpha \end{array}\right|=-2{\eta }^3-2\alpha \eta -2{\alpha }^2.$$

(21)

The outer resultant in Equation (7) is computed by

$$\begin{array}{cc}\hfill Res(Res(f_1,g;\xi ),& Res(f_2,g;\xi );\eta )\hfill \\ & =\left|\begin{array}{cccccc}3/2& 1& 1/2& 2\alpha & 0& 0\\ 0& 3/2& 1& 1/2& 2\alpha & 0\\ 0& 0& 3/2& 1& 1/2& 2\alpha \\ -2& 0& -2\alpha & -2{\alpha }^2& 0& 0\\ 0& -2& 0& -2\alpha & -2{\alpha }^2& 0\\ 0& 0& -2& 0& -2\alpha & -2{\alpha }^2\end{array}\right|\hfill \\ & =27{\alpha }^6+126{\alpha }^5+173{\alpha }^4+54{\alpha }^3+3{\alpha }^2,\hfill \end{array}$$

(22)

and

$$\begin{array}{cc}\hfill RES(\alpha ,0,0)& =\frac{Res(Res(f_1,g;\xi ),Res(f_2,g;\xi );\eta )}{{\alpha }^{1\times 2}}\hfill \\ & =27{\alpha }^4+126{\alpha }^3+173{\alpha }^2+54\alpha +3.\hfill \end{array}$$

(23)

2. 

The Sturm sequence $\left\{\mathbf{Sturm}\right(\alpha ,0,0\left)\right\}$ of $RES(\alpha ,0,0)$ can be obtained by the method of successive division. That is

$$\begin{array}{c}\hfill {\mathbf{Sturm}}^{\left(0\right)}=-0.017341-0.312139\alpha -{\alpha }^2-0.728324{\alpha }^3-0.156069{\alpha }^4,\\ {\mathbf{Sturm}}^{\left(1\right)}=-0.142857-0.915344\alpha -{\alpha }^2-0.285714{\alpha }^3,\hfill \\ {\mathbf{Sturm}}^{\left(2\right)}=-0.211034-\alpha -0.393103{\alpha }^2,\hfill \\ {\mathbf{Sturm}}^{\left(3\right)}=-0.056683+\alpha ,\hfill \\ {\mathbf{Sturm}}^{\left(4\right)}=1.\hfill \end{array}$$

(24)

3. 

By substituting $\alpha =0$ into Equation (24), we obtain the sequence $\mathbf{\Gamma }=[-0.017341,-0.142857,$$-0.211034,-0.056683,1]$. It is observed that $Num\_pre(0,0)=3$, corresponding to the sign-preserving of Γ.

Analogously, we find that $Num\_pre(0,1)=2$, $Num\_pre(1,0)=2$, and $Num\_pre(1,1)=3$, respectively. Consequently, according to Equation (8), it is concluded that the Bézier curves $\mathit{P}\left(\xi \right)$ and $\mathit{Q}\left(\eta \right)$ intersect exactly once.

4. A Fast Subdivision-Based Intersection Algorithm

In the realm of geometric modeling and computer graphics, accurately determining the intersections between two parametric curves is a fundamental challenge with numerous applications. This task is particularly demanding due to the need for computational efficiency and high precision. In this section, we propose a novel subdivision-based approach specifically designed to compute all intersections between B-spline curves efficiently.

4.1. Pseudo-Curvature of B-Spline Curves

Central to our intersection algorithm lies an effective subdivision scheme, capitalizing on the curvature characteristics of B-spline curves for more efficient computation. This scheme introduces an effective curvature-based subdivision strategy, boosting computational efficiency and accuracy when identifying intersection points.

Generally, the point-wise curvature at any given parameter $\xi$ along a parametric curve is determined by

$${\kappa }_{\mathrm{p}}\left(\mathcal{C}\left(\xi \right)\right)=\frac{\left|{\mathcal{C}}^{\prime }\left(\xi \right)\times {\mathcal{C}}^{″}\left(\xi \right)\right|}{{\left|{\mathcal{C}}^{\prime }\left(\xi \right)\right|}^3}.$$

(25)

To ascertain the total curvature energy across the curve, one integrates the point-wise curvature as follows:

$${\kappa }_{\mathrm{tradition}}=\int {\kappa }_{\mathrm{p}}\left(\mathcal{C}\left(\xi \right)\right)\mathrm{d}\xi .$$

(26)

In practice, the above integration is typically executed using numerical quadrature techniques, such as Simpson’s rule or Gauss-Legendre quadrature.

Although the above traditional methods for calculating the curvature of parametric curves are accurate, they involve intricate calculations that can be computationally expensive. To address this, we introduce the concept of pseudo-curvature, with a particular focus on B-spline curves. This approach seeks to offer a less computationally intensive alternative while maintaining sufficient accuracy for intersection determination tasks.

For a B-spline curve represented as $\mathcal{C}\left(\xi \right)={\sum }_{i=0}^n{\mathbf{P}}_i{N}_i\left(\xi \right)$, we define its pseudo-curvature as:

$$\kappa \left(\mathcal{C}\right)=\sum _{i=0}^{n-1}\frac{{∥\overline{{\mathbf{P}}_i{\mathbf{P}}_{i+1}}∥}_2}{{∥\overline{{\mathbf{P}}_0{\mathbf{P}}_n}∥}_2},$$

(27)

where $∥·∥$ stands for the Euclidean length of the vector. This formula effectively approximates the bending of the curve, as shown in Figure 4, by considering the relative lengths of segments formed by consecutive control points $\overline{{\mathbf{P}}_i{\mathbf{P}}_{i+1}}$ ($i=0,1,\dots ,n-1$) in comparison to the segment formed by the first and last control points $\overline{{\mathbf{P}}_0{\mathbf{P}}_n}$.

The principle underlying the pseudo-curvature defined in Equation (27) capitalizes on the convex hull property inherent to B-spline curves, which ensures that the B-spline curve lies entirely within the convex polygon formed by its control points. Leveraging the convex hull property, we ascertain that the bending of the B-spline curve approaches its minimum when the value of the pseudo-curvature in Equation (27) approaches 1.

The pseudo-curvature offers a significant reduction in computational complexity compared to traditional curvature metrics. By leveraging the geometric properties inherent in the control points of B-spline curves, we can obtain a useful approximation of the original curve with far less computational effort.

4.2. Subdivision-Based Potential Intersection Ranges Computation

Central to our methodology is a sophisticated subdivision scheme tailored to the unique characteristics of B-spline curves, which significantly enhances the efficiency and accuracy of intersection point identification. With the aforementioned pseudo-curvature concept established in Section 4.1, we have the essential tools for splitting curves effectively. In this section, we integrate subdivision techniques with bounding box intersection detection to compute potential intersection ranges between B-spline curves efficiently. The algorithmic framework for this process is presented in Algorithm 2.

The algorithm depicted in Algorithm 2 begins by computing the bounding boxes for the input curves. It then proceeds to check for bounding box intersections. If no intersection is detected, the algorithm concludes that the curves themselves cannot intersect, reducing unnecessary computation.

In cases where bounding boxes intersect, the algorithm evaluates the pseudo-curvature of the curves in Equation (27). If the curvature falls below a predefined threshold ${\kappa }_{\max }$, indicative of minimal bending, direct intersection checks between curve segments are performed. Conversely, if the curvature exceeds this threshold, suggesting potential complex intersections, the algorithm recursively applies itself to subdivided curve segments, iteratively refining the search space.

This iterative process continues until all potential intersection ranges have been exhaustively examined, ensuring a comprehensive identification of intersection points without omission.

In Algorithm 2, the selection of ${\kappa }_{\max }$ is critical. Setting it too low may result in excessive subdivision, thereby prolonging computation time. Conversely, setting it too high may risk overlooking complex intersections. Striking the right balance for this parameter is essential to achieving optimal performance. In our numerical experiments, we set the default value of ${\kappa }_{\max }$ to $1+5\times {10}^{-6}$.

While bounding boxes provide an initial filter for potential intersections, the precision of intersection detection relies on the subsequent subdivision and direct intersection checks. Maintaining geometric precision throughout the process is crucial. To enhance computational accuracy and efficiency, we employ the projected Gauss-Newton method in the following section.

Algorithm 2: Subdivision-based potential intersection ranges computation

4.3. Projected Gauss-Newton Method

Upon obtaining the output from Algorithm 2, we acquire pairs of potential intersection ranges where the curve segments exhibit very limited bending. According to Bézout’s theorem, as these segments closely resemble straight lines (degree-1 curves), only one intersection point is anticipated within each pair. Hence, leveraging Gauss-Newton iteration proves to be an efficient approach for intersection computation. However, it is essential to note that the computed results may extend beyond the parameter ranges of the current segments. To mitigate this issue, we implement the projected Gauss-Newton method, thereby ensuring that the computed intersections remain within the designated segments.

Consider a nonlinear system $\mathcal{F}\left(\mathit{\xi }\right)=\mathbf{0}$. The essence of the projected Gauss-Newton method is captured by a two-fold iterative process. Initially, we calculate an intermediate point ${\tilde{\mathit{\xi }}}^{k+1}$ using the classic Gauss-Newton iteration formula

$${\tilde{\mathit{\xi }}}^{k+1}={\mathit{\xi }}^k-{\mathit{J}}_k^{†}\mathit{F}\left({\mathit{\xi }}^k\right),$$

(28)

where ${\mathit{J}}_k^{†}={\left(A^{\mathrm{T}}A\right)}^{-1}{A}^{\mathrm{T}}$ is the Moore-Penrose inverse of the Jacobian matrix of the nonlinear system $\mathcal{F}$ at iteration k. The Moore-Penrose inverse ${\mathit{J}}_k^{†}$ can be computed using the singular value decomposition (SVD).

Subsequently, we refine our estimate by projecting ${\tilde{\mathit{\xi }}}^{k+1}$ onto the feasible set ${\mathbb{D}}_k=[{\xi }_{\min }^k,{\xi }_{\max }^k]\times [{\eta }_{\min }^k,{\eta }_{\max }^k]$, optimizing for proximity to the initial estimate

$${\mathit{\xi }}^{k+1}={\mathbb{P}}_{\mathit{\xi }\in {\mathbb{D}}_k}\left({\tilde{\mathit{\xi }}}^{k+1}\right),$$

(29)

where ${\mathbb{P}}_{\mathit{\xi }\in {\mathbb{D}}_k}$ denotes the projection operator. In our implementation, we adopt the following projection operator

$${\mathbb{P}}_{\mathit{\xi }\in {\mathbb{D}}_k}\left(\tilde{\mathit{\xi }}\right)=\left(min\left(max(\tilde{\xi },{\xi }_{\min }^k),{\xi }_{\max }^k\right),min\left(max\left(\tilde{\eta },{\eta }_{\min }^k\right),{\eta }_{\max }^k\right)\right).$$

(30)

This projection ensures our solution remains within the parameter ranges of the curves, an essential consideration for maintaining geometric validity.

In the case of 3D scenarios, our approach to the curve intersection problem involves addressing the following nonlinear system

$$\mathcal{F}(\xi ,\eta )={\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)=\left(\begin{array}{c}{x}_A\left(\xi \right)-x_B\left(\eta \right)\\ y_A\left(\xi \right)-y_B\left(\eta \right)\\ z_A\left(\xi \right)-z_B\left(\eta \right)\end{array}\right)=\mathbf{0}.$$

(31)

The corresponding Jacobian matrix, a $3\times 2$ matrix, can be computed as follows:

$$\mathcal{J}(\xi ,\eta )=\left[\begin{array}{cc}\frac{\mathrm{d}{\mathcal{C}}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{\mathcal{C}}_B\left(\eta \right)}{\mathrm{d}\eta }\end{array}\right]=\left[\begin{array}{cc}\frac{\mathrm{d}{x}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{x}_B\left(\eta \right)}{\mathrm{d}\eta }\\ \frac{\mathrm{d}{y}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{y}_B\left(\eta \right)}{\mathrm{d}\eta }\\ \frac{\mathrm{d}{z}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{z}_B\left(\eta \right)}{\mathrm{d}\eta }\end{array}\right].$$

(32)

For the 2D case, the nonlinear system becomes

$$\mathcal{F}(\xi ,\eta )={\mathcal{C}}_A\left(\xi \right)-{\mathcal{C}}_B\left(\eta \right)=\left(\begin{array}{c}{x}_A\left(\xi \right)-x_B\left(\eta \right)\\ y_A\left(\xi \right)-y_B\left(\eta \right)\end{array}\right)=\mathbf{0}.$$

(33)

The corresponding Jacobian matrix, a $2\times 2$ matrix, is computed by

$$\mathcal{J}(\xi ,\eta )=\left[\begin{array}{cc}\frac{\mathrm{d}{\mathcal{C}}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{\mathcal{C}}_B\left(\eta \right)}{\mathrm{d}\eta }\end{array}\right]=\left[\begin{array}{cc}\frac{\mathrm{d}{x}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{x}_B\left(\eta \right)}{\mathrm{d}\eta }\\ \frac{\mathrm{d}{y}_A\left(\xi \right)}{\mathrm{d}\xi }& -\frac{\mathrm{d}{y}_B\left(\eta \right)}{\mathrm{d}\eta }\end{array}\right].$$

(34)

Different from the 3D case, the Moore-Penrose inverse in Equation (28) degenerates into the inverse of the Jacobian matrix (34).

The comprehensive algorithm is presented in Algorithm 3.

Algorithm 3: Projected Gauss-Newton method

5. Numerical Experiments

This section presents numerical experiments conducted to evaluate the effectiveness and efficiency of the proposed algorithm across both Bézier curves and B-spline curves. We implemented the algorithm in C++ using Clang-15 compiler and conducted the experiments on a MacBook Pro (14-inch, 2021) equipped with an Apple M1 Pro CPU and 16 GB of RAM. Matlab^® 2023a and ParaView (version 5.10.1) were employed for visualization.

In all experiments, the parameter settings were standardized as follows: maximum curvature, ${\kappa }_{\max }$, was set to $1+5\times {10}^{-6}$; residual tolerance, tol_res, was defined at $1\times {10}^{-5}$; variation tolerance for $\mathit{\xi }$, denoted as ${\mathrm{tol}}_{\xi }$, was also set at $1\times {10}^{-5}$; and the limit for the maximum number of iterations, $maxIter$, was fixed at 20. Computing all the intersections between two Bézier curves is a topic that has received considerable attention in the CAGD community. The Bézier clipping algorithm, as proposed by Sederberg in 1990 [18], stands out for its second-order convergence rate, a finding further supported by Schulz in 2009 [19]. Since then, the algorithm has seen several enhancements aimed at improving its efficiency. Noteworthy among these are the quadratic and cubic clipping methods [20,22], which are based on degree reduction techniques; the geometry interval clipping algorithm [23]; and the cubic hybrid clipping method [25].

For B-spline curves, a prevalent strategy is to first convert them into Bézier segments using the Bézier extraction technique. This transformation allows the application of the previously mentioned clipping methods to each pair of Bézier segments, facilitating the computation of intersections. In this section, we compare our proposed method with these existing intersections computation techniques through detailed analysis. Specifically, different methods are conducted on four examples depicted in Figure 5, which include cubic Bézier curves with two intersections, a quadratic B-spline and a cubic B-spline each with two intersections, and two cubic B-splines with eight intersections, including a scenario where two intersection are notably proximal.

Table 1. Computation results comparison of the proposed method against existing methods. The curves information includes the numbers of control points, Bézier segments, and intersections for clearer interpretation. Data highlighted in bold indicate the best performance.

Example	Curves Info.	Methods	Time ($\mathsf{\mu }$s)	Relative Time
		Bézier Clipping [18]	210	1.84
		Quadratic Degree Reduction [20]	3293	28.89
Example 1	${\mathcal{C}}_A$: 4/1/2	Cubic Degree Reduction [22]	17,074	149.77
(Figure 5a)	${\mathcal{C}}_B$: 4/1/2	Geometry Interval Clipping [23]	2941	25.79
		Cubic Hybrid Clipping [25]	2861	25.10
		Our Method	114	1.00
		Bézier Clipping [18]	6435	38.53
		Quadratic Degree Reduction [20]	7312	43.78
Example 2	${\mathcal{C}}_A$: 29/26/2	Cubic Degree Reduction [22]	7028	42.08
(Figure 5b)	${\mathcal{C}}_B$: 29/27/2	Geometry Interval Clipping [23]	4219	25.26
		Cubic Hybrid Clipping [25]	1873	11.22
		Our Method	167	1.00
		Bézier Clipping [18]	2176	3.74
		Quadratic Degree Reduction [20]	4440	7.63
Example 3	${\mathcal{C}}_A$: 63/60/8	Cubic Degree Reduction [22]	3384	5.81
(Figure 5c)	${\mathcal{C}}_B$: 24/21/8	Geometry Interval Clipping [23]	2066	3.55
		Cubic Hybrid Clipping [25]	1667	2.86
		Our Method	582	1.00
		Bézier Clipping [18]	1863	4.05
		Quadratic Degree Reduction [20]	3948	8.58
Example 4	${\mathcal{C}}_A$: 63/60/8	Cubic Degree Reduction [22]	2970	6.46
(Figure 5d)	${\mathcal{C}}_B$: 24/21/8	Geometry Interval Clipping [23]	1558	3.39
		Cubic Hybrid Clipping [25]	1184	2.57
		Our Method	460	1.00

presents a comprehensive comparison of the computational efficiency between our proposed method and several established techniques, including Bézier Clipping [18], Quadratic Degree Reduction [20], Cubic Degree Reduction [22], Geometry Interval Clipping [23], and Cubic Hybrid Clipping [25]. The evaluation encompasses a range of scenarios, as depicted in the table and Figure 5, characterized by varying numbers of control points, Bézier segments, and intersections, thereby offering a broad perspective on performance across different curve complexities.

The findings underscore the superior performance of our method, consistently achieving the lowest computation times (highlighted in bold) across all examples. Specifically, in Example 1 (Figure 5a), our method demonstrates a significant efficiency advantage, achieving a computation time of only 114 microseconds, which is substantially lower than the next fastest method, Bézier Clipping, by a factor of 1.84. This trend of outperformance persists across the other examples, with our method maintaining the position of optimal efficiency with the lowest relative times of 1.00 in all cases.

Example 2 (Figure 5b) further illustrates the robustness of our approach, handling complex curves (29 control points, 26/27 Bézier segments, and 2 intersections) with a computation time of merely 167 microseconds, whereas the competing methods exhibit significantly higher times, up to 43.78 times slower in the case of quadratic degree reduction.

In Examples 3 (Figure 5c) and 4 (Figure 5d), which involve an even greater number of intersections and control points, our method continues to exhibit unparalleled efficiency, further evidencing its potential as a superior computational tool for the intersection of Bézier and B-spline curves.

Traditional techniques require computing intersections for each pair of Bézier segments, leading to a computational complexity of $\mathcal{O}(\#BezSe{g}^2)$, with $\#BezSeg$ denoting the number of Bézier segments. For instance, in Example 3, curve ${\mathcal{C}}_A$ comprises 60 Bézier segments, as illustrated in Figure 6, whereas curve ${\mathcal{C}}_B$ contains 21 segments. Consequently, the process of determining Bézier intersections must be executed as many as $60\times 21=1260$ times. As the number of Bézier segments grows—reflecting an increase in the number of inner knot values—the computational demand escalates.

In contrast, our method is engineered to function independently of Bézier curve segmentation. As demonstrated in Figure 7, the computation time of our method remains constant regardless of the increase in Bézier segment count. However, the computation times of the other methods rise linearly with $\#BezSe{g}^2$, which coincides with the theoretical computational complexity analysis. This design choice eliminates the necessity for intersection calculations between segment pairs, thereby significantly reducing computational overhead. This strategic departure from segment-pair dependency not only streamlines the process but also enhances efficiency, providing a robust solution for handling complex curve intersections with greater computational economy.

The comparative analysis not only highlights the computational merits of our proposed method but also underscores the significance of algorithmic efficiency in the realm of geometric computation, especially for applications requiring the rapid processing of complex B-spline curve intersections.

However, it is crucial to highlight that the selection of parameters in the algorithm presented in this paper is heuristic. A larger tolerance will result in the obtained intersection point deviating from the actual position. Conversely, a smaller tolerance will increase the number of iterations, thereby reducing the speed of the solution. The selection of the most suitable parameters is therefore imperative for users, following the specific circumstances of the problem. Furthermore, the subdivision method based on pseudo-curvature is only applicable to parametric curves formed by control points and basis functions.

6. Conclusions

In this paper, we develop a rapid and robust algorithm designed to address the intersection challenges of Bézier and B-spline curves, which are fundamental components in the realms of computer graphics and computer-aided geometric design. By incorporating Sturm’s theorem, we have established a methodical approach to determine the number of intersections within a predetermined region. To handle multiple intersections, a fast subdivision-based algorithm is proposed. The bounding box detection and the definition of pseudo-curvature for curve segmentation further refine the process, enabling the partitioning of curves into segments that closely emulate straight lines. The deployment of the projected Gauss-Newton method facilitates efficient convergence to intersection points, marking a significant advancement over traditional techniques.

For future work, the exploration of efficient and robust algorithms for the surface/surface intersection problem presents a more challenging and interesting frontier.