Improve edge handling for mapim (#2681)

* almost there

* smol changes

* getting close now

segv on the test suite, annoyingly

* argh still not quite thert

* phew

ChangeLog

@@ -1,14 +1,17 @@
 21/11/21 started 8.13
 - configure fails for requested but unmet dependencies [remicollet]
 - add support for another quantiser [DarthSim]
+- add "extend", "background" and "premultiplied" to mapim to fix edge
+  antialiasing [GavinJoyce]
 - add support for HDR HEIC and AVIF images
 - add spngsave
 - jpeg2000 load left-justifies bitdepth
 - add "password" option to pdfload
 - improve the pixel rng
-- add meson build system [tintou]
+- new meson build system [tintou]
 - improve introspection annotations [tintou]
 - add @unlimited to heifload [lovell]
+- edge antialiasing for mapim with @background, @extend and @premultiplied
 - add support for regions in C++ API [shado23]
 
 26/11/21 started 8.12.3

libvips/iofuncs/type.c

@@ -700,7 +700,7 @@ vips_blob_get( VipsBlob *blob, size_t *length )
 
 /* vips_blob_set:
  * @blob: #VipsBlob to set
- * @free_fn: (scope async) (nullable): @data will be freed with this function
+ * @free_fn: (scope async) (allow-none): @data will be freed with this function
  * @data: (array length=length) (element-type guint8) (transfer full): data to store
  * @length: number of bytes in @data
  *

libvips/resample/affine.c

@@ -374,8 +374,7 @@ vips_affine_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 					(int) iy - window_offset + 
 						window_size - 1 ) );
 
-				interpolate( affine->interpolate, 
+				interpolate( affine->interpolate, q, ir, ix, iy );
-					q, ir, ix, iy );
 			}
 			else {
 				/* Out of range: paint the background.

libvips/resample/mapim.c

@@ -9,6 +9,9 @@
  * 	- we were not offsetting pixel fetches by window_offset
  * 30/1/21 afontenot
  * 	- avoid NaN
+ * 21/12/21
+ * 	- improve edge antialiasing with "background" and "extend"
+ * 	- add "premultiplied" param
  */
 
 /*
@@ -67,6 +70,22 @@ typedef struct _VipsMapim {
 	VipsImage *index;
 	VipsInterpolate *interpolate;
 
+	/* How to generate extra edge pixels.
+	 */
+	VipsExtend extend;
+
+	/* Background colour.
+	 */
+	VipsArrayDouble *background;
+
+	/* The [double] converted to the input image format.
+	 */
+	VipsPel *ink;
+
+	/* True if the input is already premultiplied (and we don't need to).
+	 */
+	gboolean premultiplied;
+
 	/* Need an image vector for start_many / stop_many
 	 */
 	VipsImage *in_array[3];
@@ -77,11 +96,11 @@ typedef VipsResampleClass VipsMapimClass;
 
 G_DEFINE_TYPE( VipsMapim, vips_mapim, VIPS_TYPE_RESAMPLE );
 
-/* Minmax of a line of pixels. Pass in a thing to convert back to int 
+/* Minmax of a line of pixels. 
- * coordinates.
  */
-#define MINMAX( TYPE, CLIP ) { \
+#define MINMAX( TYPE ) { \
 	TYPE * restrict p1 = (TYPE *) p; \
+	\
 	TYPE t_max_x = max_x; \
 	TYPE t_min_x = min_x; \
 	TYPE t_max_y = max_y; \
@@ -114,94 +133,28 @@ G_DEFINE_TYPE( VipsMapim, vips_mapim, VIPS_TYPE_RESAMPLE );
 		p1 += 2; \
 	} \
 	\
-	min_x = CLIP( t_min_x ); \
+	min_x = t_min_x; \
-	max_x = CLIP( t_max_x ); \
+	max_x = t_max_x; \
-	min_y = CLIP( t_min_y ); \
+	min_y = t_min_y; \
-	max_y = CLIP( t_max_y ); \
+	max_y = t_max_y; \
 }
 
-/* Minmax of a line of float pixels. We have to ignore NaN. 
- */
-#define FMINMAX( TYPE ) { \
-	TYPE * restrict p1 = (TYPE *) p; \
-	TYPE t_max_x = max_x; \
-	TYPE t_min_x = min_x; \
-	TYPE t_max_y = max_y; \
-	TYPE t_min_y = min_y; \
-	\
-	for( x = 0; x < r->width; x++ ) { \
-		TYPE px = p1[0]; \
-		TYPE py = p1[1]; \
-		\
-		if( !VIPS_ISNAN( px ) && \
-			!VIPS_ISNAN( py ) ) { \
-			if( first ) { \
-				t_min_x = px; \
-				t_max_x = px; \
-				t_min_y = py; \
-				t_max_y = py; \
-				\
-				first = FALSE; \
-			} \
-			else { \
-				if( px > t_max_x ) \
-					t_max_x = px; \
-				else if( px < t_min_x ) \
-					t_min_x = px; \
-				\
-				if( py > t_max_y ) \
-					t_max_y = py; \
-				else if( py < t_min_y ) \
-					t_min_y = py; \
-			} \
-		}\
-		\
-		p1 += 2; \
-	} \
-	\
-	if( !first ) { \
-		min_x = VIPS_CLIP( 0, floor( t_min_x ), VIPS_MAX_COORD ); \
-		max_x = VIPS_CLIP( 0, floor( t_max_x ), VIPS_MAX_COORD ); \
-		min_y = VIPS_CLIP( 0, floor( t_min_y ), VIPS_MAX_COORD ); \
-		max_y = VIPS_CLIP( 0, floor( t_max_y ), VIPS_MAX_COORD ); \
-	} \
-}
-
-/* All the clippers. These vary with TYPE.
- */
-
-/* Clip a small (max val < VIPS_MAX_COORD) unsigned int type.
- */
-#define CLIP_UINT_SMALL( X ) (X)
-
-/* Clip a small (max val < VIPS_MAX_COORD) signed int type.
- */
-#define CLIP_SINT_SMALL( X ) VIPS_MAX( X, 0 );
-
-/* An unsigned int type larger than VIPS_MAX_COORD. Trim upper range.
- */
-#define CLIP_UINT_LARGE( X ) VIPS_MIN( X, VIPS_MAX_COORD ); 
-
-/* A large signed int.
- */
-#define CLIP_SINT_LARGE( X ) VIPS_CLIP( 0, X, VIPS_MAX_COORD );
-
 /* Scan a region and find min/max in the two axes.
  */
 static void
 vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 {
-	int min_x;
+	double min_x;
-	int max_x;
+	double max_x;
-	int min_y;
+	double min_y;
-	int max_y;
+	double max_y;
 	gboolean first;
 	int x, y;
 
-	min_x = 0;
+	min_x = 0.0;
-	max_x = 0;
+	max_x = 0.0;
-	min_y = 0;
+	min_y = 0.0;
-	max_y = 0;
+	max_y = 0.0;
 	first = TRUE;
 	for( y = 0; y < r->height; y++ ) {
 		VipsPel * restrict p = 
@@ -209,37 +162,37 @@ vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 
 		switch( region->im->BandFmt ) {
 		case VIPS_FORMAT_UCHAR: 	
-			MINMAX( unsigned char, CLIP_UINT_SMALL ); 
+			MINMAX( unsigned char );
 			break; 
 
 		case VIPS_FORMAT_CHAR: 	
-			MINMAX( signed char, CLIP_SINT_SMALL ); 
+			MINMAX( signed char );
 			break; 
 
 		case VIPS_FORMAT_USHORT: 
-			MINMAX( unsigned short, CLIP_UINT_SMALL ); 
+			MINMAX( unsigned short );
 			break; 
 
 		case VIPS_FORMAT_SHORT: 	
-			MINMAX( signed short, CLIP_SINT_SMALL ); 
+			MINMAX( signed short );
 			break; 
 
 		case VIPS_FORMAT_UINT: 	
-			MINMAX( unsigned int, CLIP_UINT_LARGE ); 
+			MINMAX( unsigned int );
 			break; 
 
 		case VIPS_FORMAT_INT: 	
-			MINMAX( signed int, CLIP_SINT_LARGE ); 
+			MINMAX( signed int );
 			break; 
 
 		case VIPS_FORMAT_FLOAT: 		
 		case VIPS_FORMAT_COMPLEX: 
-			FMINMAX( float );
+			MINMAX( float );
 			break; 
 
 		case VIPS_FORMAT_DOUBLE:	
 		case VIPS_FORMAT_DPCOMPLEX: 
-			FMINMAX( double );
+			MINMAX( double );
 			break;
 
 		default:
@@ -247,10 +200,26 @@ vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 		}
 	}
 
+	/* bounds is the bounding box -- we must round left/top down and round
+	 * bottom/right up.
+	 */
+	min_x = floor( min_x );
+	min_y = floor( min_y );
+	max_x = ceil( max_x );
+	max_y = ceil( max_y );
+
+	/* bounds uses ints, so we must clip the range down from double.
+	 * Coordinates can be negative for the antialias edges.
+	 */
+	min_x = VIPS_CLIP( -1, min_x, VIPS_MAX_COORD );
+	min_y = VIPS_CLIP( -1, min_y, VIPS_MAX_COORD );
+	max_x = VIPS_CLIP( -1, max_x, VIPS_MAX_COORD );
+	max_y = VIPS_CLIP( -1, max_y, VIPS_MAX_COORD );
+
 	bounds->left = min_x;
 	bounds->top = min_y;
-	bounds->width = max_x - min_x + 1;
+	bounds->width = (max_x - min_x) + 1;
-	bounds->height = max_y - min_y + 1;
+	bounds->height = (max_y - min_y) + 1;
 }
 
 /* Unsigned int types.
@@ -265,18 +234,19 @@ vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 		if( px >= clip_width || \
 			py >= clip_height ) { \
 			for( z = 0; z < ps; z++ )  \
-				q[z] = 0; \
+				q[z] = mapim->ink[z]; \
 		} \
 		else \
 			interpolate( mapim->interpolate, q, ir[0], \
-				px + window_offset, py + window_offset ); \
+				px + window_offset + 1, \
+				py + window_offset + 1 ); \
 		\
 		p1 += 2; \
 		q += ps; \
 	} \
 }
 
-/* Signed int types.
+/* Signed int types. We allow -1 for x/y to get edge antialiasing.
  */
 #define LOOKUP( TYPE ) { \
 	TYPE * restrict p1 = (TYPE *) p; \
@@ -285,23 +255,24 @@ vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 		TYPE px = p1[0]; \
 		TYPE py = p1[1]; \
 		\
-		if( px < 0 || \
+		if( px < -1 || \
 			px >= clip_width || \
-			py < 0 || \
+			py < -1 || \
 			py >= clip_height ) { \
 			for( z = 0; z < ps; z++ )  \
-				q[z] = 0; \
+				q[z] = mapim->ink[z]; \
 		} \
 		else \
 			interpolate( mapim->interpolate, q, ir[0], \
-				px + window_offset, py + window_offset ); \
+				px + window_offset + 1, \
+				py + window_offset + 1 ); \
 		\
 		p1 += 2; \
 		q += ps; \
 	} \
 }
 
-/* Float types.
+/* Float types. We allow -1 for x/y to get edge antialiasing.
  */
 #define FLOOKUP( TYPE ) { \
 	TYPE * restrict p1 = (TYPE *) p; \
@@ -312,16 +283,17 @@ vips_mapim_region_minmax( VipsRegion *region, VipsRect *r, VipsRect *bounds )
 		\
 		if( VIPS_ISNAN( px ) || \
 			VIPS_ISNAN( py ) || \
-			px < 0 || \
+			px < -1 || \
 			px >= clip_width || \
-			py < 0 || \
+			py < -1 || \
 			py >= clip_height ) { \
-			for( z = 0; z < ps; z++ )  \
+			for( z = 0; z < ps; z++ ) \
-				q[z] = 0; \
+				q[z] = mapim->ink[z]; \
 		} \
 		else \
 			interpolate( mapim->interpolate, q, ir[0], \
-				px + window_offset, py + window_offset ); \
+				px + window_offset + 1, \
+				py + window_offset + 1 ); \
 		\
 		p1 += 2; \
 		q += ps; \
@@ -335,7 +307,6 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 	VipsRegion **ir = (VipsRegion **) seq;
 	const VipsImage **in_array = (const VipsImage **) a;
 	const VipsMapim *mapim = (VipsMapim *) b; 
-	const VipsResample *resample = VIPS_RESAMPLE( mapim );
 	const VipsImage *in = in_array[0];
 	const int window_size = 
 		vips_interpolate_get_window_size( mapim->interpolate );
@@ -344,10 +315,10 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 	const VipsInterpolateMethod interpolate = 
 		vips_interpolate_get_method( mapim->interpolate );
 	const int ps = VIPS_IMAGE_SIZEOF_PEL( in );
-	const int clip_width = resample->in->Xsize;
+	const int clip_width = in->Xsize - window_size;
-	const int clip_height = resample->in->Ysize;
+	const int clip_height = in->Ysize - window_size;
 
-	VipsRect bounds, image, clipped;
+	VipsRect bounds, need, image, clipped;
 	int x, y, z;
 	
 #ifdef DEBUG_VERBOSE
@@ -359,7 +330,7 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 		r->height );
 #endif /*DEBUG_VERBOSE*/
 
-	/* Fetch the chunk of the mapim image we need, and find the max/min in
+	/* Fetch the chunk of the index image we need, and find the max/min in
 	 * x and y.
 	 */
 	if( vips_region_prepare( ir[1], r ) )
@@ -371,11 +342,15 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 
 	VIPS_GATE_STOP( "vips_mapim_gen: work" ); 
 
-	/* The bounding box of that area is what we will need from @in. Add
+	/* Enlarge by the stencil size.
-	 * enough for the interpolation stencil as well.
 	 */
-	bounds.width += window_size - 1;
+	need.width = bounds.width + window_size - 1;
-	bounds.height += window_size - 1;
+	need.height = bounds.height + window_size - 1;
+
+	/* Offset for the antialias edge we have top and left.
+	 */
+	need.left = bounds.left + 1;
+	need.top = bounds.top + 1;
 
 	/* Clip against the expanded image.
 	 */
@@ -383,7 +358,7 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 	image.top = 0;
 	image.width = in->Xsize;
 	image.height = in->Ysize;
-	vips_rect_intersectrect( &bounds, &image, &clipped );
+	vips_rect_intersectrect( &need, &image, &clipped );
 
 #ifdef DEBUG_VERBOSE
 	printf( "vips_mapim_gen: "
@@ -395,7 +370,7 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 #endif /*DEBUG_VERBOSE*/
 
 	if( vips_rect_isempty( &clipped ) ) {
-		vips_region_black( or );
+		vips_region_paint_pel( or, r, mapim->ink );
 		return( 0 );
 	}
 	if( vips_region_prepare( ir[0], &clipped ) )
@@ -403,7 +378,7 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 
 	VIPS_GATE_START( "vips_mapim_gen: work" ); 
 
-	/* Resample! x/y loop over pixels in the output image (5).
+	/* Resample! x/y loop over pixels in the output (and index) images.
 	 */
 	for( y = 0; y < r->height; y++ ) {
 		VipsPel * restrict p = 
@@ -424,12 +399,9 @@ vips_mapim_gen( VipsRegion *or, void *seq, void *a, void *b, gboolean *stop )
 			ULOOKUP( unsigned int ); break; 
 		case VIPS_FORMAT_INT: 	
 			LOOKUP( signed int ); break; 
-
 		case VIPS_FORMAT_FLOAT: 		
 		case VIPS_FORMAT_COMPLEX: 
 			FLOOKUP( float ); break; 
-			break;
-
 		case VIPS_FORMAT_DOUBLE:	
 		case VIPS_FORMAT_DPCOMPLEX: 
 			FLOOKUP( double ); break;
@@ -450,12 +422,17 @@ vips_mapim_build( VipsObject *object )
 	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
 	VipsResample *resample = VIPS_RESAMPLE( object );
 	VipsMapim *mapim = (VipsMapim *) object;
-	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 6 );
 
 	VipsImage *in;
 	int window_size;
 	int window_offset;
 
+	/* TRUE if we've premultiplied and need to unpremultiply.
+	 */
+	gboolean have_premultiplied;
+	VipsBandFormat unpremultiplied_format;
+
 	if( VIPS_OBJECT_CLASS( vips_mapim_parent_class )->build( object ) )
 		return( -1 );
 
@@ -474,30 +451,77 @@ vips_mapim_build( VipsObject *object )
 		vips_interpolate_get_window_offset( mapim->interpolate );
 
 	/* Add new pixels around the input so we can interpolate at the edges.
+	 *
+	 * We add the interpolate stencil, plus one extra pixel on all the
+	 * edges. This means when we clip in generate (above) we can be sure 
+	 * we clip outside the real pixels and don't get jaggies on edges.
+	 *
+	 * We allow for the +1 in the adjustment to window_offset in generate.
 	 */
 	if( vips_embed( in, &t[1], 
-		window_offset, window_offset, 
+		window_offset + 1, window_offset + 1, 
-		in->Xsize + window_size - 1, in->Ysize + window_size - 1,
+		in->Xsize + window_size - 1 + 2, 
-		"extend", VIPS_EXTEND_COPY,
+		in->Ysize + window_size - 1 + 2,
+		"extend", mapim->extend,
+		"background", mapim->background,
 		NULL ) )
 		return( -1 );
 	in = t[1];
 
-	if( vips_image_pipelinev( resample->out, VIPS_DEMAND_STYLE_SMALLTILE, 
+	/* If there's an alpha and we've not premultiplied, we have to 
+	 * premultiply before resampling. 
+	 */
+	have_premultiplied = FALSE;
+	if( vips_image_hasalpha( in ) &&
+		!mapim->premultiplied ) { 
+		if( vips_premultiply( in, &t[2], NULL ) ) 
+			return( -1 );
+		have_premultiplied = TRUE;
+
+		/* vips_premultiply() makes a float image. When we
+		 * vips_unpremultiply() below, we need to cast back to the
+		 * pre-premultiply format.
+		 */
+		unpremultiplied_format = in->BandFmt;
+		in = t[2];
+	}
+
+	/* Convert the background to the image's format.
+	 */
+	if( !(mapim->ink = vips__vector_to_ink( class->nickname, 
+		in,
+		VIPS_AREA( mapim->background )->data, NULL, 
+		VIPS_AREA( mapim->background )->n )) )
+		return( -1 );
+
+	t[3] = vips_image_new();
+	if( vips_image_pipelinev( t[3], VIPS_DEMAND_STYLE_SMALLTILE, 
 		in, NULL ) )
 		return( -1 );
 
-	resample->out->Xsize = mapim->index->Xsize;
+	t[3]->Xsize = mapim->index->Xsize;
-	resample->out->Ysize = mapim->index->Ysize;
+	t[3]->Ysize = mapim->index->Ysize;
 
 	mapim->in_array[0] = in;
 	mapim->in_array[1] = mapim->index;
 	mapim->in_array[2] = NULL;
-	if( vips_image_generate( resample->out, 
+	if( vips_image_generate( t[3], 
 		vips_start_many, vips_mapim_gen, vips_stop_many, 
 		mapim->in_array, mapim ) )
 		return( -1 );
 
+	in = t[3];
+
+	if( have_premultiplied ) {
+		if( vips_unpremultiply( in, &t[4], NULL ) || 
+			vips_cast( t[4], &t[5], unpremultiplied_format, NULL ) )
+			return( -1 );
+		in = t[5];
+	}
+
+	if( vips_image_write( in, resample->out ) )
+		return( -1 );
+
 	return( 0 );
 }
 
@@ -528,12 +552,35 @@ vips_mapim_class_init( VipsMapimClass *class )
 		VIPS_ARGUMENT_OPTIONAL_INPUT, 
 		G_STRUCT_OFFSET( VipsMapim, interpolate ) );
 
+	VIPS_ARG_ENUM( class, "extend", 117, 
+		_( "Extend" ), 
+		_( "How to generate the extra pixels" ),
+		VIPS_ARGUMENT_OPTIONAL_INPUT,
+		G_STRUCT_OFFSET( VipsMapim, extend ),
+		VIPS_TYPE_EXTEND, VIPS_EXTEND_BACKGROUND );
+
+	VIPS_ARG_BOXED( class, "background", 116, 
+		_( "Background" ), 
+		_( "Background value" ),
+		VIPS_ARGUMENT_OPTIONAL_INPUT,
+		G_STRUCT_OFFSET( VipsMapim, background ),
+		VIPS_TYPE_ARRAY_DOUBLE );
+
+	VIPS_ARG_BOOL( class, "premultiplied", 117,
+		_( "Premultiplied" ),
+		_( "Images have premultiplied alpha" ),
+		VIPS_ARGUMENT_OPTIONAL_INPUT,
+		G_STRUCT_OFFSET( VipsMapim, premultiplied ),
+		FALSE );
+
 }
 
 static void
 vips_mapim_init( VipsMapim *mapim )
 {
 	mapim->interpolate = vips_interpolate_new( "bilinear" );
+	mapim->extend = VIPS_EXTEND_BACKGROUND;
+	mapim->background = vips_array_double_newv( 1, 0.0 );
 }
 
 /**
@@ -546,6 +593,9 @@ vips_mapim_init( VipsMapim *mapim )
  * Optional arguments:
  *
  * * @interpolate: interpolate pixels with this
+ * * @extend: #VipsExtend how to generate new pixels 
+ * * @background: #VipsArrayDouble colour for new pixels 
+ * * @premultiplied: %gboolean, images are already premultiplied
  *
  * This operator resamples @in using @index to look up pixels. @out is
  * the same size as @index, with each pixel being fetched from that position in
@@ -557,12 +607,21 @@ vips_mapim_init( VipsMapim *mapim )
  *
  * If @index has one band, that band must be complex. Otherwise, @index must
  * have two bands of any format. 
+ *
  * Coordinates in @index are in pixels, with (0, 0) being the top-left corner 
  * of @in, and with y increasing down the image. Use vips_xyz() to build index
  * images. 
  *
  * @interpolate defaults to bilinear. 
  *
+ * By default, new pixels are filled with @background. This defaults to 
+ * zero (black). You can set other extend types with @extend. #VIPS_EXTEND_COPY 
+ * is better for image upsizing.
+ *
+ * Image are normally treated as unpremultiplied, so this operation can be used
+ * directly on PNG images. If your images have been through vips_premultiply(),
+ * set @premultiplied. 
+ *
  * This operation does not change xres or yres. The image resolution needs to
  * be updated by the application. 
  *

test/test-suite/test_resample.py

@@ -242,7 +242,7 @@ class TestResample:
         # distorted, but the rest should not be too bad
         a = r.crop(50, 0, im.width - 50, im.height).gaussblur(2)
         b = im.crop(50, 0, im.width - 50, im.height).gaussblur(2)
-        assert (a - b).abs().max() < 40
+        assert (a - b).abs().max() < 50
 
         # this was a bug at one point, strangely, if executed with debug
         # enabled