Add peak sanity check

4 files changed, 86 insertions, 20 deletions

diff --git a/src/image.h b/src/image.h
index 0479debc..85265028 100644
--- a/src/image.h
+++ b/src/image.h
@@ -67,6 +67,16 @@ struct rvec
 };
 
 
+struct reflhit {
+	signed int h;
+	signed int k;
+	signed int l;
+	double min_distance;
+	int x;
+	int y;
+};
+
+
 /* Structure describing an image */
 struct image {
 
@@ -77,6 +87,8 @@ struct image {
 	int			ncells;
 	struct detector		det;
 	const char		*filename;
+	struct reflhit		*hits;
+	int			n_hits;
 
 	int			id;
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index c21dc22f..db433369 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -67,6 +67,7 @@ struct process_args
 struct process_result
 {
 	int hit;
+	int peaks_sane;
 };
 
 
@@ -194,6 +195,10 @@ static struct image *get_simage(struct image *template, int alternate)
 	image->indexed_cell = template->indexed_cell;
 	image->f0 = template->f0;
 
+	/* Prevent muppetry */
+	image->hits = NULL;
+	image->n_hits = 0;
+
 	return image;
 }
 
@@ -252,6 +257,8 @@ static void *process_image(void *pargsv)
 	image.indexed_cell = NULL;
 	image.id = pargs->id;
 	image.filename = filename;
+	image.hits = NULL;
+	image.n_hits = 0;
 
 	STATUS("Processing '%s'\n", image.filename);
 
@@ -318,6 +325,15 @@ static void *process_image(void *pargsv)
 	/* No cell at this point?  Then we're done. */
 	if ( image.indexed_cell == NULL ) goto done;
 
+	/* Sanity check */
+	if ( !peak_sanity_check(&image, image.indexed_cell) ) {
+		STATUS("Failed peak sanity check.\n");
+		result->peaks_sane = 0;
+		goto done;
+	} else {
+		result->peaks_sane = 1;
+	}
+
 	/* Measure intensities if requested */
 	if ( config_nearbragg ) {
 		/* Use original data (temporarily) */
diff --git a/src/peaks.c b/src/peaks.c
index 20779bb3..fe5dcbea 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -44,15 +44,6 @@
 /* Degree of polarisation of X-ray beam */
 #define POL (1.0)
 
-struct reflhit {
-	signed int h;
-	signed int k;
-	signed int l;
-	double min_distance;
-	int x;
-	int y;
-};
-
 
 #define PEAK_WINDOW_SIZE (10)
 #define MAX_PEAKS (2048)
@@ -385,23 +376,17 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
 }
 
 
-void output_intensities(struct image *image, UnitCell *cell,
-                        pthread_mutex_t *mutex, int unpolar)
+static int find_projected_peaks(struct image *image, UnitCell *cell)
 {
 	int x, y;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
-	double a, b, c, al, be, ga;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	struct reflhit hits[MAX_HITS];
+	struct reflhit *hits;
 	int n_hits = 0;
-	int i;
-	int n_found;
-	int n_indclose = 0;
-	int n_foundclose = 0;
+
+	hits = malloc(sizeof(struct reflhit)*MAX_HITS);
+	if ( hits == NULL ) return 0;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
@@ -462,6 +447,58 @@ void output_intensities(struct image *image, UnitCell *cell,
 	}
 
 	STATUS("Found %i reflections\n", n_hits);
+	image->hits = hits;
+	image->n_hits = n_hits;
+
+	return n_hits;
+}
+
+
+int peak_sanity_check(struct image *image, UnitCell *cell)
+{
+	int i;
+	const int n_hits = image->n_hits;
+	const struct reflhit *hits = image->hits;
+	int n_sane = 0;
+
+	find_projected_peaks(image, cell);
+	if ( image->n_hits == 0 ) return 0;  /* Failed sanity check: no peaks */
+
+	for ( i=0; i<n_hits; i++ ) {
+
+		double d;
+		int idx;
+		struct imagefeature *f;
+
+		f = image_feature_closest(image->features, hits[i].x, hits[i].y,
+		                          &d, &idx);
+		if ( (f != NULL) && (d < PEAK_CLOSE) ) {
+			n_sane++;
+		}
+
+	}
+
+	if ( (float)n_sane / (float)n_hits < 0.8 ) return 0;
+
+	return 1;
+}
+
+
+void output_intensities(struct image *image, UnitCell *cell,
+                        pthread_mutex_t *mutex, int unpolar)
+{
+	int i;
+	int n_found;
+	int n_indclose = 0;
+	int n_foundclose = 0;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double a, b, c, al, be, ga;
+	int n_hits = image->n_hits;
+	struct reflhit *hits = image->hits;
+
+	if ( image->n_hits == 0 ) return;
 
 	/* Get exclusive access to the output stream if necessary */
 	if ( mutex != NULL ) pthread_mutex_lock(mutex);
diff --git a/src/peaks.h b/src/peaks.h
index fa912455..10cf5b0a 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -23,5 +23,6 @@ extern void search_peaks(struct image *image);
 extern void dump_peaks(struct image *image, pthread_mutex_t *mutex);
 extern void output_intensities(struct image *image, UnitCell *cell,
                                pthread_mutex_t *mutex, int unpolar);
+extern int peak_sanity_check(struct image *image, UnitCell *cell);
 
 #endif	/* PEAKS_H */