qemu/pixman/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 * ax * ax +
 113                 (-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6;
 114     }
 115     else if (ax >= 1 && ax < 2)
 116     {
 117         return ((-B - 6 * C) * ax * ax * ax +
 118                 (6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) *
 119                 ax + (8 * B + 24 * C)) / 6;
 120     }
 121     else
 122     {
 123         return 0;
 124     }
 125 }
 126
 127 static double
 128 cubic_kernel (double x)
 129 {
 130     /* This is the Mitchell-Netravali filter.
 131      *
 132      * (0.0, 0.5) would give us the Catmull-Rom spline,
 133      * but that one seems to be indistinguishable from Lanczos2.
 134      */
 135     return general_cubic (x, 1/3.0, 1/3.0);
 136 }
 137
 138 static const filter_info_t filters[] =
 139 {
 140     { PIXMAN_KERNEL_IMPULSE,            impulse_kernel,   0.0 },
 141     { PIXMAN_KERNEL_BOX,                box_kernel,       1.0 },
 142     { PIXMAN_KERNEL_LINEAR,             linear_kernel,    2.0 },
 143     { PIXMAN_KERNEL_CUBIC,              cubic_kernel,     4.0 },
 144     { PIXMAN_KERNEL_GAUSSIAN,           gaussian_kernel,  6 * SIGMA },
 145     { PIXMAN_KERNEL_LANCZOS2,           lanczos2_kernel,  4.0 },
 146     { PIXMAN_KERNEL_LANCZOS3,           lanczos3_kernel,  6.0 },
 147     { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
 148 };
 149
 150 /* This function scales @kernel2 by @scale, then
 151  * aligns @x1 in @kernel1 with @x2 in @kernel2 and
 152  * and integrates the product of the kernels across @width.
 153  *
 154  * This function assumes that the intervals are within
 155  * the kernels in question. E.g., the caller must not
 156  * try to integrate a linear kernel ouside of [-1:1]
 157  */
 158 static double
 159 integral (pixman_kernel_t kernel1, double x1,
 160           pixman_kernel_t kernel2, double scale, double x2,
 161           double width)
 162 {
 163     /* If the integration interval crosses zero, break it into
 164      * two separate integrals. This ensures that filters such
 165      * as LINEAR that are not differentiable at 0 will still
 166      * integrate properly.
 167      */
 168     if (x1 < 0 && x1 + width > 0)
 169     {
 170         return
 171             integral (kernel1, x1, kernel2, scale, x2, - x1) +
 172             integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
 173     }
 174     else if (x2 < 0 && x2 + width > 0)
 175     {
 176         return
 177             integral (kernel1, x1, kernel2, scale, x2, - x2) +
 178             integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
 179     }
 180     else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
 181     {
 182         assert (width == 0.0);
 183         return filters[kernel2].func (x2 * scale);
 184     }
 185     else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
 186     {
 187         assert (width == 0.0);
 188         return filters[kernel1].func (x1);
 189     }
 190     else
 191     {
 192         /* Integration via Simpson's rule */
 193 #define N_SEGMENTS 128
 194 #define SAMPLE(a1, a2)                                                  \
 195         (filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
 196
 197         double s = 0.0;
 198         double h = width / (double)N_SEGMENTS;
 199         int i;
 200
 201         s = SAMPLE (x1, x2);
 202
 203         for (i = 1; i < N_SEGMENTS; i += 2)
 204         {
 205             double a1 = x1 + h * i;
 206             double a2 = x2 + h * i;
 207
 208             s += 2 * SAMPLE (a1, a2);
 209
 210             if (i >= 2 && i < N_SEGMENTS - 1)
 211                 s += 4 * SAMPLE (a1, a2);
 212         }
 213
 214         s += SAMPLE (x1 + width, x2 + width);
 215
 216         return h * s * (1.0 / 3.0);
 217     }
 218 }
 219
 220 static pixman_fixed_t *
 221 create_1d_filter (int             *width,
 222                   pixman_kernel_t  reconstruct,
 223                   pixman_kernel_t  sample,
 224                   double           scale,
 225                   int              n_phases)
 226 {
 227     pixman_fixed_t *params, *p;
 228     double step;
 229     double size;
 230     int i;
 231
 232     size = scale * filters[sample].width + filters[reconstruct].width;
 233     *width = ceil (size);
 234
 235     p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t));
 236     if (!params)
 237         return NULL;
 238
 239     step = 1.0 / n_phases;
 240
 241     for (i = 0; i < n_phases; ++i)
 242     {
 243         double frac = step / 2.0 + i * step;
 244         pixman_fixed_t new_total;
 245         int x, x1, x2;
 246         double total;
 247
 248         /* Sample convolution of reconstruction and sampling
 249          * filter. See rounding.txt regarding the rounding
 250          * and sample positions.
 251          */
 252
 253         x1 = ceil (frac - *width / 2.0 - 0.5);
 254         x2 = x1 + *width;
 255
 256         total = 0;
 257         for (x = x1; x < x2; ++x)
 258         {
 259             double pos = x + 0.5 - frac;
 260             double rlow = - filters[reconstruct].width / 2.0;
 261             double rhigh = rlow + filters[reconstruct].width;
 262             double slow = pos - scale * filters[sample].width / 2.0;
 263             double shigh = slow + scale * filters[sample].width;
 264             double c = 0.0;
 265             double ilow, ihigh;
 266
 267             if (rhigh >= slow && rlow <= shigh)
 268             {
 269                 ilow = MAX (slow, rlow);
 270                 ihigh = MIN (shigh, rhigh);
 271
 272                 c = integral (reconstruct, ilow,
 273                               sample, 1.0 / scale, ilow - pos,
 274                               ihigh - ilow);
 275             }
 276
 277             total += c;
 278             *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5);
 279         }
 280
 281         /* Normalize */
 282         p -= *width;
 283         total = 1 / total;
 284         new_total = 0;
 285         for (x = x1; x < x2; ++x)
 286         {
 287             pixman_fixed_t t = (*p) * total + 0.5;
 288
 289             new_total += t;
 290             *p++ = t;
 291         }
 292
 293         if (new_total != pixman_fixed_1)
 294             *(p - *width / 2) += (pixman_fixed_1 - new_total);
 295     }
 296
 297     return params;
 298 }
 299
 300 /* Create the parameter list for a SEPARABLE_CONVOLUTION filter
 301  * with the given kernels and scale parameters
 302  */
 303 PIXMAN_EXPORT pixman_fixed_t *
 304 pixman_filter_create_separable_convolution (int             *n_values,
 305                                             pixman_fixed_t   scale_x,
 306                                             pixman_fixed_t   scale_y,
 307                                             pixman_kernel_t  reconstruct_x,
 308                                             pixman_kernel_t  reconstruct_y,
 309                                             pixman_kernel_t  sample_x,
 310                                             pixman_kernel_t  sample_y,
 311                                             int              subsample_bits_x,
 312                                             int              subsample_bits_y)
 313 {
 314     double sx = fabs (pixman_fixed_to_double (scale_x));
 315     double sy = fabs (pixman_fixed_to_double (scale_y));
 316     pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL;
 317     int subsample_x, subsample_y;
 318     int width, height;
 319
 320     subsample_x = (1 << subsample_bits_x);
 321     subsample_y = (1 << subsample_bits_y);
 322
 323     horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x);
 324     vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y);
 325
 326     if (!horz || !vert)
 327         goto out;
 328
 329     *n_values = 4 + width * subsample_x + height * subsample_y;
 330
 331     params = malloc (*n_values * sizeof (pixman_fixed_t));
 332     if (!params)
 333         goto out;
 334
 335     params[0] = pixman_int_to_fixed (width);
 336     params[1] = pixman_int_to_fixed (height);
 337     params[2] = pixman_int_to_fixed (subsample_bits_x);
 338     params[3] = pixman_int_to_fixed (subsample_bits_y);
 339
 340     memcpy (params + 4, horz,
 341             width * subsample_x * sizeof (pixman_fixed_t));
 342     memcpy (params + 4 + width * subsample_x, vert,
 343             height * subsample_y * sizeof (pixman_fixed_t));
 344
 345 out:
 346     free (horz);
 347     free (vert);
 348
 349     return params;
 350 }