diff options
| -rw-r--r-- | pixman/pixman/pixman-private.h | 9 | ||||
| -rw-r--r-- | pixman/pixman/pixman-radial-gradient.c | 454 |
2 files changed, 265 insertions, 198 deletions
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 <math.h>
#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;
|