[Pixman] [PATCH v14 10/22] pixman-filter: Correct integration with impulse filters

[spitzak at gmail.com]

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.

Signed-off-by: Bill Spitzak <spitzak at gmail.com>
---
 pixman/pixman-filter.c | 111 ++++++++++++++++++++++---------------------------
 1 file changed, 49 insertions(+), 62 deletions(-)

diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index bd5ac4c..fef2189 100644
--- a/pixman/pixman-filter.c
+++ b/pixman/pixman-filter.c
@@ -147,45 +147,46 @@ static const filter_info_t filters[] =
     { 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 scales @kernel2 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 kernel1, double x1,
-	  pixman_kernel_t kernel2, double scale, double x2,
-	  double width)
+integral (pixman_kernel_t kernel1,
+	  pixman_kernel_t kernel2, double scale, double pos,
+	  double low, double high)
 {
+    if (high < low)
+    {
+	return 0.0;
+    }
+    else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
+    {
+	return filters[kernel1].func (-pos);
+    }
+    else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
+    {
+	return filters[kernel2].func (-pos * scale);
+    }
     /* 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 (low < 0 && high > 0)
     {
 	return
-	    integral (kernel1, x1, kernel2, scale, x2, - x1) +
-	    integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
+	    integral (kernel1, kernel2, scale, pos, low, 0) +
+	    integral (kernel1, kernel2, scale, pos, 0, high);
     }
-    else if (x2 < 0 && x2 + width > 0)
+    else if (low < pos && high > pos)
     {
 	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);
+	    integral (kernel1, kernel2, scale, pos, low, pos) +
+	    integral (kernel1, kernel2, scale, pos, pos, high);
     }
     else
     {
@@ -197,30 +198,28 @@ integral (pixman_kernel_t kernel1, double x1,
 	 * filter is 6 wide.
 	 */
 #define N_SEGMENTS 12
-#define SAMPLE(a1, a2)							\
-	(filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
-	
+#define SAMPLE(a)							\
+	(filters[kernel1].func ((a)) * filters[kernel2].func (((a) - pos)))
+
 	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);
     }
@@ -235,55 +234,43 @@ create_1d_filter (int              width,
 		  pixman_fixed_t *p)
 {
     double step;
+    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;
+	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, 1.0 / size, ilow - pos,
-			      ihigh - ilow);
-	    }
-
+	for (x = 0; x < width; ++x, ++pos)
+	{
+	    double ilow = MAX (pos - swidth2, -rwidth2);
+	    double ihigh = MIN (pos + swidth2, rwidth2);
+	    double c = integral (reconstruct,
+				 sample, 1.0 / size, pos,
+				 ilow, ihigh);
 	    total += c;
-            *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5);
-        }
+	    *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5);
+	}
 
 	/* 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;
 
-- 
1.9.1