[cairo] Bug in cairo_in_stroke()
Behdad Esfahbod
behdad at behdad.org
Fri Mar 27 21:51:37 PDT 2009
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
More information about the cairo
mailing list