• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* -*- Mode: c; c-basic-offset: 4; tab-width: 8; indent-tabs-mode: t; -*- */
2 /*
3  *
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.
8  *
9  *
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.
19  *
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
27  * SOFTWARE.
28  */
29 
30 #ifdef HAVE_CONFIG_H
31 #include <config.h>
32 #endif
33 #include <stdlib.h>
34 #include <math.h>
35 #include "pixman-private.h"
36 
37 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)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)
44 {
45     /*
46      * Exact computation, assuming that the input values can
47      * be represented as pixman_fixed_16_16_t
48      */
49     return x1 * x2 + y1 * y2 + z1 * z2;
50 }
51 
52 static inline double
fdot(double x1,double y1,double z1,double x2,double y2,double z2)53 fdot (double x1,
54       double y1,
55       double z1,
56       double x2,
57       double y2,
58       double z2)
59 {
60     /*
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
64      * obvious
65      */
66     return x1 * x2 + y1 * y2 + z1 * z2;
67 }
68 
69 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)70 radial_compute_color (double                    a,
71 		      double                    b,
72 		      double                    c,
73 		      double                    inva,
74 		      double                    dr,
75 		      double                    mindr,
76 		      pixman_gradient_walker_t *walker,
77 		      pixman_repeat_t           repeat)
78 {
79     /*
80      * In this function error propagation can lead to bad results:
81      *  - discr 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 sqrtdiscr;
85      *    if discr has the wrong sign or b is very small, this can lead to bad
86      *    results
87      *
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
92      *
93      *  - the above problems are worse if a is small (as inva becomes bigger)
94      */
95     double discr;
96 
97     if (a == 0)
98     {
99 	double t;
100 
101 	if (b == 0)
102 	    return 0;
103 
104 	t = pixman_fixed_1 / 2 * c / b;
105 	if (repeat == PIXMAN_REPEAT_NONE)
106 	{
107 	    if (0 <= t && t <= pixman_fixed_1)
108 		return _pixman_gradient_walker_pixel (walker, t);
109 	}
110 	else
111 	{
112 	    if (t * dr >= mindr)
113 		return _pixman_gradient_walker_pixel (walker, t);
114 	}
115 
116 	return 0;
117     }
118 
119     discr = fdot (b, a, 0, b, -c, 0);
120     if (discr >= 0)
121     {
122 	double sqrtdiscr, t0, t1;
123 
124 	sqrtdiscr = sqrt (discr);
125 	t0 = (b + sqrtdiscr) * inva;
126 	t1 = (b - sqrtdiscr) * inva;
127 
128 	/*
129 	 * The root that must be used is the biggest one that belongs
130 	 * to the valid range ([0,1] for PIXMAN_REPEAT_NONE, any
131 	 * solution that results in a positive radius otherwise).
132 	 *
133 	 * If a > 0, t0 is the biggest solution, so if it is valid, it
134 	 * is the correct result.
135 	 *
136 	 * If a < 0, only one of the solutions can be valid, so the
137 	 * order in which they are tested is not important.
138 	 */
139 	if (repeat == PIXMAN_REPEAT_NONE)
140 	{
141 	    if (0 <= t0 && t0 <= pixman_fixed_1)
142 		return _pixman_gradient_walker_pixel (walker, t0);
143 	    else if (0 <= t1 && t1 <= pixman_fixed_1)
144 		return _pixman_gradient_walker_pixel (walker, t1);
145 	}
146 	else
147 	{
148 	    if (t0 * dr >= mindr)
149 		return _pixman_gradient_walker_pixel (walker, t0);
150 	    else if (t1 * dr >= mindr)
151 		return _pixman_gradient_walker_pixel (walker, t1);
152 	}
153     }
154 
155     return 0;
156 }
157 
158 static uint32_t *
radial_get_scanline_narrow(pixman_iter_t * iter,const uint32_t * mask)159 radial_get_scanline_narrow (pixman_iter_t *iter, const uint32_t *mask)
160 {
161     /*
162      * Implementation of radial gradients following the PDF specification.
163      * See section 8.7.4.5.4 Type 3 (Radial) Shadings of the PDF Reference
164      * Manual (PDF 32000-1:2008 at the time of this writing).
165      *
166      * In the radial gradient problem we are given two circles (c₁,r₁) and
167      * (c₂,r₂) that define the gradient itself.
168      *
169      * Mathematically the gradient can be defined as the family of circles
170      *
171      *     ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
172      *
173      * excluding those circles whose radius would be < 0. When a point
174      * belongs to more than one circle, the one with a bigger t is the only
175      * one that contributes to its color. When a point does not belong
176      * to any of the circles, it is transparent black, i.e. RGBA (0, 0, 0, 0).
177      * Further limitations on the range of values for t are imposed when
178      * the gradient is not repeated, namely t must belong to [0,1].
179      *
180      * The graphical result is the same as drawing the valid (radius > 0)
181      * circles with increasing t in [-inf, +inf] (or in [0,1] if the gradient
182      * is not repeated) using SOURCE operator composition.
183      *
184      * It looks like a cone pointing towards the viewer if the ending circle
185      * is smaller than the starting one, a cone pointing inside the page if
186      * the starting circle is the smaller one and like a cylinder if they
187      * have the same radius.
188      *
189      * What we actually do is, given the point whose color we are interested
190      * in, compute the t values for that point, solving for t in:
191      *
192      *     length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂
193      *
194      * Let's rewrite it in a simpler way, by defining some auxiliary
195      * variables:
196      *
197      *     cd = c₂ - c₁
198      *     pd = p - c₁
199      *     dr = r₂ - r₁
200      *     length(t·cd - pd) = r₁ + t·dr
201      *
202      * which actually means
203      *
204      *     hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
205      *
206      * or
207      *
208      *     ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
209      *
210      * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
211      *
212      *     (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
213      *
214      * where we can actually expand the squares and solve for t:
215      *
216      *     t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
217      *       = r₁² + 2·r₁·t·dr + t²·dr²
218      *
219      *     (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
220      *         (pdx² + pdy² - r₁²) = 0
221      *
222      *     A = cdx² + cdy² - dr²
223      *     B = pdx·cdx + pdy·cdy + r₁·dr
224      *     C = pdx² + pdy² - r₁²
225      *     At² - 2Bt + C = 0
226      *
227      * The solutions (unless the equation degenerates because of A = 0) are:
228      *
229      *     t = (B ± ⎷(B² - A·C)) / A
230      *
231      * The solution we are going to prefer is the bigger one, unless the
232      * radius associated to it is negative (or it falls outside the valid t
233      * range).
234      *
235      * Additional observations (useful for optimizations):
236      * A does not depend on p
237      *
238      * A < 0 <=> one of the two circles completely contains the other one
239      *   <=> for every p, the radiuses associated with the two t solutions
240      *       have opposite sign
241      */
242     pixman_image_t *image = iter->image;
243     int x = iter->x;
244     int y = iter->y;
245     int width = iter->width;
246     uint32_t *buffer = iter->buffer;
247 
248     gradient_t *gradient = (gradient_t *)image;
249     radial_gradient_t *radial = (radial_gradient_t *)image;
250     uint32_t *end = buffer + width;
251     pixman_gradient_walker_t walker;
252     pixman_vector_t v, unit;
253 
254     /* reference point is the center of the pixel */
255     v.vector[0] = pixman_int_to_fixed (x) + pixman_fixed_1 / 2;
256     v.vector[1] = pixman_int_to_fixed (y) + pixman_fixed_1 / 2;
257     v.vector[2] = pixman_fixed_1;
258 
259     _pixman_gradient_walker_init (&walker, gradient, image->common.repeat);
260 
261     if (image->common.transform)
262     {
263 	if (!pixman_transform_point_3d (image->common.transform, &v))
264 	    return iter->buffer;
265 
266 	unit.vector[0] = image->common.transform->matrix[0][0];
267 	unit.vector[1] = image->common.transform->matrix[1][0];
268 	unit.vector[2] = image->common.transform->matrix[2][0];
269     }
270     else
271     {
272 	unit.vector[0] = pixman_fixed_1;
273 	unit.vector[1] = 0;
274 	unit.vector[2] = 0;
275     }
276 
277     if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
278     {
279 	/*
280 	 * Given:
281 	 *
282 	 * t = (B ± ⎷(B² - A·C)) / A
283 	 *
284 	 * where
285 	 *
286 	 * A = cdx² + cdy² - dr²
287 	 * B = pdx·cdx + pdy·cdy + r₁·dr
288 	 * C = pdx² + pdy² - r₁²
289 	 * det = B² - A·C
290 	 *
291 	 * Since we have an affine transformation, we know that (pdx, pdy)
292 	 * increase linearly with each pixel,
293 	 *
294 	 * pdx = pdx₀ + n·ux,
295 	 * pdy = pdy₀ + n·uy,
296 	 *
297 	 * we can then express B, C and det through multiple differentiation.
298 	 */
299 	pixman_fixed_32_32_t b, db, c, dc, ddc;
300 
301 	/* warning: this computation may overflow */
302 	v.vector[0] -= radial->c1.x;
303 	v.vector[1] -= radial->c1.y;
304 
305 	/*
306 	 * B and C are computed and updated exactly.
307 	 * If fdot was used instead of dot, in the worst case it would
308 	 * lose 11 bits of precision in each of the multiplication and
309 	 * summing up would zero out all the bit that were preserved,
310 	 * thus making the result 0 instead of the correct one.
311 	 * This would mean a worst case of unbound relative error or
312 	 * about 2^10 absolute error
313 	 */
314 	b = dot (v.vector[0], v.vector[1], radial->c1.radius,
315 		 radial->delta.x, radial->delta.y, radial->delta.radius);
316 	db = dot (unit.vector[0], unit.vector[1], 0,
317 		  radial->delta.x, radial->delta.y, 0);
318 
319 	c = dot (v.vector[0], v.vector[1],
320 		 -((pixman_fixed_48_16_t) radial->c1.radius),
321 		 v.vector[0], v.vector[1], radial->c1.radius);
322 	dc = dot (2 * (pixman_fixed_48_16_t) v.vector[0] + unit.vector[0],
323 		  2 * (pixman_fixed_48_16_t) v.vector[1] + unit.vector[1],
324 		  0,
325 		  unit.vector[0], unit.vector[1], 0);
326 	ddc = 2 * dot (unit.vector[0], unit.vector[1], 0,
327 		       unit.vector[0], unit.vector[1], 0);
328 
329 	while (buffer < end)
330 	{
331 	    if (!mask || *mask++)
332 	    {
333 		*buffer = radial_compute_color (radial->a, b, c,
334 						radial->inva,
335 						radial->delta.radius,
336 						radial->mindr,
337 						&walker,
338 						image->common.repeat);
339 	    }
340 
341 	    b += db;
342 	    c += dc;
343 	    dc += ddc;
344 	    ++buffer;
345 	}
346     }
347     else
348     {
349 	/* projective */
350 	/* Warning:
351 	 * error propagation guarantees are much looser than in the affine case
352 	 */
353 	while (buffer < end)
354 	{
355 	    if (!mask || *mask++)
356 	    {
357 		if (v.vector[2] != 0)
358 		{
359 		    double pdx, pdy, invv2, b, c;
360 
361 		    invv2 = 1. * pixman_fixed_1 / v.vector[2];
362 
363 		    pdx = v.vector[0] * invv2 - radial->c1.x;
364 		    /*    / pixman_fixed_1 */
365 
366 		    pdy = v.vector[1] * invv2 - radial->c1.y;
367 		    /*    / pixman_fixed_1 */
368 
369 		    b = fdot (pdx, pdy, radial->c1.radius,
370 			      radial->delta.x, radial->delta.y,
371 			      radial->delta.radius);
372 		    /*  / pixman_fixed_1 / pixman_fixed_1 */
373 
374 		    c = fdot (pdx, pdy, -radial->c1.radius,
375 			      pdx, pdy, radial->c1.radius);
376 		    /*  / pixman_fixed_1 / pixman_fixed_1 */
377 
378 		    *buffer = radial_compute_color (radial->a, b, c,
379 						    radial->inva,
380 						    radial->delta.radius,
381 						    radial->mindr,
382 						    &walker,
383 						    image->common.repeat);
384 		}
385 		else
386 		{
387 		    *buffer = 0;
388 		}
389 	    }
390 
391 	    ++buffer;
392 
393 	    v.vector[0] += unit.vector[0];
394 	    v.vector[1] += unit.vector[1];
395 	    v.vector[2] += unit.vector[2];
396 	}
397     }
398 
399     iter->y++;
400     return iter->buffer;
401 }
402 
403 static uint32_t *
radial_get_scanline_wide(pixman_iter_t * iter,const uint32_t * mask)404 radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask)
405 {
406     uint32_t *buffer = radial_get_scanline_narrow (iter, NULL);
407 
408     pixman_expand_to_float (
409 	(argb_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width);
410 
411     return buffer;
412 }
413 
414 void
_pixman_radial_gradient_iter_init(pixman_image_t * image,pixman_iter_t * iter)415 _pixman_radial_gradient_iter_init (pixman_image_t *image, pixman_iter_t *iter)
416 {
417     if (iter->iter_flags & ITER_NARROW)
418 	iter->get_scanline = radial_get_scanline_narrow;
419     else
420 	iter->get_scanline = radial_get_scanline_wide;
421 }
422 
423 PIXMAN_EXPORT pixman_image_t *
pixman_image_create_radial_gradient(const pixman_point_fixed_t * inner,const pixman_point_fixed_t * outer,pixman_fixed_t inner_radius,pixman_fixed_t outer_radius,const pixman_gradient_stop_t * stops,int n_stops)424 pixman_image_create_radial_gradient (const pixman_point_fixed_t *  inner,
425                                      const pixman_point_fixed_t *  outer,
426                                      pixman_fixed_t                inner_radius,
427                                      pixman_fixed_t                outer_radius,
428                                      const pixman_gradient_stop_t *stops,
429                                      int                           n_stops)
430 {
431     pixman_image_t *image;
432     radial_gradient_t *radial;
433 
434     image = _pixman_image_allocate ();
435 
436     if (!image)
437 	return NULL;
438 
439     radial = &image->radial;
440 
441     if (!_pixman_init_gradient (&radial->common, stops, n_stops))
442     {
443 	free (image);
444 	return NULL;
445     }
446 
447     image->type = RADIAL;
448 
449     radial->c1.x = inner->x;
450     radial->c1.y = inner->y;
451     radial->c1.radius = inner_radius;
452     radial->c2.x = outer->x;
453     radial->c2.y = outer->y;
454     radial->c2.radius = outer_radius;
455 
456     /* warning: this computations may overflow */
457     radial->delta.x = radial->c2.x - radial->c1.x;
458     radial->delta.y = radial->c2.y - radial->c1.y;
459     radial->delta.radius = radial->c2.radius - radial->c1.radius;
460 
461     /* computed exactly, then cast to double -> every bit of the double
462        representation is correct (53 bits) */
463     radial->a = dot (radial->delta.x, radial->delta.y, -radial->delta.radius,
464 		     radial->delta.x, radial->delta.y, radial->delta.radius);
465     if (radial->a != 0)
466 	radial->inva = 1. * pixman_fixed_1 / radial->a;
467 
468     radial->mindr = -1. * pixman_fixed_1 * radial->c1.radius;
469 
470     return image;
471 }
472