pixman/pixman-radial-gradient.c

   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      *  - 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
  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 det;
  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     det = fdot (b, a, 0, b, -c, 0);
 120     if (det >= 0)
 121     {
 122         double sqrtdet, t0, t1;
 123
 124         sqrtdet = sqrt (det);
 125         t0 = (b + sqrtdet) * inva;
 126         t1 = (b - sqrtdet) * inva;
 127
 128         if (repeat == PIXMAN_REPEAT_NONE)
 129         {
 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);
 134         }
 135         else
 136         {
 137             if (t0 * dr > mindr)
 138                 return _pixman_gradient_walker_pixel (walker, t0);
 139             else if (t1 * dr > mindr)
 140                 return _pixman_gradient_walker_pixel (walker, t1);
 141         }
 142     }
 143
 144     return 0;
 145 }
 146
 147 static uint32_t *
 148 radial_get_scanline_narrow (pixman_iter_t *iter, const uint32_t *mask)
 149 {
 150     /*
 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).
 154      *
 155      * In the radial gradient problem we are given two circles (c₁,r₁) and
 156      * (c₂,r₂) that define the gradient itself.
 157      *
 158      * Mathematically the gradient can be defined as the family of circles
 159      *
 160      *     ((1-t)·c₁ + t·(c₂), (1-t)·r₁ + t·r₂)
 161      *
 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].
 168      *
 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.
 172      *
 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.
 177      *
 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:
 180      *
 181      *     length((1-t)·c₁ + t·(c₂) - p) = (1-t)·r₁ + t·r₂
 182      *
 183      * Let's rewrite it in a simpler way, by defining some auxiliary
 184      * variables:
 185      *
 186      *     cd = c₂ - c₁
 187      *     pd = p - c₁
 188      *     dr = r₂ - r₁
 189      *     lenght(t·cd - pd) = r₁ + t·dr
 190      *
 191      * which actually means
 192      *
 193      *     hypot(t·cdx - pdx, t·cdy - pdy) = r₁ + t·dr
 194      *
 195      * or
 196      *
 197      *     ⎷((t·cdx - pdx)² + (t·cdy - pdy)²) = r₁ + t·dr.
 198      *
 199      * If we impose (as stated earlier) that r₁ + t·dr >= 0, it becomes:
 200      *
 201      *     (t·cdx - pdx)² + (t·cdy - pdy)² = (r₁ + t·dr)²
 202      *
 203      * where we can actually expand the squares and solve for t:
 204      *
 205      *     t²cdx² - 2t·cdx·pdx + pdx² + t²cdy² - 2t·cdy·pdy + pdy² =
 206      *       = r₁² + 2·r₁·t·dr + t²·dr²
 207      *
 208      *     (cdx² + cdy² - dr²)t² - 2(cdx·pdx + cdy·pdy + r₁·dr)t +
 209      *         (pdx² + pdy² - r₁²) = 0
 210      *
 211      *     A = cdx² + cdy² - dr²
 212      *     B = pdx·cdx + pdy·cdy + r₁·dr
 213      *     C = pdx² + pdy² - r₁²
 214      *     At² - 2Bt + C = 0
 215      *
 216      * The solutions (unless the equation degenerates because of A = 0) are:
 217      *
 218      *     t = (B ± ⎷(B² - A·C)) / A
 219      *
 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
 222      * range).
 223      *
 224      * Additional observations (useful for optimizations):
 225      * A does not depend on p
 226      *
 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
 229      *       have opposite sign
 230      */
 231     pixman_image_t *image = iter->image;
 232     int x = iter->x;
 233     int y = iter->y;
 234     int width = iter->width;
 235     uint32_t *buffer = iter->buffer;
 236
 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;
 242
 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;
 247
 248     _pixman_gradient_walker_init (&walker, gradient, image->common.repeat);
 249
 250     if (image->common.transform)
 251     {
 252         if (!pixman_transform_point_3d (image->common.transform, &v))
 253             return iter->buffer;
 254
 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];
 258     }
 259     else
 260     {
 261         unit.vector[0] = pixman_fixed_1;
 262         unit.vector[1] = 0;
 263         unit.vector[2] = 0;
 264     }
 265
 266     if (unit.vector[2] == 0 && v.vector[2] == pixman_fixed_1)
 267     {
 268         /*
 269          * Given:
 270          *
 271          * t = (B ± ⎷(B² - A·C)) / A
 272          *
 273          * where
 274          *
 275          * A = cdx² + cdy² - dr²
 276          * B = pdx·cdx + pdy·cdy + r₁·dr
 277          * C = pdx² + pdy² - r₁²
 278          * det = B² - A·C
 279          *
 280          * Since we have an affine transformation, we know that (pdx, pdy)
 281          * increase linearly with each pixel,
 282          *
 283          * pdx = pdx₀ + n·ux,
 284          * pdy = pdy₀ + n·uy,
 285          *
 286          * we can then express B, C and det through multiple differentiation.
 287          */
 288         pixman_fixed_32_32_t b, db, c, dc, ddc;
 289
 290         /* warning: this computation may overflow */
 291         v.vector[0] -= radial->c1.x;
 292         v.vector[1] -= radial->c1.y;
 293
 294         /*
 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
 302          */
 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);
 307
 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],
 313                   0,
 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);
 317
 318         while (buffer < end)
 319         {
 320             if (!mask || *mask++)
 321             {
 322                 *buffer = radial_compute_color (radial->a, b, c,
 323                                                 radial->inva,
 324                                                 radial->delta.radius,
 325                                                 radial->mindr,
 326                                                 &walker,
 327                                                 image->common.repeat);
 328             }
 329
 330             b += db;
 331             c += dc;
 332             dc += ddc;
 333             ++buffer;
 334         }
 335     }
 336     else
 337     {
 338         /* projective */
 339         /* Warning:
 340          * error propagation guarantees are much looser than in the affine case
 341          */
 342         while (buffer < end)
 343         {
 344             if (!mask || *mask++)
 345             {
 346                 if (v.vector[2] != 0)
 347                 {
 348                     double pdx, pdy, invv2, b, c;
 349
 350                     invv2 = 1. * pixman_fixed_1 / v.vector[2];
 351
 352                     pdx = v.vector[0] * invv2 - radial->c1.x;
 353                     /*    / pixman_fixed_1 */
 354
 355                     pdy = v.vector[1] * invv2 - radial->c1.y;
 356                     /*    / pixman_fixed_1 */
 357
 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 */
 362
 363                     c = fdot (pdx, pdy, -radial->c1.radius,
 364                               pdx, pdy, radial->c1.radius);
 365                     /*  / pixman_fixed_1 / pixman_fixed_1 */
 366
 367                     *buffer = radial_compute_color (radial->a, b, c,
 368                                                     radial->inva,
 369                                                     radial->delta.radius,
 370                                                     radial->mindr,
 371                                                     &walker,
 372                                                     image->common.repeat);
 373                 }
 374                 else
 375                 {
 376                     *buffer = 0;
 377                 }
 378             }
 379
 380             ++buffer;
 381
 382             v.vector[0] += unit.vector[0];
 383             v.vector[1] += unit.vector[1];
 384             v.vector[2] += unit.vector[2];
 385         }
 386     }
 387
 388     iter->y++;
 389     return iter->buffer;
 390 }
 391
 392 static uint32_t *
 393 radial_get_scanline_wide (pixman_iter_t *iter, const uint32_t *mask)
 394 {
 395     uint32_t *buffer = radial_get_scanline_narrow (iter, NULL);
 396
 397     pixman_expand ((uint64_t *)buffer, buffer, PIXMAN_a8r8g8b8, iter->width);
 398
 399     return buffer;
 400 }
 401
 402 void
 403 _pixman_radial_gradient_iter_init (pixman_image_t *image,
 404                                    pixman_iter_t *iter,
 405                                    int x, int y, int width, int height,
 406                                    uint8_t *buffer, iter_flags_t flags)
 407 {
 408     if (flags & ITER_NARROW)
 409         iter->get_scanline = radial_get_scanline_narrow;
 410     else
 411         iter->get_scanline = radial_get_scanline_wide;
 412 }
 413
 414 PIXMAN_EXPORT pixman_image_t *
 415 pixman_image_create_radial_gradient (pixman_point_fixed_t *        inner,
 416                                      pixman_point_fixed_t *        outer,
 417                                      pixman_fixed_t                inner_radius,
 418                                      pixman_fixed_t                outer_radius,
 419                                      const pixman_gradient_stop_t *stops,
 420                                      int                           n_stops)
 421 {
 422     pixman_image_t *image;
 423     radial_gradient_t *radial;
 424
 425     image = _pixman_image_allocate ();
 426
 427     if (!image)
 428         return NULL;
 429
 430     radial = &image->radial;
 431
 432     if (!_pixman_init_gradient (&radial->common, stops, n_stops))
 433     {
 434         free (image);
 435         return NULL;
 436     }
 437
 438     image->type = RADIAL;
 439
 440     radial->c1.x = inner->x;
 441     radial->c1.y = inner->y;
 442     radial->c1.radius = inner_radius;
 443     radial->c2.x = outer->x;
 444     radial->c2.y = outer->y;
 445     radial->c2.radius = outer_radius;
 446
 447     /* warning: this computations may overflow */
 448     radial->delta.x = radial->c2.x - radial->c1.x;
 449     radial->delta.y = radial->c2.y - radial->c1.y;
 450     radial->delta.radius = radial->c2.radius - radial->c1.radius;
 451
 452     /* computed exactly, then cast to double -> every bit of the double
 453        representation is correct (53 bits) */
 454     radial->a = dot (radial->delta.x, radial->delta.y, -radial->delta.radius,
 455                      radial->delta.x, radial->delta.y, radial->delta.radius);
 456     if (radial->a != 0)
 457         radial->inva = 1. * pixman_fixed_1 / radial->a;
 458
 459     radial->mindr = -1. * pixman_fixed_1 * radial->c1.radius;
 460
 461     return image;
 462 }
 463