Topology Guaranteed B-Spline Surface/Surface Intersection

The surface/surface intersection technique serves as one of the most fundamental functions in modern Computer Aided Design (CAD) systems. Despite the long research history and successful applications of surface intersection algorithms in various CAD industrial software, challenges still exist in balancing computational efficiency, accuracy, as well as topology correctness. Specifically, most practical intersection algorithms fail to guarantee the correct topology of the intersection curve(s) when two surfaces are in near-critical positions, which brings instability to CAD systems. Even in one of the most successfully used commercial geometry engines ACIS, such complicated intersection topology can still be a tough nut to crack. In this paper, we present a practical topology guaranteed algorithm for computing the intersection loci of two B-spline surfaces. Our algorithm well treats all types of common and complicated intersection topology with practical efficiency, including those intersections with multiple branches or cross singularities, contacts in several isolated singular points or highorder contacts along a curve, as well as intersections along boundary curves. We present representative examples of these hard topology situations that challenge not only the open-source geometry engine OCCT but also the commercial engine ACIS. We compare our algorithm in both efficiency and topology correctness on plenty of common and complicated models with the open-source intersection package in SISL, OCCT, and the commercial engine ACIS.

The surface/surface intersection technique serves as one of the most fundamental functions in modern Computer Aided Design (CAD) systems.Despite the long research history and successful applications of surface intersection algorithms in various CAD industrial software, challenges still exist in balancing computational efficiency, accuracy, as well as topology correctness.Specifically, most practical intersection algorithms fail to guarantee the correct topology of the intersection curve(s) when two surfaces are in near-critical positions, which brings instability to CAD systems.Even in one of the most successfully used commercial geometry engines ACIS, such complicated intersection topology can still be a tough nut to crack.
In this paper, we present a practical topology guaranteed algorithm for computing the intersection loci of two B-spline surfaces.Our algorithm well treats all types of common and complicated intersection topology with practical efficiency, including those intersections with multiple branches or cross singularities, contacts in several isolated singular points or highorder contacts along a curve, as well as intersections along boundary curves.We present representative examples of these hard topology situations that challenge not only the open-source geometry engine OCCT but also the commercial engine ACIS.We compare our algorithm in both efficiency and topology correctness on plenty of common and complicated models with the open-source intersection package in SISL, OCCT, and the commercial engine ACIS.
CCS Concepts: • Computing methodologies → Computer graphics; Shape modeling; Parametric curve and surface models; Algebraic algorithms.

INTRODUCTION
Surface/Surface intersection is a fundamental task in Computer Aided Design (CAD) and geometric modeling systems.An efficient, robust, and topology guaranteed surface/surface intersection algorithm is highly desired in many modeling operations of CAD systems, such as Boolean operations, mesh generation, and rendering.Given the prevalent utilization of surface/surface intersection in CAD systems, research on this task extends beyond the single goal of efficiency.Consequently, an algorithm that not only fulfills the fundamental requirements but also exhibits robustness by adapting to diverse scenarios is highly required.Additionally, the algorithm must guarantee the topological correctness by effectively handling intricate situations, such as small loops and tangential situations.
There exists extensive research on surface intersection algorithms starting from the 1980s, during which several seminal algorithms gained significant attention, i.e. the algebraic method, the marching method, the lattice method, the subdivision method, and the hybrid method that combines two of these methods.Algebraic methods usually convert the surface intersection problem into solving algebraic equations, which can be challenging in certain complex scenarios.The marching method computes a series of starting points and then traces out the intersection locus using classic differential geometry [Bajaj et al. 1988a; Barnhill and Kersey 1990;Kriezis et al. 1992].Nevertheless, a full computation of starting points on all branches, which directly determines the topology correctness of the intersection curve, is difficult to achieve.The lattice method decomposes the surface/surface intersection problem into curve/surface intersection problems by discretizing the parametric domain of the surface [Rossignac and Requicha 1987].However, inadequate discretization may result in the loss of small loops and isolated points in the intersection curve.The subdivision method keeps dividing the surfaces into smaller patches and then converts the intersection problem into the intersection of these near-planar patches [Patrikalakis 1993], which may lead to the inaccurate calculation of tangential curves with insufficient number of subdivisions.Hybrid methods combine two of the above-discussed methods, for example, the algebraic method and the marching method, or the subdivision method and the marching method.Among all these methods, algebraic-related methods possess the advantage of preserving the intersection topology, which is difficult for the lattice or subdivision method.However, the algebraic method is generally inefficient and depends heavily on some other algebraic techniques, such as implicitization or solving polynomial systems.
Although extensive exploration has been conducted on the surface intersection problem, the existing methodologies can still be improved in terms of implementation.For example, some algebraic methods resort to surface implicitization or solving a system of polynomial equations, which involve two key steps in an intersection algorithm but are actually hard problems in algebraic geometry.Over the past decades, the implicitization and the polynomial system solving problems have been greatly developed in their own community; nevertheless, these developments were seldom brought back to the surface intersection literature [Hoffmann 1989;Sarraga 1983].Additionally, records pertaining to the timing and topology of intersection examples are few, and a limited number of works [Jia et al. 2022;Lin et al. 2013] have performed comparative analyses to evaluate the performance of different methods.
The surface/surface intersection problem has not only been widely studied in academic research but rapidly developed in the industrial field.Many open-source and commercial software packages containing intersection functions have been developed and extensively employed in both research and industrial domains.Among these packages, ACIS [2023] is a renowned 3D modeling kernel that provides powerful capabilities for creating, manipulating, and analyzing geometric models.ACIS has established itself as one of the leading 3D modeling kernels in the CAD industry due to its robustness, versatility, and wide range of functionalities.The surface intersection function in ACIS has also become one of the industry standards due to its high robustness and efficiency.Nevertheless, when confronted with scenarios involving high-order contact or other complex intersection topology, the outcomes produced by ACIS are sometimes unsatisfactory as illustrated in Fig. 1.
In this paper, we present a practical topology guaranteed algorithm for computing the intersection of two B-Spline surfaces.The proposed algorithm demonstrates topology correctness not only in all common transversal examples of B-spline surface intersections as in ACIS but also in challenging intersection topologies, such as when the two surfaces are in high-order contact at some points or along a curve.We mainly adopt an algebraic routine in the algorithm by bringing novel techniques or improvements to state-of-the-art works in the key steps of implicitization, topology determination, and redundant curve clipping.Our main contributions are as follows: (1) The fast moving plane technique is introduced to the traditional implicitization approach by Dixon matrices, which facilitates the intersection process even in the presence of base points.
(2) The topology of the intersection curve in the parametric domain is pre-computed to ensure the correctness of the adjacency relationships between critical points and boundary points, as well as considering the intersection along the boundaries of the parametric domain.
(3) A novel clipping method is proposed for the redundant part of the intersection curve caused by implicitization.The method is based on the inversion formula deduced from the Dixon matrix and successfully restricts the intersection curve to its required parametric domain.

RELATED WORK
Extensive research has been conducted on the problem of surface/surface intersection, and most of the proposed methods or approaches can be adapted to the B-Spline surface intersection problem [Bajaj et al. 1988b;Dokken 1997;Krishnan and Manocha 1997;Manocha and Canny 1991;Patrikalakis 1993;Sederberg and Meyers 1988;Shen et al. 2016;Ventura and Guedes Soares 2012].These techniques can typically be classified into five categories: the algebraic method, the marching method, the lattice method, the subdivision method, and the hybrid method.Algebraic methods.Given two B-Spline surfaces, algebraic methods typically convert one of them into its implicit equation and then substitute the parametric form of the other surface into this implicit equation.This process results in an algebraic curve with preserved algebraic characteristics [Manocha and Canny 1991;Sarraga 1983].However, the implicitization of rational surfaces using the traditional algebraic method is known to be challenging [Busé et al. 2003].Moreover, the parametric domain of the surface will be omitted after the implicitization, which leads to undesired intersection curves outside the parametric domain of surfaces during the intersection process.Sederberg and Chen [1995] were the first to discover the technique of moving surfaces for surface implicitization, where they mentioned the Dixon matrix as a possible candidate for implicitization.However, at that time, the Dixon matrix always failed to provide implicitization when the surface contained base points where the homogenous coordinates of the surface simultaneously vanish.Therefore, the technique they proposed was the first to use blending functions to construct a series of moving surfaces for implicitization, which was a different and novel technique.However, despite the disadvantage of the Dixon matrix of not being directly applicable to surfaces with base points, it remains the simplest method for implicitization, because it has much smaller matrix coefficients than other implicitization approaches.
Marching methods for surface intersection are widely used because of their generality and simplicity.Marching methods usually have two main steps.First, a series of starting points are determined on each branch of the intersection curves.Then, by using the local differential geometry of surfaces, the intersection locus is traced out from each starting point [Bajaj et al. 1988b;Barnhill and Kersey 1990;Patrikalakis 1993;Sederberg and Meyers 1988;Ventura and Guedes Soares 2012].However, this approach also has some limitations, such as the necessity of identifying starting points, which can be challenging in the presence of small loops, and the determination of the tracing direction over singular points.
Lattice methods decompose the surface/surface intersection problem into multiple curve/surface intersection problems.In this approach, one of the surfaces is disintegrated into its isoparametric curves, which are then intersected with the other surface to determine their intersection points.The ultimate intersection curve of two surfaces is produced by linking the discrete intersection points [Rossignac and Requicha 1987].Owing to the discretization of the lattice method, the inadequate sampling density of the isoparametric curves selected on the surface may result in missing isolated intersection points or small loops.
Subdivision methods are a recursive approach that involves subdividing the B-Spline surfaces into smaller patches until the intersection can be computed with sufficient accuracy.This approach consists of two iterative steps.First, each B-Spline surface is divided into a set of smaller patches by refining its control points.Then, an intersection test is performed between the patches of the two surfaces to determine possible intersections.The accuracy and efficiency of subdivision methods heavily depend on the choice of subdivision strategy and termination criteria [Lasser 1986;Sederberg and Meyers 1988].[De Figueiredo 1996] utilized affine arithmetic to calculate the bounding box of small patches and then quickly discarded the pairs of patches with no box intersection.[Lin et al. 2013] accelerated the affine arithmetic intersection algorithm using a GPU, significantly enhancing its efficiency.Although subdivision methods can obtain intersection curves at relatively high speed, they may miss small loops and isolated intersection points, and the topological correctness near singular points is not guaranteed under this approach.
Hybrid methods pertain to the integration of two or more previously mentioned methods, with the objective of harnessing their distinct advantages.A common example is to use the subdivision method to find several starting points and then trace out the whole intersection curve using the marching method in order to enhance the efficiency of the intersection process [Barnhill and Kersey 1990;Dokken 1997].Another representative approach involves utilizing algebraic methods to locate starting points, and subsequently employing marching methods to trace out the intersection curve from these calculated starting points [Krishnan and Manocha 1997].The fusion of algebraic and tracing methods is frequently employed to strike a balance between the efficiency and topological accuracy of the intersection curve.

PRELIMINARIES
Before introducing the proposed intersection algorithm, we first introduce some preliminaries about base points, moving planes, and Dixon matrices.

Dixon Matrices
Given three bivariate polynomials of degree (, ): Fig. 2. Overview of the proposed algorithm.Surface P is represented in blue and Surface Q is represented in orange.First, a fast collision detection of the bounding boxes of the two surfaces is conducted.Then, if two bounding boxes collide with each other, Surface P is implicitized to  (, , ,  ) using its corresponding Dixon matrix  (, , ,  ).Subsequently, left helping points (red), right helping points (blue) and boundary points (green) are calculated to determine the topology of the intersection curve.These helping points then serve as the starting points for locus computation.After tracing the intersection curve in the parametric domain, we lift it to 3D space and employ the Dixon matrix to clip the 3D curve within the appropriate region.
the Dixon matrix  ( , , ℎ) is a 2 × 2 matrix given by the following equation [Dixon et al. 1908]: where  and  are two variables.The Dixon resultant of  , , ℎ is defined by the determinant det ( ( , , ℎ)) of the Dixon matrix.

B-SPLINE SURFACE INTERSECTIONS
A B-Spline surface of bi-degree (, ) is defined by ( + 1) × ( + 1) control points p , and two knot vectors and Riesenfeld 1974].The parametric representation of the B-Spline surface can be written as where  , () and  , () are the spline basis functions that can be derived by the Cox-de Boor recursion formula [de Boor 1971].
Given two B-Spline surfaces P and Q of bi-degrees ( 1 ,  1 ) and ( 2 ,  2 ), with control points P, Q and knot vectors Ū1 , V1 and Ū2 , V2 : the general outline of computing the intersection of P and Q is as follows, illustrated in Fig. 2.
Several key steps in the above outline, for example, implicitization, topology determination, and clipping are hard problems themselves and have been developed in their own research field over the decades.On one hand, these developments are seldom brought back to the surface intersection algorithm; on the other hand, more exploration is needed in these steps to give a detailed and practical algorithm of intersection computation.

Fast Determination of Intersections
A fast and rough collision detection is performed between the two given spline surfaces using their bounding boxes constructed from control points.The separation of the two bounding boxes indicates the separation of the two surfaces.Otherwise, the surfaces are further split into Bézier patches and similar collision detection is taken on each pair of the Bézier patches from the two splines.
The oriented bounding box (OBB) [Assarsson and Moller 2000] of the control points of the B-Spline/Bézier surface is adopted, shown in Fig. 3(a).The OBB for a given point set is the smallest rectangular box that contains the point set with its axes not necessarily aligned with the global coordinate system.Chang et al. [2011] presented a fast algorithm for computing the OBB of the control points of B-Spline surfaces.Compared to the axis-aligned bounding box (AABB) in Fig. 3(b), OBB provides a tighter bound for the surface, hence is more efficient in discarding the non-intersection patches.

Fast Implicitization of Parametric Surfaces
Since we reduce the intersection of two B-spline surfaces into the intersections of patch pairs from the two splines, we next focus on computing the intersection of two Bezier patches.To do this, we plug the parametric form of one patch into the implicit form of the other.Next, we explore the implicitization technique of moving planes.Given two parametric surfaces P(, ) and Q(, ) of bi-degree ( 1 ,  1 ) and ( 2 ,  2 ) with the following homogeneous parametrization: we convert one of the surfaces with the lower implicit degree1 , say P(, ), into its implicit equation  (, , , ) = 0.
In order to implicitize P(, ), we define three moving planes Obviously,  , , ℎ all follow surface P(, ).We compute the Dixon matrix of polynomials  , , ℎ with respect to variables (, ): Proposition 1.When the parametrization P(, ) has no base points, rank(D) = 2m 1 n 1 , and the determinant of  gives the implicit equation of the surface [Sederberg et al. 1984].Otherwise, if the parametrization contains base points, rank(D) < 2m 1 n 1 , and the implicit equation is a factor of the maximal nonzero minors of  [Manocha and Canny 1992].
Remark 1.Since the rank  of polynomial matrix  is computationally expensive, a practical alternative strategy is to assign some generic values of X  = (  ,   ,   ,   ),  = 1, • • • ,  to matrix  and compute   = rank(D(X i )).According to the fact that set  := {X 0 = ( 0 ,  0 ,  0 ,  0 ) ∈ R 4 |rank(D(X 0 )) < r} is a zero-measure set in R 4 , most values of   equal the rank of polynomial matrix  .Once rank  is determined, the implicit equation of the parametric surface can be computed.Here a fast algorithm for computing the Dixon matrix presented in [Chionh et al. 2002] is adopted.
then point X 0 is outside the given parametric domain of the surface P(, ).

Topology Determination of the Intersection Curve in the Parametric Domain
Next we determine the topology of the intersection curve in the parametric domain of the surface Q(, ).First, an implicit expression of the intersection curve in the parametric domain is computed; Then, the characteristic points of the intersection curve are computed; Finally, the connection of branches is determined, and the topology of the intersection curve is computed.
Let  ( ) = 0 be the implicit equation of the parametric surface P(, ).To compute the intersection of P(, ) and Q(, ), we substitute  = Q(, ) into  ( ) = 0: We next determine the topology of the curve We first compute two types of characteristic points: critical points and boundary points.For each  = 1, • • • , , we record the roots of     () = 0 and  ℎ () = 0 as: Now two sets of points on the curve C are obtained: which include all s-critical points and those points on C that have the same s-coordinates as s-critical points.In the following context, we call   and   left and right helping points, respectively.According to [Jin and Cheng 2021], once these left and right helping points are computed, the topology of the intersection curve in the parametric domain (, ) is determined.Refer to Fig. 14 and the appendix for the detailed example.

Tracing and Clipping the Intersection Curve
The advantage of the previous step in topology determination lies in that it greatly simplifies the locus tracing of the intersection curve  (, ) = 0 in the parametric domain (, ), especially in the neighborhood of singular points, where the standard tracing scheme [Barnhill et al. 1987] tends to lose its direction.Starting from our topology graph, tracing from a singular point, such as in Fig. 14(a) in the appendix, is decomposed by tracing from some left helping points or right helping points in a -neighborhood of this singular point along its related curve branches.For the choice of the tracing direction and step length, one can refer to [Chen et al. 1997] and the appendix.
In the following, we clip the obtained 3D intersection curve C (Fig. 7) such that for every point X 0 on C, its corresponding parameter values on P(, ) and Q(, ) are within (2) Clip C: For each split branch   ,  = 1, 2, • • • , , a random representative point q  on   is selected as shown in Fig. 9(c).We decide whether to discard or keep   based on the behavior of q  : Compute the parameters (, ) for q  by solving  (q  )U = 0 for (, ).According to Remark 2, if there is one solution in the branch   is retained; otherwise, it is discarded.Then C is clipped within the (, ) domain, see Fig. 9(d).
It is worth mentioning that after we split the original branches at the surface boundaries, a representative point of each split branch lies on the surface if and only if all points on the split branch lie on the surface.Hence, we focus only on those representative points.

EVALUATIONS AND COMPARISON
We implemented our algorithm and performed comparisons with the discrete intersection algorithm for mesh Booleans [Cherchi et al. 2022], the affine arithmetic-based intersection method with GPU acceleration [Lin et al. 2013  Table 1.The compared methods and their parameter settings in experiments.

Topology correctness
We first focus on the topology correctness of the intersection of two general B-spline surfaces, and compare our results with the open-source geometry kernel OCCT and the commercial engine ACIS.Plenty of experiments show that when the two surfaces intersect transversely, all three algorithms give the correct results.
However, when the two surfaces are in certain tangential or special intersection topology, our algorithm performs most robustly in these challenging situations, while OCCT and even ACIS sometimes fail.Fig. 1 shows representative examples with different intersection topologies.
In Fig. 1(a), the intersection curve contains two loops that are crossing at two common points.Our algorithm, OCCT and ACIS all give the correct intersection result; in (b), the two surfaces intersect at two separate curve branches and are in contact at an isolated point.OCCT loses this isolated contact point and gives only the two separate intersection curves, while our algorithm and ACIS give the correct result.The omission of tangent points in the OCCT intersection algorithm stems from its reliance on a discrete method for the initial point search.By contrast, our proposed algorithm utilizes polynomial root solving to guarantee the inclusion of all tangent points; in (c), the two surfaces are very close to each other in the triangular domain, and the intersection curve is a curved triangle, two edges of which are the boundary curves of the two surfaces and the left edge is a regular intersection curve.OCCT loses the two common boundary curves, while both our algorithm and ACIS give the complete result; in (d), the two surfaces are again very close to each other, and their intersection is even more complex than (c), that is, the intersection contains two quadrilaterals sharing one edge.The two surfaces are in contact along the common edge of the two quadrilaterals, and intersect transversely at the other six curved edges of the two quadrilaterals, which are also the boundaries of the surfaces.OCCT returns three of the six edges of the two quadrilaterals, but has lost the contact curve in the middle and the other three boundary curves, while both our algorithm and ACIS give the complete result.Examples (c) and (d) illustrate the limitations of the OCCT intersection algorithm when computing intersection curves at surface boundaries.The omission of intersection curves along the boundaries might result from improper termination criteria during the tracing process.Our proposed algorithm addresses this issue by incorporating specific handling for boundary intersection curves, thereby preventing the occurrence of this problem; in (e), the intersection curve contains two loops that are contact at one common point.OCCT fails to return any intersection point, while both our algorithm and ACIS give the correct result; in (f), the two surfaces are in high-order contact along a line.Neither OCCT nor ACIS returns any intersection point, while ours correctly returns the contact line.The failures of OCCT in examples (e) and (f), as well as ACIS in example (f), demonstrate that when two surfaces are tangent or in high-order contact, their normals at the tangent curves are the same.Existing algorithms exhibit deficiencies in calculating the marching direction during the tracing process.In contrast, our proposed algorithm performs curve tracing on planar polynomial curves in the parametric domain, thereby enabling the easy determination of the tracing direction; in (g), the two surfaces are in contact at three isolated points.Neither OCCT nor ACIS provides a complete set of contact points, while our algorithm ensures the detection of all three points of contact.The reason for the missing intersection points in the OCCT intersection algorithm is analogous to example (b), both are attributed to the discrete algorithm of the initial point detection.As for ACIS, the cause of losing one of the intersection points may be linked to the fact that this particular point also serves as a self-intersection point for the surfaces.The ACIS algorithm demonstrates shortcomings in handling such intricate intersection points.
Overall speaking, our algorithm performs most robustly in all these challenging situations in topology.ACIS outperforms OCCT, and gives the correct results in most situations, but fails at examples with high-order contact curves or multiple isolated contact points.The topological robustness of our algorithm lies in our algebraic treatment of the implicitization and the topology analysis of the intersection curve.On one hand, although traditional implicitization methods such as resultants and Gröbner basis computation are not very practical in both efficiency and numerical computations, the technique of moving planes, which we have adopted, has overcome these difficulties and paved the way for the application of implicitization to many geometry problems.On the other hand, the topological analysis and clipping method ensure topological correctness, especially for those tangential intersections or isolated contact points, which are really difficult to capture through the discrete techniques adopted by most software.

Comprehensive performance comparison
We have performed plenty of experiments on bi-cubic or bi-quadratic B-Spline surfaces, and parts of them are shown in Table 2. Since for practical use in CAD systems, efficiency is one of the most important factors in evaluating the intersection algorithm, we compare our algorithm with efficient and practical state-of-the-art intersection methods, including the affine arithmetic-based intersection method with GPU acceleration (AA-GPU), the mesh intersection approach (Mesh Booleans), the open source intersection packages in OCCT and SISL, and the commercial 3D ACIS modeler (ACIS), from the perspective of the computation accuracy, topology correctness and efficiency.We claim that our algorithm has the most robust topology correctness and comparable practical efficiency with the existing methods or packages.
In Table 2, six examples of B-Spline surface intersections are provided, in which our method always gives the correct topology, while the other five algorithms either fail or give the undesired intersection topology.The results are illustrated for every method, where 'NULL' means that the computation result returns as no intersections.The intersection results show that tangential intersections or isolated singular points are very challenging for all competitors.Our algorithm is particularly designed based on topology determination and inversion formulas, so tangential intersections can always be treated well.The Mesh Booleans method can obtain relatively smooth intersection curves when dealing with the transversal intersection of surfaces.However, this method tends to give a 2D grid region instead of a 1D curve along the tangential intersection curve or the isolated point.A 2D grid-like intersection is generated because the triangular meshes corresponding to the two surfaces cannot guarantee to pass the tangential curve or isolated point simultaneously, resulting in significant errors in the intersection lines of the meshes near these tangential curves or isolated points.AA-GPU fails to get correct results in examples (1∼5) when the intersection is tangential or near boundaries, and gives partial results in example (6) when the intersection has an isolated point.Since the two surfaces are very close to each other near the tangential curves, the overestimation of the bounding boxes with the affine arithmetic may contribute to the very thick strips of intersection grids on parametric regions.Thus, the undesired results may be attributed to the phenomenon in which the required accuracy near the boundary and tangential curves exceed the maximum subdivision accuracy during the subdivision of the parametric regions of B-Spline surfaces.OCCT fails to give the intersection result in examples (1)(2)(4), and gives partial intersection results but loses a whole branch in examples (3)(5)(6); SISL fails on all the six examples; even ACIS only returns several intersection points in examples (1)(2), returns no result in example (4), loses a whole branch in examples (3)(5), and loses the isolated singular point in example (6).The loss of the tangential curve in OCCT, SISL, and ACIS may be caused by the fact that the normal vectors of the two surfaces are parallel at every point on the tangential curve, which results in the inability of the 3D tracing step to track the entire intersection curve.Meanwhile, the omission of the isolated point may be caused by the use of some discrete methods in the search for the starting point, which results in the isolated point being ignored during the discretization process.
Table 2. Intersections of two B-Spline Surfaces.The first column shows the surfaces and their bi-degrees.The remaining columns represent the intersection results and the consuming times (in seconds) of different algorithms.AA-GPU refers to the affine arithmetic-based B-Spline surface intersection method with GPU acceleration.Mesh Booleans refers to the interactive and robust mesh boolean algorithm.OCCT refers to the intersection algorithm in the Open CASCADE Technology.SISL refers to the intersection algorithm in the SINTEF Spline Library.ACIS refers to the intersection algorithm in 3D ACIS Modeler.In the first column, the blue and orange surfaces depict the two surfaces to be intersected, and the red curve represents the intersection result of our method.The second column displays two surfaces in blue and red, with the intersection result of the AA method represented in green.N/A indicates that the program reported an error in this example, and the 'NULL' in the image indicates that the result obtained by this method is no intersection.

Intersection of Surfaces with Special Geometry
According to our method, given two parametric surfaces or spline surfaces, we transform one of them into its implicit expression.However, for specific surface types, such as quadrics, tori, cyclides, and surfaces of revolution, the implicit equations can be directly derived from their shape parameters, eliminating the need for the implicitization step.Consequently, our algorithm simplifies the intersection of these surfaces by directly deriving the implicit equation from the shape parameters.This streamlined approach replaces the conventional implicitization step in our implementation.
Table 3 shows the comparison of our method with ACIS in the intersection of an ellipsoid and a torus, where the implicit equation of the ellipsoid is directly written down (note that here we prefer to select the implicit equation of the surface with lower implicit degree).A comparison is made in the geometry of the intersection curve and the time cost.The shape parameters of the torus and ellipsoid of the corresponding examples are listed in Table 4.
It is worth mentioning that when dealing with surface intersections, ACIS classifies all surface types into planes, spheres, cones, tori, and splines.Ellipsoids are converted to spline surfaces in ACIS before the intersection computation.
Our extensive testing of intersections between ellipsoids and tori indicates that our method always achieves the same correct intersection geometry as ACIS.Table 3 shows only several representative examples.In terms of efficiency, if the intersection curve contains cross singular points, which are marked with red dots, our method Table 4.The parameters of torus and ellipsoid.In the Torus column, 'Center' represents the center of the torus, 'Normal' the normal of torus, 'R' the major radius and 'r' the minor radius.In the Ellipsoid column, 'Center' represents the center of the ellipsoid, 'x-axis' and 'y-axis' define a Cartesian coordinates frame, 'a', 'b' and 'c' are the length of the semi-axes.

Accelerating Tracing with TBB
During the process of tracing, the computation of each local branch (Fig. 7) is independent and uniquely determined by the starting point, the ending point, and the implicit equation, which have been calculated in the implicitization and topology determination steps.Consequently, in order to accelerate the tracing process, we employed the threading building blocks (TBB) library function, parallel_for, in the tracing step of local branches.Table 5 illustrates the comparison of the average time cost of the above-mentioned examples before and after TBB acceleration.Based on the obtained results, it can be inferred that when the intersection curve has only a single branch within the parametric domain (e.g.examples (1)(2)(4) in B-Spline surfaces), the employment of TBB doesn't improve the efficiency of our algorithm.Nevertheless, in cases where the intersection curve has multiple local branches within the parametric domain, the utilization of TBB enables simultaneous tracing of these local branches.This capability leads to a substantial reduction in algorithmic computation time, up to 2-3 times faster than the original calculation time.See Fig. 10 for an illustration of the local branches in the parametric domain.

Special Cases in Topology Determination
In

Applications
The intersection algorithm is crucial in the Boolean operations of solids.By computing the intersection curves of B-Spline patches in B-Rep solids, we can get different Boolean operation results such as union, intersection and subtraction.Fig. 12 shows the Boolean operation results of two hummingbirds using our intersection algorithm.We apply our algorithm to the Boolean operations of complex CAD models, and perform comparisons with OCCT and ACIS.When the adjacent surfaces in the CAD models intersect transversely, our algorithm gives the same Boolean operation results as OCCT and ACIS; nevertheless, when the adjacent surfaces are in high-order contact or intersect along their parametric boundaries, OCCT and ACIS sometimes provide incorrect results.Fig. 13 shows two such examples.
In the first example, we test the Boolean operations on two Cardan joint models.The intersection curves of OCCT are represented in green, where part of the intersection branch in the lower-left region and a tangential line in the middle region are missing.The miscalculation of the intersection curves leads to an error in the Boolean intersection performed by OCCT.The incorrect parts are highlighted within a blue box.Meanwhile, the incorrect Boolean subtraction results in OCCT arise from the inaccurate intersection calculation and the erroneous discarding of surface portions of solid A. The intersection curves of ACIS are represented in yellow.ACIS returns the transversal intersecting curves, but misses the tangential intersection curve, which is highlighted in blue in our corresponding intersection result.The missing tangential curve leads to the wrong omission of surfaces in the Boolean intersection results and the failure to subtract the portions of solid A that lie within solid B in Boolean subtraction results.The missing surface in the Boolean intersection is shown in the blue box, while the surface that is not properly trimmed in the Boolean subtraction is represented in green.
In the second example, we conduct Boolean operations on an additional pair of mechanical component models.OCCT mistakenly calculates two extra intersecting lines near the tangential intersection, which leads to inaccuracies of surfaces near the tangent region in both Boolean union and Boolean intersection results.The surplus intersections also lead to an erroneous trimming of a surface in solid A in the Boolean subtraction outcome.The intersection result of ACIS fails to capture a tangent intersection, subsequently impacting the surface classification process within the Boolean union and Boolean intersection outcomes, leading to a partial omission in the Boolean union and an excessive portion in the Boolean intersection.The Boolean subtraction result of ACIS also leads to errors due to the same reason.Contrary to OCCT and ACIS, our algorithm accurately calculates the intersection curves of the tangential and boundary regions of surfaces, resulting in the correct Boolean union, intersection, and subtraction outcomes.

CONCLUSIONS AND FUTURE WORK
We present a topology guaranteed algorithm for computing the intersection of two B-Spline surfaces.The algorithm begins by performing fast collision detection using the oriented bounding boxes to eliminate non-intersection patch pairs.Subsequently, the implicitization for one Bézier surface within each Bézier patch pair is computed efficiently using the Dixon matrix of moving planes.Following this, a topology determination strategy is applied to the intersection curve in the parametric domain and the helping points on the curve are calculated for topology guaranteed tracing.Finally, a clipping method based on the inversion formula of the Dixon matrix is demonstrated to restrict the intersection curve within the parametric domain of two surfaces.We have illustrated the topological correctness of our method with various examples of B-Spline surface intersection with different topologies.Furthermore, a detailed comparison of the efficiency and topology correctness of our algorithm with the affine arithmetic-based method, mesh-to-mesh intersection algorithm and intersection commands in OCCT, SISL and ACIS is presented, which indicates that our method can maintain topology correctness even when two surfaces are under special relative position.
In the proposed method, we derive the implicit equation of the parametric surface from its Dixon matrix and then use it to determine the critical points of the intersection curves.However, the Dixon matrix contains all the necessary surface information required for performing intersection computations.Thus, our future research will focus on simplifying this process by directly computing the critical points of the intersection using the Dixon matrix, thereby enhancing the efficiency of our algorithm.Moreover, in our current work, we convert two B-Spline surfaces into multiple pairs of parametric surfaces to perform surface intersection within these pairs.However, this conversion step introduces an additional computational overhead.To address this limitation, future studies should explore ways to exploit the geometric properties of B-Spline surfaces.By making full use of these properties, we can further accelerate our algorithm and improve its performance.Table 6.Comparison of the intersections of two B-Spline Surfaces.The first column shows the surfaces and their bi-degrees.The remaining columns represent the intersection results and the consuming times (in seconds) of different algorithms.AA-GPU refers to the affine arithmetic-based B-Spline surface intersection method with GPU acceleration.Mesh Booleans refers to the interactive and robust mesh boolean algorithm.OCCT refers to the intersection algorithm in the Open CASCADE Technology.SISL refers to the intersection algorithm in the SINTEF Spline Library.ACIS refers to the intersection algorithm in 3D ACIS Modeler.N/A indicates that the program reported an error in this example, and the 'NULL' indicates that the result obtained by this method is no intersection.4), and it obtains unsatisfactory result near cross point in example (3); besides, it returns errors in the cases that two surfaces intersect at their boundaries in examples (5)(6).ACIS can accurately calculate the intersections in most quadratic cases, but still miss the tangential line in example (1).Among all experiments, our method can guarantee the topology of intersections and calculate them within an acceptable time, which is comparable to the computation time of ACIS in correct cases.

Fig. 1 .
Fig.1.We proposed a topology guaranteed B-Spline surface intersection method that is robust toward various intersection topology, including (a) cross intersections, (b) isolated contacts, (c)(d) intersections along boundaries, (e) contact in different intersection branches, (f) high-order contact along a whole curve and (g) multiple isolated contacts.The intersection loci of our method are presented in red curves, while the results of ACIS and OCCT are presented in yellow and green curves respectively.'NULL' indicates that the result is an empty curve.The boxes with dashed boundary lines indicate wrong topology.Our method provides correct intersection topology in all cases, while ACIS and OCCT sometimes fail.

Fig. 3 .
Fig. 3.The OBB (a) and AABB (b) of a spline surface.The control points are marked in red.The OBB is tighter than the AABB.

Fig. 5 .
Fig. 5. Examples of different points on the same implicitized surface  (, , ,  ) = 0. Figure (a) presents points  1 and  2 in Example 1; figure (b) illustrates the detail of the relative position of the points and the surface.The figures show that the yellow point is within the required domain of parametric surface P(, ), while the red point lies outside the required domain of the parametric surface.

Fig. 7 .
Fig. 7. Intersection curve lifting from 2D parametric domain (s,t) to 3D space.Left is the intersection curve C with different local branches.Right is the lifted intersection curve C with different viewpoints.Each local branch in C is lifted from the local branch in C with the corresponding number.
The strategy consists of two parts: 1. Split the so far obtained C in Fig. 8(b) at its intersection points with the boundary curve of the surface P(, ), as is presented in Fig. 8(c) and Fig. 8(d); 2. Clip C and select the part within (, ) ∈ [ 0 ,  1 ] × [ 0 ,  1 ], as is presented in Fig. 8(e).The details are as follows.(1) Split C: We sample (, ) on the four boundary lines of surface P in [ 0 ,  1 ] × [ 0 ,  1 ] with a step size of 1/ , as is shown in Fig. 9(a), and then obtain a series of approximate intersection points p  ,  = 1, • • • ,  of C with these four boundary lines.These points are called splitting points, which are marked with red asterisks in Fig. 9(b).Split C at p  ,  = 1, • • • ,  into branches  1 ,  2 , • • • ,   .See Fig. 8(d).

Fig. 8 .
Fig. 8. Clip a 3D intersection curve to adapt the required parametric domain.(a) is the 3D intersection curve C before clipping; (b) presents different intersection branches in the parametric domain (,  ) by distinct colors; (c) illustrates the sample points of boundary curves of the surface P; (d) shows the splitting points in red asterisks; (e) is the selected branches after clipping; (f) depicts the final intersection curve of two surfaces in the required parametric domain.
Fig. 9. Curve clipping detail.(a) marks the sample points on boundary lines in magenta; (b) shows the splitting points with red asterisks, which are composed of the points on the intersection branches that are close enough to the sample points in (a); (c) records the selected points (black point and red triangle) on each split branch, which are used to determine whether the branch needs to be abandoned.The black ones are the branches to be remained while the red ones are to be abandoned.(d) represents the final selected branches.

Table 3 .
Intersections of torus and ellipsoid.The first column provides the parameters of the torus and the ellipsoid for each example; the second column displays two surfaces; the third column represents the intersection results; the remaining columns illustrate the time cost (in seconds) for both ACIS and our method.Given the absence of the ellipsoid class in ACIS, we represent the ellipsoids as piece-wise spline surfaces and intersect them with the torus.Bold fonts indicate better performance.

Fig. 10 .Fig. 11 .
Fig. 10.Two examples of the intersection curve of B-Spline surfaces in the parametric domain (s,t).The figure on the left represents example (1), while the figure on the right represents example (5).Example (1) has only one local branch while example (5) has several local branches.
the topology determination step of our algorithm, it is possible to encounter situations where the left polynomial     () and the right polynomial  ℎ () corresponding to  have no roots in the parametric domain  ∈ [ 0 ,  1 ], which indicates there are no left helping points and right helping points corresponding to the root .As is shown in Fig. 11, we categorize this situation into four distinct cases and discuss each of them individually.(a) If the roots of     () and  ℎ () are outside the parametric domain  ∈ [ 0 ,  1 ] but the curve traverses the region [ 0 ,  1 ] × [ 0 ,  1 ], we can get the curve through boundary point tracing.(b) If the roots of     () and  ℎ () are outside the parametric domain  ∈ [ 0 ,  1 ] and the curve does not fall within the parametric domain, we just ignore this curve.(c) If the curve is a vertical line  = , we can detect it by substituting  =  into the polynomial  (, ) in Eq. 12 and determining if  (, ) is a zero polynomial.(d) If the intersection is an isolated point, we can find the point by substituting  =  into the polynomial  (, ) and solve  (, ) = 0 in the parametric domain  ∈ [ 0 ,  1 ].To sum up, when     () and  ℎ () have no roots in the parametric domain  ∈ [ 0 ,  1 ], we substitute  =  into the polynomial  (, ) to get polynomial  (, ).If  (, ) is a zero polynomial, then we uniformly sample  =  in  ∈ [ 0 ,  1 ] to get the vertical intersection line in the parametric domain.If  (, ) is a nonzero polynomial, we solve  (, ) = 0 in  ∈ [ 0 ,  1 ] to get the possible isolated intersection point.Our algorithm also provides dedicated treatment for the intersection topology along the boundaries of the parametric domain.We first detect whether the boundary lines of the domain (, ) ∈ [ 0 ,  1 ] × [ 0 ,  1 ] are intersection curves of the two surfaces as follows.We substitute  =  0 into the polynomial  (, ) to get   0 ().If the   0 () is a zero polynomial, it indicates that the boundary line   0 := {(, ) ∈ R 2 |  =  0 ,  ∈ [ 0 ,  1 ]} corresponds to an intersection of the two surfaces.Similarly for the other boundary lines  =  1 ,  =  0 , and  =  1 .Then we uniformly sample the boundary intersection lines and map them to the 3D space to get the corresponding boundary intersection curves in R 3 .

Fig. 12 .Fig. 13 .
Fig. 12. Boolean operation results of two hummingbird models.The hummingbird model consists of 312 bi-cubic B-Spline patches, with each B-Spline patch containing 6×6 control points.Our method has detected 99 pairs of intersecting patches.Union, intersection and subtraction are performed after computing the intersection curves of two surfaces.
boundary intersections in example (5); additionally, it fails to address the tangential conditions in examples (1)(2).SISL gets a bunch of intersection points near tangential line and tangential point in examples (1)(