"Retry" mechanism for indexing

This increases the indexing rate a bit in situations where there are
lots of weak peaks which, although they may be real, don't help
indexing.  The weakest 10% of peaks get cut out and the indexing re-run.

This also allows multiple hits to be indexed, using the "inelegant peak
subtraction method", by retrying indexing in the same way after deleting
peaks which are accounted for by the lattice just found.

3 files changed, 153 insertions, 1 deletions

diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 3dfccea3..cd7047dd 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -95,6 +95,45 @@ ImageFeatureList *image_feature_list_new()
 }
 
 
+static int comp(const void *a, const void *b)
+{
+	const struct imagefeature *ap = a;
+	const struct imagefeature *bp = b;
+
+	return ap->intensity < bp->intensity;
+}
+
+
+/* Strongest first.  Returned list is guaranteed not to have any holes
+ * (feature->valid = 0) */
+ImageFeatureList *sort_peaks(ImageFeatureList *flist)
+{
+	ImageFeatureList *n = image_feature_list_new();
+	int nf, i;
+
+	if ( n == NULL ) return NULL;
+
+	n->features = malloc(flist->n_features*sizeof(struct imagefeature));
+	if ( n->features == NULL ) {
+		free(n);
+		return NULL;
+	}
+
+	nf = 0;
+	for ( i=0; i<flist->n_features; i++ ) {
+		struct imagefeature *f;
+		f = image_get_feature(flist, i);
+		if ( f == NULL ) continue;
+		n->features[nf++] = flist->features[i];
+	}
+	n->n_features = nf;
+
+	qsort(n->features, nf, sizeof(struct imagefeature), comp);
+
+	return n;
+}
+
+
 void image_feature_list_free(ImageFeatureList *flist)
 {
 	if ( !flist ) return;
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 37904308..e9cd76fc 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -224,6 +224,7 @@ extern Reflection *image_reflection_closest(RefList *rlist,
 
 extern int image_feature_count(ImageFeatureList *flist);
 extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx);
+extern ImageFeatureList *sort_peaks(ImageFeatureList *flist);
 
 extern void image_add_crystal(struct image *image, Crystal *cryst);
 extern void free_all_crystals(struct image *image);
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 99e40053..2e54859d 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -13,6 +13,7 @@
  *   2012      Lorenzo Galli
  *   2013      Cornelius Gati <cornelius.gati@cfel.de>
  *   2015      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *   2014      Takanori Nakane <nakane.t@gmail.com>
  *
  * This file is part of CrystFEL.
  *
@@ -40,6 +41,7 @@
 #include <math.h>
 #include <string.h>
 #include <assert.h>
+#include <fenv.h>
 
 #include "image.h"
 #include "utils.h"
@@ -265,10 +267,91 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 }
 
 
+static int delete_weakest_peaks(ImageFeatureList *peaks)
+{
+	int i;
+	int np, ndel, n;
+
+	np = image_feature_count(peaks);
+	n = 0;
+	for ( i=0; i<np; i++ ) {
+		if ( image_get_feature(peaks, i) != NULL ) n++;
+	}
+
+	if ( n < 6 ) return 1;
+
+	/* Delete 10% weakest peaks */
+	ndel = n/10;
+	if ( ndel == 0 ) return 1;
+
+	for ( i=n-1; i>n-ndel-1; i-- ) {
+		image_remove_feature(peaks, i);
+	}
+
+	return 0;
+}
+
+
+static int delete_explained_peaks(struct image *image, Crystal *cr)
+{
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	const double min_dist = 0.25;
+	int i, nspots = 0, nindexed = 0;
+
+	/* Round towards nearest */
+	fesetround(1);
+
+	/* Cell basis vectors for this image */
+	cell_get_cartesian(crystal_get_cell(cr), &ax, &ay, &az,
+	                   &bx, &by, &bz, &cx, &cy, &cz);
+
+	/* Loop over peaks, checking proximity to nearest reflection */
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+		struct imagefeature *f;
+		struct rvec q;
+		double h, k, l, hd, kd, ld;
+		double dsq;
+
+		f = image_get_feature(image->features, i);
+		if ( f == NULL ) continue;
+		nspots++;
+
+		/* Reciprocal space position of found peak */
+		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
+
+		/* Decimal and fractional Miller indices of nearest
+		 * reciprocal lattice point */
+		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);
+
+		/* Check distance */
+		dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
+
+		if ( sqrt(dsq) < min_dist ) {
+			image_remove_feature(image->features, i);
+			nindexed++;
+		}
+	}
+
+	/* Return TRUE if not enough peaks to continue */
+	return (nspots - nindexed) < 5;
+
+}
+
+
 void index_pattern(struct image *image,
                    IndexingMethod *indms, IndexingPrivate **iprivs)
 {
 	int n = 0;
+	ImageFeatureList *peaks;
+	ImageFeatureList *orig;
 
 	if ( indms == NULL ) return;
 
@@ -276,14 +359,43 @@ void index_pattern(struct image *image,
 	image->crystals = NULL;
 	image->n_crystals = 0;
 
+	peaks = sort_peaks(image->features);
+	orig = image->features;
+	image->features = peaks;
+
 	while ( indms[n] != INDEXING_NONE ) {
 
-		if ( try_indexer(image, indms[n], iprivs[n]) ) break;
+		int done = 0;
+		int r;
+		int ntry = 0;
+
+		do {
+
+			r = try_indexer(image, indms[n], iprivs[n]);
+			ntry++;
+
+			if ( r == 0 ) {
+				/* Retry with fewer peaks */
+				done = delete_weakest_peaks(image->features);
+			} else {
+				/* Remove "used" spots and try for
+				 * another lattice */
+				Crystal *cr;
+				cr = image->crystals[image->n_crystals-1];
+				done = delete_explained_peaks(image, cr);
+			}
+
+			if ( ntry > 5 ) done = 1;
+
+		} while ( !done );
+
 		n++;
 
 	}
 
 	image->indexed_by = indms[n];
+	image->features = orig;
+	image_feature_list_free(peaks);
 }