Indexing stuff

11 files changed, 279 insertions, 91 deletions

diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 160df719..0b2debda 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -53,6 +53,7 @@
 #include "dirax.h"
 #include "utils.h"
 #include "peaks.h"
+#include "cell-utils.h"
 
 
 #define DIRAX_VERBOSE 0
@@ -68,6 +69,13 @@ typedef enum {
 } DirAxInputType;
 
 
+struct dirax_private {
+	IndexingMethod          indm;
+	float                   *ltl;
+	UnitCell                *template;
+};
+
+
 struct dirax_data {
 
 	/* DirAx auto-indexing low-level stuff */
@@ -85,10 +93,56 @@ struct dirax_data {
 	int                     best_acl_nh;
 	int                     acls_tried[MAX_CELL_CANDIDATES];
 	int                     n_acls_tried;
+	UnitCell                *cur_cell;
+	int                     done;
+	int                     success;
+
+	struct dirax_private    *dp;
 
 };
 
 
+static int check_cell(struct dirax_private *dp, struct image *image,
+                      UnitCell *cell)
+{
+	UnitCell *out;
+	Crystal *cr;
+
+	if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) {
+
+		out = match_cell(cell, dp->template, 0, dp->ltl, 1);
+		if ( out == NULL ) return 0;
+
+	} else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
+
+		out = match_cell(cell, dp->template, 0, dp->ltl, 0);
+		if ( out == NULL ) return 0;
+
+	} else {
+		out = cell_new_from_cell(cell);
+	}
+
+	cr = crystal_new();
+	if ( cr == NULL ) {
+		ERROR("Failed to allocate crystal.\n");
+		return 0;
+	}
+
+	crystal_set_cell(cr, cell);
+
+	if ( dp->indm & INDEXING_CHECK_PEAKS ) {
+		if ( !peak_sanity_check(image, &cr, 1) ) {
+			crystal_free(cr);  /* Frees the cell as well */
+			return 0;
+		}
+	}
+
+	image_add_crystal(image, cr);
+
+	return 1;
+}
+
+
 static void dirax_parseline(const char *line, struct image *image,
                             struct dirax_data *dirax)
 {
@@ -119,7 +173,7 @@ static void dirax_parseline(const char *line, struct image *image,
 		if ( line[i] == 'R' ) rf = 1;
 		if ( (line[i] == 'D') && rf ) {
 			dirax->read_cell = 1;
-			image->candidate_cells[image->ncells] = cell_new();
+			dirax->cur_cell = cell_new();
 			return;
 		}
 		i++;
@@ -127,6 +181,7 @@ static void dirax_parseline(const char *line, struct image *image,
 
 	/* Parse unit cell vectors as appropriate */
 	if ( dirax->read_cell == 1 ) {
+
 		/* First row of unit cell values */
 		float ax, ay, az;
 		int r;
@@ -136,14 +191,16 @@ static void dirax_parseline(const char *line, struct image *image,
 			ERROR("Couldn't understand cell line:\n");
 			ERROR("'%s'\n", line);
 			dirax->read_cell = 0;
-			cell_free(image->candidate_cells[image->ncells]);
+			cell_free(dirax->cur_cell);
 			return;
 		}
-		cell_set_cartesian_a(image->candidate_cells[image->ncells],
+		cell_set_cartesian_a(dirax->cur_cell,
 		                     ax*1e-10, ay*1e-10, az*1e-10);
 		dirax->read_cell++;
 		return;
+
 	} else if ( dirax->read_cell == 2 ) {
+
 		/* Second row of unit cell values */
 		float bx, by, bz;
 		int r;
@@ -153,14 +210,16 @@ static void dirax_parseline(const char *line, struct image *image,
 			ERROR("Couldn't understand cell line:\n");
 			ERROR("'%s'\n", line);
 			dirax->read_cell = 0;
-			cell_free(image->candidate_cells[image->ncells]);
+			cell_free(dirax->cur_cell);
 			return;
 		}
-		cell_set_cartesian_b(image->candidate_cells[image->ncells],
+		cell_set_cartesian_b(dirax->cur_cell,
 		                     bx*1e-10, by*1e-10, bz*1e-10);
 		dirax->read_cell++;
 		return;
+
 	} else if ( dirax->read_cell == 3 ) {
+
 		/* Third row of unit cell values */
 		float cx, cy, cz;
 		int r;
@@ -170,13 +229,21 @@ static void dirax_parseline(const char *line, struct image *image,
 			ERROR("Couldn't understand cell line:\n");
 			ERROR("'%s'\n", line);
 			dirax->read_cell = 0;
-			cell_free(image->candidate_cells[image->ncells]);
+			cell_free(dirax->cur_cell);
 			return;
 		}
-		cell_set_cartesian_c(image->candidate_cells[image->ncells++],
+		cell_set_cartesian_c(dirax->cur_cell,
 		                     cx*1e-10, cy*1e-10, cz*1e-10);
 		dirax->read_cell = 0;
+
+		/* Finished reading a cell.  Time to check it... */
+		if ( check_cell(dirax->dp, image, dirax->cur_cell) ) {
+			dirax->done = 1;
+			dirax->success = 1;
+		}
+
 		return;
+
 	}
 
 	dirax->read_cell = 0;
@@ -351,8 +418,7 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax)
 
 			switch ( type ) {
 
-			case DIRAX_INPUT_LINE :
-
+				case DIRAX_INPUT_LINE :
 				block_buffer = malloc(i+1);
 				memcpy(block_buffer, dirax->rbuffer, i);
 				block_buffer[i] = '\0';
@@ -366,20 +432,17 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax)
 				endbit_length = i+2;
 				break;
 
-			case DIRAX_INPUT_PROMPT :
-
+				case DIRAX_INPUT_PROMPT :
 				dirax_send_next(image, dirax);
 				endbit_length = i+7;
 				break;
 
-			case DIRAX_INPUT_ACL :
-
+				case DIRAX_INPUT_ACL :
 				dirax_send_next(image,dirax );
 				endbit_length = i+10;
 				break;
 
-			default :
-
+				default :
 				/* Obviously, this never happens :) */
 				ERROR("Unrecognised DirAx input mode! "
 				      "I don't know how to understand DirAx\n");
@@ -456,7 +519,7 @@ static void write_drx(struct image *image)
 }
 
 
-void run_dirax(struct image *image)
+int run_dirax(struct image *image, IndexingPrivate *ipriv)
 {
 	unsigned int opts;
 	int status;
@@ -468,13 +531,13 @@ void run_dirax(struct image *image)
 	dirax = malloc(sizeof(struct dirax_data));
 	if ( dirax == NULL ) {
 		ERROR("Couldn't allocate memory for DirAx data.\n");
-		return;
+		return 0;
 	}
 
 	dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL);
 	if ( dirax->pid == -1 ) {
 		ERROR("Failed to fork for DirAx: %s\n", strerror(errno));
-		return;
+		return 0;
 	}
 	if ( dirax->pid == 0 ) {
 
@@ -505,6 +568,9 @@ void run_dirax(struct image *image)
 	dirax->read_cell = 0;
 	dirax->n_acls_tried = 0;
 	dirax->best_acl_nh = 0;
+	dirax->done = 0;
+	dirax->success = 0;
+	dirax->dp = (struct dirax_private *)ipriv;
 
 	do {
 
@@ -543,7 +609,7 @@ void run_dirax(struct image *image)
 			rval = 1;
 		}
 
-	} while ( !rval );
+	} while ( !rval && !dirax->success );
 
 	close(dirax->pty);
 	free(dirax->rbuffer);
@@ -553,5 +619,24 @@ void run_dirax(struct image *image)
 		ERROR("DirAx doesn't seem to be working properly.\n");
 	}
 
+	rval = dirax->success;
 	free(dirax);
+	return rval;
+}
+
+
+IndexingPrivate *dirax_prepare(IndexingMethod indm, UnitCell *cell,
+                               const char *filename, struct detector *det,
+                               struct beam_params *beam, float *ltl)
+{
+	struct dirax_private *dp;
+
+	dp = malloc(sizeof(struct dirax_private));
+	if ( dp == NULL ) return NULL;
+
+	dp->ltl = ltl;
+	dp->template = cell;
+	dp->indm = indm;
+
+	return (IndexingPrivate *)dp;
 }
diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h
index ce3cd4b0..09a794b0 100644
--- a/libcrystfel/src/dirax.h
+++ b/libcrystfel/src/dirax.h
@@ -3,11 +3,11 @@
  *
  * Invoke the DirAx auto-indexing program
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010,2012 Thomas White <taw@physics.org>
+ *   2010,2012-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -33,10 +33,13 @@
 #include <config.h>
 #endif
 
-#include "utils.h"
+#include "index.h"
 
+extern int run_dirax(struct image *image, IndexingPrivate *ipriv);
 
-extern void run_dirax(struct image *image);
-
+extern IndexingPrivate *dirax_prepare(IndexingMethod indm,
+                                      UnitCell *cell, const char *filename,
+                                      struct detector *det,
+                                      struct beam_params *beam, float *ltl);
 
 #endif	/* DIRAX_H */
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index d54cf43f..7cbe387a 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -175,4 +175,6 @@ extern struct imagefeature *image_feature_closest(ImageFeatureList *flist,
 extern int image_feature_count(ImageFeatureList *flist);
 extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx);
 
+extern void image_add_crystal(struct image *image, Crystal *cryst);
+
 #endif	/* IMAGE_H */
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index fd4384d7..2370e26e 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -61,7 +61,7 @@ static const char *maybes(int n)
 
 IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
                                    const char *filename, struct detector *det,
-                                   double nominal_photon_energy, float *ltl)
+                                   struct beam_params *beam, float *ltl)
 {
 	int n;
 	int nm = 0;
@@ -76,7 +76,8 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
 		switch ( indm[n] & INDEXING_METHOD_MASK ) {
 
 			case INDEXING_DIRAX :
-			iprivs[n] = NULL;
+			iprivs[n] = dirax_prepare(indm[nm], cell, filename, det,
+			                          beam, ltl);
 			break;
 
 			case INDEXING_MOSFLM :
@@ -84,7 +85,7 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
 			break;
 
 			case INDEXING_REAX :
-			iprivs[n] = reax_prepare(cell, filename, det);
+			iprivs[n] = reax_prepare(cell, filename, det, beam);
 			break;
 
 			default :
@@ -101,35 +102,30 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
 }
 
 
-void cleanup_indexing(IndexingPrivate **priv)
+void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
 {
 	int n = 0;
 
-	if ( priv == NULL ) return;  /* Nothing to do */
+	if ( indms == NULL ) return;  /* Nothing to do */
+	if ( privs == NULL ) return;  /* Nothing to do */
 
-	while ( priv[n] != NULL ) {
+	while ( indms[n] != INDEXING_NONE ) {
 
-		switch ( priv[n]->indm & INDEXING_METHOD_MASK ) {
+		switch ( indms[n] & INDEXING_METHOD_MASK ) {
 
 			case INDEXING_NONE :
-			free(priv[n]);
-			break;
-
 			case INDEXING_DIRAX :
-			free(priv[n]);
-			break;
-
 			case INDEXING_MOSFLM :
-			free(priv[n]);
+			/* No cleanup */
 			break;
 
 			case INDEXING_REAX :
-			reax_cleanup(priv[n]);
+			reax_cleanup(privs[n]);
 			break;
 
 			default :
 			ERROR("Don't know how to clean up indexing method %i\n",
-			      priv[n]->indm);
+			      indms[n]);
 			break;
 
 		}
@@ -162,19 +158,43 @@ void map_all_peaks(struct image *image)
 }
 
 
-static void try_indexer(struct image *image, IndexingMethod indm,
+/* Return non-zero for "success" */
+static int try_indexer(struct image *image, IndexingMethod indm,
                         IndexingPrivate *ipriv)
 {
+	switch ( indm & INDEXING_METHOD_MASK ) {
+
+		case INDEXING_NONE :
+		break;
+
+		case INDEXING_DIRAX :
+		return run_dirax(image, ipriv);
+		break;
+
+		case INDEXING_MOSFLM :
+		return run_mosflm(image, ipriv);
+		break;
+
+		case INDEXING_REAX :
+		return reax_index(image, ipriv);
+		break;
+
+		default :
+		ERROR("Unrecognised indexing method: %i\n", indm);
+		break;
+
+	}
+
+	return 0;
 }
 
 
 void index_pattern(struct image *image,
                    IndexingMethod *indms, IndexingPrivate **iprivs)
 {
-	int i;
 	int n = 0;
 
-	if ( indm == NULL ) return;
+	if ( indms == NULL ) return;
 
 	map_all_peaks(image);
 	image->crystals = NULL;
@@ -182,7 +202,7 @@ void index_pattern(struct image *image,
 
 	while ( indms[n] != INDEXING_NONE ) {
 
-		if ( try_indexer(image, indms[n], iprivs[i]) ) break;
+		if ( try_indexer(image, indms[n], iprivs[n]) ) break;
 		n++;
 
 	}
@@ -216,7 +236,8 @@ static IndexingMethod set_bad(IndexingMethod a)
 /* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */
 static IndexingMethod set_axes(IndexingMethod a)
 {
-	return (a & ~INDEXING_CHECK_COMBINATIONS) | INDEXING_CHECK_CELL_AXES;
+	return (a & ~INDEXING_CHECK_CELL_COMBINATIONS)
+	          | INDEXING_CHECK_CELL_AXES;
 }
 
 
@@ -234,7 +255,7 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell)
 	n = assplode(str, ",-", &methods, ASSPLODE_NONE);
 	list = malloc((n+1)*sizeof(IndexingMethod));
 
-	*nmeth = -1;  /* So that the first method is #0 */
+	nmeth = -1;  /* So that the first method is #0 */
 	for ( i=0; i<n; i++ ) {
 
 		if ( strcmp(methods[i], "dirax") == 0) {
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index d5b01b43..e0dcb8e0 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -87,15 +87,11 @@ extern IndexingMethod *build_indexer_list(const char *str, int *need_cell);
 extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
                                           const char *filename,
                                           struct detector *det,
-                                          double nominal_photon_energy,
-                                          float *ltl);
+                                          struct beam_params *beam, float *ltl);
 
 extern void index_pattern(struct image *image,
                           IndexingMethod *indms, IndexingPrivate **iprivs);
 
-extern void cleanup_indexing(IndexingPrivate **priv);
-
-/* For indexing methods to call */
-extern IndexingPrivate *indexing_private(IndexingMethod indm);
+extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs);
 
 #endif	/* INDEX_H */
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 548e3834..8b300012 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -80,6 +80,7 @@
 #include "mosflm.h"
 #include "utils.h"
 #include "peaks.h"
+#include "cell-utils.h"
 
 
 #define MOSFLM_VERBOSE 0
@@ -92,6 +93,14 @@ typedef enum {
 } MOSFLMInputType;
 
 
+
+struct mosflm_private {
+	IndexingMethod          indm;
+	float                   *ltl;
+	UnitCell                *template;
+};
+
+
 struct mosflm_data {
 
 	/* MOSFLM auto-indexing low-level stuff */
@@ -107,10 +116,53 @@ struct mosflm_data {
 	char                    sptfile[128];
 	int                     step;
 	int                     finished_ok;
-	UnitCell                *target_cell;
+	int                     done;
+	int                     success;
+
+	struct mosflm_private  *mp;
 
 };
 
+static int check_cell(struct mosflm_private *mp, struct image *image,
+                      UnitCell *cell)
+{
+	UnitCell *out;
+	Crystal *cr;
+
+	if ( mp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) {
+
+		out = match_cell(cell, mp->template, 0, mp->ltl, 1);
+		if ( out == NULL ) return 0;
+
+	} else if ( mp->indm & INDEXING_CHECK_CELL_AXES ) {
+
+		out = match_cell(cell, mp->template, 0, mp->ltl, 0);
+		if ( out == NULL ) return 0;
+
+	} else {
+		out = cell_new_from_cell(cell);
+	}
+
+	cr = crystal_new();
+	if ( cr == NULL ) {
+		ERROR("Failed to allocate crystal.\n");
+		return 0;
+	}
+
+	crystal_set_cell(cr, cell);
+
+	if ( mp->indm & INDEXING_CHECK_PEAKS ) {
+		if ( !peak_sanity_check(image, &cr, 1) ) {
+			crystal_free(cr);  /* Frees the cell as well */
+			return 0;
+		}
+	}
+
+	image_add_crystal(image, cr);
+
+	return 1;
+}
+
 
 static void mosflm_parseline(const char *line, struct image *image,
                              struct mosflm_data *dirax)
@@ -130,7 +182,8 @@ static void mosflm_parseline(const char *line, struct image *image,
 }
 
 
-static int read_newmat(const char *filename, struct image *image)
+static int read_newmat(struct mosflm_data *mosflm, const char *filename,
+                       struct image *image)
 {
 	FILE * fh;
 	float asx, asy, asz;
@@ -138,6 +191,7 @@ static int read_newmat(const char *filename, struct image *image)
 	float csx, csy, csz;
 	int n;
 	double c;
+	UnitCell *cell;
 
 	fh = fopen(filename, "r");
 	if ( fh == NULL ) {
@@ -155,18 +209,21 @@ static int read_newmat(const char *filename, struct image *image)
 	/* MOSFLM "A" matrix is multiplied by lambda, so fix this */
 	c = 1.0/image->lambda;
 
-	image->candidate_cells[0] = cell_new();
+	cell = cell_new();
 
 	/* The relationship between the coordinates in the spot file and the
 	 * resulting matrix is diabolically complicated.  This transformation
 	 * seems to work, but is not derived by working through all the
 	 * transformations. */
-	cell_set_reciprocal(image->candidate_cells[0],
+	cell_set_reciprocal(cell,
 	                    -asy*c, -asz*c, asx*c,
 	                    -bsy*c, -bsz*c, bsx*c,
 	                    -csy*c, -csz*c, csx*c);
 
-        image->ncells = 1;
+	if ( check_cell(mosflm->mp, image, cell) ) {
+		mosflm->success = 1;
+		mosflm->done = 1;
+	}
 
         return 0;
 }
@@ -352,9 +409,9 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 		break;
 
 		case 2 :
-		if ( mosflm->target_cell != NULL ) {
+		if ( mosflm->mp->template != NULL ) {
 			const char *symm;
-			symm = spacegroup_for_lattice(mosflm->target_cell);
+			symm = spacegroup_for_lattice(mosflm->mp->template);
 			snprintf(tmp, 255, "SYMM %s\n", symm);
 			mosflm_sendline(tmp, mosflm);
 		} else {
@@ -527,7 +584,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
 }
 
 
-void run_mosflm(struct image *image, UnitCell *cell)
+int run_mosflm(struct image *image, IndexingPrivate *ipriv)
 {
 	struct mosflm_data *mosflm;
 	unsigned int opts;
@@ -537,11 +594,9 @@ void run_mosflm(struct image *image, UnitCell *cell)
 	mosflm = malloc(sizeof(struct mosflm_data));
 	if ( mosflm == NULL ) {
 		ERROR("Couldn't allocate memory for MOSFLM data.\n");
-		return;
+		return 0;
 	}
 
-	mosflm->target_cell = cell;
-
 	snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id);
 	write_img(image, mosflm->imagefile); /* Dummy image */
 
@@ -556,7 +611,7 @@ void run_mosflm(struct image *image, UnitCell *cell)
 	if ( mosflm->pid == -1 ) {
 		ERROR("Failed to fork for MOSFLM: %s\n", strerror(errno));
 		free(mosflm);
-		return;
+		return 0;
 	}
 	if ( mosflm->pid == 0 ) {
 
@@ -584,6 +639,9 @@ void run_mosflm(struct image *image, UnitCell *cell)
 
 	mosflm->step = 1;	/* This starts the "initialisation" procedure */
 	mosflm->finished_ok = 0;
+	mosflm->mp = (struct mosflm_private *)ipriv;
+	mosflm->done = 0;
+	mosflm->success = 0;
 
 	do {
 
@@ -632,8 +690,27 @@ void run_mosflm(struct image *image, UnitCell *cell)
 		ERROR("MOSFLM doesn't seem to be working properly.\n");
 	} else {
 		/* Read the mosflm NEWMAT file and get cell if found */
-		read_newmat(mosflm->newmatfile, image);
+		read_newmat(mosflm, mosflm->newmatfile, image);
 	}
 
+	rval = mosflm->success;
 	free(mosflm);
+	return rval;
+}
+
+
+IndexingPrivate *mosflm_prepare(IndexingMethod indm, UnitCell *cell,
+                               const char *filename, struct detector *det,
+                               struct beam_params *beam, float *ltl)
+{
+	struct mosflm_private *mp;
+
+	mp = malloc(sizeof(struct mosflm_private));
+	if ( mp == NULL ) return NULL;
+
+	mp->ltl = ltl;
+	mp->template = cell;
+	mp->indm = indm;
+
+	return (IndexingPrivate *)mp;
 }
diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h
index ba1eb73d..1d296cb7 100644
--- a/libcrystfel/src/mosflm.h
+++ b/libcrystfel/src/mosflm.h
@@ -35,10 +35,10 @@
 #include <config.h>
 #endif
 
-#include "utils.h"
+#include "index.h"
 
 
-extern void run_mosflm(struct image *image, UnitCell *cell);
+extern int run_mosflm(struct image *image, IndexingPrivate *ipriv);
 
 
 #endif	/* MOSFLM_H */
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index c350388b..08595ef5 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -594,7 +594,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 }
 
 
-int peak_sanity_check(struct image *image)
+int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst)
 {
 	int n_feat = 0;
 	int n_sane = 0;
@@ -616,13 +616,13 @@ int peak_sanity_check(struct image *image)
 		/* Reciprocal space position of found peak */
 		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
 
-		for ( j=0; j<image->n_crystals; j++ ) {
+		for ( j=0; j<n_cryst; j++ ) {
 
 			double ax, ay, az;
 			double bx, by, bz;
 			double cx, cy, cz;
 
-			cell_get_cartesian(crystal_get_cell(image->crystals[j]),
+			cell_get_cartesian(crystal_get_cell(crystals[j]),
 			                   &ax, &ay, &az,
 			                   &bx, &by, &bz,
 			                   &cx, &cy, &cz);
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 40136eca..b00edf00 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -47,7 +47,8 @@ extern void integrate_reflections(struct image *image,
                                   double ir_inn, double ir_mid, double ir_out,
                                   int integrate_saturated);
 
-extern int peak_sanity_check(struct image *image);
+extern int peak_sanity_check(struct image *image, Crystal **crystals,
+                             int n_cryst);
 
 extern void estimate_resolution(RefList *list, UnitCell *cell,
                                 double *min, double *max);
diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h
index 01dfe649..af383e29 100644
--- a/libcrystfel/src/reax.h
+++ b/libcrystfel/src/reax.h
@@ -34,16 +34,24 @@
 #endif
 
 #include "cell.h"
+#include "beam-parameters.h"
+#include "detector.h"
 
 #ifdef HAVE_FFTW
 
-extern IndexingPrivate *reax_prepare(void);
+extern IndexingPrivate *reax_prepare(UnitCell *cell, const char *filename,
+                                     struct detector *det,
+                                     struct beam_params *beam);
+
 extern void reax_cleanup(IndexingPrivate *pp);
-extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell);
+
+extern int reax_index(struct image *image, IndexingPrivate *p);
 
 #else /* HAVE_FFTW */
 
-static IndexingPrivate *reax_prepare()
+static IndexingPrivate *reax_prepare(UnitCell *cell, const char *filename,
+                                     struct detector *det,
+                                     struct beam_params *beam)
 {
 	return NULL;
 }
@@ -52,7 +60,7 @@ static void reax_cleanup(IndexingPrivate *pp)
 {
 }
 
-static void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell)
+static void reax_index(struct image *image, IndexingPrivate *p)
 {
 }
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index accfa788..01377de3 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -243,7 +243,6 @@ int main(int argc, char *argv[])
 	struct index_args iargs;
 	struct beam_params *beam = NULL;
 	char *element = NULL;
-	double nominal_photon_energy;
 	int stream_flags = STREAM_INTEGRATED;
 	char *hdf5_peak_path = NULL;
 	struct copy_hdf5_field *copyme;
@@ -556,7 +555,14 @@ int main(int argc, char *argv[])
 	}
 
 	if ( geometry == NULL ) {
-		ERROR("You need to specify a geometry file with --geometry\n");
+		ERROR("You need to provide a geometry file (please read the"
+		      " manual for more details).\n");
+		return 1;
+	}
+
+	if ( beam == NULL ) {
+		ERROR("You need to provide a beam parameters file (please read"
+		      " the manual for more details).\n");
 		return 1;
 	}
 
@@ -582,17 +588,6 @@ int main(int argc, char *argv[])
 
 	write_stream_header(ofh, argc, argv);
 
-	if ( beam != NULL ) {
-		nominal_photon_energy = beam->photon_energy;
-	} else {
-		STATUS("No beam parameters file was given, so I'm taking the"
-		       " nominal photon energy to be 2 keV.\n");
-		ERROR("I'm also going to assume 1 ADU per photon, which is");
-		ERROR(" almost certainly wrong.  Peak sigmas will be"
-		      " incorrect.\n");
-		nominal_photon_energy = 2000.0;
-	}
-
 	/* Get first filename and use it to set up the indexing */
 	prepare_line = malloc(1024);
 	rval = fgets(prepare_line, 1023, fh);
@@ -613,8 +608,8 @@ int main(int argc, char *argv[])
 
 	/* Prepare the indexer */
 	if ( indm != NULL ) {
-		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
-		                         nominal_photon_energy);
+		ipriv = prepare_indexing(indm, cell, prepare_filename,
+		                         det, beam);
 		if ( ipriv == NULL ) {
 			ERROR("Failed to prepare indexing.\n");
 			return 1;