Change the sanity check so that it does what we say it does

4 files changed, 30 insertions, 55 deletions

diff --git a/src/index.c b/src/index.c
index 862f5145..ced8bbb8 100644
--- a/src/index.c
+++ b/src/index.c
@@ -30,6 +30,7 @@
 #include "index.h"
 #include "index-priv.h"
 #include "reax.h"
+#include "geometry.h"
 
 
 /* Base class constructor for unspecialised indexing private data */
@@ -204,8 +205,12 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
 			if ( new_cell == NULL ) continue;
 
 			/* Sanity check */
+			image->reflections = find_intersections(image,
+			                                        new_cell);
 			if ( !config_insane &&
-			     !peak_sanity_check(image, new_cell, 0, 0.1) ) {
+			     !peak_sanity_check(image->reflections,
+			                        image->features) )
+			{
 				cell_free(new_cell);
 				continue;
 			}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 384195ac..d27ff4d4 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -314,20 +314,16 @@ static void process_image(void *pp, int cookie)
 	image.data = data_for_measurement;
 
 	/* Calculate orientation matrix (by magic) */
+	image.div = beam->divergence;
+	image.bw = beam->bandwidth;
+	image.profile_radius = 0.0001e9;
 	index_pattern(&image, cell, indm, pargs->static_args.cellr,
 		      config_verbose, pargs->static_args.ipriv,
 		      pargs->static_args.config_insane);
 
-	if ( image.indexed_cell != NULL ) pargs->indexable = 1;
-
-	/* Measure intensities */
 	if ( image.indexed_cell != NULL ) {
 
-		image.div = beam->divergence;
-		image.bw = beam->bandwidth;
-		image.profile_radius = 0.0001e9;
-		image.reflections = find_intersections(&image,
-			                               image.indexed_cell);
+		pargs->indexable = 1;
 
 		if ( image.reflections != NULL ) {
 			integrate_reflections(&image, config_polar,
diff --git a/src/peaks.c b/src/peaks.c
index 67b8d9c4..1c3564d1 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -31,6 +31,7 @@
 #include "detector.h"
 #include "filters.h"
 #include "diffraction.h"
+#include "reflist-utils.h"
 
 
 /* How close a peak must be to an indexed position to be considered "close"
@@ -449,62 +450,36 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
 }
 
 
-int peak_sanity_check(struct image *image, UnitCell *cell,
-                      int circular_domain, double domain_r)
+int peak_sanity_check(RefList *rlist, ImageFeatureList *flist)
 {
 	int i;
 	int n_feat = 0;
 	int n_sane = 0;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	double aslen, bslen, cslen;
 
-	/* "Borrow" direction values to get reciprocal lengths */
-	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	aslen = modulus(ax, ay, az);
-	bslen = modulus(bx, by, bz);
-	cslen = modulus(cx, cy, cz);
-
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
-	fesetround(1);  /* Round towards nearest */
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
+	for ( i=0; i<image_feature_count(flist); i++ ) {
 
 		double dist;
-		struct rvec q;
 		struct imagefeature *f;
-		double hd, kd, ld;
-		signed int h, k, l;
-		double dh, dk, dl;
+		Reflection *refl;
+		RefListIterator *iter;
 
-		f = image_get_feature(image->features, i);
+		f = image_get_feature(flist, i);
 		if ( f == NULL ) continue;
 		n_feat++;
 
-		/* Get closest hkl */
-		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
-
-		hd = q.u * ax + q.v * ay + q.w * az;
-		kd = q.u * bx + q.v * by + q.w * bz;
-		ld = q.u * cx + q.v * cy + q.w * cz;
-
-		h = lrint(hd);  k = lrint(kd);  l = lrint(ld);
-
-		dh = hd - h;  dk = kd - k;  dl = ld - l;
-
-		if ( circular_domain ) {
-
-			/* Circular integration domain */
-			dist = sqrt(pow(dh*aslen, 2.0) + pow(dk*bslen, 2.0)
-			                              + pow(dl*cslen, 2.0));
-			if ( dist <= domain_r ) n_sane++;
-
-		} else {
-
-			/* "Crystallographic" integration domain */
-			dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
-			if ( dist <= domain_r ) n_sane++;
+		/* Find closest predicted peak */
+
+		for ( refl = first_refl(rlist, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			double fs, ss;
+			get_detector_pos(refl, &fs, &ss);
+			dist = sqrt(pow(fs-f->fs, 2.0) + pow(ss-f->ss, 2.0));
+			if ( dist < PEAK_CLOSE ) {
+				n_sane++;
+				continue;
+			}
 		}
 
 	}
diff --git a/src/peaks.h b/src/peaks.h
index 3d21bbde..07c257e6 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -27,8 +27,7 @@ extern void search_peaks(struct image *image, float threshold,
 extern void integrate_reflections(struct image *image,
                                   int polar, int use_closer, int bgsub);
 
-extern int peak_sanity_check(struct image *image, UnitCell *cell,
-                             int circular_domain, double domain_r);
+extern int peak_sanity_check(RefList *rlist, ImageFeatureList *flist);
 
 /* Exported so it can be poked by integration_check */
 extern int integrate_peak(struct image *image, int cfs, int css,