zneak zneak - 28 days ago 13
C++ Question

Checking if two cubic Bézier curves intersect

For a personal project, I'd need to find out if two cubic Bézier curves intersect. I don't need to know where: I just need to know if they do. However, I'd need to do it fast.

I've been scavenging the place and I found several resources. Mostly, there's this question here that had a promising answer.

So after I figured what is a Sylvester matrix, what is a determinant, what is a resultant and why it's useful, I thought I figured how the solution works. However, reality begs to differ, and it doesn't work so well.




Messing Around



I've used my graphing calculator to draw two Bézier splines (that we'll call B0 and B1) that intersect. Their coordinates are as follow (P0, P1, P2, P3):

(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)


The result is the following, B0 being the "horizontal" curve and B1 the other one:

Two cubic Bézier curves that intersect

Following directions from the aforementioned question's top-voted answer, I've subtracted B0 to B1. It left me with two equations (the X and the Y axes) that, according to my calculator, are:

x = 9t^3 - 9t^2 - 3t + 2
y = 9t^3 - 9t^2 - 6t + 4





The Sylvester Matrix



And from that I've built the following Sylvester matrix:

9 -9 -3 2 0 0
0 9 -9 -3 2 0
0 0 9 -9 -3 2
9 -9 -6 4 0 0
0 9 -9 -6 4 0
0 0 9 -9 -6 4


After that, I've made a C++ function to calculate determinants of matrices using Laplace expansion:

template<int size>
float determinant(float* matrix)
{
float total = 0;
float sign = 1;
float temporaryMatrix[(size - 1) * (size - 1)];
for (int i = 0; i < size; i++)
{
if (matrix[i] != 0)
{
for (int j = 1; j < size; j++)
{
float* targetOffset = temporaryMatrix + (j - 1) * (size - 1);
float* sourceOffset = matrix + j * size;
int firstCopySize = i * sizeof *matrix;
int secondCopySize = (size - i - 1) * sizeof *matrix;
memcpy(targetOffset, sourceOffset, firstCopySize);
memcpy(targetOffset + i, sourceOffset + i + 1, secondCopySize);
}
float subdeterminant = determinant<size - 1>(temporaryMatrix);
total += matrix[i] * subdeterminant * sign;
}
sign *= -1;
}
return total;
}

template<>
float determinant<1>(float* matrix)
{
return matrix[0];
}


It seems to work pretty well on relatively small matrices (2x2, 3x3 and 4x4), so I'd expect it to work on 6x6 matrices too. I didn't conduct extensive tests however, and there's a possibility that it's broken.




The Problem



If I understood correctly the answer from the other question, the determinant should be 0 since the curves intersect. However, feeding my program the Sylvester matrix I made above, it's -2916.

Is it a mistake on my end or on their end? What's the correct way to find out if two cubic Bézier curves intersect?

Answer

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

Although they don't explain the algorithm they actually use there, except to say that it's from p. 137 of "The Metafont Book", it appears that the key to it is two important properties of Bezier curves (which are explained elsewhere on that site though I can't find 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 and an algorithm for intersecting polygons, you 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 will be fast if the curves don't intersect -- is that the usual case?

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

Comments