pixman-filter.c (13306B)
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 <pixman-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 *pstart, 241 pixman_fixed_t *pend 242 ) 243 { 244 pixman_fixed_t *p = pstart; 245 double step; 246 int i; 247 if(width <= 0) return; 248 step = 1.0 / n_phases; 249 250 for (i = 0; i < n_phases; ++i) 251 { 252 double frac = step / 2.0 + i * step; 253 pixman_fixed_t new_total; 254 int x, x1, x2; 255 double total, e; 256 257 /* Sample convolution of reconstruction and sampling 258 * filter. See rounding.txt regarding the rounding 259 * and sample positions. 260 */ 261 262 x1 = ceil (frac - width / 2.0 - 0.5); 263 x2 = x1 + width; 264 assert( p >= pstart && p + (x2 - x1) <= pend ); /* assert validity of the following loop */ 265 total = 0; 266 for (x = x1; x < x2; ++x) 267 { 268 double pos = x + 0.5 - frac; 269 double rlow = - filters[reconstruct].width / 2.0; 270 double rhigh = rlow + filters[reconstruct].width; 271 double slow = pos - scale * filters[sample].width / 2.0; 272 double shigh = slow + scale * filters[sample].width; 273 double c = 0.0; 274 double ilow, ihigh; 275 276 if (rhigh >= slow && rlow <= shigh) 277 { 278 ilow = MAX (slow, rlow); 279 ihigh = MIN (shigh, rhigh); 280 281 c = integral (reconstruct, ilow, 282 sample, 1.0 / scale, ilow - pos, 283 ihigh - ilow); 284 } 285 286 *p = (pixman_fixed_t)floor (c * 65536.0 + 0.5); 287 total += *p; 288 p++; 289 } 290 291 /* Normalize, with error diffusion */ 292 p -= width; 293 assert(p >= pstart && p + (x2 - x1) <= pend); /* assert validity of the following loop */ 294 295 total = 65536.0 / total; 296 new_total = 0; 297 e = 0.0; 298 for (x = x1; x < x2; ++x) 299 { 300 double v = (*p) * total + e; 301 pixman_fixed_t t = floor (v + 0.5); 302 303 e = v - t; 304 new_total += t; 305 *p++ = t; 306 } 307 308 /* pixman_fixed_e's worth of error may remain; put it 309 * at the first sample, since that is the only one that 310 * hasn't had any error diffused into it. 311 */ 312 313 assert(p - width >= pstart && p - width < pend); /* assert... */ 314 *(p - width) += pixman_fixed_1 - new_total; 315 } 316 } 317 318 319 static int 320 filter_width (pixman_kernel_t reconstruct, pixman_kernel_t sample, double size) 321 { 322 return ceil (filters[reconstruct].width + size * filters[sample].width); 323 } 324 325 #ifdef PIXMAN_GNUPLOT 326 327 /* If enable-gnuplot is configured, then you can pipe the output of a 328 * pixman-using program to gnuplot and get a continuously-updated plot 329 * of the horizontal filter. This works well with demos/scale to test 330 * the filter generation. 331 * 332 * The plot is all the different subposition filters shuffled 333 * together. This is misleading in a few cases: 334 * 335 * IMPULSE.BOX - goes up and down as the subfilters have different 336 * numbers of non-zero samples 337 * IMPULSE.TRIANGLE - somewhat crooked for the same reason 338 * 1-wide filters - looks triangular, but a 1-wide box would be more 339 * accurate 340 */ 341 static void 342 gnuplot_filter (int width, int n_phases, const pixman_fixed_t* p) 343 { 344 double step; 345 int i, j; 346 int first; 347 348 step = 1.0 / n_phases; 349 350 printf ("set style line 1 lc rgb '#0060ad' lt 1 lw 0.5 pt 7 pi 1 ps 0.5\n"); 351 printf ("plot [x=%g:%g] '-' with linespoints ls 1\n", -width*0.5, width*0.5); 352 /* Print a point at the origin so that y==0 line is included: */ 353 printf ("0 0\n\n"); 354 355 /* The position of the first sample of the phase corresponding to 356 * frac is given by: 357 * 358 * ceil (frac - width / 2.0 - 0.5) + 0.5 - frac 359 * 360 * We have to find the frac that minimizes this expression. 361 * 362 * For odd widths, we have 363 * 364 * ceil (frac - width / 2.0 - 0.5) + 0.5 - frac 365 * = ceil (frac) + K - frac 366 * = 1 + K - frac 367 * 368 * for some K, so this is minimized when frac is maximized and 369 * strictly growing with frac. So for odd widths, we can simply 370 * start at the last phase and go backwards. 371 * 372 * For even widths, we have 373 * 374 * ceil (frac - width / 2.0 - 0.5) + 0.5 - frac 375 * = ceil (frac - 0.5) + K - frac 376 * 377 * The graph for this function (ignoring K) looks like this: 378 * 379 * 0.5 380 * | |\ 381 * | | \ 382 * | | \ 383 * 0 | | \ 384 * |\ | 385 * | \ | 386 * | \ | 387 * -0.5 | \| 388 * --------------------------------- 389 * 0 0.5 1 390 * 391 * So in this case we need to start with the phase whose frac is 392 * less than, but as close as possible to 0.5, then go backwards 393 * until we hit the first phase, then wrap around to the last 394 * phase and continue backwards. 395 * 396 * Which phase is as close as possible 0.5? The locations of the 397 * sampling point corresponding to the kth phase is given by 398 * 1/(2 * n_phases) + k / n_phases: 399 * 400 * 1/(2 * n_phases) + k / n_phases = 0.5 401 * 402 * from which it follows that 403 * 404 * k = (n_phases - 1) / 2 405 * 406 * rounded down is the phase in question. 407 */ 408 if (width & 1) 409 first = n_phases - 1; 410 else 411 first = (n_phases - 1) / 2; 412 413 for (j = 0; j < width; ++j) 414 { 415 for (i = 0; i < n_phases; ++i) 416 { 417 int phase = first - i; 418 double frac, pos; 419 420 if (phase < 0) 421 phase = n_phases + phase; 422 423 frac = step / 2.0 + phase * step; 424 pos = ceil (frac - width / 2.0 - 0.5) + 0.5 - frac + j; 425 426 printf ("%g %g\n", 427 pos, 428 pixman_fixed_to_double (*(p + phase * width + j))); 429 } 430 } 431 432 printf ("e\n"); 433 fflush (stdout); 434 } 435 436 #endif 437 438 /* Create the parameter list for a SEPARABLE_CONVOLUTION filter 439 * with the given kernels and scale parameters 440 */ 441 PIXMAN_EXPORT pixman_fixed_t * 442 pixman_filter_create_separable_convolution (int *n_values, 443 pixman_fixed_t scale_x, 444 pixman_fixed_t scale_y, 445 pixman_kernel_t reconstruct_x, 446 pixman_kernel_t reconstruct_y, 447 pixman_kernel_t sample_x, 448 pixman_kernel_t sample_y, 449 int subsample_bits_x, 450 int subsample_bits_y) 451 { 452 double sx = fabs (pixman_fixed_to_double (scale_x)); 453 double sy = fabs (pixman_fixed_to_double (scale_y)); 454 pixman_fixed_t *params; 455 int subsample_x, subsample_y; 456 int width, height; 457 458 width = filter_width (reconstruct_x, sample_x, sx); 459 subsample_x = (1 << subsample_bits_x); 460 461 height = filter_width (reconstruct_y, sample_y, sy); 462 subsample_y = (1 << subsample_bits_y); 463 464 *n_values = 4 + width * subsample_x + height * subsample_y; 465 466 params = malloc (*n_values * sizeof (pixman_fixed_t)); 467 if (!params) 468 return NULL; 469 470 params[0] = pixman_int_to_fixed (width); 471 params[1] = pixman_int_to_fixed (height); 472 params[2] = pixman_int_to_fixed (subsample_bits_x); 473 params[3] = pixman_int_to_fixed (subsample_bits_y); 474 475 { 476 pixman_fixed_t 477 *xparams = params+4, 478 *yparams = xparams + width*subsample_x, 479 *endparams = params + *n_values; 480 create_1d_filter(width, reconstruct_x, sample_x, sx, subsample_x, 481 xparams, yparams); 482 create_1d_filter(height, reconstruct_y, sample_y, sy, subsample_y, 483 yparams, endparams); 484 } 485 486 #ifdef PIXMAN_GNUPLOT 487 gnuplot_filter(width, subsample_x, params + 4); 488 #endif 489 490 return params; 491 }