Uses filter evaluators that know the size of a pixel, and thus can integrate with reconstruction filters directly. The result is equivalent to having the box and impulse sampling filters use a box reconstruction, and all others use box for sizes smaller than their Nyquist rate and impulse for larger sizes.
For back-compatibility the "reconstruction" filter argument is used at a size of 1 if it is not box or impulse and the filter size is less than 1. IMHO this preserves all the useful results of the previous code, is a good deal simpler, and the caller no longer has to alter the arguments depending on whether the size is less than one or not. --- pixman/pixman-filter.c | 392 ++++++++++++++++++++++++------------------------ 1 file changed, 194 insertions(+), 198 deletions(-) diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c index b2bf53f..1087f6d 100644 --- a/pixman/pixman-filter.c +++ b/pixman/pixman-filter.c @@ -33,42 +33,83 @@ #endif #include "pixman-private.h" -typedef double (* kernel_func_t) (double x); +/* Produce contribution of a filter of size r for pixel centered on x. + * For a typical low-pass function this evaluates the function at x/r. + * If the frequency is higher than 1/2, such as when r is less than 1, + * this may need to integrate several samples, see cubic for examples. + */ +typedef double (* kernel_func_t) (double x, double r); typedef struct { pixman_kernel_t kernel; kernel_func_t func; - double width; + double width; /* note special handling if <= 1 */ } filter_info_t; +/* PIXMAN_KERNEL_IMPULSE: Returns pixel nearest the center. This + * matches PIXMAN_FILTER_NEAREST. This is useful if you wish to + * combine the result of nearest in one direction with another filter + * in the other. + */ static double -impulse_kernel (double x) +impulse_kernel (double x, double r) { - return (x == 0.0)? 1.0 : 0.0; + return 1; } +/* PIXMAN_KERNEL_BOX: Intersection of a box of width r with square + * pixels. This is the smallest possible filter such that the output + * image contains an equal contribution from all the input + * pixels. Lots of software uses this. The function is a trapazoid of + * width r+1, not a box. + * + * When r == 1.0, PIXMAN_KERNEL_BOX and PIXMAN_KERNEL_LINEAR produce + * the same filter, allowing them to be exchanged at this point. + */ static double -box_kernel (double x) +box_kernel (double x, double r) { - return 1; + return MAX (0.0, MIN (MIN (r, 1.0), + MIN ((r + 1) / 2 - x, (r + 1) / 2 + x))); } +/* PIXMAN_KERNEL_LINEAR: Triangle of width 2r. Lots of software uses + * this as a "better" filter, twice the size of a box but smaller than + * a cubic. + * + * When r == 1.0, PIXMAN_KERNEL_BOX and PIXMAN_KERNEL_LINEAR produce + * the same filter, allowing them to be exchanged at this point. + */ static double -linear_kernel (double x) +linear_kernel (double x, double r) { - return 1 - fabs (x); + if (r < 1.0) + return box_kernel(x, r); + else + return MAX (1.0 - fabs(x / r), 0.0); } +/* PIXMAN_KERNEL_GAUSSIAN: Gaussian function, cut off somewhere near + +/-2.5. Output looks like it has been blurred by a gaussian filter + of size 1. This may be useful if the output is to be blurred + further as that blur can be reduced by 1. +*/ static double -gaussian_kernel (double x) +gaussian_kernel (double x, double r) { -#define SQRT2 (1.4142135623730950488016887242096980785696718753769480) -#define SIGMA (SQRT2 / 2.0) - - return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI)); + if (r < 0.5) + return + gaussian_kernel (x * 2 - .5, r * 2) + + gaussian_kernel (x * 2 + .5, r * 2); + return exp (- x * x) / sqrt (M_PI); } +/* PIXMAN_KERNEL_LANCZOS2: lanczos windowed sinc function from -2 to + * +2. This is very close to the cubic filter that is often used + * by other software. Should be slightly higher quality. + */ + static double sinc (double x) { @@ -79,60 +120,91 @@ sinc (double x) } static double -lanczos (double x, int n) +lanczos (double x, double r, double n) { - return sinc (x) * sinc (x * (1.0 / n)); + if (r < 1.0) + return + lanczos (x * 2 - .5, r * 2, n) + + lanczos (x * 2 + .5, r * 2, n); + x /= r; + return fabs (x) < n ? sinc (x) * sinc (x * (1.0 / n)) : 0.0; } static double -lanczos2_kernel (double x) +lanczos2_kernel (double x, double r) { - return lanczos (x, 2); + return lanczos (x, r, 2.0); } +/* PIXMAN_KERNEL_LANCZOS3: lanczos windowed sinc function from -3 to + * +3. Very popular with high-end software though I think any + * advantage over cubics is hidden by quantization and programming + * mistakes. You will see LANCZOS5 or even 7 sometimes. + */ static double -lanczos3_kernel (double x) +lanczos3_kernel (double x, double r) { - return lanczos (x, 3); + return lanczos (x, r, 3.0); } +/* PIXMAN_KERNEL_LANCZOS3_STRETCHED - The LANCZOS3 kernel widened by + * 4/3. Recommended by Jim Blinn + * http://graphics.cs.cmu.edu/nsp/course/15-462/Fall07/462/papers/jaggy.pdf + */ static double -nice_kernel (double x) +nice_kernel (double x, double r) { - return lanczos3_kernel (x * 0.75); + return lanczos3_kernel (x, r * (4.0/3)); } +/* Cubic functions described in the Mitchell-Netravali paper. + * http://mentallandscape.com/Papers_siggraph88.pdf. This describes + * all possible cubic functions that can be used for sampling. + */ static double -general_cubic (double x, double B, double C) +general_cubic (double x, double r, double B, double C) { - double ax = fabs(x); + double ax; + if (r < 1.0) + return + general_cubic(x * 2 - .5, r * 2, B, C) + + general_cubic(x * 2 + .5, r * 2, B, C); + + ax = fabs (x / r); if (ax < 1) { - return ((12 - 9 * B - 6 * C) * ax * ax * ax + - (-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6; + return (((12 - 9 * B - 6 * C) * ax + + (-18 + 12 * B + 6 * C)) * ax * ax + + (6 - 2 * B)) / 6; } - else if (ax >= 1 && ax < 2) + else if (ax < 2) { - return ((-B - 6 * C) * ax * ax * ax + - (6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) * - ax + (8 * B + 24 * C)) / 6; + return ((((-B - 6 * C) * ax + + (6 * B + 30 * C)) * ax + + (-12 * B - 48 * C)) * ax + + (8 * B + 24 * C)) / 6; } else { - return 0; + return 0.0; } } +/* PIXMAN_KERNEL_CUBIC: Cubic recommended by the Mitchell-Netravali + * paper. This has negative values and because the values at +/-1 are + * not zero it does not interpolate the pixels, meaning it will change + * an image even if there is no translation. + * + * Warning: this is not what a lot of other software calls "cubic". + * That often refers to the Catmull-Rom spline, or another interpolating + * function. (0.0, 0.5) would give us the Catmull-Rom spline, but + * that one seems indistinguishable from Lanczos2. + */ static double -cubic_kernel (double x) +cubic_kernel (double x, double r) { - /* This is the Mitchell-Netravali filter. - * - * (0.0, 0.5) would give us the Catmull-Rom spline, - * but that one seems to be indistinguishable from Lanczos2. - */ - return general_cubic (x, 1/3.0, 1/3.0); + return general_cubic (x, r, 1/3.0, 1/3.0); } static const filter_info_t filters[] = @@ -140,163 +212,62 @@ static const filter_info_t filters[] = { PIXMAN_KERNEL_IMPULSE, impulse_kernel, 0.0 }, { PIXMAN_KERNEL_BOX, box_kernel, 1.0 }, { PIXMAN_KERNEL_LINEAR, linear_kernel, 2.0 }, - { PIXMAN_KERNEL_CUBIC, cubic_kernel, 4.0 }, - { PIXMAN_KERNEL_GAUSSIAN, gaussian_kernel, 6 * SIGMA }, + { PIXMAN_KERNEL_CUBIC, cubic_kernel, 4.0 }, + { PIXMAN_KERNEL_GAUSSIAN, gaussian_kernel, 5.0 }, { PIXMAN_KERNEL_LANCZOS2, lanczos2_kernel, 4.0 }, { PIXMAN_KERNEL_LANCZOS3, lanczos3_kernel, 6.0 }, { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel, 8.0 }, }; -/* This function scales @kernel2 by @scale, then - * aligns @x1 in @kernel1 with @x2 in @kernel2 and - * and integrates the product of the kernels across @width. - * - * This function assumes that the intervals are within - * the kernels in question. E.g., the caller must not - * try to integrate a linear kernel ouside of [-1:1] - */ -static double -integral (pixman_kernel_t kernel1, double x1, - pixman_kernel_t kernel2, double scale, double x2, - double width) -{ - /* If the integration interval crosses zero, break it into - * two separate integrals. This ensures that filters such - * as LINEAR that are not differentiable at 0 will still - * integrate properly. - */ - if (x1 < 0 && x1 + width > 0) - { - return - integral (kernel1, x1, kernel2, scale, x2, - x1) + - integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1); - } - else if (x2 < 0 && x2 + width > 0) - { - return - integral (kernel1, x1, kernel2, scale, x2, - x2) + - integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2); - } - else if (kernel1 == PIXMAN_KERNEL_IMPULSE) - { - assert (width == 0.0); - return filters[kernel2].func (x2 * scale); - } - else if (kernel2 == PIXMAN_KERNEL_IMPULSE) - { - assert (width == 0.0); - return filters[kernel1].func (x1); - } - else - { - /* Integration via Simpson's rule */ -#define N_SEGMENTS 128 -#define SAMPLE(a1, a2) \ - (filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale)) - - double s = 0.0; - double h = width / (double)N_SEGMENTS; - int i; - - s = SAMPLE (x1, x2); - - for (i = 1; i < N_SEGMENTS; i += 2) - { - double a1 = x1 + h * i; - double a2 = x2 + h * i; - s += 2 * SAMPLE (a1, a2); - - if (i >= 2 && i < N_SEGMENTS - 1) - s += 4 * SAMPLE (a1, a2); - } - - s += SAMPLE (x1 + width, x2 + width); - - return h * s * (1.0 / 3.0); - } -} - -static pixman_fixed_t * -create_1d_filter (int *width, - pixman_kernel_t reconstruct, - pixman_kernel_t sample, - double scale, - int n_phases) +/* Fills in one dimension of the filter array */ +static void get_filter(pixman_kernel_t filter, double r, + int width, int subsample, + pixman_fixed_t* out) { - pixman_fixed_t *params, *p; - double step; - double size; int i; + pixman_fixed_t *p = out; + int n_phases = 1 << subsample; + double step = 1.0 / n_phases; + kernel_func_t func = filters[filter].func; - size = scale * filters[sample].width + filters[reconstruct].width; - *width = ceil (size); - - p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t)); - if (!params) - return NULL; - - step = 1.0 / n_phases; + /* special-case the impulse filter: */ + if (width <= 1) + { + for (i = 0; i < n_phases; ++i) + *p++ = pixman_fixed_1; + return; + } for (i = 0; i < n_phases; ++i) { - double frac = step / 2.0 + i * step; - pixman_fixed_t new_total; - int x, x1, x2; - double total; - - /* Sample convolution of reconstruction and sampling - * filter. See rounding.txt regarding the rounding - * and sample positions. - */ - - x1 = ceil (frac - *width / 2.0 - 0.5); - x2 = x1 + *width; - - total = 0; - for (x = x1; x < x2; ++x) - { - double pos = x + 0.5 - frac; - double rlow = - filters[reconstruct].width / 2.0; - double rhigh = rlow + filters[reconstruct].width; - double slow = pos - scale * filters[sample].width / 2.0; - double shigh = slow + scale * filters[sample].width; - double c = 0.0; - double ilow, ihigh; - - if (rhigh >= slow && rlow <= shigh) - { - ilow = MAX (slow, rlow); - ihigh = MIN (shigh, rhigh); - - c = integral (reconstruct, ilow, - sample, 1.0 / scale, ilow - pos, - ihigh - ilow); - } - - total += c; - *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5); - } + double frac = (i + .5) * step; + /* Center of left-most pixel: */ + double x1 = ceil (frac - width / 2.0 - 0.5) - frac + 0.5; + double total = 0; + pixman_fixed_t new_total = 0; + int j; + + for (j = 0; j < width; ++j) + { + double v = func(x1 + j, r); + total += v; + p[j] = pixman_double_to_fixed (v); + } /* Normalize */ - p -= *width; total = 1 / total; - new_total = 0; - for (x = x1; x < x2; ++x) - { - pixman_fixed_t t = (*p) * total + 0.5; + for (j = 0; j < width; ++j) + new_total += (p[j] *= total); - new_total += t; - *p++ = t; - } + /* Put any error on center pixel */ + p[width / 2] += (pixman_fixed_1 - new_total); - if (new_total != pixman_fixed_1) - *(p - *width / 2) += (pixman_fixed_1 - new_total); + p += width; } - - return params; } + /* Create the parameter list for a SEPARABLE_CONVOLUTION filter * with the given kernels and scale parameters */ @@ -313,38 +284,63 @@ pixman_filter_create_separable_convolution (int *n_values, { double sx = fabs (pixman_fixed_to_double (scale_x)); double sy = fabs (pixman_fixed_to_double (scale_y)); - pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL; - int subsample_x, subsample_y; - int width, height; + int width_x, width_y, size_x, size_y; + pixman_fixed_t *params; - subsample_x = (1 << subsample_bits_x); - subsample_y = (1 << subsample_bits_y); + /* Simulate older version of this function: + * + * Old version attempted to convolve two filters, the "reconstruct" + * filter at a scale of 1 and the main filter at the given scale. + * + * These filter generators simulate a BOX or IMPULSE + * reconstruction filter depending on the scale and other + * factors. If another filter is specified as the reconstruction + * filter, then it is used at scale = 1 if the scale <= 1 to + * simulate the previous results. Since these filters are + * supposed to be low pass filters, the convolution is achieved by + * just choosing the larger one. + */ + if (reconstruct_x > PIXMAN_KERNEL_BOX && sx < 1.0) + { + sample_x = reconstruct_x; + sx = 1.0; + } + if (reconstruct_y > PIXMAN_KERNEL_BOX && sy < 1.0) + { + sample_y = reconstruct_y; + sy = 1.0; + } + + width_x = filters[sample_x] . width; + if (width_x < 1.0) + { width_x = 1.0; subsample_bits_x = 0; } + else if (width_x < 2.0) + width_x = sx < 1.0 ? 2 : ceil (sx + width_x); + else + width_x = MAX (2, ceil (sx * width_x)); + + width_y = filters[sample_y] . width; + if (width_y < 1.0) + { width_y = 1.0; subsample_bits_y = 0; } + else if (width_y < 2.0) + width_y = sy < 1.0 ? 2 : ceil (sy + width_y); + else + width_y = MAX (2, ceil (sy * width_y)); - horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x); - vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y); + size_x = (1 << subsample_bits_x) * width_x; + size_y = (1 << subsample_bits_y) * width_y; + *n_values = 4 + size_x + size_y; - if (!horz || !vert) - goto out; - - *n_values = 4 + width * subsample_x + height * subsample_y; - params = malloc (*n_values * sizeof (pixman_fixed_t)); - if (!params) - goto out; + if (!params) return 0; - params[0] = pixman_int_to_fixed (width); - params[1] = pixman_int_to_fixed (height); + params[0] = pixman_int_to_fixed (width_x); + params[1] = pixman_int_to_fixed (width_y); params[2] = pixman_int_to_fixed (subsample_bits_x); params[3] = pixman_int_to_fixed (subsample_bits_y); - memcpy (params + 4, horz, - width * subsample_x * sizeof (pixman_fixed_t)); - memcpy (params + 4 + width * subsample_x, vert, - height * subsample_y * sizeof (pixman_fixed_t)); - -out: - free (horz); - free (vert); + get_filter(sample_x, sx, width_x, subsample_bits_x, params + 4); + get_filter(sample_y, sy, width_y, subsample_bits_y, params + 4 + size_x); return params; } -- 1.7.9.5 _______________________________________________ Pixman mailing list Pixman@lists.freedesktop.org http://lists.freedesktop.org/mailman/listinfo/pixman