[cairo] Bug in cairo_in_stroke()

[Bertram Felgenhauer]

Bertram Felgenhauer wrote:
> The fault was actually in the extents calculation. _cairo_spline_bound
> had a wrong sign:
[snip] 
> I'll write more about the math here later.

1. The current code.
--------------------

The code in question is the FIND_EXTREMES(a, b, c) macro in
cairo-spline.c. It finds the roots of the quadratic polynomial
p(t) = a*t² + 2*b*t + c. The interesting case is when the p(t)
has two distinct roots, i.e. when a != 0 and b² - a*c > 0.
(b² - a*c is the discriminant of the polynomial).

The code for that case looks like this:

    double b2 = b * b;
    double delta = b2 - a * c;

    /* Note that at this point, a != 0 and delta > 0 */
    
    cairo_bool_t feasible;
    double _2ab = 2 * a * b;
    /* We are only interested in solutions t that satisfy 0<t<1
     * here.  We do some checks to avoid sqrt if the solutions
     * are not in that range.  The checks can be derived from:
     *
     *   0 < (-b±√delta)/a < 1
     */
    if (_2ab >= 0)
        feasible = delta > b2 && delta < a*a + b2 + _2ab;
    else if (-b / a >= 1)
        feasible = delta < b2 && delta > a*a + b2 + _2ab;
    else
        feasible = delta < b2 || delta < a*a + b2 + _2ab;
   
    if (unlikely (feasible)) {
        double sqrt_delta = sqrt (delta);
        ADD ((-b - sqrt_delta) / a);
        ADD ((-b + sqrt_delta) / a);
    }

The last part just applies the quadratic formula. ADD() will check
that the roots that were found are indeed in (0, 1).

But let's check the calculation of 'feasible'. That variable must be
set if p(t) has a root in (0, 1).

I will use the notation x # y # z to signify that x < y < z or
x > y > z, i.e. that y lies between x and y. Note that 0 # y # z
implies 0 < x² < y².

Start with

  0 < (-b ± √delta)/a < 1.

We have a != 0, so we can multiply by a. This may flip the comparisons.

  0 # -b ± √delta # a
  b # ± √delta # a+b                                   (1)

Now there are three cases:

  a) 0 # b # a+b or 0 = b. In that case, squaring (1) gives
         b² < delta < (a+b)²
     We can detect this case by checking whether a*b >= 0.

  b) b # a+b # 0 or a+b = 0. In that case,
         b² > delta > (a+b)²
     We can detect this case using a*b < 0 and |b| >= |a|, or -b/a >= 1.

  c) b # 0 # a+b. In that case, (1) is equivalent to
         b # ±√delta # 0  or  0 # ±√delta # a+b
     Squaring gives
         delta < b²  or  delta < (a+b)²

These checks agree with the code.
(NB: I think (a+b)*(a+b) would be clearer that a*a + b2 + _2ab.)


2. Alternatives for calculation of 'feasible'
---------------------------------------------

0)  the current code, as discussed above.

1)  feasible = 1;

    This is admissible, because ADD() does its own checks, but defeats
    the purpose of avoiding sqrt calculations.

2)  feasible = fabs(3*a + 8*b + 8*c) < fabs(a) + fabs(4*a + 8*b);

    This is based on writing 

        p(t) = a*t² + 2*b*t + c
             = a*(t-1/2)² + (2*b+a)(t-1/2) + a/4 + b + c
             = a*((t-1/2)²-1/8) + (2*b+a)(t-1/2) + 3/8 a + b + c

    If p(t) = 0, and 0 < t < 1, then
        |3/8 a + b + c| = |a*((t-1/2)²-1/8) + (2*b+a)(t-1/2)|
                       <= |a*((t-1/2)²-1/8)| + |(2*b+a)(t-1/2)|
                       <  |a/8| + |(2b+a)/2|

    which, after multiplying by 8 gives the test above.

    The test may miss cases when there are no roots in (0,1),
    but it's quite cheap.

3)  double p1 = a + 2*b + c;
    if (c < 0 && p1 > 0)
        feasible = TRUE;
    else if (c > 0 && p1 < 0)
        feasible = TRUE;
    else if (a < 0)
        feasible = b > 0 && b < -a && c < 0;
    else
        feasible = b < 0 && b > -a && c > 0;

    We have p(0) = c and p(1) = a + 2*b + c. By the intermediate value
    theorem, we have a root in (0,1) if the signs of p(0) and p(1)
    differ. This case is handled by the first two ifs.

    If p(0) and p(1) have the same sign, then there can only be a
    root if the extremum of p(t) lies in (0,1), i.e. if
    0 < -b/a < 1, or 0 # b # -a.

    We already know that p has two roots, so p(-b/a) and a have
    different signs. There is a root in (0,1) if the sign of the
    extremum and the signs of p(0) differ.

    If a < 0, the condition become 0 < b < -a and p(0) = c < 0.

    For a > 0, the comparisons get turned around.

    The huge penalty on mispredicted jumps makes this version
    unattractive on most modern processors. But it may be useful
    on embedded hardware with very poor FPUs.

Experimentally, using Soeren's curve in various orientation as the
test case, variant 2) generates about 1% false positives. A similar
test with circles of radius 50 produced no false positives. I don't
know how representative these two examples are.


3. Conclusion
-------------

The current code is correct and works. It could be a bit more readable;
replacing a*a + b2 + _2ab by (a + b)*(a + b) would help.

Personally, I'm fond of variant 2) above, because it's quite effective
and easier to understand than the current code.

Bertram

P.S. I have the beginnings of a test case for these cairo_in_foo bugs;
    see attachment. Is there an established way of incorporating non-
    visual tests into the testsuite?
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cairo-in-foo-test.c
Type: text/x-csrc
Size: 3039 bytes
Desc: not available
Url : http://lists.cairographics.org/archives/cairo/attachments/20090328/33d9132c/attachment.c