Why splitting bezier for length calculation?

[waruyama]

I have been reading your blog and part of you code to learn more about arc length computation for cubic bezier curves and found that in cubicBezierLength() you are splitting the curve before applying the Gauss quadrature to calculate the length of each piece of curve. As I understand it this should not be necessary, because you can choose t0 and t1 (or a and b in your case) to calculate the length of each piece without splitting.

Is this a possible optimization for your code?

[tdewolff]

I think it gave slightly better results by splitting the Béziers, but not too much. It really depends on the curve, as some curves are very hard to calculate accurately. Your mileage may vary, but it should work fine without splitting. This needs better testing to see if it really improves accuracy.

[waruyama]

Thanks. I have also tried all kind of things to improve accuracy for the worst case curves. Someone recommended to use Romberg's method and I actually implemented it, but results were worse (by a factor 2) than Gauss quadrature.
One thing that actually improves accuracy quite a bit is splitting the curve (or setting intermediate t values) at "fake curvature extrema", where x'*x'' + y'*y'' = 0. These roots are not the real curvature extrema, but often very close and if there is a the cusp they include it. Even though calculating the roots of this is not exactly slow, it takes as long as the whole Gauss quadrature with n=16. This shows how fast and awesome Gauss quadrature is, but I wish there was a simple way to improve accuracy for edge cases.
Anyway, thanks for your great blog post, it helped me a lot for understanding how to approach arc length approximation.

[devshgraphicsprogramming]

We're porting a huge CAD application to Vulkan and they use a lot of striped "pens", and this is what we can tell you:

1. quadratic beziers are 2x or more faster than cubics to deal with even after you account for the fact that you need multiple quadratics to fit a single cubic
2. actual analytical solution to quad arclength is 2^-19 accurate or more, and comparatively fast to LG
3. Newton Raphson to invert arclength is stupid fast and accurate, 8 iterations get you to 2^-24 of the arclength you want to invert
4. beziers need to be clipped to screen to avoid IEEE754 issues on a 4k screen when determining feathering/coverage and stipple pattern from arc-length
5. drawing beziers "straight from SDF" instead of turning them into poly-lines is dope and 100% pixel perfect

[tdewolff]

Hi @devshgraphicsprogramming , thanks for the advice. Is this general advice or did you stumble upon a problem of accuracy / performance?

So to sum up: you advice to convert cubics to quadratics and use the analytical implementation for the length that I already have implemented (https://github.com/tdewolff/canvas/blob/master/path_util.go#L412)?

I don't understand point 4. And how do you mean straight from SDF, is that using OpenGL shaders? Or something built-in the software rasterizer?

[devshgraphicsprogramming]

Hi @devshgraphicsprogramming , thanks for the advice. Is this general advice or did you stumble upon a problem of accuracy / performance?

We've measured the per-pixel math for just the SDF computation to be 2x slower to perform than for 2 quad beziers that specify a similar curve + overlap. The arclength functions are also much faster if you want accurate arclength and not LG quadrature with 6 or 8 coefficients.

Furthermore when dealing with beziers the SDF is distance to closest point. For a cubic the SDF is numerical root finding of 5 roots, with either mediocre accuracy or lots of ALU waste, also most root finders' convergence is sensitive to input function characteristics like the magnitude of the slope (Newton) or the slope not being 0 (Halley). Quadratic closest point only requires a bit of special care to account for when the control points are close to becoming a straight like.

The landscape has changed a lot since 2000s when a lot was written about beziers and their rendering, we now have stupidly fast inverse-trig, trig and special functions, even on GPUs!

Accuracy of arclength (for stipple patterns) and SDF (if you want pixel-perfect AA/coverage and feathering) becomes an issue when you zoom in a lot, especially relevant with a 4k or 8k display.

Anyway because of FP32 precision and 4k displays, you'll need to clip your beziers' control points dynamically at runtime to the screen.

Solving a quadratic equation for the clipping is easy and you always get 0-2 separate curve segments as the result of clipping against an axis-slab, clipping a cubic requires solving a cubic and can produce 0-3 curve segments.

So to sum up: you advice to convert cubics to quadratics and use the analytical implementation for the length that I already have implemented (https://github.com/tdewolff/canvas/blob/master/path_util.go#L412)?

Yes, basically this: https://www.youtube.com/watch?v=8l1SJAx3yoM
Except that we do it with adaptive number of quadratics not always 2, as we fit not only cubics, but also elliptical arcs, etc. and we have an error bound whereas Yuksel does not.

You implementation has a numerical stability issue, you don't check for b*b-4ac being close to zero (also happens when bezier is close to folding on itself) which will also cover the a=0 case as well as many others. Here's our code

float calcArcLength_Analytic(in float t, in vec2 p0, in vec2 p1, in vec2 p2)
{
    vec2 A = p0 - 2.0*p1 + p2;
    vec2 B = 2.0*(p1 - p0);
    vec2 C = p0;
    float lenA2 = dot(A,A);
    float AdotB = dot(A,B);
        
    float a = 4.0 * lenA2;
    float b = 4.0 * AdotB;
    float c = dot(B,B);
    
    float b_over_4a = AdotB/a;
    
    float lenTan = sqrt(t*(a*t+b)+c);
    float retval = 0.5f*t*lenTan;
    // we skip this because when |a| -> we have += 0/0 * 0 here resulting in NaN
    if (lenA2>=exp2(-23.f))
        retval += b_over_4a*(lenTan-sqrt(c));
    
    // sin2 multiplied by length of A and B
    float det_over_16 = AdotB*AdotB-lenA2*c;
    // because `b` is linearly dependent on `a` this will also ensure `b_over_4a` is not NaN, ergo `a` has a minimum value
    if (det_over_16>=exp2(-23.f))
    {
    // TODO: finish
        //retval += det_over_16*...;
    }
    return retval;
}

I can ping you when the improved version is in our repository.

I don't understand point 4. And how do you mean straight from SDF, is that using OpenGL shaders? Or something built-in the software rasterizer?

Well that would be point (5), point (4) is about knowing the arclength and manipulating closest points using inverse arclength per-pixel.

Yes pixel-shaders, you use the SDF value at a pixel minus the half-width of the stroke to determine your pixel's alpha.

This means you only construct very crude bounding polygons around the beziers.

How do you stop "overblending" neighbouring quad beziers on a piecewise line like this, you may ask?
https://i.ibb.co/pW4xwHN/image.png

We use a very weird algorithm that leverages fragment shader interlock (which AMD does not have). If someone needs to care about platforms that don't have this, they need to very diligently construct mitres on the cages, also shear them off in order to stop "exploding mitres" to infinity, also account of cases when curvature of any bezier is soo great that it exceeds the quarter-width of the stroke. We can just fall back to using an OBB in the last case and not implement mitres because the ends of our beziers can overlap without overblending.

Here's a giv that shows (WIP) the bounding polygons and the IEEE754 accuracy issues if you don't handle special cases in the arclength calculation.
https://cdn.discordapp.com/attachments/1115628593922445342/1139116684670668821/Animation.gif