Reliable test for intersection of two Bezier curves

Question

Is there a test which reliably decides whether two Bezier curves intersect or not?

I don't need to know how many intersections there are or at what parameters they appear at. I just would like to know if there is some intersection or not. By "reliably" I mean that there should be no false positives or false negatives: for instance, checking if the bounding boxes overlap is not "sufficient" as they might overlap even if there is no curve intersection.

Alternatively, if it helps, I would be also interested in a weaker test: it would answer "yes" if the curves intersect (one or more times). If the curves don't intersect this test's answer would be "yes" or "no", meaning that false positives would be allowed.

An example of what I'm after: according to Rational polynomial parametric/rational polynomial parametric curve intersection, "If the two wedges do not overlap, the curves cannot intersect more than once". This would answer my question, were it "more than zero" instead.

Answer 1

I'd test the convex hulls of the control points for overlap.

• If there is no overlap, then the answer is definitely "no" because the curves stay inside the convex hull.
• If there is overlap in such a way that the endpoints of the first curve are not inside the convex hull for the second curve, and vice versa, and the direct connections of the endpoints intersect, then the answer is "yes"
• Otherwise (i.e., there is overlap, but one curve might sneak around the other) the asnwer is "maybe" and we can recurse: Subdivide both of the curves (or just one? But which one?) and check each combination. This looks like it might become exponential but in most real life situations only one of the four sub-pairs should ba a "maybe" case again.

It may be advisable to stop recursion at a certain depth (and answer "yes"). Otherwise, touching curves might cause "stack overflow". Also, if you limit yourself to a fixed number (three, say) of recursions, you can both be sure to bound the complexity while at the same time producing false positives at most (and presumably only in rare situations. The heuristic is that after a few subdivisions, the subcurves are so little bent (i.e., the convex hulls are very "thin" compared to their lengths) that we may expect either a "no" or a clear "yes".

Answer 2

As far as I know, there is no certified and complete method that can determine arbitrary intersections without transformations and implicitization. However, the bounding box checks, including repeated subdivision of the curves to some application-dependent resolution, are excellent filters for curves in general position.

As you already noticed, a purely numerical method (including subdivision when you only use fixed-precision floating point arithmetics) will not be able to reliably handle degeneracies like (nearly) tangential intersections. The argument is the same as in the univariate case, just that the precision bounds are even worse in higher dimension.

For the bivariate solving, the main methods to analyze and intersect (pairs of) algebraic curves in monomial basis are RUR (rational univariate representation) by Fabrice Rouillier et al. and CAD (cylindrical algebraic decomposition). The former involves some not-so-lightweight algebraic machinery like Gröbner bases, the latter can be done with GCD and resultant computations only. Note that I'm oh-so-slightly biased both as an author of a state-of-the-art variant of the latter method and as a coauthor of Fabrice. Obviously, algebraic curve analysis is a well-studied topic, and I refer to the small selection of a vast number of references that is mentioned in that paper.

I also met and worked with Chee Yap who's the author of the paper you mentioned, but didn't know this work. From a very, very quick glance, he needs a genericity assumption that ensures that sufficient refinement eventually reveals only ordinary intersections. In particular, the method does not seem to handle tangential intersections.

Any complete approach with both variants will require exact or bigfloat arithmetic, but that alone will not give you a complete algorithm unless you compute with a priori worst-case precision bounds. More elegant formulations use adaptive precision handling and combine symbolic and numeric computations. The crucial ingredients are, again, univariate polynomial root isolation as well as GCD, (sub-)resultants and/or Gröbner basis computations, depending on your favorite method.

I vaguely remember rumours about resultant computations in Bernstein basis, but IIRC, none of them has been investigated and developed to the point of a complete solver. From a theoretical point of view, there is some evidence and conjectures – but no proof! – that conversion to implicit form is not a bottleneck step (and that, indeed, the previously mentioned method might be asymptotically optimal). Again, the conjecture basically boils down to univariate root isolation complexity and precision demand bounds. In practice, the situation might be different; but I guess that the difference is negligible for the degenerate or "almost degenerate" situations that remain after some filter steps (essentially, intersection of convex hulls of control point polygons).

BTW, for implementations, have a look at the CGAL project; it features a library for 2D arrangements according to the exact geometric computation paradigm. (Again, I'm biased because... You know. Just one of many, many wheels here in all those developments.)

Answer 3

I think you're asking for something impossible, so you'll have to compromise somewhere.

I assume that you will implement the algorithm using floating point arithmetic, and that's one of the big sources of trouble. Two curves can be arbitrarily close without intersecting, so you have to define what you mean by "intersection" in floating point arithmetic. I suggest that a reasonable definition of an intersection is a pair of points, one lying on each curve, whose distance apart is less than some small number $\varepsilon$. And, of course, $\varepsilon$ will need to be an order of magitude (or two) bigger than the error $\delta$ involved in calculating points on the curve. Using double precision arithmetic, and assuming your curves are of unit size, $\delta$ will be around $10^{-14}$, so $\varepsilon$ will need to be around $10^{-12}$, say.

Depending on your application, it might make sense to use a much larger value for $\varepsilon$. For engineering or manufacturing, $10^{-7}$ is probably good enough, and for graphics $10^{-5}$ might suffice. Once you have chosen $\varepsilon$, you know when you can stop the various subdivision processes that compute intersections.

There is a pretty thorough analysis of intersection algorithms in Tom Sederberg's notes. But most of the techniques are based on bounding volume tests and subdivision, which you don't seem to like.

In this chapter of his notes, and in a paper in Computer-Aided Design Volume 18, Issue 1, January–February 1986, Pages 58–63, Sederberg discusses intersection calculation by implicitization. This is much faster than subdivision-based techniques for curves of low degree (less that 5 or so). The idea is to obtain an implicit equation $F(x,y)=0$ for one of the curves. Then if the other curve has parametric equation $C(t)$, the intersection points are the roots of the polynomial equation $\phi(t) = F(C(t))=0$. In floating point arithmetic, $\phi(t)$ will never be zero, of course, so you still have to deal with the $\varepsilon$ problem I outlined above. And, in fact, in this scenario, the $\varepsilon$ problem is harder because the value of $\phi$ has no clear geometric meaning, so it's hard to decide how small a value indicates an "intersection".

Answer 4

If you are talking about quadratic bezier curves like from a TrueType font definition etc, then (except for the horizontal line case) you can solve for $x$ as a function of $y$ (something you'd already have code for, to rasterize the interior, etc.) and there are 2 cases (for simplicity, consider the full unbounded curves, for unbounded $t$):

1. non-degenerate case: $x = Ay + B \pm\sqrt{Cy + D}$

2. degenerate: $x = Ay^2 + By + C$

Set two of the above kinds of functions to be equal, and square away the square-roots, and you get a quartic formula that can be solved in closed form. (It seems you always get a quartic if one of the curves is non-degenerate, so there's no way to rotate or skew $x$ and $y$ to make it into something lower-order, unfortunately).

Be sure to back-test the answer because extending to unbounded $t$ introduces false positives. To back-test, chop the curve into sections that are monotonic in both $x$ and $y$ (something you'd do anyway for initial bounding-box tests. This also gives you a single square-root signum to work on at a time.) Then you get bounding boxes for the two sections you're working on and you can easily check whether your answer is within both boxes.

There is also the case where one of the "bezier curves" is a horizontal line, but if you happen to be using this for banded union/intersection operations the horizontal line case already falls away. Finally, as others have noted, floating-point arithmetic will have problems when the curves are "near intersection" so you have to accept some inaccuracy there. Oh and the quartic formula solver code is complicated and can easily end up with hard-to-find bugs.