6 ms·

I love bezier curves. I'm working on rebuilding a project I started in college, a vector graphics editor. The function for determining intersections between two

I love bezier curves. I'm working on rebuilding a project I started in college, a vector graphics editor. The function for determining intersections between two cubic bezier curves, which I found on a professor's website, is absolutely insane: https://github.com/artursapek/mondrian/blob/master/src/coffee/geometry/cubic-bezier-line-segment.coffee#L397 https://github.com/artursapek/mondrian/blob/master/src/coffe... It's a wall of this: c11.x*c12.y*c13.x*c13.y*c23.x*c23.y - c11.y*c12.x*c13.x*c13.y*c23.x*c23.y + 6*c21.x*c22.x*c13y3*c23.x + 3*c11.x*c12.x*c13.x*c13.y*c23y2 + 6*c10.x*c13.x*c13y2*c23.x*c23.y - 3*c11.x*c12.x*c13y2*c23.x*c23.y - 3*c11.y*c12.y*c13.x*c13.y*c23x2 - 6*c10.y*c13x2*c13.y*c23.x*c23.y - 6*c20.x*c13.x*c13y2*c23.x*c23.y + 3*c11.y*c12.y*c13x2*c23.x*c23.y - 2*c12.x*c12y2*c13.x*c23.x*c23.y - 6*c21.x*c13.x*c22.x*c13y2*c23.y - 6*c21.x*c13.x*c13y2*c22.y*c23.x - 6*c13.x*c21.y*c22.x*c13y2*c23.x + 6*c21.x*c13x2*c13.y*c22.y*c23.y + Maybe there's a way to refactor it, but I wouldn't know because I still don't understand it; I just know that it works.

nkoren 10y ago
What the heck...!?! I would really love for someone to explain how something like this happens.
- artursapek 10y ago
  Me too. It's constructing some polynomial, but the rest is lost on me.
  - T-hawk 10y ago
    Looks like what you would get if you wrote out all the operations for matrix multiplication? It must be related to expressing the coordinates of the nodes as matrices. I would guess the factors of 3 and 6 arise from derivatives of the cubics: x^3 derives to 3x^2 derives to 6x.
    - wnissen 10y ago
      The intersection of two NURBS (Non-uniform Rational B-Spline) surfaces has a degree into the hundreds. Those are the surface types used in virtually all CAD modelers. For that reason the intersection is not calculated exactly but approximated to a geometric tolerance. It is harder, mathematically and algorithm-wise to calculate the approximation, but much faster to evaluate afterwards. It's amazing what magic lurks below the surface of even something as simple as "find the intersection of two polynomial functions".
      - surement 10y ago
        > NURBS (Non-uniform Rational B-Spline) I believe NURBS stands for Nobody Understands Rational B-Splines [0] :) [0] http://www.springer.com/productFlyer_978-3-540-61545-3.pdf?SGWID=0-0-1297-1491853-0 http://www.springer.com/productFlyer_978-3-540-61545-3.pdf?S... (pdf)
- phkahler 10y ago
  It happens because the curves are not polynomials in the normal sense. The parameter t that ranges from [0..1] is not the same between the two curves, so it is not equivalent to setting two polynomials equal to each other and solving for x, or t in this case. A solution to the intersection is equivalent to solving for the t parameter of both curves simultaneously and IIRC that ends up having an order that is the product of the orders of the curves (9th order for 2 cubics?). On second thought it's not the product, but it is definitely much higher order. That is the reason people use the matrix representations, once you get a feel for that the complexity tends to get hidden by representing an entire matrix of Nth order polynomials by a letter like B. If anyone is really interested in this stuff and wants to make a major contribution to open source software... First study these and then move on to NURBS, surfaces, and then trimmed NURBS surfaces. These are the core of all CAD software and a clean library for trimmed NURBS is sorely lacking. There are two that I'm aware of, the one in OpenCascade which is a huge ball of mud (but works) and the one in SolveSpace which doesn't scale well and has some other limitations - that is the one I feel is most worthy of someones time.
  - Klapaucius 10y ago
    There is GoTools (https://github.com/SINTEF-Geometry/GoTools https://github.com/SINTEF-Geometry/GoTools) and SISL (https://github.com/SINTEF-Geometry/SISL https://github.com/SINTEF-Geometry/SISL). These are open-source geometry libraries for parametric surface modeling, that include NURBS surfaces and intersections and trimmings of such. Intersections are computed approximatively within some prescribed geometrical tolerance, and intersection curves are "marched out".
    - phkahler 10y ago
      Thanks. The links have been forwarded to people who may be interested. BTW they don't explicitly say "trimmed NURBS" but doing boolean operations implies that is the case ;-)
  - blt 10y ago
    NLib is the gold standard in closed-source, right?
TheRealPomax 10y ago
This kind of feels like an export of a WolframAlpha query, exported and written out as code =D
- santaclaus 10y ago
  Haha yea that is some straight up Mathematica -> C codegen (damn you //cform). A couple of years ago I inherited a project with a 90 meg .c file that was auto generated by Mathematica... the thing caused GCC to allocate ~4 gigs of ram. Simplifying the math by hand we got it down to a few hundred lines of code, and reduced the build time by three orders of magnitude.
  - TheRealPomax 10y ago
    hahaha that is amazing!
  - bhaak 10y ago
    Please, tell me you simplified the math in Mathematica and that you didn't go through the C file manually? If it's the latter, how did you keep your sanity?
    - santaclaus 10y ago
      We got a couple of big ass whiteboards and carefully went through the math. Lots of big expressions simplified to zero that Mathematica couldn't see, and we were able to factor out a lot of shared expressions.
vladdanilov 10y ago
It's not even a complete solution which is roots for a ninth-degree polynomial. Approximation with quadratics is what makes it feasible with some caveats [1]. [1] https://youtu.be/OmfliNQsk88?t=19m7s https://youtu.be/OmfliNQsk88?t=19m7s
vanderZwan 10y ago
I had a similar experience looking up algorithms for space-filling curves: https://www.reddit.com/r/ProgrammerHumor/comments/4xzi9a/oh_great_these_mathematicians_actually_provided/ https://www.reddit.com/r/ProgrammerHumor/comments/4xzi9a/oh_...
- bigger_cheese 10y ago
  Same thing when I was looking for Reinsch Spline code it's a mess. http://stackoverflow.com/questions/5072758/spline-smoothing-algorithm-in-c http://stackoverflow.com/questions/5072758/spline-smoothing-... I converted Algol code from a paper into C because I couldn't find any clean C code on the web but it's still ugly eg. here is a chunk of it: for (i=m1; i < m2; i++) //Choleksy Decomposition {(R^T)R of [(Q^T)(D^2)Q + pT]} { r1[i-1] = f*r[i-1]; r2[i-2] = g*r[i-2]; r[i] = 1/(p*b[i]+t[i]-f*r1[i-1]-g*r2[i-2]); u[i] = a[i]-r1[i-1]*u[i-1]-r2[i-2]*u[i-2]; f = p*c[i]+t1[i]-h*r1[i-1]; g = h; h = d[i]*p; }
BenFrantzDale 10y ago
I'm not sure what that wall does, but the easiest way to find intersections, roots, etc. of Bézier curves is to use their convex-hull property and recursively bisect them. The convex hull property is that the curve is within the convex hull of the control points. So given two Bézier curves, if their convex hulls don't overlap, then they don't intersect. If they do, then split them both in half and see which of those Béziers' convex hulls overlap. Recur until you reach a tolerance, then just approximate the curves by a straight line and find the intersection point (or basically equivalently, just use the average position of the control points, since at that point the line is pretty much straight).

nkoren 10y ago

What the heck...!?! I would really love for someone to explain how something like this happens.

artursapek 10y ago
Me too. It's constructing some polynomial, but the rest is lost on me.
- T-hawk 10y ago
  Looks like what you would get if you wrote out all the operations for matrix multiplication? It must be related to expressing the coordinates of the nodes as matrices. I would guess the factors of 3 and 6 arise from derivatives of the cubics: x^3 derives to 3x^2 derives to 6x.
  - wnissen 10y ago
    The intersection of two NURBS (Non-uniform Rational B-Spline) surfaces has a degree into the hundreds. Those are the surface types used in virtually all CAD modelers. For that reason the intersection is not calculated exactly but approximated to a geometric tolerance. It is harder, mathematically and algorithm-wise to calculate the approximation, but much faster to evaluate afterwards. It's amazing what magic lurks below the surface of even something as simple as "find the intersection of two polynomial functions".
    - surement 10y ago
      > NURBS (Non-uniform Rational B-Spline) I believe NURBS stands for Nobody Understands Rational B-Splines [0] :) [0] http://www.springer.com/productFlyer_978-3-540-61545-3.pdf?SGWID=0-0-1297-1491853-0 http://www.springer.com/productFlyer_978-3-540-61545-3.pdf?S... (pdf)
phkahler 10y ago
It happens because the curves are not polynomials in the normal sense. The parameter t that ranges from [0..1] is not the same between the two curves, so it is not equivalent to setting two polynomials equal to each other and solving for x, or t in this case. A solution to the intersection is equivalent to solving for the t parameter of both curves simultaneously and IIRC that ends up having an order that is the product of the orders of the curves (9th order for 2 cubics?). On second thought it's not the product, but it is definitely much higher order. That is the reason people use the matrix representations, once you get a feel for that the complexity tends to get hidden by representing an entire matrix of Nth order polynomials by a letter like B. If anyone is really interested in this stuff and wants to make a major contribution to open source software... First study these and then move on to NURBS, surfaces, and then trimmed NURBS surfaces. These are the core of all CAD software and a clean library for trimmed NURBS is sorely lacking. There are two that I'm aware of, the one in OpenCascade which is a huge ball of mud (but works) and the one in SolveSpace which doesn't scale well and has some other limitations - that is the one I feel is most worthy of someones time.
- Klapaucius 10y ago
  There is GoTools (https://github.com/SINTEF-Geometry/GoTools https://github.com/SINTEF-Geometry/GoTools) and SISL (https://github.com/SINTEF-Geometry/SISL https://github.com/SINTEF-Geometry/SISL). These are open-source geometry libraries for parametric surface modeling, that include NURBS surfaces and intersections and trimmings of such. Intersections are computed approximatively within some prescribed geometrical tolerance, and intersection curves are "marched out".
  - phkahler 10y ago
    Thanks. The links have been forwarded to people who may be interested. BTW they don't explicitly say "trimmed NURBS" but doing boolean operations implies that is the case ;-)
- blt 10y ago
  NLib is the gold standard in closed-source, right?

TheRealPomax 10y ago

This kind of feels like an export of a WolframAlpha query, exported and written out as code =D

santaclaus 10y ago
Haha yea that is some straight up Mathematica -> C codegen (damn you //cform). A couple of years ago I inherited a project with a 90 meg .c file that was auto generated by Mathematica... the thing caused GCC to allocate ~4 gigs of ram. Simplifying the math by hand we got it down to a few hundred lines of code, and reduced the build time by three orders of magnitude.
- TheRealPomax 10y ago
  hahaha that is amazing!
- bhaak 10y ago
  Please, tell me you simplified the math in Mathematica and that you didn't go through the C file manually? If it's the latter, how did you keep your sanity?
  - santaclaus 10y ago
    We got a couple of big ass whiteboards and carefully went through the math. Lots of big expressions simplified to zero that Mathematica couldn't see, and we were able to factor out a lot of shared expressions.

vladdanilov 10y ago

It's not even a complete solution which is roots for a ninth-degree polynomial. Approximation with quadratics is what makes it feasible with some caveats [1]. [1] https://youtu.be/OmfliNQsk88?t=19m7s https://youtu.be/OmfliNQsk88?t=19m7s

vanderZwan 10y ago

I had a similar experience looking up algorithms for space-filling curves: https://www.reddit.com/r/ProgrammerHumor/comments/4xzi9a/oh_great_these_mathematicians_actually_provided/ https://www.reddit.com/r/ProgrammerHumor/comments/4xzi9a/oh_...

bigger_cheese 10y ago
Same thing when I was looking for Reinsch Spline code it's a mess. http://stackoverflow.com/questions/5072758/spline-smoothing-algorithm-in-c http://stackoverflow.com/questions/5072758/spline-smoothing-... I converted Algol code from a paper into C because I couldn't find any clean C code on the web but it's still ugly eg. here is a chunk of it: for (i=m1; i < m2; i++) //Choleksy Decomposition {(R^T)R of [(Q^T)(D^2)Q + pT]} { r1[i-1] = f*r[i-1]; r2[i-2] = g*r[i-2]; r[i] = 1/(p*b[i]+t[i]-f*r1[i-1]-g*r2[i-2]); u[i] = a[i]-r1[i-1]*u[i-1]-r2[i-2]*u[i-2]; f = p*c[i]+t1[i]-h*r1[i-1]; g = h; h = d[i]*p; }

BenFrantzDale 10y ago

I'm not sure what that wall does, but the easiest way to find intersections, roots, etc. of Bézier curves is to use their convex-hull property and recursively bisect them. The convex hull property is that the curve is within the convex hull of the control points. So given two Bézier curves, if their convex hulls don't overlap, then they don't intersect. If they do, then split them both in half and see which of those Béziers' convex hulls overlap. Recur until you reach a tolerance, then just approximate the curves by a straight line and find the intersection point (or basically equivalently, just use the average position of the control points, since at that point the line is pretty much straight).