[Pixman] [PATCH v14 09/22] pixman-filter: Correct Simpsons integration

[spitzak at gmail.com]

From: Bill Spitzak <spitzak at gmail.com>

Simpsons uses cubic curve fitting, with 3 samples defining each cubic. This
makes the weights of the samples be in a pattern of 1,4,2,4,2...4,1, and then
dividing the result by 3.

The previous code was using weights of 1,2,0,6,0,6...,2,1.

With this fix the integration is accurate enough that the number of samples
could be reduced a lot. Multiples of 12 seem to work best.

v7: Merged with patch to reduce from 128 samples to 16
v9: Changed samples from 16 to 12
v10: Fixed rebase error that made it not compile
v11: minor whitespace change
v14: more whitespace changes

Signed-off-by: Bill Spitzak <spitzak at gmail.com>
Reviewed-by: Oded Gabbay <oded.gabbay at gmail.com>
---
 pixman/pixman-filter.c | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index c03a7f6..bd5ac4c 100644
--- a/pixman/pixman-filter.c
+++ b/pixman/pixman-filter.c
@@ -189,13 +189,19 @@ integral (pixman_kernel_t kernel1, double x1,
     }
     else
     {
-	/* Integration via Simpson's rule */
-#define N_SEGMENTS 128
+	/* Integration via Simpson's rule
+	 * See http://www.intmath.com/integration/6-simpsons-rule.php
+	 * 12 segments (6 cubic approximations) seems to produce best
+	 * result for lanczos3.linear, which was the combination that
+	 * showed the most errors.  This makes sense as the lanczos3
+	 * filter is 6 wide.
+	 */
+#define N_SEGMENTS 12
 #define SAMPLE(a1, a2)							\
 	(filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
 	
 	double s = 0.0;
-	double h = width / (double)N_SEGMENTS;
+	double h = width / N_SEGMENTS;
 	int i;
 
 	s = SAMPLE (x1, x2);
@@ -204,11 +210,14 @@ integral (pixman_kernel_t kernel1, double x1,
 	{
 	    double a1 = x1 + h * i;
 	    double a2 = x2 + h * i;
+	    s += 4 * SAMPLE (a1, a2);
+	}
 
+	for (i = 2; 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);
-- 
1.9.1