[Pixman] [PATCH 13/15] pixman-filter: Fix several issues related to normalization

[spitzak at gmail.com]

From: Søren Sandmann Pedersen <soren.sandmann at gmail.com>

There are a few bugs in the current normalization code

(1) The normalization is based on the sum of the *floating point*
    values generated by integral(). But in order to get the sum to be
    close to pixman_fixed_1, the sum of the rounded fixed point values
    should be used.

(2) The multiplications in the normalization loops often round the
    same way, so the residual error can fairly large.

(3) The residual error is added to the sample located at index
    (width - width / 2), which is not the midpoint for odd widths (and
    for width 1 is in fact outside the array).

This patch fixes these issues by (1) using the sum of the fixed point
values as the total to divide by, (2) doing error diffusion in the
normalization loop.

Changes in this version:

Keep the total as a pixman_fixed_t.

Merge the .5 offset used for rounding into the e value. This produces the same
results but matches how I have seen error diffusion done.

The possible error range at the end contains only one integer, zero. Floating
point rounding errors will not change this unless there are more than 2^20
samples in the filter. So there is never any possible way to add the residual
error, so the code for this is removed.

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

diff --git a/pixman/pixman-filter.c b/pixman/pixman-filter.c
index 11e7d0e..b3c36b8 100644
--- a/pixman/pixman-filter.c
+++ b/pixman/pixman-filter.c
@@ -245,9 +245,9 @@ create_1d_filter (int              width,
     for (i = 0; i < n_phases; ++i)
     {
         double frac = step / 2.0 + i * step;
-	pixman_fixed_t new_total;
+	pixman_fixed_t total;
         int x, x1, x2;
-	double total;
+	double e;
 
 	/* Sample convolution of reconstruction and sampling
 	 * filter. See rounding.txt regarding the rounding
@@ -278,24 +278,28 @@ create_1d_filter (int              width,
 			      ihigh - ilow);
 	    }
 
-	    total += c;
-            *p++ = (pixman_fixed_t)(c * 65536.0 + 0.5);
-        }
+	    *p = (pixman_fixed_t)floor (c * 65536.0 + 0.5);
+	    total += *p;
+	    p++;
+	}
 
-	/* Normalize */
+	/* Normalize, with error diffusion */
 	p -= width;
-        total = 1 / total;
-        new_total = 0;
+	e = 0.5;
 	for (x = x1; x < x2; ++x)
 	{
-	    pixman_fixed_t t = (*p) * total + 0.5;
-
-	    new_total += t;
+	    double v = (*p) * 65536.0 / total + e;
+	    pixman_fixed_t t = floor (v);
+	    e = v - t;
 	    *p++ = t;
 	}
-
-	if (new_total != pixman_fixed_1)
-	    *(p - width / 2) += (pixman_fixed_1 - new_total);
+	/* The sum of all the t values calculated above is 65536.5 - e
+	 * (where e is the final value of e). This can be seen because
+	 * total == sum(*p).  Since 0 <= e <= 1 the only integer in
+	 * range is 65536.  So there can be no remaining error. For a
+	 * few dozen entries in which can vary by no more than 2^16
+	 * double precision error does not effect this result.
+	 */
     }
 }
 
-- 
1.9.1