Perform prediction refinement straight after indexing

This allows indexing to be attempted again (either a new method or with
"retry") if the prediction refinement fails, increasing overall indexing
rate.

Side-effect: there are some hoops which would need to be jumped through
to store the profile radius before refinement and hence enable
scripts/plot-predict-refine to work.  For now, we'll ignore this as it's
clear that the prediction refinement is working.

6 files changed, 69 insertions, 138 deletions

diff --git a/Makefile.am b/Makefile.am
index 1940affc..ae4b2913 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -179,7 +179,7 @@ script_DATA = scripts/alternate-stream scripts/cell-please \
               scripts/gen-sfs-expand scripts/add-beam-params \
               scripts/find-pairs scripts/plot-cc-and-scale.R \
               scripts/ave-resolution scripts/crystal-frame-number \
-              scripts/plot-radius-resolution scripts/plot-predict-refine \
+              scripts/plot-radius-resolution \
               scripts/detector-shift scripts/turbo-index
 
 EXTRA_DIST += $(script_DATA)
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index cd7047dd..82ad8f7e 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -280,6 +280,26 @@ void image_add_crystal(struct image *image, Crystal *cryst)
 }
 
 
+void remove_flagged_crystals(struct image *image)
+{
+	int i;
+
+	for ( i=0; i<image->n_crystals; i++ ) {
+		if ( crystal_get_user_flag(image->crystals[i]) ) {
+			int j;
+			Crystal *deleteme = image->crystals[i];
+			cell_free(crystal_get_cell(deleteme));
+			crystal_free(deleteme);
+			for ( j=i; j<image->n_crystals-1; j++ ) {
+				image->crystals[j] = image->crystals[j+1];
+			}
+			image->n_crystals--;
+		}
+	}
+
+}
+
+
 /* Free all crystals, including their RefLists and UnitCells */
 void free_all_crystals(struct image *image)
 {
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index e9cd76fc..ddbfc9dc 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -227,6 +227,7 @@ 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 remove_flagged_crystals(struct image *image);
 extern void free_all_crystals(struct image *image);
 
 #ifdef __cplusplus
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 2b84a621..6a74a958 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -55,6 +55,7 @@
 #include "geometry.h"
 #include "cell-utils.h"
 #include "felix.h"
+#include "predict-refine.h"
 
 
 static int debug_index(struct image *image)
@@ -234,36 +235,61 @@ void map_all_peaks(struct image *image)
 static int try_indexer(struct image *image, IndexingMethod indm,
                        IndexingPrivate *ipriv)
 {
+	int i, r, n_bad;
+
 	switch ( indm & INDEXING_METHOD_MASK ) {
 
 		case INDEXING_NONE :
 		return 0;
 
 		case INDEXING_DIRAX :
-		return run_dirax(image, ipriv);
+		r = run_dirax(image, ipriv);
+		break;
 
 		case INDEXING_ASDF :
-		return run_asdf(image, ipriv);
+		r = run_asdf(image, ipriv);
+		break;
 
 		case INDEXING_MOSFLM :
-		return run_mosflm(image, ipriv);
+		r = run_mosflm(image, ipriv);
+		break;
 
 		case INDEXING_XDS :
-		return run_xds(image, ipriv);
+		r = run_xds(image, ipriv);
+		break;
 
 		case INDEXING_DEBUG :
-		return debug_index(image);
+		r = debug_index(image);
+		break;
 
 		case INDEXING_FELIX :
-		return felix_index(image, ipriv);
+		r = felix_index(image, ipriv);
+		break;
 
 		default :
 		ERROR("Unrecognised indexing method: %i\n", indm);
-		break;
+		return 0;
 
 	}
 
-	return 0;
+	/* Attempt prediction refinement */
+	n_bad = 0;
+	for ( i=0; i<r; i++ ) {
+		Crystal *cr = image->crystals[image->n_crystals-i-1];
+		crystal_set_image(cr, image);
+		crystal_set_user_flag(cr, 0);
+		crystal_set_profile_radius(cr, 0.02e9);
+		crystal_set_mosaicity(cr, 0.0);
+		if ( refine_prediction(image, cr) != 0 ) {
+			crystal_set_user_flag(cr, 1);
+			n_bad++;
+		}
+	}
+
+	remove_flagged_crystals(image);
+
+	if ( n_bad == r ) return 0;
+	return r;
 }
 
 
diff --git a/scripts/plot-predict-refine b/scripts/plot-predict-refine
deleted file mode 100755
index 04f6ffa8..00000000
--- a/scripts/plot-predict-refine
+++ /dev/null
@@ -1,64 +0,0 @@
-#!/usr/bin/env python
-
-#
-# Visualise behaviour of prediction refinement
-#
-# Copyright (c) 2015 Deutsches Elektronen-Synchrotron DESY,
-#                    a research centre of the Helmholtz Association.
-#
-# Author:
-#    2015 Thomas White <taw@physics.org>
-#
-
-import sys
-import os
-import re
-import matplotlib.pyplot as plt
-
-
-def show_frac(ratios, v):
-    n = len([x for x in ratios if x < v])
-    print "%4i (%.2f%%) new/old ratios are below %.2f" % (n, 100*float(n)/len(ratios), v)
-
-
-f = open(sys.argv[1], 'r')
-
-oldR = []
-newR = []
-
-prog1 = re.compile("^predict_refine/R\sold\s=\s([0-9\.\-]+)\snew\s=\s([0-9\.\-]+)\s")
-
-while True:
-
-    fline = f.readline()
-    if not fline:
-        break
-
-    match = prog1.match(fline)
-    if match:
-        old = float(match.group(1))
-        new = float(match.group(2))
-        oldR.append(old)
-        newR.append(new)
-
-f.close()
-
-mean_oldR = sum(oldR) / len(oldR)
-mean_newR = sum(newR) / len(newR)
-ratios = [new/old for new,old in zip(newR, oldR)];
-print 'Mean profile radius before: %.2e, after: %.2e nm^-1' % (mean_oldR,mean_newR)
-show_frac(ratios, 1.2)
-show_frac(ratios, 1.1)
-show_frac(ratios, 1.0)
-show_frac(ratios, 0.9)
-show_frac(ratios, 0.8)
-
-#plt.plot(oldR, newR, 'rx')
-plt.hist(ratios, 50)
-#plt.axis([-2,2,-2,2])
-plt.title('Profile radius before and after refinement')
-plt.xlabel('x shift / mm')
-plt.ylabel('y shift / mm')
-plt.grid(True)
-plt.show()
-
diff --git a/src/process_image.c b/src/process_image.c
index 5d4212e2..9b0d79d8 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -54,29 +54,6 @@
 #include "im-sandbox.h"
 
 
-static void try_refine_autoR(struct image *image, Crystal *cr)
-{
-	double old_R, new_R;
-	char notes[1024];
-
-	refine_radius(cr, image);
-	old_R = crystal_get_profile_radius(cr);
-
-	if ( refine_prediction(image, cr) ) {
-		crystal_set_user_flag(cr, 1);
-		return;
-	}
-
-	/* Estimate radius again with better geometry */
-	refine_radius(cr, image);
-	new_R = crystal_get_profile_radius(cr);
-
-	snprintf(notes, 1024, "predict_refine/R old = %.5f new = %.5f nm^-1",
-	                      old_R/1e9, new_R/1e9);
-	crystal_add_notes(cr, notes);
-}
-
-
 static float **backup_image_data(float **dp, struct detector *det)
 {
 	float **bu;
@@ -130,7 +107,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	int r;
 	int ret;
 	char *rn;
-	int n_crystals_left;
 	float **prefilter;
 	int any_crystals;
 
@@ -218,6 +194,18 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		return;
 	}
 
+	/* Set beam parameters */
+	if ( iargs->fix_divergence >= 0.0 ) {
+		image.div = iargs->fix_divergence;
+	} else {
+		image.div = 0.0;
+	}
+	if ( iargs->fix_bandwidth >= 0.0 ) {
+		image.bw = iargs->fix_bandwidth;
+	} else {
+		image.bw = 0.00000001;
+	}
+
 	/* Index the pattern */
 	index_pattern(&image, iargs->indm, iargs->ipriv);
 
@@ -229,22 +217,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	}
 	free(rn);
 
-	for ( i=0; i<image.n_crystals; i++ ) {
-		crystal_set_image(image.crystals[i], &image);
-		crystal_set_user_flag(image.crystals[i], 0);
-	}
-
 	/* Set beam/crystal parameters */
-	if ( iargs->fix_divergence >= 0.0 ) {
-		image.div = iargs->fix_divergence;
-	} else {
-		image.div = 0.0;
-	}
-	if ( iargs->fix_bandwidth >= 0.0 ) {
-		image.bw = iargs->fix_bandwidth;
-	} else {
-		image.bw = 0.00000001;
-	}
 	if ( iargs->fix_profile_r >= 0.0 ) {
 		for ( i=0; i<image.n_crystals; i++ ) {
 			crystal_set_profile_radius(image.crystals[i],
@@ -259,34 +232,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	}
 
 	if ( iargs->fix_profile_r < 0.0 ) {
-
 		for ( i=0; i<image.n_crystals; i++ ) {
-			if ( iargs->predict_refine ) {
-				try_refine_autoR(&image, image.crystals[i]);
-			} else {
-				refine_radius(image.crystals[i], &image);
-			}
+			refine_radius(image.crystals[i], &image);
 		}
-
-	} else {
-
-		for ( i=0; i<image.n_crystals; i++ ) {
-			if ( iargs->predict_refine ) {
-				refine_prediction(&image, image.crystals[i]);
-			}
-		}
-
-	}
-
-	/* If there are no crystals left, set the indexing flag back to zero */
-	n_crystals_left = 0;
-	for ( i=0; i<image.n_crystals; i++ ) {
-		if ( crystal_get_user_flag(image.crystals[i]) == 0 ) {
-			n_crystals_left++;
-		}
-	}
-	if ( n_crystals_left == 0 ) {
-		image.indexed_by = INDEXING_NONE;
 	}
 
 	/* Integrate! */