qemu/pixman/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      *  - 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₁ + 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 *
 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 }