Modeling: Part

The next generation geometric modeling systems have higher requirements on the reliability and the efficiency of Boolean operations. The two parts of these technical re ports describe the details of the nonlinear intersection algorithms developed at EDRC. The prototype of the algorithms has been shown to be robust, accurate, and efficient, without the overhead of subdivision. The second part of this report describes the intersection routines for two dimen sional entities.


Introduction
In the second part of this technical report, we will describe the intersection algorithms involving two dimensional entities (i.e., surfaces and planes). These problems are more complicated (than the ones described in part I) because the search spaces are two dimensional and more sophisticated algorithms are required to find the intersection. In this report, we will limit the discussion to transversal intersection [3] only. In a nut shell, transversal intersection excludes the cases of tangential contacts, which mathematically correspond higher order singularities. Non-transversal intersection is still under research and will be discussed in a later report. The (transversal) intersection in SPI and SSI can be classified into two categories: visible and invisible curves (see figure 1). For the former, the intersection curve starts from and/or ends at the boundary of the domain. In other words, the intersection curve is visible from the boundary and can be obtained once the boundary intersection point is found. For the latter, the intersection curve is totally embedded in the domain, thus invisible from the domain boundary. In the literature, this problem is sometimes referred as loop detection ( [14], [12]).
As mentioned in the part one of this report, here we only discuss the intersection algorithms for untrimmed entities. For example, a planar face, which can be a nonuniform multiply-connected polygon, is generalized to an infinite plane for the SPI operation. All the free-form surfaces are assumed to have the unit square as its domain. Post processing will be employed to adjust the intersection result to the appropriate domain. A more efficient way is to consider the effect of the trimmed domain in the preprocessing. The algorithms described in this paper can be easily modified if desired.

2*1 General Description
The problem is to find the intersection between a parametric surface and an infinite plane, ax + by + cz + d = 0. From the convex hull property of B-spline surfaces, we can easily exclude the non-intersecting case by inspecting the location of the control points of the surface. This inspection can be further simplified by applying a transformation to the surface and the plane such that the plane equation is equal to z = 0. As will be shown later, this simplification also unifies the treatment of rational and non-rational surfaces.

2.2.1
Transforming a plane to the nominal position (z=0): This operation contains two transformations: translating the plane such that it passes the origin, and orienting the plane such that the normal is aligned with the Z-axis. This is exactly the same operation as in CPI. Please refer to part I of the report. Note that if the surface is rational, in which case W(u,v) is no longer a constant function, the formulation is still the same. That is, Z-component is the only effective term.

Derivation of the intersection curve f(u,v)
From the above, we know that SPI has a close form solution f(u^v) = 0. However, in geometric modeling, it is more important to know the trace of the intersection curve than the exact algebraic representation. Since this is a problem with one equation and two variables, continuation method [11] can be applied to trace the solution set /(ti,v) = 0.
The evaluations of function, Jacobian, and Hessian of the function f(u,v) come from the Z-components of surface evaluation, first derivatives, and second derivatives respectively. However, there are some important issues: • There may be more than one components in f(u,v) = 0. How to make sure that enough starting points have been located such that no component is missing.
• How to realize whether a starting point has been traced by other branches such that no repetitive effort is spent on that point.
For the former issue, our approach can find all the necessary points of the visible curves, thanks to the robust CPI solver. Starting points for invisible curves are also found using a robust algorithm, which will be discussed later (section 2.2.4). For the second issue, the intersection curves are separated at the characteristic points. The bookkeeping of tracing record at these points is the key of our method.

Numerical Computation of Characteristic Points:
Since the intersection branches start from characteristic points, accurate computation of the characteristic points is crucial to curve tracing. The accuracy is achieved by numerical computation at these points. In the following, we will describe the definition and the formulation used in Newton iteration at each type of points.
border points : the intersection points at the boundary of the surface. Numerical computation of these points is actually an one-dimensional problem. These points are calculated accurately using the CPI routine.
u-turning points : A turning point is a point at which the tangent vector is parallel to either axis in the parametric space. U-turning points (with tangent parallel to v-axis) are the key to our curve tracing processing (see section 2.2.4). Vturning points (with tangent parallel to u-axis) are only useful in the case of the complication of border points (section 4). The following two equations are satisfied at a u-turning point: With two variables and two equations, we can compute a turning point accurately using the Newton iteration, provided that a good initial guess is given. Note that if we use Newton-Raphson iteration, the function, and therefore the surface is assumed to be at least twice differentiate (C 2 ).
Continuity in general will not be an issue, unless the point happens to be exactly located at the boundary of a component subpatch where the continuity condition is not guaranteed. In that case, we need to break the surface into its component patches such that the continuity condition is satisfied l . ular points (nodes) : Nodes are the points at which the intersection curve crosses itself. At these points, three conditions are satisfied: / = f u = f v = 0. We formulate the problem to be the gradient map of the objective function f 2 + fu 2 + fv 2 -That is, In this formulation, we have the same concern of continuity. Actually, this formulation requires one order higher of continuity (C 3 ).
^ote that in general continuity should also be a concern for curve tracing. However, in continuation method, the requirement is only C 1 (Jacobian is needed for correctors in curve tracing). Based on the properties of B-spline surface, we can have a simple preprocessing routine to detect the tangential discontinuous / fu + fv, fuu + fv fuv -/ fv "t" fu fuv 4* fv fvv = features.

2.2.4
Finding all U-turning points: As was previous mentioned, intersection curves are categorized into visible and invisible. Visible curves are easily discovered by solving the border points. For invisible curves, we utilize the following simple fact: On every invisible curve, there is at least one u-turntng point or one singular point.
Thus, the solution strategy is to find all the u-turning points and singular points. That is, find all the points that satisfy / = f v = 0. The nature of this problem is to find a methodology to solve all the roots in a system of equations. We have reviewed several methods: homotopy, interval arithmetic, and trajectory methods ([2], [5], [7], [9]). In the current implementation, we have experimented the following package: CONSOL developed by Dr. Alexander Morgan in General Motors Research Laboratory [9]. The package uses the homotopy concept and perfects its implementation in polynomial system. As is well known, finding the appropriate homotopy formulation is the most difficult problem. He has developed a systematic way of finding the homotopy equation that will solve the polynomial system.
To use this package, we need to formulate the problem as a polynomial system. If the surface is not a Bezier surface, we need to subdivide the surface into Bezier subpatches, which are by definition polynomials. Although the subdivision might create a lot more subpatches to be considered, many subpatches can be excluded because either they have no intersection with the plane, or it is possible to have invisible curves in them. The time and storage spent on the overhead of subdivision is moderate, but the robustness is greatly enhanced.
The system is formulated as follows: This is a two-variable, two-equation problem. The equations are of the form: = 0 = 0 By expanding the first one into monomial form 2 , the first equation can be express and K is (TO + 1) x (n + 1) coefficient matrix.
The efficiency of the homotopy method in this case is directly proportional to the number of paths explored, which is proportional to the total degree of the system. The following preprocessing routines are employed to reduce the degree, and thus to increase the efficiency.
• Delete the terms that are negligible because of small coefficients.
• Factor out any variable in each equation. Make su:-e that in every equation, the term of lowest degree in each variable is a constant.
• Introduce a new variable x$ = uv, if this substitution can reduce the total degree. Notice that this substitution will increase the system size to three, but at this point, we assume that the total is the dominant factor of execution time. 3 Interpretation of the results: CONSOL will calculate all the roots (real and imaginary) in the system. In general, we are only interested in the real roots that are inside the domain. However, the tolerance used forjudging whether a root is real might vary depending on the type of the solution. For a regular simple solution (e.g. the u-turning point in a circle), the imaginary part is quite small (10~1 5 ). As the degree of singularity gets larger (e.g, nodes and cusps), the imaginary part is also larger. It might be as large as 10~3. Thus, the choice of the threshold of real roots needs to take the type of roots into consideration.
Infinite number of roots: CONSOL is designed to solve the system with (geometrically) isolated roots. Thus, we need to make an effort to exclude the case of 2 We found that the numerical precision is not a big problem for this expansion if we use double precision arithmetics. 3 Although this is the conclusion we drew from several experiments, this conjecture needs to confirmed.

Figure 3: It is impossible to have a hidden straight segment in SPI.
non-isolated roots as best as we can 4 . In the current context, this effort is translated to extracting the presence of every non-isolated u-turning points. We claim that this task can be achieved if we can extract all the vertical straight component in the solution set. Furthermore, the extraction can be done by checking the location of the border points. Here are the explanations: • The intersection between a Bezier surface and a plane is an algebraic curve. Thus, it is impossible to have a straight vertical segment connected to some curve branches, as shown in figure 3.

• If we ever discover some u-turning point belonging to part of vertical intersection line (with zero curvature), this point must belong to a visible vertical component of the intersection curve. And this intersection component must be of the form u = ti*, where u* is a border point we discovered in CPI.
The extraction is implemented numerically. Our experience showed that with double precision arithmetic, it is quite reliable.

2.2.5
More on curve tracing: There are two issues that we need to consider. First, usually we get more than enough starting points to complete all components of the intersection curve. In order to avoid tracing the same branch many times, we need to do some bookkeeping on the starting points (i.e., the characteristic points). That is, we need to turn off the appropriate / ^ Divergent Corrector Inappropriate predictor size For nodes, the choice is determined by comparing the incoming direction of the branch with the four branches emanating from the node. The choice for cusps is the most difficult. This will be discussed in a later section. Second, there are some features that might hinder the efficiency of curve tracing. The curve tracing process employs the predictor-corrector method to follow the curve. If we use a big predictor to start off a small loop, the corrector step in continuation might fail to converge. In that case, we need to reduce che size of the predictor until the corrector behaves satisfactorily. This iteration process can be costly if we are dealing with high order surfaces. A better solution is to start off with the predictor of appropriate size determined from on the local geometry, e.g. the magnitude of the curvature at that point.

Current Development
In addition to CONSOL , we are experimenting another root solution package using interval arithmetic, GENBIS [7]. The initial experiments show that GENBIS is in general a lot more efficient. However, as reported in [7], GENBIS is poor in handling singular solutions. In our experience, for the configurations with nodes and cusps, the solution process will not stop until reaching the prescribed number of evaluations. Moreover, the robustness (in terms of the capability of solving every root) is sometimes worse than CONSOL . Since interval arithmetic method is used, it is especially vulnerable to the cases where there might be infinite number of roots. Thus, extracting non-isolated roots becomes crucial in applying this method. Because the efficiency gain of using GENBIS for the cases that worked is quite big, we are contemplating writing a high level decision routine to choose the appropriate solution methodology.

An Application -Slicing
Slicing is an operation to calculate a set of cross sections of a solid object. It is useful in many applications. For example, these cross sections are used in rapid prototyping technologies for constructing the layers of the model. The criteria of a good slicing algorithm are twofold: The cross sections calculated must be accurate, and the operation must be reasonably efficient. Most of the parts that we want to deal with are described in free form geometry. We have two options of handling this kind of object: linearize the part (i.e., approximate the free form surface by linear facets) and employ linear intersection algorithms, or calculate the intersection directly for the free form geometry. For the former, the efficiency of intersection algorithm is in general better, but the accuracy depends on the approximation of the linear facets. To generate a good linear approximation of a free form object is quite costly, both in terms of creating and handling of the data. Thus, we have decided to tackle the problem by the latter option. With the SPI algorithm described in the previous section, we believe we can accomplish an accurate slicing result with reasonable efficiency.
There are two steps to generate a slice : Calculate the intersection curves between the plane and all the surfaces in the model, and then connect the intersection curves to form a loop 5 . In the following, we will discuss only the computation of intersection curves. Since we have to intersect the same surface with a set of parallel planes, we need to make some modification to the general SPI algorithm. Without the loss of generality, we assume that the plane and the surface have been transformed such that the plane is parallel to XY plane. That is, the cutting planes are of the equations z = c t , where c t is a distinct constant for each plane.
Plugging in the component functions of the surface r(u, v), we get: The critical points are the points at which the following equations are satisfied: f u = f v = 0. It is clear that within each invisible loop, there must be at least one critical point . If r(w,v) is a non-rational surface, W(u, v) is a constant. Thus, the critical points for every slice will be the same (independent of c t ) . Thus, the timeconsuming critical point calculation will only be necessary for the first slice. The slicing operation can be described in the following: • Solve the visible curves using the same routines as in SPI. For rational surfaces where the homogeneous coordinate term W(u,v) is a function of u and v, the critical points for each slice will be different. The general SPI algorithm has to be used for each slice.

General Description
There are three requirements for a surface-surface intersection (SSI) algorithm: accuracy, efficiency, and robustness. Accuracy is measured by the deviation between the derived intersection curve and the actual curve lying on the participating surfaces. Efficiency is measured by the time and space needed for the operation. A robust algorithm should correctly identify every segment of the intersection curve. In the following, we review several types of SSI algorithms, and compare their achievements in these three requirements.
• Analytic approaches seek analytic expressions for the intersection of two algebraic surfaces [13]. Although it can reliably provide a closed form representation of the intersection, the complexity of the representation makes it impractical for the application of moderately high order parametric polynomial patches.

• Subdivision methods [6] decompose the surfaces into polygons or triangles and an approximation of the intersection is calculated from intersections of the linear facets. These methods may fail to compute all loops and branches for the intersections, due to the finite level of subdivision. Further increase of the subdivision level does not assure correct topology of the intersection, and can rather result in excessive computation. One advantage of the subdivision method is that bounding boxes for subpatches can be built on the basis of control polyhedra [15].
Further subdivision of subpatches can be eliminated if the bounding boxes do not intersect.

• The basic idea of the lattice method is to reduce the dimensions of the operand of the intersection ([1], [16]). A set of spatial curves lying on one surface (e.g., isoparametric curves) is defined and their intersection points with the other sur-
face can be computed. These points are then sorted and connected to form an approximation of the intersection curve. The SSI problem is transformed into a series of curve-surface intersections, for which more robust algorithms are available. However, because of the discrete nature of the method, it is difficult to capture all small loops and singularities.

• By formulating the problem as a system of nonlinear equations, one can employ marching methods to obtain the intersection curve [4]. The critical issue is to provide all necessary starting points for every branch of the intersection curves. This is commonly done by intersecting linear approximations of the two surfaces.
These methods are usually sensitive to singularities. Small numerical perturbation might result in the change of topology of the intersection curve. Nevertheless, since every point on the curve can be computed accurately by numerical techniques, marching methods are more accurate.
In this paper, we will describe another method we developed -the Augmented Lattice method. The basic idea is to discover all invisible intersection curves by locating the critical points of the distance function between the two surfaces. Since the distance function is quite complex, the powerful techniques used in SPI can no longer be applied. The strategy we use for exploring the nonlinear distance function is based on our previous work [17]. The lattice evaluation is added to increase the robustness. Also, the path tracker is changed to a more accurate continuation solver [11]. In summary, the solution procedures for solving SSI is the following: 1. Locate all visible curves -solve the CSI subproblems. Since the operands are assumed to be untrimmed surfaces, eight CSI problems need to be solved. Bounding boxes are used to increase the efficiency of the operation.

Augmented Lattice:
Step 1 : Generate a set of sampling lines on one surface. Explore the distance function along these sampling lines. Record the positions where the gradient vectors (of the distance function) are perpendicular to the sampling lines. The sampling lines are the isoparametric lines at the node position of the more complicated surface.
Step 2 : For each interval between the sampling lines, trace the U-characteristic lines. The goal of this tracing is to find the u-turning points and singular points of the intersection curve.

Trace the invisible curves starting from the new u-turning points.
Comparing this approach with the criteria of a good SSI algorithm, we conclude: Accuracy It is important to make sure that every intersection point is accurately computed, and that the deviation of the interpolation segment is within the tolerance. In the current approach, the interpolation points are calculated using continuation method. Thus, the accuracy is guaranteed to be satisfied (unlike the ODE tracing method). The deviation can be estimated using heuristic rules.
Efficiency The performance of our prototype implementation is similar to the data presented in [8]. It is our opinion that the execution time is not a good comparison, because the performance depends on the actual implementation. However, the algorithm does not use subdivision, which is always a benefit in saving the extra calculation and storage. Moreover, since the CSI subproblems (the border point computation and the sampling lines exploration) are done independently, this algorithm is suitable for parallel implementation. [8]. Moreover, because of the modification of additional sampling lines, the current implementation is more robust than our earlier work.

Robustness We have successfully tested the implementations of the previous [17] and the current works over all the cases listed in
Misc. Since we work directly on the distance function, it can be applied to the intersection of offset surfaces, without much modification.
The details of the implementation will appear in a later report. Here we only discuss the technical details of this problem.
This is a more complicated way of formulating a three-variable, four-equation problem (as compared to formulation 1). In this formulation, the parameters (u, v,s,t) are regarded as independent variables. Note that eq(ll) gives the notion of distance function as described in the second formulation, and that eqs. (12) and (13) are the same orthogonal conditions as in the second formulation.

Continuation Formulation:
Among the formulations described above, it is important to choose the appropriate one for different context.

• For tracing SSI curves, it is advantageous to choose the least expensive formulation, equation (5).
• In the exploration of the distance functions, the second and third formulations are more appropriate.

Numerical Computation of Characteristic Points:
As described in the SPI section, characteristic points are used to facilitate the bookkeeping of the curve tracing process. Notice that we only need to define the characteristic points on the first of the two operands in the Boolean operation, since the connectivity information can be processed with just one set of points. In the following, we will describe the numerical formulations for different types of points in the SSI context.

u-turning points : The points on the intersection curve with vertical tangent vectors.
From the second formulation (6), the following equations are satisfied at a uturning point: subjected to (p-q)-q. = 0 (16) (r -q) • q t = 0 (17) Note that nested Newton iterations is needed in this case. That is, the outer iteration is designed to converge to a u-turning point by adjusting the (ix,i?) (eqs (14) and (15)), while the inner iteration is designed to maintain the orthogonal conditions, (eqs (16) and (17)). Another set of simpler but equivalent conditions can be derived using the third formulation (eqs(ll) to (13)):

T(U,V)
-q(s,t) = 0 (18) n(j,t).r tt (ti,t;) = 0 Here the problem is a system of four variables. No nested Newton is needed but the size of the system is larger than the previous one. Note that equation (19) , is designed to converge to a u-turning point (see eqs. 15 and 10), and the first three (19) will guarantee to stay on the intersection curve (thus satisfying eqs (15), (16) and (17)).

Singular points (nodes)
There is one more condition to be satisfied than the uturning point, <f>u = 0. This extra equation will make the system over determined.
In the SPI case, we formulate the system to be the gradient map of a functional: f 2 + fu + fv-I n the current case, we can get away in a cheaper way by some not-very stringent assumptions.
We know the following equations are satisfied at a node:  To determine whether there are valid branches out of the turning point, we employ the local parameterization and the curvature information as follows.
We will explain the case of V-turning points in the following. The results of Uturning points can be derived in a similar fashion. In a planar parametric curve, the principle normal direction is defined to be 90 degrees counterclockwise from the tangent direction [3]. In the current case, there is no particular direction associated with the intersection curve. We choose, without loss of generality, that the curve goes in increasing u direction, as shown in figure 7.
The curvature for a planar parameterized curve 7(0( = ^(0>^(0) * s k(t) = uv -uv It is sufficient to inspect the sign of curvature only. Thus,

Sign(k) = Sign(uv -uv)
Since we take the local parameterization as t = u, and ^ = l,i = l,ii = 0. Thus,

Sign(k) = Sign(v)
Along the intersection curve,

5.2.1
Calculate the branch points: To start the curve tracing from a cusp, we need to obtain the starting points on every branch of the cusp. We develop a numerical method to calculate the branch point of a cusp as follows: • Perturb along the common tangent of the cusp. Unless the cusp point is on the boundary of the domain, this perturbation is performed on both orientations of the tangent (the possibility of being a double cusp). It is theoretically difficult to distinguish a double cusp from a single cusp.
• For each orientation, starting from the perturbed point, search for the branches on both sides of the tangent along a direction perpendicular to the tangent. The branch points are the points that satisfy f(u,v) = 0, where f(u,v), as defined in previous sections, is either the Z-component functions in SPI, or the distance function <f> in SSI. This is a one dimensional search problem, since u and v are related by the search direction. The search implementation is similar to the one in CPI (see part I of the report). Note that the search is terminated once two branches are found.
The success of the implementation depends on the numerical values chosen for the perturbation along the tangent, and the step size used in the one dimensional search. The perturbation has to be significantly large such that the branches are numerically distinguishable. A large step size in the one dimensional search is prone to missing the branch point, especially for the case of cusp of second kind.

5.2.2
Determine the incoming branch: As mentioned in section 2.2.5, in order to avoid duplicate tracing, we need to do the bookkeeping on starting points. The key to this bookkeeping is the comparison between 20 F(u,v) -epsilon = 0 F(u,v) = 0 A point whose tangent is perpenticular to the common tangent of the cusp Figure 10: Determine the incoming branches of a cusp the incoming direction and the branch direction. However, this is not applicable for single cusps, since both branches of a single cusp have the same tangent orientation 7 . Initial investigation suggests that the information of the curvature can be used to identify the branches for the cusp, but the usage is limited to cusps of the first kind. A more robust method is developed in the following such that the handling of both kinds of cusp is uniform. The basic idea to handle the branch classification is the following fact: For every cusp, there exists a separation line that passes through the cusp point, and separates the branches on both sides of the line.
Let us discuss the concept more carefully. The intersection curve that contains the cusp is an instance of the level set of the function f(u y v) = 0. Observe the trace of the perturbed function, that is, the curve of the function f{u,v) -6 = 0. It is always possible to find a point whose tangent is perpendicular to the common tangent of the cusp in the neighborhood of the cusp, as shown in figure 10. The branch on one side of the separation line differs from that of the other side by a sign, which denotes the orientation of the vector formed by the cross product of the tangent of the cusp and the local gradient vector. Thus, the sign for each branch can be determined by the branch points calculated in the last section. The incoming direction is determined by the sign of the last points in the curve tracing.

Cusp Branch Handling:
In section 2.2.5, we have mentioned the issue of branch management. The key idea is to avoid duplicate tracing by turning off the starting direction at the end of the branch. This method poses a problem for the cusps. The issue is that cusps are less desirable points to start off the curve tracing. Because the branches from a cusp are close, any over-sized predictor might go to another branch. On the other hand, to determine the incoming direction of a cusp is difficult. Although we have developed a robust method However, we can use this to distinguish the orientation of a double cusp.

21
in the last section, the continuation points might have already wandered to the other branch before we actually do this testing. There are remedies for both problems: If we want to start from a cusp, we need to pay extra caution on the size of the predictor at the beginning, until we are reasonably far away from the singularity. After we have decided the branch is to be terminated at a cusp, we need to check the incoming continuation points and decide which is actually the incoming branch that most of the points are following.