pixman-filter: Corrections to the integral() function

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).

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.

v7: This is a merge of 4 patches and lots of new code cleanup and fixes
 determined by examining the gnuplot output

v9: Restored the recursion splitting at zero for linear filter

v10: Small change from here moved to previous Simpsons patch so it compiles
     Merged patch to get correct subsample positions when subsample_bits==0

v11: Whitespace fixes

Signed-off-by: Bill Spitzak <spitzak@gmail.com>
Reviewed-by: Oded Gabbay <oded.gabbay@gmail.com>
Acked-by: Oded Gabbay <oded.gabbay@gmail.com>

1 files changed, 79 insertions, 75 deletions

diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index 8b8fb82..e82871f 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,51 @@ static const filter_info_t filters[] =
     { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
 };
 
-/* This function scales the @sample filter by @size, then
- * aligns @x1 in @reconstruct with @x2 in @sample and
- * and integrates the product of the kernels across @width.
+/* This function scales the @sample filter by @size, 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 size, double x2,
-	  double width)
+integral (pixman_kernel_t reconstruct,
+	  pixman_kernel_t sample, double size, double pos,
+	  double low, double high)
 {
+    if (high < low)
+    {
+	return 0.0;
+    }
+    else if (sample == PIXMAN_KERNEL_IMPULSE)
+    {
+	return filters[reconstruct].func (-pos);
+    }
+    else if (reconstruct == PIXMAN_KERNEL_IMPULSE)
+    {
+	return filters[sample].func (-pos / size) / size;
+    }
+    else if (reconstruct == PIXMAN_KERNEL_BOX && sample == PIXMAN_KERNEL_BOX)
+    {
+	assert (high <= low + 1.0);
+	return (high - low) / size;
+    }
     /* 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)
+    else if (reconstruct == PIXMAN_KERNEL_LINEAR && low < 0 && high > 0)
     {
 	return
-	    integral (reconstruct, x1, sample, size, x2, - x1) +
-	    integral (reconstruct, 0, sample, size, x2 - x1, width + x1);
+	    integral (reconstruct, sample, size, pos, low, 0) +
+	    integral (reconstruct, sample, size, pos, 0, high);
     }
-    else if (x2 < 0 && x2 + width > 0)
+    else if (sample == PIXMAN_KERNEL_LINEAR && low < pos && high > pos)
     {
 	return
-	    integral (reconstruct, x1, sample, size, x2, - x2) +
-	    integral (reconstruct, x1 - x2, sample, size, 0, width + x2);
-    }
-    else if (reconstruct == PIXMAN_KERNEL_IMPULSE)
-    {
-	assert (width == 0.0);
-	return filters[sample].func (x2 / size);
-    }
-    else if (sample == PIXMAN_KERNEL_IMPULSE)
-    {
-	assert (width == 0.0);
-	return filters[reconstruct].func (x1);
+	    integral (reconstruct, sample, size, pos, low, pos) +
+	    integral (reconstruct, sample, size, pos, pos, high);
     }
     else
     {
@@ -197,32 +203,30 @@ integral (pixman_kernel_t reconstruct, double x1,
 	 * filter is 6 wide.
 	 */
 #define N_SEGMENTS 12
-#define SAMPLE(a1, a2)							\
-	(filters[reconstruct].func ((a1)) * filters[sample].func ((a2) / size))
-	
+#define SAMPLE(a)							\
+	(filters[reconstruct].func ((a)) * filters[sample].func (((a) - pos) / size))
+
 	double s = 0.0;
-	double h = width / N_SEGMENTS;
+	double h = (high - low) / 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) / size;
     }
 }
 
@@ -234,65 +238,65 @@ 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 = size * 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;
+	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;
+	if (n_phases & 1)
+	    pos = frac - width / 2.0;
+	else
+	    pos = ceil (frac - width / 2.0 - 0.5) + 0.5 - frac;
 
 	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 - size * filters[sample].width / 2.0;
-	    double shigh = slow + size * 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, size, ilow - pos,
-			      ihigh - ilow);
-	    }
-
+	for (x = 0; x < width; ++x)
+	{
+	    double ilow = MAX (pos - swidth2, -rwidth2);
+	    double ihigh = MIN (pos + swidth2, rwidth2);
+	    double c = integral (reconstruct,
+				 sample, size, 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)
+	total = 1 / total;
+	new_total = 0;
+	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;
     }
 }