Intersection Algorithms (Geometric Modeling) Part 2

A Curve Newton-Raphson Method

To find the intersection of two smooth curves parameterized by functions p(u) and q(u) amounts to solving an equation of the form

tmp973d-35_thumb[2]

We describe an approach discussed in [HosL93]. Lettmp973d-36_thumb[2]

tmp973d-37_thumb[2]

(1)    Pick a pointtmp973d-38_thumb[2]on the curve p(u).


(2)    Find the intersection of the tangent linetmp973d-39_thumb[2]with the curve q(u).

Since this tangent line is defined bytmp973d-40_thumb[2]and    has    equation

tmp973d-46_thumb[2]

the equation we need to solve is

tmp973d-47_thumb[2]

 

or

tmp973d-48_thumb[2]

 

 

Curve intersections with the Newton-Raphson method.

Figure 13.4. Curve intersections with the Newton-Raphson method.

Lettmp973d-50_thumb[2]be a solution. Now find the intersection of the tangent linetmp973d-51_thumb[2]at

tmp973d-52_thumb[2]with the curve p(u). Let this point betmp973d-53_thumb[2]

(3) Continue this process, getting sequencestmp973d-54_thumb[2]Stop    whentmp973d-55_thumb[2]is sufficiently small.

See Figure 13.4(a). The method works as long as the tangents keep intersecting the curves. In Figure 13.4(b) this corresponds to starting at points along the boundary of the shaded region.

Numerical techniques for tracing curves work pretty well if one does not run into any singularities. The problem at singularities is that the approximation that one is using for the curve fails to be sufficiently accurate at such points.

Curve Recursive Subdivision Methods

We can extend the line-curve recursive subdivision (divide-and-conquer) approaches described earlier. The only difference is that we keep subdividing both curves and checking for intersections of bounding boxes. For curves that satisfy the variation diminishing property we get the algorithm from [LanR80]. The algorithm in [KopM83] applies to arbitrary curves.

The paper [SedN90] describes a better algorithm for planar Bezier curves. The authors refer to their approach as Bizier clipping. It is also a bounding box type approach.

Definition. A fat line is the region between two parallel lines.

Note that a fat line can be thought of as special case of the slabs defined in Section 6.2. The reader should compare the computations we will make for fat lines with those we made for slabs. Some other related concepts are the “fat arcs” in [SeWZ89], used as a criterion for curve flatness, and what could be called “fat planes”, used in [Carl82] for an intersection algorithm for Bezier surfaces.

The first thing we want to do is associate a fat line to an arbitrary Bezier curve p(u). Lettmp973d-56_thumb[2]be the control points of p(u). See Figure 13.5. Let

L be the line throughtmp973d-57_thumb[2]Represent    the    line    by    an    equation    of the form

The fat line for a Bezier curve

Figure 13.5. The fat line for a Bezier curve.

tmp973d-67_thumb[2]

and define

tmp973d-68_thumb[2]

It is easy to see that d(x,y) is just the signed distance from a point (x,y) to L because n = (a,b) is a unit normal vector fortmp973d-69_thumb[2]and

set

tmp973d-71_thumb[2]tmp973d-72_thumb[2]

Quadratic Curve.

tmp973d-73_thumb[2]

Cubic Curve. Although precise values can be computed in this case also, the formulas are complicated and the extra computational effort is not justified. Instead the following approximations are suggested in [SedN90]:

If ■tmp973d-74_thumb[2]then    let

tmp973d-76_thumb[2]

The convex hull property of B6zier curves implies that the curve is contained in the fat line defined by the two lines Li and L2, which are parallel to L and a distance of dmin and dmax from L, respectively. Actually, as the example in Figure 13.5 indicates, we can do better. The curve is contained in a thinner fat line and the lines Li and L2 can be moved closer to L. The thinnest fat line parallel to L is really determined by the minimum and maximum of d(u), the signed distance from p(u) to L. The polynomial d(u) can be explicitly computed in the quadratic and cubic curve case and one can use the following improved values for dmjn and dmax: tmp973d-77_thumb[2]

Sederberg-Nishita Bezier clipping.

 Sederberg-Nishita Bezier clipping.

Figure 13.6. Sederberg-Nishita Bezier clipping.

If tmp973d-79_thumb[2] then let

tmp973d-80_thumb[2]

After these preliminaries, we are ready to describe the Bezier clipping algorithm. Let p(u) and q(u) be the two Bezier curves with domain [0,1] whose intersection we want to determine. Let X be a fat line that contains the curve q(u). See Figure 13.6(a). In that example, X is the region between the parallel lines L1 and L2. We want to determine the values of u for which p(u) lies outside of X. Express p(u) in the form

tmp973d-81_thumb[2]

wheretmp973d-82_thumb[2]are the Bernstein polynomialsandtmp973d-83_thumb[2]are the control points of  p(u). Let equation (13.8) be the equation of the line throughtmp973d-84_thumb[2]and let d(u) be the signed distance from p(u) to L. It is easy to show that

tmp973d-88_thumb[2]

Finally, define the parametric curve

tmp973d-89_thumb[2]

wheretmp973d-90_thumb[2]The    fact    that    the    x-coordinate    of    D(u) is u  follows  from the identities

tmp973d-92_thumb[2]

Figure 13.6(b) shows the graph of D(u) for the example in Figure 13.6(a). Clearly, if

tmp973d-93_thumb[2] then p(u) will not lie on the curve q(u). Therefore one can trim away those parts of p(u). Actually, we shall only trim away parts at the ends of the curve. Lettmp973d-94_thumb[2]and tmp973d-95_thumb[2]be the largest and smallest u so that the convex hull of the D; misses the fat line Y betweentmp973d-96_thumb[2]over the intervalstmp973d-97_thumb[2]respectively.

For the example in Figure 13.6(b) it turned out thattmp973d-98_thumb[2]We

now subdivide the Bezier curvetmp973d-99_thumb[2]using    the de Casteljau algorithm and only keep the part overtmp973d-100_thumb[2]That    part    is    then    reparameterized    to    [0,1].

Note that we are being conservative here and we could have trimmed a little more. As Figure 13.6(b) shows, the curve D(u) meets the fat line Y over an interval smaller thantmp973d-101_thumb[2]

After clipping p(u) against q(u), we repeat this process, this time clipping q(u) against the fat line associated to the clipped p(u). We alternate back and forth between clipping against p(u) and q(u) until we arrive at our intersection point. This is the basic idea behind the Bezier clipping algorithm for finding curve intersections, but there are some important details that need to be added.

The first observation is that we can use any fat line for the clipping. The choice above is usually effective, but another choice that is sometimes better is to use a fat line orthogonal to the line containing p0 and pn. An example of that is shown in Figure 13.7. The fat line between L1 and L2 is a better choice when clipping q(u) against p(u) than using one that is parallel to the line containing the endpoints A and B of the curve p(u). The extra time used in checking which choice is better seems to be worth it.

Another point that has not been addressed is whether or not the above algorithm actually converges. In fact, it may not if we only clip little if anything at each stage.

An alternative fat line.

Figure 13.7. An alternative fat line.

Multiple-curve intersections.

Figure 13.8. Multiple-curve intersections.

Figure 13.8 shows an example of multiple intersections, which lead to such situations. Therefore, the suggested solution for the problem is the following recursive subdivision step: If Bezier clipping does not reduce the parameter range of each of the curves by at least 20%, then subdivide the curve with the largest remaining parameter interval and do a Bezier clip of the other curve with the two halves. Although this solves our convergence problem, if curves are almost tangent or two intersections are very close, then the algorithm reduces to a divide-and-conquer binary search type algorithm in the neighborhood of those points. There is a very efficient way to handle that situation, but we refer the reader to [SedN90] for the details. This finishes our discussion of the Sederberg-Nishita Bezier clipping algorithm. Timing comparisons in [SedN90] show that for curves of degree less than five, the implicitization algorithm in [SedP86] runs typically between 1 and 3 times faster than the Bezier clip algorithm, which in turn is between 2 and 10 times faster than the algorithms from [LanR80] and [KopM83]. For higher-degree curves, the Bezier clipping algorithm seems to win. [Rock90] also has a comparison of various divide-and-conquer methods for curves.

Curve Algebraic Methods

[SedP86] describes an algorithm for computing the intersection of two rational curves. By implicitizing one of the curves and substituting the parameterization of the other into the equation for the first one ends up having to solve a set of homogeneous linear equations. Another approach combining elimination theory and matrix computations is described in [ManD94] and [ManK97]. An algorithm for finding the intersection of two algebraic curves is described in [Sede89].

In general, depending on how smooth curves are presented, earlier comments imply that the intersection problem can be reduced to finding the roots of an equation of one or two variables. The interesting case is the two-variable case that involves solving a polynomial equation of the form

tmp973d-112_thumb[2]

like equation (13.6). This is the type of problem that algebraic geometry can help solve and is really part of the more general problem of computing implicitly defined objects about we have more to say in the next topic.

Sometimes one wants not only the points in the plane that constitute the intersection of two parameterized curves but the parameter values at the points of intersection. If one used an implicitization approach to find the intersection, then the parameter variables would have disappeared. This means that after finding an intersection point q, one is left with another problem, namely, that of finding the parameter u so that p(u) = q.

Special Surface Intersections

Ray-Surface Intersections

Problem 10.2.2.6 solved the ray-surface intersection problem for the faceted surface case. Now assume that S is a parameterized surface with parameterization p(u,v) and that X is a ray from a point p in a direction v. Let L be the line through p with direction vector v. Using a rigid motion we can transform this situation into one where p is the origin and v is e3, so that L is the z-axis. Like in the ray-curve case, there is little cost associated to assuming that this has been done.

If p(u,v) = (p1(u,v),p2(u,v),p3(u,v)), then finding the point where the z-axis pierces S is equivalent to solving the equations

tmp973d-113_thumb[2]

along with the condition thattmp973d-114_thumb[2]These    equations can be solved using a

Newton-Raphson method but one needs a starting guess at a solution. One way to get one is to subdivide the u-v dom in into small rectangles with cornerstmp973d-115_thumb[2]and  then approximate the surface by the grid of points p(ui,vj). If we find a quadrilateral pierced by the ray, then any one of its corners can be used as an initial guess for the Newton-Raphson iteration.

Again, like in all Newton-Raphson method approaches, one has to be ready to deal with at least two problems. The first is that the partial derivatives of the functions p1(u,v) and p2(u,v) may not be linearly independent at some points causing the usual problems with a Newton-Raphson method. The other more complicated problem is caused by the fact that we have a bounded domain and the iteration may want to cross one of the boundary lines. When this happens one should check if the solution does not in fact lie on the boundary by applying a Newton-Raphson iteration along the boundary, which is a one-dimensional problem. There will be additional problems at the corners of the patch. If at any time the solution seems to want to move back inside the patch, revert to the two-dimensional search problem.

Curve-Surface Intersections

We only consider the case where both the curve and surface are smooth subsets of R3 and are defined by parameterizations. Assuming that the curve is parameterized by q(t) and the surface by p(u,v), we need to solve the equation

tmp973d-118_thumb[2]

We describe a Newton-Raphson approach. The idea is to start with a guesstmp973d-119_thumb[2] for a solution to equation (13.11) and then define a sequencetmp973d-120_thumb[2]that converges

to an actual solution. Assume that we have already definedtmp973d-121_thumb[2]Using the linear function

tmp973d-125_thumb[2]

as an approximation to the functiontmp973d-126_thumb[2]in    a    neighborhood    oftmp973d-127_thumb[2]we    get

the next iteratetmp973d-128_thumb[2]by    solvingtmp973d-129_thumb[2]and

tmp973d-130_thumb[2]It is easy to show that so that we need to solve

tmp973d-137_thumb[2]tmp973d-138_thumb[2]tmp973d-139_thumb[2]

Buttmp973d-140_thumb[2]is orthogonal totmp973d-141_thumb[2]and

tmp973d-144_thumb[2]

Finally, lettmp973d-145_thumb[2]be    the normal to the surface attmp973d-146_thumb[2]

Then equation (13.13) can be solved for v. A similar argument can be applied to the variables u and t. What we finally get are equations

tmp973d-149_thumb[2]

The sequencetmp973d-150_thumb[2]defined by  equations (13.14) will converge to  a  point in the inter section of our curve and surface provided that we do not run into the usual problems associated to Newton-Raphson methods. Getting all the intersection points hinges on being able to come up with enough initial guesses.

Next post:

Previous post: