revise unpremultiply, again

We were not detecting division by zero carefully enough, nor clipping
the alpha range sufficiently in unpremultiply.

see https://github.com/libvips/libvips/issues/1941

also see https://github.com/libvips/libvips/pull/1675 for another
difficult test case

ChangeLog

@@ -7,6 +7,7 @@
 - fix includes of glib headers in C++ [lovell]
 - fix build with more modern librsvg [lovell]
 - fix a possible segv with very wide images [f1ac]
+- revise premultiply, again [jjonesrs]
 
 18/12/20 started 8.10.5
 - fix potential /0 in animated webp load [lovell]

libvips/conversion/unpremultiply.c

@@ -7,6 +7,8 @@
  * 	- max_alpha defaults to 65535 for RGB16/GREY16
  * 24/11/17 lovell
  * 	- match normalised alpha to output type
+ * 27/2/21 jjonesrs
+ * 	- revise range clipping and 1/x, again
  */
 
 /*
@@ -78,18 +80,12 @@ G_DEFINE_TYPE( VipsUnpremultiply, vips_unpremultiply, VIPS_TYPE_CONVERSION );
 	\
 	for( x = 0; x < width; x++ ) { \
 		IN alpha = p[alpha_band]; \
-		\
+		IN clip_alpha = VIPS_CLIP( 0.0, alpha, max_alpha ); \
-		if( alpha != 0 ) { \
+		OUT factor = max_alpha * vips_recip( clip_alpha ); \
-			OUT factor = max_alpha / alpha; \
 		\
 		for( i = 0; i < alpha_band; i++ ) \
 			q[i] = factor * p[i]; \
-			q[alpha_band] = alpha; \
+		q[alpha_band] = clip_alpha; \
-		} \
-		else \
-			for( i = 0; i < alpha_band + 1; i++ ) \
-				q[i] = 0; \
-		\
 		for( i = alpha_band + 1; i < bands; i++ ) \
 			q[i] = p[i]; \
 		\
@@ -106,21 +102,13 @@ G_DEFINE_TYPE( VipsUnpremultiply, vips_unpremultiply, VIPS_TYPE_CONVERSION );
 	\
 	for( x = 0; x < width; x++ ) { \
 		IN alpha = p[3]; \
-		\
+		IN clip_alpha = VIPS_CLIP( 0.0, alpha, max_alpha ); \
-		if( alpha != 0 ) { \
+		OUT factor = max_alpha * vips_recip( clip_alpha ); \
-			OUT factor = max_alpha / alpha; \
 		\
 		q[0] = factor * p[0]; \
 		q[1] = factor * p[1]; \
 		q[2] = factor * p[2]; \
-			q[3] = alpha; \
+		q[3] = clip_alpha; \
-		} \
-		else { \
-			q[0] = 0; \
-			q[1] = 0; \
-			q[2] = 0; \
-			q[3] = 0; \
-		} \
 		\
 		p += 4; \
 		q += 4; \

libvips/include/vips/util.h

@@ -84,6 +84,20 @@ extern "C" {
 #define VIPS_FMIN( A, B ) VIPS_MIN( A, B )
 #endif
 
+/* Calculate 1.0 / x, avoiding division by zero when near zero.
+ */
+static inline double
+vips_recip( double x )
+{
+	static const double epsilon = 1e-10;
+
+	int sign = x < 0.0 ? -1.0 : 1.0;
+
+	return( sign * x >= epsilon ?
+		1.0 / x :
+		sign / epsilon );
+}
+
 /* Testing status before the function call saves a lot of time.
  */
 #define VIPS_ONCE( ONCE, FUNC, CLIENT ) \