diff options
Diffstat (limited to 'pixman')
| -rw-r--r-- | pixman/pixman/pixman-radial-gradient.c | 931 |
1 files changed, 471 insertions, 460 deletions
diff --git a/pixman/pixman/pixman-radial-gradient.c b/pixman/pixman/pixman-radial-gradient.c index 63c712cc2..ecbca6f63 100644 --- a/pixman/pixman/pixman-radial-gradient.c +++ b/pixman/pixman/pixman-radial-gradient.c @@ -1,460 +1,471 @@ -/* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */
-/*
- *
- * Copyright © 2000 Keith Packard, member of The XFree86 Project, Inc.
- * Copyright © 2000 SuSE, Inc.
- * 2005 Lars Knoll & Zack Rusin, Trolltech
- * Copyright © 2007 Red Hat, Inc.
- *
- *
- * Permission to use, copy, modify, distribute, and sell this software and its
- * documentation for any purpose is hereby granted without fee, provided that
- * the above copyright notice appear in all copies and that both that
- * copyright notice and this permission notice appear in supporting
- * documentation, and that the name of Keith Packard not be used in
- * advertising or publicity pertaining to distribution of the software without
- * specific, written prior permission. Keith Packard makes no
- * representations about the suitability of this software for any purpose. It
- * is provided "as is" without express or implied warranty.
- *
- * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS
- * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
- * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
- * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN
- * AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING
- * OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS
- * SOFTWARE.
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <stdlib.h>
-#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)
- {
- double t;
-
- if (b == 0)
- return 0;
-
- t = pixman_fixed_1 / 2 * c / b;
- if (repeat == PIXMAN_REPEAT_NONE)
- {
- if (0 <= t && t <= pixman_fixed_1)
- return _pixman_gradient_walker_pixel (walker, t);
- }
- else
- {
- if (t * dr > mindr)
- return _pixman_gradient_walker_pixel (walker, t);
- }
-
- return 0;
- }
-
- 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 uint32_t *
-radial_get_scanline_narrow (pixman_iter_t *iter, 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.
- *
- * Mathematically the gradient can be defined as the family of circles
- *
- * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
- *
- * 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].
- *
- * 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.
- *
- * 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.
- *
- * 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:
- *
- * 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:
- *
- * cd = c₂ - c₁
- * pd = p - c₁
- * dr = r₂ - r₁
- * lenght(t·cd - pd) = r₁ + t·dr
- *
- * which actually means
- *
- * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
- *
- * or
- *
- * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
- *
- * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
- *
- * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
- *
- * where we can actually expand the squares and solve for t:
- *
- * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
- * = r₁² + 2·r₁·t·dr + t²·dr²
- *
- * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
- * (pdx² + pdy² - r₁²) = 0
- *
- * 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:
- *
- * t = (B ± ⎷(B² - A·C)) / A
- *
- * 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).
- *
- * Additional observations (useful for optimizations):
- * A does not depend on p
- *
- * 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
- */
- pixman_image_t *image = iter->image;
- int x = iter->x;
- int y = iter->y;
- int width = iter->width;
- uint32_t *buffer = iter->buffer;
-
- gradient_t *gradient = (gradient_t *)image;
- radial_gradient_t *radial = (radial_gradient_t *)image;
- uint32_t *end = buffer + width;
- pixman_gradient_walker_t walker;
- 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, image->common.repeat);
-
- if (image->common.transform)
- {
- if (!pixman_transform_point_3d (image->common.transform, &v))
- return iter->buffer;
-
- unit.vector[0] = image->common.transform->matrix[0][0];
- unit.vector[1] = image->common.transform->matrix[1][0];
- unit.vector[2] = image->common.transform->matrix[2][0];
- }
- else
- {
- unit.vector[0] = pixman_fixed_1;
- unit.vector[1] = 0;
- unit.vector[2] = 0;
- }
-
- if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
- {
- /*
- * Given:
- *
- * t = (B ± ⎷(B² - A·C)) / A
- *
- * where
- *
- * 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·ux,
- * pdy = pdy₀ + n·uy,
- *
- * 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
- */
- 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],
- -((pixman_fixed_48_16_t) radial->c1.radius),
- v.vector[0], v.vector[1], radial->c1.radius);
- dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0],
- 2 * (pixman_fixed_48_16_t) 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++)
- {
- *buffer = radial_compute_color (radial->a, b, c,
- radial->inva,
- radial->delta.radius,
- radial->mindr,
- &walker,
- image->common.repeat);
- }
-
- 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++)
- {
- if (v.vector[2] != 0)
- {
- double pdx, pdy, invv2, b, c;
-
- invv2 = 1. * pixman_fixed_1 / v.vector[2];
-
- pdx = v.vector[0] * invv2 - radial->c1.x;
- /* / pixman_fixed_1 */
-
- pdy = v.vector[1] * invv2 - radial->c1.y;
- /* / pixman_fixed_1 */
-
- b = fdot (pdx, pdy, radial->c1.radius,
- radial->delta.x, radial->delta.y,
- radial->delta.radius);
- /* / pixman_fixed_1 / pixman_fixed_1 */
-
- 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,
- image->common.repeat);
- }
- else
- {
- *buffer = 0;
- }
- }
-
- ++buffer;
-
- v.vector[0] += unit.vector[0];
- v.vector[1] += unit.vector[1];
- v.vector[2] += unit.vector[2];
- }
- }
-
- iter->y++;
- return iter->buffer;
-}
-
-static uint32_t *
-radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask)
-{
- uint32_t *buffer = radial_get_scanline_narrow (iter, NULL);
-
- pixman_expand ((uint64_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width);
-
- return buffer;
-}
-
-void
-_pixman_radial_gradient_iter_init (pixman_image_t *image, pixman_iter_t *iter)
-{
- if (iter->flags & ITER_NARROW)
- iter->get_scanline = radial_get_scanline_narrow;
- else
- iter->get_scanline = radial_get_scanline_wide;
-}
-
-PIXMAN_EXPORT pixman_image_t *
-pixman_image_create_radial_gradient (pixman_point_fixed_t * inner,
- pixman_point_fixed_t * outer,
- pixman_fixed_t inner_radius,
- pixman_fixed_t outer_radius,
- const pixman_gradient_stop_t *stops,
- int n_stops)
-{
- pixman_image_t *image;
- radial_gradient_t *radial;
-
- image = _pixman_image_allocate ();
-
- if (!image)
- return NULL;
-
- radial = &image->radial;
-
- if (!_pixman_init_gradient (&radial->common, stops, n_stops))
- {
- free (image);
- return NULL;
- }
-
- image->type = RADIAL;
-
- radial->c1.x = inner->x;
- radial->c1.y = inner->y;
- radial->c1.radius = inner_radius;
- radial->c2.x = outer->x;
- radial->c2.y = outer->y;
- radial->c2.radius = outer_radius;
-
- /* 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;
-
- return image;
-}
-
+/* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */ +/* + * + * Copyright © 2000 Keith Packard, member of The XFree86 Project, Inc. + * Copyright © 2000 SuSE, Inc. + * 2005 Lars Knoll & Zack Rusin, Trolltech + * Copyright © 2007 Red Hat, Inc. + * + * + * Permission to use, copy, modify, distribute, and sell this software and its + * documentation for any purpose is hereby granted without fee, provided that + * the above copyright notice appear in all copies and that both that + * copyright notice and this permission notice appear in supporting + * documentation, and that the name of Keith Packard not be used in + * advertising or publicity pertaining to distribution of the software without + * specific, written prior permission. Keith Packard makes no + * representations about the suitability of this software for any purpose. It + * is provided "as is" without express or implied warranty. + * + * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS + * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND + * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY + * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN + * AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING + * OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS + * SOFTWARE. + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif +#include <stdlib.h> +#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: + * - discr 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 sqrtdiscr; + * if discr 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 discr; + + if (a == 0) + { + double t; + + if (b == 0) + return 0; + + t = pixman_fixed_1 / 2 * c / b; + if (repeat == PIXMAN_REPEAT_NONE) + { + if (0 <= t && t <= pixman_fixed_1) + return _pixman_gradient_walker_pixel (walker, t); + } + else + { + if (t * dr > mindr) + return _pixman_gradient_walker_pixel (walker, t); + } + + return 0; + } + + discr = fdot (b, a, 0, b, -c, 0); + if (discr >= 0) + { + double sqrtdiscr, t0, t1; + + sqrtdiscr = sqrt (discr); + t0 = (b + sqrtdiscr) * inva; + t1 = (b - sqrtdiscr) * inva; + + /* + * The root that must be used if the biggest one that belongs + * to the valid range ([0,1] for PIXMAN_REPEAT_NONE, any + * solution that results in a positive radius otherwise). + * + * If a > 0, t0 is the biggest solution, so if it is valid, it + * is the correct result. + * + * If a < 0, only one of the solutions can be valid, so the + * order in which they are tested is not important. + */ + 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 uint32_t * +radial_get_scanline_narrow (pixman_iter_t *iter, 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. + * + * Mathematically the gradient can be defined as the family of circles + * + * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂) + * + * 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]. + * + * 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. + * + * 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. + * + * 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: + * + * 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: + * + * cd = c₂ - c₁ + * pd = p - c₁ + * dr = r₂ - r₁ + * lenght(t·cd - pd) = r₁ + t·dr + * + * which actually means + * + * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr + * + * or + * + * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr. + * + * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes: + * + * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)² + * + * where we can actually expand the squares and solve for t: + * + * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² = + * = r₁² + 2·r₁·t·dr + t²·dr² + * + * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t + + * (pdx² + pdy² - r₁²) = 0 + * + * 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: + * + * t = (B ± ⎷(B² - A·C)) / A + * + * 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). + * + * Additional observations (useful for optimizations): + * A does not depend on p + * + * 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 + */ + pixman_image_t *image = iter->image; + int x = iter->x; + int y = iter->y; + int width = iter->width; + uint32_t *buffer = iter->buffer; + + gradient_t *gradient = (gradient_t *)image; + radial_gradient_t *radial = (radial_gradient_t *)image; + uint32_t *end = buffer + width; + pixman_gradient_walker_t walker; + 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, image->common.repeat); + + if (image->common.transform) + { + if (!pixman_transform_point_3d (image->common.transform, &v)) + return iter->buffer; + + unit.vector[0] = image->common.transform->matrix[0][0]; + unit.vector[1] = image->common.transform->matrix[1][0]; + unit.vector[2] = image->common.transform->matrix[2][0]; + } + else + { + unit.vector[0] = pixman_fixed_1; + unit.vector[1] = 0; + unit.vector[2] = 0; + } + + if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1) + { + /* + * Given: + * + * t = (B ± ⎷(B² - A·C)) / A + * + * where + * + * 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·ux, + * pdy = pdy₀ + n·uy, + * + * 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 + */ + 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], + -((pixman_fixed_48_16_t) radial->c1.radius), + v.vector[0], v.vector[1], radial->c1.radius); + dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0], + 2 * (pixman_fixed_48_16_t) 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++) + { + *buffer = radial_compute_color (radial->a, b, c, + radial->inva, + radial->delta.radius, + radial->mindr, + &walker, + image->common.repeat); + } + + 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++) + { + if (v.vector[2] != 0) + { + double pdx, pdy, invv2, b, c; + + invv2 = 1. * pixman_fixed_1 / v.vector[2]; + + pdx = v.vector[0] * invv2 - radial->c1.x; + /* / pixman_fixed_1 */ + + pdy = v.vector[1] * invv2 - radial->c1.y; + /* / pixman_fixed_1 */ + + b = fdot (pdx, pdy, radial->c1.radius, + radial->delta.x, radial->delta.y, + radial->delta.radius); + /* / pixman_fixed_1 / pixman_fixed_1 */ + + 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, + image->common.repeat); + } + else + { + *buffer = 0; + } + } + + ++buffer; + + v.vector[0] += unit.vector[0]; + v.vector[1] += unit.vector[1]; + v.vector[2] += unit.vector[2]; + } + } + + iter->y++; + return iter->buffer; +} + +static uint32_t * +radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask) +{ + uint32_t *buffer = radial_get_scanline_narrow (iter, NULL); + + pixman_expand ((uint64_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width); + + return buffer; +} + +void +_pixman_radial_gradient_iter_init (pixman_image_t *image, pixman_iter_t *iter) +{ + if (iter->flags & ITER_NARROW) + iter->get_scanline = radial_get_scanline_narrow; + else + iter->get_scanline = radial_get_scanline_wide; +} + +PIXMAN_EXPORT pixman_image_t * +pixman_image_create_radial_gradient (pixman_point_fixed_t * inner, + pixman_point_fixed_t * outer, + pixman_fixed_t inner_radius, + pixman_fixed_t outer_radius, + const pixman_gradient_stop_t *stops, + int n_stops) +{ + pixman_image_t *image; + radial_gradient_t *radial; + + image = _pixman_image_allocate (); + + if (!image) + return NULL; + + radial = &image->radial; + + if (!_pixman_init_gradient (&radial->common, stops, n_stops)) + { + free (image); + return NULL; + } + + image->type = RADIAL; + + radial->c1.x = inner->x; + radial->c1.y = inner->y; + radial->c1.radius = inner_radius; + radial->c2.x = outer->x; + radial->c2.y = outer->y; + radial->c2.radius = outer_radius; + + /* 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; + + return image; +} + |