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
