# Computing Intersections Between a Cubic Bezier Curve and a Line

Posted on August 10th, 2013
Previous Article :: Next Article

## Introduction

Last week I was running a Starfish simulation of a molecular transport through a vent. Starfish was designed to support both linear and cubic spline representation of surfaces. But while testing the code by tracing particles, I noticed that the particle hits with curved surfaces were not being computed correctly. Digging deeper I found the issue: I never implemented the algorithm for finding intersection between a cubic and a line! Instead, the code was using the line-line intersection by connecting the two end points of the curve. Oops!

Before continuing, let me just say that I am really starting to enjoy algorithm development using interactive HTML technologies! In the past, I would write the code in Java or C++ and have it output the splines and intersection points to a file which I would subsequently visualize using some plotting program. Or using Matlab, I could do the plotting right from the program. But I would still be confined to testing just a single case. With HTML, we can do something much better: we can interactively manipulate the curves and see the algorithm respond in real time!

## Demo

You can try this out in the demo above. If everything loaded fine, you should see a blue cubic Bezier curve and a red line. You will also see two white circles, these are the two control points $$\mathbf{P}_1$$ and $$\mathbf{P}_2$$ defining the cubic. As you change the curves by dragging the large circles, you should see a small black dot track the intersection point. You will see additional points appear if you orient the curve such that there are multiple intersections. This is shown below in Figure 1.

## Math Background

So how does this code work? The visualization is similar to the article on smooth splines through prescribed points. This mathematical algorithm is based on this answer. One way to represent an infinitely-long line is as follows
$$\displaystyle \frac{(x-x_1)}{(x_2-x_1)} = \frac{(y-y_1)}{(y_2-y_1)}$$
which can be rewritten as
$$x(y_2-y_1) + y(x_1-x_2) + x_1(y_1-y_2)+y_1(x_2-x_1)=0$$
or
$$Ax + By + C = 0 \qquad (1)$$

We next constrain the $$(x,y)$$ pairs to those located on the cubic curve. A cubic Bezier curve is given by
$$\mathbf{r}(t)=(1-t)^3\mathbf{P}_0 + 3(1-t)^2t\mathbf{P}_1 + 3(1-t)t^2\mathbf{P}_2 + t^3\mathbf{P}_3, \quad t\in [0,1]$$
where $$\mathbf{r}(t)$$ is the position vector. If we substitute these $$(x,y)$$ components into equation (1), we obtain a cubic equation in $$t$$. Finding the intersection points is then a “simple” matter of finding the roots of the cubic equation.

### Cubic Roots

One way to find a single root is using Newton’s method. Unfortunately, a cubic can have up to 3 roots. This is because, as shown in Figure 1, a line can intersect a cubic spline in up to 3 locations. Since we are using this algorithm for particle tracing, we are interested in the first intersection along the line. There is no guarantee that the Newton’s method will converge to this root. As such, we need to find all existing roots and sort them. Finding additional roots with Newton’s method is possible but not trivial. Third-order polynomials also have an analytical solution for their roots. But unlike the well known quadratic formula, there are multiple equations for cubic roots. In the end, I ended up using algorithm from Stephen Schmitt’s site:

/*based on http://mysite.verizon.net/res148h4j/javascript/script_exact_cubic.html#the%20source%20code*/
function cubicRoots(P)
{
var a=P[0];
var b=P[1];
var c=P[2];
var d=P[3];

var A=b/a;
var B=c/a;
var C=d/a;

var Q, R, D, S, T, Im;

var Q = (3*B - Math.pow(A, 2))/9;
var R = (9*A*B - 27*C - 2*Math.pow(A, 3))/54;
var D = Math.pow(Q, 3) + Math.pow(R, 2);    // polynomial discriminant

var t=Array();

if (D >= 0)                                 // complex or duplicate roots
{
var S = sgn(R + Math.sqrt(D))*Math.pow(Math.abs(R + Math.sqrt(D)),(1/3));
var T = sgn(R - Math.sqrt(D))*Math.pow(Math.abs(R - Math.sqrt(D)),(1/3));

t[0] = -A/3 + (S + T);                    // real root
t[1] = -A/3 - (S + T)/2;                  // real part of complex root
t[2] = -A/3 - (S + T)/2;                  // real part of complex root
Im = Math.abs(Math.sqrt(3)*(S - T)/2);    // complex part of root pair

if (Im!=0)
{
t[1]=-1;
t[2]=-1;
}

}
else                                          // distinct real roots
{
var th = Math.acos(R/Math.sqrt(-Math.pow(Q, 3)));

t[0] = 2*Math.sqrt(-Q)*Math.cos(th/3) - A/3;
t[1] = 2*Math.sqrt(-Q)*Math.cos((th + 2*Math.PI)/3) - A/3;
t[2] = 2*Math.sqrt(-Q)*Math.cos((th + 4*Math.PI)/3) - A/3;
Im = 0.0;
}

/*discard out of spec roots*/
for (var i=0;i<3;i++)
if (t[i]<0 || t[i]>1.0) t[i]=-1;

/*sort but place -1 at the end*/
t=sortSpecial(t);

console.log(t[0]+" "+t[1]+" "+t[2]);
return t;
}

This algorithm returns an array of parametric intersection locations along the cubic, with -1 indicating an out-of-bounds intersection (before or after the end point or in the imaginary plane). We also need to verify that the intersections are within the limits of the linear segment. This is done by the following code:

/*computes intersection between a cubic spline and a line segment*/
function computeIntersections(px,py,lx,ly)
{
var X=Array();

var A=ly[1]-ly[0];	    //A=y2-y1
var B=lx[0]-lx[1];	    //B=x1-x2
var C=lx[0]*(ly[0]-ly[1]) +
ly[0]*(lx[1]-lx[0]);	//C=x1*(y1-y2)+y1*(x2-x1)

var bx = bezierCoeffs(px[0],px[1],px[2],px[3]);
var by = bezierCoeffs(py[0],py[1],py[2],py[3]);

var P = Array();
P[0] = A*bx[0]+B*by[0];		/*t^3*/
P[1] = A*bx[1]+B*by[1];		/*t^2*/
P[2] = A*bx[2]+B*by[2];		/*t*/
P[3] = A*bx[3]+B*by[3] + C;	/*1*/

var r=cubicRoots(P);

/*verify the roots are in bounds of the linear segment*/
for (var i=0;i<3;i++)
{
t=r[i];

X[0]=bx[0]*t*t*t+bx[1]*t*t+bx[2]*t+bx[3];
X[1]=by[0]*t*t*t+by[1]*t*t+by[2]*t+by[3];

/*above is intersection point assuming infinitely long line segment,
make sure we are also in bounds of the line*/
var s;
if ((lx[1]-lx[0])!=0)           /*if not vertical line*/
s=(X[0]-lx[0])/(lx[1]-lx[0]);
else
s=(X[1]-ly[0])/(ly[1]-ly[0]);

/*in bounds?*/
if (t<0 || t>1.0 || s<0 || s>1.0)
{
X[0]=-100;  /*move off screen*/
X[1]=-100;
}

/*move intersection point*/
I[i].setAttributeNS(null,"cx",X[0]);
I[i].setAttributeNS(null,"cy",X[1]);
}
}

As you can see, we are always plotting 3 intersection locations, but the out-of-bounds intersections are moved off screen to location (-100,-100). The above code also does not sort the intersections along the line, but this change is ease to implement by storing the s parametric positions in array.

### Source Code

And that’s it. You can download the code by right clicking and selecting “save as” on this link: cubic-line.svg

Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.

### 24 comments to “Computing Intersections Between a Cubic Bezier Curve and a Line”

1. Oli
January 13, 2014 at 8:40 pm

Hello,

first I want to say that I really like the blog and topics discussed here. I remember finding the Boris-push post years ago and it helped my a lot at that time.

But on the topic: How do you check if particles are inside your domain? In my problems I just have an arbitrary shaped boundary. I classify each cell as either inside, outside or boundary cell. For the boundary cells I use a (kind of reduced) point-in-polygon check. But how would one do it for cubic splines? Just always calculate if there is an intersection point?

I hope you can clear things up like in the past :).

Cheers

Oli

• August 28, 2014 at 4:43 pm

Hey Oli,

as a quick rejection of collision with bezier, you can check for collision with the hull created by the control points of the bezier. If it doesn’t collide with the hull, it can’t collide with the curve.

from wikipedia:
Bézier curves are widely used in computer graphics to model smooth curves. As the curve is completely contained in the convex hull of its control points, the points can be graphically displayed and used to manipulate the curve intuitively

2. August 12, 2014 at 10:47 am

Hi there,

I guess you may feel interested about my research on curve-to-curve intersection too, https://sites.google.com/site/curvesintersection/

Enjoy my show! Cheers,

Hunt Chang

• August 14, 2014 at 6:08 am

Thanks for sharing!

3. November 18, 2015 at 1:23 am

Hey! Just wanted to say this was incredibly helpful to me in a project I’m working on. Thanks a lot!

• November 18, 2015 at 4:01 pm

Thanks and glad to hear that!

4. esseci
March 11, 2016 at 4:48 am

Your code saved really a lot of time!
Many thanks!!

5. Luca
June 14, 2016 at 11:55 am

Thanks for sharing this, it really helped!

6. Ben
January 18, 2017 at 1:38 am

It’s good to use generic line equation.

7. Joshua Pilkington
May 23, 2017 at 10:47 pm

Thanks for the helpful code!

I did find one issue where the calculated coefficients could end up having zero values for either or both of the cube or square with certain configurations (e.g. for A*x^3 + B*x^2 + C*x + D, A and/or B could end up being zero). This causes a divide-by-zero error when attempting to calculate the cubic roots (resulting in the line segment going through the bezier curve but NOT showing an intersection point). I simply added quadratic and linear cases to handle such scenarios.

Here is a scenario that causes the zero coefficients:

	/*create control points for the cubic*/
P[0] = createMarkerV(100,50,"blue");
P[1] = createMarkerK(200,250);
P[2] = createMarkerK(300,250);
P[3] = createMarkerV(400,450,"blue");

/*create control points for the line*/
P[4] = createMarkerV(200,50,"red");
P[5] = createMarkerV(200,350,"red");


You can see it happen as you move one of the line segment endpoints around from being vertical to non-vertical.

To fix the issue, I added this to the cubicRoots function:

   if( P[0] == 0 )
{
if( P[1] == 0 )
{
console.log("Linear formula detected");

var t = Array();
t[0] = -1 * ( P[3] / P[2] );
t[1] = -1;
t[2] = -1;

/*discard out of spec roots*/
for (var i=0;i<1;i++)
if (t[i]1.0) t[i]=-1;

/*sort but place -1 at the end*/
t=sortSpecial(t);

console.log(t[0]);
return t;
}

var DQ = Math.pow(P[2], 2) + 4*P[1]*P[3]; // quadratic discriminant
if( DQ >= 0 )
{
DQ = Math.sqrt(DQ);
var t = Array();
t[0] = -1 * ( ( DQ + P[2] ) / ( 2 * P[1] ) );
t[1] = ( ( DQ - P[2] ) / ( 2 * P[1] ) );
t[2] = -1;

/*discard out of spec roots*/
for (var i=0;i<2;i++)
if (t[i]1.0) t[i]=-1;

/*sort but place -1 at the end*/
t=sortSpecial(t);

console.log(t[0]+" "+t[1]);
return t;
}
}


Again, awesome code!

• NounoursPanda
May 2, 2018 at 10:52 am

Hey man, thank a lot for this improvement ! But I found a syntax error. Here is the final code I think :
function cubicRoots(P)
{
if( P[0] == 0 )
{
if( P[1] == 0 )
{
console.log(“Linear formula detected”);

var t = Array();
t[0] = -1 * ( P[3] / P[2] );
t[1] = -1;
t[2] = -1;

/*discard out of spec roots*/
for (var i=0;i<1;i++)
if (t[i]1.0) t[i]=-1;

/*sort but place -1 at the end*/
t=sortSpecial(t);

console.log(t[0]);
return t;
}

var DQ = Math.pow(P[2], 2) – 4*P[1]*P[3]; // quadratic discriminant
if( DQ >= 0 )
{
DQ = Math.sqrt(DQ);
var t = Array();
t[0] = -1 * ( ( DQ + P[2] ) / ( 2 * P[1] ) );
t[1] = ( ( DQ – P[2] ) / ( 2 * P[1] ) );
t[2] = -1;

/*discard out of spec roots*/
for (var i=0;i<2;i++)
if (t[i]1.0) t[i]=-1;

/*sort but place -1 at the end*/
t=sortSpecial(t);

console.log(t[0]+” “+t[1]);
return t;
}
}

• NounoursPanda
May 2, 2018 at 10:56 am

well…
I understand now why there was an syntax error. This appear when we copy past in this comment section the code near “/*discard out of spec roots*/”

You just need to copy the original post for this part, don’t copy past the code I sent.

• Mark
December 20, 2019 at 1:46 am

Should be: if (t[i]<0 || t[i]>1.0) t[i]=-1;

When typing the symbols for less-than (<) and greater-than (>), the browser will first try to interpret them as HTML tags. Only when it fails to do so will it display them as symbols.

• May 19, 2018 at 5:54 pm

Thanks!

8. Joshua Pilkington
May 23, 2017 at 10:57 pm

By the way,

var DQ = Math.pow(P[2], 2) + 4*P[1]*P[3]; // quadratic discriminant


should have been

var DQ = Math.pow(P[2], 2) - 4*P[1]*P[3]; // quadratic discriminant


I accidentally used a ‘+’ in the place of a ‘-‘.

9. shengpei zheng
July 4, 2017 at 4:48 pm

Thank u a lot.it’s interesting.

10. Matt H.
January 18, 2018 at 6:20 pm

This is a GREAT demo. Any updates to this? It is very useful sample. I’m working with it trying to make it work for multiple lines but still allowing for multiple intersections per line, and I have also run into the issue Joshua ran into. If you have any updates or info to share that would be very appreciated.

• January 18, 2018 at 6:34 pm

I wasn’t really planning on making any updates to this, but I guess I should implement Joshua’s fixes. One of these days. I can help you with your code offline (email me directly) if needed. Truly, I am quite surprised by how much interest this article generated. This is really an off-shoot of my main work (plasma modeling) but it seems this is the most interesting article on the entire blog!

• Jesse
July 19, 2018 at 5:02 pm

That’s definitely because bezier curves are more widespread than plasma modeling

11. Sergiy
March 13, 2018 at 12:37 am

Hi,
Thanx for your code, it works perfectly.

However, it seems I cannot find the explanation for that part of code:

    var P = Array();
P[0] = A*bx[0]+B*by[0];		/*t^3*/
P[1] = A*bx[1]+B*by[1];		/*t^2*/
P[2] = A*bx[2]+B*by[2];		/*t*/
P[3] = A*bx[3]+B*by[3] + C;	        /*1*/


What does it do, what math is behind it?
It looks like this is the key to the whole thing of plugging two equations together and it’s not explained or I could not comprehend it.

Regards, Sergiy

• March 14, 2018 at 10:53 am

Hi Sergiy,

I probably should have included more of the intermediate math. Basically, we have the equation Ax+By+C=0. We can however rewrite both x and y as functions of t. This will then give you A*x(t) + B*y(t) + C = 0, which is now an equation in only a single unknown (t) that we can find root(s) for. A Bezier curve can be converted to form x = bx[0]*t^3 + bx[1]*t^2 + bx[2]*t + bx[3] (ignore the coefficient ordering, in hindsight they should be numbered in reverse). Similar equation exists for y. You can plug these into the A*x(t) + B*y(t) + C = 0 equation. You will then obtain the following equation (A*bx[0]+B*by[0])*t^3 + (A*bx[1]+B*by[1])*t^2 + (A*bx[2]+B*by[2])*t + (A*bx[3]+B*by[3]+C) = 0, or P0*t^3 + P1*t^2 + P2*t + P3 = 0.

12. NounoursPanda
May 2, 2018 at 10:50 am

Hello, thank’s a lot, this was very helpful for a project !!!!

13. March 3, 2019 at 7:35 pm

Thanks for this – very useful! I was able to translate this into NodeBox 3 (a visual language).

One problem. I found erratic behavior in some cases which turned out to be rounding errors when some of the intermediate values became extremely large or small. This happened in some cases when control point values had repeating decimals (e.g. 133.3333333333)

I was able to solve this by rounding the X and Y values of the initial four bezier points to two digits of precision (e.g. 133.33333333333 becomes 133.33). This made no visible impact on the lines, curves, or intersection points and solved the problem.

14. Ben Zamos
June 20, 2019 at 1:39 pm

Hi! How would I use this method in a C++ program?