Intersection Algorithms (Geometric Modeling) Part 4

Surface Homotopy Method

The problem of finding intersections of implicit surfaces can be thought of as a special case of the more general problem of solving systems of polynomial equations. One much-studied technique that is used is called the homotopy method. It has also been called the homotopy continuation method or simply a continuation or embedding method. The idea is to find the solutions to a related set of simpler equations first and then to deform these equations and their solutions into the given ones. A simple example, presented in [PraG86], considers the problem of finding the solutions to the equations

We note that the following related equations

are easily solved and the solutions to the parameterized equations

define a homotopy between the solutions to equations (13.20) and (13.21). Therefore, to solve the equations (13.20) we compute the incremental changes to the solutions to (13.21) as t changes from 0 to 1. One surface intersection algorithm that uses the homotopy method can be found in [AbdY96].

Here is an overview of the general homotopy method. A good survey can be found in [AllG90] and [Wats86]. Some other helpful papers are [GarZ79], [Morg83], [Wrig85], and [AllG93], where one can also find many additional references.

A system of m polynomial equations in n variables corresponds to a polynomial mapand an equation f (x) = 0.    (13.23)

We choose another system of equations

defined by a mapwhose  zeros  are known  and  consider  the  homotopy  h: between g(x) and f(x) defined by

Let x0 be a zero of g(x). The object is to find a curve g(t) in the zero set of h that starts atand ends at a pointwhere   x1 is a zero of f. More precisely, we look for a curve

with the property that

One way to find g(t) is to differentiate the equation in condition (2). This shows that g(t) must satisfy

where J(x,t) is the Jacobian matrix

If the rank of J(x,t) is n for all values t, then the implicit function theorem will guarantee that g(t) exists. We are left with the problem of solving the system of ordinary differential equations defined by (13.26) and can use any of the well-known ways to solve such equations.

A second way to find g(t) is to use a standard root finding method like the Newton-Raphson method to the equation in condition (2) above.

We mention two potential problems with using the homotopy method. In general terms, what sets the homotopy method apart from previous iterative schemes is the fact that it replaces a “local convergence” approach with a “global convergence” one. With the homotopy method we already have some knowledge about certain initial zeros (if our initial guess for a zero is poor when using a standard Newton-Raphson method then we may have extreme difficulty in converging to a zero) and we can use them to iterate to all of the zeros of our function. Nevertheless, we are still dealing with an iterative process and one big problem with the type of iterative methods that are used is that when relevant Jacobian matrices have singularities, then it is very hard to guarantee the convergence to a desired solution. It took a lot of work to overcome the singular Jacobian matrix problem and make the homotopy method work and be reasonably efficient.

Another problem that caused great inefficiencies in the homotopy method initially was the choice of a start function g(x). The original approach was to use Bezout’s theorem and choose a polynomial function g(x) of total degree d that was the product of the degrees of the polynomials in f(x). In practice, however, the number of zeros of f(x) was often much smaller than d and so a lot of computation effort was wasted in generating curves g(t) that diverged toIt    was    important to find a better bound on the number of zeros of f(x). This can be done in the case of sparse polynomials, that is, polynomials that have a relatively small number of monomial terms. One gets a more efficient homotopy method when f(x) is a sparse polynomial. See, for example, the paper [VeVC94], which concentrates on sparse polynomial systems.

Surface Recursive Subdivision Methods

The general idea of using recursive subdivision (divide-and-conquer) in computer graphics is an old one. Display algorithms in early papers, such as [Catm74] and [LCWB80], made use of it. The first suggestion that it might be useful for intersection problems can be found in [LanR80]. Bezier and B-spline surfaces with their control points are good candidates for this approach because their subdivisions have associated control nets. In [LanR80] the idea is to subdivide the surfaces enough so that the individual pieces are essentially flat (see also [Clar79] and [LCWB80]). One then finds the intersection of the approximating planar patches to get a polygonal approximation to the actual intersection. Because the surfaces satisfy the convex hull property, one can use bounding boxes to make the algorithm more efficient. This is the analog of the recursive subdivision for curve intersections in Sections 13.3.1 and 13.3.4.

One problem associated to recursive subdivision is the amount of data that could be generated potentially. The way to mitigate this problem is by not doing the subdivision down to a certain level blindly. One should keep subdividing an object only at those places where there is a possibility of it intersecting the other object. The typical way that such an adaptive approach is carried out is by using some sort of bounding box approach. Min-max boxes, which are obtained from the minimum and maximum values of the coordinates of the points of an object, are common choices because they are so simple. Other bounding objects are slabs (see Section 6.2), convex hulls in the case of Bezier or B-spline surfaces, or “fat” planes ([Carl82]) for Bezier surfaces. If the bounding objects do not intersect, then the objects will not intersect. Checking for intersections of the bounding objects is much easier than checking for intersections between the objects themselves. One subdivides the objects only at those places where the bounding boxes intersect. This approach implies that one has a termination criterion. The usual such is based on a flatness test. There are many ways to check for flatness, but one simple one is to use the normal vectors. The normal vectors will not change much over regions that are close to flat. Of course, no matter how simplified a test one uses, this will still take some extra time and so the more one can avoid having to use such tests, the better. One might be able to say in advance how much subdivision is needed.

The algorithm in [HEFS85] is an example of a recursive subdivision algorithm for C1 parametric surfaces. It proceeds in four steps: subdivision, intersection, sorting, and refinement. The subdivision is done in stages. First, the surfaces are subdivided until the subpatches satisfy an initial flatness and edge linearity criterion. A stack of possibly intersecting pairs of subpatches is created. The test for intersection here is based on whether associated “oriented” bounding boxes (bounding boxes that roughly line up with their patches) intersect. The leftover subpatches are successively subdivided further to higher flatness tolerances. A tree data structure is used to maintain the data generated. When subpatches are flat enough, they are approximated by two triangles and the intersections of these triangles are used as approximations to the intersection of the patches. This produces a collection of segments that the sorting phase must then combine to get the polygonal curve segments that become the approximation of the intersection curves of the surfaces. Since the algorithm is based on adaptive subdivision there may be gaps between adjacent triangles, but the algorithm checks for this and redefines subpieces appropriately to prevent the problem from occurring. Finally, the triangle intersections are only approximations, therefore a refinement step iteratively tries to find points closer to the surfaces. We described a slightly improved version of this refinement in the Barnhill-Kersey algorithm in Section 13.5.2. An algorithm similar to the one in [HEFS85] is described in [FiMM86]. There bounding boxes are defined using derivative bounds.

On a related topic, [Nasr87] discusses finding intersections of recursive subdivision surfaces.

Surface Algebraic Methods

When one tries to use algebraic methods to solve intersection problems, one looks first for special cases that can be handled by special techniques. For example, plane/plane or plane/quadric intersections produce lines and conics, respectively, and can be solved for explicitly. In general, because the implicit/parametric case is relatively easy, attempts have been made to reduce other cases to this one.

An implicit surface S can be parameterized, but not necessarily by rational polynomial functions. If S is defined by linear or quadratic polynomials, then S can be parameterized by rational polynomial functions. If S is defined by higher-degree polynomials, then it may not admit such a parameterization. See Section 10.15 in [AgoM05] for more details. Parameterized surfaces, on the other hand, can always be represented implicitly by rational polynomial functions if the parameterization was also of that form. The only problem is that the standard implicitization algorithms may produce very complicated equations. See Sections 10.9 and 10.10 in [AgoM05]. For example, it can be shown that a bicubic patch is equivalent to an algebraic surface of degree 18 whose equation contains 1330 terms. For these reasons, algebraic geometry methods seem to be impractical currently, but there is a lot of ongoing research.

At any rate, because it is known that every algebraic curve in R3 can be mapped to an algebraic curve in R2 (although the latter may be more complicated than the former), one general algebraic approach to solving the surface intersection problem is:

(1) Map the intersection curve in R3 to a plane curve defined by an equation

(2)  Solve equation (13.23).

(3) Map the solution back to R3.

Step (2) is considered in more detail in Sections 14.5.1 and 14.6. Here we describe two approaches to (1) and (3). One is based on substitutions and the other, on projections. Let S1 and S2 be surfaces in R3.

The Substitution Approach. Suppose that surface S1 is defined implicitly by an equation

and surface S2 is defined via a parameterization

Substituting into equation (13.24) gives

If we solve equation (13.25) in the u-v plane, then we can map the solution back to R3 using g.

If both surfaces are defined parametrically, we can implicitize one of them using the Grobner basis approach and thus reduce the problem to the previous case. We could also use the resultant to implicitize a surface, but that method has the disadvantage of introducing extraneous factors. See Section 10.9 in [AgoM05].

If both surfaces are defined implicitly, we would like to parameterize one of them to again reduce the problem to the case solved above. Unfortunately, as was pointed out earlier, it is not always possible to parameterize an implicitly defined surface by rational functions. On the other hand, it can be shown that the intersection of two implicitly defined surfaces always lies on a parameterizable surface X. Here is an algorithm that will produce such a surface. Assume that S1 and S2 are the zero sets of functions f(x,y,z) and g(x,y,z), respectively.

Step 1. Homogenize the function f(x,y,z) and g(x,y,z) to get homogeneous functions F(x,y,z,w) and G(x,y,z,w), respectively. The intersection of the surfaces S1 and S2 corresponds to the nonideal points of the intersection of the homogeneous hypersurfaces defined by F and G.

Step 2. Choose one of the variables appearing in F or G and express F and G as polynomials in that variable (with coefficients that are polynomial in the other variables). If we suppose that w was chosen, then we write

where ai and bj are polynomials in x, y, and z. Assume without loss of generality that We can assume thatbecause otherwise f and g would be homogeneous polynomials and that special case will not be considered here. Define

We can think of F1 and G1 as having been derived from F and G by removing their highest, respectively, lowest degree terms. Note that both are linear combinations of F and G and hence the intersection of the hypersurfaces defined by them contains the intersection of the hypersurfaces defined by F and G.

Step 3. Since the degree of F1 and G1 is less than n, we repeat Step 2 with F1 and G1 replacing F and G, thereby generating a sequence of polynomials Fi and Gi, until we finally end up with an Fk or Gk, which has degree 1. (The case where Fi = Gi for some i and where we go from a degree larger than 1 to a degree 0 in one step is a special case not dealt with in our algorithm.) Using this linear polynomial in w we see that our (homogeneous) intersection lies on a hypersurface defined by an equation of the form

By induction, the polynomial H is a linear combination of F and G.

Step 4. The hypersurface defined by equation (13.28) can be parameterized (with homogeneous coordinates) by

Substituting these functions into the formula for G, we get a homogeneous plane curve. Dehomogenizing this curve gives us an affine plane curve that is now solved.

Example. Consider the ellipsoid S1 centered at the origin and the sphere S2 centered at (1,0,0) defined by

and

respectively. Figure 13.19 shows the x-z plane cross-section of the two surfaces. We show how to use Steps 1-4 above to map the intersection of S1 and S2 to a planar curve. (Of course, because the equations are so simple, we could have done this directly without following any fancy steps, but this is beside the point.)

Solution. Step 1 produces

Step 2, based on the variable x, leads to

Figure 13.19. The ellipsoids of Example 13.5.5.1.

Choosing the simpler equation F1 = 0 to solve for x, we get

Substituting this expression for x in G and dehomogenizing, that is, setting w to 1, gives us the equation

Equation (13.32) defines the circle

in the y-z plane. (To see this, letin    equation (13.32) and solve for u.)

Notice that we can run into serious problems if the wrong variable is chosen at Step 2. For example, choosing y we would get

and we are not able to solve for y. The reason that x worked and y does not is that the intersection projects nicely in a one-to-one fashion to the y-z plane using an orthographic projection parallel to the x-axis whereas projecting parallel to the y-axis collapses the intersection to a closed segment.

One problem with the above approach is that it has the potential to introduce extraneous factors and solutions. The intersection of X and S1 may be larger than the intersection of S1 and S2. See [Hoff89] for examples of this.

The Projection Approach. The idea here is to use a central projection from some point p to project the space curve to a plane curve. The only problem is choosing p correctly. We do not want the projection to introduce any singularities so that the map cannot be inverted. [Hoff89] describes a method that works with a high probability.

Step 1. Transform the surface equations by a linear transformation defined by a matrix with symbolic coefficients.

Step 2. Use the resultant to project the intersection.

Step 3. Substitute random values for the symbolic coefficients and check that the projection has the desired properties.

Example. We redo Example 13.5.5.1 using the projection approach.

Solution. We can skip Steps 1 and 3 this time because there are no singularities when resolving on x. Using the resultant which eliminates x from f and g, we get

This again defines a circle in the y-z plane. If we had tried to eliminate the y variable we would have a problem because

We need Step 1 and 3 in general to make sure that there are no singularities in the projection. We were lucky when we chose x originally.

A major problem with algebraic approaches is that algorithms for finding solutions to high-degree polynomial equations are numerically unstable. It is also computationally expensive to compute resultants. By combining elimination theory with matrix computations the authors of [ManD94] tried to avoid these problems. Rather than finding roots of polynomials they had to find eigenvalues of matrices. Algorithms for finding eigenvalues are known to be stable. The algorithm in [ManK97] builds on this approach. Because one only needs the eigenvalues in a certain range, the new algorithm saves time by not computing those that are not needed.

Summary

In this topic we have seen some of the difficulties involved in finding the intersection of objects. For that reason, many special cases have been studied extensively to get the best possible results. We have looked at some of these. Intersections of lines, conics, and quadrics are other special cases about which much is known. Here are a few of those and references to optimized algorithms for them:

Rectangular solids and convex polyhedra:	[Gree94]
Line and cylinder:	[Shen94]
Ray and cylinder:	[CycW94]
Line and cone:	[Shen95]
Ray and quadric surface:	[CycW92]
Parametric cubics:	[Klas94]

Some more references for intersections of quadrics are [OckS84], [GolM87], [Mill87], [OweR87], [SheJ87], [FaNO89], and [Gold83] (for quadrics of revolution). The special case of ruled surfaces is considered in [HeKE99]. [HMPY97] describes a robust interval analysis approach. [FaNO89] studies a class of degenerate quadric intersections that are rather common cases. Because the intersection curves in these cases admitted rational parameterizations, algebraic methods could be applied. Specifically, the approach described in [FaNO89] made use of both the Segre characteristic of a quadratic form (which is something determined from the multiplicities of roots to an associated polynomial) and the projecting cone of a space curve (which is the ruled surface determined by the pencil of lines through the origin and points on the curve) with multivariate polynomial factorization algorithms. An equivalent approach was described in [WilM93] that was somewhat simpler and facilitated the display of the intersection curves.

For cyclide intersections see [MaPS86] and [John93]. For formulas for the intersection of a ray with various objects see [Hain89] and [Hanr89]. More references for the general intersection problem can be found in [AbdY97] and [HosL93]. [MarM89] and [LuMM95] address the difficult tangential surface-surface intersection problem.

Note that nothing has been said about set operations on CSG objects. The reason is that there is nothing new to say in that regard. The fact is that if we apply set operations to two CSG objects, then we get another CSG object, so that the problem is dealt with in the context of the boundary evaluation algorithm for CSG. See Section 5.9.

Methods, such as marching methods, which use Newton-Raphson iteration somewhere along the line, are the most common. [BarK90] has a comparison of some marching methods. [DoSY89] compares marching and recursive subdivision methods. Pure marching techniques tend to be faster, because the Newton-Raphson method has a quadratic rate of convergence, but are liable to failure because they are very sensitive to local singularities. Pure recursive subdivision techniques need no starting points and tend to be more robust but are less efficient. They can produce an excessive amount of data for a fixed tolerance when singularities are present. [DoSY89] describes a method that is a combination of the two where one uses recursion to discover the geometry and iteration for accuracy. [AzBB90] compares the two methods for Bezier surfaces and shows that iteration is more accurate. An approach that combines marching and algebraic methods can be found in [KriM97]. It uses a matrix representation for the intersection curve.

The homotopy method described in Section 13.5.3 has not been used much in geometric modeling. Marching methods can actually be considered to be a kind of special case. The papers [LamM95] and [LamM96] compare the two methods, and the authors argue that the homotopy method has advantages over methods based on standard Newton-Raphson iteration. It also needs a “starting solution” but often has better convergence properties. Basically, the problem is that the basins of attraction for the Newton-Raphson method (the points to which one converges) tend to be fractal-like whereas the corresponding sets for homotopy methods are smoother. They are semialgebraic sets when one is solving algebraic systems. It is the chaotic nature of the Newton-Raphson method that causes its problems. The negative aspect of homotopy methods is that according to [LamM96] they are roughly 10-20 times slower than Newton-Raphson methods (assuming that the latter converge). The paper describes how homotopy curves are followed.

As mentioned above, getting the wrong connectivity of the intersection curve is a big problem. Three aspects of this problem are loop detection, jumping components, and incorrectly ordered components.

Loop Detection. The reason why it is good to have criteria for the existence of loops in the intersection is that if there are none, then the intersection must start in the boundary of the surfaces and one could then start with the simpler problem of intersecting a curve with a surface. Unfortunately, the flatness conditions in subdivision algorithms have a hard time dealing with loops. For example, consider the case of a large sphere and a plane. As the sphere approaches the plane, one passes from no intersection to the case where the intersection is a point and finally to where it is a circle. Locally the sphere is very flat and it becomes very hard to detect the circles when they are still small.

Jumping Components. In the case of tracing algorithms, one is trying to generate sequences of points, each of which describe one component of the intersection curve. It is not easy choosing a step size that will prevent the algorithm from jumping from one component of the intersection curve to another if they are very close. This would produce the wrong connectivity in the result.

Incorrectly Ordered Component. Similarly, one might accidentally jump from one part of the curve to another part of the same component if those parts passed close by each other. This would generate an incorrectly ordered sequence of points for that component.

[Hohm91] describes a loop detection criterion that can also be used to guarantee that a tracing algorithm does not jump components or generate points out of order. The criterion is based on properties of the Gauss map for the two surfaces. It isolates intersections with no loops quite well if the surfaces are relatively flat. An approach in case the loop criterion is not satisfied can be found in [KriM97]. [SedM88] describes a criterion that ensures that all branches of the intersection curve will be detected. An improved criterion is proved in [SeCK89]. [MaLe98] has an overview of approaches for loop detection with additional references and gives a topological criterion for the existence of loops in the intersection of two C2 parametric surfaces.

[Wang92] also describes a way to deal with the topological inconsistency problems.

Finally, in this topic we have been satisfied with simply getting the curve that is the intersection of two surfaces. In most algorithms one gets a polygonal approximation. There are times when one might want to get a higher-order approximation to the intersection. In that case it would be helpful to know the derivatives of the curve up to some arbitrary order k. Formulas for these and other geometric invariants such as tangents, curvature, and torsion, are derived in [YeMa99].