Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  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. }
  472.