1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/cairo/libpixman/src/pixman-radial-gradient.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,727 @@
1.4 +/* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */
1.5 +/*
1.6 + *
1.7 + * Copyright © 2000 Keith Packard, member of The XFree86 Project, Inc.
1.8 + * Copyright © 2000 SuSE, Inc.
1.9 + * 2005 Lars Knoll & Zack Rusin, Trolltech
1.10 + * Copyright © 2007 Red Hat, Inc.
1.11 + *
1.12 + *
1.13 + * Permission to use, copy, modify, distribute, and sell this software and its
1.14 + * documentation for any purpose is hereby granted without fee, provided that
1.15 + * the above copyright notice appear in all copies and that both that
1.16 + * copyright notice and this permission notice appear in supporting
1.17 + * documentation, and that the name of Keith Packard not be used in
1.18 + * advertising or publicity pertaining to distribution of the software without
1.19 + * specific, written prior permission. Keith Packard makes no
1.20 + * representations about the suitability of this software for any purpose. It
1.21 + * is provided "as is" without express or implied warranty.
1.22 + *
1.23 + * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS
1.24 + * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
1.25 + * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
1.26 + * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
1.27 + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN
1.28 + * AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING
1.29 + * OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS
1.30 + * SOFTWARE.
1.31 + */
1.32 +
1.33 +#ifdef HAVE_CONFIG_H
1.34 +#include <config.h>
1.35 +#endif
1.36 +#include <stdlib.h>
1.37 +#include <math.h>
1.38 +#include "pixman-private.h"
1.39 +
1.40 +#include "pixman-dither.h"
1.41 +
1.42 +static inline pixman_fixed_32_32_t
1.43 +dot (pixman_fixed_48_16_t x1,
1.44 + pixman_fixed_48_16_t y1,
1.45 + pixman_fixed_48_16_t z1,
1.46 + pixman_fixed_48_16_t x2,
1.47 + pixman_fixed_48_16_t y2,
1.48 + pixman_fixed_48_16_t z2)
1.49 +{
1.50 + /*
1.51 + * Exact computation, assuming that the input values can
1.52 + * be represented as pixman_fixed_16_16_t
1.53 + */
1.54 + return x1 * x2 + y1 * y2 + z1 * z2;
1.55 +}
1.56 +
1.57 +static inline double
1.58 +fdot (double x1,
1.59 + double y1,
1.60 + double z1,
1.61 + double x2,
1.62 + double y2,
1.63 + double z2)
1.64 +{
1.65 + /*
1.66 + * Error can be unbound in some special cases.
1.67 + * Using clever dot product algorithms (for example compensated
1.68 + * dot product) would improve this but make the code much less
1.69 + * obvious
1.70 + */
1.71 + return x1 * x2 + y1 * y2 + z1 * z2;
1.72 +}
1.73 +
1.74 +static uint32_t
1.75 +radial_compute_color (double a,
1.76 + double b,
1.77 + double c,
1.78 + double inva,
1.79 + double dr,
1.80 + double mindr,
1.81 + pixman_gradient_walker_t *walker,
1.82 + pixman_repeat_t repeat)
1.83 +{
1.84 + /*
1.85 + * In this function error propagation can lead to bad results:
1.86 + * - discr can have an unbound error (if b*b-a*c is very small),
1.87 + * potentially making it the opposite sign of what it should have been
1.88 + * (thus clearing a pixel that would have been colored or vice-versa)
1.89 + * or propagating the error to sqrtdiscr;
1.90 + * if discr has the wrong sign or b is very small, this can lead to bad
1.91 + * results
1.92 + *
1.93 + * - the algorithm used to compute the solutions of the quadratic
1.94 + * equation is not numerically stable (but saves one division compared
1.95 + * to the numerically stable one);
1.96 + * this can be a problem if a*c is much smaller than b*b
1.97 + *
1.98 + * - the above problems are worse if a is small (as inva becomes bigger)
1.99 + */
1.100 + double discr;
1.101 +
1.102 + if (a == 0)
1.103 + {
1.104 + double t;
1.105 +
1.106 + if (b == 0)
1.107 + return 0;
1.108 +
1.109 + t = pixman_fixed_1 / 2 * c / b;
1.110 + if (repeat == PIXMAN_REPEAT_NONE)
1.111 + {
1.112 + if (0 <= t && t <= pixman_fixed_1)
1.113 + return _pixman_gradient_walker_pixel (walker, t);
1.114 + }
1.115 + else
1.116 + {
1.117 + if (t * dr >= mindr)
1.118 + return _pixman_gradient_walker_pixel (walker, t);
1.119 + }
1.120 +
1.121 + return 0;
1.122 + }
1.123 +
1.124 + discr = fdot (b, a, 0, b, -c, 0);
1.125 + if (discr >= 0)
1.126 + {
1.127 + double sqrtdiscr, t0, t1;
1.128 +
1.129 + sqrtdiscr = sqrt (discr);
1.130 + t0 = (b + sqrtdiscr) * inva;
1.131 + t1 = (b - sqrtdiscr) * inva;
1.132 +
1.133 + /*
1.134 + * The root that must be used is the biggest one that belongs
1.135 + * to the valid range ([0,1] for PIXMAN_REPEAT_NONE, any
1.136 + * solution that results in a positive radius otherwise).
1.137 + *
1.138 + * If a > 0, t0 is the biggest solution, so if it is valid, it
1.139 + * is the correct result.
1.140 + *
1.141 + * If a < 0, only one of the solutions can be valid, so the
1.142 + * order in which they are tested is not important.
1.143 + */
1.144 + if (repeat == PIXMAN_REPEAT_NONE)
1.145 + {
1.146 + if (0 <= t0 && t0 <= pixman_fixed_1)
1.147 + return _pixman_gradient_walker_pixel (walker, t0);
1.148 + else if (0 <= t1 && t1 <= pixman_fixed_1)
1.149 + return _pixman_gradient_walker_pixel (walker, t1);
1.150 + }
1.151 + else
1.152 + {
1.153 + if (t0 * dr >= mindr)
1.154 + return _pixman_gradient_walker_pixel (walker, t0);
1.155 + else if (t1 * dr >= mindr)
1.156 + return _pixman_gradient_walker_pixel (walker, t1);
1.157 + }
1.158 + }
1.159 +
1.160 + return 0;
1.161 +}
1.162 +
1.163 +static uint32_t *
1.164 +radial_get_scanline_narrow (pixman_iter_t *iter, const uint32_t *mask)
1.165 +{
1.166 + /*
1.167 + * Implementation of radial gradients following the PDF specification.
1.168 + * See section 8.7.4.5.4 Type 3 (Radial) Shadings of the PDF Reference
1.169 + * Manual (PDF 32000-1:2008 at the time of this writing).
1.170 + *
1.171 + * In the radial gradient problem we are given two circles (c₁,r₁) and
1.172 + * (c₂,r₂) that define the gradient itself.
1.173 + *
1.174 + * Mathematically the gradient can be defined as the family of circles
1.175 + *
1.176 + * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
1.177 + *
1.178 + * excluding those circles whose radius would be < 0. When a point
1.179 + * belongs to more than one circle, the one with a bigger t is the only
1.180 + * one that contributes to its color. When a point does not belong
1.181 + * to any of the circles, it is transparent black, i.e. RGBA (0, 0, 0, 0).
1.182 + * Further limitations on the range of values for t are imposed when
1.183 + * the gradient is not repeated, namely t must belong to [0,1].
1.184 + *
1.185 + * The graphical result is the same as drawing the valid (radius > 0)
1.186 + * circles with increasing t in [-inf, +inf] (or in [0,1] if the gradient
1.187 + * is not repeated) using SOURCE operator composition.
1.188 + *
1.189 + * It looks like a cone pointing towards the viewer if the ending circle
1.190 + * is smaller than the starting one, a cone pointing inside the page if
1.191 + * the starting circle is the smaller one and like a cylinder if they
1.192 + * have the same radius.
1.193 + *
1.194 + * What we actually do is, given the point whose color we are interested
1.195 + * in, compute the t values for that point, solving for t in:
1.196 + *
1.197 + * length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂
1.198 + *
1.199 + * Let's rewrite it in a simpler way, by defining some auxiliary
1.200 + * variables:
1.201 + *
1.202 + * cd = c₂ - c₁
1.203 + * pd = p - c₁
1.204 + * dr = r₂ - r₁
1.205 + * length(t·cd - pd) = r₁ + t·dr
1.206 + *
1.207 + * which actually means
1.208 + *
1.209 + * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
1.210 + *
1.211 + * or
1.212 + *
1.213 + * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
1.214 + *
1.215 + * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
1.216 + *
1.217 + * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
1.218 + *
1.219 + * where we can actually expand the squares and solve for t:
1.220 + *
1.221 + * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
1.222 + * = r₁² + 2·r₁·t·dr + t²·dr²
1.223 + *
1.224 + * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
1.225 + * (pdx² + pdy² - r₁²) = 0
1.226 + *
1.227 + * A = cdx² + cdy² - dr²
1.228 + * B = pdx·cdx + pdy·cdy + r₁·dr
1.229 + * C = pdx² + pdy² - r₁²
1.230 + * At² - 2Bt + C = 0
1.231 + *
1.232 + * The solutions (unless the equation degenerates because of A = 0) are:
1.233 + *
1.234 + * t = (B ± ⎷(B² - A·C)) / A
1.235 + *
1.236 + * The solution we are going to prefer is the bigger one, unless the
1.237 + * radius associated to it is negative (or it falls outside the valid t
1.238 + * range).
1.239 + *
1.240 + * Additional observations (useful for optimizations):
1.241 + * A does not depend on p
1.242 + *
1.243 + * A < 0 <=> one of the two circles completely contains the other one
1.244 + * <=> for every p, the radiuses associated with the two t solutions
1.245 + * have opposite sign
1.246 + */
1.247 + pixman_image_t *image = iter->image;
1.248 + int x = iter->x;
1.249 + int y = iter->y;
1.250 + int width = iter->width;
1.251 + uint32_t *buffer = iter->buffer;
1.252 +
1.253 + gradient_t *gradient = (gradient_t *)image;
1.254 + radial_gradient_t *radial = (radial_gradient_t *)image;
1.255 + uint32_t *end = buffer + width;
1.256 + pixman_gradient_walker_t walker;
1.257 + pixman_vector_t v, unit;
1.258 +
1.259 + /* reference point is the center of the pixel */
1.260 + v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2;
1.261 + v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2;
1.262 + v.vector[2] = pixman_fixed_1;
1.263 +
1.264 + _pixman_gradient_walker_init (&walker, gradient, image->common.repeat);
1.265 +
1.266 + if (image->common.transform)
1.267 + {
1.268 + if (!pixman_transform_point_3d (image->common.transform, &v))
1.269 + return iter->buffer;
1.270 +
1.271 + unit.vector[0] = image->common.transform->matrix[0][0];
1.272 + unit.vector[1] = image->common.transform->matrix[1][0];
1.273 + unit.vector[2] = image->common.transform->matrix[2][0];
1.274 + }
1.275 + else
1.276 + {
1.277 + unit.vector[0] = pixman_fixed_1;
1.278 + unit.vector[1] = 0;
1.279 + unit.vector[2] = 0;
1.280 + }
1.281 +
1.282 + if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
1.283 + {
1.284 + /*
1.285 + * Given:
1.286 + *
1.287 + * t = (B ± ⎷(B² - A·C)) / A
1.288 + *
1.289 + * where
1.290 + *
1.291 + * A = cdx² + cdy² - dr²
1.292 + * B = pdx·cdx + pdy·cdy + r₁·dr
1.293 + * C = pdx² + pdy² - r₁²
1.294 + * det = B² - A·C
1.295 + *
1.296 + * Since we have an affine transformation, we know that (pdx, pdy)
1.297 + * increase linearly with each pixel,
1.298 + *
1.299 + * pdx = pdx₀ + n·ux,
1.300 + * pdy = pdy₀ + n·uy,
1.301 + *
1.302 + * we can then express B, C and det through multiple differentiation.
1.303 + */
1.304 + pixman_fixed_32_32_t b, db, c, dc, ddc;
1.305 +
1.306 + /* warning: this computation may overflow */
1.307 + v.vector[0] -= radial->c1.x;
1.308 + v.vector[1] -= radial->c1.y;
1.309 +
1.310 + /*
1.311 + * B and C are computed and updated exactly.
1.312 + * If fdot was used instead of dot, in the worst case it would
1.313 + * lose 11 bits of precision in each of the multiplication and
1.314 + * summing up would zero out all the bit that were preserved,
1.315 + * thus making the result 0 instead of the correct one.
1.316 + * This would mean a worst case of unbound relative error or
1.317 + * about 2^10 absolute error
1.318 + */
1.319 + b = dot (v.vector[0], v.vector[1], radial->c1.radius,
1.320 + radial->delta.x, radial->delta.y, radial->delta.radius);
1.321 + db = dot (unit.vector[0], unit.vector[1], 0,
1.322 + radial->delta.x, radial->delta.y, 0);
1.323 +
1.324 + c = dot (v.vector[0], v.vector[1],
1.325 + -((pixman_fixed_48_16_t) radial->c1.radius),
1.326 + v.vector[0], v.vector[1], radial->c1.radius);
1.327 + dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0],
1.328 + 2 * (pixman_fixed_48_16_t) v.vector[1] + unit.vector[1],
1.329 + 0,
1.330 + unit.vector[0], unit.vector[1], 0);
1.331 + ddc = 2 * dot (unit.vector[0], unit.vector[1], 0,
1.332 + unit.vector[0], unit.vector[1], 0);
1.333 +
1.334 + while (buffer < end)
1.335 + {
1.336 + if (!mask || *mask++)
1.337 + {
1.338 + *buffer = radial_compute_color (radial->a, b, c,
1.339 + radial->inva,
1.340 + radial->delta.radius,
1.341 + radial->mindr,
1.342 + &walker,
1.343 + image->common.repeat);
1.344 + }
1.345 +
1.346 + b += db;
1.347 + c += dc;
1.348 + dc += ddc;
1.349 + ++buffer;
1.350 + }
1.351 + }
1.352 + else
1.353 + {
1.354 + /* projective */
1.355 + /* Warning:
1.356 + * error propagation guarantees are much looser than in the affine case
1.357 + */
1.358 + while (buffer < end)
1.359 + {
1.360 + if (!mask || *mask++)
1.361 + {
1.362 + if (v.vector[2] != 0)
1.363 + {
1.364 + double pdx, pdy, invv2, b, c;
1.365 +
1.366 + invv2 = 1. * pixman_fixed_1 / v.vector[2];
1.367 +
1.368 + pdx = v.vector[0] * invv2 - radial->c1.x;
1.369 + /* / pixman_fixed_1 */
1.370 +
1.371 + pdy = v.vector[1] * invv2 - radial->c1.y;
1.372 + /* / pixman_fixed_1 */
1.373 +
1.374 + b = fdot (pdx, pdy, radial->c1.radius,
1.375 + radial->delta.x, radial->delta.y,
1.376 + radial->delta.radius);
1.377 + /* / pixman_fixed_1 / pixman_fixed_1 */
1.378 +
1.379 + c = fdot (pdx, pdy, -radial->c1.radius,
1.380 + pdx, pdy, radial->c1.radius);
1.381 + /* / pixman_fixed_1 / pixman_fixed_1 */
1.382 +
1.383 + *buffer = radial_compute_color (radial->a, b, c,
1.384 + radial->inva,
1.385 + radial->delta.radius,
1.386 + radial->mindr,
1.387 + &walker,
1.388 + image->common.repeat);
1.389 + }
1.390 + else
1.391 + {
1.392 + *buffer = 0;
1.393 + }
1.394 + }
1.395 +
1.396 + ++buffer;
1.397 +
1.398 + v.vector[0] += unit.vector[0];
1.399 + v.vector[1] += unit.vector[1];
1.400 + v.vector[2] += unit.vector[2];
1.401 + }
1.402 + }
1.403 +
1.404 + iter->y++;
1.405 + return iter->buffer;
1.406 +}
1.407 +
1.408 +static uint32_t *
1.409 +radial_get_scanline_16 (pixman_iter_t *iter, const uint32_t *mask)
1.410 +{
1.411 + /*
1.412 + * Implementation of radial gradients following the PDF specification.
1.413 + * See section 8.7.4.5.4 Type 3 (Radial) Shadings of the PDF Reference
1.414 + * Manual (PDF 32000-1:2008 at the time of this writing).
1.415 + *
1.416 + * In the radial gradient problem we are given two circles (c₁,r₁) and
1.417 + * (c₂,r₂) that define the gradient itself.
1.418 + *
1.419 + * Mathematically the gradient can be defined as the family of circles
1.420 + *
1.421 + * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
1.422 + *
1.423 + * excluding those circles whose radius would be < 0. When a point
1.424 + * belongs to more than one circle, the one with a bigger t is the only
1.425 + * one that contributes to its color. When a point does not belong
1.426 + * to any of the circles, it is transparent black, i.e. RGBA (0, 0, 0, 0).
1.427 + * Further limitations on the range of values for t are imposed when
1.428 + * the gradient is not repeated, namely t must belong to [0,1].
1.429 + *
1.430 + * The graphical result is the same as drawing the valid (radius > 0)
1.431 + * circles with increasing t in [-inf, +inf] (or in [0,1] if the gradient
1.432 + * is not repeated) using SOURCE operator composition.
1.433 + *
1.434 + * It looks like a cone pointing towards the viewer if the ending circle
1.435 + * is smaller than the starting one, a cone pointing inside the page if
1.436 + * the starting circle is the smaller one and like a cylinder if they
1.437 + * have the same radius.
1.438 + *
1.439 + * What we actually do is, given the point whose color we are interested
1.440 + * in, compute the t values for that point, solving for t in:
1.441 + *
1.442 + * length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂
1.443 + *
1.444 + * Let's rewrite it in a simpler way, by defining some auxiliary
1.445 + * variables:
1.446 + *
1.447 + * cd = c₂ - c₁
1.448 + * pd = p - c₁
1.449 + * dr = r₂ - r₁
1.450 + * length(t·cd - pd) = r₁ + t·dr
1.451 + *
1.452 + * which actually means
1.453 + *
1.454 + * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
1.455 + *
1.456 + * or
1.457 + *
1.458 + * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
1.459 + *
1.460 + * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
1.461 + *
1.462 + * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
1.463 + *
1.464 + * where we can actually expand the squares and solve for t:
1.465 + *
1.466 + * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
1.467 + * = r₁² + 2·r₁·t·dr + t²·dr²
1.468 + *
1.469 + * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
1.470 + * (pdx² + pdy² - r₁²) = 0
1.471 + *
1.472 + * A = cdx² + cdy² - dr²
1.473 + * B = pdx·cdx + pdy·cdy + r₁·dr
1.474 + * C = pdx² + pdy² - r₁²
1.475 + * At² - 2Bt + C = 0
1.476 + *
1.477 + * The solutions (unless the equation degenerates because of A = 0) are:
1.478 + *
1.479 + * t = (B ± ⎷(B² - A·C)) / A
1.480 + *
1.481 + * The solution we are going to prefer is the bigger one, unless the
1.482 + * radius associated to it is negative (or it falls outside the valid t
1.483 + * range).
1.484 + *
1.485 + * Additional observations (useful for optimizations):
1.486 + * A does not depend on p
1.487 + *
1.488 + * A < 0 <=> one of the two circles completely contains the other one
1.489 + * <=> for every p, the radiuses associated with the two t solutions
1.490 + * have opposite sign
1.491 + */
1.492 + pixman_image_t *image = iter->image;
1.493 + int x = iter->x;
1.494 + int y = iter->y;
1.495 + int width = iter->width;
1.496 + uint16_t *buffer = iter->buffer;
1.497 + pixman_bool_t toggle = ((x ^ y) & 1);
1.498 +
1.499 + gradient_t *gradient = (gradient_t *)image;
1.500 + radial_gradient_t *radial = (radial_gradient_t *)image;
1.501 + uint16_t *end = buffer + width;
1.502 + pixman_gradient_walker_t walker;
1.503 + pixman_vector_t v, unit;
1.504 +
1.505 + /* reference point is the center of the pixel */
1.506 + v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2;
1.507 + v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2;
1.508 + v.vector[2] = pixman_fixed_1;
1.509 +
1.510 + _pixman_gradient_walker_init (&walker, gradient, image->common.repeat);
1.511 +
1.512 + if (image->common.transform)
1.513 + {
1.514 + if (!pixman_transform_point_3d (image->common.transform, &v))
1.515 + return iter->buffer;
1.516 +
1.517 + unit.vector[0] = image->common.transform->matrix[0][0];
1.518 + unit.vector[1] = image->common.transform->matrix[1][0];
1.519 + unit.vector[2] = image->common.transform->matrix[2][0];
1.520 + }
1.521 + else
1.522 + {
1.523 + unit.vector[0] = pixman_fixed_1;
1.524 + unit.vector[1] = 0;
1.525 + unit.vector[2] = 0;
1.526 + }
1.527 +
1.528 + if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
1.529 + {
1.530 + /*
1.531 + * Given:
1.532 + *
1.533 + * t = (B ± ⎷(B² - A·C)) / A
1.534 + *
1.535 + * where
1.536 + *
1.537 + * A = cdx² + cdy² - dr²
1.538 + * B = pdx·cdx + pdy·cdy + r₁·dr
1.539 + * C = pdx² + pdy² - r₁²
1.540 + * det = B² - A·C
1.541 + *
1.542 + * Since we have an affine transformation, we know that (pdx, pdy)
1.543 + * increase linearly with each pixel,
1.544 + *
1.545 + * pdx = pdx₀ + n·ux,
1.546 + * pdy = pdy₀ + n·uy,
1.547 + *
1.548 + * we can then express B, C and det through multiple differentiation.
1.549 + */
1.550 + pixman_fixed_32_32_t b, db, c, dc, ddc;
1.551 +
1.552 + /* warning: this computation may overflow */
1.553 + v.vector[0] -= radial->c1.x;
1.554 + v.vector[1] -= radial->c1.y;
1.555 +
1.556 + /*
1.557 + * B and C are computed and updated exactly.
1.558 + * If fdot was used instead of dot, in the worst case it would
1.559 + * lose 11 bits of precision in each of the multiplication and
1.560 + * summing up would zero out all the bit that were preserved,
1.561 + * thus making the result 0 instead of the correct one.
1.562 + * This would mean a worst case of unbound relative error or
1.563 + * about 2^10 absolute error
1.564 + */
1.565 + b = dot (v.vector[0], v.vector[1], radial->c1.radius,
1.566 + radial->delta.x, radial->delta.y, radial->delta.radius);
1.567 + db = dot (unit.vector[0], unit.vector[1], 0,
1.568 + radial->delta.x, radial->delta.y, 0);
1.569 +
1.570 + c = dot (v.vector[0], v.vector[1],
1.571 + -((pixman_fixed_48_16_t) radial->c1.radius),
1.572 + v.vector[0], v.vector[1], radial->c1.radius);
1.573 + dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0],
1.574 + 2 * (pixman_fixed_48_16_t) v.vector[1] + unit.vector[1],
1.575 + 0,
1.576 + unit.vector[0], unit.vector[1], 0);
1.577 + ddc = 2 * dot (unit.vector[0], unit.vector[1], 0,
1.578 + unit.vector[0], unit.vector[1], 0);
1.579 +
1.580 + while (buffer < end)
1.581 + {
1.582 + if (!mask || *mask++)
1.583 + {
1.584 + *buffer = dither_8888_to_0565(
1.585 + radial_compute_color (radial->a, b, c,
1.586 + radial->inva,
1.587 + radial->delta.radius,
1.588 + radial->mindr,
1.589 + &walker,
1.590 + image->common.repeat),
1.591 + toggle);
1.592 + }
1.593 +
1.594 + toggle ^= 1;
1.595 + b += db;
1.596 + c += dc;
1.597 + dc += ddc;
1.598 + ++buffer;
1.599 + }
1.600 + }
1.601 + else
1.602 + {
1.603 + /* projective */
1.604 + /* Warning:
1.605 + * error propagation guarantees are much looser than in the affine case
1.606 + */
1.607 + while (buffer < end)
1.608 + {
1.609 + if (!mask || *mask++)
1.610 + {
1.611 + if (v.vector[2] != 0)
1.612 + {
1.613 + double pdx, pdy, invv2, b, c;
1.614 +
1.615 + invv2 = 1. * pixman_fixed_1 / v.vector[2];
1.616 +
1.617 + pdx = v.vector[0] * invv2 - radial->c1.x;
1.618 + /* / pixman_fixed_1 */
1.619 +
1.620 + pdy = v.vector[1] * invv2 - radial->c1.y;
1.621 + /* / pixman_fixed_1 */
1.622 +
1.623 + b = fdot (pdx, pdy, radial->c1.radius,
1.624 + radial->delta.x, radial->delta.y,
1.625 + radial->delta.radius);
1.626 + /* / pixman_fixed_1 / pixman_fixed_1 */
1.627 +
1.628 + c = fdot (pdx, pdy, -radial->c1.radius,
1.629 + pdx, pdy, radial->c1.radius);
1.630 + /* / pixman_fixed_1 / pixman_fixed_1 */
1.631 +
1.632 + *buffer = dither_8888_to_0565 (
1.633 + radial_compute_color (radial->a, b, c,
1.634 + radial->inva,
1.635 + radial->delta.radius,
1.636 + radial->mindr,
1.637 + &walker,
1.638 + image->common.repeat),
1.639 + toggle);
1.640 + }
1.641 + else
1.642 + {
1.643 + *buffer = 0;
1.644 + }
1.645 + }
1.646 +
1.647 + ++buffer;
1.648 + toggle ^= 1;
1.649 +
1.650 + v.vector[0] += unit.vector[0];
1.651 + v.vector[1] += unit.vector[1];
1.652 + v.vector[2] += unit.vector[2];
1.653 + }
1.654 + }
1.655 +
1.656 + iter->y++;
1.657 + return iter->buffer;
1.658 +}
1.659 +static uint32_t *
1.660 +radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask)
1.661 +{
1.662 + uint32_t *buffer = radial_get_scanline_narrow (iter, NULL);
1.663 +
1.664 + pixman_expand_to_float (
1.665 + (argb_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width);
1.666 +
1.667 + return buffer;
1.668 +}
1.669 +
1.670 +void
1.671 +_pixman_radial_gradient_iter_init (pixman_image_t *image, pixman_iter_t *iter)
1.672 +{
1.673 + if (iter->iter_flags & ITER_16)
1.674 + iter->get_scanline = radial_get_scanline_16;
1.675 + else if (iter->iter_flags & ITER_NARROW)
1.676 + iter->get_scanline = radial_get_scanline_narrow;
1.677 + else
1.678 + iter->get_scanline = radial_get_scanline_wide;
1.679 +}
1.680 +
1.681 +
1.682 +PIXMAN_EXPORT pixman_image_t *
1.683 +pixman_image_create_radial_gradient (const pixman_point_fixed_t * inner,
1.684 + const pixman_point_fixed_t * outer,
1.685 + pixman_fixed_t inner_radius,
1.686 + pixman_fixed_t outer_radius,
1.687 + const pixman_gradient_stop_t *stops,
1.688 + int n_stops)
1.689 +{
1.690 + pixman_image_t *image;
1.691 + radial_gradient_t *radial;
1.692 +
1.693 + image = _pixman_image_allocate ();
1.694 +
1.695 + if (!image)
1.696 + return NULL;
1.697 +
1.698 + radial = &image->radial;
1.699 +
1.700 + if (!_pixman_init_gradient (&radial->common, stops, n_stops))
1.701 + {
1.702 + free (image);
1.703 + return NULL;
1.704 + }
1.705 +
1.706 + image->type = RADIAL;
1.707 +
1.708 + radial->c1.x = inner->x;
1.709 + radial->c1.y = inner->y;
1.710 + radial->c1.radius = inner_radius;
1.711 + radial->c2.x = outer->x;
1.712 + radial->c2.y = outer->y;
1.713 + radial->c2.radius = outer_radius;
1.714 +
1.715 + /* warning: this computations may overflow */
1.716 + radial->delta.x = radial->c2.x - radial->c1.x;
1.717 + radial->delta.y = radial->c2.y - radial->c1.y;
1.718 + radial->delta.radius = radial->c2.radius - radial->c1.radius;
1.719 +
1.720 + /* computed exactly, then cast to double -> every bit of the double
1.721 + representation is correct (53 bits) */
1.722 + radial->a = dot (radial->delta.x, radial->delta.y, -radial->delta.radius,
1.723 + radial->delta.x, radial->delta.y, radial->delta.radius);
1.724 + if (radial->a != 0)
1.725 + radial->inva = 1. * pixman_fixed_1 / radial->a;
1.726 +
1.727 + radial->mindr = -1. * pixman_fixed_1 * radial->c1.radius;
1.728 +
1.729 + return image;
1.730 +}
