From 84e04e18d90bc9f5ee56b736cfa6a6fd421c8d85 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 6 Jan 2010 17:50:29 +0100
Subject: Hit finding "improvements"

---
 src/process_images.c | 48 +++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 41 insertions(+), 7 deletions(-)

diff --git a/src/process_images.c b/src/process_images.c
index 8effed1e..cd06a97e 100644
--- a/src/process_images.c
+++ b/src/process_images.c
@@ -43,32 +43,66 @@ static void show_help(const char *s)
 static int sum_of_peaks(struct image *image)
 {
 	int x, y;
-	int integr = 0;
+	int integr, n;
+	float fintegr, mean, sd, th;
 
-	for ( x=0; x<400; x++ ) {
+	/* Measure mean */
+	integr = 0;
+	n = 0;
+	for ( x=0; x<1024; x++ ) {
 	for ( y=600; y<1024; y++ ) {
 
 		int val;
-
+		if ( (x>400) && (x<600) ) continue;
 		val = image->data[x+image->height*y];
+		if ( val < 0 ) continue;
+		integr += val;
+		n++;
+
+	}
+	}
+	mean = (float)integr / n;  /* As integer to keep maths fast */
+
+	/* Standard deviation */
+	integr = 0;
+	for ( x=0; x<1024; x++ ) {
+	for ( y=600; y<1024; y++ ) {
 
-		if ( val > 1000 ) integr+=val;
+		float val;
+
+		if ( (x>400) && (x<600) ) continue;
+		val = (float)image->data[x+image->height*y];
+		if ( val < 0 ) continue;
+
+		val -= mean;
+		val = powf(val, 2.0);
+		fintegr += val;
 
 	}
 	}
-	for ( x=700; x<1023; x++ ) {
+	sd = sqrtf(fintegr / n);
+
+	/* Threshold */
+	th = mean + 5*sd;
+	STATUS("mean=%f ,sd=%f, th=%f\n", mean, sd, th);
+
+	/* Sum intensity above threshold */
+	integr = 0;
+	for ( x=0; x<1024; x++ ) {
 	for ( y=600; y<1024; y++ ) {
 
 		int val;
 
+		/* Chop out streaky region */
+		if ( (x>400) && (x<600) ) continue;
+
 		val = image->data[x+image->height*y];
 
-		if ( val > 1000 ) integr+=val;
+		if ( val > th ) integr+=val;
 
 	}
 	}
 
-
 	return integr;
 }
 
-- 
cgit v1.2.3