1 /* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */
4 * Copyright © 2000 Keith Packard, member of The XFree86 Project, Inc.
5 * Copyright © 2000 SuSE, Inc.
6 * 2005 Lars Knoll & Zack Rusin, Trolltech
7 * Copyright © 2007 Red Hat, Inc.
10 * Permission to use, copy, modify, distribute, and sell this software and its
11 * documentation for any purpose is hereby granted without fee, provided that
12 * the above copyright notice appear in all copies and that both that
13 * copyright notice and this permission notice appear in supporting
14 * documentation, and that the name of Keith Packard not be used in
15 * advertising or publicity pertaining to distribution of the software without
16 * specific, written prior permission. Keith Packard makes no
17 * representations about the suitability of this software for any purpose. It
18 * is provided "as is" without express or implied warranty.
20 * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS
21 * SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
22 * FITNESS, IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
23 * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN
25 * AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING
26 * OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS
35 #include "pixman-private.h"
37 static inline pixman_fixed_32_32_t
38 dot (pixman_fixed_48_16_t x1,
39 pixman_fixed_48_16_t y1,
40 pixman_fixed_48_16_t z1,
41 pixman_fixed_48_16_t x2,
42 pixman_fixed_48_16_t y2,
43 pixman_fixed_48_16_t z2)
46 * Exact computation, assuming that the input values can
47 * be represented as pixman_fixed_16_16_t
49 return x1 * x2 + y1 * y2 + z1 * z2;
61 * Error can be unbound in some special cases.
62 * Using clever dot product algorithms (for example compensated
63 * dot product) would improve this but make the code much less
66 return x1 * x2 + y1 * y2 + z1 * z2;
70 radial_compute_color (double a,
76 pixman_gradient_walker_t *walker,
77 pixman_repeat_t repeat)
80 * In this function error propagation can lead to bad results:
81 * - det can have an unbound error (if b*b-a*c is very small),
82 * potentially making it the opposite sign of what it should have been
83 * (thus clearing a pixel that would have been colored or vice-versa)
84 * or propagating the error to sqrtdet;
85 * if det has the wrong sign or b is very small, this can lead to bad
88 * - the algorithm used to compute the solutions of the quadratic
89 * equation is not numerically stable (but saves one division compared
90 * to the numerically stable one);
91 * this can be a problem if a*c is much smaller than b*b
93 * - the above problems are worse if a is small (as inva becomes bigger)
104 t = pixman_fixed_1 / 2 * c / b;
105 if (repeat == PIXMAN_REPEAT_NONE)
107 if (0 <= t && t <= pixman_fixed_1)
108 return _pixman_gradient_walker_pixel (walker, t);
113 return _pixman_gradient_walker_pixel (walker, t);
119 det = fdot (b, a, 0, b, -c, 0);
122 double sqrtdet, t0, t1;
124 sqrtdet = sqrt (det);
125 t0 = (b + sqrtdet) * inva;
126 t1 = (b - sqrtdet) * inva;
128 if (repeat == PIXMAN_REPEAT_NONE)
130 if (0 <= t0 && t0 <= pixman_fixed_1)
131 return _pixman_gradient_walker_pixel (walker, t0);
132 else if (0 <= t1 && t1 <= pixman_fixed_1)
133 return _pixman_gradient_walker_pixel (walker, t1);
138 return _pixman_gradient_walker_pixel (walker, t0);
139 else if (t1 * dr > mindr)
140 return _pixman_gradient_walker_pixel (walker, t1);
148 radial_get_scanline_narrow (pixman_iter_t *iter, const uint32_t *mask)
151 * Implementation of radial gradients following the PDF specification.
152 * See section 8.7.4.5.4 Type 3 (Radial) Shadings of the PDF Reference
153 * Manual (PDF 32000-1:2008 at the time of this writing).
155 * In the radial gradient problem we are given two circles (c₁,r₁) and
156 * (c₂,r₂) that define the gradient itself.
158 * Mathematically the gradient can be defined as the family of circles
160 * ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
162 * excluding those circles whose radius would be < 0. When a point
163 * belongs to more than one circle, the one with a bigger t is the only
164 * one that contributes to its color. When a point does not belong
165 * to any of the circles, it is transparent black, i.e. RGBA (0, 0, 0, 0).
166 * Further limitations on the range of values for t are imposed when
167 * the gradient is not repeated, namely t must belong to [0,1].
169 * The graphical result is the same as drawing the valid (radius > 0)
170 * circles with increasing t in [-inf, +inf] (or in [0,1] if the gradient
171 * is not repeated) using SOURCE operatior composition.
173 * It looks like a cone pointing towards the viewer if the ending circle
174 * is smaller than the starting one, a cone pointing inside the page if
175 * the starting circle is the smaller one and like a cylinder if they
176 * have the same radius.
178 * What we actually do is, given the point whose color we are interested
179 * in, compute the t values for that point, solving for t in:
181 * length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂
183 * Let's rewrite it in a simpler way, by defining some auxiliary
189 * lenght(t·cd - pd) = r₁ + t·dr
191 * which actually means
193 * hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
197 * ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
199 * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
201 * (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
203 * where we can actually expand the squares and solve for t:
205 * t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
206 * = r₁² + 2·r₁·t·dr + t²·dr²
208 * (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
209 * (pdx² + pdy² - r₁²) = 0
211 * A = cdx² + cdy² - dr²
212 * B = pdx·cdx + pdy·cdy + r₁·dr
213 * C = pdx² + pdy² - r₁²
216 * The solutions (unless the equation degenerates because of A = 0) are:
218 * t = (B ± ⎷(B² - A·C)) / A
220 * The solution we are going to prefer is the bigger one, unless the
221 * radius associated to it is negative (or it falls outside the valid t
224 * Additional observations (useful for optimizations):
225 * A does not depend on p
227 * A < 0 <=> one of the two circles completely contains the other one
228 * <=> for every p, the radiuses associated with the two t solutions
231 pixman_image_t *image = iter->image;
234 int width = iter->width;
235 uint32_t *buffer = iter->buffer;
237 gradient_t *gradient = (gradient_t *)image;
238 radial_gradient_t *radial = (radial_gradient_t *)image;
239 uint32_t *end = buffer + width;
240 pixman_gradient_walker_t walker;
241 pixman_vector_t v, unit;
243 /* reference point is the center of the pixel */
244 v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2;
245 v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2;
246 v.vector[2] = pixman_fixed_1;
248 _pixman_gradient_walker_init (&walker, gradient, image->common.repeat);
250 if (image->common.transform)
252 if (!pixman_transform_point_3d (image->common.transform, &v))
255 unit.vector[0] = image->common.transform->matrix[0][0];
256 unit.vector[1] = image->common.transform->matrix[1][0];
257 unit.vector[2] = image->common.transform->matrix[2][0];
261 unit.vector[0] = pixman_fixed_1;
266 if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
271 * t = (B ± ⎷(B² - A·C)) / A
275 * A = cdx² + cdy² - dr²
276 * B = pdx·cdx + pdy·cdy + r₁·dr
277 * C = pdx² + pdy² - r₁²
280 * Since we have an affine transformation, we know that (pdx, pdy)
281 * increase linearly with each pixel,
286 * we can then express B, C and det through multiple differentiation.
288 pixman_fixed_32_32_t b, db, c, dc, ddc;
290 /* warning: this computation may overflow */
291 v.vector[0] -= radial->c1.x;
292 v.vector[1] -= radial->c1.y;
295 * B and C are computed and updated exactly.
296 * If fdot was used instead of dot, in the worst case it would
297 * lose 11 bits of precision in each of the multiplication and
298 * summing up would zero out all the bit that were preserved,
299 * thus making the result 0 instead of the correct one.
300 * This would mean a worst case of unbound relative error or
301 * about 2^10 absolute error
303 b = dot (v.vector[0], v.vector[1], radial->c1.radius,
304 radial->delta.x, radial->delta.y, radial->delta.radius);
305 db = dot (unit.vector[0], unit.vector[1], 0,
306 radial->delta.x, radial->delta.y, 0);
308 c = dot (v.vector[0], v.vector[1],
309 -((pixman_fixed_48_16_t) radial->c1.radius),
310 v.vector[0], v.vector[1], radial->c1.radius);
311 dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0],
312 2 * (pixman_fixed_48_16_t) v.vector[1] + unit.vector[1],
314 unit.vector[0], unit.vector[1], 0);
315 ddc = 2 * dot (unit.vector[0], unit.vector[1], 0,
316 unit.vector[0], unit.vector[1], 0);
320 if (!mask || *mask++)
322 *buffer = radial_compute_color (radial->a, b, c,
324 radial->delta.radius,
327 image->common.repeat);
340 * error propagation guarantees are much looser than in the affine case
344 if (!mask || *mask++)
346 if (v.vector[2] != 0)
348 double pdx, pdy, invv2, b, c;
350 invv2 = 1. * pixman_fixed_1 / v.vector[2];
352 pdx = v.vector[0] * invv2 - radial->c1.x;
353 /* / pixman_fixed_1 */
355 pdy = v.vector[1] * invv2 - radial->c1.y;
356 /* / pixman_fixed_1 */
358 b = fdot (pdx, pdy, radial->c1.radius,
359 radial->delta.x, radial->delta.y,
360 radial->delta.radius);
361 /* / pixman_fixed_1 / pixman_fixed_1 */
363 c = fdot (pdx, pdy, -radial->c1.radius,
364 pdx, pdy, radial->c1.radius);
365 /* / pixman_fixed_1 / pixman_fixed_1 */
367 *buffer = radial_compute_color (radial->a, b, c,
369 radial->delta.radius,
372 image->common.repeat);
382 v.vector[0] += unit.vector[0];
383 v.vector[1] += unit.vector[1];
384 v.vector[2] += unit.vector[2];
393 radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask)
395 uint32_t *buffer = radial_get_scanline_narrow (iter, NULL);
397 pixman_expand ((uint64_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width);
403 _pixman_radial_gradient_iter_init (pixman_image_t *image, pixman_iter_t *iter)
405 if (iter->flags & ITER_NARROW)
406 iter->get_scanline = radial_get_scanline_narrow;
408 iter->get_scanline = radial_get_scanline_wide;
411 PIXMAN_EXPORT pixman_image_t *
412 pixman_image_create_radial_gradient (pixman_point_fixed_t * inner,
413 pixman_point_fixed_t * outer,
414 pixman_fixed_t inner_radius,
415 pixman_fixed_t outer_radius,
416 const pixman_gradient_stop_t *stops,
419 pixman_image_t *image;
420 radial_gradient_t *radial;
422 image = _pixman_image_allocate ();
427 radial = &image->radial;
429 if (!_pixman_init_gradient (&radial->common, stops, n_stops))
435 image->type = RADIAL;
437 radial->c1.x = inner->x;
438 radial->c1.y = inner->y;
439 radial->c1.radius = inner_radius;
440 radial->c2.x = outer->x;
441 radial->c2.y = outer->y;
442 radial->c2.radius = outer_radius;
444 /* warning: this computations may overflow */
445 radial->delta.x = radial->c2.x - radial->c1.x;
446 radial->delta.y = radial->c2.y - radial->c1.y;
447 radial->delta.radius = radial->c2.radius - radial->c1.radius;
449 /* computed exactly, then cast to double -> every bit of the double
450 representation is correct (53 bits) */
451 radial->a = dot (radial->delta.x, radial->delta.y, -radial->delta.radius,
452 radial->delta.x, radial->delta.y, radial->delta.radius);
454 radial->inva = 1. * pixman_fixed_1 / radial->a;
456 radial->mindr = -1. * pixman_fixed_1 * radial->c1.radius;