Exploit the fact that all our filters are separable to simplify the

Tue Apr  1 15:33:51 2003  Owen Taylor  <otaylor@redhat.com>

        * pixops/pixops.c (make_weights): Exploit the fact that all
        our filters are separable to simplify the calculation of
        the weight tables. (Based on a patch from Brian Cameron.)

        * pixops/Makefile.am pixbuf-transform-math.ltx:
        Add a short article describing how the math in pixops.c
        works.

3 files changed, 395 insertions, 346 deletions

diff --git a/gdk-pixbuf/pixops/Makefile.am b/gdk-pixbuf/pixops/Makefile.am
index 16bc8d4b7b..f666d9b0b6 100644
--- a/gdk-pixbuf/pixops/Makefile.am
+++ b/gdk-pixbuf/pixops/Makefile.am
@@ -29,4 +29,5 @@ libpixops_la_SOURCES =  		\
 
 EXTRA_DIST =				\
 	DETAILS				\
+	pixbuf-transform-math.ltx	\
 	makefile.msc
diff --git a/gdk-pixbuf/pixops/pixbuf-transform-math.ltx b/gdk-pixbuf/pixops/pixbuf-transform-math.ltx
new file mode 100644
index 0000000000..19e231308d
--- /dev/null
+++ b/gdk-pixbuf/pixops/pixbuf-transform-math.ltx
@@ -0,0 +1,112 @@
+\documentclass{article}
+
+\begin{document}
+
+\title{Some image transform math}
+\author{Owen Taylor}
+\date{18 February 2003}
+\maketitle
+
+\section{Basics}
+
+The transform process is composed of three steps;
+first we reconstruct a continuous image from the 
+source data \(A_{i,j}\):
+\[a(u,v) = \sum_{i = -\infty}^{\infty} \sum_{j = -\infty}^{\infty} A_{i,j}F\left( {u - i \atop v - j} \right) \]
+Then we transform from destination coordinates to source coordinates:
+\[b(x,y) = a\left(u(x,y) \atop v(x,y)\right)
+         = a\left(t_{00}x + t_{01}y + t_{02} \atop t_{10}x + t_{11}y + t_{12} \right)\]
+Finally, we resample using a sampling function \(G\):
+\[B_{x_0,y_0} = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} b(x,y)G\left( {x - x_0 \atop y - y_0} \right) dxdy\]
+Putting all of these together:
+\[B_{x_0,y_0} = 
+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
+\sum_{i = -\infty}^{\infty} \sum_{j = -\infty}^{\infty} A_{i,j}
+F\left( {u(x,y) - i \atop v(x,y) - j} \right)
+G\left( {x - x_0 \atop y - y_0} \right) dxdy\]
+We can reverse the order of the integrals and the sums:
+\[B_{x_0,y_0} = 
+\sum_{i = -\infty}^{\infty} \sum_{j = -\infty}^{\infty} A_{i,j}
+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
+F\left( {u(x,y) - i \atop v(x,y) - j} \right)
+G\left( {x - x_0 \atop y - y_0} \right) dxdy\]
+Which shows that the destination pixel values are a linear combination of the 
+source pixel values. But the coefficents depend on \(x_0\) and \(y_0\). 
+To simplify this a bit, define:
+\[i_0 = \lfloor u(x_0,y_0) \rfloor = \lfloor {t_{00}x_0 + t_{01}y_0 + t_{02}} \rfloor \]
+\[j_0 = \lfloor v(x_0,y_0) \rfloor = \lfloor {t_{10}x_0 + t_{11}y_0 + t_{12}} \rfloor \]
+\[\Delta_u = u(x_0,y_0) - i_0 = t_{00}x_0 + t_{01}y_0 + t_{02} - \lfloor {t_{00}x_0 + t_{01}y_0 + t_{02}} \rfloor \]
+\[\Delta_v = v(x_0,y_0) - j_0 = t_{10}x_0 + t_{11}y_0 + t_{12} - \lfloor {t_{10}x_0 + t_{11}y_0 + t_{12}} \rfloor \]
+Then making the transforms \(x' = x - x_0\), \(y' = y - x_0\), \(i' = i - i_0\), \(j' = j - x_0\)
+\begin{eqnarray*}
+F(u,v) & = & F\left( {t_{00}x + t_{01}y + t_{02} - i \atop t_{10}x + t_{11}y + t_{12} - j} \right)\\
+       & = & F\left( {t_{00}(x'+x_0) + t_{01}(y'+y_0) + t_{02} - (i'+i_0) \atop 
+                      t_{10}(x'+x_0) + t_{11}(y'+y_0) + t_{12} - (j'+j_0)} \right) \\
+       & = & F\left( {\Delta_u + t_{00}x' + t_{01}y' - i' \atop 
+                      \Delta_v + t_{10}x' + t_{11}y' - j'} \right)
+\end{eqnarray*}
+Using that, we can then reparameterize the sums and integrals and
+define coefficients that depend only on \((\Delta_u,\Delta_v)\),
+which we'll call the \emph{phase} at the point \((x_0,y_0)\):
+\[
+B_{x_0,y_0} = 
+\sum_{i = -\infty}^{\infty} \sum_{j = -\infty}^{\infty} A_{i_0+i,j_0+j} C_{i,j}(\Delta_u,\Delta_v)
+\]
+\[
+C_{i,j}(\Delta_u,\Delta_v) =
+\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
+F\left( {\Delta_u + t_{00}x + t_{01}y - i \atop 
+         \Delta_v + t_{10}x + t_{11}y - j} \right)
+G\left( {x \atop y} \right) dxdy
+\]
+\section{Separability}
+A frequent special case is when the reconstruction and sampling functions
+are of the form:
+\[F(u,v) = f(u)f(v)\]
+\[G(x,y) = g(x)g(y)\]
+If we also have a transform that is purely a scale and translation;
+(\(t_{10} = 0\), \(t_{01} = 0\)), then we can separate 
+\(C_{i,j}(\Delta_u,\Delta_v)\) into the product of a \(x\) portion
+and a \(y\) portion:
+\[C_{i,j}(\Delta_u,\Delta_v) = c_{i}(\Delta_u) c_{j}(\Delta_v)\]
+\[c_{i}(\Delta_u) = \int_{-\infty}^{\infty} f(\Delta_u + t_{00}x - i)g(x)dx\]
+\[c_{j}(\Delta_v) = \int_{-\infty}^{\infty} f(\Delta_v + t_{11}y - j)g(y)dy\]
+
+\section{Some filters}
+gdk-pixbuf provides 4 standard filters for scaling, under the names ``NEAREST'',
+``TILES'', ``BILINEAR'', and ``HYPER''. All of turn out to be separable
+as discussed in the previous section. 
+For ``NEAREST'' filter, the reconstruction function is simple replication
+and the sampling function is a delta function\footnote{A delta function is an infinitely narrow spike, such that:
+\[\int_{-\infty}^{\infty}\delta(x)f(x) = f(0)\]}:
+\[f(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                0, & otherwise}\]
+\[g(t) = \delta(t - 0.5)\]
+For ``TILES'', the reconstruction function is again replication, but we
+replace the delta-function for sampling with a box filter:
+\[f(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                0, & otherwise}\]
+\[g(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                0, & otherwise}\]
+The ``HYPER'' filter (in practice, it was originally intended to be 
+something else) uses bilinear interpolation for reconstruction and
+a box filter for sampling:
+\[f(t) = \cases{1 - |t - 0.5|, & if \(-0.5 \le t \le 1.5\); \cr 
+                0, & otherwise}\]
+\[g(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                0, & otherwise}\]
+The ``BILINEAR'' filter is defined in a somewhat more complicated way;
+the definition depends on the scale factor in the transform (\(t_{00}\)
+or \(t_{01}]\). In the \(x\) direction, for \(t_{00} < 1\), it is
+the same as for ``TILES'':
+\[f_u(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                  0, & otherwise}\]
+\[g_u(t) = \cases{1, & if \(0 \le t \le 1\); \cr 
+                  0, & otherwise}\]
+but for \(t_{10} > 1\), we use bilinear reconstruction and delta-function
+sampling:
+\[f_u(t) = \cases{1 - |t - 0.5|, & if \(-0.5 \le t \le 1.5\); \cr 
+                  0, & otherwise}\]
+\[g_u(t) = \delta(t - 0.5)\]
+The behavior in the \(y\) direction depends in the same way on \(t_{11}\).
+\end{document}
\ No newline at end of file
diff --git a/gdk-pixbuf/pixops/pixops.c b/gdk-pixbuf/pixops/pixops.c
index 1a6c06118f..04a429d296 100644
--- a/gdk-pixbuf/pixops/pixops.c
+++ b/gdk-pixbuf/pixops/pixops.c
@@ -11,14 +11,20 @@
 #define SCALE_SHIFT 16
 
 typedef struct _PixopsFilter PixopsFilter;
+typedef struct _PixopsFilterDimension PixopsFilterDimension;
+
+struct _PixopsFilterDimension
+{
+  int n;
+  double offset;
+  double *weights;
+};
 
 struct _PixopsFilter
 {
-  int *weights;
-  int n_x;
-  int n_y;
-  double x_offset;
-  double y_offset;
+  PixopsFilterDimension x;
+  PixopsFilterDimension y;
+  double overall_alpha;
 }; 
 
 typedef guchar *(*PixopsLineFunc) (int *weights, int n_x, int n_y,
@@ -938,6 +944,63 @@ process_pixel (int *weights, int n_x, int n_y,
   (*pixel_func) (dest, dest_x, dest_channels, dest_has_alpha, src_has_alpha, check_size, color1, color2, r, g, b, a);
 }
 
+static void 
+correct_total (int    *weights, 
+               int    n_x, 
+               int    n_y,
+               int    total, 
+               double overall_alpha)
+{
+  int correction = (int)(0.5 + 65536 * overall_alpha) - total;
+
+  if (correction != 0)
+    {
+      int i;
+      for (i = n_x * n_y - 1; i >= 0; i--) 
+        {
+          if (*(weights + i) + correction >= 0) 
+            {
+              *(weights + i) += correction;
+              break;
+            }
+        }
+    }
+}
+
+static int *
+make_filter_table (PixopsFilter *filter)
+{
+  int i_offset, j_offset;
+  int n_x = filter->x.n;
+  int n_y = filter->y.n;
+  int *weights = g_new (int, SUBSAMPLE * SUBSAMPLE * n_x * n_y);
+
+  for (i_offset=0; i_offset < SUBSAMPLE; i_offset++)
+    for (j_offset=0; j_offset < SUBSAMPLE; j_offset++)
+      {
+        double weight;
+        int *pixel_weights = weights + ((i_offset*SUBSAMPLE) + j_offset) * n_x * n_y;
+        int total = 0;
+        int i, j;
+
+        for (i=0; i < n_y; i++)
+          for (j=0; j < n_x; j++)
+            {
+              weight = filter->x.weights[(j_offset * n_x) + j] *
+                       filter->y.weights[(i_offset * n_y) + i] *
+                       filter->overall_alpha * 65536 + 0.5;
+
+              total += (int)weight;
+
+              *(pixel_weights + n_x * i + j) = weight;
+            }
+
+        correct_total (pixel_weights, n_x, n_y, total, filter->overall_alpha);
+      }
+
+  return weights;
+}
+
 static void
 pixops_process (guchar         *dest_buf,
 		int             render_x0,
@@ -966,38 +1029,42 @@ pixops_process (guchar         *dest_buf,
 {
   int i, j;
   int x, y;			/* X and Y position in source (fixed_point) */
-  guchar **line_bufs = g_new (guchar *, filter->n_y);
+  
+  guchar **line_bufs = g_new (guchar *, filter->y.n);
+  int *filter_weights = make_filter_table (filter);
 
   int x_step = (1 << SCALE_SHIFT) / scale_x; /* X step in source (fixed point) */
   int y_step = (1 << SCALE_SHIFT) / scale_y; /* Y step in source (fixed point) */
 
   int check_shift = check_size ? get_check_shift (check_size) : 0;
 
-  int scaled_x_offset = floor (filter->x_offset * (1 << SCALE_SHIFT));
+  int scaled_x_offset = floor (filter->x.offset * (1 << SCALE_SHIFT));
 
   /* Compute the index where we run off the end of the source buffer. The furthest
    * source pixel we access at index i is:
    *
-   *  ((render_x0 + i) * x_step + scaled_x_offset) >> SCALE_SHIFT + filter->n_x - 1
+   *  ((render_x0 + i) * x_step + scaled_x_offset) >> SCALE_SHIFT + filter->x.n - 1
    *
    * So, run_end_index is the smallest i for which this pixel is src_width, i.e, for which:
    *
-   *  (i + render_x0) * x_step >= ((src_width - filter->n_x + 1) << SCALE_SHIFT) - scaled_x_offset
+   *  (i + render_x0) * x_step >= ((src_width - filter->x.n + 1) << SCALE_SHIFT) - scaled_x_offset
    *
    */
 #define MYDIV(a,b) ((a) > 0 ? (a) / (b) : ((a) - (b) + 1) / (b))    /* Division so that -1/5 = -1 */
   
-  int run_end_x = (((src_width - filter->n_x + 1) << SCALE_SHIFT) - scaled_x_offset);
+  int run_end_x = (((src_width - filter->x.n + 1) << SCALE_SHIFT) - scaled_x_offset);
   int run_end_index = MYDIV (run_end_x + x_step - 1, x_step) - render_x0;
   run_end_index = MIN (run_end_index, render_x1 - render_x0);
 
-  y = render_y0 * y_step + floor (filter->y_offset * (1 << SCALE_SHIFT));
+  y = render_y0 * y_step + floor (filter->y.offset * (1 << SCALE_SHIFT));
   for (i = 0; i < (render_y1 - render_y0); i++)
     {
       int dest_x;
       int y_start = y >> SCALE_SHIFT;
       int x_start;
-      int *run_weights = filter->weights + ((y >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) * filter->n_x * filter->n_y * SUBSAMPLE;
+      int *run_weights = filter_weights +
+                         ((y >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) *
+                         filter->x.n * filter->y.n * SUBSAMPLE;
       guchar *new_outbuf;
       guint32 tcolor1, tcolor2;
       
@@ -1015,7 +1082,7 @@ pixops_process (guchar         *dest_buf,
 	  tcolor2 = color2;
 	}
 
-      for (j=0; j<filter->n_y; j++)
+      for (j=0; j<filter->y.n; j++)
 	{
 	  if (y_start <  0)
 	    line_bufs[j] = (guchar *)src_buf;
@@ -1033,7 +1100,7 @@ pixops_process (guchar         *dest_buf,
 
       while (x_start < 0 && outbuf < outbuf_end)
 	{
-	  process_pixel (run_weights + ((x >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) * (filter->n_x * filter->n_y), filter->n_x, filter->n_y,
+	  process_pixel (run_weights + ((x >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) * (filter->x.n * filter->y.n), filter->x.n, filter->y.n,
 			 outbuf, dest_x, dest_channels, dest_has_alpha,
 			 line_bufs, src_channels, src_has_alpha,
 			 x >> SCALE_SHIFT, src_width,
@@ -1045,7 +1112,7 @@ pixops_process (guchar         *dest_buf,
 	  outbuf += dest_channels;
 	}
 
-      new_outbuf = (*line_func) (run_weights, filter->n_x, filter->n_y,
+      new_outbuf = (*line_func) (run_weights, filter->x.n, filter->y.n,
 				 outbuf, dest_x,
 				 dest_buf + dest_rowstride * i + run_end_index * dest_channels,
 				 dest_channels, dest_has_alpha,
@@ -1059,7 +1126,7 @@ pixops_process (guchar         *dest_buf,
 
       while (outbuf < outbuf_end)
 	{
-	  process_pixel (run_weights + ((x >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) * (filter->n_x * filter->n_y), filter->n_x, filter->n_y,
+	  process_pixel (run_weights + ((x >> (SCALE_SHIFT - SUBSAMPLE_BITS)) & SUBSAMPLE_MASK) * (filter->x.n * filter->y.n), filter->x.n, filter->y.n,
 			 outbuf, dest_x, dest_channels, dest_has_alpha,
 			 line_bufs, src_channels, src_has_alpha,
 			 x >> SCALE_SHIFT, src_width,
@@ -1074,320 +1141,216 @@ pixops_process (guchar         *dest_buf,
     }
 
   g_free (line_bufs);
+  g_free (filter_weights);
 }
 
-static void 
-correct_total (int    *weights, 
-	       int    n_x, 
-	       int    n_y,
-	       int    total, 
-	       double overall_alpha)
-{
-  int correction = (int)(0.5 + 65536 * overall_alpha) - total;
-  int i;
-  for (i = n_x * n_y - 1; i >= 0; i--) 
-    {
-      if (*(weights + i) + correction >= 0) 
-	{
-	  *(weights + i) += correction;
-	  break;
-	}
-    }
-}  
-
+/* Compute weights for reconstruction by replication followed by
+ * sampling with a box filter
+ */
 static void
-tile_make_weights (PixopsFilter *filter, double x_scale, double y_scale, double overall_alpha)
+tile_make_weights (PixopsFilterDimension *dim,
+		   double                 scale)
 {
-  int i_offset, j_offset;
+  int n = ceil (1 / scale + 1);
+  double *pixel_weights = g_new (double, SUBSAMPLE * n);
+  int offset;
+  int i;
 
-  int n_x = ceil(1/x_scale + 1);
-  int n_y = ceil(1/y_scale + 1);
+  dim->n = n;
+  dim->offset = 0;
+  dim->weights = pixel_weights;
 
-  filter->x_offset = 0;
-  filter->y_offset = 0;
-  filter->n_x = n_x;
-  filter->n_y = n_y;
-  filter->weights = g_new (int, SUBSAMPLE * SUBSAMPLE * n_x * n_y);
+  for (offset = 0; offset < SUBSAMPLE; offset++)
+    {
+      double x = (double)offset / SUBSAMPLE;
+      double a = x + 1 / scale;
 
-  for (i_offset=0; i_offset<SUBSAMPLE; i_offset++)
-    for (j_offset=0; j_offset<SUBSAMPLE; j_offset++)
-      {
-	int *pixel_weights = filter->weights + ((i_offset*SUBSAMPLE) + j_offset) * n_x * n_y;
-	double x = (double)j_offset / SUBSAMPLE;
-	double y = (double)i_offset / SUBSAMPLE;
-	int i,j;
-	int total = 0;
-	  
-	for (i = 0; i < n_y; i++)
-	  {
-	    double tw, th;
-		
-	    if (i < y)
-	      {
-
-		if (i + 1 > y)
-		  th = MIN(i+1, y + 1/y_scale) - y;
-		else
-		  th = 0;
-	      }
-	    else
-	      {
-		if (y + 1/y_scale > i)
-		  th = MIN(i+1, y + 1/y_scale) - i;
-		else
-		  th = 0;
-	      }
-		
-	    for (j = 0; j < n_x; j++)
-	      {
-		int weight;
-		
-		if (j < x)
-		  {
-		    if (j + 1 > x)
-		      tw = MIN(j+1, x + 1/x_scale) - x;
-		    else
-		      tw = 0;
-		  }
-		else
-		  {
-		    if (x + 1/x_scale > j)
-		      tw = MIN(j+1, x + 1/x_scale) - j;
-		    else
-		      tw = 0;
-		  }
-
-		weight = 65536 * tw * x_scale * th * y_scale * overall_alpha + 0.5;
-		total += weight;
-		*(pixel_weights + n_x * i + j) = weight;
-	      }
-	  }
-	
-	correct_total (pixel_weights, n_x, n_y, total, overall_alpha);
-      }
+      for (i = 0; i < n; i++)
+        {
+          if (i < x)
+            {
+              if (i + 1 > x)
+                *(pixel_weights++)  = (MIN (i + 1, a) - x) * scale;
+              else
+                *(pixel_weights++) = 0;
+            }
+          else
+            {
+              if (a > i)
+                *(pixel_weights++)  = (MIN (i + 1, a) - i) * scale;
+              else
+                *(pixel_weights++) = 0;
+            }
+       }
+    }
 }
 
+/* Compute weights for a filter that, for minification
+ * is the same as 'tiles', and for magnification, is bilinear
+ * reconstruction followed by a sampling with a delta function.
+ */
 static void
-bilinear_make_fast_weights (PixopsFilter *filter, double x_scale, double y_scale, double overall_alpha)
+bilinear_magnify_make_weights (PixopsFilterDimension *dim,
+			       double                 scale)
 {
-  int i_offset, j_offset;
-  double *x_weights, *y_weights;
-  int n_x, n_y;
-
-  if (x_scale > 1.0)		/* Bilinear */
-    {
-      n_x = 2;
-      filter->x_offset = 0.5 * (1/x_scale - 1);
-    }
-  else				/* Tile */
-    {
-      n_x = ceil(1.0 + 1.0/x_scale);
-      filter->x_offset = 0.0;
-    }
+  double *pixel_weights;
+  int n;
+  int offset;
+  int i;
 
-  if (y_scale > 1.0)		/* Bilinear */
+  if (scale > 1.0)            /* Linear */
     {
-      n_y = 2;
-      filter->y_offset = 0.5 * (1/y_scale - 1);
+      n = 2;
+      dim->offset = 0.5 * (1 / scale - 1);
     }
-  else				/* Tile */
+  else                          /* Tile */
     {
-      n_y = ceil(1.0 + 1.0/y_scale);
-      filter->y_offset = 0.0;
+      n = ceil (1.0 + 1.0 / scale);
+      dim->offset = 0.0;
     }
 
-  filter->n_y = n_y;
-  filter->n_x = n_x;
-  filter->weights = g_new (int, SUBSAMPLE * SUBSAMPLE * n_x * n_y);
+  dim->n = n;
+  dim->weights = g_new (double, SUBSAMPLE * n);
 
-  x_weights = g_new (double, n_x);
-  y_weights = g_new (double, n_y);
+  pixel_weights = dim->weights;
 
-  for (i_offset=0; i_offset<SUBSAMPLE; i_offset++)
-    for (j_offset=0; j_offset<SUBSAMPLE; j_offset++)
-      {
-	int *pixel_weights = filter->weights + ((i_offset*SUBSAMPLE) + j_offset) * n_x * n_y;
-	double x = (double)j_offset / SUBSAMPLE;
-	double y = (double)i_offset / SUBSAMPLE;
-	int i,j;
-	int total = 0;
-
-	if (x_scale > 1.0)	/* Bilinear */
-	  {
-	    for (i = 0; i < n_x; i++)
-	      {
-		x_weights[i] = ((i == 0) ? (1 - x) : x) / x_scale;
-	      }
-	  }
-	else			/* Tile */
-	  {
-	    /*           x
-	     * ---------|--.-|----|--.-|-------  SRC
-	     * ------------|---------|---------  DEST
-	     */
-	    for (i = 0; i < n_x; i++)
-	      {
-		if (i < x)
-		  {
-		    if (i + 1 > x)
-		      x_weights[i] = MIN(i+1, x + 1/x_scale) - x;
-		    else
-		      x_weights[i] = 0;
-		  }
-		else
-		  {
-		    if (x + 1/x_scale > i)
-		      x_weights[i] = MIN(i+1, x + 1/x_scale) - i;
-		    else
-		      x_weights[i] = 0;
-		  }
-	      }
-	  }
-
-	if (y_scale > 1.0)	/* Bilinear */
-	  {
-	    for (i = 0; i < n_y; i++)
-	      {
-		y_weights[i] = ((i == 0) ? (1 - y) : y) / y_scale;
-	      }
-	  }
-	else			/* Tile */
-	  {
-	    /*           y
-	     * ---------|--.-|----|--.-|-------  SRC
-	     * ------------|---------|---------  DEST
-	     */
-	    for (i = 0; i < n_y; i++)
-	      {
-		if (i < y)
-		  {
-		    if (i + 1 > y)
-		      y_weights[i] = MIN(i+1, y + 1/y_scale) - y;
-		    else
-		      y_weights[i] = 0;
-		  }
-		else
-		  {
-		    if (y + 1/y_scale > i)
-		      y_weights[i] = MIN(i+1, y + 1/y_scale) - i;
-		    else
-		      y_weights[i] = 0;
-		  }
-	      }
-	  }
-
-	for (i = 0; i < n_y; i++)
-	  for (j = 0; j < n_x; j++)
-	    {
-	      int weight = 65536 * x_weights[j] * x_scale * y_weights[i] * y_scale * overall_alpha + 0.5;
-	      *(pixel_weights + n_x * i + j) = weight;
-	      total += weight;
-	    }
-
-	correct_total (pixel_weights, n_x, n_y, total, overall_alpha);
-      }
-
-  g_free (x_weights);
-  g_free (y_weights);
+  for (offset=0; offset < SUBSAMPLE; offset++)
+    {
+      double x = (double)offset / SUBSAMPLE;
+
+      if (scale > 1.0)      /* Linear */
+        {
+          for (i = 0; i < n; i++)
+            *(pixel_weights++) = (((i == 0) ? (1 - x) : x) / scale) * scale;
+        }
+      else                  /* Tile */
+        {
+          double a = x + 1 / scale;
+
+          /*           x
+           * ---------|--.-|----|--.-|-------  SRC
+           * ------------|---------|---------  DEST
+           */
+          for (i = 0; i < n; i++)
+            {
+              if (i < x)
+                {
+                  if (i + 1 > x)
+                    *(pixel_weights++) = (MIN (i + 1, a) - x) * scale;
+                  else
+                    *(pixel_weights++) = 0;
+                }
+              else
+                {
+                  if (a > i)
+                    *(pixel_weights++) = (MIN (i + 1, a) - i) * scale;
+                  else
+                    *(pixel_weights++) = 0;
+                }
+            }
+        }
+    }
 }
 
+/* Computes the integral from b0 to b1 of
+ *
+ * f(x) = x; 0 <= x < 1
+ * f(x) = 0; otherwise
+ *
+ * We combine two of these to compute the convolution of
+ * a box filter with a triangular spike.
+ */
 static double
-bilinear_quadrant (double bx0, double bx1, double by0, double by1)
+linear_box_half (double b0, double b1)
 {
-  double ax0, ax1, ay0, ay1;
-  double x0, x1, y0, y1;
+  double a0, a1;
+  double x0, x1;
 
-  ax0 = 0.;
-  ax1 = 1.;
-  ay0 = 0.;
-  ay1 = 1.;
-
-  if (ax0 < bx0)
-    {
-      if (ax1 > bx0)
-	{
-	  x0 = bx0;
-	  x1 = MIN (ax1, bx1);
-	}
-      else
-	return 0;
-    }
-  else
-    {
-      if (bx1 > ax0)
-	{
-	  x0 = ax0;
-	  x1 = MIN (ax1, bx1);
-	}
-      else
-	return 0;
-    }
+  a0 = 0.;
+  a1 = 1.;
 
-  if (ay0 < by0)
+  if (a0 < b0)
     {
-      if (ay1 > by0)
-	{
-	  y0 = by0;
-	  y1 = MIN (ay1, by1);
-	}
+      if (a1 > b0)
+        {
+          x0 = b0;
+          x1 = MIN (a1, b1);
+        }
       else
-	return 0;
+        return 0;
     }
   else
     {
-      if (by1 > ay0)
-	{
-	  y0 = ay0;
-	  y1 = MIN (ay1, by1);
-	}
+      if (b1 > a0)
+        {
+          x0 = a0;
+          x1 = MIN (a1, b1);
+        }
       else
-	return 0;
+        return 0;
     }
 
-  return 0.25 * (x1*x1 - x0*x0) * (y1*y1 - y0*y0);
+  return 0.5 * (x1*x1 - x0*x0);
 }
 
+/* Compute weights for reconstructing with bilinear
+ * interpolation, then sampling with a box filter
+ */
 static void
-bilinear_make_weights (PixopsFilter *filter, double x_scale, double y_scale, double overall_alpha)
+bilinear_box_make_weights (PixopsFilterDimension *dim,
+			   double                 scale)
 {
-  int i_offset, j_offset;
+  int n = ceil (1/scale + 2.0);
+  double *pixel_weights = g_new (double, SUBSAMPLE * n);
+  double w;
+  int offset, i;
 
-  int n_x = ceil(1/x_scale + 2.0);
-  int n_y = ceil(1/y_scale + 2.0);
+  dim->offset = -1.0;
+  dim->n = n;
+  dim->weights = pixel_weights;
 
-  filter->x_offset = -1.0;
-  filter->y_offset = -1.0;
-  filter->n_x = n_x;
-  filter->n_y = n_y;
-  
-  filter->weights = g_new (int, SUBSAMPLE * SUBSAMPLE * n_x * n_y);
+  for (offset =0 ; offset < SUBSAMPLE; offset++)
+    {
+      double x = (double)offset / SUBSAMPLE;
+      double a = x + 1 / scale;
 
-  for (i_offset=0; i_offset<SUBSAMPLE; i_offset++)
-    for (j_offset=0; j_offset<SUBSAMPLE; j_offset++)
-      {
-	int *pixel_weights = filter->weights + ((i_offset*SUBSAMPLE) + j_offset) * n_x * n_y;
-	double x = (double)j_offset / SUBSAMPLE;
-	double y = (double)i_offset / SUBSAMPLE;
-	int i,j;
-	int total = 0;
-
-	for (i = 0; i < n_y; i++)
-	  for (j = 0; j < n_x; j++)
-	    {
-	      double w;
-	      int weight;
-	      
-	      w = bilinear_quadrant  (0.5 + j - (x + 1 / x_scale), 0.5 + j - x, 0.5 + i - (y + 1 / y_scale), 0.5 + i - y);
-	      w += bilinear_quadrant (1.5 + x - j, 1.5 + (x + 1 / x_scale) - j, 0.5 + i - (y + 1 / y_scale), 0.5 + i - y);
-	      w += bilinear_quadrant (0.5 + j - (x + 1 / x_scale), 0.5 + j - x, 1.5 + y - i, 1.5 + (y + 1 / y_scale) - i);
-	      w += bilinear_quadrant (1.5 + x - j, 1.5 + (x + 1 / x_scale) - j, 1.5 + y - i, 1.5 + (y + 1 / y_scale) - i);
-	      weight = 65536 * w * x_scale * y_scale * overall_alpha + 0.5;
-	      *(pixel_weights + n_x * i + j) = weight;
-	      total += weight;
-	    }
-	
-	correct_total (pixel_weights, n_x, n_y, total, overall_alpha);
-      }
+      for (i = 0; i < n; i++)
+        {
+          w  = linear_box_half (0.5 + i - a, 0.5 + i - x);
+          w += linear_box_half (1.5 + x - i, 1.5 + a - i);
+      
+          *(pixel_weights++) = w * scale;
+        }
+    }
+}
+
+static void
+make_weights (PixopsFilter     *filter,
+	      PixopsInterpType  interp_type,	      
+	      double            scale_x,
+	      double            scale_y)
+{
+  switch (interp_type)
+    {
+    case PIXOPS_INTERP_NEAREST:
+      g_assert_not_reached ();
+      break;
+
+    case PIXOPS_INTERP_TILES:
+      tile_make_weights (&filter->x, scale_x);
+      tile_make_weights (&filter->y, scale_y);
+      break;
+      
+    case PIXOPS_INTERP_BILINEAR:
+      bilinear_magnify_make_weights (&filter->x, scale_x);
+      bilinear_magnify_make_weights (&filter->y, scale_y);
+      break;
+      
+    case PIXOPS_INTERP_HYPER:
+      bilinear_box_make_weights (&filter->x, scale_x);
+      bilinear_box_make_weights (&filter->y, scale_y);
+      break;
+    }
 }
 
 void
@@ -1430,38 +1393,28 @@ pixops_composite_color (guchar         *dest_buf,
 
   if (!src_has_alpha && overall_alpha == 255)
     {
-    pixops_scale (dest_buf, render_x0, render_y0, render_x1, render_y1,
-		  dest_rowstride, dest_channels, dest_has_alpha,
-		  src_buf, src_width, src_height, src_rowstride, src_channels,
-		  src_has_alpha, scale_x, scale_y, interp_type);
+      pixops_scale (dest_buf, render_x0, render_y0, render_x1, render_y1,
+		    dest_rowstride, dest_channels, dest_has_alpha,
+		    src_buf, src_width, src_height, src_rowstride, src_channels,
+		    src_has_alpha, scale_x, scale_y, interp_type);
       return;
     }
 
-  switch (interp_type)
+  if (interp_type == PIXOPS_INTERP_NEAREST)
     {
-    case PIXOPS_INTERP_NEAREST:
       pixops_composite_color_nearest (dest_buf, render_x0, render_y0, render_x1, render_y1,
 				      dest_rowstride, dest_channels, dest_has_alpha,
 				      src_buf, src_width, src_height, src_rowstride, src_channels, src_has_alpha,
 				      scale_x, scale_y, overall_alpha,
 				      check_x, check_y, check_size, color1, color2);
       return;
-
-    case PIXOPS_INTERP_TILES:
-      tile_make_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
-      
-    case PIXOPS_INTERP_BILINEAR:
-      bilinear_make_fast_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
-      
-    case PIXOPS_INTERP_HYPER:
-      bilinear_make_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
     }
+  
+  filter.overall_alpha = overall_alpha / 255.;
+  make_weights (&filter, interp_type, scale_x, scale_y);
 
 #ifdef USE_MMX
-  if (filter.n_x == 2 && filter.n_y == 2 &&
+  if (filter.x.n == 2 && filter.y.n == 2 &&
       dest_channels == 4 && src_channels == 4 && src_has_alpha && !dest_has_alpha && found_mmx)
     line_func = composite_line_color_22_4a4_mmx_stub;
   else
@@ -1474,7 +1427,8 @@ pixops_composite_color (guchar         *dest_buf,
 		  src_has_alpha, scale_x, scale_y, check_x, check_y, check_size, color1, color2,
 		  &filter, line_func, composite_pixel_color);
 
-  g_free (filter.weights);
+  g_free (filter.x.weights);
+  g_free (filter.y.weights);
 }
 
 /**
@@ -1536,36 +1490,26 @@ pixops_composite (guchar        *dest_buf,
 
   if (!src_has_alpha && overall_alpha == 255)
     {
-    pixops_scale (dest_buf, render_x0, render_y0, render_x1, render_y1,
-		  dest_rowstride, dest_channels, dest_has_alpha,
-		  src_buf, src_width, src_height, src_rowstride, src_channels,
-		  src_has_alpha, scale_x, scale_y, interp_type);
+      pixops_scale (dest_buf, render_x0, render_y0, render_x1, render_y1,
+		    dest_rowstride, dest_channels, dest_has_alpha,
+		    src_buf, src_width, src_height, src_rowstride, src_channels,
+		    src_has_alpha, scale_x, scale_y, interp_type);
       return;
     }
 
-  switch (interp_type)
+  if (interp_type == PIXOPS_INTERP_NEAREST)
     {
-    case PIXOPS_INTERP_NEAREST:
       pixops_composite_nearest (dest_buf, render_x0, render_y0, render_x1, render_y1,
 				dest_rowstride, dest_channels, dest_has_alpha,
 				src_buf, src_width, src_height, src_rowstride, src_channels,
 				src_has_alpha, scale_x, scale_y, overall_alpha);
       return;
-
-    case PIXOPS_INTERP_TILES:
-      tile_make_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
-      
-    case PIXOPS_INTERP_BILINEAR:
-      bilinear_make_fast_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
-      
-    case PIXOPS_INTERP_HYPER:
-      bilinear_make_weights (&filter, scale_x, scale_y, overall_alpha / 255.);
-      break;
     }
+  
+  filter.overall_alpha = overall_alpha / 255.;
+  make_weights (&filter, interp_type, scale_x, scale_y);
 
-  if (filter.n_x == 2 && filter.n_y == 2 &&
+  if (filter.x.n == 2 && filter.y.n == 2 &&
       dest_channels == 4 && src_channels == 4 && src_has_alpha && !dest_has_alpha)
     {
 #ifdef USE_MMX
@@ -1584,7 +1528,8 @@ pixops_composite (guchar        *dest_buf,
 		  src_has_alpha, scale_x, scale_y, 0, 0, 0, 0, 0, 
 		  &filter, line_func, composite_pixel);
 
-  g_free (filter.weights);
+  g_free (filter.x.weights);
+  g_free (filter.y.weights);
 }
 
 void
@@ -1620,29 +1565,19 @@ pixops_scale (guchar        *dest_buf,
   if (scale_x == 0 || scale_y == 0)
     return;
 
-  switch (interp_type)
+  if (interp_type == PIXOPS_INTERP_NEAREST)
     {
-    case PIXOPS_INTERP_NEAREST:
       pixops_scale_nearest (dest_buf, render_x0, render_y0, render_x1, render_y1,
 			    dest_rowstride, dest_channels, dest_has_alpha,
 			    src_buf, src_width, src_height, src_rowstride, src_channels, src_has_alpha,
 			    scale_x, scale_y);
       return;
-
-    case PIXOPS_INTERP_TILES:
-      tile_make_weights (&filter, scale_x, scale_y, 1.0);
-      break;
-      
-    case PIXOPS_INTERP_BILINEAR:
-      bilinear_make_fast_weights (&filter, scale_x, scale_y, 1.0);
-      break;
-      
-    case PIXOPS_INTERP_HYPER:
-      bilinear_make_weights (&filter, scale_x, scale_y, 1.0);
-      break;
     }
+  
+  filter.overall_alpha = 1.0;
+  make_weights (&filter, interp_type, scale_x, scale_y);
 
-  if (filter.n_x == 2 && filter.n_y == 2 && dest_channels == 3 && src_channels == 3)
+  if (filter.x.n == 2 && filter.y.n == 2 && dest_channels == 3 && src_channels == 3)
     {
 #ifdef USE_MMX
       if (found_mmx)
@@ -1660,6 +1595,7 @@ pixops_scale (guchar        *dest_buf,
 		  src_has_alpha, scale_x, scale_y, 0, 0, 0, 0, 0,
 		  &filter, line_func, scale_pixel);
 
-  g_free (filter.weights);
+  g_free (filter.x.weights);
+  g_free (filter.y.weights);
 }