Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Parallel curves of cubic Béziers #80

Closed
raphlinus opened this issue Sep 2, 2022 · 9 comments
Closed

Parallel curves of cubic Béziers #80

raphlinus opened this issue Sep 2, 2022 · 9 comments

Comments

@raphlinus
Copy link
Owner

I have a better approach to making parallel (offset) curves of cubic Béziers. This is a followup to cleaner parallel curves and also moves the curve fitting work forward.

The core concept is to curve-fit the parallel curve of the source Bézier, ie find the cubic Bézier that most closely matches it. This possibility was raised in the curve fitting post, but not followed through until now. An array of numerical techniques supports this basic idea:

  • Green's theorem to compute area and moment of the offset curve
  • Gauss-Legendre quadrature to do the numerical integration
  • Quartic solving (Orellana and De Michele)
  • Arc length parametrized error estimation
  • ITP method for solving inverse arc length

The algorithm does work better when split at cusps, but this is not absolutely necessary; it can just render them.

There is a bit of refinement to the curve fitting algorithm. For the "near misses" identified in the blog post, we take the real value of complex roots. Also the d0 = 0 and d1 = 0 cases are added.

There is a JavaScript interactive implementation. It's not perfect, but is arguably good enough for a blog post.

Loose ends:

  • The current error metric is a bit slow (lots of inverse arc length)
    • Total arc length of the curve can very likely be used to prune bad candidates
    • There's a draft implementation of a faster but less robust technique; maybe can be rescued?
  • I have the idea but not the finished implementation for robust cusp finding based on intervals.
  • The numerical integration is not driven by error bounds (these are probably analogous to the ones for [arc length].
  • Real performance numbers
  • Quantitative comparison with existing approaches, including Tiller-Hanson.

Compare to Bézier primer section on offsetting. Also worth referencing: Comparing offset curve approximation methods.

@waruyama
Copy link

What precision do you need for your arc length parametrization? Is an average error of 0.25% and a maximum error of 3.5% sufficient? If that is good enough I have a quite fast solution for you.

@raphlinus
Copy link
Owner Author

Well, I'm not at present using an arc length parametrization, I've adopted the ray intersection approach from Tiller and Hanson, which generally seems to be much faster and just as accurate (though there was an issue I had to fix when the ray fails to intersect the candidate curve). That said, I am quite interested in your idea, those numbers may be good enough for at least some uses.

@waruyama
Copy link

waruyama commented Sep 13, 2022

I was more thinking of the "inverse arc length" that you mentioned above.

As you know, wt calculate the arc length, we need the integral of a nasty function that looks like this:

f(x) = sqrt(((((k4)*t+k3)*t+k2)*t+k1)*t+k0)

This is usually done with Gauss quadrature, which is reasonably fast, but for reverse arc length this is slow because it is done with bisection and repeated Gauss quadratures.

I ran into a blog post by Taco de Wolff at https://tacodewolff.nl/posts/20190525-arc-length/
He tried to create an arc length parametrization using a cubic polynomial. This looked like a good idea for inverse path length, too, until I realized that a cubic polynomial cannot easily be inverted and is not necessarily monotonic, so we can get cases with the funny effect of decreasing curve length while t is increasing.

But I found another nice solution. I split the length function at its inflection points, which can be found quite easily and each segment is approximated with a (drommroll) quadratic Bézier curve. The reason for this is that we can make sure that the curve is monotonic and mapping from t to as and from s to t both require only solving one quadratic equation, which is blazing fast. The end points of the quadratic Bézier curves [ti, si] are given by the t values and arc lengths at the inflection points. The slope at these points (the function f(x) above is the derivative of the length function, so that is the slope). This is all we need to find the quadratics' control points.

In other words, arc length and inverse arc length require only solving one quadratic equation. The whole initial setup is solving one cubic equation to get the inflection points, n+1 Gauss quadratures (n is the number of inflection points), and then finding the quadratic Bézier curves, which is just some addition, multiplication and division.

I will try to post some (pseudo) code.

@waruyama
Copy link

waruyama commented Sep 13, 2022

Here we go:

https://github.com/waruyama/Cuve-Length-Parametrization/blob/main/CurveLengthParametrization.js

At the top you can see the import of some functions that I cannot extract right now, but they are fairly simple.

solveQuadratic and solveCubic do just that, solve a quadratic or cubic equation. The have the same signature as these functions in paper.js

The integrate(fn, n, ...ts) function performs a Gauss integration for the function fn, using n quadrature points, for the specified t values and returns an array of integration intervals [0,t1], [t1,t2] ... [tm,1].

@raphlinus
Copy link
Owner Author

Fun! This might be good enough for dashing and so on.

One minor point I'd make - bisection is a pretty crude approach for inverting arc length. Much better to use ITP method. Even so, solving a quadratic Bézier will be much faster and has a pleasing symmetry between forward and inverse solve.

@waruyama
Copy link

Yes, it is good for dashing, but also it can also be used to find evenly spaced sample points for estimating the error of an offset curve, especially with increasing number of sample points.

I did not know the ITP method. I would have thought that in those very basic fields of mathematics everything has been settled for decades or centuries (Regula Falsi was my method of choice, invented 2000 years ago), but it is amazing that for stuff like quartic root finding and bracketed root approximation they still find ways to improve (even a lot in some cases).

@waruyama
Copy link

BTW, I think there is a bug in your ITP code (and in the Wikipedia article). In the current form it fails in many cases badly, but after applying the fix below it works all the time, at least in my approx. 10000 test runs.

Your code:

    if (yitp > 0) {
      b = xitp;
      yb = yitp;
    } else if (yitp < 0) {
      a = xitp;
      ya = yitp;
    } else {
      return xitp;
    }

Fixed:

    if (yitp * yb > 0) {
      b = xitp;
      yb = yitp;
    } else if (yitp * yb < 0) {
      a = xitp;
      ya = yitp;
    } else {
      return xitp;
    }

You can see in this ITP implementation that it should probably be done this way:

https://github.com/paulnorthrop/itp/blob/main/src/itp_c.cpp#L161

@waruyama
Copy link

I asked the author of the ITP implementation and promptly got an answer:

paulnorthrop/itp#2

In short, the Wikipedia article is correct, but it requires ya < 0 < yb.

@raphlinus
Copy link
Owner Author

@waruyama This is documented, and something of a style choice. Being able to handle either sign would be more friendly, but slightly slower. For many of the applications (inverse arc length is the primary use in kurbo at the moment), I know the signs, so this wouldn't help.

Thanks for the interest and the link to Paul's project. I know of few other implementations.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants