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.

Introduction

Given an interval [α, β] ⊂ ℝ, a Bézier curve P(t), t ∈ [α, β] is defined aswhere are Bernstein polynomials in [α, β], and are the control points [1]. Given two Bézier curves P(t), t ∈ [α, β] and Q(s), s ∈ [ξ, η], the problem considered in the present study is how to compute all (t∗, s∗), t∗ ∈ [α, β], s∗ ∈ [ξ, η] such that P(t∗) = Q(s∗). Various methods have been developed to solve this task, such as a subdivision-based approach [2], binary subdivision approach [3], implicitization [4] and Bézier clipping [5].

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 k-th iteration interval is denoted as [α, β] × [ξ, η]. Let h = β − α and d = η − ξ. If a constant γ exists such thatwhere C are constants independent of k and the curves, and thus γ is the convergence rate of the sequence intervals {[α, β]} and {[ξ, η]}. The key problem is to find an algorithm for which γ 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 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 second-order 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.

Methods

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.

Fig. 1

A cubic Bézier curve expressed as a quadratic hybrid curve

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 at t = 0.5, is first evaluated and the resulting point is then used to evaluate , 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 , there exists an equivalent cubic hybrid curve  with two fixed control points  and two moving control points . The two moving control points  are Bézier curves of degree n − 3 with control points  and , respectively, where

a =  − (n − i)(n − i − 1)(n − i − 2), b = n(n − 1)(n − 2), c =  − i(i − 1)(i − 2), and i ∈ {1, …, n − 1}.

Proof. The degree n × m tensor product Bézier surface patch [1] is defined aswhere is the product of the two Bernstein bases in [0, 1], and Q, 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 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.,

If m = 3, then j ∈ {0, 1, 2, 3} and (j, k) ∈ {(0, i), (1, i − 1), (2, i − 2), (3, i − 3)}. Expanding the summation and rearranging the terms, the following is obtained:

If the control points P of degree n = d + 3 diagonal curve P(t) are known, Q0,  = P0 and Q3,  = P can be set. Thus,

Simplifying the above formulas, the following is achieved:for i ∈ {1, …, n − 1}. Setting a =  − (n − i)(n − i − 1)(n − i − 2), b = n(n − 1)(n − 2), c =  − i(i − 1)(i − 2), and observing that a + b + c = 3i(n − i)(n − 2), the following occur:

Because Q0,  = P0 and Q3,  = P, the s = t diagonal curve of Q(s, t) can be evaluated using the following formula:where are the degree n − 3 Bézier curves comprising the control points and , respectively, where

This is a cubic hybrid curve with two moving control points , and fixed control points P0, P.

From Theorem 1, if i = 1 or i = n − 1, the first control point of and the last control point of 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 and the second moving control point is denoted by , then the control points of 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 and the first moving control point is denoted by , the control points of can be obtained from Eq. (4), i.e.,

Case 3

In general, and are moving to control points. Because they are equivalent to P(t),

Hence, λ ∈ [0, 1] exists such that

Through a simple approach,and obtainwhere i ∈ {0, …, n − 3}; in addition, and 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.

Fig. 2

Intersection of cubic strip [, ] of and fat line of

In each step, one curve is bounded by two lines, called fat lines, which were first introduced in ref. [5]. Let be a line passing through P0 and P of a degree n Bézier curve, P(t), and suppose has an implicit equation:

The fat line of P is defined as a regionwhere [dmin, dmax] = [min0 ≤ d(P), max0 ≤ d(P)], and d(P) = ax + by + c, P = (x, y).

The steps of Algorithm 1 are described in more detail in the following:

In line 1, if the intervals are within the specified level of accuracy, the parameter intervals corresponding to the intersection in line 21 can be directly obtained.

In line 2, the curve with a larger priority parameter interval is always clipped.

In line 5, is a cubic hybrid polynomial in the Bernstein basis,

where

In line 6, to obtain the cubic lower and upper bounds, defining as the intervals containing the coefficients and of , respectively, is bound using an interval Bernstein polynomial [16]:

The lower and upper bounds of are defined through cubic polynomials in a simple manner:

Hence, . To obtain a tighter bound , the following optimization function is used:where , and j ∈ {0, 1, …, n − 3}. In Eq. (22), if the second moving point of is a fixed point denoted by , and the first moving point is denoted by , the following is obtained:

If the first moving point of is a fixed point denoted by , and the second moving point is denoted by , 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 . If a ≥ b, λ = 1 is set in Eq. (28). Otherwise, λ = 0. The tighter cubic bounds of can then be obtained using Eqs. (23) and (24).

In lines 7–10, the intervals of t for which P(t) lies outside of L correspond to regions where , as shown in Fig. 3. The values t ∈ [0, 1], for which and cross dmin and dmax, correspond to the roots of

Fig. 3

Clipping computed from

Because and are cubic polynomials, these roots can be solved directly using the cubic formula.

In line 7, if the intersection is empty, no intersection exists between the two curves. In line 11, if the lengths of these intervals are sufficiently small compared to the previous intervals [α, β], HybridClip is applied to line 14. Otherwise, the curve is subdivided into two subsegments and applies HybridClip to the two halves (line 12).

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 passes through the end control points P0 and P 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, is represented using the implicit equation

The fat plane containing curve P(t) and its control points are defined aswhere [dmin, dmax] = [min0 ≤ d(P), max0 ≤ d(P)], and d(P) = ax + by + cz + e, P = (x, y, z). 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.

Results

Proof for the convergence rate

Although Yuan’s method [14] is based on cubic HybClip, 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:

Lemma 1. For any given polynomial P, there exists two constants C and D depending solely on P, such that for all intervals [α, β] ⊆ [0, 1] the lower bound m and the upper bound M generated in line 6 of Algorithm 1 satisfywhere .

Proof. According to Eqs. (22) and (24), P(α) = m(α), P(β) = m(β), and P(t) ≥ m(t), ∀ t ∈ [α, β], and thuswhere P1(t) is a continuous function of degree n − 2, and m1(t) is a linear function. Let g(t) = b0(β − t) + b1(t − α) be a line passing through the lowest control point and parallel to the line connecting the end points of P1(t), such that P1(t) − g(t) ≥ 0, ∀ t ∈ [α, β], and thuswhere the constant C depends solely on P.where are the control points of g after the degree elevation [1], , and P3(t) = (a − c)(t − α) ≥ 0.

Let t1, t2 be the minimum values of P2(t), P3(t) in [α, β], respectively, i.e.,where s1, s2 ∈ [α, β]. Hence,

From Eqs. (34), (35), and (38),

Similarly, |M(t) − P(t)| < Dh4.

Lemma 2. For any given polynomial P, there exist constantsdepending solely on P, such that for all intervals [α, β] ⊆ [0, 1] the lower bound m and upper bound M generated in line 6 of Algorithm 1 for ∀i ∈ {0, 1, 2, 3} satisfywhere.

Proof A new norm in [α, β] is introduced as

According to the equivalence of norms in a finite-dimensional real linear space, there exists a constant C such thatwhere the constant C does not depend on the intervals [α, β], again owing to the affine invariance. Using arguments similar to those in the previous proof, let r = P − m,

Similarly, .

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 ([α, β]), which converges to t∗, is a sequence of intervals generated by Algorithm 1, with lengths h = β − α. 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 [α, β] based on line 6 of Algorithm 1. Because P′ is continuous, and owing to Lemma 2, the following inequalitieshold for all but a finite number of values of i. These three inequalities above imply thatand hence

From Lemma 1, the vertical height δ = δmin + δmax of m and M is bounded by . Thus, the length h of the intervals satisfiesfor all but a finite number of values of i (Fig. 4).

Fig. 4

Illustration of Theorem 2

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.

Fig. 5

Illustration of Theorem 3

Theorem 3. (Double root) If the polynomial P has a root t∗in [α, β], and provided that this root has multiplicity 2, 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 ([α, β]) is analyzed with lengths h = β − α 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 [α, β]. Because P″ is continuous, and based on Lemma 2, the inequalitieshold for all but a finite number of values of i. These two inequalities imply thatand thus . Letting τ = t − t∗, and based on

. From Lemmas 1 and 2,

Letting t1, t2 be the roots of m, t∗ ∈ [t1, t2], and τ2 = t2 − t∗ > 0, τ1 = t1 − t∗ < 0, the following is obtained:

Because τ1 ≤ h and h → 0, D3|τ1| → 0,for a sufficiently large i. Therefore, , and hence

Similarly, the following bound for t2 is obtained:

Because τ1 < 0, τ2 > 0,

Hence, the sequence (h) has a convergence rate of 2.

From Theorems 2 and 3, it can be seen that the new algorithm has a higher convergence rate compared with geometry interval clipping [11] and quadratic clipping [8] when computing all roots of a univariate polynomial equation. The following theorem provides the convergence rate for the curve/curve intersection problems.

Theorem 4. Suppose f(t), g(s) have a transversal intersection (f′(t∗) × g′(s∗) ≠ 0) at p∗ = f(t∗) = g(s∗). Furthermore, supposing that [α, β]is the sequence of generated intervals that contain t∗, and [ξ, η]is the corresponding sequence of generated intervals that contain s∗, there then exist constants C1, C2, C3, C4 depending solely on f and g, such that

Proof From line 11 of Algorithm 1, it can be seen that the length of intervals [ξ, η] tends toward zero as i tends toward infinity, that is, the interval [ξ, η] tends toward s∗.

Let be the line or plane that passes through the endpoints b0, b of g in [ξ, η]. Denote n as the unit normal vector of . Then, the distance function from f(t) to is defined as

Denote as the tangent line of f at t∗. Let be the angle between and , and be the angles between and b0b. As h = [ξ, η] tends toward 0, the line or plane converges at b0b, and angle φ converges at θ. Thus, for a sufficiently small h, , and thus .

The angle ρ between f′(t∗) and n is either or . Using this, the derivative of the distance function can be bound at the intersection as

Because d′(t∗) ≠ 0, and for convenience, w = d′(t∗) > 0 is denoted (otherwise, the vector −n can be considered instead of n).

Because d′(t) is continuous, the inequalityholds for all but a finite number of values of i. Hence,

From line 6 of Algorithm 1, the cubic polynomial bound [d(t), d(t)] of the distance function d(t) can be obtained. Based on Lemma 2,and by Eq. (61), the following is obtained:

From Fig. 6, the bound for h is obtained as

Fig. 6

Illustration of Eq. (64)

Based on Lemma 1, the vertical heights δ of d and d are bounded as follows:

Let t1 and t2 be the roots of d and d respectively. From Eq. (63),

From Eq. (64), the above inequality implies that

This thus implies the first inequality in Eq. (57) from (see the proof of Theorem 6 in ref. [7]). Similarly, in the next iteration step, the following is obtained:where C′ is solely dependent on f, and C′ 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 l1 and l3. Therefore, a transversal intersection is required in the proof. From Theorem 4, the two sequences {[α, β]} and {[ξ, η]} 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].

For convenience, denote Bézier clipping as BezClip [5]; quadratic clipping and cubic clipping based on a degree reduction as 2-DegClip [8] and 3-DegClip [10], respectively; geometry interval clipping as 2-HybClip [11]; and cubic HybClip based on hybrid curves in ref. [14] as 3-HybClip*. In addition, the proposed cubic HybClip algorithm is denoted as 3-HybClip.

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 a 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.

Table 1

Transversal intersections

(n, m)	ε	10−6					10−10
		BezClip	2-HybClip	3-HybClip	2-DegClip	3-DegClip	BezClip	2-HybClip	3-HybClip	2-DegClip	3-DegClip
(4,4)	[NP, NQ]	[4,3]	[3,3]	[3,2]	[3,3]	[3,2]	[16,15]	[4,3]	[3,3]	[4,3]	[3,3]
	t/μs	125	120	120	150	150	90	90	100	125	130
(8,4)	[NP, NQ]	[4,3]	[3,3]	[3,2]	[3,3]	[3,2]	[15,14]	[4,3]	[3,2]	[4,3]	[3,3]
	t/μs	160	150	150	190	190	100	120	100	140	130
(8,8)	[NP, NQ]	[4,4]	[4,4]	[3,3]	[3,3]	[3,3]	[16,16]	[5,4]	[3,3]	[4,4]	[3,3]
	t/μs	190	180	180	250	250	190	130	140	180	170

Number of iterations [N, N] and computing time t (in microseconds) of accuracy ε. In addition, (n, m) are degrees of , , respectively

Fig. 7

Computing time t in milliseconds vs accuracy ε of curve pairs with degrees (4, 4), (8, 4), and (8, 8) from the top to down. a Transversal intersections; b Tangent intersections

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.

Table 2

Tangent intersections

(n, m)	ε	10−6					10−10
		BezClip	2-HybClip	3-HybClip	2-DegClip	3-DegClip	BezClip	2-HybClip	3-HybClip	2-DegClip	3-DegClip
(4,4)	[NP, NQ]	[12,10]	[6,5]	[6,5]	[6,6]	[5,5]	[20,20]	[6,6]	[6,6]	[7,6]	[6,5]
	t/μs	138	117	93	119	138	164	129	85	124	136
(8,4)	[NP, NQ]	[20,20]	[8,8]	[10,10]	[6,6]	[5,6]	[34,34]	[9,8]	[11,10]	[7,6]	[5,7]
	t/μs	160	175	170	190	200	160	175	185	232	246
(8,8)	[NP, NQ]	[18,17]	[7,7]	[10,9]	[7,7]	[4,6]	[32,30]	[8,7]	[10,10]	[8,7]	[5,6]
	t/μs	170	185	180	226	263	177	157	133	272	282

Number of iterations [N, N] and computing time t (in microseconds) of accuracy ε. In addition, (n, m) are degrees of , , respectively

To compare these six algorithms numerically, statistics are generated on 40,000 pairs of randomly generated polynomial curves of degree 4–10 for single and multiple intersections. Figure 8 shows computing time needed to achieve the given accuracy of the five algorithms. The relative computing iterations and computing time for these tests are listed in Table 3. As shown in Table 3, 3-HybClip requires 2% fewer iterations and 8% less time than 3-HybClip* [14]. In addition, 3-HybClip has 2% fewer computing iterations than 3-DegClip, and at least 10% fewer iterations than 2-HybClip and 2-DegClip. With respect to the computing time, 3-HybClip is at least 60% faster than 3-DegClip and 2-DegClip, and at least 30% faster than 2-HybClip.

Fig. 8

Statistical comparisons: Computing time t in seconds vs accuracy ε. a Single intersections; b Multiple intersections

Table 3

Relative computing iterations N and computing time t

ε		BezClip	2-HybClip	3-HybClip*	3-HybClip	2-DegClip	3-DegClip
10−6	N	1.36	1.14	1.02	1	1.18	1.02
	t	1.85	1.39	1.08	1	1.95	1.71
10−10	N	1.79	1.16	1.03	1	1.18	1.03
	t	2.60	1.37	1.09	1	1.88	1.66

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 polynomials. Then the intersection of “fat planes” and the two cubic polynomials is the strip domain containing the intersections. The details of the algorithm and comparisons with previous approaches are left for future work.