[ACCEPTED]-Checking if two cubic Bézier curves intersect-bezier

Accepted answer
Score: 22

Intersection of Bezier curves is done by 14 the (very cool) Asymptote vector graphics language: look 13 for intersect() here.

Although they don't explain the algorithm 12 they actually use there, except to say that 11 it's from p. 137 of "The Metafont Book", it 10 appears that the key to it is two important 9 properties of Bezier curves (which are explained 8 elsewhere on that site though I can't find 7 the page right now):

  • A Bezier curve is always contained within the bounding box defined by its 4 control points
  • A Bezier curve can always be subdivided at an arbitrary t value into 2 sub-Bezier curves

With these two properties 6 and an algorithm for intersecting polygons, you 5 can recurse to arbitrary precision:

bezInt(B1, B2):

  1. Does bbox(B1) intersect bbox(B2)?
    • No: Return false.
    • Yes: Continue.
  2. Is area(bbox(B1)) + area(bbox(B2)) < threshold?
    • Yes: Return true.
    • No: Continue.
  3. Split B1 into B1a and B1b at t = 0.5
  4. Split B2 into B2a and B2b at t = 0.5
  5. Return bezInt(B1a, B2a) || bezInt(B1a, B2b) || bezInt(B1b, B2a) || bezInt(B1b, B2b).

This 4 will be fast if the curves don't intersect 3 -- is that the usual case?

[EDIT] It looks like 2 the algorithm for splitting a Bezier curve 1 in two is called de Casteljau's algorithm.

Score: 9

If you're doing this for production code, I'd 13 suggest the Bezier clipping algorithm. It's 12 explained well in section 7.7 of this free online CAGD text (pdf), works for any 11 degree of Bezier curve, and is fast and 10 robust.

While using standard rootfinders 9 or matrices might be more straightforward 8 from a mathematical point of view, Bezier 7 clipping is relatively easy to implement 6 and debug, and will actually have less floating 5 point error. This is because whenever it's 4 creating new numbers, it's doing weighted 3 averages (convex combinations) so there's 2 no chance of extrapolating based on noisy 1 inputs.

Score: 3

Is it a mistake on my end or on their end?

Are 18 you basing your interpretation of the determinant 17 on the 4th comment attached to this answer? If so, I 16 believe that's where the mistake lies. Reproducing 15 the comment here:

If the determinant is zero 14 there is is a root in X and Y at *exactly 13 the same value of t, so there is an intersection 12 of the two curves. (the t may not be in 11 the interval 0..1 though).

I don't see 10 any problems with this part, but the author 9 goes on to say:

If the determinant is <> zero 8 you can be sure that the curves don't intersect 7 anywhere.

I don't think that's correct. It's 6 perfectly possible for the two curves to 5 intersect in a location where the t values 4 differ, and in that case, there will be 3 an intersection even though the matrix has 2 a non-zero determinant. I believe this is 1 what's happening in your case.

Score: 2

I don't know how fast it will be, but if 8 you have two curves C1(t) and C2(k) they 7 intersect if C1(t) == C2(k). So you have 6 two equations (per x and per y) for two 5 variables (t, k). You can solve the system 4 using numerical methods with enough for 3 you accuracy. When you've found t,k parameters 2 you should check if there is a parameter 1 on [0, 1]. If it is they intersects on [0, 1].

Score: 2

This is a hard problem. I would split each 2 of the 2 Bezier curves into say 5-10 discrete 1 line segments, and just do line-line intersections.

enter image description here

foreach SampledLineSegment line1 in Bezier1
    foreach SampledLineSegment line2 in Bezier2
        if( line1 intersects line2 )
            then Bezier1 intersects Bezier2
Score: 1

I'm by no way an expert on this kind of 7 thing, but I follow a nice blog that talks a 6 lot about curves. He has link to two nice 5 articles talking about your problem (the 4 second link has an interactive demonstration 3 and some source code). Other people may 2 have much better insight into the problem 1 but I hope this helps!

http://cagd.cs.byu.edu/~557/text/ch7.pdf (archived copy)

Score: 0

I would say that the easiest and likely 99 fastest answer is to subdivide them into 98 very small lines and find the points where 97 the curves intersect, if they actually do.

public static void towardsCubic(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
    y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;
    xy[0] = x;
    xy[1] = y;
}

public static void towardsQuad(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * x0 + 2 * (1 - t) * t * x1 + t * t * x2;
    y = (1 - t) * (1 - t) * y0 + 2 * (1 - t) * t * y1 + t * t * y2;
    xy[0] = x;
    xy[1] = y;
}

Obviously 96 the brute force answer is bad See bo{4}'s 95 answer, and there's a lot of linear geometry 94 and collision detection that will actually 93 help quite a lot.


Pick the number of slices 92 you want for the curves. Something like 91 100 should be great.

We take all the segments 90 and sort them based on the largest value 89 of y they have. We then, add a second flag 88 in the list for the smallest value of y 87 for that segment.

We keep a list of active 86 edges.

We iterate through the y-sorted list 85 of segment, when we encounter a leading 84 segment we add it to the active list. When 83 we hit the small-y flagged value, we remove 82 that segment from the active list.

We then 81 can simply iterate through the entire set 80 of segments with what amounts to a scan 79 line, increasing our y monotonically as 78 we simply iterate the list. We iterate through 77 the values in our sorted list, which will 76 typically remove one segment and add a new 75 segment (or for split and merge nodes, add 74 two segments or remove two segments). Thereby 73 keeping an active list of relevant segments.

We 72 run the fast fail intersect check as we 71 add a new active segment to the list of 70 active segments, against only that segment 69 and and the currently active segments.

So 68 at all times we know exactly which line 67 segments are relevant, as we iterate through 66 the sampled segments of our curves. We know 65 these segments share overlap in the y-coords. We 64 can fast fail any new segment that does 63 not overlap with regard to its x-coords. In 62 the rare case that they overlap in the x-coords, we 61 then check if these segments intersect.

This 60 will likely reduce the number of line intersection 59 checks to basically the number of intersections.

foreach(segment in sortedSegmentList) {
    if (segment.isLeading()) {
        checkAgainstActives(segment);
        actives.add(segment);
    }
    else actives.remove(segment)
}

checkAgainstActive() would 58 simply check if this segment and any segment 57 in the active list overlap their x-coords, if 56 they do, then run the line intersect check 55 on them, and take the appropriate action.


Also 54 note, this will work for any number of curves, any 53 type of curves, any order of curves, in 52 any mixture. And if we iterate through the 51 entire list of segments it will find every 50 intersection. It will find every point a 49 Bezier intersects a circle or every intersection 48 that a dozen quadratic Bezier curves have 47 with each other (or themselves), and all 46 in the same split second.

The oft cited Chapter 7 document notes, for 45 the subdivision algorithm:

"Once a pair 44 of curves has been subdivided enough that 43 they can each be approximated by a line 42 segment to within a tolerance"

We can 41 literally just skip the middleman. We can 40 do this fast enough so to simply compare 39 line segments with a tolerable amount of 38 error from the curve. In the end, that's 37 what the typical answer does.


Secondly, note, that 36 the bulk of the speed increase for the collision 35 detection here, namely the ordered list 34 of segments sorted based on their highest 33 y to add to the active list, and lowest 32 y to remove from the active list, can likewise 31 be done for the hull polygons of the Bezier 30 curve directly. Our line segment is simply 29 a polygon of order 2, but we could do any 28 number of points just as trivially, and 27 checking the bounding box of all the control 26 points in whatever order of curve we're 25 dealing with. So rather than a list of active 24 segments, we would have a list of active 23 hull polygons points. We could simply use 22 De Casteljau's algorithm to split the curves, thereby 21 sampling as these as subdivided curves rather 20 than line segments. So rather than 2 points 19 we'd have 4 or 7 or whatnot, and run the 18 same routine being quite apt towards fast 17 failing.

Adding the relevant group of points 16 at max y, removing it at min y, and comparing 15 only the active list. We can thus quickly 14 and better implement the Bezier subdivision 13 algorithm. By simply finding bounding box 12 overlap, and then subdividing those curves 11 which overlapped, and removing those that 10 did not. As, again, we can do any number 9 of curves, even ones subdivided from curves 8 in the previous iteration. And like our 7 segment approximation quickly solve for 6 every intersection position between hundreds 5 of different curves (even of different orders) very 4 quickly. Just check all curves to see if 3 the bounding boxes overlap, if they do, we 2 subdivide those, until our curves are small 1 enough or we ran out of them.

Score: 0

Yes, I know, this thread is accepted and 35 closed for a long time, but...

First, I'd 34 like to thank you, zneak, for an inspiration. Your 33 effort allows to find the right way.

Second, I 32 was not quite happy with the accepted solution. This 31 kind is used in my favorite freeware IPE, and 30 its bugzilla is plenty of complains for 29 a low accuracy and reliability about an 28 intersection issue, my among them.

The missing 27 trick which allows to solve the problem 26 in a manner proposed by zneak : It is enough 25 to shorten one of curves by a factor k>0, then 24 the determinant of Sylvester matrix will 23 be equal to zero. It is obvious, if a shortened 22 curve will intersect, then original will 21 too. Now, the task is changed into the search 20 for a suitable value for k factor. This has 19 led to the problem of solving an univariate 18 polynomial of 9 degree. A real and positive 17 roots of this polynomial are responsible 16 for potential points of intersection. (This shouldn't 15 be a surprise, two cubic Bezier curves can 14 intersect up to 9 times.) The final selection 13 is performed to find only those k factors, which 12 provide 0<=t<=1 for both curves.

Now 11 the Maxima code, where the starting point 10 is set of 8 points provided by zneak :

p0x:1; p0y:1;
p1x:2; p1y:4;
p2x:3; p2y:4;
p3x:4; p3y:3;

q0x:3; q0y:5;
q1x:3; q1y:6;
q2x:0; q2y:1;
q3x:3; q3y:1;

c0x:p0x;
c0y:p0y;
c1x:3*(p1x-p0x);
c1y:3*(p1y-p0y);
c2x:3*(p2x+p0x)-6*p1x;
c2y:3*(p2y+p0y)-6*p1y;
c3x:3*(p1x-p2x)+p3x-p0x;
c3y:3*(p1y-p2y)+p3y-p0y;

d0x:q0x;
d0y:q0y;
d1x:3*(q1x-q0x);
d1y:3*(q1y-q0y);
d2x:3*(q2x+q0x)-6*q1x;
d2y:3*(q2y+q0y)-6*q1y;
d3x:3*(q1x-q2x)+q3x-q0x;
d3y:3*(q1y-q2y)+q3y-q0y;

x:c0x-d0x + (c1x-d1x*k)*t+ (c2x-d2x*k^2)*t^2+ (c3x-d3x*k^3)*t^3;
y:c0y-d0y + (c1y-d1y*k)*t+ (c2y-d2y*k^2)*t^2+ (c3y-d3y*k^3)*t^3;

z:resultant(x,y,t);

At this 9 point, Maxima answered:

(%o35)−2*(1004*k^9−5049*k^8+5940*k^7−1689*k^6+10584*k^5−8235*k^4−2307*k^3+1026*k^2+108*k+76)

Let Maxima solve 8 this equation:

rr: float( realroots(z,1e-20))  

The answer is:

(%o36) [k=−0.40256438624399,k=0.43261490325108,k=0.84718739982868,k=2.643321910825066,k=2.71772491293651]

Now the code 7 to select a right value of k:

for item in rr do ( 
  evk:ev(k,item),
  if evk>0 then (  
    /*print("k=",evk),*/ 
    xx:ev(x,item),  rx:float( realroots(xx,1e-20)),/*print("x(t)=",xx,"   roots: ",rx),*/
    yy:ev(y,item),  ry:float( realroots(yy,1e-20)),/*print("y(t)=",yy,"   roots: ",ry),*/
    for it1 in rx do (  t1:ev(t,it1),
    for it2 in ry do (  t2:ev(t,it2),
       dt:abs(t1-t2),
       if dt<1e-10 then (
         /*print("Common root=",t1,"  delta t=",dt),*/
         if (t1>0) and (t1<=1) then ( t2:t1*evk,
         if (t2>0) and (t2<=1) then (                           
                 x1:c0x + c1x*t1+ c2x*t1^2+ c3x*t1^3,
                 y1:c0y + c1y*t1+ c2y*t1^2+ c3y*t1^3,
                 print("Intersection point:     x=",x1, "      y=",y1)
)))))/*,disp ("-----")*/
));

Maxima's answer:

"Intersection point:     x="1.693201254437358"      y="2.62375005067273
(%o37) done

Theres 6 not only a honey, though. The presented 5 method is unusable, if:

  • P0=Q0, or more generally, if P0 lies on the second curve (or its extension). One can try to swap the curves.
  • an extremely rare cases, when both curves belong to one K-family (eg. their infinite extensions are the same).
  • keep an eyes on (sqr(c3x)+sqr(c3y))=0 or (sqr(d3x)+sqr(3y))=0 cases, here a quadratic are pretending to be a cubic Bezier curves.

One can ask, why 4 a shortening is performed only once. It's 3 enough because of reverse-inverse law, which 2 was discovered en passant, but this is an another 1 story.

More Related questions