4

You can find the intersection points of two Bezier curves fully analytically using a technique called implicitization. It's nearly fully analytical, as you have to find the roots of a 9th-degree polynomial, but let's not be picky. The best resource I found on this subject are the Computer Aided Geometric Design notes by Thomas W. Sederberg.

I implemented this technique. I found out that it works really well (also performance-wise) in most cases, including cases one might call tough (e.g., nine points of intersection). The following screenshots come from the custom GUI helper tool using that implementation.

Two Bezier curves with nine intersections

Interestingly, some cases appear extremely simple, yet pose a great challenge to the implicitization-based intersection finding algorithm. These two curves have a single prominent intersection:

Two Bezier curves with a single intersection

The simple approximation intersection finding algorithms, like subdivision-based intersection finding, don't have any issues with this case (it isn't even challenging in any way to them). For implicitization, though, this is kryptonite.

If we minimally adjust a single control point, the implicitization-based algorithm magically starts finding the solution:

Two Bezier curves with a single intersection, minimally adjusted

I picked an example that's not trying to hide the underlying problem. These are two simple Bezier curves that were "extracted" from a single self-intersecting curve.

Another example of such "problematic" curves, this time using screenshots from a vector computer graphics application:

Pair of intersecting Bezier curves

Now, highlighting their underlying curve:

A loop Bezier curve

Their intersection polynomial is a constant 0 polynomial. It has infinitely many solutions, as the implicitization-based approach can't help but analyze the curves in their whole t-value domain (not just in the primary 0..1 range).

One might argue that this is actually a correct answer (as both curves actually overlap in a certain sense), but it's not. The curves we started with are defined in the 0..1 range, and it's a mathematical nuance how they "behave" outside of that range. In fact, this is a fundamental weakness of the implicitization-based intersection finding algorithm when applied to the standard Bezier curves defined in the 0..1 range.

How to analytically find the intersection of two cubic Bezier curves (restricted to the 0..1 range) in general? By saying "in general", I mean that the method should work in all cases that appear trivial at first glance, including the corner case I mentioned. If that helps, it might not work well (or at all) in obviously degenerate cases, like Bezier curves degenerating to a line or a point.

While a fully or mostly analytical solution is preferred (as I've come this far sticking to as analytical methods as possible), some ideas for approximated solutions described from a mathematical perspective are also welcome. For example, one of my ideas is to detect that the two curves overlap outside of their 0..1 range and then just... "shake" them a bit (pseudo-randomly move their control points by an epsilon in any direction). This should give a nearly perfect solution by "tricking" the intersection polynomial to become non-zero. Sounds silly, but might work.

Ideally, the solution wouldn't add computational overhead in the general case (99.999999.. % of cases) to handle the pathological cases. One of my ideas is to pick any point on the first curve and then pass it to the second curve's inverse polynomial to find the second curve's t-value. From that t-value, we can easily get the corresponding point by evaluating its basis function. If that's not the same point that we started with (the one from the first curve), then we know that this is not the pathological case and both curves are truly different (and it's safe to solve the intersection polynomial). Otherwise, we can proceed to some more sophisticated checks, as we already know that this either is a pathological case (and it's fine that it requires more computations) or we're just very unlucky (and it's also fine to add more computations in rare unlucky cases).

cubuspl42
  • 145
  • 1
    Can't you just compute the resultant or the Gröbner basis? – Sugai Ryomei Jul 12 '25 at 19:01
  • @SugaiRyomei "Can't just compute the resultant" What resultant exactly? The implicitization technique involves computing a resultant, but the issue is that some curves, although apparently distinct, actually fully overlap if you consider their whole domain of t-values. The "kryptonite" pair of curves I gave as an example has this property. There are infinitely many (t1, t2) pairs that correspond to the same (x, y) point. Are we on the same page in this matter? If we are, I might be misunderstanding your suggestion. – cubuspl42 Jul 13 '25 at 09:42
  • Oh sorry, I was misunderstanding your example. – Sugai Ryomei Jul 13 '25 at 11:05
  • Can you clarify what you mean with the curves overlapping? It is hard to understand from the pictures. – Jap88 Jul 13 '25 at 11:30
  • I added some extra explanation and more pictures, now it should be clearer – cubuspl42 Jul 13 '25 at 11:30
  • @Jap88 Do the new pictures I added a few seconds ago help? – cubuspl42 Jul 13 '25 at 11:30
  • By "overlapping" I mean two Bezier curves that share exactly the same set of points. They might have different control points, but one "completely covers" another. In the corner cases I mention, this takes place only in the real range, while the curves appear innocent when narrowed to 0..1 range. – cubuspl42 Jul 13 '25 at 11:33
  • Yes now i understand. The complete curve is called the Zariski closure and what you are looking for in this case is a self intersection. A numerically stable way to handle this is a bit complicated and best made by solving a a generalized eigenvalue or svd problem involving resultant matrices. I suggest you search for: "solving for Bezier curve intersections using resultant matrices and eigenvalues including rank deficiencies ". This is will give you what is needed. You cannot form the actual implicit form directly but needed to keep it on resultant matrix form in order to handle robustly. – Jap88 Jul 13 '25 at 11:40
  • I inserted that phrase to Google and the results aren't super satisfying... Is there any paper or any specific source (like a link) you have in mind? Does the solution you have in mind that I'm not informed up-front whether I'm dealing with the general case or a corner case? This is the reason I phrased this question the way I did - I intentionally didn't say "I have this corner case, how to solve it?" as it's not my problem. I'm given two cubic curves described by four points, and I have to figure out what to do. – cubuspl42 Jul 13 '25 at 11:48
  • While I might be missing something, verifying analytically that the two curves are problematic isn't trivial. What I said about the zero polynomial is my conceptual understanding of the problem. In reality, the intersection polynomial just can't be trusted and starts misbehaving numerically when you try to find the intersections of two curves that fully overlap (outside of their 0..1 range). I have some heuristics in mind, but I'd still like to know what's the general mathematical algorithm for finding the intersections of two arbitrary Bezier curves. – cubuspl42 Jul 13 '25 at 11:52
  • @cubuspl42 I provided references as an answer below. The details are very complicated and boils down to writing a lot of code with lots of case handling. – Jap88 Jul 14 '25 at 19:37

2 Answers2

2

These papers and related papers in the references of these papers cover this:

A numerically stable way to handle this is typically handled by solving a a generalized eigenvalue or svd problem involving resultant matrices.

The trick in order to handle this in general is that you cannot form the actual implicit form directly but needed to keep it on resultant matrix form in order to handle robustly.

The papers will show methods based on numerical linear algebra (eigenvalue and svd computation) to solve.

In general, degenerate cases such as self intersection (as in OPs example) are indicated by rank deficiency of the resultant matrices. By manipulating these matrices one can also find those self-intersections as eigenvalues of singular values. The methods are complicated and the papers merely provide sketches of methods.

Another example of rank deficiency is when a 3rd degree Bezier curve actually is a 2nd degree curve when implicitizing. This too shows up as rank deficiency in the resultant matrices. Similarly if you place all control points on a line. This corresponds to overlapping lines and is rank deficient.

Jap88
  • 1,553
  • Thank you for these references! I'd want to point out that the problem that I originally described is not related to numerical stability only. As I understand it, numerical stability issues arise once you start implementing a mathematical method using finite precision arithmetic. The corner case I mentioned causes the implicitization method fail even in theory (unless I made a serious logical mistake somewhere). – cubuspl42 Jul 15 '25 at 09:50
  • The more complex methods you mentioned might indeed give a solution to the problem I encountered (I'm not making any argument against that), but it's somewhat interesting that when the more complex methods are introduced in these papers, the implicitization method is quoted as prior work, but it's not said that the newer methods work around a conceptual/mathematical limitation of the implitiziation-based approach; merely the subject of numerical stability is taken into consideration. – cubuspl42 Jul 15 '25 at 09:55
  • My whole mathematical reasoning is based on an assumption that if you build an intersection polynomial of B(t) and B(t - m) via implicitization, in theory you should get a zero polynomial, and it's not happening in software exactly for numerical stability reasons. Yet, if you did get a zero polynomial (as the theoretical math indicates), it still wouldn't give you the answer you are looking for (unless you amended your method with another supporting method), as you won't find the (quasi-self-) intersection by finding its roots (as it has infinite useless roots). – cubuspl42 Jul 15 '25 at 10:00
  • But my reasoning might be completely biased by the fact that I'm an engineer, and I want to build a fully working software solution. In mathematics, a limitation might possibly not be considered an actual limitation if there's a well-known solution for the limitation (even if the problem and its solutions aren't explicitly acknowledged in the mathematical papers). – cubuspl42 Jul 15 '25 at 10:08
  • The matrix-based implicitation methods are heavy-duty in this context and if one wants to simplify one can instead use Bezier polygon refinement-based methods. The nice thing with the matrix-based methods is that they are probably optimal with respect to robustness. They can be used to handle tangential interscections etc. However, one needs to add lots of "if statements" to handle all cases. – Jap88 Jul 16 '25 at 22:58
  • Related to this is detecting overlapping curves A and B.. An alternative approach is to sample a random point on curve A and see if it lies on curve B (incidence). This is a point-curve inversion problem and can be solved with some standard Newton iteration method. – Jap88 Jul 16 '25 at 23:01
  • Note that you shouldn't compute the intersection between B(t) and a shifted B(t-m) since the difference is zero at the point of intersection. This is just a reparameterization of the same curve. When you implicitize (however you do it) you will get the curve's algebraic representation which is insensitive to reparameterizations. Instead you consider the curve when intersected with itself as two curves B(s) and B(t) and implicitize as if they were different curves. – Jap88 Jul 16 '25 at 23:05
1

I read the coordinates of the control points from the example image and calculated the coefficients of the curves. Both curves are quadratic. Below are the Sage Math notebook and the execution results.

p0x=689; p1x=438; p2x=248; p3x=119
p0y=817; p1y=564; p2y=338; p3y=139
p4x=699; p5x=563; p6x=356; p7x=78
p4y=93; p5y=273; p6y=473; p7y=693
a0=-p0x+3*p1x-3*p2x+p3x
a1=3*p0x-6*p1x+3*p2x
a2=-3*p0x+3*p1x
a3=p0x
b0=-p0y+3*p1y-3*p2y+p3y
b1=3*p0y-6*p1y+3*p2y
b2=-3*p0y+3*p1y
b3=p0y
c0=-p4x+3*p5x-3*p6x+p7x
c1=3*p4x-6*p5x+3*p6x
c2=-3*p4x+3*p5x
c3=p4x
d0=-p4y+3*p5y-3*p6y+p7y
d1=3*p4y-6*p5y+3*p6y
d2=-3*p4y+3*p5y
d3=p4y
a0,a1,a2,a3
b0,b1,b2,b3
c0,c1,c2,c3
d0,d1,d2,d3

(0, 183, -753, 689)
(0, 81, -759, 817)
(0, -213, -408, 699)
(0, 60, 540, 93)

When I applied the implicit equation for a cubic curve, all the coefficients became $0$, as you described.

var('t u')
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
F0=a3^3*b0^3-a2*a3^2*b0^2*b1+a1*a3^2*b0*b1^2-a0*a3^2*b1^3+a2^2*a3*b0^2*b2-2*a1*a3^2*b0^2*b2-a1*a2*a3*b0*b1*b2+3*a0*a3^2*b0*b1*b2+a0*a2*a3*b1^2*b2+a1^2*a3*b0*b2^2-2*a0*a2*a3*b0*b2^2-a0*a1*a3*b1*b2^2+a0^2*a3*b2^3-a2^3*b0^2*b3+3*a1*a2*a3*b0^2*b3-3*a0*a3^2*b0^2*b3+a1*a2^2*b0*b1*b3-2*a1^2*a3*b0*b1*b3-a0*a2*a3*b0*b1*b3-a0*a2^2*b1^2*b3+2*a0*a1*a3*b1^2*b3-a1^2*a2*b0*b2*b3+2*a0*a2^2*b0*b2*b3+a0*a1*a3*b0*b2*b3+a0*a1*a2*b1*b2*b3-3*a0^2*a3*b1*b2*b3-a0^2*a2*b2^2*b3+a1^3*b0*b3^2-3*a0*a1*a2*b0*b3^2+3*a0^2*a3*b0*b3^2-a0*a1^2*b1*b3^2+2*a0^2*a2*b1*b3^2+a0^2*a1*b2*b3^2-a0^3*b3^3+(-3*a3^2*b0^3+2*a2*a3*b0^2*b1-2*a1*a3*b0*b1^2+2*a0*a3*b1^3-a2^2*b0^2*b2+4*a1*a3*b0^2*b2+a1*a2*b0*b1*b2-6*a0*a3*b0*b1*b2-a0*a2*b1^2*b2-a1^2*b0*b2^2+2*a0*a2*b0*b2^2+a0*a1*b1*b2^2-a0^2*b2^3-3*a1*a2*b0^2*b3+6*a0*a3*b0^2*b3+2*a1^2*b0*b1*b3+a0*a2*b0*b1*b3-2*a0*a1*b1^2*b3-a0*a1*b0*b2*b3+3*a0^2*b1*b2*b3-3*a0^2*b0*b3^2)*x+(a2^3*b0^2-3*a1*a2*a3*b0^2+3*a0*a3^2*b0^2-a1*a2^2*b0*b1+2*a1^2*a3*b0*b1+a0*a2*a3*b0*b1+a0*a2^2*b1^2-2*a0*a1*a3*b1^2+a1^2*a2*b0*b2-2*a0*a2^2*b0*b2-a0*a1*a3*b0*b2-a0*a1*a2*b1*b2+3*a0^2*a3*b1*b2+a0^2*a2*b2^2-2*a1^3*b0*b3+6*a0*a1*a2*b0*b3-6*a0^2*a3*b0*b3+2*a0*a1^2*b1*b3-4*a0^2*a2*b1*b3-2*a0^2*a1*b2*b3+3*a0^3*b3^2)*y+(3*a3*b0^3-a2*b0^2*b1+a1*b0*b1^2-a0*b1^3-2*a1*b0^2*b2+3*a0*b0*b1*b2-3*a0*b0^2*b3)*x^2+(3*a1*a2*b0^2-6*a0*a3*b0^2-2*a1^2*b0*b1-a0*a2*b0*b1+2*a0*a1*b1^2+a0*a1*b0*b2-3*a0^2*b1*b2+6*a0^2*b0*b3)*x*y+(a1^3*b0-3*a0*a1*a2*b0+3*a0^2*a3*b0-a0*a1^2*b1+2*a0^2*a2*b1+a0^2*a1*b2-3*a0^3*b3)*y^2-b0^3*x^3+3*a0*b0^2*x^2*y-3*a0^2*b0*x*y^2+a0^3*y^3
F0.simplify_full()
x=a0*t^3+a1*t^2+a2*t+a3
y=b0*t^3+b1*t^2+b2*t+b3
F1=c3^3*d0^3-c2*c3^2*d0^2*d1+c1*c3^2*d0*d1^2-c0*c3^2*d1^3+c2^2*c3*d0^2*d2-2*c1*c3^2*d0^2*d2-c1*c2*c3*d0*d1*d2+3*c0*c3^2*d0*d1*d2+c0*c2*c3*d1^2*d2+c1^2*c3*d0*d2^2-2*c0*c2*c3*d0*d2^2-c0*c1*c3*d1*d2^2+c0^2*c3*d2^3-c2^3*d0^2*d3+3*c1*c2*c3*d0^2*d3-3*c0*c3^2*d0^2*d3+c1*c2^2*d0*d1*d3-2*c1^2*c3*d0*d1*d3-c0*c2*c3*d0*d1*d3-c0*c2^2*d1^2*d3+2*c0*c1*c3*d1^2*d3-c1^2*c2*d0*d2*d3+2*c0*c2^2*d0*d2*d3+c0*c1*c3*d0*d2*d3+c0*c1*c2*d1*d2*d3-3*c0^2*c3*d1*d2*d3-c0^2*c2*d2^2*d3+c1^3*d0*d3^2-3*c0*c1*c2*d0*d3^2+3*c0^2*c3*d0*d3^2-c0*c1^2*d1*d3^2+2*c0^2*c2*d1*d3^2+c0^2*c1*d2*d3^2-c0^3*d3^3+(-3*c3^2*d0^3+2*c2*c3*d0^2*d1-2*c1*c3*d0*d1^2+2*c0*c3*d1^3-c2^2*d0^2*d2+4*c1*c3*d0^2*d2+c1*c2*d0*d1*d2-6*c0*c3*d0*d1*d2-c0*c2*d1^2*d2-c1^2*d0*d2^2+2*c0*c2*d0*d2^2+c0*c1*d1*d2^2-c0^2*d2^3-3*c1*c2*d0^2*d3+6*c0*c3*d0^2*d3+2*c1^2*d0*d1*d3+c0*c2*d0*d1*d3-2*c0*c1*d1^2*d3-c0*c1*d0*d2*d3+3*c0^2*d1*d2*d3-3*c0^2*d0*d3^2)*x+(c2^3*d0^2-3*c1*c2*c3*d0^2+3*c0*c3^2*d0^2-c1*c2^2*d0*d1+2*c1^2*c3*d0*d1+c0*c2*c3*d0*d1+c0*c2^2*d1^2-2*c0*c1*c3*d1^2+c1^2*c2*d0*d2-2*c0*c2^2*d0*d2-c0*c1*c3*d0*d2-c0*c1*c2*d1*d2+3*c0^2*c3*d1*d2+c0^2*c2*d2^2-2*c1^3*d0*d3+6*c0*c1*c2*d0*d3-6*c0^2*c3*d0*d3+2*c0*c1^2*d1*d3-4*c0^2*c2*d1*d3-2*c0^2*c1*d2*d3+3*c0^3*d3^2)*y+(3*c3*d0^3-c2*d0^2*d1+c1*d0*d1^2-c0*d1^3-2*c1*d0^2*d2+3*c0*d0*d1*d2-3*c0*d0^2*d3)*x^2+(3*c1*c2*d0^2-6*c0*c3*d0^2-2*c1^2*d0*d1-c0*c2*d0*d1+2*c0*c1*d1^2+c0*c1*d0*d2-3*c0^2*d1*d2+6*c0^2*d0*d3)*x*y+(c1^3*d0-3*c0*c1*c2*d0+3*c0^2*c3*d0-c0*c1^2*d1+2*c0^2*c2*d1+c0^2*c1*d2-3*c0^3*d3)*y^2-d0^3*x^3+3*c0*d0^2*x^2*y-3*c0^2*d0*x*y^2+c0^3*y^3
F1.simplify_full()

0
0

The roots can be found by applying the implicit equation for a quadratic curve.

var('t u')
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
F0=a3^2*b1^2-a2*a3*b1*b2+a1*a3*b2^2+a2^2*b1*b3-2*a1*a3*b1*b3-a1*a2*b2*b3+a1^2*b3^2+(-2*a3*b1^2+a2*b1*b2-a1*b2^2+2*a1*b1*b3)*x+(-a2^2*b1+2*a1*a3*b1+a1*a2*b2-2*a1^2*b3)*y+b1^2*x^2-2*a1*b1*x*y+a1^2*y^2
F0.roots(u,ring=CC)
x=a0*t^3+a1*t^2+a2*t+a3
y=b0*t^3+b1*t^2+b2*t+b3
F1=c3^2*d1^2-c2*c3*d1*d2+c1*c3*d2^2+c2^2*d1*d3-2*c1*c3*d1*d3-c1*c2*d2*d3+c1^2*d3^2+(-2*c3*d1^2+c2*d1*d2-c1*d2^2+2*c1*d1*d3)*x+(-c2^2*d1+2*c1*c3*d1+c1*c2*d2-2*c1^2*d3)*y+d1^2*x^2-2*c1*d1*x*y+c1^2*y^2
F1.roots(t,ring=CC)

[(-2.48777459512159, 1),
 (0.629823165976091, 1),
 (-3.74172878016745 - 2.50017946124494*I, 1),
 (-3.74172878016745 + 2.50017946124494*I, 1)]
[(0.501246287552726, 1),
 (3.67921018009561, 1),
 (5.23619828833076 - 2.54859050511897*I, 1),
 (5.23619828833076 + 2.54859050511897*I, 1)]

I selected the real root in the range $[0,1]$ and calculated the coordinates, they matched, so I think the intersection was found correctly.

UR=F0.roots(u,ring=CC)
TR=F1.roots(t,ring=CC)
u=UR[1][0]
t=TR[0][0]
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
x,y
x=a0*t^3+a1*t^2+a2*t+a3
y=b0*t^3+b1*t^2+b2*t+b3
x,y

(357.539900336524, 456.905142851098)
(357.539900336524, 456.905142851098)

Please also refer to my post https://qiita.com/kumanomi1000/items/bf69f2c32d52cc13b93b. This post includes C++ source code for an implementation example. Since it is switching between cubic, quadratic, and linear using "if" statement, I'm not sure if it meets your expectations, i.e., whether it can be said to be numerically stable.

kumanomi1000
  • 11
  • 1
  • I think that you misunderstood the question. I didn't ask about the specific coordinates of the intersection point of the example curves, but I described a conceptual challenge of the implicitization-based intersection finding algorithm and asked for possible general solutions. – cubuspl42 Dec 27 '25 at 11:25