pixman/pixman-filter.c

   1 /*
   2  * Copyright 2012, Red Hat, Inc.
   3  * Copyright 2012, Soren Sandmann
   4  *
   5  * Permission is hereby granted, free of charge, to any person obtaining a
   6  * copy of this software and associated documentation files (the "Software"),
   7  * to deal in the Software without restriction, including without limitation
   8  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9  * and/or sell copies of the Software, and to permit persons to whom the
  10  * Software is furnished to do so, subject to the following conditions:
  11  *
  12  * The above copyright notice and this permission notice (including the next
  13  * paragraph) shall be included in all copies or substantial portions of the
  14  * Software.
  15  *
  16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  19  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  21  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  22  * DEALINGS IN THE SOFTWARE.
  23  *
  24  * Author: Soren Sandmann <soren.sandmann@gmail.com>
  25  */
  26 #include <string.h>
  27 #include <stdlib.h>
  28 #include <stdio.h>
  29 #include <math.h>
  30 #include <assert.h>
  31 #ifdef HAVE_CONFIG_H
  32 #include <config.h>
  33 #endif
  34 #include "pixman-private.h"
  35
  36 typedef double (* kernel_func_t) (double x);
  37
  38 typedef struct
  39 {
  40     pixman_kernel_t     kernel;
  41     kernel_func_t       func;
  42     double              width;
  43 } filter_info_t;
  44
  45 static double
  46 impulse_kernel (double x)
  47 {
  48     return (x == 0.0)? 1.0 : 0.0;
  49 }
  50
  51 static double
  52 box_kernel (double x)
  53 {
  54     return 1;
  55 }
  56
  57 static double
  58 linear_kernel (double x)
  59 {
  60     return 1 - fabs (x);
  61 }
  62
  63 static double
  64 gaussian_kernel (double x)
  65 {
  66 #define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
  67 #define SIGMA (SQRT2 / 2.0)
  68
  69     return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
  70 }
  71
  72 static double
  73 sinc (double x)
  74 {
  75     if (x == 0.0)
  76         return 1.0;
  77     else
  78         return sin (M_PI * x) / (M_PI * x);
  79 }
  80
  81 static double
  82 lanczos (double x, int n)
  83 {
  84     return sinc (x) * sinc (x * (1.0 / n));
  85 }
  86
  87 static double
  88 lanczos2_kernel (double x)
  89 {
  90     return lanczos (x, 2);
  91 }
  92
  93 static double
  94 lanczos3_kernel (double x)
  95 {
  96     return lanczos (x, 3);
  97 }
  98
  99 static double
 100 nice_kernel (double x)
 101 {
 102     return lanczos3_kernel (x * 0.75);
 103 }
 104
 105 static double
 106 general_cubic (double x, double B, double C)
 107 {
 108     double ax = fabs(x);
 109
 110     if (ax < 1)
 111     {
 112         return (((12 - 9 * B - 6 * C) * ax +
 113                  (-18 + 12 * B + 6 * C)) * ax * ax +
 114                 (6 - 2 * B)) / 6;
 115     }
 116     else if (ax < 2)
 117     {
 118         return ((((-B - 6 * C) * ax +
 119                   (6 * B + 30 * C)) * ax +
 120                  (-12 * B - 48 * C)) * ax +
 121                 (8 * B + 24 * C)) / 6;
 122     }
 123     else
 124     {
 125         return 0;
 126     }
 127 }
 128
 129 static double
 130 cubic_kernel (double x)
 131 {
 132     /* This is the Mitchell-Netravali filter.
 133      *
 134      * (0.0, 0.5) would give us the Catmull-Rom spline,
 135      * but that one seems to be indistinguishable from Lanczos2.
 136      */
 137     return general_cubic (x, 1/3.0, 1/3.0);
 138 }
 139
 140 static const filter_info_t filters[] =
 141 {
 142     { PIXMAN_KERNEL_IMPULSE,            impulse_kernel,   0.0 },
 143     { PIXMAN_KERNEL_BOX,                box_kernel,       1.0 },
 144     { PIXMAN_KERNEL_LINEAR,             linear_kernel,    2.0 },
 145     { PIXMAN_KERNEL_CUBIC,              cubic_kernel,     4.0 },
 146     { PIXMAN_KERNEL_GAUSSIAN,           gaussian_kernel,  5.0 },
 147     { PIXMAN_KERNEL_LANCZOS2,           lanczos2_kernel,  4.0 },
 148     { PIXMAN_KERNEL_LANCZOS3,           lanczos3_kernel,  6.0 },
 149     { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
 150 };
 151
 152 /* This function scales @kernel2 by @scale, then
 153  * aligns @x1 in @kernel1 with @x2 in @kernel2 and
 154  * and integrates the product of the kernels across @width.
 155  *
 156  * This function assumes that the intervals are within
 157  * the kernels in question. E.g., the caller must not
 158  * try to integrate a linear kernel ouside of [-1:1]
 159  */
 160 static double
 161 integral (pixman_kernel_t kernel1, double x1,
 162           pixman_kernel_t kernel2, double scale, double x2,
 163           double width)
 164 {
 165     if (kernel1 == PIXMAN_KERNEL_BOX && kernel2 == PIXMAN_KERNEL_BOX)
 166     {
 167         return width;
 168     }
 169     /* The LINEAR filter is not differentiable at 0, so if the
 170      * integration interval crosses zero, break it into two
 171      * separate integrals.
 172      */
 173     else if (kernel1 == PIXMAN_KERNEL_LINEAR && x1 < 0 && x1 + width > 0)
 174     {
 175         return
 176             integral (kernel1, x1, kernel2, scale, x2, - x1) +
 177             integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
 178     }
 179     else if (kernel2 == PIXMAN_KERNEL_LINEAR && x2 < 0 && x2 + width > 0)
 180     {
 181         return
 182             integral (kernel1, x1, kernel2, scale, x2, - x2) +
 183             integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
 184     }
 185     else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
 186     {
 187         assert (width == 0.0);
 188         return filters[kernel2].func (x2 * scale);
 189     }
 190     else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
 191     {
 192         assert (width == 0.0);
 193         return filters[kernel1].func (x1);
 194     }
 195     else
 196     {
 197         /* Integration via Simpson's rule
 198          * See http://www.intmath.com/integration/6-simpsons-rule.php
 199          * 12 segments (6 cubic approximations) seems to produce best
 200          * result for lanczos3.linear, which was the combination that
 201          * showed the most errors.  This makes sense as the lanczos3
 202          * filter is 6 wide.
 203          */
 204 #define N_SEGMENTS 12
 205 #define SAMPLE(a1, a2)                                                  \
 206         (filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
 207
 208         double s = 0.0;
 209         double h = width / N_SEGMENTS;
 210         int i;
 211
 212         s = SAMPLE (x1, x2);
 213
 214         for (i = 1; i < N_SEGMENTS; i += 2)
 215         {
 216             double a1 = x1 + h * i;
 217             double a2 = x2 + h * i;
 218             s += 4 * SAMPLE (a1, a2);
 219         }
 220
 221         for (i = 2; i < N_SEGMENTS; i += 2)
 222         {
 223             double a1 = x1 + h * i;
 224             double a2 = x2 + h * i;
 225             s += 2 * SAMPLE (a1, a2);
 226         }
 227
 228         s += SAMPLE (x1 + width, x2 + width);
 229
 230         return h * s * (1.0 / 3.0);
 231     }
 232 }
 233
 234 static void
 235 create_1d_filter (int              width,
 236                   pixman_kernel_t  reconstruct,
 237                   pixman_kernel_t  sample,
 238                   double           scale,
 239                   int              n_phases,
 240                   pixman_fixed_t *p)
 241 {
 242     double step;
 243     int i;
 244
 245     step = 1.0 / n_phases;
 246
 247     for (i = 0; i < n_phases; ++i)
 248     {
 249         double frac = step / 2.0 + i * step;
 250         pixman_fixed_t new_total;
 251         int x, x1, x2;
 252         double total, e;
 253
 254         /* Sample convolution of reconstruction and sampling
 255          * filter. See rounding.txt regarding the rounding
 256          * and sample positions.
 257          */
 258
 259         x1 = ceil (frac - width / 2.0 - 0.5);
 260         x2 = x1 + width;
 261
 262         total = 0;
 263         for (x = x1; x < x2; ++x)
 264         {
 265             double pos = x + 0.5 - frac;
 266             double rlow = - filters[reconstruct].width / 2.0;
 267             double rhigh = rlow + filters[reconstruct].width;
 268             double slow = pos - scale * filters[sample].width / 2.0;
 269             double shigh = slow + scale * filters[sample].width;
 270             double c = 0.0;
 271             double ilow, ihigh;
 272
 273             if (rhigh >= slow && rlow <= shigh)
 274             {
 275                 ilow = MAX (slow, rlow);
 276                 ihigh = MIN (shigh, rhigh);
 277
 278                 c = integral (reconstruct, ilow,
 279                               sample, 1.0 / scale, ilow - pos,
 280                               ihigh - ilow);
 281             }
 282
 283             *p = (pixman_fixed_t)floor (c * 65536.0 + 0.5);
 284             total += *p;
 285             p++;
 286         }
 287
 288         /* Normalize, with error diffusion */
 289         p -= width;
 290         total = 65536.0 / total;
 291         new_total = 0;
 292         e = 0.0;
 293         for (x = x1; x < x2; ++x)
 294         {
 295             double v = (*p) * total + e;
 296             pixman_fixed_t t = floor (v + 0.5);
 297
 298             e = v - t;
 299             new_total += t;
 300             *p++ = t;
 301         }
 302
 303         /* pixman_fixed_e's worth of error may remain; put it
 304          * at the first sample, since that is the only one that
 305          * hasn't had any error diffused into it.
 306          */
 307         *(p - width) += pixman_fixed_1 - new_total;
 308     }
 309 }
 310
 311
 312 static int
 313 filter_width (pixman_kernel_t reconstruct, pixman_kernel_t sample, double size)
 314 {
 315     return ceil (filters[reconstruct].width + size * filters[sample].width);
 316 }
 317
 318 #ifdef PIXMAN_GNUPLOT
 319
 320 /* If enable-gnuplot is configured, then you can pipe the output of a
 321  * pixman-using program to gnuplot and get a continuously-updated plot
 322  * of the horizontal filter. This works well with demos/scale to test
 323  * the filter generation.
 324  *
 325  * The plot is all the different subposition filters shuffled
 326  * together. This is misleading in a few cases:
 327  *
 328  *  IMPULSE.BOX - goes up and down as the subfilters have different
 329  *                numbers of non-zero samples
 330  *  IMPULSE.TRIANGLE - somewhat crooked for the same reason
 331  *  1-wide filters - looks triangular, but a 1-wide box would be more
 332  *                   accurate
 333  */
 334 static void
 335 gnuplot_filter (int width, int n_phases, const pixman_fixed_t* p)
 336 {
 337     double step;
 338     int i, j;
 339     int first;
 340
 341     step = 1.0 / n_phases;
 342
 343     printf ("set style line 1 lc rgb '#0060ad' lt 1 lw 0.5 pt 7 pi 1 ps 0.5\n");
 344     printf ("plot [x=%g:%g] '-' with linespoints ls 1\n", -width*0.5, width*0.5);
 345     /* Print a point at the origin so that y==0 line is included: */
 346     printf ("0 0\n\n");
 347
 348     /* The position of the first sample of the phase corresponding to
 349      * frac is given by:
 350      *
 351      *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
 352      *
 353      * We have to find the frac that minimizes this expression.
 354      *
 355      * For odd widths, we have
 356      *
 357      *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
 358      *   = ceil (frac) + K - frac
 359      *   = 1 + K - frac
 360      *
 361      * for some K, so this is minimized when frac is maximized and
 362      * strictly growing with frac. So for odd widths, we can simply
 363      * start at the last phase and go backwards.
 364      *
 365      * For even widths, we have
 366      *
 367      *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
 368      *   = ceil (frac - 0.5) + K - frac
 369      *
 370      * The graph for this function (ignoring K) looks like this:
 371      *
 372      *        0.5
 373      *           |    |\
 374      *           |    | \
 375      *           |    |  \
 376      *         0 |    |   \
 377      *           |\   |
 378      *           | \  |
 379      *           |  \ |
 380      *      -0.5 |   \|
 381      *   ---------------------------------
 382      *           0    0.5   1
 383      *
 384      * So in this case we need to start with the phase whose frac is
 385      * less than, but as close as possible to 0.5, then go backwards
 386      * until we hit the first phase, then wrap around to the last
 387      * phase and continue backwards.
 388      *
 389      * Which phase is as close as possible 0.5? The locations of the
 390      * sampling point corresponding to the kth phase is given by
 391      * 1/(2 * n_phases) + k / n_phases:
 392      *
 393      *         1/(2 * n_phases) + k / n_phases = 0.5
 394      *
 395      * from which it follows that
 396      *
 397      *         k = (n_phases - 1) / 2
 398      *
 399      * rounded down is the phase in question.
 400      */
 401     if (width & 1)
 402         first = n_phases - 1;
 403     else
 404         first = (n_phases - 1) / 2;
 405
 406     for (j = 0; j < width; ++j)
 407     {
 408         for (i = 0; i < n_phases; ++i)
 409         {
 410             int phase = first - i;
 411             double frac, pos;
 412
 413             if (phase < 0)
 414                 phase = n_phases + phase;
 415
 416             frac = step / 2.0 + phase * step;
 417             pos = ceil (frac - width / 2.0 - 0.5) + 0.5 - frac + j;
 418
 419             printf ("%g %g\n",
 420                     pos,
 421                     pixman_fixed_to_double (*(p + phase * width + j)));
 422         }
 423     }
 424
 425     printf ("e\n");
 426     fflush (stdout);
 427 }
 428
 429 #endif
 430
 431 /* Create the parameter list for a SEPARABLE_CONVOLUTION filter
 432  * with the given kernels and scale parameters
 433  */
 434 PIXMAN_EXPORT pixman_fixed_t *
 435 pixman_filter_create_separable_convolution (int             *n_values,
 436                                             pixman_fixed_t   scale_x,
 437                                             pixman_fixed_t   scale_y,
 438                                             pixman_kernel_t  reconstruct_x,
 439                                             pixman_kernel_t  reconstruct_y,
 440                                             pixman_kernel_t  sample_x,
 441                                             pixman_kernel_t  sample_y,
 442                                             int              subsample_bits_x,
 443                                             int              subsample_bits_y)
 444 {
 445     double sx = fabs (pixman_fixed_to_double (scale_x));
 446     double sy = fabs (pixman_fixed_to_double (scale_y));
 447     pixman_fixed_t *params;
 448     int subsample_x, subsample_y;
 449     int width, height;
 450
 451     width = filter_width (reconstruct_x, sample_x, sx);
 452     subsample_x = (1 << subsample_bits_x);
 453
 454     height = filter_width (reconstruct_y, sample_y, sy);
 455     subsample_y = (1 << subsample_bits_y);
 456
 457     *n_values = 4 + width * subsample_x + height * subsample_y;
 458
 459     params = malloc (*n_values * sizeof (pixman_fixed_t));
 460     if (!params)
 461         return NULL;
 462
 463     params[0] = pixman_int_to_fixed (width);
 464     params[1] = pixman_int_to_fixed (height);
 465     params[2] = pixman_int_to_fixed (subsample_bits_x);
 466     params[3] = pixman_int_to_fixed (subsample_bits_y);
 467
 468     create_1d_filter (width, reconstruct_x, sample_x, sx, subsample_x,
 469                       params + 4);
 470     create_1d_filter (height, reconstruct_y, sample_y, sy, subsample_y,
 471                       params + 4 + width * subsample_x);
 472
 473 #ifdef PIXMAN_GNUPLOT
 474     gnuplot_filter(width, subsample_x, params + 4);
 475 #endif
 476
 477     return params;
 478 }