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
We describe an approach discussed in [HosL93]. Let
(1) Pick a pointon the curve p(u).
(2) Find the intersection of the tangent linewith the curve q(u).
Since this tangent line is defined byand has equation
the equation we need to solve is
or
Figure 13.4. Curve intersections with the Newton-Raphson method.
Letbe a solution. Now find the intersection of the tangent lineat
with the curve p(u). Let this point be
(3) Continue this process, getting sequencesStop whenis 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). Letbe the control points of p(u). See Figure 13.5. Let
L be the line throughRepresent the line by an equation of the form
Figure 13.5. The fat line for a Bezier curve.
and define
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 forand
set
Quadratic Curve.
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]:
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:
Sederberg-Nishita Bezier clipping.
Figure 13.6. Sederberg-Nishita Bezier clipping.
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
whereare the Bernstein polynomialsandare the control points of p(u). Let equation (13.8) be the equation of the line throughand let d(u) be the signed distance from p(u) to L. It is easy to show that
Finally, define the parametric curve
whereThe fact that the x-coordinate of D(u) is u follows from the identities
Figure 13.6(b) shows the graph of D(u) for the example in Figure 13.6(a). Clearly, if
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. Letand be the largest and smallest u so that the convex hull of the D; misses the fat line Y betweenover the intervalsrespectively.
For the example in Figure 13.6(b) it turned out thatWe
now subdivide the Bezier curveusing the de Casteljau algorithm and only keep the part overThat 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 than
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.
Figure 13.7. An alternative fat line.
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
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
along with the condition thatThese 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 cornersand 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
We describe a Newton-Raphson approach. The idea is to start with a guess for a solution to equation (13.11) and then define a sequencethat converges
to an actual solution. Assume that we have already definedUsing the linear function
as an approximation to the functionin a neighborhood ofwe get
It is easy to show that so that we need to solve
Finally, letbe the normal to the surface at
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
The sequencedefined 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.