diff --git a/libvips/resample/bicubic.cpp b/libvips/resample/bicubic.cpp
index 7388ab06..97e4a403 100644
--- a/libvips/resample/bicubic.cpp
+++ b/libvips/resample/bicubic.cpp
@@ -79,15 +79,16 @@ typedef VipsInterpolate VipsInterpolateBicubic;
 typedef VipsInterpolateClass VipsInterpolateBicubicClass;
 
 /* Precalculated interpolation matrices. int (used for pel
- * sizes up to short), and double (for all others).
+ * sizes up to short), and double (for all others). We go to
+ * scale + 1 so we can round-to-nearest safely.
  */
 
 /* We could keep a large set of 2d 4x4 matricies, but this actually
  * works out slower since for many resizes the thing will no longer
  * fit in L1.
  */
-static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE][4];
-static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE][4];
+static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE + 1][4];
+static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE + 1][4];
 
 /* We need C linkage for this.
  */
@@ -497,13 +498,19 @@ static void
 vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
 	void *out, VipsRegion *in, double x, double y )
 {
-	/* Find the mask index.
+	/* Find the mask index. We round-to-nearest, so we need to generate 
+	 * indexes in 0 to VIPS_TRANSFORM_SCALE, 2^n + 1 values. We multiply 
+	 * by 2 more than we need to, add one, mask, then shift down again to 
+	 * get the extra range.
 	 */
-	const int sx = x * VIPS_TRANSFORM_SCALE;
-	const int sy = y * VIPS_TRANSFORM_SCALE;
+	const int sx = x * VIPS_TRANSFORM_SCALE * 2;
+	const int sy = y * VIPS_TRANSFORM_SCALE * 2;
 
-	const int tx = sx & (VIPS_TRANSFORM_SCALE - 1);
-	const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
+	const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
+	const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
+
+	const int tx = (six + 1) >> 1;
+	const int ty = (siy + 1) >> 1;
 
 	/* We know x/y are always positive, so we can just (int) them. 
 	 */
@@ -636,7 +643,7 @@ vips_interpolate_bicubic_class_init( VipsInterpolateBicubicClass *iclass )
 
 	/* Build the tables of pre-computed coefficients.
 	 */
-	for( int x = 0; x < VIPS_TRANSFORM_SCALE; x++ ) {
+	for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
 		calculate_coefficients_catmull( vips_bicubic_matrixf[x], 
 			(float) x / VIPS_TRANSFORM_SCALE ); 
 
diff --git a/libvips/resample/reduceh.cpp b/libvips/resample/reduceh.cpp
index cc829c01..3a88f4ef 100644
--- a/libvips/resample/reduceh.cpp
+++ b/libvips/resample/reduceh.cpp
@@ -11,7 +11,6 @@
  * 6/6/20 kleisauke
  * 	- deprecate @centre option, it's now always on
  * 	- fix pixel shift
- * 	- remove unnecessary round-to-nearest behaviour
  */
 
 /*
@@ -79,10 +78,11 @@ typedef struct _VipsReduceh {
 	double hoffset;
 
 	/* Precalculated interpolation matrices. int (used for pel
-	 * sizes up to short), and double (for all others).
+	 * sizes up to short), and double (for all others). We go to
+	 * scale + 1 so we can round-to-nearest safely.
 	 */
-	int *matrixi[VIPS_TRANSFORM_SCALE];
-	double *matrixf[VIPS_TRANSFORM_SCALE];
+	int *matrixi[VIPS_TRANSFORM_SCALE + 1];
+	double *matrixf[VIPS_TRANSFORM_SCALE + 1];
 
 	/* Deprecated.
 	 */
@@ -346,8 +346,9 @@ vips_reduceh_gen( VipsRegion *out_region, void *seq,
 		for( int x = 0; x < r->width; x++ ) {
 			const int ix = (int) X;
 			VipsPel *p = p0 + ix * ps;
-			const int sx = X * VIPS_TRANSFORM_SCALE;
-			const int tx = sx & (VIPS_TRANSFORM_SCALE - 1);
+			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 = reduceh->matrixi[tx];
 			const double *cxf = reduceh->matrixf[tx];
 
@@ -479,7 +480,7 @@ vips_reduceh_build( VipsObject *object )
 
 	/* Build the tables of pre-computed coefficients.
 	 */
-	for( int x = 0; x < VIPS_TRANSFORM_SCALE; x++ ) {
+	for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
 		reduceh->matrixf[x] = 
 			VIPS_ARRAY( object, reduceh->n_point, double ); 
 		reduceh->matrixi[x] = 
diff --git a/libvips/resample/reducev.cpp b/libvips/resample/reducev.cpp
index 396069f2..847f79ec 100644
--- a/libvips/resample/reducev.cpp
+++ b/libvips/resample/reducev.cpp
@@ -21,7 +21,6 @@
  * 	- deprecate @centre option, it's now always on
  * 	- fix pixel shift
  * 	- speed up the mask construction for uchar/ushort images
- * 	- remove unnecessary round-to-nearest behaviour
  */
 
 /*
@@ -118,14 +117,15 @@ typedef struct _VipsReducev {
 	double voffset;
 
 	/* Precalculated interpolation matrices. int (used for pel
-	 * sizes up to short), and double (for all others).
+	 * sizes up to short), and double (for all others). We go to
+	 * scale + 1 so we can round-to-nearest safely.
 	 */
-	int *matrixi[VIPS_TRANSFORM_SCALE];
-	double *matrixf[VIPS_TRANSFORM_SCALE];
+	int *matrixi[VIPS_TRANSFORM_SCALE + 1];
+	double *matrixf[VIPS_TRANSFORM_SCALE + 1];
 
 	/* And another set for orc: we want 2.6 precision.
 	 */
-	int *matrixo[VIPS_TRANSFORM_SCALE];
+	int *matrixo[VIPS_TRANSFORM_SCALE + 1];
 
 	/* The passes we generate for this mask.
 	 */
@@ -154,7 +154,7 @@ vips_reducev_finalize( GObject *gobject )
 	for( int i = 0; i < reducev->n_pass; i++ )
 		VIPS_FREEF( vips_vector_free, reducev->pass[i].vector );
 	reducev->n_pass = 0;
-	for( int i = 0; i < VIPS_TRANSFORM_SCALE; i++ ) {
+	for( int i = 0; i < VIPS_TRANSFORM_SCALE + 1; i++ ) {
 		VIPS_FREE( reducev->matrixf[i] );
 		VIPS_FREE( reducev->matrixi[i] );
 		VIPS_FREE( reducev->matrixo[i] );
@@ -555,8 +555,9 @@ vips_reducev_gen( VipsRegion *out_region, void *vseq,
 			VIPS_REGION_ADDR( out_region, r->left, r->top + y );
 		const int py = (int) Y;
 		VipsPel *p = VIPS_REGION_ADDR( ir, r->left, py );
-		const int sy = Y * VIPS_TRANSFORM_SCALE;
-		const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
+		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 = reducev->matrixi[ty];
 		const double *cyf = reducev->matrixf[ty];
 		const int lskip = VIPS_REGION_LSKIP( ir );
@@ -677,8 +678,9 @@ vips_reducev_vector_gen( VipsRegion *out_region, void *vseq,
 		VipsPel *q = 
 			VIPS_REGION_ADDR( out_region, r->left, r->top + y );
 		const int py = (int) Y;
-		const int sy = Y * VIPS_TRANSFORM_SCALE;
-		const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
+		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 *cyo = reducev->matrixo[ty];
 
 #ifdef DEBUG_PIXELS
@@ -739,7 +741,7 @@ vips_reducev_raw( VipsReducev *reducev, VipsImage *in, VipsImage **out )
 	 */
 	if( in->BandFmt == VIPS_FORMAT_UCHAR &&
 		vips_vector_isenabled() ) 
-		for( int y = 0; y < VIPS_TRANSFORM_SCALE; y++ ) {
+		for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
 			reducev->matrixo[y] = 
 				VIPS_ARRAY( NULL, reducev->n_point, int ); 
 			if( !reducev->matrixo[y] )
@@ -851,7 +853,7 @@ vips_reducev_build( VipsObject *object )
 
 	/* Build the tables of pre-computed coefficients.
 	 */
-	for( int y = 0; y < VIPS_TRANSFORM_SCALE; y++ ) {
+	for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
 		reducev->matrixf[y] = 
 			VIPS_ARRAY( NULL, reducev->n_point, double ); 
 		reducev->matrixi[y] =