From 2fb0f434a2823b0b2a9512ba2c71a784ebc68468 Mon Sep 17 00:00:00 2001
From: John Cupitt <jcupitt@gmail.com>
Date: Wed, 8 Apr 2015 09:40:31 +0100
Subject: [PATCH] avoid /0 in some cases in spcor

see https://github.com/jcupitt/libvips/issues/252 and others
---
 libvips/convolution/spcor.c | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/libvips/convolution/spcor.c b/libvips/convolution/spcor.c
index 20a65a79..c95b304a 100644
--- a/libvips/convolution/spcor.c
+++ b/libvips/convolution/spcor.c
@@ -36,6 +36,8 @@
  * 	- cleanups
  * 7/11/13
  * 	- redone as a class
+ * 8/4/15
+ * 	- avoid /0 for constant reference or zero image
  */
 
 /*
@@ -176,20 +178,20 @@ vips_spcor_pre_generate( VipsCorrelation *correlation )
 	sum3 = 0.0; \
 	for( j = 0; j < ref->Ysize; j++ ) { \
 		for( i = 0; i < sz; i += bands ) { \
-			/* Reference pel, and input pel. \
+			/* Reference pel and input pel. \
 			 */ \
 			IN ip = p1a[i]; \
 			IN rp = r1a[i]; \
 			\
-			/* Accumulate sum-of-squares-of- \
-			 * differences for input image. \
+			/* Accumulate sum-of-squares-of-differences for \
+			 * input image. \
 			 */ \
 			double t = ip - imean; \
 			sum2 += t * t; \
 			\
 			/* Accumulate product-of-difference from mean. \
 			 */ \
-			sum3 += (rp - spcor->rmean[b]) * (ip - imean); \
+			sum3 += (rp - spcor->rmean[b]) * t; \
 		} \
 		\
 		p1a += in_lsk; \
@@ -265,8 +267,15 @@ vips_spcor_correlation( VipsCorrelation *correlation,
 					return; 
 				}
 
-				c2 = sqrt( sum2 );
-				cc = sum3 / (spcor->c1[b] * c2);
+				c2 = spcor->c1[b] * sqrt( sum2 );
+
+				if( c2 == 0.0 )
+					/* Something like constant ref.
+					 * We regard this as uncorrelated.
+					 */
+					cc = 0.0;
+				else
+					cc = sum3 / c2;
 
 				*q++ = cc;
 			}