[Pixman] [PATCH 07/13] pixman-filter: Corrections to the integral() function
spitzak at gmail.com
spitzak at gmail.com
Sun Jan 3 19:12:11 PST 2016
From: Bill Spitzak <spitzak at gmail.com>
The IMPULSE special-cases did not sample the center of the of the region.
This caused it to sample the filters outside their range, and produce
assymetric filters and other errors. Fixing this required changing the
arguments to integral() so the correct point could be determined.
I fixed the nice filter and the integration to directly produce normalized
values. Re-normalization is still needed for impulse.box or impulse.triangle
so I did not remove it.
Distribute fixed error over all filter samples, to remove a high-frequency
bit of noise in the center of some filters (lancoz at large scale value).
I removed the recursion where it splits at zero. This would only be necessary
if the curve was not smooth at this point, but it is for all integrations
(even linear.linear is smooth). Plots show that removing this makes no difference.
box.box, which I expect will be very common as it is the proposed "good" filter,
was made a lot faster and more accurate. This is easy as the caller already
intersected the two boxes, so the width is the integral.
---
pixman/pixman-filter.c | 126 +++++++++++++++++++++----------------------------
1 file changed, 55 insertions(+), 71 deletions(-)
diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index 5677431..0749e51 100644
--- a/pixman/pixman-filter.c
+++ b/pixman/pixman-filter.c
@@ -79,7 +79,7 @@ sinc (double x)
}
static double
-lanczos (double x, int n)
+lanczos (double x, double n)
{
return sinc (x) * sinc (x * (1.0 / n));
}
@@ -99,7 +99,7 @@ lanczos3_kernel (double x)
static double
nice_kernel (double x)
{
- return lanczos3_kernel (x * 0.75);
+ return lanczos3_kernel (x * 0.75) * 0.75;
}
static double
@@ -147,45 +147,34 @@ static const filter_info_t filters[] =
{ PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel, 8.0 },
};
-/* This function scales @sample by @scale, then
- * aligns @x1 in @reconstruct with @x2 in @sample and
- * and integrates the product of the kernels across @width.
+/* This function scales @sample by @scale, shifts it by @pos,
+ * and then integrates the product of the kernels across low..high
*
* 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 reconstruct, double x1,
- pixman_kernel_t sample, double scale, double x2,
- double width)
+integral (pixman_kernel_t reconstruct,
+ pixman_kernel_t sample, double scale, double pos,
+ double low, double high)
{
- /* 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)
+ if (high < low)
{
- return
- integral (reconstruct, x1, sample, scale, x2, - x1) +
- integral (reconstruct, 0, sample, scale, x2 - x1, width + x1);
+ return 0.0;
}
- else if (x2 < 0 && x2 + width > 0)
+ else if (sample == PIXMAN_KERNEL_IMPULSE)
{
- return
- integral (reconstruct, x1, sample, scale, x2, - x2) +
- integral (reconstruct, x1 - x2, sample, scale, 0, width + x2);
+ return filters[reconstruct].func (-pos);
}
else if (reconstruct == PIXMAN_KERNEL_IMPULSE)
{
- assert (width == 0.0);
- return filters[sample].func (x2 / scale);
+ return filters[sample].func (-pos / scale) / scale;
}
- else if (sample == PIXMAN_KERNEL_IMPULSE)
+ else if (reconstruct == PIXMAN_KERNEL_BOX && sample == PIXMAN_KERNEL_BOX)
{
- assert (width == 0.0);
- return filters[reconstruct].func (x1);
+ assert (high <= low + 1.0);
+ return (high - low) / scale;
}
else
{
@@ -193,32 +182,30 @@ integral (pixman_kernel_t reconstruct, double x1,
* See http://www.intmath.com/integration/6-simpsons-rule.php
*/
#define N_SEGMENTS 16
-#define SAMPLE(a1, a2) \
- (filters[reconstruct].func ((a1)) * filters[sample].func ((a2) / scale))
+#define SAMPLE(a) \
+ (filters[reconstruct].func ((a)) * filters[sample].func (((a) - pos) / scale))
double s = 0.0;
- double h = width / (double)N_SEGMENTS;
+ double h = (high - low) / (double)N_SEGMENTS;
int i;
- s = SAMPLE (x1, x2);
+ s = SAMPLE (low);
for (i = 1; i < N_SEGMENTS; i += 2)
{
- double a1 = x1 + h * i;
- double a2 = x2 + h * i;
- s += 4 * SAMPLE(a1, a2);
+ double a1 = low + h * i;
+ s += 4 * SAMPLE(a1);
}
for (i = 2; i < N_SEGMENTS; i += 2)
{
- double a1 = x1 + h * i;
- double a2 = x2 + h * i;
- s += 2 * SAMPLE(a1, a2);
+ double a1 = low + h * i;
+ s += 2 * SAMPLE(a1);
}
- s += SAMPLE (x1 + width, x2 + width);
+ s += SAMPLE (high);
- return h * s * (1.0 / 3.0);
+ return h * s * (1.0 / 3.0) / scale;
}
}
@@ -238,65 +225,62 @@ create_1d_filter (int width,
int n_phases,
pixman_fixed_t *p)
{
- double step;
+ double step = 1.0 / n_phases;
+ double rwidth2 = filters[reconstruct].width / 2.0;
+ double swidth2 = scale * filters[sample].width / 2.0;
int i;
- step = 1.0 / n_phases;
-
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;
+ int x;
+ double pos, 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;
+ pos = ceil (frac - width / 2.0 - 0.5) + 0.5 - frac;
total = 0;
- for (x = x1; x < x2; ++x)
+ for (x = 0; x < width; ++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, scale, ilow - pos,
- ihigh - ilow);
- }
-
+ double ilow = MAX (pos - swidth2, -rwidth2);
+ double ihigh = MIN (pos + swidth2, rwidth2);
+ double c = integral (reconstruct,
+ sample, scale, pos,
+ ilow, ihigh);
total += c;
- *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5);
+ p[x] = (pixman_fixed_t)(c * 65536.0 + 0.5);
+ pos++;
}
/* Normalize */
- p -= width;
total = 1 / total;
new_total = 0;
- for (x = x1; x < x2; ++x)
+ for (x = 0; x < width; ++x)
{
- pixman_fixed_t t = (*p) * total + 0.5;
-
+ pixman_fixed_t t = p[x] * total + 0.5;
new_total += t;
- *p++ = t;
+ p[x] = t;
}
+ /* Distribute any remaining error over all samples */
if (new_total != pixman_fixed_1)
- *(p - width / 2) += (pixman_fixed_1 - new_total);
+ {
+ pixman_fixed_t delta = new_total - pixman_fixed_1;
+ pixman_fixed_t t = 0;
+ for (x = 0; x < width; ++x)
+ {
+ pixman_fixed_t new_t = delta * (x + 1) / width;
+ p[x] += new_t - t;
+ t = new_t;
+ }
+ }
+
+ p += width;
}
}
--
1.9.1
More information about the Pixman
mailing list