Fix floating point error in bezier bounding box calculation

[bouhek]

Description

During my recent work with this library, I encountered an issue related to the calculation of bounding boxes for certain cubic Bezier curves. This problem appeared to stem from a floating point error that occurs during the computation of the denominator, which then affects subsequent conditions.

To resolve this, I have implemented a solution in this PR: the introduction of a FLOAT_EPSILON constant. This constant serves as a threshold for acceptable error in floating point calculations, effectively addressing the issue at hand.

To demonstrate the efficacy of this fix, I've added the zero_denominator_bezier_curve test. Included below are before-and-after images to illustrate the improvement. In these images, the Bezier curve is depicted by a black line, while the red rectangle represents the bounding box as calculated. This visual comparison clearly shows the enhancement brought about by this PR.

Before the fix

After the fix

Other test cases added

Test start_end_bbox_bezier_curve

Test general_bezier_curve

[bouhek]

Hey, @mathandy, could you have a look at this when you have time? Thanks!

[Brishen]

I found this error as well. I would propose this method instead to avoid setting a constant.

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) < ulp(a[0]) and abs(denom) > 0:
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)

[bouhek]

I found this error as well. I would propose this method instead to avoid setting a constant.

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) < ulp(a[0]) and abs(denom) > 0:
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)

Hey @Brishen, thanks for the suggestion! I like the idea of getting rid of the constant, but I feel like your solution wouldn't work for the arguments where denominator actually equals to 0 without the floating point error due to the and abs(denom) > 0 in your condition is bigger than ulp(a[0]). I'll try to revisit this tomorrow and add a test case for this too.

[Brishen]

@bouhek I think I got my conditions slightly wrong. It should probably be abs(denom) > ulp(a[0]) (and a[x] should probably be chosen to be the largest value).

from math import ulp

def bezier_real_minmax(p):
    """returns the minimum and maximum for any real cubic bezier"""
    local_extremizers = [0, 1]
    if len(p) == 4:  # cubic case
        a = [p.real for p in p]
        denom = a[0] - 3*a[1] + 3*a[2] - a[3]
        if abs(denom) > ulp(a[0]):
            delta = a[1]**2 - (a[0] + a[1])*a[2] + a[2]**2 + (a[0] - a[1])*a[3]
            if delta >= 0:  # otherwise no local extrema
                sqdelta = sqrt(delta)
                tau = a[0] - 2*a[1] + a[2]
                r1 = (tau + sqdelta)/denom
                r2 = (tau - sqdelta)/denom
                if 0 < r1 < 1:
                    local_extremizers.append(r1)
                if 0 < r2 < 1:
                    local_extremizers.append(r2)
            local_extrema = [bezier_point(a, t) for t in local_extremizers]
            return min(local_extrema), max(local_extrema)

    # find reverse standard coefficients of the derivative
    dcoeffs = bezier2polynomial(a, return_poly1d=True).deriv().coeffs

    # find real roots, r, such that 0 <= r <= 1
    local_extremizers += polyroots01(dcoeffs)
    local_extrema = [bezier_point(a, t) for t in local_extremizers]
    return min(local_extrema), max(local_extrema)