From 2131cbeb6f10b52a48d2ce25cfe27417e173afa5 Mon Sep 17 00:00:00 2001 From: marha Date: Wed, 13 Oct 2010 05:49:35 +0000 Subject: pixman update 13/10/2010 --- pixman/pixman/pixman-private.h | 9 +- pixman/pixman/pixman-radial-gradient.c | 454 +++++++++++++++++++-------------- 2 files changed, 265 insertions(+), 198 deletions(-) (limited to 'pixman') diff --git a/pixman/pixman/pixman-private.h b/pixman/pixman/pixman-private.h index 8250d1f7a..e756bdbed 100644 --- a/pixman/pixman/pixman-private.h +++ b/pixman/pixman/pixman-private.h @@ -152,10 +152,11 @@ struct radial_gradient circle_t c1; circle_t c2; - double cdx; - double cdy; - double dr; - double A; + + circle_t delta; + double a; + double inva; + double mindr; }; struct conical_gradient diff --git a/pixman/pixman/pixman-radial-gradient.c b/pixman/pixman/pixman-radial-gradient.c index 6f00c4113..bc4a13412 100644 --- a/pixman/pixman/pixman-radial-gradient.c +++ b/pixman/pixman/pixman-radial-gradient.c @@ -1,3 +1,4 @@ +/* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */ /* * * Copyright © 2000 Keith Packard, member of The XFree86 Project, Inc. @@ -33,6 +34,100 @@ #include #include "pixman-private.h" +static inline pixman_fixed_32_32_t +dot (pixman_fixed_48_16_t x1, + pixman_fixed_48_16_t y1, + pixman_fixed_48_16_t z1, + pixman_fixed_48_16_t x2, + pixman_fixed_48_16_t y2, + pixman_fixed_48_16_t z2) +{ + /* + * Exact computation, assuming that the input values can + * be represented as pixman_fixed_16_16_t + */ + return x1 * x2 + y1 * y2 + z1 * z2; +} + +static inline double +fdot (double x1, + double y1, + double z1, + double x2, + double y2, + double z2) +{ + /* + * Error can be unbound in some special cases. + * Using clever dot product algorithms (for example compensated + * dot product) would improve this but make the code much less + * obvious + */ + return x1 * x2 + y1 * y2 + z1 * z2; +} + +static uint32_t +radial_compute_color (double a, + double b, + double c, + double inva, + double dr, + double mindr, + pixman_gradient_walker_t *walker, + pixman_repeat_t repeat) +{ + /* + * In this function error propagation can lead to bad results: + * - det can have an unbound error (if b*b-a*c is very small), + * potentially making it the opposite sign of what it should have been + * (thus clearing a pixel that would have been colored or vice-versa) + * or propagating the error to sqrtdet; + * if det has the wrong sign or b is very small, this can lead to bad + * results + * + * - the algorithm used to compute the solutions of the quadratic + * equation is not numerically stable (but saves one division compared + * to the numerically stable one); + * this can be a problem if a*c is much smaller than b*b + * + * - the above problems are worse if a is small (as inva becomes bigger) + */ + double det; + + if (a == 0) + { + return _pixman_gradient_walker_pixel (walker, + pixman_fixed_1 / 2 * c / b); + } + + det = fdot (b, a, 0, b, -c, 0); + if (det >= 0) + { + double sqrtdet, t0, t1; + + sqrtdet = sqrt (det); + t0 = (b + sqrtdet) * inva; + t1 = (b - sqrtdet) * inva; + + if (repeat == PIXMAN_REPEAT_NONE) + { + if (0 <= t0 && t0 <= pixman_fixed_1) + return _pixman_gradient_walker_pixel (walker, t0); + else if (0 <= t1 && t1 <= pixman_fixed_1) + return _pixman_gradient_walker_pixel (walker, t1); + } + else + { + if (t0 * dr > mindr) + return _pixman_gradient_walker_pixel (walker, t0); + else if (t1 * dr > mindr) + return _pixman_gradient_walker_pixel (walker, t1); + } + } + + return 0; +} + static void radial_gradient_get_scanline_32 (pixman_image_t *image, int x, @@ -42,118 +137,85 @@ radial_gradient_get_scanline_32 (pixman_image_t *image, const uint32_t *mask) { /* + * Implementation of radial gradients following the PDF specification. + * See section 8.7.4.5.4 Type 3 (Radial) Shadings of the PDF Reference + * Manual (PDF 32000-1:2008 at the time of this writing). + * * In the radial gradient problem we are given two circles (c₁,r₁) and - * (c₂,r₂) that define the gradient itself. Then, for any point p, we - * must compute the value(s) of t within [0.0, 1.0] representing the - * circle(s) that would color the point. - * - * There are potentially two values of t since the point p can be - * colored by both sides of the circle, (which happens whenever one - * circle is not entirely contained within the other). - * - * If we solve for a value of t that is outside of [0.0, 1.0] then we - * use the extend mode (NONE, REPEAT, REFLECT, or PAD) to map to a - * value within [0.0, 1.0]. - * - * Here is an illustration of the problem: - * - * p₂ - * p • - * • ╲ - * · ╲r₂ - * p₁ · ╲ - * • θ╲ - * ╲ ╌╌• - * ╲r₁ · c₂ - * θ╲ · - * ╌╌• - * c₁ + * (c₂,r₂) that define the gradient itself. * - * Given (c₁,r₁), (c₂,r₂) and p, we must find an angle θ such that two - * points p₁ and p₂ on the two circles are collinear with p. Then, the - * desired value of t is the ratio of the length of p₁p to the length - * of p₁p₂. + * Mathematically the gradient can be defined as the family of circles * - * So, we have six unknown values: (p₁x, p₁y), (p₂x, p₂y), θ and t. - * We can also write six equations that constrain the problem: + * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂) * - * Point p₁ is a distance r₁ from c₁ at an angle of θ: + * excluding those circles whose radius would be < 0. When a point + * belongs to more than one circle, the one with a bigger t is the only + * one that contributes to its color. When a point does not belong + * to any of the circles, it is transparent black, i.e. RGBA (0, 0, 0, 0). + * Further limitations on the range of values for t are imposed when + * the gradient is not repeated, namely t must belong to [0,1]. * - * 1. p₁x = c₁x + r₁·cos θ - * 2. p₁y = c₁y + r₁·sin θ + * The graphical result is the same as drawing the valid (radius > 0) + * circles with increasing t in [-inf, +inf] (or in [0,1] if the gradient + * is not repeated) using SOURCE operatior composition. * - * Point p₂ is a distance r₂ from c₂ at an angle of θ: + * It looks like a cone pointing towards the viewer if the ending circle + * is smaller than the starting one, a cone pointing inside the page if + * the starting circle is the smaller one and like a cylinder if they + * have the same radius. * - * 3. p₂x = c₂x + r2·cos θ - * 4. p₂y = c₂y + r2·sin θ + * What we actually do is, given the point whose color we are interested + * in, compute the t values for that point, solving for t in: * - * Point p lies at a fraction t along the line segment p₁p₂: + * length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂ + * + * Let's rewrite it in a simpler way, by defining some auxiliary + * variables: * - * 5. px = t·p₂x + (1-t)·p₁x - * 6. py = t·p₂y + (1-t)·p₁y + * cd = c₂ - c₁ + * pd = p - c₁ + * dr = r₂ - r₁ + * lenght(t·cd - pd) = r₁ + t·dr * - * To solve, first subtitute 1-4 into 5 and 6: + * which actually means * - * px = t·(c₂x + r₂·cos θ) + (1-t)·(c₁x + r₁·cos θ) - * py = t·(c₂y + r₂·sin θ) + (1-t)·(c₁y + r₁·sin θ) + * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr * - * Then solve each for cos θ and sin θ expressed as a function of t: + * or * - * cos θ = (-(c₂x - c₁x)·t + (px - c₁x)) / ((r₂-r₁)·t + r₁) - * sin θ = (-(c₂y - c₁y)·t + (py - c₁y)) / ((r₂-r₁)·t + r₁) + * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr. * - * To simplify this a bit, we define new variables for several of the - * common terms as shown below: + * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes: * - * p₂ - * p • - * • ╲ - * · ┆ ╲r₂ - * p₁ · ┆ ╲ - * • pdy┆ ╲ - * ╲ ┆ •c₂ - * ╲r₁ ┆ · ┆ - * ╲ ·┆ ┆cdy - * •╌╌╌╌┴╌╌╌╌╌╌╌┘ - * c₁ pdx cdx + * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)² * - * cdx = (c₂x - c₁x) - * cdy = (c₂y - c₁y) - * dr = r₂-r₁ - * pdx = px - c₁x - * pdy = py - c₁y + * where we can actually expand the squares and solve for t: * - * Note that cdx, cdy, and dr do not depend on point p at all, so can - * be pre-computed for the entire gradient. The simplifed equations - * are now: + * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² = + * = r₁² + 2·r₁·t·dr + t²·dr² * - * cos θ = (-cdx·t + pdx) / (dr·t + r₁) - * sin θ = (-cdy·t + pdy) / (dr·t + r₁) + * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t + + * (pdx² + pdy² - r₁²) = 0 * - * Finally, to get a single function of t and eliminate the last - * unknown θ, we use the identity sin²θ + cos²θ = 1. First, square - * each equation, (we knew a quadratic was coming since it must be - * possible to obtain two solutions in some cases): + * A = cdx² + cdy² - dr² + * B = pdx·cdx + pdy·cdy + r₁·dr + * C = pdx² + pdy² - r₁² + * At² - 2Bt + C = 0 + * + * The solutions (unless the equation degenerates because of A = 0) are: * - * cos²θ = (cdx²t² - 2·cdx·pdx·t + pdx²) / (dr²·t² + 2·r₁·dr·t + r₁²) - * sin²θ = (cdy²t² - 2·cdy·pdy·t + pdy²) / (dr²·t² + 2·r₁·dr·t + r₁²) + * t = (B ± ⎷(B² - A·C)) / A * - * Then add both together, set the result equal to 1, and express as a - * standard quadratic equation in t of the form At² + Bt + C = 0 + * The solution we are going to prefer is the bigger one, unless the + * radius associated to it is negative (or it falls outside the valid t + * range). * - * (cdx² + cdy² - dr²)·t² - 2·(cdx·pdx + cdy·pdy + r₁·dr)·t + (pdx² + pdy² - r₁²) = 0 + * Additional observations (useful for optimizations): + * A does not depend on p * - * In other words: - * - * A = cdx² + cdy² - dr² - * B = -2·(pdx·cdx + pdy·cdy + r₁·dr) - * C = pdx² + pdy² - r₁² - * - * And again, notice that A does not depend on p, so can be - * precomputed. From here we just use the quadratic formula to solve - * for t: - * - * t = (-2·B ± ⎷(B² - 4·A·C)) / 2·A + * A < 0 <=> one of the two circles completely contains the other one + * <=> for every p, the radiuses associated with the two t solutions + * have opposite sign */ gradient_t *gradient = (gradient_t *)image; @@ -161,153 +223,149 @@ radial_gradient_get_scanline_32 (pixman_image_t *image, radial_gradient_t *radial = (radial_gradient_t *)image; uint32_t *end = buffer + width; pixman_gradient_walker_t walker; - pixman_bool_t affine = TRUE; - double cx = 1.; - double cy = 0.; - double cz = 0.; - double rx = x + 0.5; - double ry = y + 0.5; - double rz = 1.; + pixman_vector_t v, unit; + + /* reference point is the center of the pixel */ + v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2; + v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2; + v.vector[2] = pixman_fixed_1; _pixman_gradient_walker_init (&walker, gradient, source->common.repeat); if (source->common.transform) { - pixman_vector_t v; - /* reference point is the center of the pixel */ - v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2; - v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2; - v.vector[2] = pixman_fixed_1; - if (!pixman_transform_point_3d (source->common.transform, &v)) return; - - cx = source->common.transform->matrix[0][0] / 65536.; - cy = source->common.transform->matrix[1][0] / 65536.; - cz = source->common.transform->matrix[2][0] / 65536.; - rx = v.vector[0] / 65536.; - ry = v.vector[1] / 65536.; - rz = v.vector[2] / 65536.; - - affine = - source->common.transform->matrix[2][0] == 0 && - v.vector[2] == pixman_fixed_1; + unit.vector[0] = source->common.transform->matrix[0][0]; + unit.vector[1] = source->common.transform->matrix[1][0]; + unit.vector[2] = source->common.transform->matrix[2][0]; + } + else + { + unit.vector[0] = pixman_fixed_1; + unit.vector[1] = 0; + unit.vector[2] = 0; } - if (affine) + if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1) { - /* When computing t over a scanline, we notice that some expressions - * are constant so we can compute them just once. Given: + /* + * Given: * - * t = (-2·B ± ⎷(B² - 4·A·C)) / 2·A + * t = (B ± ⎷(B² - A·C)) / A * * where * - * A = cdx² + cdy² - dr² [precomputed as radial->A] - * B = -2·(pdx·cdx + pdy·cdy + r₁·dr) + * A = cdx² + cdy² - dr² + * B = pdx·cdx + pdy·cdy + r₁·dr * C = pdx² + pdy² - r₁² + * det = B² - A·C * * Since we have an affine transformation, we know that (pdx, pdy) * increase linearly with each pixel, * - * pdx = pdx₀ + n·cx, - * pdy = pdy₀ + n·cy, - * - * we can then express B in terms of an linear increment along - * the scanline: + * pdx = pdx₀ + n·ux, + * pdy = pdy₀ + n·uy, * - * B = B₀ + n·cB, with - * B₀ = -2·(pdx₀·cdx + pdy₀·cdy + r₁·dr) and - * cB = -2·(cx·cdx + cy·cdy) - * - * Thus we can replace the full evaluation of B per-pixel (4 multiplies, - * 2 additions) with a single addition. + * we can then express B, C and det through multiple differentiation. + */ + pixman_fixed_32_32_t b, db, c, dc, ddc; + + /* warning: this computation may overflow */ + v.vector[0] -= radial->c1.x; + v.vector[1] -= radial->c1.y; + + /* + * B and C are computed and updated exactly. + * If fdot was used instead of dot, in the worst case it would + * lose 11 bits of precision in each of the multiplication and + * summing up would zero out all the bit that were preserved, + * thus making the result 0 instead of the correct one. + * This would mean a worst case of unbound relative error or + * about 2^10 absolute error */ - double r1 = radial->c1.radius / 65536.; - double r1sq = r1 * r1; - double pdx = rx - radial->c1.x / 65536.; - double pdy = ry - radial->c1.y / 65536.; - double A = radial->A; - double invA = -65536. / (2. * A); - double A4 = -4. * A; - double B = -2. * (pdx*radial->cdx + pdy*radial->cdy + r1*radial->dr); - double cB = -2. * (cx*radial->cdx + cy*radial->cdy); - pixman_bool_t invert = A * radial->dr < 0; + b = dot (v.vector[0], v.vector[1], radial->c1.radius, + radial->delta.x, radial->delta.y, radial->delta.radius); + db = dot (unit.vector[0], unit.vector[1], 0, + radial->delta.x, radial->delta.y, 0); + + c = dot (v.vector[0], v.vector[1], -radial->c1.radius, + v.vector[0], v.vector[1], radial->c1.radius); + dc = dot (2 * v.vector[0] + unit.vector[0], + 2 * v.vector[1] + unit.vector[1], + 0, + unit.vector[0], unit.vector[1], 0); + ddc = 2 * dot (unit.vector[0], unit.vector[1], 0, + unit.vector[0], unit.vector[1], 0); while (buffer < end) { if (!mask || *mask++) { - pixman_fixed_48_16_t t; - double det = B * B + A4 * (pdx * pdx + pdy * pdy - r1sq); - if (det <= 0.) - t = (pixman_fixed_48_16_t) (B * invA); - else if (invert) - t = (pixman_fixed_48_16_t) ((B + sqrt (det)) * invA); - else - t = (pixman_fixed_48_16_t) ((B - sqrt (det)) * invA); - - *buffer = _pixman_gradient_walker_pixel (&walker, t); + *buffer = radial_compute_color (radial->a, b, c, + radial->inva, + radial->delta.radius, + radial->mindr, + &walker, + source->common.repeat); } - ++buffer; - pdx += cx; - pdy += cy; - B += cB; + b += db; + c += dc; + dc += ddc; + ++buffer; } } else { /* projective */ + /* Warning: + * error propagation guarantees are much looser than in the affine case + */ while (buffer < end) { if (!mask || *mask++) { - double pdx, pdy; - double B, C; - double det; - double c1x = radial->c1.x / 65536.0; - double c1y = radial->c1.y / 65536.0; - double r1 = radial->c1.radius / 65536.0; - pixman_fixed_48_16_t t; - double x, y; - - if (rz != 0) + if (v.vector[2] != 0) { - x = rx / rz; - y = ry / rz; - } - else - { - x = y = 0.; - } + double pdx, pdy, invv2, b, c; - pdx = x - c1x; - pdy = y - c1y; + invv2 = 1. * pixman_fixed_1 / v.vector[2]; - B = -2 * (pdx * radial->cdx + - pdy * radial->cdy + - r1 * radial->dr); - C = (pdx * pdx + pdy * pdy - r1 * r1); + pdx = v.vector[0] * invv2 - radial->c1.x; + /* / pixman_fixed_1 */ - det = (B * B) - (4 * radial->A * C); - if (det < 0.0) - det = 0.0; + pdy = v.vector[1] * invv2 - radial->c1.y; + /* / pixman_fixed_1 */ - if (radial->A * radial->dr < 0) - t = (pixman_fixed_48_16_t) ((-B - sqrt (det)) / (2.0 * radial->A) * 65536); - else - t = (pixman_fixed_48_16_t) ((-B + sqrt (det)) / (2.0 * radial->A) * 65536); + b = fdot (pdx, pdy, radial->c1.radius, + radial->delta.x, radial->delta.y, + radial->delta.radius); + /* / pixman_fixed_1 / pixman_fixed_1 */ - *buffer = _pixman_gradient_walker_pixel (&walker, t); + c = fdot (pdx, pdy, -radial->c1.radius, + pdx, pdy, radial->c1.radius); + /* / pixman_fixed_1 / pixman_fixed_1 */ + + *buffer = radial_compute_color (radial->a, b, c, + radial->inva, + radial->delta.radius, + radial->mindr, + &walker, + source->common.repeat); + } + else + { + *buffer = 0; + } } ++buffer; - rx += cx; - ry += cy; - rz += cz; + v.vector[0] += unit.vector[0]; + v.vector[1] += unit.vector[1]; + v.vector[2] += unit.vector[2]; } } } @@ -351,12 +409,20 @@ pixman_image_create_radial_gradient (pixman_point_fixed_t * inner, radial->c2.x = outer->x; radial->c2.y = outer->y; radial->c2.radius = outer_radius; - radial->cdx = pixman_fixed_to_double (radial->c2.x - radial->c1.x); - radial->cdy = pixman_fixed_to_double (radial->c2.y - radial->c1.y); - radial->dr = pixman_fixed_to_double (radial->c2.radius - radial->c1.radius); - radial->A = (radial->cdx * radial->cdx + - radial->cdy * radial->cdy - - radial->dr * radial->dr); + + /* warning: this computations may overflow */ + radial->delta.x = radial->c2.x - radial->c1.x; + radial->delta.y = radial->c2.y - radial->c1.y; + radial->delta.radius = radial->c2.radius - radial->c1.radius; + + /* computed exactly, then cast to double -> every bit of the double + representation is correct (53 bits) */ + radial->a = dot (radial->delta.x, radial->delta.y, -radial->delta.radius, + radial->delta.x, radial->delta.y, radial->delta.radius); + if (radial->a != 0) + radial->inva = 1. * pixman_fixed_1 / radial->a; + + radial->mindr = -1. * pixman_fixed_1 * radial->c1.radius; image->common.property_changed = radial_gradient_property_changed; -- cgit v1.2.3