2D intersection algorithm for cubic bezier curves... who has one?
Pomax
does anyone know where I might find an algorithm for determining the
intersections (numerical, not analytical) for two 2D cubic bezier curves,
worked out as source code in a relatively common programming language?
I've spent a few hours reading through several papers that discuss how one
might determine the intersections for bezier curves, without them actually
giving any concrete code that does what the papers talk about. If anyone has
example code that does this trick (I'm sure plenty of people on this
newsgroup have such code lying around, or know where to get it) in a
relatively common and understandable programming language (but please, no
VB), I would be most appreciative if I could have a look at it.
Ideally I'm looking for something like the code that's used by Tomoyuki
Nishita's java applet for determining insections (he uses Bezier Clipping,
but after a few hours of googling for "bezier clipping" I'm not really
anywhere nearer to having found an actual implementation that tells me how
the code works, rather than what the motivation behind using the procedure
that results in whatever code that ended up being used!). Please, no "this is
the math behind it, now write it yourself" solutions, it's a problem that's
been tackled enough times over the course of the last three decades to have
some kind of standard source code solution by now (plus, I'm not good enough
at analytical math to convert a theoretical solution for 2nd order parametric
equalities into concrete source code, anyway).
Many thanks to whoever feels like helping.
- Mike "Pomax" Kamermans
nihongoresources.com
Jongware
Pomax, would a conversion to straight line segments work for you?
There are a couple of different takes on that; one is divide the curve
into /n/ segments, regardless of its length & complexity -- very fast
--, another is to sub-divide the curve until you approach a reasonable
calculation limit per segment -- for the screen that'd typically be 2
pixels.
I can provide C code for both.
[Jw]
Angus Johnson
> does anyone know where I might find an algorithm for determining the
> intersections (numerical, not analytical) for two 2D cubic bezier curves,
> worked out as source code in a relatively common programming language?
I have a pending update (with published beta code*) to my polygon
clipping library 'Clipper' that does pretty much what you ask. While the
main Clipper library (written in Delphi, C++ and C#) only does polygon
clipping (see http://sourceforge.net/projects/polyclipping/), the update
also performs line (and bezier and spline) clipping. This update
however is currently written in Delphi (object pascal) only.
* http://www.angusj.com/delphi/clipper.php#download
Intersections are determined by first flattening the bezier curves into
series of straight lines (using the de Casteljau algorithm). Points of
intersection are determined by clipping one set of curves against
another (ie subjects and clips). However, the clipping library doesn't
return the points of intersection per se, nor the clipped flattened
curves, what is returned is a series of reconstructed 'clipped' beziers.
If you're only interested in the points of intersection, the clipping
code could be modified to trigger an 'event' whenever intersections are
detected.
For other possible solutions, the following link may help too:
http://stackoverflow.com/questions/109364/bezier-clipping
Pomax
> would a conversion to straight line segments work for you?
>
> There are a couple of different takes on that; one is divide the curve
> into /n/ segments, regardless of its length & complexity -- very fast
> --, another is to sub-divide the curve until you approach a reasonable
> calculation limit per segment -- for the screen that'd typically be 2
> pixels.
That would probably work quite well, yes. and "fast" is always appealing =)
Pomax
> Intersections are determined by first flattening the bezier curves into
> series of straight lines (using the de Casteljau algorithm). Points of
> intersection are determined by clipping one set of curves against
> another (ie subjects and clips). However, the clipping library doesn't
> return the points of intersection per se, nor the clipped flattened
> curves, what is returned is a series of reconstructed 'clipped' beziers.
> If you're only interested in the points of intersection, the clipping
> code could be modified to trigger an 'event' whenever intersections are
> detected.
While a library is a nice resource, it's typically too much code to easily
port to other languages. I had seen the discussion on the stackoverflow post,
but the lack of a link to the finished product made it simply too theoretical
for me, hence the request for an algorithm that already worked, rather than
hints on how to implement such an algorithm by oneself.
It is rather essential to get the actual intersection coordinates too
(because they will be used for visual feedback), so I hope Jongware's example
code will be exactly what I'm looking for =)
Thanks for your reply though!
Jongware
> Hi Jw,
>
>
>> would a conversion to straight line segments work for you?
>>
>> There are a couple of different takes on that; one is divide the curve
>> --, another is to sub-divide the curve until you approach a reasonable
>> calculation limit per segment -- for the screen that'd typically be 2
>> pixels.
>
> That would probably work quite well, yes. and "fast" is always appealing =)
>
> - Mike "Pomax" Kamermans
> nihongoresources.com
Okay. Don't have the routines with me but I'll look'em up for you to-nite.
[Jw]
Pomax
> Okay. Don't have the routines with me but I'll look'em up for you to-nite.
>
> [Jw]
>
Looking forward to them! And thank you.
- Mike
[Jongware]
wrote:
> On 11/1/2010 1:01:29 PM, Jongware wrote:
>> Okay. Don't have the routines with me but I'll look'em up for you
>> to-nite.
What a misery. My Win machine crashed a while ago, & I'm still figgering
out what's there and what's gone. Fortunately, I could find a couple o'old
sources with this: Bezier Subdivision. It's two routines, you're supposed
to call only the main one yourself. I see here that for speed, I didn't
test the length, but used a global (ugh) variable to count the number of
nestings, and that's because unlimited nesting may cause run time problems.
void beziersub (int ax,int ay, int bx,int by, int cx,int cy, int dx,int
dy, unsigned short color)
{
static int recur = 0;
int abx,aby;
int abcx,abcy;
int abcdx,abcdy;
int bcdx,bcdy;
int cdx,cdy;
recur++;
if (recur < 4)
{
abx = (ax + bx)>>1;
aby = (ay + by)>>1;
abcx = ((ax+cx)>>2) + (bx>>1);
abcy = ((ay+cy)>>2) + (by>>1);
abcdx = ((ax+dx)>>3) + ((bx+cx)>>3)+((bx+cx)>>2);
abcdy = ((ay+dy)>>3) + ((by+cy)>>3)+((by+cy)>>2);
bcdx = ((bx+dx)>>2) + (cx>>1);
bcdy = ((by+dy)>>2) + (cy>>1);
cdx = (cx + dx)>>1;
cdy = (cy + dy)>>1;
beziersub (ax,ay, abx,aby, abcx,abcy, abcdx,abcdy, color);
beziersub (abcdx,abcdy, bcdx,bcdy, cdx,cdy, dx,dy, color);
} else
{
ax >>= 8;
ay >>= 8;
dx >>= 8;
dy >>= 8;
line (ax,ay, dx,dy, color);
}
recur--;
}
void bezier (int ax,int ay, int bx,int by, int cx,int cy, int dx,int dy,
unsigned short color)
{
int abx,aby;
int abcx,abcy;
int abcdx,abcdy;
int bcdx,bcdy;
int cdx,cdy;
if (ax == bx && ay == by && dx == cx && dy == cy)
{
line (ax,ay, dx,dy, color);
return;
}
ax <<= 8; ay <<= 8;
bx <<= 8; by <<= 8;
cx <<= 8; cy <<= 8;
dx <<= 8; dy <<= 8;
abx = (ax + bx)>>1;
aby = (ay + by)>>1;
abcx = ((ax+cx)>>2) + (bx>>1);
abcy = ((ay+cy)>>2) + (by>>1);
abcdx = ((ax+dx)>>3) + ((bx+cx)>>3)+((bx+cx)>>2);
abcdy = ((ay+dy)>>3) + ((by+cy)>>3)+((by+cy)>>2);
bcdx = ((bx+dx)>>2) + (cx>>1);
bcdy = ((by+dy)>>2) + (cy>>1);
cdx = (cx + dx)>>1;
cdy = (cy + dy)>>1;
beziersub (ax,ay, abx,aby, abcx,abcy, abcdx,abcdy, color);
beziersub (abcdx,abcdy, bcdx,bcdy, cdx,cdy, dx,dy, color);
}
Here is the other one, adapted from Maxim Shemanarev's Antigrain code.
This is not C, it's Javascript, but the general principle should be the
same. nSteps, by the way, can be initialized to any large enough value you
care to test for.
// Code adapted from Maxim Shemanarev's AntiGrain
// http://www.antigrain.com/research/bezier_interpolation/
function curve4(x1, y1, //Anchor1
x2, y2, //Control1
x3, y3, //Control2
x4, y4, //Anchor2
pointList ) // Destination
{
var dx1 = x2 - x1;
var dy1 = y2 - y1;
var dx2 = x3 - x2;
var dy2 = y3 - y2;
var dx3 = x4 - x3;
var dy3 = y4 - y3;
var subdiv_step = 1.0 / (nSteps + 1);
var subdiv_step2 = subdiv_step*subdiv_step;
var subdiv_step3 = subdiv_step*subdiv_step*subdiv_step;
var pre1 = 3.0 * subdiv_step;
var pre2 = 3.0 * subdiv_step2;
var pre4 = 6.0 * subdiv_step2;
var pre5 = 6.0 * subdiv_step3;
var tmp1x = x1 - x2 * 2.0 + x3;
var tmp1y = y1 - y2 * 2.0 + y3;
var tmp2x = (x2 - x3)*3.0 - x1 + x4;
var tmp2y = (y2 - y3)*3.0 - y1 + y4;
var fx = x1;
var fy = y1;
var dfx = (x2 - x1)*pre1 + tmp1x*pre2 + tmp2x*subdiv_step3;
var dfy = (y2 - y1)*pre1 + tmp1y*pre2 + tmp2y*subdiv_step3;
var ddfx = tmp1x*pre4 + tmp2x*pre5;
var ddfy = tmp1y*pre4 + tmp2y*pre5;
var dddfx = tmp2x*pre5;
var dddfy = tmp2y*pre5;
var step = nSteps;
pointList.push ([x1, y1]); // Start Here
while(step--)
{
fx += dfx;
fy += dfy;
dfx += ddfx;
dfy += ddfy;
ddfx += dddfx;
ddfy += dddfy;
pointList.push ([fx, fy]);
}
pointList.push ([x4, y4]); // Last step must go exactly to x4, y4
}
Enjoy.
[Jw]
Angus Johnson
> I didn't test the length, but used a global (ugh) variable to count
> the number of nestings, and that's because unlimited nesting may cause
> run time problems.
> Here is the other one, adapted from Maxim Shemanarev's Antigrain code.
> This is not C, it's Javascript, but the general principle should be
> the same. nSteps, by the way, can be initialized to any large enough
> value you care to test for.
The issue of appropriate levels of recursion/nesting was discussed in
some detail in this newsgroup back in 2005 (and Maxim was part of that
discussion):
http://groups.google.com/group/comp.graphics.algorithms/tree/browse_frm/thread/d85ca902fdbd746e
The thread was started by Maxim's excellent article on this subject:
http://antigrain.com/research/adaptive_bezier/index.html
Finally here is my Delphi code for cubic beziers ...
function GetCBezierPoints(const control_points: array of TFloatPoint):
TArrayOfFloatPoint;
var
i, j, arrayLen, resultCnt: integer;
ctrlPts: array [ 0..3] of TFloatPoint;
const
cbezier_tolerance = 0.5; //customize to desired pixel precision
half = 0.5;
procedure RecursiveCBezier(const p1, p2, p3, p4: TFloatPoint);
var
p12, p23, p34, p123, p234, p1234: TFloatPoint;
begin
if abs(p1.x + p3.x - 2*p2.x) + abs(p2.x + p4.x - 2*p3.x) +
abs(p1.y + p3.y - 2*p2.y) + abs(p2.y + p4.y - 2*p3.y) <
cbezier_tolerance then
begin
if resultCnt = length(result) then
setLength(result, length(result) +128);
result[resultCnt] := p4;
inc(resultCnt);
end else
begin
p12.X := (p1.X + p2.X) *half;
p12.Y := (p1.Y + p2.Y) *half;
p23.X := (p2.X + p3.X) *half;
p23.Y := (p2.Y + p3.Y) *half;
p34.X := (p3.X + p4.X) *half;
p34.Y := (p3.Y + p4.Y) *half;
p123.X := (p12.X + p23.X) *half;
p123.Y := (p12.Y + p23.Y) *half;
p234.X := (p23.X + p34.X) *half;
p234.Y := (p23.Y + p34.Y) *half;
p1234.X := (p123.X + p234.X) *half;
p1234.Y := (p123.Y + p234.Y) *half;
RecursiveCBezier(p1, p12, p123, p1234);
RecursiveCBezier(p1234, p234, p34, p4);
end;
end;
begin
//first check that the 'control_points' count is valid ...
arrayLen := length(control_points);
if (arrayLen < 4) or ((arrayLen -1) mod 3 <> 0) then exit;
//start with a reasonable length return array to minimize
//unnecessary cpu cycles with repeated resizing ...
setLength(result, 128); //will shorten or length as needed
result[0] := control_points[0];
resultCnt := 1;
for i := 0 to (arrayLen div 3)-1 do
begin
for j := 0 to 3 do
ctrlPts[j] := control_points[i*3 +j];
RecursiveCBezier(ctrlPts[0], ctrlPts[1], ctrlPts[2], ctrlPts[3]);
end;
SetLength(result,resultCnt); //resize array to final length
end;
[Jongware]
wrote:
> The issue of appropriate levels of recursion/nesting was discussed in
> some detail in this newsgroup back in 2005 (and Maxim was part of that
> discussion):
Yeah, my motto is if you're going to borrow code, borrow from the best :-)
Do you know how deep your recursive code typically recurses? I can imagine
some worst-case scenario's ... I had to clip off my C version because of
some pretty bad performance (and that was with integer multiplication &
division -- hence, Pomax, the shifts by 8 in my C variant).
All said & done: while this code is great to draw beziers from, how would
the actual performance be wrt Pomax' original question -- determining the
intersection point(s) for two of 'em? Your (and mine -- properly adjusted)
code will yield a list of straight lines, which have to be calculated and
stored *twice* (once for each bezier curve); then each straight segment of
one of the curves has to be tested with a fairly expensive test against
each of the other one's straight segments.
(*Fairly* expensive; a simple bounding box test will disregard most of the
segments right away. But still.)
It'd be worthwhile to take a look at the available literature for a
mathematical solution. It'd also have the advantage of being *exact*,
whereas converting to straight lines always is an approximation.
[Jw]
--
Using Opera's revolutionary email client: http://www.opera.com/mail/
Angus Johnson
> Do you know how deep your recursive code typically recurses? I can
> imagine some worst-case scenario's ... I had to clip off my C version
> because of some pretty bad performance (and that was with integer
> multiplication & division -- hence, Pomax, the shifts by 8 in my C
> variant).
No, I don't but I've never had any problems. Given that it's a binary
division of each curve segment the recursion can't really be all that deep.
Fernando Cacciola
> It'd be worthwhile to take a look at the available literature for a mathematical
> solution.
There is no closed-form solution, and the numerical methods that have been
published all boil down to recusively subdividing the curve. The difference
being in that you don't effectively construct a polyline as a separate first
step (as in this thread) but you search for the intersection *along* the recursion.
> It'd also have the advantage of being *exact*, whereas converting to
> straight lines always is an approximation.
>
Unfortunately not. Known methods still apply a succesive approximation.
P.S.: I haven't post any reference to the published algorithms because the ones
I know of are papers or books, but not code, which the OP needs.
But I would look around in David's site though:
http://www.geometrictools.com/index.html
HTH
--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com
Pomax
<< There is no closed-form solution, and the numerical methods that have been
published all boil down to recusively subdividing the curve. The difference
being in that you don't effectively construct a polyline as a separate first
step (as in this thread) but you search for the intersection *along* the
recursion.>>
I wrote an email to Richard Kinch, whose page on finding all intersections
for two bezier curves is basically the top hit on google, asking him where I
might find his efficient "unpublished algorithm" --at the time his page was
written, in 2007-- and he simply replied with "Flatten one curve, then the
intersection of the other is a cubic polynomical which has a closed form
solution", so basically it's half the work because you're only flattening one
of the curves, it's "closed form", but I don't think you can call it "exact",
since it's still numerical in nature due to the curve flattening.
Granted, the precision would be arbitrary (pick subpixel precision and for
computer graphics it's a done deal) but if you want mathematically sound
immediate computation... let's just say that wolfram alpha returns a VERY
elaborate formula for the generalised solution:
http://www.wolframalpha.com/input/?i=(A*(1-t)^3+%2B+B*t*(1-t)^2+%2B+C*(1-t)*t^2+%2B+D*t^3)-(E*(1-t)^3+%2B+F*t*(1-t)^2+%2B+G*(1-t)
*t^2+%2B+H*t^3)+%3D+0 -- specifically, the "solutions for variable t". It's
not the fastes of computations to execute!
- Mike
Fernando Cacciola
> On 11/2/2010 4:58:56 AM, Fernando Cacciola wrote:
>
> << There is no closed-form solution, and the numerical methods that have been
> published all boil down to recusively subdividing the curve. The difference
> being in that you don't effectively construct a polyline as a separate first
> step (as in this thread) but you search for the intersection *along* the
> recursion.>>
>
> I wrote an email to Richard Kinch, whose page on finding all intersections
> for two bezier curves is basically the top hit on google, asking him where I
> might find his efficient "unpublished algorithm" --at the time his page was
> written, in 2007-- and he simply replied with "Flatten one curve, then the
> intersection of the other is a cubic polynomical which has a closed form
> solution",
Ah, now that is progress. In fact, twice as good as the previous results.
> so basically it's half the work because you're only flattening one
> of the curves, it's "closed form",
Only the second part is closed form of course. The solution to the entire
problem is not.
> but I don't think you can call it "exact",
> since it's still numerical in nature due to the curve flattening.
>
Rigth.
> Granted, the precision would be arbitrary (pick subpixel precision and for
> computer graphics it's a done deal) but if you want mathematically sound
> immediate computation... let's just say that wolfram alpha returns a VERY
> elaborate formula for the generalised solution:
>
>
> http://www.wolframalpha.com/input/?i=(A*(1-t)^3+%2B+B*t*(1-t)^2+%2B+C*(1-t)*t^2+%2B+D*t^3)-(E*(1-t)^3+%2B+F*t*(1-t)^2+%2B+G*(1-t)
> *t^2+%2B+H*t^3)+%3D+0 -- specifically, the "solutions for variable t". It's
> not the fastes of computations to execute!
>
Definitely not.
But still an improvement.
Pomax
thank you very much for the code examples! I shall pour over them this
evening and try to make them work for my application.
It's actually not that much of a problem than one is in C and the other in
javascript: I'm working in the "Processing" language, which is a
visualisation programming language based on java syntax. It's wonderful for
its use of 2D drawing primitived (want a bezier curve? just call
bezier(x1,y1, x1c,y1c, x2c,y2c, x2,y2) and there's a bezier curve on your
screen, without any fiddling on lower level graphic context updates etc) but
it completely lacks any predefine polygon/polybezier shape operations like
unions, intersections and complements.
While it supports external libraries to give 2D shape algebra, one of the
super benefits of Processing is that there's also a "processing.js"
initiative, which is a project to run Processing code directly on a webpage,
by virtue of a javascript translation parser that turns Processing source
into javascript source, and then runs it inside a <canvas> element. Of
course, if you want to write Processing source that runs without problems on
a webpage too, that means you can't rely on libraries... but the upside is
it's a ridiculously simple programming language, so while I'll need to write
the polygon/polybezier shape algebra myself in order for things to keep
working on a webpage, it shouldn't be very hard with your code as example.
This should get me well on my way, thanks again!
Fernando Cacciola
> I'm working in the "Processing" language, which is a
> visualisation programming language based on java syntax.
I didn't know about it. I've been just looking at it (as well as processing,js,
spde and such) and I already like it a lot!
Cheers
Scott Hemphill
> On 2/11/2010 8:12 AM, [Jongware] wrote:
>> I didn't test the length, but used a global (ugh) variable to count
>> the number of nestings, and that's because unlimited nesting may
>> cause run time problems.
> ...
>> Here is the other one, adapted from Maxim Shemanarev's Antigrain
>> code. This is not C, it's Javascript, but the general principle
>> should be the same. nSteps, by the way, can be initialized to any
>> large enough value you care to test for.
>
>
> The issue of appropriate levels of recursion/nesting was discussed in
> some detail in this newsgroup back in 2005 (and Maxim was part of that
> discussion):
> http://groups.google.com/group/comp.graphics.algorithms/tree/browse_frm/thread/d85ca902fdbd746e
>
> The thread was started by Maxim's excellent article on this subject:
> http://antigrain.com/research/adaptive_bezier/index.html
[snip]
As far as I can tell, the highest performing intersection finders are
the ones which are the ones which use implicitization for one curve and
substitution for the other to create a 9th degree polynomial. They then
find the roots of the polynomial and check that the parameters are in
range. Here's a good paper on the topic:
http://tom.cs.byu.edu/~tom/papers/cubiciii.pdf
Scott
--
Scott Hemphill hemp...@alumni.caltech.edu
"This isn't flying. This is falling, with style." -- Buzz Lightyear
Pomax
> As far as I can tell, the highest performing intersection finders are
> the ones which are the ones which use implicitization for one curve and
> substitution for the other to create a 9th degree polynomial. They then
> find the roots of the polynomial and check that the parameters are in
> range.
Ah, another great moment in "The Sciences Raping English"... "implicitize"
has made it firmly onto my monthly atrocious words list =)
That aside, given that there is no generalised solution for root finding of
any polynomial of degree 5 or higher, I am slightly skeptical that this is
really the highest performing intersection finding algorithm for low level
polynomial intersections. I'm looking at the pdf you linked to, but it's a
bit more mathematics than I can read through in a casual sitting, and appears
to be from 25 years ago. While of course anyone with a pen and a highly
mathematically trained brain can come up with a 'best way' to determine the
intersection of two polynomial curves, I find it hard to believe that this
approach has not either been improved upon by the masses of computer
engineers that graduated and found jobs in the world of CAD/CAM software
producing companies, or even completely rewritten because there are far more
efficient ways to do it.
Would you happen to know of a more recent work that either confirms this is
the best approach, or refers to it while offering a more efficient approach?
- Mike
Scott Hemphill
Here are some course notes (by the same author) dated 2009.
(See section 7.4.1):
http://tom.cs.byu.edu/~557/text/ch7.pdf