Curve intersection based on cubic hybrid clipping

This study presents a novel approach to computing all intersections between two Bézier curves using cubic hybrid clipping. Each intersection is represented by two strip intervals that contain an intersection. In each step, one curve is bounded by two fat lines, and the other is bounded by two cubic Bézier curves, clipping away the domain that does not contain the intersections. By selecting the moving control points of the cubic hybrid curves, better cubic polynomial bounds are obtained to make the proposed method more efficient. It was proved that the two strip intervals have second- and fourth-order convergence rates for transversal intersections. Experimental results show that the new algorithm is the most efficient among all existing curve/curve intersection approaches.

The most common approach consists of clipping away the regions of the curves that are guaranteed to not intersect. Each intersection parameter pair (s * , t * ) is replaced with an interval that is iteratively computed. The The key problem is to find an algorithm for which γ i is as large as possible with as few computations as possible during each iteration. The above problem plays an important role in many engineering fields, such as computer-aided design and manufacturing (CAD/CAM), collision detection, and geometric modeling [1], and is a basic operation in solid modeling. In geometric processing, the intersections and intersection curves in a solid model are extremely important for the visualization, analysis, and manufacturing of the model [6]. With the continuous development of computer-aided geometric design and CAD/CAM, as well as the continuous progress made in science and technology, the numbers of calculations and data to be processed for intersection problems are increasing. It is therefore important to develop efficient and stable methods for dealing with intersection problems.
To solve such problems, the Bézier clipping algorithm introduced in ref. [5] is a widely used, fast, and (2) Page 2 of 13 Wu and Li Visual Computing for Industry, Biomedicine, and Art (2022) 5:17 robust method. The Bézier clipping algorithm was proven to have a second-order convergence rate [7]. Subsequently, several different approaches have been proposed to improve the Bézier clipping algorithm. Bartoň and Jüttler [8] and Liu et al. [9] developed quadratic and cubic clipping techniques based on a degree reduction to compute all roots of a univariate polynomial equation. Lou and Liu [10] extended the approach in ref. [9] to curve/curve intersection problems and proved that the algorithm achieves at least a secondorder convergence rate. In addition, North [11] developed a geometry interval clipping algorithm based on quadratic hybrid curves [12] for use with curve/curve intersection problems, Liu and Li [13] proved that the algorithm achieves a quadratic convergence rate. Moreover, Yuan [14] recently developed a cubic hybrid clipping (HybClip) based on hybrid curves to compute all roots of a univariate polynomial equation with a numerically verified fourth-order convergence rate. In this study, the approach in ref. [14] is extended to handle curve/curve intersection problems. Unlike the approach in ref. [14], a better bound is chosen for cubic HybClip, and thus the algorithm requires 8% less time than a method that directly uses the cubic hybrid curve [14]. In addition, it is proved that the two sequences in the new clipping algorithm have second-and fourth-order convergence rates. Subsequently, a complete comparison is provided with all existing curve/curve intersection algorithms based on subdivisions on a random 40,000 curve/curve intersection database. The new algorithm requires 30% less time than the geometry interval clipping algorithm [11] and 60% less time than the cubic clipping algorithm [10].
The remainder of this paper is organized as follows. In Methods section, the cubic hybrid curves are presented with two moving control points, and the details of the curve/curve intersection algorithms are described when applied in both 2D and 3D. In Results section, a proof of the convergence rate of the new intersection algorithm is provided, and the six techniques are compared from various perspectives. Finally, some concluding remarks are provided in Conclusions section and areas of future work are discussed in Discussion section.

Hybrid curve
A hybrid curve refers to a curve with at least one moving control point, which is itself a parametric curve and shares one parameter with the hybrid. Sederberg and Kakimoto [12] originated the idea of using hybrid polynomial Bézier curves to approximate rational Bézier curves. Later, North [11] transformed all polynomial Bézier curves of degree d ≥ 2 into equivalent quadratic hybrid curves with a single moving control point and fixed endpoints. As an illustration, a simple quadratic hybrid curve was constructed with a single moving control point, equivalent to a cubic Bézier curve, as shown in Fig. 1.
To evaluate a point on a hybrid curve, the locations of all moving control points are first determined at the given parameter value t. Once the moving control points are determined, the hybrid curve can be  evaluated as a common curve. For example, to evaluate P (t) at t = 0.5, P 1 (0.5) is first evaluated and the resulting point is then used to evaluate P (0.5) , as shown in Fig. 1.
Using the same principles as in ref. [11], by properly selecting the moving control points, the hybrid curve can produce any traditional Bézier curve. In this study, a hybrid curve having the following form is focused on: Theorem 1. Given a degree n ≥ 3 Bézier curve P(t) with control points {P i } n i=0 , there exists an equivalent cubic hybrid curve P (t) with two fixed control points P 0 = P 0 ,P 3 = P n and two moving control points P 1 (t),P 2 (t). The two moving control points P 1 (t),P 2 (t) are Bézier curves of degree n − 3 with Proof. The degree n × m tensor product Bézier surface patch [1] is defined as where B n i (s)B m j (t), 0 ≤ s, t ≤ 1 is the product of the two Bernstein bases in [0, 1], and Q i, j , i = 0, …, n; j = 0, …, m are the control points of Q(s, t).
From ref. [15], a degree m + d Bézier curve P(t) with control points P i can be described as the diagonal curve P(t) = Q(t, t) of a degree m × d Bézier surface Q(s, t), i.e., Expanding the summation and rearranging the terms, the following is obtained: If the control points P i of degree n = d + 3 diagonal curve P(t) are known, Q 0, i = P 0 and Q 3, i − 3 = P n can be set. Thus, Simplifying the above formulas, the following is achieved: , and observing that a i + b i + c i = 3i(n − i)(n − 2), the following occur: Because Q 0, i = P 0 and Q 3, i − 3 = P n , the s = t diagonal curve of Q(s, t) can be evaluated using the following formula: where P 1 (t) andP 2 (t) are the degree n − 3 Bézier curves comprising the control points P 1,i−1 and P 2,i−2 , respectively, where This is a cubic hybrid curve with two moving control points P 1 (t),P 2 (t) , and fixed control points P 0 , P n .
From Theorem 1, if i = 1 or i = n − 1, the first control point of P 1 (t) and the last control point of P 2 (t) are fixed as follows: Theorem 1 indicates that the two moving control points are relevant, for which three cases are discussed:

Case 1
If the first moving control point is a fixed point denoted by Q 1 and the second moving control point is denoted by Q 2 (t) , then the control points Q 2,i n−3 i=0 of Q 2 (t) can be calculated from Eq. (4), as indicated by Yuan [14]:

Case 2
If the second moving control point is a fixed point denoted by R 2 and the first moving control point is can be obtained from Eq. (4), i.e.,

Case 3
In general, P 1 (t) and P 2 (t) are moving to control points. Because they are equivalent to P(t), Through a simple approach, and obtain where i ∈ {0, …, n − 3}; in addition, Q 1 ,Q 2,i and R 1,i ,R 2 are known from the two cases above, and based on Eqs. (4) and (13), the two moving control points depend on the value of λ.

Curve/curve intersection based on cubic HybClip
Given two Bézier curves P(t), t ∈ [α, β] and Q(s), s ∈ [ξ, η], in this section, a cubic hybrid clipping algorithm is proposed for computing all intersections.

2D curve/curve intersection
The algorithm for two planar Bézier curves is first discussed. This algorithm is presented in Algorithm 1, and illustrated in Fig. 2.
In each step, one curve is bounded by two lines, called fat lines, which were first introduced in ref. [5]. Let L be a line passing through P 0 and P n of a degree n Bézier curve, P(t), and suppose L has an implicit equation: The fat line of P is defined as a region The steps of Algorithm 1 are described in more detail in the following: The lower and upper bounds of d (t) are defined through cubic polynomials in a simple manner: To obtain a tighter bound d (t) , the following optimization function is used: denoted by q 2 , and the first moving point is denoted by q 1 (t) , the following is obtained: If the first moving point of d (t) is a fixed point denoted by r 1 , and the second moving point is denoted by r 2 (t) , the following is obtained: There exists λ ∈ [0, 1] such that From Eqs. (19) and (25), the problem becomes linear, i.e., where j ∈ {0, 1, …, n − 3}. Let = max j q 1,j − min j q 1,j and b = max j r 2,j − min j r 2,j . If a ≥ b, λ = 1 is set in Eq.

3D curve/curve intersection
The above algorithm can be naturally generalized to handle 3D Bézier curve/curve intersection problems. In the 3D case, "fat lines" in 2D are replaced with several bounding planes, which are called "fat planes. " Plane L passes through the end control points P 0 and P n of a degree n Bézier curve, P(t). Because a plane consists of three points that are not collinear, an arbitrary control point is simply chosen that is not on the endpoint line. Here, L is represented using the implicit equation The fat plane containing curve P(t) and its control points are defined as where [d min , d max ] = [min 0 ≤ i ≤ n d(P i ), max 0 ≤ i ≤ n d(P i )], and d(P i ) = ax i + by i + cz i + e, P i = (x i , y i , z i ). The distance from one curve in a cubic hybrid form is then bound to the fat plane using two cubic polynomials and a strip domain containing the intersections is computed, which is similar to Algorithm 1 described in 2D curve/curve intersection section.

Proof for the convergence rate
Although Yuan's method [14] is based on cubic Hyb-Clip, it is mainly used to solve univariate polynomial root problems. However, a theoretical convergence rate or proof is not provided.
In this section, the theoretical results are provided on the convergence rate of the new curve/curve intersection algorithm. This begins with two technical lemmas: where P 1 (t) is a continuous function of degree n − 2, and m 1 (t) is a linear function. Let g(t) = b 0 (β − t) + b 1 (t − α) be a line passing through the lowest control point and parallel to the line connecting the end points of P 1 (t), such that P 1 (t) − g(t) ≥ 0, ∀ t ∈ [α, β], and thus where the constant C depends solely on P.
where {c i } n−2 i=0 are the control points of g after the degree elevation [1], and ∞ ≤ D P h 4 . From the above lemmas, the convergence rate can be analyzed using the HybClip algorithm. In Algorithm 1, if Q = 0, the curve/curve intersection problem P(t) = Q(s) becomes a root-finding problem P(t) = 0; that is, the cubic HybClip technique may be used to compute the roots of the polynomials and intersections of the two curves. These two cases are discussed separately.
Theorem 2. (Single root) If polynomial P has a root t * in [α, β], and provided that this root has multiplicity 1, the sequence of the lengths of the intervals generated through cubic HybClip containing that root has the convergence rate d = 4.
Proof. Suppose that ([α i , β i ]) i = 0, 1, 2, … , which converges to t * , is a sequence of intervals generated by Algorithm 1, with lengths h i = β i − α i . It is assumed that the first derivative satisfies P ′ (t * ) > 0 (otherwise, the polynomial −P can be considered instead of P).
Two cubic Bernstein polynomials m and M can be obtained as the lower and upper bounds of P in [α i , β i ] based on line 6 of Algorithm 1. Because P ′ is continuous, and owing to Lemma 2, the following inequalities hold for all but a finite number of values of i. These three inequalities above imply that and hence From Lemma 1, the vertical height δ = δ min + δ max of m and M is bounded by C P h 4 i . Thus, the length h i of the intervals satisfies for all but a finite number of values of i (Fig. 4).
For other clipping techniques [8,9], multiple roots reduce the convergence rate. The convergence rate of cubic HybClip is now discussed in the case of double roots, as illustrated in Fig. 5.
Theorem 3. (Double root) If the polynomial P has a root t * in [α, β], and provided that this root has multiplicity 2, i the sequence of the lengths of the intervals generated by cubic HybClip containing that root has the convergence rate d = 2.
Proof Similar to the proof of the previous theorem, the sequence of intervals ([α i , β i ]) i = 0, 1, 2, … is analyzed with lengths h i = β i − α i generated by Algorithm 1, which contains the double root. It is assumed that the second derivative satisfies P ″ > 0. Otherwise, polynomial −P can be considered instead of P.
Again, two cubic Bernstein polynomials m and M can be obtained as the lower and upper bounds of P in [α i , β i ]. Because P ″ is continuous, and based on Lemma 2, the inequalities hold for all but a finite number of values of i. These two inequalities imply that and thus m ′′ (t) Letting t 1 , t 2 be the roots of m, t * ∈ [t 1 , t 2 ], and τ 2 = t 2 − t * > 0, τ 1 = t 1 − t * < 0, the following is obtained: Page 9 of 13 Wu and Li Visual Computing for Industry, Biomedicine, and Art This thus implies the first inequality in Eq. (57) from d max − d min < h 2 i,g C g (see the proof of Theorem 6 in ref. [7]). Similarly, in the next iteration step, the following is obtained: where C′ f is solely dependent on f, and C′ g is solely dependent on g. Based on the first inequality of Eq. (57), the following is obtained: which implies the second inequality.
Note that the property of w being nonzero is key to binding l 1 and l 3 . Therefore, a transversal intersection is required in the proof. From Theorem 4, the two sequences {[α i , β i ]} i and {[ξ i , η i ]} i of the new intersection algorithm have second-and fourth-order convergence rates, respectively, and the 3D curve intersection problem yields the same results.

Experimental results
In this section, all six algorithms are compared based on three criteria: the amount of time per iteration step, the number of iterations, and the computing time required to achieve a certain accuracy. All algorithms were implemented in C++ on a PC with an 2.60-GHz Intel (R) Core (TM) i7-9750H CPU and 16.0 GB of RAM. In all experiments, both curves P(t) and Q(s) have a parameter domain [0, 1].
To analyze the relationship between the computational effort and the desired accuracy, two examples representing polynomials with transversal and tangent intersections are discussed. The five algorithms are first applied     to three pairs of Bézier curves with a transversal intersection. Table 1 reports the number of pairs of iterations and the computing time in microseconds of the desired accuracy for computing the transversal intersections between the three curve pairs with various degrees. Figure 7a shows the relationship between the computing time and desired accuracy, and indicates that 3-HybClip based on cubic hybrid curves is significantly improved in comparison with BezClip, 2HybClip, 2-DegClip, and 3-DegClip.
(70)  The five algorithms are applied to three pairs of Bézier curves with tangent intersections. Table 2 and Fig. 7b report the number of pairs of iterations and the computing time in milliseconds of various desired accuracies ε for computing the tangent intersections between the three curve pairs with various degrees. Experimental results show that the quadratic and cubic clipping techniques are better than Bézier clipping; however, compared with quadratic clipping based on hybrid curves or a degree reduction, the cubic clipping techniques show no substantial improvements. This is due to the fact that all clipping algorithms require a large number of binary subdivisions for tangent intersections.

Conclusions
In this study, an algorithm called 3-HybClip was derived for computing all intersections between two Bézier curves within a given domain. By selecting the moving control points, better bounds were obtained than those in ref. [14]. It was proved that the two sequences of bounded intervals for intersections have second-and fourth-order convergence rates for transversal intersections. The experimental results show that the newly proposed 3-HybClip algorithm requires 2% fewer iterations and 8% less time than 3-HybClip* from ref. [14], 10% fewer iterations than 2-HybClip and 2-DegClip, and at least 30% less time than other techniques such as BezClip, 2-HybClip, 2-DegClip, and 3-DegClip.

Discussion
As discussed in 3D curve/curve intersection section, for 3D curve/curve intersection problems, the "fat planes" are computed to bound a 3D Bézier curve. The distance from one curve in a cubic hybrid form to the fat plane is bound by two cubic polynomials, and a strip domain containing the intersections is then computed. Similarly, in curve/surface intersection problems, "fat planes" can also be used to bind a Bézier surface, and then the distance from the curve to the fat plane is bound by two cubic