Home | History | Annotate | Download | only in pixman
      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 (at) gmail.com>
     25  */
     26 #include <string.h>
     27 #include <stdlib.h>
     28 #include <stdio.h>
     29 #include <math.h>
     30 #include <assert.h>
     31 #include <config.h>
     32 #include "pixman-private.h"
     33 
     34 typedef double (* kernel_func_t) (double x);
     35 
     36 typedef struct
     37 {
     38     pixman_kernel_t	kernel;
     39     kernel_func_t	func;
     40     double		width;
     41 } filter_info_t;
     42 
     43 static double
     44 impulse_kernel (double x)
     45 {
     46     return (x == 0.0)? 1.0 : 0.0;
     47 }
     48 
     49 static double
     50 box_kernel (double x)
     51 {
     52     return 1;
     53 }
     54 
     55 static double
     56 linear_kernel (double x)
     57 {
     58     return 1 - fabs (x);
     59 }
     60 
     61 static double
     62 gaussian_kernel (double x)
     63 {
     64 #define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
     65 #define SIGMA (SQRT2 / 2.0)
     66 
     67     return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
     68 }
     69 
     70 static double
     71 sinc (double x)
     72 {
     73     if (x == 0.0)
     74 	return 1.0;
     75     else
     76 	return sin (M_PI * x) / (M_PI * x);
     77 }
     78 
     79 static double
     80 lanczos (double x, int n)
     81 {
     82     return sinc (x) * sinc (x * (1.0 / n));
     83 }
     84 
     85 static double
     86 lanczos2_kernel (double x)
     87 {
     88     return lanczos (x, 2);
     89 }
     90 
     91 static double
     92 lanczos3_kernel (double x)
     93 {
     94     return lanczos (x, 3);
     95 }
     96 
     97 static double
     98 nice_kernel (double x)
     99 {
    100     return lanczos3_kernel (x * 0.75);
    101 }
    102 
    103 static double
    104 general_cubic (double x, double B, double C)
    105 {
    106     double ax = fabs(x);
    107 
    108     if (ax < 1)
    109     {
    110 	return ((12 - 9 * B - 6 * C) * ax * ax * ax +
    111 		(-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6;
    112     }
    113     else if (ax >= 1 && ax < 2)
    114     {
    115 	return ((-B - 6 * C) * ax * ax * ax +
    116 		(6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) *
    117 		ax + (8 * B + 24 * C)) / 6;
    118     }
    119     else
    120     {
    121 	return 0;
    122     }
    123 }
    124 
    125 static double
    126 cubic_kernel (double x)
    127 {
    128     /* This is the Mitchell-Netravali filter.
    129      *
    130      * (0.0, 0.5) would give us the Catmull-Rom spline,
    131      * but that one seems to be indistinguishable from Lanczos2.
    132      */
    133     return general_cubic (x, 1/3.0, 1/3.0);
    134 }
    135 
    136 static const filter_info_t filters[] =
    137 {
    138     { PIXMAN_KERNEL_IMPULSE,	        impulse_kernel,   0.0 },
    139     { PIXMAN_KERNEL_BOX,	        box_kernel,       1.0 },
    140     { PIXMAN_KERNEL_LINEAR,	        linear_kernel,    2.0 },
    141     { PIXMAN_KERNEL_CUBIC,		cubic_kernel,     4.0 },
    142     { PIXMAN_KERNEL_GAUSSIAN,	        gaussian_kernel,  6 * SIGMA },
    143     { PIXMAN_KERNEL_LANCZOS2,	        lanczos2_kernel,  4.0 },
    144     { PIXMAN_KERNEL_LANCZOS3,	        lanczos3_kernel,  6.0 },
    145     { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
    146 };
    147 
    148 /* This function scales @kernel2 by @scale, then
    149  * aligns @x1 in @kernel1 with @x2 in @kernel2 and
    150  * and integrates the product of the kernels across @width.
    151  *
    152  * This function assumes that the intervals are within
    153  * the kernels in question. E.g., the caller must not
    154  * try to integrate a linear kernel ouside of [-1:1]
    155  */
    156 static double
    157 integral (pixman_kernel_t kernel1, double x1,
    158 	  pixman_kernel_t kernel2, double scale, double x2,
    159 	  double width)
    160 {
    161     /* If the integration interval crosses zero, break it into
    162      * two separate integrals. This ensures that filters such
    163      * as LINEAR that are not differentiable at 0 will still
    164      * integrate properly.
    165      */
    166     if (x1 < 0 && x1 + width > 0)
    167     {
    168 	return
    169 	    integral (kernel1, x1, kernel2, scale, x2, - x1) +
    170 	    integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
    171     }
    172     else if (x2 < 0 && x2 + width > 0)
    173     {
    174 	return
    175 	    integral (kernel1, x1, kernel2, scale, x2, - x2) +
    176 	    integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
    177     }
    178     else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
    179     {
    180 	assert (width == 0.0);
    181 	return filters[kernel2].func (x2 * scale);
    182     }
    183     else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
    184     {
    185 	assert (width == 0.0);
    186 	return filters[kernel1].func (x1);
    187     }
    188     else
    189     {
    190 	/* Integration via Simpson's rule */
    191 #define N_SEGMENTS 128
    192 #define SAMPLE(a1, a2)							\
    193 	(filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
    194 
    195 	double s = 0.0;
    196 	double h = width / (double)N_SEGMENTS;
    197 	int i;
    198 
    199 	s = SAMPLE (x1, x2);
    200 
    201 	for (i = 1; i < N_SEGMENTS; i += 2)
    202 	{
    203 	    double a1 = x1 + h * i;
    204 	    double a2 = x2 + h * i;
    205 
    206 	    s += 2 * SAMPLE (a1, a2);
    207 
    208 	    if (i >= 2 && i < N_SEGMENTS - 1)
    209 		s += 4 * SAMPLE (a1, a2);
    210 	}
    211 
    212 	s += SAMPLE (x1 + width, x2 + width);
    213 
    214 	return h * s * (1.0 / 3.0);
    215     }
    216 }
    217 
    218 static pixman_fixed_t *
    219 create_1d_filter (int             *width,
    220 		  pixman_kernel_t  reconstruct,
    221 		  pixman_kernel_t  sample,
    222 		  double           scale,
    223 		  int              n_phases)
    224 {
    225     pixman_fixed_t *params, *p;
    226     double step;
    227     double size;
    228     int i;
    229 
    230     size = scale * filters[sample].width + filters[reconstruct].width;
    231     *width = ceil (size);
    232 
    233     p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t));
    234     if (!params)
    235         return NULL;
    236 
    237     step = 1.0 / n_phases;
    238 
    239     for (i = 0; i < n_phases; ++i)
    240     {
    241         double frac = step / 2.0 + i * step;
    242 	pixman_fixed_t new_total;
    243         int x, x1, x2;
    244 	double total;
    245 
    246 	/* Sample convolution of reconstruction and sampling
    247 	 * filter. See rounding.txt regarding the rounding
    248 	 * and sample positions.
    249 	 */
    250 
    251 	x1 = ceil (frac - *width / 2.0 - 0.5);
    252         x2 = x1 + *width;
    253 
    254 	total = 0;
    255         for (x = x1; x < x2; ++x)
    256         {
    257 	    double pos = x + 0.5 - frac;
    258 	    double rlow = - filters[reconstruct].width / 2.0;
    259 	    double rhigh = rlow + filters[reconstruct].width;
    260 	    double slow = pos - scale * filters[sample].width / 2.0;
    261 	    double shigh = slow + scale * filters[sample].width;
    262 	    double c = 0.0;
    263 	    double ilow, ihigh;
    264 
    265 	    if (rhigh >= slow && rlow <= shigh)
    266 	    {
    267 		ilow = MAX (slow, rlow);
    268 		ihigh = MIN (shigh, rhigh);
    269 
    270 		c = integral (reconstruct, ilow,
    271 			      sample, 1.0 / scale, ilow - pos,
    272 			      ihigh - ilow);
    273 	    }
    274 
    275 	    total += c;
    276             *p++ = (pixman_fixed_t)(c * 65535.0 + 0.5);
    277         }
    278 
    279 	/* Normalize */
    280 	p -= *width;
    281         total = 1 / total;
    282         new_total = 0;
    283 	for (x = x1; x < x2; ++x)
    284 	{
    285 	    pixman_fixed_t t = (*p) * total + 0.5;
    286 
    287 	    new_total += t;
    288 	    *p++ = t;
    289 	}
    290 
    291 	if (new_total != pixman_fixed_1)
    292 	    *(p - *width / 2) += (pixman_fixed_1 - new_total);
    293     }
    294 
    295     return params;
    296 }
    297 
    298 /* Create the parameter list for a SEPARABLE_CONVOLUTION filter
    299  * with the given kernels and scale parameters
    300  */
    301 PIXMAN_EXPORT pixman_fixed_t *
    302 pixman_filter_create_separable_convolution (int             *n_values,
    303 					    pixman_fixed_t   scale_x,
    304 					    pixman_fixed_t   scale_y,
    305 					    pixman_kernel_t  reconstruct_x,
    306 					    pixman_kernel_t  reconstruct_y,
    307 					    pixman_kernel_t  sample_x,
    308 					    pixman_kernel_t  sample_y,
    309 					    int              subsample_bits_x,
    310 					    int	             subsample_bits_y)
    311 {
    312     double sx = fabs (pixman_fixed_to_double (scale_x));
    313     double sy = fabs (pixman_fixed_to_double (scale_y));
    314     pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL;
    315     int subsample_x, subsample_y;
    316     int width, height;
    317 
    318     subsample_x = (1 << subsample_bits_x);
    319     subsample_y = (1 << subsample_bits_y);
    320 
    321     horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x);
    322     vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y);
    323 
    324     if (!horz || !vert)
    325         goto out;
    326 
    327     *n_values = 4 + width * subsample_x + height * subsample_y;
    328 
    329     params = malloc (*n_values * sizeof (pixman_fixed_t));
    330     if (!params)
    331         goto out;
    332 
    333     params[0] = pixman_int_to_fixed (width);
    334     params[1] = pixman_int_to_fixed (height);
    335     params[2] = pixman_int_to_fixed (subsample_bits_x);
    336     params[3] = pixman_int_to_fixed (subsample_bits_y);
    337 
    338     memcpy (params + 4, horz,
    339 	    width * subsample_x * sizeof (pixman_fixed_t));
    340     memcpy (params + 4 + width * subsample_x, vert,
    341 	    height * subsample_y * sizeof (pixman_fixed_t));
    342 
    343 out:
    344     free (horz);
    345     free (vert);
    346 
    347     return params;
    348 }
    349