Merge branch 'add-reduce'

ChangeLog

@@ -1,3 +1,8 @@
+29/1/16 started 8.3
+- add vips_reduce*() ... a fast path for bicubic downsize
+- vips_resize() and vips_similarity use it when they can
+- bicubic is better on 32-bit int images
+
 27/1/16 started 8.2.3
 - fix a crash with SPARC byte-order labq vips images
 

TODO

@@ -1,15 +1,22 @@
-- look for FCLIP / FABS use ... VIPS_ABS() on float type
+- new vips_reduce:
 
-- write a shrinker that does two 1D shrinks, vertical and horizontal, with
+	affine
-  bicubic interpolation
+	$ time vipsthumbnail wtc.tif -o x.tif -s 7000
+	real	0m2.325s
+	user	0m8.316s
+	sys	0m0.256s
 
-  should be much quicker than affine + bicubic
+	reduce
-
+	$ time vipsthumbnail wtc.tif -o x.tif -s 7000
-  vips_resize() can use this, for the default path at least
+	real	0m1.677s
+	user	0m5.996s
+	sys	0m0.300s
 
 - get some brightly coloured spots with nohalo / vsqbs on wobble.ws ... very
   odd, the interpolation shouldn't change between bands
 
+  probably over/underflow in interpolation ... clipping problems?
+
 - still not happy about float->int mask conversion in im_vips2mask.c
 
 - need to follow up on

configure.ac

@@ -2,7 +2,7 @@
 
 # also update the version number in the m4 macros below
 
-AC_INIT([vips], [8.2.3], [vipsip@jiscmail.ac.uk])
+AC_INIT([vips], [8.3.0], [vipsip@jiscmail.ac.uk])
 # required for gobject-introspection
 AC_PREREQ(2.62)
 
@@ -17,8 +17,8 @@ AC_CONFIG_MACRO_DIR([m4])
 
 # user-visible library versioning
 m4_define([vips_major_version], [8])
-m4_define([vips_minor_version], [2])
+m4_define([vips_minor_version], [3])
-m4_define([vips_micro_version], [3])
+m4_define([vips_micro_version], [0])
 m4_define([vips_version],
           [vips_major_version.vips_minor_version.vips_micro_version])
 
@@ -37,9 +37,9 @@ VIPS_VERSION_STRING=$VIPS_VERSION-`date`
 #   binary interface changes backwards compatible?: increment age
 #   binary interface changes not backwards compatible?: reset age to 0
 
-LIBRARY_CURRENT=45
+LIBRARY_CURRENT=46
-LIBRARY_REVISION=3
+LIBRARY_REVISION=0
-LIBRARY_AGE=3
+LIBRARY_AGE=4
 
 # patched into include/vips/version.h
 AC_SUBST(VIPS_VERSION)

cplusplus/include/vips/vips-operators.h

@@ -1,5 +1,5 @@
 // headers for vips operations
-// Sat Jan  9 15:06:40 GMT 2016
+// Sun Feb  7 16:27:14 GMT 2016
 // this file is generated automatically, do not edit!
 
 static void system( char * cmd_format , VOption *options = 0 );
@@ -146,7 +146,9 @@ VImage mapim( VImage index , VOption *options = 0 );
 VImage shrink( double xshrink , double yshrink , VOption *options = 0 );
 VImage shrinkh( int xshrink , VOption *options = 0 );
 VImage shrinkv( int yshrink , VOption *options = 0 );
-VImage shrink2( double xshrink , double yshrink , VOption *options = 0 );
+VImage reduceh( double xshrink , VOption *options = 0 );
+VImage reducev( double yshrink , VOption *options = 0 );
+VImage reduce( double xshrink , double yshrink , VOption *options = 0 );
 VImage quadratic( VImage coeff , VOption *options = 0 );
 VImage affine( std::vector<double> matrix , VOption *options = 0 );
 VImage similarity( VOption *options = 0 );

cplusplus/vips-operators.cpp

@@ -1,5 +1,5 @@
 // bodies for vips operations
-// Sat Jan  9 15:05:58 GMT 2016
+// Sun Feb  7 16:26:59 GMT 2016
 // this file is generated automatically, do not edit!
 
 void VImage::system( char * cmd_format , VOption *options )
@@ -1783,11 +1783,37 @@ VImage VImage::shrinkv( int yshrink , VOption *options )
     return( out );
 }
 
-VImage VImage::shrink2( double xshrink , double yshrink , VOption *options )
+VImage VImage::reduceh( double xshrink , VOption *options )
 {
     VImage out;
 
-    call( "shrink2" ,
+    call( "reduceh" ,
+        (options ? options : VImage::option()) ->
+            set( "in", *this ) ->
+            set( "out", &out ) ->
+            set( "xshrink", xshrink ) );
+
+    return( out );
+}
+
+VImage VImage::reducev( double yshrink , VOption *options )
+{
+    VImage out;
+
+    call( "reducev" ,
+        (options ? options : VImage::option()) ->
+            set( "in", *this ) ->
+            set( "out", &out ) ->
+            set( "yshrink", yshrink ) );
+
+    return( out );
+}
+
+VImage VImage::reduce( double xshrink , double yshrink , VOption *options )
+{
+    VImage out;
+
+    call( "reduce" ,
         (options ? options : VImage::option()) ->
             set( "in", *this ) ->
             set( "out", &out ) ->

doc/libvips-docs.xml

@@ -9,7 +9,7 @@
   <bookinfo>
     <title>VIPS Reference Manual</title>
     <releaseinfo>
-      For VIPS 8.2.2.
+      For VIPS 8.3.0.
       The latest version of this documentation can be found on the 
       <ulink role="online-location"
         url="http://www.vips.ecs.soton.ac.uk/index.php?title=Documentation">VIPS website</ulink>.

libvips/include/vips/resample.h

@@ -43,8 +43,11 @@ int vips_shrink( VipsImage *in, VipsImage **out,
 	__attribute__((sentinel));
 int vips_shrinkh( VipsImage *in, VipsImage **out, int xshrink, ... );
 int vips_shrinkv( VipsImage *in, VipsImage **out, int yshrink, ... );
-int vips_shrink2( VipsImage *in, VipsImage **out, 
+
+int vips_reduce( VipsImage *in, VipsImage **out, 
 	double xshrink, double yshrink, ... );
+int vips_reduceh( VipsImage *in, VipsImage **out, double xshrink, ... );
+int vips_reducev( VipsImage *in, VipsImage **out, double yshrink, ... );
 
 int vips_similarity( VipsImage *in, VipsImage **out, ... )
 	__attribute__((sentinel));

libvips/resample/Makefile.am

@@ -9,7 +9,9 @@ libresample_la_SOURCES = \
 	shrink.c \
 	shrinkh.c \
 	shrinkv.c \
-	shrink2.c \
+	reduce.c \
+	reduceh.cpp \
+	reducev.cpp \
 	interpolate.c \
 	transform.c \
 	bicubic.cpp \

libvips/resample/bicubic.cpp

@@ -2,6 +2,8 @@
  *
  * 12/8/10
  * 	- revise window_size / window_offset stuff again
+ * 7/2/16
+ * 	- double intermediate for 32-bit int types
  */
 
 /*
@@ -231,7 +233,7 @@ bicubic_signed_int_tab( void *pout, const VipsPel *pin,
 	}
 }
 
-/* Floating-point version, for int/float types.
+/* Floating-point version.
  */
 template <typename T>
 static void inline
@@ -294,6 +296,135 @@ bicubic_float_tab( void *pout, const VipsPel *pin,
 	}
 }
 
+/* uint32 version needs a double intermediate.
+ */
+
+template <typename T, int max_value>
+static void inline
+bicubic_unsigned_int32_tab( void *pout, const VipsPel *pin,
+	const int bands, const int lskip,
+	const double *cx, const double *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+        const int l1_plus_b1 = l1 + b1;
+        const int l1_plus_b2 = l1 + b2;
+        const int l1_plus_b3 = l1 + b3;
+        const int l2_plus_b1 = l2 + b1;
+        const int l2_plus_b2 = l2 + b2;
+        const int l2_plus_b3 = l2 + b3;
+        const int l3_plus_b1 = l3 + b1;
+        const int l3_plus_b2 = l3 + b2;
+        const int l3_plus_b3 = l3 + b3;
+
+	for( int z = 0; z < bands; z++ ) {
+		const T uno_one = in[0];
+		const T uno_two = in[b1];
+		const T uno_thr = in[b2];
+		const T uno_fou = in[b3];
+
+		const T dos_one = in[l1];
+		const T dos_two = in[l1_plus_b1];
+		const T dos_thr = in[l1_plus_b2];
+		const T dos_fou = in[l1_plus_b3];
+
+		const T tre_one = in[l2];
+		const T tre_two = in[l2_plus_b1];
+		const T tre_thr = in[l2_plus_b2];
+		const T tre_fou = in[l2_plus_b3];
+
+		const T qua_one = in[l3];
+		const T qua_two = in[l3_plus_b1];
+		const T qua_thr = in[l3_plus_b2];
+		const T qua_fou = in[l3_plus_b3];
+
+		double bicubic = bicubic_float<double>(
+			uno_one, uno_two, uno_thr, uno_fou,
+			dos_one, dos_two, dos_thr, dos_fou,
+			tre_one, tre_two, tre_thr, tre_fou,
+			qua_one, qua_two, qua_thr, qua_fou,
+			cx, cy );
+
+		bicubic = VIPS_CLIP( 0, bicubic, max_value ); 
+
+		out[z] = bicubic;
+
+		in += 1;
+	}
+}
+
+template <typename T, int min_value, int max_value>
+static void inline
+bicubic_signed_int32_tab( void *pout, const VipsPel *pin,
+	const int bands, const int lskip,
+	const double *cx, const double *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+        const int l1_plus_b1 = l1 + b1;
+        const int l1_plus_b2 = l1 + b2;
+        const int l1_plus_b3 = l1 + b3;
+        const int l2_plus_b1 = l2 + b1;
+        const int l2_plus_b2 = l2 + b2;
+        const int l2_plus_b3 = l2 + b3;
+        const int l3_plus_b1 = l3 + b1;
+        const int l3_plus_b2 = l3 + b2;
+        const int l3_plus_b3 = l3 + b3;
+
+	for( int z = 0; z < bands; z++ ) {
+		const T uno_one = in[0];
+		const T uno_two = in[b1];
+		const T uno_thr = in[b2];
+		const T uno_fou = in[b3];
+
+		const T dos_one = in[l1];
+		const T dos_two = in[l1_plus_b1];
+		const T dos_thr = in[l1_plus_b2];
+		const T dos_fou = in[l1_plus_b3];
+
+		const T tre_one = in[l2];
+		const T tre_two = in[l2_plus_b1];
+		const T tre_thr = in[l2_plus_b2];
+		const T tre_fou = in[l2_plus_b3];
+
+		const T qua_one = in[l3];
+		const T qua_two = in[l3_plus_b1];
+		const T qua_thr = in[l3_plus_b2];
+		const T qua_fou = in[l3_plus_b3];
+
+		double bicubic = bicubic_float<double>(
+			uno_one, uno_two, uno_thr, uno_fou,
+			dos_one, dos_two, dos_thr, dos_fou,
+			tre_one, tre_two, tre_thr, tre_fou,
+			qua_one, qua_two, qua_thr, qua_fou,
+			cx, cy );
+
+		bicubic = VIPS_CLIP( min_value, bicubic, max_value ); 
+
+		out[z] = bicubic;
+
+		in += 1;
+	}
+}
+
 /* Ultra-high-quality version for double images.
  */
 template <typename T>
@@ -461,12 +592,14 @@ vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
 		break;
 
 	case VIPS_FORMAT_UINT:
-		bicubic_float_tab<unsigned int>( out, p, bands, lskip,
+		bicubic_unsigned_int32_tab<unsigned int, INT_MAX>( 
+			out, p, bands, lskip,
 			cxf, cyf );
 		break;
 
 	case VIPS_FORMAT_INT:
-		bicubic_float_tab<signed int>( out, p, bands, lskip,
+		bicubic_signed_int32_tab<signed int, INT_MIN, INT_MAX>( 
+			out, p, bands, lskip,
 			cxf, cyf );
 		break;
 

libvips/resample/reduce.c

@@ -0,0 +1,156 @@
+/* 2D reduce ... call reduceh and reducev
+ *
+ * 27/1/16
+ * 	- from shrink.c 
+ */
+
+/*
+
+    This file is part of VIPS.
+    
+    VIPS is free software; you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+    02110-1301  USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+/*
+#define DEBUG
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <vips/vips.h>
+#include <vips/debug.h>
+#include <vips/internal.h>
+
+#include "presample.h"
+
+typedef struct _VipsReduce {
+	VipsResample parent_instance;
+
+	double xshrink;		/* Shrink factors */
+	double yshrink;
+
+} VipsReduce;
+
+typedef VipsResampleClass VipsReduceClass;
+
+G_DEFINE_TYPE( VipsReduce, vips_reduce, VIPS_TYPE_RESAMPLE );
+
+static int
+vips_reduce_build( VipsObject *object )
+{
+	VipsResample *resample = VIPS_RESAMPLE( object );
+	VipsReduce *reduce = (VipsReduce *) object;
+	VipsImage **t = (VipsImage **) 
+		vips_object_local_array( object, 3 );
+
+	if( VIPS_OBJECT_CLASS( vips_reduce_parent_class )->build( object ) )
+		return( -1 );
+
+	if( vips_reducev( resample->in, &t[0], reduce->yshrink, NULL ) ||
+		vips_reduceh( t[0], &t[1], reduce->xshrink, NULL ) ||
+		vips_image_write( t[1], resample->out ) )
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_reduce_class_init( VipsReduceClass *class )
+{
+	GObjectClass *gobject_class = G_OBJECT_CLASS( class );
+	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class );
+	VipsOperationClass *operation_class = VIPS_OPERATION_CLASS( class );
+
+	VIPS_DEBUG_MSG( "vips_reduce_class_init\n" );
+
+	gobject_class->set_property = vips_object_set_property;
+	gobject_class->get_property = vips_object_get_property;
+
+	vobject_class->nickname = "reduce";
+	vobject_class->description = _( "reduce an image" );
+	vobject_class->build = vips_reduce_build;
+
+	operation_class->flags = VIPS_OPERATION_SEQUENTIAL;
+
+	VIPS_ARG_DOUBLE( class, "xshrink", 8, 
+		_( "Xshrink" ), 
+		_( "Horizontal shrink factor" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsReduce, xshrink ),
+		1.0, 1000000.0, 1.0 );
+
+	VIPS_ARG_DOUBLE( class, "yshrink", 9, 
+		_( "Yshrink" ), 
+		_( "Vertical shrink factor" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsReduce, yshrink ),
+		1.0, 1000000.0, 1.0 );
+
+}
+
+static void
+vips_reduce_init( VipsReduce *reduce )
+{
+}
+
+/**
+ * vips_reduce:
+ * @in: input image
+ * @out: output image
+ * @xshrink: horizontal shrink
+ * @shrinke: vertical shrink
+ * @...: %NULL-terminated list of optional named arguments
+ *
+ * Reduce @in by a pair of factors with a pair of 1D cubic interpolators. This 
+ * will not work well for shrink factors greater than two.
+ *
+ * This is a very low-level operation: see vips_resize() for a more
+ * convenient way to resize images. 
+ *
+ * This operation does not change xres or yres. The image resolution needs to
+ * be updated by the application. 
+ *
+ * See also: vips_resize(), vips_affine().
+ *
+ * Returns: 0 on success, -1 on error
+ */
+int
+vips_reduce( VipsImage *in, VipsImage **out, 
+	double xshrink, double yshrink, ... )
+{
+	va_list ap;
+	int result;
+
+	va_start( ap, yshrink );
+	result = vips_call_split( "reduce", ap, in, out, xshrink, yshrink );
+	va_end( ap );
+
+	return( result );
+}

libvips/resample/reduceh.cpp

@@ -0,0 +1,530 @@
+/* horizontal reduce by a float factor
+ *
+ * 29/1/16
+ * 	- from shrinkh.c
+ */
+
+/*
+
+    This file is part of VIPS.
+    
+    VIPS is free software; you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+    02110-1301  USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+/*
+#define DEBUG
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <vips/vips.h>
+#include <vips/debug.h>
+#include <vips/internal.h>
+
+#include "presample.h"
+#include "templates.h"
+
+typedef struct _VipsReduceh {
+	VipsResample parent_instance;
+
+	double xshrink;		/* Reduce factor */
+
+} VipsReduceh;
+
+typedef VipsResampleClass VipsReducehClass;
+
+/* Precalculated interpolation matrices. int (used for pel
+ * sizes up to short), and double (for all others). We go to
+ * scale + 1 so we can round-to-nearest safely.
+ */
+
+static int vips_reduceh_matrixi[VIPS_TRANSFORM_SCALE + 1][4];
+static double vips_reduceh_matrixf[VIPS_TRANSFORM_SCALE + 1][4];
+
+/* We need C linkage for this.
+ */
+extern "C" {
+G_DEFINE_TYPE( VipsReduceh, vips_reduceh, VIPS_TYPE_RESAMPLE );
+}
+
+template <typename T, int max_value>
+static void inline
+reduceh_unsigned_int_tab( VipsPel *pout, const VipsPel *pin,
+	const int bands, const int *cx )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const int c0 = cx[0];
+	const int c1 = cx[1];
+	const int c2 = cx[2];
+	const int c3 = cx[3];
+
+	for( int z = 0; z < bands; z++ ) {
+		int cubich = unsigned_fixed_round( 
+			c0 * in[0] +
+			c1 * in[b1] +
+			c2 * in[b2] +
+			c3 * in[b3] ); 
+
+		cubich = VIPS_CLIP( 0, cubich, max_value ); 
+
+		out[z] = cubich;
+
+		in += 1;
+	}
+}
+
+template <typename T, int min_value, int max_value>
+static void inline
+reduceh_signed_int_tab( VipsPel *pout, const VipsPel *pin,
+	const int bands, const int *cx )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const int c0 = cx[0];
+	const int c1 = cx[1];
+	const int c2 = cx[2];
+	const int c3 = cx[3];
+
+	for( int z = 0; z < bands; z++ ) {
+		int cubich = signed_fixed_round( 
+			c0 * in[0] +
+			c1 * in[b1] +
+			c2 * in[b2] +
+			c3 * in[b3] ); 
+
+		cubich = VIPS_CLIP( min_value, cubich, max_value ); 
+
+		out[z] = cubich;
+
+		in += 1;
+	}
+}
+
+/* Floating-point version, for int/float types.
+ */
+template <typename T>
+static void inline
+reduceh_float_tab( VipsPel *pout, const VipsPel *pin,
+	const int bands, const double *cx )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const double c0 = cx[0];
+	const double c1 = cx[1];
+	const double c2 = cx[2];
+	const double c3 = cx[3];
+
+	for( int z = 0; z < bands; z++ ) {
+		out[z] = 
+			c0 * in[0] +
+			c1 * in[b1] +
+			c2 * in[b2] +
+			c3 * in[b3];
+
+		in += 1;
+	}
+}
+
+/* 32-bit output needs a double intermediate.
+ */
+
+template <typename T, int max_value>
+static void inline
+reduceh_unsigned_int32_tab( VipsPel *pout, const VipsPel *pin,
+	const int bands, const double *cx )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const double c0 = cx[0];
+	const double c1 = cx[1];
+	const double c2 = cx[2];
+	const double c3 = cx[3];
+
+	for( int z = 0; z < bands; z++ ) {
+		double cubich = 
+			c0 * in[0] +
+			c1 * in[b1] +
+			c2 * in[b2] +
+			c3 * in[b3]; 
+
+		cubich = VIPS_CLIP( 0, cubich, max_value ); 
+
+		out[z] = cubich;
+
+		in += 1;
+	}
+}
+
+template <typename T, int min_value, int max_value>
+static void inline
+reduceh_signed_int32_tab( VipsPel *pout, const VipsPel *pin,
+	const int bands, const double *cx )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	const double c0 = cx[0];
+	const double c1 = cx[1];
+	const double c2 = cx[2];
+	const double c3 = cx[3];
+
+	for( int z = 0; z < bands; z++ ) {
+		double cubich = 
+			c0 * in[0] +
+			c1 * in[b1] +
+			c2 * in[b2] +
+			c3 * in[b3]; 
+
+		cubich = VIPS_CLIP( min_value, cubich, max_value ); 
+
+		out[z] = cubich;
+
+		in += 1;
+	}
+}
+
+/* Ultra-high-quality version for double images.
+ */
+template <typename T>
+static void inline
+reduceh_notab( VipsPel *pout, const VipsPel *pin,
+	const int bands, double x )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int b1 = bands;
+	const int b2 = b1 + b1;
+	const int b3 = b1 + b2;
+
+	double cx[4];
+
+	calculate_coefficients_catmull( x, cx );
+
+	for( int z = 0; z < bands; z++ ) {
+		const T one = in[0];
+		const T two = in[b1];
+		const T thr = in[b2];
+		const T fou = in[b3];
+
+		out[z] = cubic_float<T>( one, two, thr, fou, cx );
+
+		in += 1;
+	}
+}
+
+static int
+vips_reduceh_gen( VipsRegion *out_region, void *seq, 
+	void *a, void *b, gboolean *stop )
+{
+	VipsImage *in = (VipsImage *) a;
+	VipsReduceh *reduceh = (VipsReduceh *) b;
+	const int ps = VIPS_IMAGE_SIZEOF_PEL( in );
+	VipsRegion *ir = (VipsRegion *) seq;
+	VipsRect *r = &out_region->valid;
+
+	/* Double bands for complex.
+	 */
+	const int bands = in->Bands * 
+		(vips_band_format_iscomplex( in->BandFmt ) ?  2 : 1);
+
+	VipsRect s;
+
+#ifdef DEBUG
+	printf( "vips_reduceh_gen: generating %d x %d at %d x %d\n",
+		r->width, r->height, r->left, r->top ); 
+#endif /*DEBUG*/
+
+	s.left = r->left * reduceh->xshrink;
+	s.top = r->top;
+	s.width = r->width * reduceh->xshrink + 4;
+	s.height = r->height;
+	if( vips_region_prepare( ir, &s ) )
+		return( -1 );
+
+	VIPS_GATE_START( "vips_reduceh_gen: work" ); 
+
+	for( int y = 0; y < r->height; y ++ ) { 
+		VipsPel *q;
+		double X;
+
+		q = VIPS_REGION_ADDR( out_region, r->left, r->top + y );
+		X = r->left * reduceh->xshrink;
+
+		for( int x = 0; x < r->width; x++ ) {
+			int ix = (int) X;
+			VipsPel *p = VIPS_REGION_ADDR( ir, ix, r->top + y );
+			const int sx = X * VIPS_TRANSFORM_SCALE * 2;
+			const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
+			const int tx = (six + 1) >> 1;
+			const int *cxi = vips_reduceh_matrixi[tx];
+			const double *cxf = vips_reduceh_matrixf[tx];
+
+			switch( in->BandFmt ) {
+			case VIPS_FORMAT_UCHAR:
+				reduceh_unsigned_int_tab
+					<unsigned char, UCHAR_MAX>(
+					q, p, bands, cxi );
+				break;
+
+			case VIPS_FORMAT_CHAR:
+				reduceh_signed_int_tab
+					<signed char, SCHAR_MIN, SCHAR_MAX>(
+					q, p, bands, cxi );
+				break;
+
+			case VIPS_FORMAT_USHORT:
+				reduceh_unsigned_int_tab
+					<unsigned short, USHRT_MAX>(
+					q, p, bands, cxi );
+				break;
+
+			case VIPS_FORMAT_SHORT:
+				reduceh_signed_int_tab
+					<signed short, SHRT_MIN, SHRT_MAX>(
+					q, p, bands, cxi );
+				break;
+
+			case VIPS_FORMAT_UINT:
+				reduceh_unsigned_int32_tab
+					<unsigned int, INT_MAX>(
+					q, p, bands, cxf );
+				break;
+
+			case VIPS_FORMAT_INT:
+				reduceh_signed_int32_tab
+					<signed int, INT_MIN, INT_MAX>(
+					q, p, bands, cxf );
+				break;
+
+			case VIPS_FORMAT_FLOAT:
+			case VIPS_FORMAT_COMPLEX:
+				reduceh_float_tab<float>( q, p, bands, cxf );
+				break;
+
+			case VIPS_FORMAT_DOUBLE:
+			case VIPS_FORMAT_DPCOMPLEX:
+				reduceh_notab<double>( q, p, bands, X - ix );
+				break;
+
+			default:
+				g_assert_not_reached();
+				break;
+			}
+
+			X += reduceh->xshrink;
+			q += ps;
+		}
+	}
+
+	VIPS_GATE_STOP( "vips_reduceh_gen: work" ); 
+
+	return( 0 );
+}
+
+static int
+vips_reduceh_build( VipsObject *object )
+{
+	VipsObjectClass *object_class = VIPS_OBJECT_GET_CLASS( object );
+	VipsResample *resample = VIPS_RESAMPLE( object );
+	VipsReduceh *reduceh = (VipsReduceh *) object;
+	VipsImage **t = (VipsImage **) 
+		vips_object_local_array( object, 2 );
+
+	VipsImage *in;
+
+	if( VIPS_OBJECT_CLASS( vips_reduceh_parent_class )->build( object ) )
+		return( -1 );
+
+	in = resample->in; 
+
+	if( reduceh->xshrink < 1 ) { 
+		vips_error( object_class->nickname, 
+			"%s", _( "reduce factors should be >= 1" ) );
+		return( -1 );
+	}
+	if( reduceh->xshrink > 3 )  
+		vips_warn( object_class->nickname, 
+			"%s", _( "reduce factor greater than 3" ) );
+
+	if( reduceh->xshrink == 1 ) 
+		return( vips_image_write( in, resample->out ) );
+
+	/* Unpack for processing.
+	 */
+	if( vips_image_decode( in, &t[0] ) )
+		return( -1 );
+	in = t[0];
+
+	/* Add new pixels around the input so we can interpolate at the edges.
+	 */
+	if( vips_embed( in, &t[1], 
+		1, 0, 
+		in->Xsize + 3, in->Ysize,
+		"extend", VIPS_EXTEND_COPY,
+		NULL ) )
+		return( -1 );
+	in = t[1];
+
+	if( vips_image_pipelinev( resample->out, 
+		VIPS_DEMAND_STYLE_THINSTRIP, in, NULL ) )
+		return( -1 );
+
+	/* Size output. Note: we round the output width down!
+	 *
+	 * Don't change xres/yres, leave that to the application layer. For
+	 * example, vipsthumbnail knows the true reduce factor (including the
+	 * fractional part), we just see the integer part here.
+	 */
+	resample->out->Xsize = (in->Xsize - 3) / reduceh->xshrink;
+	if( resample->out->Xsize <= 0 ) { 
+		vips_error( object_class->nickname, 
+			"%s", _( "image has shrunk to nothing" ) );
+		return( -1 );
+	}
+
+#ifdef DEBUG
+	printf( "vips_reduceh_build: reducing %d x %d image to %d x %d\n", 
+		in->Xsize, in->Ysize, 
+		resample->out->Xsize, resample->out->Ysize );  
+#endif /*DEBUG*/
+
+	if( vips_image_generate( resample->out,
+		vips_start_one, vips_reduceh_gen, vips_stop_one, 
+		in, reduceh ) )
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_reduceh_class_init( VipsReducehClass *reduceh_class )
+{
+	GObjectClass *gobject_class = G_OBJECT_CLASS( reduceh_class );
+	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( reduceh_class );
+	VipsOperationClass *operation_class = 
+		VIPS_OPERATION_CLASS( reduceh_class );
+
+	VIPS_DEBUG_MSG( "vips_reduceh_class_init\n" );
+
+	gobject_class->set_property = vips_object_set_property;
+	gobject_class->get_property = vips_object_get_property;
+
+	vobject_class->nickname = "reduceh";
+	vobject_class->description = _( "shrink an image horizontally" );
+	vobject_class->build = vips_reduceh_build;
+
+	operation_class->flags = VIPS_OPERATION_SEQUENTIAL_UNBUFFERED;
+
+	VIPS_ARG_DOUBLE( reduceh_class, "xshrink", 3, 
+		_( "Xshrink" ), 
+		_( "Horizontal shrink factor" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsReduceh, xshrink ),
+		1, 1000000, 1 );
+
+	/* Build the tables of pre-computed coefficients.
+	 */
+	for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
+		calculate_coefficients_catmull(
+			(float) x / VIPS_TRANSFORM_SCALE,
+			vips_reduceh_matrixf[x] );
+
+		for( int i = 0; i < 4; i++ )
+			vips_reduceh_matrixi[x][i] =
+				vips_reduceh_matrixf[x][i] * 
+				VIPS_INTERPOLATE_SCALE;
+	}
+
+}
+
+static void
+vips_reduceh_init( VipsReduceh *reduceh )
+{
+}
+
+/**
+ * vips_reduceh:
+ * @in: input image
+ * @out: output image
+ * @xshrink: horizontal reduce
+ * @...: %NULL-terminated list of optional named arguments
+ *
+ * Reduce @in horizontally by a float factor. The pixels in @out are
+ * interpolated with a 1D cubic mask. This operation will not work well for
+ * a reduction of more than a factor of two.
+ *
+ * This is a very low-level operation: see vips_resize() for a more
+ * convenient way to resize images. 
+ *
+ * This operation does not change xres or yres. The image resolution needs to
+ * be updated by the application. 
+ *
+ * See also: vips_shrink(), vips_resize(), vips_affine().
+ *
+ * Returns: 0 on success, -1 on error
+ */
+int
+vips_reduceh( VipsImage *in, VipsImage **out, double xshrink, ... )
+{
+	va_list ap;
+	int result;
+
+	va_start( ap, xshrink );
+	result = vips_call_split( "reduceh", ap, in, out, xshrink );
+	va_end( ap );
+
+	return( result );
+}

libvips/resample/reducev.cpp

@@ -0,0 +1,531 @@
+/* vertical reduce by a float factor
+ *
+ * 29/1/16
+ * 	- from shrinkv.c
+ */
+
+/*
+
+    This file is part of VIPS.
+    
+    VIPS is free software; you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+    02110-1301  USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+/*
+#define DEBUG
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <vips/vips.h>
+#include <vips/debug.h>
+#include <vips/internal.h>
+
+#include "presample.h"
+#include "templates.h"
+
+typedef struct _VipsReducev {
+	VipsResample parent_instance;
+
+	double yshrink;		/* Shrink factor */
+
+} VipsReducev;
+
+typedef VipsResampleClass VipsReducevClass;
+
+/* Precalculated interpolation matrices. int (used for pel
+ * sizes up to short), and double (for all others). We go to
+ * scale + 1 so we can round-to-nearest safely.
+ */
+
+static int vips_reducev_matrixi[VIPS_TRANSFORM_SCALE + 1][4];
+static double vips_reducev_matrixf[VIPS_TRANSFORM_SCALE + 1][4];
+
+/* We need C linkage for this.
+ */
+extern "C" {
+G_DEFINE_TYPE( VipsReducev, vips_reducev, VIPS_TYPE_RESAMPLE );
+}
+
+template <typename T, int max_value>
+static void inline
+reducev_unsigned_int_tab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	const int *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	const int c0 = cy[0];
+	const int c1 = cy[1];
+	const int c2 = cy[2];
+	const int c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		int cubicv = unsigned_fixed_round( 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3] ); 
+
+		cubicv = VIPS_CLIP( 0, cubicv, max_value ); 
+
+		out[z] = cubicv;
+
+		in += 1;
+	}
+}
+
+template <typename T, int min_value, int max_value>
+static void inline
+reducev_signed_int_tab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	const int *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	const int c0 = cy[0];
+	const int c1 = cy[1];
+	const int c2 = cy[2];
+	const int c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		int cubicv = signed_fixed_round( 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3] ); 
+
+		cubicv = VIPS_CLIP( min_value, cubicv, max_value ); 
+
+		out[z] = cubicv;
+
+		in += 1;
+	}
+}
+
+/* Floating-point version.
+ */
+
+template <typename T>
+static void inline
+reducev_float_tab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	const double *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	const double c0 = cy[0];
+	const double c1 = cy[1];
+	const double c2 = cy[2];
+	const double c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		out[z] = 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3]; 
+
+		in += 1;
+	}
+}
+
+/* 32-bit int version needs a double intermediate.
+ */
+
+template <typename T, int max_value>
+static void inline
+reducev_unsigned_int32_tab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	const double *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	const double c0 = cy[0];
+	const double c1 = cy[1];
+	const double c2 = cy[2];
+	const double c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		double cubicv = 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3]; 
+
+		cubicv = VIPS_CLIP( 0, cubicv, max_value ); 
+
+		out[z] = cubicv;
+
+		in += 1;
+	}
+}
+
+template <typename T, int min_value, int max_value>
+static void inline
+reducev_signed_int32_tab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	const double *cy )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	const double c0 = cy[0];
+	const double c1 = cy[1];
+	const double c2 = cy[2];
+	const double c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		double cubicv = 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3]; 
+
+		cubicv = VIPS_CLIP( min_value, cubicv, max_value ); 
+
+		out[z] = cubicv;
+
+		in += 1;
+	}
+}
+
+/* Ultra-high-quality version for double images.
+ */
+template <typename T>
+static void inline
+reducev_notab( VipsPel *pout, const VipsPel *pin,
+	const int ne, const int lskip,
+	double y )
+{
+	T* restrict out = (T *) pout;
+	const T* restrict in = (T *) pin;
+
+	const int l1 = lskip / sizeof( T );
+	const int l2 = l1 + l1;
+	const int l3 = l1 + l2;
+
+	double cy[4];
+
+	calculate_coefficients_catmull( y, cy );
+
+	const double c0 = cy[0];
+	const double c1 = cy[1];
+	const double c2 = cy[2];
+	const double c3 = cy[3];
+
+	for( int z = 0; z < ne; z++ ) {
+		out[z] = 
+			c0 * in[0] +
+			c1 * in[l1] +
+			c2 * in[l2] +
+			c3 * in[l3]; 
+
+		in += 1;
+	}
+}
+
+static int
+vips_reducev_gen( VipsRegion *out_region, void *seq, 
+	void *a, void *b, gboolean *stop )
+{
+	VipsImage *in = (VipsImage *) a;
+	VipsReducev *reducev = (VipsReducev *) b;
+	VipsRegion *ir = (VipsRegion *) seq;
+	VipsRect *r = &out_region->valid;
+
+	/* Double bands for complex.
+	 */
+	const int bands = in->Bands * 
+		(vips_band_format_iscomplex( in->BandFmt ) ?  2 : 1);
+	int ne = r->width * bands;
+
+	VipsRect s;
+
+#ifdef DEBUG
+	printf( "vips_reducev_gen: generating %d x %d at %d x %d\n",
+		r->width, r->height, r->left, r->top ); 
+#endif /*DEBUG*/
+
+	s.left = r->left;
+	s.top = r->top * reducev->yshrink;
+	s.width = r->width;
+	s.height = r->height * reducev->yshrink + 4;
+	if( vips_region_prepare( ir, &s ) )
+		return( -1 );
+
+	VIPS_GATE_START( "vips_reducev_gen: work" ); 
+
+	for( int y = 0; y < r->height; y ++ ) { 
+		VipsPel *q = VIPS_REGION_ADDR( out_region, r->left, r->top + y );
+		const double Y = (r->top + y) * reducev->yshrink; 
+		VipsPel *p = VIPS_REGION_ADDR( ir, r->left, (int) Y ); 
+		const int sy = Y * VIPS_TRANSFORM_SCALE * 2;
+		const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
+		const int ty = (siy + 1) >> 1;
+		const int *cyi = vips_reducev_matrixi[ty];
+		const double *cyf = vips_reducev_matrixf[ty];
+		const int lskip = VIPS_REGION_LSKIP( ir );
+
+		switch( in->BandFmt ) {
+		case VIPS_FORMAT_UCHAR:
+			reducev_unsigned_int_tab
+				<unsigned char, UCHAR_MAX>(
+				q, p, ne, lskip, cyi );
+			break;
+
+		case VIPS_FORMAT_CHAR:
+			reducev_signed_int_tab
+				<signed char, SCHAR_MIN, SCHAR_MAX>(
+				q, p, ne, lskip, cyi );
+			break;
+
+		case VIPS_FORMAT_USHORT:
+			reducev_unsigned_int_tab
+				<unsigned short, USHRT_MAX>(
+				q, p, ne, lskip, cyi );
+			break;
+
+		case VIPS_FORMAT_SHORT:
+			reducev_signed_int_tab
+				<signed short, SHRT_MIN, SHRT_MAX>(
+				q, p, ne, lskip, cyi );
+			break;
+
+		case VIPS_FORMAT_UINT:
+			reducev_unsigned_int32_tab
+				<unsigned int, INT_MAX>(
+				q, p, ne, lskip, cyf );
+			break;
+
+		case VIPS_FORMAT_INT:
+			reducev_signed_int32_tab
+				<signed int, INT_MIN, INT_MAX>(
+				q, p, ne, lskip, cyf );
+			break;
+
+		case VIPS_FORMAT_FLOAT:
+		case VIPS_FORMAT_COMPLEX:
+			reducev_float_tab<float>( q, p, ne, lskip, cyf );
+			break;
+
+		case VIPS_FORMAT_DPCOMPLEX:
+		case VIPS_FORMAT_DOUBLE:
+			reducev_notab<double>( q, p, ne, lskip, Y - (int) Y );
+			break;
+
+		default:
+			g_assert_not_reached();
+			break;
+		}
+	}
+
+	VIPS_GATE_STOP( "vips_reducev_gen: work" ); 
+
+	return( 0 );
+}
+
+static int
+vips_reducev_build( VipsObject *object )
+{
+	VipsObjectClass *object_class = VIPS_OBJECT_GET_CLASS( object );
+	VipsResample *resample = VIPS_RESAMPLE( object );
+	VipsReducev *reducev = (VipsReducev *) object;
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 2 );
+
+	VipsImage *in;
+
+	if( VIPS_OBJECT_CLASS( vips_reducev_parent_class )->build( object ) )
+		return( -1 );
+
+	in = resample->in; 
+
+	if( reducev->yshrink < 1 ) { 
+		vips_error( object_class->nickname, 
+			"%s", _( "reduce factors should be >= 1" ) );
+		return( -1 );
+	}
+	if( reducev->yshrink > 3 )  
+		vips_warn( object_class->nickname, 
+			"%s", _( "reduce factor greater than 3" ) );
+
+	if( reducev->yshrink == 1 ) 
+		return( vips_image_write( in, resample->out ) );
+
+	/* Unpack for processing.
+	 */
+	if( vips_image_decode( in, &t[0] ) )
+		return( -1 );
+	in = t[0];
+
+	/* Add new pixels around the input so we can interpolate at the edges.
+	 */
+	if( vips_embed( in, &t[1], 
+		0, 1, 
+		in->Xsize, in->Ysize + 3, 
+		"extend", VIPS_EXTEND_COPY,
+		NULL ) )
+		return( -1 );
+	in = t[1];
+
+	if( vips_image_pipelinev( resample->out, 
+		VIPS_DEMAND_STYLE_SMALLTILE, in, NULL ) )
+		return( -1 );
+
+	/* Size output. Note: we round the output width down!
+	 *
+	 * Don't change xres/yres, leave that to the application layer. For
+	 * example, vipsthumbnail knows the true reduce factor (including the
+	 * fractional part), we just see the integer part here.
+	 */
+	resample->out->Ysize = (in->Ysize - 3) / reducev->yshrink;
+	if( resample->out->Ysize <= 0 ) { 
+		vips_error( object_class->nickname, 
+			"%s", _( "image has shrunk to nothing" ) );
+		return( -1 );
+	}
+
+#ifdef DEBUG
+	printf( "vips_reducev_build: reducing %d x %d image to %d x %d\n", 
+		in->Xsize, in->Ysize, 
+		resample->out->Xsize, resample->out->Ysize );  
+#endif /*DEBUG*/
+
+	if( vips_image_generate( resample->out,
+		vips_start_one, vips_reducev_gen, vips_stop_one, 
+		in, reducev ) )
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_reducev_class_init( VipsReducevClass *reducev_class )
+{
+	GObjectClass *gobject_class = G_OBJECT_CLASS( reducev_class );
+	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( reducev_class );
+	VipsOperationClass *operation_class = 
+		VIPS_OPERATION_CLASS( reducev_class );
+
+	VIPS_DEBUG_MSG( "vips_reducev_class_init\n" );
+
+	gobject_class->set_property = vips_object_set_property;
+	gobject_class->get_property = vips_object_get_property;
+
+	vobject_class->nickname = "reducev";
+	vobject_class->description = _( "shrink an image vertically" );
+	vobject_class->build = vips_reducev_build;
+
+	operation_class->flags = VIPS_OPERATION_SEQUENTIAL_UNBUFFERED;
+
+	VIPS_ARG_DOUBLE( reducev_class, "yshrink", 3, 
+		_( "Xshrink" ), 
+		_( "Vertical shrink factor" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsReducev, yshrink ),
+		1, 1000000, 1 );
+
+	/* Build the tables of pre-computed coefficients.
+	 */
+	for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
+		calculate_coefficients_catmull(
+			(float) y / VIPS_TRANSFORM_SCALE,
+			vips_reducev_matrixf[y] );
+
+		for( int i = 0; i < 4; i++ )
+			vips_reducev_matrixi[y][i] =
+				vips_reducev_matrixf[y][i] * 
+				VIPS_INTERPOLATE_SCALE;
+	}
+
+}
+
+static void
+vips_reducev_init( VipsReducev *reducev )
+{
+}
+
+/**
+ * vips_reducev:
+ * @in: input image
+ * @out: output image
+ * @yshrink: horizontal reduce
+ * @...: %NULL-terminated list of optional named arguments
+ *
+ * Reduce @in vertically by a float factor. The pixels in @out are
+ * interpolated with a 1D cubic mask. This operation will not work well for
+ * a reduction of more than a factor of two.
+ *
+ * This is a very low-level operation: see vips_resize() for a more
+ * convenient way to resize images. 
+ *
+ * This operation does not change xres or yres. The image resolution needs to
+ * be updated by the application. 
+ *
+ * See also: vips_shrink(), vips_resize(), vips_affine().
+ *
+ * Returns: 0 on success, -1 on error
+ */
+int
+vips_reducev( VipsImage *in, VipsImage **out, double yshrink, ... )
+{
+	va_list ap;
+	int result;
+
+	va_start( ap, yshrink );
+	result = vips_call_split( "reducev", ap, in, out, yshrink );
+	va_end( ap );
+
+	return( result );
+}

libvips/resample/resample.c

@@ -52,6 +52,28 @@
 
 #include "presample.h"
 
+/**
+ * SECTION: resample
+ * @short_description: resample images in various ways
+ * @stability: Stable
+ * @include: vips/vips.h
+ *
+ * There are three types of operation in this section.
+ *
+ * First, vips_affine() applies an affine transform to an image. This is any
+ * sort of 2D transform which preserves straight lines; so any combination of 
+ * stretch, sheer, rotate and translate. You supply an interpolator for it to
+ * use to generate pixels, see vips_interpolate_new(). It will not produce
+ * good results for very large shrinks.
+ *
+ * Next, vips_resize() specialises in the common task of image reduce and 
+ * enlarge. It strings together combinations of vips_shrink(), vips_reduce(),
+ * vips_affine() and others to implement a general, high-quality image
+ * resizer.
+ *
+ * Finally, vips_mapim() can apply arbitrary 2D image transforms to an image.
+ */
+
 G_DEFINE_ABSTRACT_TYPE( VipsResample, vips_resample, VIPS_TYPE_OPERATION );
 
 static int
@@ -111,11 +133,13 @@ vips_resample_init( VipsResample *resample )
 void
 vips_resample_operation_init( void )
 {
-	extern GType vips_shrink_get_type( void ); 
 	extern GType vips_mapim_get_type( void ); 
+	extern GType vips_shrink_get_type( void ); 
 	extern GType vips_shrinkh_get_type( void ); 
 	extern GType vips_shrinkv_get_type( void ); 
-	extern GType vips_shrink2_get_type( void ); 
+	extern GType vips_reduce_get_type( void ); 
+	extern GType vips_reduceh_get_type( void ); 
+	extern GType vips_reducev_get_type( void ); 
 	extern GType vips_quadratic_get_type( void ); 
 	extern GType vips_affine_get_type( void ); 
 	extern GType vips_similarity_get_type( void ); 
@@ -125,7 +149,9 @@ vips_resample_operation_init( void )
 	vips_shrink_get_type(); 
 	vips_shrinkh_get_type(); 
 	vips_shrinkv_get_type(); 
-	vips_shrink2_get_type(); 
+	vips_reduceh_get_type(); 
+	vips_reducev_get_type(); 
+	vips_reduce_get_type(); 
 	vips_quadratic_get_type(); 
 	vips_affine_get_type(); 
 	vips_similarity_get_type(); 

libvips/resample/resize.c

@@ -230,12 +230,26 @@ vips_resize_build( VipsObject *object )
 		vips_info( class->nickname, "%s interpolation", nickname );
 	}
 
+	/* We have a special path for bicubic, idx/idy == 0.
+	 */
+	if( resize->interpolate &&
+		strcmp( VIPS_OBJECT_GET_CLASS( resize->interpolate )->nickname,
+			"bicubic" ) == 0 &&
+		resize->idx == 0.0 &&
+		resize->idy == 0.0 ) {
+		vips_info( class->nickname, "using fast path for residual" );
+		if( vips_reduce( in, &t[3], 
+			1.0 / hresidual, 1.0 / vresidual, NULL ) )  
+			return( -1 );
+	}
+	else {
 		if( vips_affine( in, &t[3], hresidual, 0, 0, vresidual, 
 			"interpolate", resize->interpolate,
 			"idx", resize->idx,
 			"idy", resize->idy,
 			NULL ) )  
 			return( -1 );
+	}
 	in = t[3];
 
 	/* If we are upsampling, don't sharpen. Also don't sharpen if we

libvips/resample/shrink2.c

@@ -1,461 +0,0 @@
-/* shrink with a box filter
- *
- * Copyright: 1990, N. Dessipris.
- *
- * Authors: Nicos Dessipris and Kirk Martinez
- * Written on: 29/04/1991
- * Modified on: 2/11/92, 22/2/93 Kirk Martinez - Xres Yres & cleanup 
- incredibly inefficient for box filters as LUTs are used instead of + 
- Needs converting to a smoother filter: eg Gaussian!  KM
- * 15/7/93 JC
- *	- rewritten for partial v2
- *	- ANSIfied
- *	- now shrinks any non-complex type
- *	- no longer cloned from im_convsub()
- *	- could be much better! see km comments above
- * 3/8/93 JC
- *	- rounding bug fixed
- * 11/1/94 JC
- *	- problems with .000001 and round up/down ignored! Try shrink 3738
- *	  pixel image by 9.345000000001
- * 7/10/94 JC
- *	- IM_NEW and IM_ARRAY added
- *	- more typedef
- * 3/7/95 JC
- *	- IM_CODING_LABQ handling added here
- * 20/12/08
- * 	- fall back to im_copy() for 1/1 shrink
- * 2/2/11
- * 	- gtk-doc
- * 10/2/12
- * 	- shrink in chunks to reduce peak memuse for large shrinks
- * 	- simpler
- * 12/6/12
- * 	- redone as a class
- * 	- warn about non-int shrinks
- * 	- some tuning .. tried an int coordinate path, not worthwhile
- * 16/11/12
- * 	- don't change xres/yres, see comment below
- * 8/4/13
- * 	- oops demand_hint was incorrect, thanks Jan
- * 6/6/13
- * 	- don't chunk horizontally, fixes seq problems with large shrink
- * 	  factors
- */
-
-/*
-
-    This file is part of VIPS.
-    
-    VIPS is free software; you can redistribute it and/or modify
-    it under the terms of the GNU Lesser General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU Lesser General Public License for more details.
-
-    You should have received a copy of the GNU Lesser General Public License
-    along with this program; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-    02110-1301  USA
-
- */
-
-/*
-
-    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
-
- */
-
-/*
-#define DEBUG
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif /*HAVE_CONFIG_H*/
-#include <vips/intl.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-#include <vips/vips.h>
-#include <vips/debug.h>
-#include <vips/internal.h>
-
-#include "presample.h"
-
-typedef struct _VipsShrink2 {
-	VipsResample parent_instance;
-
-	double xshrink;		/* Shrink factors */
-	double yshrink;
-
-	int mw;			/* Size of area we average */
-	int mh;
-
-	int np;			/* Number of pels we average */
-
-} VipsShrink2;
-
-typedef VipsResampleClass VipsShrink2Class;
-
-G_DEFINE_TYPE( VipsShrink2, vips_shrink2, VIPS_TYPE_RESAMPLE );
-
-/* Our per-sequence parameter struct. Somewhere to sum band elements.
- */
-typedef struct {
-	VipsRegion *ir;
-
-	VipsPel *sum;
-} VipsShrink2Sequence;
-
-/* Free a sequence value.
- */
-static int
-vips_shrink2_stop( void *vseq, void *a, void *b )
-{
-	VipsShrink2Sequence *seq = (VipsShrink2Sequence *) vseq;
-
-	VIPS_FREEF( g_object_unref, seq->ir );
-
-	return( 0 );
-}
-
-/* Make a sequence value.
- */
-static void *
-vips_shrink2_start( VipsImage *out, void *a, void *b )
-{
-	VipsImage *in = (VipsImage *) a;
-	VipsShrink2 *shrink = (VipsShrink2 *) b;
-	VipsShrink2Sequence *seq;
-
-	if( !(seq = VIPS_NEW( out, VipsShrink2Sequence )) )
-		return( NULL );
-
-	seq->ir = vips_region_new( in );
-	if( !(seq->sum = (VipsPel *) VIPS_ARRAY( out, in->Bands, double )) ) {
-		vips_shrink2_stop( seq, in, shrink );
-		return( NULL );
-	}
-
-	return( (void *) seq );
-}
-
-/* Integer shrink. 
- */
-#define ISHRINK( TYPE ) { \
-	int *sum = (int *) seq->sum; \
-	TYPE *p = (TYPE *) in; \
-	TYPE *q = (TYPE *) out; \
-	\
-	for( b = 0; b < bands; b++ ) \
-		sum[b] = 0; \
-	\
-	for( y1 = 0; y1 < shrink->mh; y1++ ) { \
-		for( i = 0, x1 = 0; x1 < shrink->mw; x1++ ) \
-			for( b = 0; b < bands; b++, i++ ) \
-				sum[b] += p[i]; \
-		\
-		p += ls; \
-	} \
-	\
-	for( b = 0; b < bands; b++ ) \
-		q[b] = (sum[b] + shrink->np / 2) / shrink->np; \
-} 
-
-/* Float shrink. 
- */
-#define FSHRINK( TYPE ) { \
-	double *sum = (double *) seq->sum; \
-	TYPE *p = (TYPE *) in; \
-	TYPE *q = (TYPE *) out; \
-	\
-	for( b = 0; b < bands; b++ ) \
-		sum[b] = 0.0; \
-	\
-	for( y1 = 0; y1 < shrink->mh; y1++ ) { \
-		for( i = 0, x1 = 0; x1 < shrink->mw; x1++ ) \
-			for( b = 0; b < bands; b++, i++ ) \
-				sum[b] += p[i]; \
-		\
-		p += ls; \
-	} \
-	\
-	for( b = 0; b < bands; b++ ) \
-		q[b] = sum[b] / shrink->np; \
-} 
-
-/* Generate an area of @or. @ir is large enough.
- */
-static void
-vips_shrink2_gen2( VipsShrink2 *shrink, VipsShrink2Sequence *seq,
-	VipsRegion *or, VipsRegion *ir,
-	int left, int top, int width, int height )
-{
-	VipsResample *resample = VIPS_RESAMPLE( shrink );
-	const int bands = resample->in->Bands;
-	const int sizeof_pixel = VIPS_IMAGE_SIZEOF_PEL( resample->in );
-	const int ls = VIPS_REGION_LSKIP( ir ) / 
-		VIPS_IMAGE_SIZEOF_ELEMENT( resample->in );
-
-	int x, y, i;
-	int x1, y1, b;
-
-	for( y = 0; y < height; y++ ) { 
-		VipsPel *out = VIPS_REGION_ADDR( or, left, top + y ); 
-
-		for( x = 0; x < width; x++ ) { 
-			int ix = (left + x) * shrink->xshrink; 
-			int iy = (top + y) * shrink->yshrink; 
-			VipsPel *in = VIPS_REGION_ADDR( ir, ix, iy ); 
-
-			switch( resample->in->BandFmt ) {
-			case VIPS_FORMAT_UCHAR: 	
-				ISHRINK( unsigned char ); break;
-			case VIPS_FORMAT_CHAR: 	
-				ISHRINK( char ); break; 
-			case VIPS_FORMAT_USHORT: 
-				ISHRINK( unsigned short ); break;
-			case VIPS_FORMAT_SHORT: 	
-				ISHRINK( short ); break; 
-			case VIPS_FORMAT_UINT: 	
-				ISHRINK( unsigned int ); break; 
-			case VIPS_FORMAT_INT: 	
-				ISHRINK( int );  break; 
-			case VIPS_FORMAT_FLOAT: 	
-				FSHRINK( float ); break; 
-			case VIPS_FORMAT_DOUBLE:	
-				FSHRINK( double ); break;
-
-			default:
-				g_assert_not_reached(); 
-			}
-
-			out += sizeof_pixel;
-		}
-	}
-}
-
-static int
-vips_shrink2_gen( VipsRegion *or, void *vseq, void *a, void *b, gboolean *stop )
-{
-	VipsShrink2Sequence *seq = (VipsShrink2Sequence *) vseq;
-	VipsShrink2 *shrink = (VipsShrink2 *) b;
-	VipsRegion *ir = seq->ir;
-	VipsRect *r = &or->valid;
-
-	/* How do we chunk up the image? We don't want to prepare the whole of
-	 * the input region corresponding to *r since it could be huge. 
-	 *
-	 * Each pixel of *r will depend on roughly mw x mh
-	 * pixels, so we walk *r in chunks which map to the tile size.
-	 *
-	 * Make sure we can't ask for a zero step.
-	 *
-	 * We don't chunk horizontally. We want "vips shrink x.jpg b.jpg 100
-	 * 100" to run sequentially. If we chunk horizontally, we will fetch
-	 * 100x100 lines from the top of the image, then 100x100 100 lines
-	 * down, etc. for each thread, then when they've finished, fetch
-	 * 100x100, 100 pixels across from the top of the image. This will
-	 * break sequentiality. 
-	 */
-	int ystep = shrink->mh > VIPS__TILE_HEIGHT ? 
-		1 : VIPS__TILE_HEIGHT / shrink->mh;
-
-	int y;
-
-#ifdef DEBUG
-	printf( "vips_shrink2_gen: generating %d x %d at %d x %d\n",
-		r->width, r->height, r->left, r->top ); 
-#endif /*DEBUG*/
-
-	for( y = 0; y < r->height; y += ystep ) {
-		/* Clip the this rect against the demand size.
-		 */
-		int height = VIPS_MIN( ystep, r->height - y );
-
-		VipsRect s;
-
-		s.left = r->left * shrink->xshrink;
-		s.top = (r->top + y) * shrink->yshrink;
-		s.width = VIPS_CEIL( r->width * shrink->xshrink );
-		s.height = VIPS_CEIL( height * shrink->yshrink );
-#ifdef DEBUG
-		printf( "shrink_gen: requesting %d x %d at %d x %d\n",
-			s.width, s.height, s.left, s.top ); 
-#endif /*DEBUG*/
-		if( vips_region_prepare( ir, &s ) )
-			return( -1 );
-
-		VIPS_GATE_START( "vips_shrink2_gen: work" ); 
-
-		vips_shrink2_gen2( shrink, seq, 
-			or, ir, 
-			r->left, r->top + y, r->width, height );
-
-		VIPS_GATE_STOP( "vips_shrink2_gen: work" ); 
-	}
-
-	return( 0 );
-}
-
-static int
-vips_shrink2_build( VipsObject *object )
-{
-	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
-	VipsResample *resample = VIPS_RESAMPLE( object );
-	VipsShrink2 *shrink = (VipsShrink2 *) object;
-	VipsImage **t = (VipsImage **) 
-		vips_object_local_array( object, 1 );
-
-	VipsImage *in;
-
-	if( VIPS_OBJECT_CLASS( vips_shrink2_parent_class )->build( object ) )
-		return( -1 );
-
-	shrink->mw = VIPS_CEIL( shrink->xshrink );
-	shrink->mh = VIPS_CEIL( shrink->yshrink );
-	shrink->np = shrink->mw * shrink->mh;
-
-	in = resample->in; 
-
-	if( vips_check_noncomplex( class->nickname, in ) )
-		return( -1 );
-
-	if( shrink->xshrink < 1.0 || 
-		shrink->yshrink < 1.0 ) {
-		vips_error( class->nickname, 
-			"%s", _( "shrink factors should be >= 1" ) );
-		return( -1 );
-	}
-
-	if( (int) shrink->xshrink != shrink->xshrink || 
-		(int) shrink->yshrink != shrink->yshrink ) 
-		vips_warn( class->nickname, 
-			"%s", _( "not integer shrink factors, "
-				"expect poor results" ) ); 
-
-	if( shrink->xshrink == 1.0 &&
-		shrink->yshrink == 1.0 )
-		return( vips_image_write( in, resample->out ) );
-
-	/* Unpack for processing.
-	 */
-	if( vips_image_decode( in, &t[0] ) )
-		return( -1 );
-	in = t[0];
-
-	/* THINSTRIP will work, anything else will break seq mode. If you 
-	 * combine shrink with conv you'll need to use a line cache to maintain
-	 * sequentiality.
-	 */
-	if( vips_image_pipelinev( resample->out, 
-		VIPS_DEMAND_STYLE_THINSTRIP, in, NULL ) )
-		return( -1 );
-
-	/* Size output. Note: we round the output width down!
-	 *
-	 * Don't change xres/yres, leave that to the application layer. For
-	 * example, vipsthumbnail knows the true shrink factor (including the
-	 * fractional part), we just see the integer part here.
-	 */
-	resample->out->Xsize = in->Xsize / shrink->xshrink;
-	resample->out->Ysize = in->Ysize / shrink->yshrink;
-	if( resample->out->Xsize <= 0 || 
-		resample->out->Ysize <= 0 ) {
-		vips_error( class->nickname, 
-			"%s", _( "image has shrunk to nothing" ) );
-		return( -1 );
-	}
-
-#ifdef DEBUG
-	printf( "vips_shrink2_build: shrinking %d x %d image to %d x %d\n", 
-		in->Xsize, in->Ysize, 
-		resample->out->Xsize, resample->out->Ysize );  
-	printf( "vips_shrink2_build: %d x %d block average\n", 
-		shrink->mw, shrink->mh ); 
-#endif /*DEBUG*/
-
-	if( vips_image_generate( resample->out,
-		vips_shrink2_start, vips_shrink2_gen, vips_shrink2_stop, 
-		in, shrink ) )
-		return( -1 );
-
-	return( 0 );
-}
-
-static void
-vips_shrink2_class_init( VipsShrink2Class *class )
-{
-	GObjectClass *gobject_class = G_OBJECT_CLASS( class );
-	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class );
-	VipsOperationClass *operation_class = VIPS_OPERATION_CLASS( class );
-
-	VIPS_DEBUG_MSG( "vips_shrink2_class_init\n" );
-
-	gobject_class->set_property = vips_object_set_property;
-	gobject_class->get_property = vips_object_get_property;
-
-	vobject_class->nickname = "shrink2";
-	vobject_class->description = _( "shrink an image" );
-	vobject_class->build = vips_shrink2_build;
-
-	operation_class->flags = VIPS_OPERATION_SEQUENTIAL;
-
-	VIPS_ARG_DOUBLE( class, "xshrink", 8, 
-		_( "Xshrink" ), 
-		_( "Horizontal shrink factor" ),
-		VIPS_ARGUMENT_REQUIRED_INPUT,
-		G_STRUCT_OFFSET( VipsShrink2, xshrink ),
-		1.0, 1000000, 1 );
-
-	VIPS_ARG_DOUBLE( class, "yshrink", 9, 
-		_( "Yshrink" ), 
-		_( "Vertical shrink factor" ),
-		VIPS_ARGUMENT_REQUIRED_INPUT,
-		G_STRUCT_OFFSET( VipsShrink2, yshrink ),
-		1.0, 1000000, 1 );
-
-}
-
-static void
-vips_shrink2_init( VipsShrink2 *shrink )
-{
-}
-
-/**
- * vips_shrink2:
- * @in: input image
- * @out: output image
- * @xshrink: horizontal shrink
- * @yshrink: vertical shrink
- * @...: %NULL-terminated list of optional named arguments
- *
- * Shrink @in by a pair of factors with a simple box filter. 
- *
- * This is an old version of vips_shrink() kept around for testing. Use
- * vips_shrink() in preference. 
- *
- * See also: vips_resize(), vips_affine().
- *
- * Returns: 0 on success, -1 on error
- */
-int
-vips_shrink2( VipsImage *in, VipsImage **out, 
-	double xshrink, double yshrink, ... )
-{
-	va_list ap;
-	int result;
-
-	va_start( ap, yshrink );
-	result = vips_call_split( "shrink2", ap, in, out, xshrink, yshrink );
-	va_end( ap );
-
-	return( result );
-}

libvips/resample/similarity.c

@@ -8,6 +8,8 @@
  * 	  rest of vips
  * 13/8/14
  * 	- oops, missing scale from b, thanks Topochicho
+ * 7/2/16
+ * 	- use vips_reduce(), if we can
  */
 
 /*
@@ -47,6 +49,7 @@
 #include <vips/intl.h>
 
 #include <math.h>
+#include <string.h>
 
 #include <vips/vips.h>
 
@@ -78,15 +81,31 @@ vips_similarity_build( VipsObject *object )
 	VipsImage **t = (VipsImage **) 
 		vips_object_local_array( object, 4 );
 
-	double a, b, c, d; 
-
 	if( VIPS_OBJECT_CLASS( vips_similarity_parent_class )->build( object ) )
 		return( -1 );
 
-	a = similarity->scale * cos( VIPS_RAD( similarity->angle ) ); 
+	/* Use vips_reduce(), if we can.
-	b = similarity->scale * -sin( VIPS_RAD( similarity->angle ) );
+	 */
-	c = -b;
+	if( similarity->interpolate && 
-	d = a;
+		strcmp( VIPS_OBJECT_GET_CLASS( similarity->interpolate )->
+			nickname, "bicubic" ) == 0 &&
+		similarity->angle == 0.0 &&
+		similarity->idx == 0.0 &&
+		similarity->idy == 0.0 &&
+		similarity->odx == 0.0 &&
+		similarity->ody == 0.0 ) {
+		if( vips_reduce( resample->in, &t[0], 
+			1.0 / similarity->scale, 
+			1.0 / similarity->scale, NULL ) )  
+			return( -1 );
+	}
+	else {
+		double a = similarity->scale * 
+			cos( VIPS_RAD( similarity->angle ) ); 
+		double b = similarity->scale * 
+			-sin( VIPS_RAD( similarity->angle ) );
+		double c = -b;
+		double d = a;
 
 		if( vips_affine( resample->in, &t[0], a, b, c, d, 
 			"interpolate", similarity->interpolate,
@@ -94,8 +113,11 @@ vips_similarity_build( VipsObject *object )
 			"ody", similarity->ody,
 			"idx", similarity->idx,
 			"idy", similarity->idy,
-		NULL ) ||
+			NULL ) )
-		vips_image_write( t[0], resample->out ) )
+			return( -1 );
+	}
+
+	if( vips_image_write( t[0], resample->out ) )
 		return( -1 ); 
 
 	return( 0 );

libvips/resample/templates.h

@@ -164,29 +164,31 @@ bicubic_unsigned_int(
 	const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
 	const int* restrict cx, const int* restrict cy )
 {
+	const int c0 = cx[0];
+	const int c1 = cx[1];
+	const int c2 = cx[2];
+	const int c3 = cx[3];
+
 	const int r0 = unsigned_fixed_round( 
-		cx[0] * uno_one +
+		c0 * uno_one +
-		cx[1] * uno_two +
+		c1 * uno_two +
-		cx[2] * uno_thr +
+		c2 * uno_thr +
-		cx[3] * uno_fou ); 
+		c3 * uno_fou ); 
-
 	const int r1 = unsigned_fixed_round( 
-		cx[0] * dos_one +
+		c0 * dos_one +
-		cx[1] * dos_two +
+		c1 * dos_two +
-		cx[2] * dos_thr +
+		c2 * dos_thr +
-		cx[3] * dos_fou );
+		c3 * dos_fou ); 
-
 	const int r2 = unsigned_fixed_round( 
-		cx[0] * tre_one +
+		c0 * tre_one +
-		cx[1] * tre_two +
+		c1 * tre_two +
-		cx[2] * tre_thr +
+		c2 * tre_thr +
-		cx[3] * tre_fou );
+		c3 * tre_fou ); 
-
 	const int r3 = unsigned_fixed_round( 
-		cx[0] * qua_one +
+		c0 * qua_one +
-		cx[1] * qua_two +
+		c1 * qua_two +
-		cx[2] * qua_thr +
+		c2 * qua_thr +
-		cx[3] * qua_fou );
+		c3 * qua_fou ); 
 
 	return( unsigned_fixed_round( 
 		cy[0] * r0 +
@@ -214,29 +216,31 @@ bicubic_signed_int(
 	const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
 	const int* restrict cx, const int* restrict cy )
 {
+	const int c0 = cx[0];
+	const int c1 = cx[1];
+	const int c2 = cx[2];
+	const int c3 = cx[3];
+
 	const int r0 = signed_fixed_round( 
-		cx[0] * uno_one +
+		c0 * uno_one +
-		cx[1] * uno_two +
+		c1 * uno_two +
-		cx[2] * uno_thr +
+		c2 * uno_thr +
-		cx[3] * uno_fou ); 
+		c3 * uno_fou ); 
-
 	const int r1 = signed_fixed_round( 
-		cx[0] * dos_one +
+		c0 * dos_one +
-		cx[1] * dos_two +
+		c1 * dos_two +
-		cx[2] * dos_thr +
+		c2 * dos_thr +
-		cx[3] * dos_fou );
+		c3 * dos_fou ); 
-
 	const int r2 = signed_fixed_round( 
-		cx[0] * tre_one +
+		c0 * tre_one +
-		cx[1] * tre_two +
+		c1 * tre_two +
-		cx[2] * tre_thr +
+		c2 * tre_thr +
-		cx[3] * tre_fou );
+		c3 * tre_fou ); 
-
 	const int r3 = signed_fixed_round( 
-		cx[0] * qua_one +
+		c0 * qua_one +
-		cx[1] * qua_two +
+		c1 * qua_two +
-		cx[2] * qua_thr +
+		c2 * qua_thr +
-		cx[3] * qua_fou );
+		c3 * qua_fou ); 
 
 	return( signed_fixed_round( 
 		cy[0] * r0 +
@@ -245,6 +249,17 @@ bicubic_signed_int(
 		cy[3] * r3 ) ); 
 }
 
+template <typename T> static T inline
+cubic_float(
+	const T one, const T two, const T thr, const T fou,
+	const double* restrict cx )
+{
+	return( cx[0] * one +
+		 cx[1] * two +
+		 cx[2] * thr +
+		 cx[3] * fou );
+}
+
 /* Floating-point bicubic, used for int/float/double types.
  */
 template <typename T> static T inline
@@ -255,26 +270,16 @@ bicubic_float(
 	const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
 	const double* restrict cx, const double* restrict cy )
 {
-	return(
+	const double r0 = cubic_float<T>( 
-		cy[0] * (cx[0] * uno_one +
+		uno_one, uno_two, uno_thr, uno_fou, cx ); 
-			 cx[1] * uno_two +
+	const double r1 = cubic_float<T>( 
-			 cx[2] * uno_thr +
+		dos_one, dos_two, dos_thr, dos_fou, cx ); 
-			 cx[3] * uno_fou)
+	const double r2 = cubic_float<T>( 
-                +
+		tre_one, tre_two, tre_thr, tre_fou, cx ); 
-		cy[1] * (cx[0] * dos_one +
+	const double r3 = cubic_float<T>( 
-			 cx[1] * dos_two +
+		qua_one, qua_two, qua_thr, qua_fou, cx ); 
-			 cx[2] * dos_thr +
+
-			 cx[3] * dos_fou)
+	return( cubic_float<T>( r0, r1, r2, r3, cy ) ); 
-                +
-		cy[2] * (cx[0] * tre_one +
-			 cx[1] * tre_two +
-			 cx[2] * tre_thr +
-			 cx[3] * tre_fou)
-                +
-		cy[3] * (cx[0] * qua_one +
-			 cx[1] * qua_two +
-			 cx[2] * qua_thr +
-			 cx[3] * qua_fou) );
 }
 
 /* Given an offset in [0,1] (we can have x == 1 when building tables),

test/test_arithmetic.py

@@ -11,7 +11,6 @@ from gi.repository import Vips
 
 Vips.leak_set(True)
 
-
 unsigned_formats = [Vips.BandFormat.UCHAR, 
                     Vips.BandFormat.USHORT, 
                     Vips.BandFormat.UINT] 

test/test_resample.py

@@ -10,6 +10,20 @@ from gi.repository import Vips
 
 Vips.leak_set(True)
 
+unsigned_formats = [Vips.BandFormat.UCHAR, 
+                    Vips.BandFormat.USHORT, 
+                    Vips.BandFormat.UINT] 
+signed_formats = [Vips.BandFormat.CHAR, 
+                  Vips.BandFormat.SHORT, 
+                  Vips.BandFormat.INT] 
+float_formats = [Vips.BandFormat.FLOAT, 
+                 Vips.BandFormat.DOUBLE]
+complex_formats = [Vips.BandFormat.COMPLEX, 
+                   Vips.BandFormat.DPCOMPLEX] 
+int_formats = unsigned_formats + signed_formats
+noncomplex_formats = int_formats + float_formats
+all_formats = int_formats + float_formats + complex_formats
+
 # Run a function expecting a complex image on a two-band image
 def run_cmplx(fn, image):
     if image.format == Vips.BandFormat.FLOAT:
@@ -100,6 +114,20 @@ class TestResample(unittest.TestCase):
  
             self.assertEqual((x - im).abs().max(), 0)
 
+    def test_reduce(self):
+        im = Vips.Image.new_from_file("images/IMG_4618.jpg")
+        bicubic = Vips.Interpolate.new("bicubic")
+        
+        for fac in [1, 1.1, 1.5, 1.999]:
+            for fmt in all_formats:
+                x = im.cast(fmt)
+                r = x.reduce(fac, fac)
+                a = x.affine([1.0 / fac, 0, 0, 1.0 / fac], 
+                             interpolate = bicubic,
+                             oarea = [0, 0, x.width / fac, x.height / fac])
+                d = (r - a).abs().max()
+                self.assertLess(d, 5)
+
     def test_resize(self):
         im = Vips.Image.new_from_file("images/IMG_4618.jpg")
         im2 = im.resize(0.25)