pixman-image: Implement PIXMAN_FILTER_GOOD/BEST as separable convolutions

In my opinion the low level code has to be in control of the filtering. This
is necessary to allow hardware implementations, to allow developement of new
filtering algorithms, and to allow optimization based on knowledge of the
exact filter being used.

This implements the same GOOD/BEST as Cairo, with minor improvements.

The GOOD will produce exactly the same BILINEAR result as before for any scale
greater than 3/4, or for scale of 1/2 with no rotation and an integer
translation. This means the output is unchanged for most current users of
GOOD.

GOOD uses:

 scale < 1/16 : BOX.BOX at size 16
 scale < 3/4 : BOX.BOX at size 1/scale
 larger : BOX.BOX at size 1

 If both directions have a scale >= 3/4 or a scale of 1/2 and an integer
 translation, the faster PIXMAN_FILTER_BILINEAR code is used. This is
 compatable at these scales with older versions of pixman where bilinear
 was always used for GOOD.

BEST uses:

 scale < 1/24 : BOX.BOX at size 24
 scale < 1/16 : BOX.BOX at size 1/scale
 scale < 1 : IMPULSE.LANCZOS2 at size 1/scale
 scale < 2.333 : IMPULSE.LANCZOS2 at size 1
 scale < 128 : BOX.LANCZOS2 at size 1/(scale-1) (antialiased square pixels)
 larger : BOX.LANCZOS2 at size 1/127 (antialias blur gets thicker)

v8: Cutoff in BEST between IMPULSE.LANCZOS2 and BOX.LANCZOS2 adjusted for
    a better match between the filters.
v9: Fixed divide-by-zero from all-zero matrix found by stress-test
v11: Whitespace and formatting fixes
     Moved demo changes to a later patch
v12: Whitespace and formatting fixes
v14: Compute subsample bits here
     Works for non-affine and when fast paths are disabled
     Fixed big memory leak of array for BEST

Signed-off-by: Bill Spitzak <spitzak@gmail.com>

3 files changed, 206 insertions, 7 deletions

diff --git a/pixman/pixman-bits-image.c b/pixman/pixman-bits-image.c
index dcdcc699..8c25c5e4 100644
--- a/pixman/pixman-bits-image.c
+++ b/pixman/pixman-bits-image.c
@@ -320,21 +320,20 @@ bits_image_fetch_pixel_filtered (bits_image_t *image,
     case PIXMAN_FILTER_NEAREST:
     case PIXMAN_FILTER_FAST:
 	return bits_image_fetch_pixel_nearest (image, x, y, get_pixel);
-	break;
 
-    case PIXMAN_FILTER_BILINEAR:
     case PIXMAN_FILTER_GOOD:
     case PIXMAN_FILTER_BEST:
+	if (image->common.flags & FAST_PATH_SEPARABLE_CONVOLUTION_FILTER)
+            return bits_image_fetch_pixel_separable_convolution (image, x, y, get_pixel);
+	/* else fall through to BILINEAR */
+    case PIXMAN_FILTER_BILINEAR:
 	return bits_image_fetch_pixel_bilinear (image, x, y, get_pixel);
-	break;
 
     case PIXMAN_FILTER_CONVOLUTION:
 	return bits_image_fetch_pixel_convolution (image, x, y, get_pixel);
-	break;
 
     case PIXMAN_FILTER_SEPARABLE_CONVOLUTION:
         return bits_image_fetch_pixel_separable_convolution (image, x, y, get_pixel);
-        break;
 
     default:
         break;
diff --git a/pixman/pixman-image.c b/pixman/pixman-image.c
index 8ad28912..c6aabfb9 100644
--- a/pixman/pixman-image.c
+++ b/pixman/pixman-image.c
@@ -28,6 +28,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <assert.h>
+#include <math.h>
 
 #include "pixman-private.h"
 
@@ -261,13 +262,34 @@ pixman_disable_out_of_bounds_workaround (void)
 {
 }
 
+/* Compute subsample bits for a given size so there are approximately
+ * the same number of samples in each filter size.
+ */
+static int
+subsample_bits(double size)
+{
+    int desired_samples = 14;
+    if (size >= 1.0)
+	desired_samples *= 2.0 / (size + 1.0);
+    else
+	desired_samples *= (size + 1.0) / (2.0 * size);
+    if (desired_samples <= 1.0)
+	return 0;
+    else if (desired_samples >= 256.0)
+	return 8;
+    else
+	return (int) ceil (log2(desired_samples) - .01);
+}
+
 static void
 compute_image_info (pixman_image_t *image)
 {
     pixman_format_code_t code;
     uint32_t flags = 0;
     int nearest_ok = FALSE;
+    int bilinear_ok = FALSE;
     pixman_fixed_t (*m)[3];
+    double dx, dy;
 
     /* Transform */
     if (!image->common.transform)
@@ -354,14 +376,182 @@ compute_image_info (pixman_image_t *image)
 	break;
 
     case PIXMAN_FILTER_BILINEAR:
-    case PIXMAN_FILTER_GOOD:
-    case PIXMAN_FILTER_BEST:
 	if (nearest_ok)
 	    flags |= (FAST_PATH_NEAREST_FILTER | FAST_PATH_NO_CONVOLUTION_FILTER);
 	else
 	    flags |= (FAST_PATH_BILINEAR_FILTER | FAST_PATH_NO_CONVOLUTION_FILTER);
 	break;
 
+    case PIXMAN_FILTER_GOOD:
+	if (nearest_ok)
+	{
+	    flags |= (FAST_PATH_NEAREST_FILTER | FAST_PATH_NO_CONVOLUTION_FILTER);
+	    break;
+	}
+
+	/* Compute filter sizes. This is the bounding box of a
+	 * diameter=1 circle transformed by the matrix. Scaling
+	 * down produces values greater than 1. See comment in
+	 * ../demos/scale.c for proof hypot is correct.
+	 *
+	 * For non-affine the circle is centered on one of the 4
+	 * points 1,1 away from the origin. Which one depends on
+	 * the signs of the values in the last row of the matrix,
+	 * chosen to avoid dividing by zero.
+	 */
+	/* This division factor both accounts for the w component
+	 * and converts from fixed to float.
+	 */
+	dy = abs(m[2][0]) + abs(m[2][1]) + abs(m[2][2]);
+	if (dy)
+	    dy = 1.0 / dy;
+	/* There are some signs that hypot is faster with numbers near 1
+	 * so the division is done first. Mathematically it should work
+	 * to divide afterwards.
+	 */
+	dx = hypot (m[0][0] * dy, m[0][1] * dy);
+	dy = hypot (m[1][0] * dy, m[1][1] * dy);
+
+	/* scale < 1/16 : BOX.BOX at size 16
+	 * scale < 3/4 : BOX.BOX at size 1/scale
+	 * larger : BOX.BOX at size 1
+	 *
+	 * If both directions have a scale >= 3/4 or a scale of
+	 * 1/2 and an integer translation, the faster
+	 * PIXMAN_FILTER_BILINEAR code is used.
+	 *
+	 * Filter size is clamped to 16 to prevent extreme slowness.
+	 */
+	if (dx <= 4.0 / 3)
+	{
+	    dx = 1.0;
+	    bilinear_ok = TRUE;
+	}
+	else if (dx > 16.0)
+	{
+	    dx = 16.0;
+	}
+	else if (dx > 1.999 && dx < 2.001 &&
+		 abs(m[0][0] * m[0][1]) < 4 &&
+		 abs(pixman_fixed_frac(m[0][2]) < 2))
+	{
+	    bilinear_ok = TRUE;
+	}
+
+	if (dy <= 4.0 / 3)
+	{
+	    dy = 1.0;
+	}
+	else if (dy > 16.0)
+	{
+	    dy = 16.0;
+	    bilinear_ok = FALSE;
+	}
+	else if (bilinear_ok)
+	{
+	    bilinear_ok =
+		(dy > 1.999 && dy < 2.001 &&
+		 abs(m[1][0] * m[1][1]) < 4 &&
+		 abs(pixman_fixed_frac(m[1][2]) < 2));
+	}
+
+	if (bilinear_ok)
+	{
+	    flags |= (FAST_PATH_BILINEAR_FILTER |
+		      FAST_PATH_NO_CONVOLUTION_FILTER);
+	    break;
+	}
+
+	if (image->common.filter_params)
+	    free (image->common.filter_params);
+
+	image->common.filter_params =
+	    pixman_filter_create_separable_convolution
+	    ( & image->common.n_filter_params,
+	      pixman_double_to_fixed(dx),
+	      pixman_double_to_fixed(dy),
+	      PIXMAN_KERNEL_BOX,
+	      PIXMAN_KERNEL_BOX,
+	      PIXMAN_KERNEL_BOX,
+	      PIXMAN_KERNEL_BOX,
+	      subsample_bits(dx), subsample_bits(dy));
+
+	flags |= FAST_PATH_SEPARABLE_CONVOLUTION_FILTER;
+	break;
+
+    case PIXMAN_FILTER_BEST:
+	if (nearest_ok)
+	{
+	    flags |= (FAST_PATH_NEAREST_FILTER | FAST_PATH_NO_CONVOLUTION_FILTER);
+	    break;
+	}
+	/* See notes above about filter sizes */
+	dy = abs(m[2][0]) + abs(m[2][1]) + abs(m[2][2]);
+	if (dy)
+	    dy = 1.0 / dy;
+	dx = hypot (m[0][0] * dy, m[0][1] * dy);
+	dy = hypot (m[1][0] * dy, m[1][1] * dy);
+
+	/* scale < 1/24 : BOX.BOX at size 24
+	 * scale < 1/16 : BOX.BOX at size 1/scale
+	 * scale < 1 : IMPULSE.LANCZOS2 at size 1/scale
+	 * scale < 2.333 : IMPULSE.LANCZOS2 at size 1
+	 * scale < 128 : BOX.LANCZOS2 at size 1/(scale-1)
+	 * larger : BOX.LANCZOS2 at size 1/127
+	 *
+	 * Filter switches to box and then clamps at 24 to prevent
+	 * extreme slowness.
+	 *
+	 * When enlarging this produces square pixels with an
+	 * anti-aliased border between them.  At scales larger
+	 * than 128x the antialias blur is increased to avoid
+	 * making lots of subsamples.
+	 */
+	if (dx > 24.0)
+	{
+	    dx = 24.0;
+	}
+	else if (dx < 1.0)
+	{
+	    if (dx >= 3.0/7)
+		dx = 1.0;
+	    else if (dx > 1.0/128)
+		dx /= 1.0 - dx;
+	    else
+		dx = 1.0/127;
+	}
+
+	if (dy > 24.0)
+	{
+	    dy = 24.0;
+	}
+	else if (dy < 1.0)
+	{
+	    if (dy >= 3.0/7)
+		dy = 1.0;
+	    else if (dy > 1.0/128)
+		dy /= 1.0 - dy;
+	    else
+		dy = 1.0/127;
+	}
+
+	if (image->common.filter_params)
+	    free (image->common.filter_params);
+
+	image->common.filter_params =
+	    pixman_filter_create_separable_convolution
+	    ( & image->common.n_filter_params,
+	      pixman_double_to_fixed(dx),
+	      pixman_double_to_fixed(dy),
+	      dx >= 1.0 && dx < 16.0 ? PIXMAN_KERNEL_IMPULSE : PIXMAN_KERNEL_BOX,
+	      dy >= 1.0 && dy < 16.0 ? PIXMAN_KERNEL_IMPULSE : PIXMAN_KERNEL_BOX,
+	      dx < 16.0 ? PIXMAN_KERNEL_LANCZOS2 : PIXMAN_KERNEL_BOX,
+	      dy < 16.0 ? PIXMAN_KERNEL_LANCZOS2 : PIXMAN_KERNEL_BOX,
+	      subsample_bits(dx), subsample_bits(dy));
+
+	flags |= FAST_PATH_SEPARABLE_CONVOLUTION_FILTER;
+	break;
+
     case PIXMAN_FILTER_CONVOLUTION:
 	break;
 
diff --git a/pixman/pixman.c b/pixman/pixman.c
index f932eac3..1dce700d 100644
--- a/pixman/pixman.c
+++ b/pixman/pixman.c
@@ -467,6 +467,16 @@ analyze_extent (pixman_image_t       *image,
 	    
 	case PIXMAN_FILTER_GOOD:
 	case PIXMAN_FILTER_BEST:
+	    if (image->common.flags & FAST_PATH_SEPARABLE_CONVOLUTION_FILTER)
+	    {
+	        params = image->common.filter_params;
+	        x_off = - pixman_fixed_e - ((params[0] - pixman_fixed_1) >> 1);
+	        y_off = - pixman_fixed_e - ((params[1] - pixman_fixed_1) >> 1);
+	        width = params[0];
+	        height = params[1];
+	        break;
+	    }
+	    /* else fall through to BILINEAR */
 	case PIXMAN_FILTER_BILINEAR:
 	    x_off = - pixman_fixed_1 / 2;
 	    y_off = - pixman_fixed_1 / 2;