bezier shape-fill: problem with finding roots for intersection test

Question

I adapted this shadertoy into a Unity shader and have noticed that, while it generally works as expected, there are certain artifacts that appear, as I will further explain. I believe they are related to the calculation of roots in the Bezier curves.

Below is a visualization of the artifact. First picture shows the Bezier shape with artifact, second picture shows the Bezier shape in earlier animation frames, working correctly:

I seemed to have narrowed down the origin of the artifact to the function that calculates the Bezier roots. When I only visualize the shape outline, there are no artifacts. Artifacts are introduced when applying the method that fills in the shape using the even-odd intersection test.

The Bezier shape is animated: In other words, the Vector2 positions that define the control points of the Bezier curve change over time (they linearly interpolate between two defined positions). At most frames behavior is expected, however, at certain frames (or positions) the artifacts appear.

The fill-rule basically counts intersections between the Bezier curve and horizontal lines at each uv.y coordinate, where a given horizontal line acts as the x-axis for the Bezier function and the roots of the function represents points that intersect with this x-axis. However, some of the roots will be complex numbers (i.e. they have a real and imaginary component, a + bi), which I think might cause the issue. However, the function that calculates the roots (starting at line 74) evidently discards complex roots, because I read the annotations where they link their sources and one of the sources is clearly where the shadertoy author got this function from: the cited article describes how this function discards complex roots. See "Implementing Cardano's algorithm for finding all real roots" about 1/4 the way down.

I am scratching my head. Does anyone reading have experience with Bezier functions? Ever seen artifacts like this? Would really appreciate any insight.

Thanks!

I have already tried addressing "near misses" and double-counting of segments at segments' beginning and ending, since this is where segments may overlap (e.g. the point that defines the end of one Bezier segment is the exact same position that defines the beginning of the next segment). This did not solve the problem.

Answer 1

The problem is definitely the computation of the real roots of the cubic polynomial.

One first misunderstanding is the scale of the intermediate values. In the segment where the errors are produced, the coefficients of the expanded and normalized polynomial are in the range of 3*1500. For instance at y = 4.84400272369 the values are

a = -1117.05571565802, b = 1706.85975024015, c = -1828.09125840538

for the polynomial t^3 + 3*a*t^2 + 3*b*t + c.

This results in coefficients in the million and billion range for the coefficients of the reduced polynomial, in my version s^3 - 3*p*s +q, in the example

p = 1246106.61213401, q = -2782036197.45853

The discriminant then gets an even larger value, for the example -178627039216.052. This is for an y level close to a double root, meaning this is a "small" discriminant value.

The three real roots here are

0: 3349.63861387448, 1: 0.568446220187752, 2: 0.960086879399569

This combination of large and small roots will have an influence on the precision of the roots. The values cited above are for a computation in double precision.

---

The implementation I used for the values above employs the typical changes to compute more accurate solutions, to avoid catastrophic cancellation.

• use the atan2 method to compute angles (overloaded atan in the shader math functions)
• stable computation of solutions of quadratic equations,
• stable combination of the cubic roots, use binomial theorem A^3+B^3=(A+B)*(A^2-A*B+B^2) if A*B<0.

This removes the errors in the shape of the object. The occasional flicker of one or two scan lines can be eliminated by applying one Newton step to the roots in and close to the unit interval.

---

Better still would be to avoid this escalation in the scale of the numbers.

To that end transform the input polynomial

A*(1-t)^3 + 3*B*(1-t)^2*t + 3*C*(1-t)*t^2 + D*t^3

in the case abs(A) < abs(D) to

A/D + 3*(B/D)*s + 3*(C/D)*s^2 + s^3  with  s = t/(1-t),  t = s/(1+s)

and in the opposite case to

s^3 + 3*(B/A)*s^2 + 3*(C/A)*s + D/A  with  s = (t-1)/t,  t = 1/(1+s)

---

In code, for the call in the counting procedure

    vec3 roots = vec3(-2,-2,-2);
 ...
    if(abs(sA.y) < abs(sD.y)) {
        cubicRoots(sC.y/sD.y, sB.y/sD.y, sA.y/sD.y, roots);
        roots = roots/(1.0+roots);
    } else {
        cubicRoots(sB.y/sA.y, sC.y/sA.y, sD.y/sA.y, roots);
        roots = 1.0/(1.0+roots);
    }            

and the modified Cardano procedure itself. Note that there are no preparatory coefficient manipulations now.

void cubicRoots(float a, float b, float c, inout vec3 roots){

    // input stands for  t^3 + 3*a*t^2 + 3*b*t + c = 0
    // insert t = s-a
    
    //  s^3-3as^2+3a^2s-a^3 + 3a*(s^2-2as+a^2) + 3b*(s-a) + c
    //  s^3 + 3*(-a^2+b)*s + (2*a^3-3*a*b+c)
    // Precompute a^2 and other powers
    float a2 = a * a;

    float p = a2-b;
    float p3 = p * p * p;

    float q = (2.0*a2*a - 3.0 * a * b + c);
    float q2 = 0.5 * q;

    // equation now s^3 - 3*p*s + q = 0
    // Cardano trick s = u+v
    // u^3+v^3+q + 3*(u+v)*(u*v-p) = 0
    // u^6+p^3+q*u^3 = 0
    // (u^3+q2)^2 = q2^2-p3 = discriminant
    
    // q^2 - 4*p^3 = (4a^6+9a^2b^2+c^2-12a^4b-6abc+4a^3c) - 4*(a^6-3a^4b+3a^2b^2-b^3)
    float discriminant = 0.25*(-3.0*a2*b*b+c*c-6.0*a*b*c)+a2*a*c+b*b*b;

    if(discriminant <= 0.0){
        // Three possible real roots s=2*r*cos(phi) 
        // => r^3 * (cos(3*phi)+3*cos(phi)) - 3*p*r*cos(phi) + q2 = 0
        // r^2=p,  cos(3*phi)=-q2/sqrt(p3), sin(3*phi) = sqrt(-discriminant)/sqrt(p3)

        float r = sqrt(p);
        float phi = atan(sqrt(-discriminant), -q2);
     
        roots.x = phi;
        roots.y = phi + TWO_PI;
        roots.z = phi - TWO_PI;
        
        roots = 2.0*r*cos(roots*O3) - a;      

    } else if(discriminant < -1e-6) {
        // Three real roots, but two of them are equal
        // s^3 - 3*p*s + q = (s+u)^2 * (s-2*u) = s^3 - 3*u^2*s - 2*u^3               
        float u1 = -cuberoot(q2);
        roots.x = 2.0 * u1;
        roots.y = roots.z= -u1;
        roots -= a;
    } else {  // discriminant > 0, one real root, two complex roots
        float sd = sqrt(discriminant);
        if(q2<0.0) sd = -sd;
        float u = -cuberoot(q2+sd);
        float v = p/u;
        roots.x = (p>=0.0)?(u+v-a):(-q/(u*u-p+v*v)-a);
    }
}