From d3ccadf212977c732e8237565e7bc046948a8bb6 Mon Sep 17 00:00:00 2001
From: John Cupitt <jcupitt@gmail.com>
Date: Sat, 27 Feb 2021 15:16:25 +0000
Subject: [PATCH] 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                          |  1 +
 libvips/conversion/unpremultiply.c | 38 ++++++++++--------------------
 libvips/include/vips/util.h        | 14 +++++++++++
 3 files changed, 28 insertions(+), 25 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 0b6981a0..91c231ba 100644
--- a/ChangeLog
+++ b/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]
diff --git a/libvips/conversion/unpremultiply.c b/libvips/conversion/unpremultiply.c
index 219fac28..5efd7696 100644
--- a/libvips/conversion/unpremultiply.c
+++ b/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 ); \
+		OUT factor = max_alpha * vips_recip( clip_alpha ); \
 		\
-		if( alpha != 0 ) { \
-			OUT factor = max_alpha / alpha; \
-			\
-			for( i = 0; i < alpha_band; i++ ) \
-				q[i] = factor * p[i]; \
-			q[alpha_band] = alpha; \
-		} \
-		else \
-			for( i = 0; i < alpha_band + 1; i++ ) \
-				q[i] = 0; \
-		\
+		for( i = 0; i < alpha_band; i++ ) \
+			q[i] = factor * p[i]; \
+		q[alpha_band] = clip_alpha; \
 		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 ); \
+		OUT factor = max_alpha * vips_recip( clip_alpha ); \
 		\
-		if( alpha != 0 ) { \
-			OUT factor = max_alpha / alpha; \
-			\
-			q[0] = factor * p[0]; \
-			q[1] = factor * p[1]; \
-			q[2] = factor * p[2]; \
-			q[3] = alpha; \
-		} \
-		else { \
-			q[0] = 0; \
-			q[1] = 0; \
-			q[2] = 0; \
-			q[3] = 0; \
-		} \
+		q[0] = factor * p[0]; \
+		q[1] = factor * p[1]; \
+		q[2] = factor * p[2]; \
+		q[3] = clip_alpha; \
 		\
 		p += 4; \
 		q += 4; \
diff --git a/libvips/include/vips/util.h b/libvips/include/vips/util.h
index feb49f72..2b3e8726 100644
--- a/libvips/include/vips/util.h
+++ b/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 ) \