From 5901a7ca87963c57318fbbad28d49bcca0768ff4 Mon Sep 17 00:00:00 2001
From: John Cupitt <jcupitt@gmail.com>
Date: Tue, 8 Sep 2009 11:25:38 +0000
Subject: [PATCH] im_stats() and im_deviate() done

---
 libvips/arithmetic/im_deviate.c | 17 +++++++-----
 libvips/arithmetic/im_stats.c   | 47 +++------------------------------
 2 files changed, 14 insertions(+), 50 deletions(-)

diff --git a/libvips/arithmetic/im_deviate.c b/libvips/arithmetic/im_deviate.c
index dcb14da3..d914a87b 100644
--- a/libvips/arithmetic/im_deviate.c
+++ b/libvips/arithmetic/im_deviate.c
@@ -169,8 +169,8 @@ deviate_scan( void *in, int n, void *seq, void *a, void *b )
 		g_assert( 0 );
 	}
 
-	ss2[0] += s;
-	ss2[1] += s2;
+	ss2[0] = s;
+	ss2[1] = s2;
 
 	return( 0 );
 }
@@ -194,7 +194,8 @@ int
 im_deviate( IMAGE *in, double *out )
 {
 	double global_ss2[2];
-	gint64 N;
+	double s, s2;
+	gint64 vals;
 
 	/* Check our args. 
 	 */
@@ -208,7 +209,8 @@ im_deviate( IMAGE *in, double *out )
 	global_ss2[0] = 0.0;
 	global_ss2[1] = 0.0;
 	if( im__wrapscan( in, 
-		deviate_start, deviate_scan, deviate_stop, in, &global_ss2 ) ) 
+		deviate_start, deviate_scan, deviate_stop, in, global_ss2 ) ) 
+		return( -1 );
 
 	/*
 	  
@@ -220,9 +222,10 @@ im_deviate( IMAGE *in, double *out )
 
 	/* Calculate and return deviation. Add a fabs to stop sqrt(<=0).
 	 */
-	N = (gint64) in->Xsize * in->Ysize * in->Bands;
-	*out = sqrt( fabs( global_ss2[1] - 
-		(global_ss2[0] * global_ss2[0] / N) ) / (N - 1) );
+	vals = (gint64) in->Xsize * (gint64) in->Ysize * (gint64) in->Bands;
+	s = global_ss2[0];
+	s2 = global_ss2[1];
+	*out = sqrt( fabs( s2 - (s * s / vals) ) / (vals - 1) );
 
 	return( 0 );
 }
diff --git a/libvips/arithmetic/im_stats.c b/libvips/arithmetic/im_stats.c
index fcf5aca4..ab7172ff 100644
--- a/libvips/arithmetic/im_stats.c
+++ b/libvips/arithmetic/im_stats.c
@@ -164,46 +164,7 @@ stats_scan( void *in, int n, void *seq, void *a, void *b )
 	/* Now generate code for all types. 
 	 */
 	switch( im->BandFmt ) {
-	case IM_BANDFMT_UCHAR:	
-		/*
-		LOOP( unsigned char ); break; 
-		 */
-
-	for( z = 0; z < im->Bands; z++ ) { 
-		PEL *q = (PEL *) in + z; 
-		double *row = stats + z * 4; 
-		PEL small, big; 
-		double sum, sum2; 
-
-		printf( "z = %d\n", z );
-		
-		small = row[0]; 
-		big = row[1]; 
-		sum = row[2]; 
-		sum2 = row[3]; 
-		
-		for( x = 0; x < n; x++ ) { 
-			PEL value = *q; 
-
-			printf( "value = %d\n", value );
-			
-			sum += value;
-			sum2 += (double) value * (double) value;
-			if( value > big ) 
-				big = value; 
-			else if( value < small ) 
-				small = value;
-			
-			q += im->Bands; 
-		}
-		
-		row[0] = small; 
-		row[1] = big; 
-		row[2] = sum; 
-		row[3] = sum2; 
-	}
-		break;
-
+	case IM_BANDFMT_UCHAR:	LOOP( unsigned char ); break;
 	case IM_BANDFMT_CHAR:	LOOP( signed char ); break; 
 	case IM_BANDFMT_USHORT:	LOOP( unsigned short ); break; 
 	case IM_BANDFMT_SHORT:	LOOP( signed short ); break; 
@@ -230,7 +191,7 @@ im_stats( IMAGE *im )
 {	
 	DOUBLEMASK *out;
 	double *row;
-	gint64 pels, vals, z;
+	gint64 pels, vals;
 	double *global_stats;
 	int i, j;
 	double value;
@@ -255,7 +216,7 @@ im_stats( IMAGE *im )
 	 * separately.
 	 */
 	if( im__wrapscan( im, stats_start, stats_scan, stats_stop, 
-		im, &global_stats ) ) 
+		im, global_stats ) ) 
 		return( NULL );
 
 	/* Calculate mean, deviation, plus overall stats.
@@ -274,7 +235,7 @@ im_stats( IMAGE *im )
 	vals = pels * im->Bands;
 
 	for( i = 0; i < im->Bands; i++ ) {
-		row = out->coeff + (z + 1) * 6;
+		row = out->coeff + (i + 1) * 6;
 		for( j = 0; j < 4; j++ )
 			row[j] = global_stats[i * 4 + j];