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 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 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 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 * 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 + tr) 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 + tr 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(tcd - pd) = r + tdr 201 * 202 * which actually means 203 * 204 * hypot(tcdx - pdx, tcdy - pdy) = r + tdr 205 * 206 * or 207 * 208 * ((tcdx - pdx) + (tcdy - pdy)) = r + tdr. 209 * 210 * If we impose (as stated earlier) that r + tdr >= 0, it becomes: 211 * 212 * (tcdx - pdx) + (tcdy - pdy) = (r + tdr) 213 * 214 * where we can actually expand the squares and solve for t: 215 * 216 * tcdx - 2tcdxpdx + pdx + tcdy - 2tcdypdy + pdy = 217 * = r + 2rtdr + tdr 218 * 219 * (cdx + cdy - dr)t - 2(cdxpdx + cdypdy + rdr)t + 220 * (pdx + pdy - r) = 0 221 * 222 * A = cdx + cdy - dr 223 * B = pdxcdx + pdycdy + rdr 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 - AC)) / 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 - AC)) / A 283 * 284 * where 285 * 286 * A = cdx + cdy - dr 287 * B = pdxcdx + pdycdy + rdr 288 * C = pdx + pdy - r 289 * det = B - AC 290 * 291 * Since we have an affine transformation, we know that (pdx, pdy) 292 * increase linearly with each pixel, 293 * 294 * pdx = pdx + nux, 295 * pdy = pdy + nuy, 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 * 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 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 * 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
