[cairo] Bug in cairo_in_stroke()

[Behdad Esfahbod]

On 03/27/2009 10:09 PM, Bertram Felgenhauer wrote:
> 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.

Neat stuff.


> 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.

Yep.  I also checked them after your patch yesterday and they matched.


> (NB: I think (a+b)*(a+b) would be clearer that a*a + b2 + _2ab.)

Makes sense.  Feel free to make that change.


> 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.

It can be improved on modern processors by using constructions like
(c < 0) ^ (p1 > 0).


> 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.

Well, for me variants 1) and 3) are easier to understand :).  At least before 
reading three lines of equations.

At any rate, feel free to replace it or add your analysis in the comments.


> 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?

There are some.  See font-options.c test for example.

behdad