pixman-filter: Correct Simpsons integration

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@gmail.com>
Reviewed-by: Oded Gabbay <oded.gabbay@gmail.com>
Reviewed-by: Søren Sandmann <soren.sandmann@gmail.com>

diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index c018a67637f559ade6446f6924abf958c251dfca..ecff852b5c5ff073f4a230f75205f10ba6c77f36 100644 (file)
--- 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);