Merge branch 'master' into taketwo

33 files changed, 2928 insertions, 602 deletions

diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 7bb08685..fef2e11f 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -11,7 +11,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
                          src/cell-utils.c src/integer_matrix.c src/crystal.c \
                          src/xds.c src/integration.c src/predict-refine.c \
                          src/histogram.c src/events.c src/felix.c \
-                         src/taketwo.c
+                         src/taketwo.c src/peakfinder8.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/asdf.c
@@ -32,7 +32,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \
                                  src/xds.h src/predict-refine.h \
                                  src/integration.h src/histogram.h \
                                  src/events.h src/asdf.h src/felix.h \
-                                 src/taketwo.h
+                                 src/taketwo.h src/peakfinder8.h
 
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
 AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index 49917ad9..0c43fe7a 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -4,12 +4,12 @@
  * Alexandra's Superior Direction Finder, or
  * Algorithm Similar to DirAx, FFT-based
  *
- * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
- *   2015      Thomas White <taw@physics.org>
+ *   2015,2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -28,8 +28,6 @@
  *
  */
 
-#define _ISOC99_SOURCE
-
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
@@ -1103,7 +1101,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
 }
 
 
-int run_asdf(struct image *image, IndexingPrivate *ipriv)
+int run_asdf(struct image *image, void *ipriv)
 {
 	int i, j;
 
@@ -1204,8 +1202,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
 }
 
 
-IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
-                               struct detector *det, float *ltl)
+void *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
+                   struct detector *det, float *ltl)
 {
 	struct asdf_private *dp;
 	int need_cell = 0;
@@ -1234,11 +1232,11 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
 
 	dp->fftw = fftw_vars_new();
 
-	return (IndexingPrivate *)dp;
+	return (void *)dp;
 }
 
 
-void asdf_cleanup(IndexingPrivate *pp)
+void asdf_cleanup(void *pp)
 {
 	struct asdf_private *p;
 	p = (struct asdf_private *)pp;
diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h
index 1e492a6f..f130d63d 100644
--- a/libcrystfel/src/asdf.h
+++ b/libcrystfel/src/asdf.h
@@ -4,12 +4,12 @@
  * Alexandra's Superior Direction Finder, or
  * Algorithm Similar to DirAx, FFT-based
  *
- * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
- *   2015      Thomas White <taw@physics.org>
+ *   2015,2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -43,33 +43,33 @@ extern "C" {
 
 #ifdef HAVE_FFTW
 
-extern int run_asdf(struct image *image, IndexingPrivate *ipriv);
+extern int run_asdf(struct image *image, void *ipriv);
 
-extern IndexingPrivate *asdf_prepare(IndexingMethod *indm,
-                                     UnitCell *cell, struct detector *det,
-                                     float *ltl);
+extern void *asdf_prepare(IndexingMethod *indm,
+                          UnitCell *cell, struct detector *det,
+                          float *ltl);
 
-extern void asdf_cleanup(IndexingPrivate *pp);
+extern void asdf_cleanup(void *pp);
 
 #else /* HAVE_FFTW */
 
-int run_asdf(struct image *image, IndexingPrivate *ipriv)
+int run_asdf(struct image *image, void *ipriv)
 {
 	ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
 	return 0;
 }
 
 
-IndexingPrivate *asdf_prepare(IndexingMethod *indm,
-                              UnitCell *cell, struct detector *det,
-                              float *ltl)
+void *asdf_prepare(IndexingMethod *indm,
+                   UnitCell *cell, struct detector *det,
+                   float *ltl)
 {
 	ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
 	ERROR("To use asdf indexing, recompile with FFTW.\n");
 	return NULL;
 }
 
-void asdf_cleanup(IndexingPrivate *pp)
+void asdf_cleanup(void *pp)
 {
 }
 
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 8fa3b894..6610abc3 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -272,6 +272,12 @@ void cell_print(UnitCell *cell)
 		STATUS("c* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n",
 		       csx, csy, csz, modulus(csx, csy, csz));
 
+		STATUS("alpha* = %6.2f deg, beta* = %6.2f deg, "
+		       "gamma* = %6.2f deg\n",
+		       rad2deg(angle_between(bsx, bsy, bsz, csx, csy, csz)),
+		       rad2deg(angle_between(asx, asy, asz, csx, csy, csz)),
+		       rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz)));
+
 		STATUS("Cell representation is %s.\n", cell_rep(cell));
 
 	} else {
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index ec591e24..3de61073 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -460,10 +460,12 @@ int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
 		case CELL_REP_RECIP:
 		/* Convert reciprocal -> crystallographic.
                  * Start by converting reciprocal -> cartesian */
-		cell_invert(cell->axs, cell->ays, cell->azs,
-		            cell->bxs, cell->bys, cell->bzs,
-		            cell->cxs, cell->cys, cell->czs,
-		            &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+		if ( cell_invert(cell->axs, cell->ays, cell->azs,
+		                 cell->bxs, cell->bys, cell->bzs,
+		                         cell->cxs, cell->cys, cell->czs,
+		                         &ax, &ay, &az,
+		                         &bx, &by, &bz,
+		                         &cx, &cy, &cz) ) return 1;
 
 		/* Now convert cartesian -> crystallographic */
 		*a = modulus(ax, ay, az);
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 42ffbe40..e3b351ff 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -31,8 +31,6 @@
  *
  */
 
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
 #include <stdlib.h>
 #include <math.h>
 #include <stdio.h>
@@ -530,37 +528,15 @@ int panel_number(struct detector *det, struct panel *p)
 }
 
 
-void fill_in_values(struct detector *det, struct hdfile *f, struct event* ev)
+void adjust_centering_for_rail(struct panel *p)
 {
-	int i;
-
-	for ( i=0; i<det->n_panels; i++ ) {
-
-		double offs;
-		struct panel *p = &det->panels[i];
-
-		if ( p->clen_from != NULL ) {
-
-			double val;
-			int r;
+	double offs;
 
-			r = hdfile_get_value(f, p->clen_from, ev, &val,
-			                     H5T_NATIVE_DOUBLE);
-			if ( r ) {
-				ERROR("Failed to read '%s'\n", p->clen_from);
-			} else {
-				p->clen = val * 1.0e-3;
-			}
-
-		}
-
-                /* Offset in +z direction from calibrated clen to actual */
-		offs = p->clen - p->clen_for_centering;
-		p->cnx += p->rail_x * offs;
-		p->cny += p->rail_y * offs;
-		p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs;
-
-	}
+	/* Offset in +z direction from calibrated clen to actual */
+	offs = p->clen - p->clen_for_centering;
+	p->cnx += p->rail_x * offs;
+	p->cny += p->rail_y * offs;
+	p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs;
 }
 
 
@@ -911,7 +887,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
 		panel->adu_per_eV = atof(val);
 	} else if ( strcmp(key, "adu_per_photon") == 0 ) {
 		panel->adu_per_photon = atof(val);
-		STATUS("got adu per photon: %s\n", val);
 	} else if ( strcmp(key, "rigid_group") == 0 ) {
 		add_to_rigid_group(find_or_add_rg(det, val), panel);
 	} else if ( strcmp(key, "clen") == 0 ) {
@@ -1079,12 +1054,15 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam,
 
 	} else if ( strcmp(key, "photon_energy") == 0 ) {
 		if ( beam != NULL ) {
-			if ( strncmp(val, "/", 1) == 0 ) {
+			double v;
+			char *end;
+			v = strtod(val, &end);
+			if ( (val[0] != '\0') && (end[0] == '\0') ) {
+				beam->photon_energy = v;
+				beam->photon_energy_from = NULL;
+			} else {
 				beam->photon_energy = 0.0;
 				beam->photon_energy_from = strdup(val);
-			} else {
-				beam->photon_energy = atof(val);
-				beam->photon_energy_from = NULL;
 			}
 		}
 
@@ -1204,7 +1182,7 @@ struct detector *get_detector_geometry_2(const char *filename,
 	int i;
 	int rgi, rgci;
 	int reject = 0;
-	int path_dim;
+	int path_dim, mask_path_dim;
 	int dim_dim;
 	int dim_reject = 0;
 	int dim_dim_reject = 0;
@@ -1389,6 +1367,7 @@ struct detector *get_detector_geometry_2(const char *filename,
 
 	}
 
+	mask_path_dim = -1;
 	for ( i=0; i<det->n_panels; i++ ) {
 
 		int panel_mask_dim = 0;
@@ -1398,8 +1377,7 @@ struct detector *get_detector_geometry_2(const char *filename,
 
 			next_instance = det->panels[i].mask;
 
-			while(next_instance)
-			{
+			while ( next_instance ) {
 				next_instance = strstr(next_instance, "%");
 				if ( next_instance != NULL ) {
 					next_instance += 1*sizeof(char);
@@ -1407,18 +1385,29 @@ struct detector *get_detector_geometry_2(const char *filename,
 				}
 			}
 
-			if ( panel_mask_dim != path_dim ) {
-				dim_reject = 1;
+			if ( mask_path_dim == -1 ) {
+				mask_path_dim = panel_mask_dim;
+			} else {
+				if ( panel_mask_dim != mask_path_dim ) {
+					dim_reject = 1;
+				}
 			}
+
 		}
 	}
 
-	if ( dim_reject ==  1) {
+	if ( dim_reject ==  1 ) {
 		ERROR("All panels' data and mask entries must have the same "
 		      "number of placeholders\n");
 		reject = 1;
 	}
 
+	if ( mask_path_dim > path_dim ) {
+		ERROR("Number of placeholders in mask cannot be larger than "
+		      "for data\n");
+		reject = 1;
+	}
+
 	det->path_dim = path_dim;
 
 	dim_dim_reject = 0;
@@ -1993,6 +1982,11 @@ double largest_q(struct image *image)
 	struct rvec q;
 	double tt;
 
+	if ( image->det == NULL ) {
+		ERROR("No detector geometry. assuming detector is infinite!\n");
+		return INFINITY;
+	}
+
 	q = get_q_for_panel(image->det->furthest_out_panel,
 	                    image->det->furthest_out_fs,
 	                    image->det->furthest_out_ss,
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index 04e3c1ec..a2be2b47 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -3,12 +3,12 @@
  *
  * Detector properties
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2009-2016 Thomas White <taw@physics.org>
+ *   2009-2017 Thomas White <taw@physics.org>
  *   2011-2012 Richard Kirian <rkirian@asu.edu>
  *   2014      Valerio Mariani
  *   2011      Andrew Aquila
@@ -42,7 +42,6 @@ struct rg_collection;
 struct detector;
 struct panel;
 struct badregion;
-struct detector;
 struct beam_params;
 struct hdfile;
 struct event;
@@ -80,6 +79,47 @@ struct rg_collection
 };
 
 
+/**
+ * panel:
+ *  @name: Text name for the panel (fixed length array)
+ *  @cnx: Location of corner, in pixels, x coordinate
+ *  @cny: Location of corner, in pixels, y coordinate
+ *  @coffset: The offset to be applied from @clen (which may come from elsewhere)
+ *  @clen: The distance from the interaction point to the corner of the first pixel
+ *  @clen_from: Location to get @clen from, e.g. from HDF5 file
+ *  @mask: Location of mask data
+ *  @mask_file: Filename for mask data
+ *  @satmap: Location of per-pixel saturation map
+ *  @satmap_file: Filename for saturation map
+ *  @res: Resolution of panel in pixels per metre
+ *  @badrow: Readout direction (for filtering out clusters of peaks)
+ *  @no_index: Non-zero if panel is entirely "bad"
+ *  @adu_per_photon: Number of detector intensity units per photon
+ *  @adu_per_eV: Number of detector intensity units per eV of photon energy
+ *  @max_adu: Saturation value
+ *  @dim_structure: Dimension structure
+ *  @fsx: Real-space x-direction of data fast-scan direction
+ *  @fsy: Real-space y-direction of data fast-scan direction
+ *  @fsz: Real-space z-direction of data fast-scan direction
+ *  @ssx: Real-space x-direction of data slow-scan direction
+ *  @ssy: Real-space y-direction of data slow-scan direction
+ *  @ssz: Real-space z-direction of data slow-scan direction
+ *  @rail_x: x direction of camera length "rail"
+ *  @rail_y: y direction of camera length "rail"
+ *  @rail_z: z direction of camera length "rail"
+ *  @clen_for_centering: Value of clen (without coffset) at which beam is centered
+ *  @xfs: Data fast-scan direction of real-space x-direction
+ *  @yfs: Data fast-scan direction of real-space y-direction
+ *  @xss: Data slow-scan direction of real-space x-direction
+ *  @yss: Data slow-scan direction of real-space y-direction
+ *  @orig_min_fs: Minimum fs coordinate of data in file
+ *  @orig_max_fs: Maximum fs coordinate of data in file
+ *  @orig_min_ss: Minimum ss coordinate of data in file (inclusive)
+ *  @orig_max_ss: Maximum ss coordinate of data in file (inclusive)
+ *  @data: Location of data in file
+ *  @w: Width of panel
+ *  @h: Height of panel
+ */
 struct panel
 {
         char     name[1024];  /* Name for this panel */
@@ -224,9 +264,8 @@ extern void get_pixel_extents(struct detector *det,
                               double *min_x, double *min_y,
                               double *max_x, double *max_y);
 
-extern void fill_in_values(struct detector *det, struct hdfile *f,
-                           struct event* ev);
 extern void fill_in_adu(struct image *image);
+extern void adjust_centering_for_rail(struct panel *p);
 
 extern int panel_is_in_rigid_group(const struct rigid_group *rg,
                                    struct panel *p);
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 19f35696..e9466a24 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -3,11 +3,11 @@
  *
  * Invoke the DirAx auto-indexing program
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2014,2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -532,7 +532,7 @@ static void write_drx(struct image *image)
 }
 
 
-int run_dirax(struct image *image, IndexingPrivate *ipriv)
+int run_dirax(struct image *image, void *ipriv)
 {
 	unsigned int opts;
 	int status;
@@ -638,8 +638,8 @@ int run_dirax(struct image *image, IndexingPrivate *ipriv)
 }
 
 
-IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
-                               struct detector *det, float *ltl)
+void *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
+                    struct detector *det, float *ltl)
 {
 	struct dirax_private *dp;
 	int need_cell = 0;
@@ -670,7 +670,7 @@ IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
 }
 
 
-void dirax_cleanup(IndexingPrivate *pp)
+void dirax_cleanup(void *pp)
 {
 	struct dirax_private *p;
 	p = (struct dirax_private *)pp;
diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h
index 9776f3f0..96ba6dbc 100644
--- a/libcrystfel/src/dirax.h
+++ b/libcrystfel/src/dirax.h
@@ -3,11 +3,11 @@
  *
  * Invoke the DirAx auto-indexing program
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010,2012-2014 Thomas White <taw@physics.org>
+ *   2010,2012-2014,2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -39,13 +39,12 @@
 extern "C" {
 #endif
 
-extern int run_dirax(struct image *image, IndexingPrivate *ipriv);
+extern int run_dirax(struct image *image, void *ipriv);
 
-extern IndexingPrivate *dirax_prepare(IndexingMethod *indm,
-                                      UnitCell *cell, struct detector *det,
-                                      float *ltl);
+extern void *dirax_prepare(IndexingMethod *indm,
+                           UnitCell *cell, struct detector *det, float *ltl);
 
-extern void dirax_cleanup(IndexingPrivate *pp);
+extern void dirax_cleanup(void *pp);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c
index 25771a69..8e4eb861 100644
--- a/libcrystfel/src/events.c
+++ b/libcrystfel/src/events.c
@@ -3,10 +3,11 @@
  *
  * Event properties
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
+ *   2017      Thomas White
  *   2014      Valerio Mariani
  *
  * This file is part of CrystFEL.
@@ -26,10 +27,8 @@
  *
  */
 
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
-
 #include "events.h"
+#include "utils.h"
 
 #include <hdf5.h>
 #include <string.h>
@@ -585,46 +584,38 @@ char *retrieve_full_path(struct event *ev, const char *data)
 {
 	int ei ;
 	char *return_value;
+	char *pholder;
 
 	return_value = strdup(data);
+	pholder = strstr(return_value,"%");
+	ei = 0;
 
-	for ( ei=0; ei<ev->path_length; ei++ ) {
+	while ( pholder != NULL ) {
 
 		char *tmp;
-		tmp = event_path_placeholder_subst(ev->path_entries[ei],
-		                                   return_value);
 
-		free(return_value);
-		return_value = tmp;
-
-	}
-
-	return return_value;
-
-}
-
-
-char *partial_event_substitution(struct event *ev, const char *data)
-{
-	int ei ;
-	char *return_value;
-	char *pholder;
-
-	return_value = strdup(data);
-	pholder = strstr(return_value,"%");
-	ei = 0;
+		/* Check we have enough things to put in the placeholders */
+		if ( ei >= ev->path_length ) {
+			ERROR("Too many placeholders ('%%') in location.\n");
+			return NULL;
+		}
 
-	while( pholder != NULL) {
+		/* Substitute one placeholder */
+		tmp = event_path_placeholder_subst(ev->path_entries[ei++],
+		                                         return_value);
 
-		char *tmp_subst_data;
+		if ( tmp == NULL ) {
+			ERROR("Couldn't substitute placeholder\n");
+			return NULL;
+		}
 
-		tmp_subst_data = event_path_placeholder_subst(ev->path_entries[ei],
-                                                      return_value);
+		/* Next time, we will substitute the next part of the path into
+		 * the partially substituted string */
 		free(return_value);
-		return_value = strdup(tmp_subst_data);
-		free(tmp_subst_data);
+		return_value = tmp;
+
 		pholder = strstr(return_value, "%");
-		ei += 1;
+
 	}
 
 	return return_value;
diff --git a/libcrystfel/src/events.h b/libcrystfel/src/events.h
index 7f9c6731..4f717209 100644
--- a/libcrystfel/src/events.h
+++ b/libcrystfel/src/events.h
@@ -78,7 +78,6 @@ extern char *get_event_string(struct event *ev);
 extern struct event *get_event_from_event_string(const char *ev_string);
 extern char *event_path_placeholder_subst(const char *ev_name,
                                           const char *data);
-extern char *partial_event_substitution(struct event *ev, const char *data);
 extern char *retrieve_full_path(struct event *ev, const char *data);
 
 
diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c
index 02d523c6..b531e3af 100644
--- a/libcrystfel/src/felix.c
+++ b/libcrystfel/src/felix.c
@@ -3,8 +3,8 @@
  *
  * Invoke Felix for multi-crystal autoindexing.
  *
- * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2015-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2015 Thomas White <taw@physics.org>
@@ -658,9 +658,8 @@ static void parse_options(const char *options, struct felix_private *gp)
 	free(option);
 }
 
-IndexingPrivate *felix_prepare(IndexingMethod *indm, UnitCell *cell,
-                                      struct detector *det, float *ltl,
-                                      const char *options)
+void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+                    struct detector *det, float *ltl, const char *options)
 {
 	struct felix_private *gp;
 
diff --git a/libcrystfel/src/felix.h b/libcrystfel/src/felix.h
index 40771933..c06e963a 100644
--- a/libcrystfel/src/felix.h
+++ b/libcrystfel/src/felix.h
@@ -3,12 +3,12 @@
  *
  * Invoke Felix for multi-crystal autoindexing
  *
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
- *   2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *   2010-2013,2017 Thomas White <taw@physics.org>
+ *   2013-2014      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -36,11 +36,9 @@
 
 #include "cell.h"
 
-extern IndexingPrivate *felix_prepare(IndexingMethod *indm,
-                                             UnitCell *cell,
-                                             struct detector *det,
-                                             float *ltl,
-                                             const char *options);
+extern void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+                           struct detector *det, float *ltl,
+                           const char *options);
 
 extern void felix_cleanup(IndexingPrivate *pp);
 
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index 6542b360..13339506 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -44,6 +44,19 @@
 #include "utils.h"
 
 
+/**
+ * SECTION:hdf5-file
+ * @short_description: HDF5 file handling
+ * @title: HDF5 file handling
+ * @section_id:
+ * @see_also:
+ * @include: "hdf5-file.h"
+ * @Image:
+ *
+ * Routines for accessing HDF5 files.
+ **/
+
+
 struct hdf5_write_location {
 
 	const char      *location;
@@ -121,8 +134,7 @@ struct hdfile *hdfile_open(const char *filename)
 }
 
 
-int hdfile_set_image(struct hdfile *f, const char *path,
-                     struct panel *p)
+int hdfile_set_image(struct hdfile *f, const char *path)
 {
 	f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT);
 	if ( f->dh < 0 ) {
@@ -320,9 +332,34 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line)
 }
 
 
-/* Get peaks from HDF5, in "CXI format" (as in "CXIDB") */
-int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
-                  struct filename_plus_event *fpe)
+/**
+ * get_peaks_cxi_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5, in "CXI format" (as in "CXIDB").  The data should be in
+ * a set of arrays under @p.  The number of peaks should be in a 1D array at
+ * @p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at
+ * @p/peakXPosRaw and @p/peakYPosRaw respectively (sorry about the naming).  The
+ * first dimension of these arrays should be the event number (as given by
+ * @fpe).  The intensities are expected to be at @p/peakTotalIntensity in a
+ * similar 2D array.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry').  The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array.  In this case, the peak coordinates should
+ * have 0.5 added to them.  This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+                    struct filename_plus_event *fpe, int half_pixel_shift)
 {
 	char path_n[1024];
 	char path_x[1024];
@@ -338,6 +375,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 	float *buf_y;
 	float *buf_i;
 
+	double peak_offset = half_pixel_shift ? 0.5 : 0.0;
+
 	if ( (fpe != NULL) && (fpe->ev != NULL)
 	  && (fpe->ev->dim_entries != NULL) )
 	{
@@ -375,8 +414,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 		float fs, ss, val;
 		struct panel *p;
 
-		fs = buf_x[pk];
-		ss = buf_y[pk];
+		fs = buf_x[pk] + peak_offset;
+		ss = buf_y[pk] + peak_offset;
 		val = buf_i[pk];
 
 		p = find_orig_panel(image->det, fs, ss);
@@ -395,7 +434,53 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 }
 
 
-int get_peaks(struct image *image, struct hdfile *f, const char *p)
+/**
+ * get_peaks_cxi:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions.  Use get_peaks_cxi_2() instead.
+ *
+ * This function is equivalent to get_peaks_cxi_2(@image, @f, @p, @fpe, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
+                  struct filename_plus_event *fpe)
+{
+	return get_peaks_cxi_2(image, f, p, fpe, 1);
+}
+
+
+/**
+ * get_peaks_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5.  The peak list should be located at @p in the HDF5 file,
+ * a 2D array where the first dimension equals the number of peaks and second
+ * dimension is three.  The first two columns contain the fast scan and slow
+ * scan coordinates, respectively, of the peaks.  The third column contains the
+ * intensities.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry').  The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array.  In this case, the peak coordinates should
+ * have 0.5 added to them.  This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+                int half_pixel_shift)
 {
 	hid_t dh, sh;
 	hsize_t size[2];
@@ -405,6 +490,7 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 	herr_t r;
 	int tw;
 	char *np;
+	double peak_offset = half_pixel_shift ? 0.5 : 0.0;
 
 	if ( image->event != NULL ) {
 		np = retrieve_full_path(image->event, p);
@@ -473,8 +559,8 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 		float fs, ss, val;
 		struct panel *p;
 
-		fs = buf[tw*i+0];
-		ss = buf[tw*i+1];
+		fs = buf[tw*i+0] + peak_offset;
+		ss = buf[tw*i+1] + peak_offset;
 		val = buf[tw*i+2];
 
 		p = find_orig_panel(image->det, fs, ss);
@@ -499,6 +585,26 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 }
 
 
+/**
+ * get_peaks:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions.  Use get_peaks_2() instead.
+ *
+ * This function is equivalent to get_peaks_2(@image, @f, @p, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks(struct image *image, struct hdfile *f, const char *p)
+{
+	return get_peaks_2(image, f, p, 1);
+}
+
+
 static void cleanup(hid_t fh)
 {
 	int n_ids, i;
@@ -1127,6 +1233,9 @@ static int get_scalar_value(struct hdfile *f, const char *name, void *val,
 		return 1;
 	}
 
+	H5Tclose(type);
+	H5Dclose(dh);
+
 	return 0;
 }
 
@@ -1154,7 +1263,7 @@ static int get_ev_based_value(struct hdfile *f, const char *name,
 	char *subst_name = NULL;
 
 	if ( ev->path_length != 0 ) {
-		subst_name = partial_event_substitution(ev, name);
+		subst_name = retrieve_full_path(ev, name);
 	} else {
 		subst_name = strdup(name);
 	}
@@ -1284,8 +1393,10 @@ int hdfile_get_value(struct hdfile *f, const char *name,
 }
 
 
-void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
-                             struct event *ev, struct image *image)
+static void hdfile_fill_in_beam_parameters(struct beam_params *beam,
+                                           struct hdfile *f,
+                                           struct event *ev,
+                                           struct image *image)
 {
 	double eV;
 
@@ -1298,8 +1409,8 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
 
 		int r;
 
-		r = hdfile_get_value(f, beam->photon_energy_from, ev, &eV,
-		                     H5T_NATIVE_DOUBLE);
+		r = hdfile_get_value(f, beam->photon_energy_from,
+		                     ev, &eV, H5T_NATIVE_DOUBLE);
 		if ( r ) {
 			ERROR("Failed to read '%s'\n",
 			      beam->photon_energy_from);
@@ -1311,6 +1422,36 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
 }
 
 
+static void hdfile_fill_in_clen(struct detector *det, struct hdfile *f,
+                                struct event *ev)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		struct panel *p = &det->panels[i];
+
+		if ( p->clen_from != NULL ) {
+
+			double val;
+			int r;
+
+			r = hdfile_get_value(f, p->clen_from, ev, &val,
+			                     H5T_NATIVE_DOUBLE);
+			if ( r ) {
+				ERROR("Failed to read '%s'\n", p->clen_from);
+			} else {
+				p->clen = val * 1.0e-3;
+			}
+
+		}
+
+		adjust_centering_for_rail(p);
+
+	}
+}
+
+
 int hdf5_read(struct hdfile *f, struct image *image, const char *element,
               int satcorr)
 {
@@ -1326,7 +1467,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
 	if ( element == NULL ) {
 		fail = hdfile_set_first_image(f, "/");
 	} else {
-		fail = hdfile_set_image(f, element, NULL);
+		fail = hdfile_set_image(f, element);
 	}
 
 	if ( fail ) {
@@ -1375,7 +1516,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
 
 	if ( image->beam != NULL ) {
 
-		fill_in_beam_parameters(image->beam, f, NULL, image);
+		hdfile_fill_in_beam_parameters(image->beam, f, NULL, image);
 
 		if ( image->lambda > 1000 ) {
 			/* Error message covers a silly value in the beam file
@@ -1634,10 +1775,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 			if ( !exists ) {
 				ERROR("Cannot find data for panel %s\n",
 				      p->name);
+				free(image->dp);
+				free(image->bad);
+				free(image->sat);
 				return 1;
 			}
 
-			fail = hdfile_set_image(f, panel_full_path, p);
+			fail = hdfile_set_image(f, panel_full_path);
 
 			free(panel_full_path);
 
@@ -1654,9 +1798,12 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 				if ( !exists ) {
 					ERROR("Cannot find data for panel %s\n",
 					      p->name);
+					free(image->dp);
+					free(image->bad);
+					free(image->sat);
 					return 1;
 				}
-				fail = hdfile_set_image(f, p->data, p);
+				fail = hdfile_set_image(f, p->data);
 
 			}
 
@@ -1664,6 +1811,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 		if ( fail ) {
 			ERROR("Couldn't select path for panel %s\n",
 			      p->name);
+			free(image->dp);
+			free(image->bad);
+			free(image->sat);
 			return 1;
 		}
 
@@ -1673,6 +1823,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 		f_count = malloc(hsd->num_dims*sizeof(hsize_t));
 		if ( (f_offset == NULL) || (f_count == NULL ) ) {
 			ERROR("Failed to allocate offset or count.\n");
+			free(image->dp);
+			free(image->bad);
+			free(image->sat);
 			return 1;
 		}
 		for ( hsi=0; hsi<hsd->num_dims; hsi++ ) {
@@ -1700,6 +1853,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 		if ( check < 0 ) {
 			ERROR("Error selecting file dataspace for panel %s\n",
 			      p->name);
+			free(image->dp);
+			free(image->bad);
+			free(image->sat);
 			return 1;
 		}
 
@@ -1713,6 +1869,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 			ERROR("Failed to allocate panel %s\n", p->name);
 			free(f_offset);
 			free(f_count);
+			free(image->dp);
+			free(image->bad);
+			free(image->sat);
 			return 1;
 		}
 		for ( i=0; i<p->w*p->h; i++ ) image->sat[pi][i] = INFINITY;
@@ -1724,6 +1883,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 			      p->name);
 			free(f_offset);
 			free(f_count);
+			for ( i=0; i<=pi; i++ ) {
+				free(image->dp[i]);
+				free(image->sat[i]);
+			}
+			free(image->dp);
+			free(image->bad);
+			free(image->sat);
 			return 1;
 		}
 
@@ -1740,7 +1906,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 			if ( load_satmap(f, ev, p, f_offset, f_count, hsd,
 			                 image->sat[pi]) )
 			{
-				ERROR("Failed to laod sat map for panel %s\n",
+				ERROR("Failed to load sat map for panel %s\n",
 				      p->name);
 			}
 		}
@@ -1753,13 +1919,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 
 	H5Dclose(f->dh);
 	f->data_open = 0;
-	fill_in_values(image->det, f, ev);
+	hdfile_fill_in_clen(image->det, f, ev);
 
 	if ( satcorr ) debodge_saturation(f, image);
 
 	if ( image->beam != NULL ) {
 
-		fill_in_beam_parameters(image->beam, f, ev, image);
+		hdfile_fill_in_beam_parameters(image->beam, f, ev, image);
 
 		if ( (image->lambda > 1.0) || (image->lambda < 1e-20) ) {
 
@@ -1966,7 +2132,7 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name,
 	char *tmp = NULL, *subst_name = NULL;
 
 	if (ev != NULL && ev->path_length != 0 ) {
-		subst_name = partial_event_substitution(ev, name);
+		subst_name = retrieve_full_path(ev, name);
 	} else {
 		subst_name = strdup(name);
 	}
@@ -2166,7 +2332,7 @@ int hdfile_set_first_image(struct hdfile *f, const char *group)
 	for ( i=0; i<n; i++ ) {
 
 		if ( is_image[i] ) {
-			hdfile_set_image(f, names[i], NULL);
+			hdfile_set_image(f, names[i]);
 			for ( j=0; j<n; j++ ) free(names[j]);
 			free(is_image);
 			free(is_group);
@@ -2484,6 +2650,9 @@ static int check_dims(struct hdfile *hdfile, struct panel *p, struct event *ev,
 		return 1;
 	}
 
+	H5Sclose(sh);
+	H5Dclose(dh);
+
 	return 0;
 }
 
diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h
index 8c89eb93..ab53dd2e 100644
--- a/libcrystfel/src/hdf5-file.h
+++ b/libcrystfel/src/hdf5-file.h
@@ -3,11 +3,11 @@
  *
  * Read/write HDF5 data files
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2015 Thomas White <taw@physics.org>
+ *   2009-2017 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
 
  *
@@ -67,8 +67,8 @@ extern int hdf5_read2(struct hdfile *f, struct image *image,
 extern int check_path_existence(hid_t fh, const char *path);
 
 extern struct hdfile *hdfile_open(const char *filename);
-int hdfile_set_image(struct hdfile *f, const char *path,
-                     struct panel *p);
+int hdfile_set_image(struct hdfile *f, const char *path);
+
 extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile,
                                          int binning, int16_t *maxp);
 extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
@@ -78,9 +78,16 @@ extern void hdfile_close(struct hdfile *f);
 
 extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
 
+extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+                       int half_pixel_shift);
+
 extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
                          struct filename_plus_event *fpe);
 
+extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+                           struct filename_plus_event *fpe,
+                           int half_pixel_shift);
+
 extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
 extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
 
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 7bb4cede..aad5017c 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -3,12 +3,12 @@
  *
  * Handle images and image features
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2014      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
- *   2011-2016 Thomas White <taw@physics.org>
+ *   2011-2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -27,14 +27,28 @@
  *
  */
 
-#define _ISOC99_SOURCE
+#include <config.h>
+
 #include <stdlib.h>
 #include <assert.h>
 #include <math.h>
 #include <stdio.h>
+#include <hdf5.h>
+
+#ifdef HAVE_CBFLIB
+#ifdef HAVE_CBFLIB_CBF_H
+#include <cbflib/cbf.h>
+#endif
+#ifdef HAVE_CBF_CBF_H
+#include <cbf/cbf.h>
+#endif
+#endif
 
 #include "image.h"
 #include "utils.h"
+#include "events.h"
+#include "hdf5-file.h"
+#include "detector.h"
 
 /**
  * SECTION:image
@@ -51,6 +65,14 @@
  */
 
 
+struct imagefile
+{
+	enum imagefile_type type;
+	char *filename;
+	struct hdfile *hdfile;
+};
+
+
 struct _imagefeaturelist
 {
 	struct imagefeature	*features;
@@ -300,3 +322,673 @@ void free_all_crystals(struct image *image)
 	free(image->crystals);
 	image->n_crystals = 0;
 }
+
+
+/**************************** Image field lists *******************************/
+
+struct imagefile_field_list
+{
+	char **fields;
+	int n_fields;
+	int max_fields;
+};
+
+
+struct imagefile_field_list *new_imagefile_field_list()
+{
+	struct imagefile_field_list *n;
+
+	n = calloc(1, sizeof(struct imagefile_field_list));
+	if ( n == NULL ) return NULL;
+
+	n->max_fields = 32;
+	n->fields = malloc(n->max_fields*sizeof(char *));
+	if ( n->fields == NULL ) {
+		free(n);
+		return NULL;
+	}
+
+	return n;
+}
+
+
+void free_imagefile_field_list(struct imagefile_field_list *n)
+{
+	int i;
+	for ( i=0; i<n->n_fields; i++ ) {
+		free(n->fields[i]);
+	}
+	free(n->fields);
+	free(n);
+}
+
+
+void add_imagefile_field(struct imagefile_field_list *copyme, const char *name)
+{
+	int i;
+
+	/* Already on the list?   Don't re-add if so. */
+	for ( i=0; i<copyme->n_fields; i++ ) {
+		if ( strcmp(copyme->fields[i], name) == 0 ) return;
+	}
+
+	/* Need more space? */
+	if ( copyme->n_fields == copyme->max_fields ) {
+
+		char **nfields;
+		int nmax = copyme->max_fields + 32;
+
+		nfields = realloc(copyme->fields, nmax*sizeof(char *));
+		if ( nfields == NULL ) {
+			ERROR("Failed to allocate space for new HDF5 field.\n");
+			return;
+		}
+
+		copyme->max_fields = nmax;
+		copyme->fields = nfields;
+
+	}
+
+	copyme->fields[copyme->n_fields] = strdup(name);
+	if ( copyme->fields[copyme->n_fields] == NULL ) {
+		ERROR("Failed to add field for copying '%s'\n", name);
+		return;
+	}
+
+	copyme->n_fields++;
+}
+
+
+/******************************* CBF files ************************************/
+
+#ifdef HAVE_CBFLIB
+
+static char *cbf_strerr(int e)
+{
+	char *err;
+
+	err = malloc(1024);
+	if ( err == NULL ) return NULL;
+
+	err[0] = '\0';
+
+	/* NB Sum of lengths of all strings must be less than 1024 */
+	if ( e & CBF_FORMAT ) strcat(err, "Invalid format");
+	if ( e & CBF_ALLOC ) strcat(err, "Memory allocation failed");
+	if ( e & CBF_ARGUMENT ) strcat(err, "Invalid argument");
+	if ( e & CBF_ASCII ) strcat(err, "Value is ASCII");
+	if ( e & CBF_BINARY ) strcat(err, "Value is binary");
+	if ( e & CBF_BITCOUNT ) strcat(err, "Wrong number of bits");
+	if ( e & CBF_ENDOFDATA ) strcat(err, "End of data");
+	if ( e & CBF_FILECLOSE ) strcat(err, "File close error");
+	if ( e & CBF_FILEOPEN ) strcat(err, "File open error");
+	if ( e & CBF_FILEREAD ) strcat(err, "File read error");
+	if ( e & CBF_FILETELL ) strcat(err, "File tell error");
+	if ( e & CBF_FILEWRITE ) strcat(err, "File write error");
+	if ( e & CBF_IDENTICAL ) strcat(err, "Name already exists");
+	if ( e & CBF_NOTFOUND ) strcat(err, "Not found");
+	if ( e & CBF_OVERFLOW ) strcat(err, "Overflow");
+	if ( e & CBF_UNDEFINED ) strcat(err, "Number undefined");
+	if ( e & CBF_NOTIMPLEMENTED ) strcat(err, "Not implemented");
+
+	return err;
+}
+
+
+static int unpack_panels(struct image *image, signed int *data, int data_width,
+                         int data_height)
+{
+	int pi;
+
+	/* FIXME: Load these masks from an HDF5 file, if filenames are
+	 * given in the geometry file */
+	uint16_t *flags = NULL;
+	float *sat = NULL;
+
+	image->dp = malloc(image->det->n_panels * sizeof(float *));
+	image->bad = malloc(image->det->n_panels * sizeof(int *));
+	image->sat = malloc(image->det->n_panels * sizeof(float *));
+	if ( (image->dp == NULL) || (image->bad == NULL)
+	  || (image->sat == NULL) )
+	{
+		ERROR("Failed to allocate panels.\n");
+		return 1;
+	}
+
+	for ( pi=0; pi<image->det->n_panels; pi++ ) {
+
+		struct panel *p;
+		int fs, ss;
+
+		p = &image->det->panels[pi];
+		image->dp[pi] = malloc(p->w*p->h*sizeof(float));
+		image->bad[pi] = calloc(p->w*p->h, sizeof(int));
+		image->sat[pi] = malloc(p->w*p->h*sizeof(float));
+		if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL)
+		  || (image->sat[pi] == NULL) )
+		{
+			ERROR("Failed to allocate panel\n");
+			return 1;
+		}
+
+		if ( p->mask != NULL ) {
+			ERROR("WARNING: Bad pixel masks do not currently work "
+			      "with CBF files\n");
+			ERROR(" (bad pixel regions specified in the geometry "
+			      "file will be used, however)\n");
+		}
+
+		if ( p->satmap != NULL ) {
+			ERROR("WARNING: Saturation maps do not currently work "
+			      "with CBF files\n");
+		}
+
+		if ( (p->orig_min_fs + p->w > data_width)
+		  || (p->orig_min_ss + p->h > data_height) )
+		{
+			ERROR("Panel %s is outside range of data in CBF file\n",
+			      p->name);
+			return 1;
+		}
+
+		for ( ss=0; ss<p->h; ss++ ) {
+		for ( fs=0; fs<p->w; fs++ ) {
+
+			int idx;
+			int cfs, css;
+			int bad = 0;
+
+			cfs = fs+p->orig_min_fs;
+			css = ss+p->orig_min_ss;
+			idx = cfs + css*data_width;
+
+			image->dp[pi][fs+p->w*ss] = data[idx];
+
+			if ( sat != NULL ) {
+				image->sat[pi][fs+p->w*ss] = sat[idx];
+			} else {
+				image->sat[pi][fs+p->w*ss] = INFINITY;
+			}
+
+			if ( p->no_index ) bad = 1;
+
+			if ( in_bad_region(image->det, p, cfs, css) ) {
+				bad = 1;
+			}
+
+			if ( flags != NULL ) {
+
+				int f;
+
+				f = flags[idx];
+
+				/* Bad if it's missing any of the "good" bits */
+				if ( (f & image->det->mask_good)
+				       != image->det->mask_good ) bad = 1;
+
+				/* Bad if it has any of the "bad" bits. */
+				if ( f & image->det->mask_bad ) bad = 1;
+
+			}
+			image->bad[pi][fs+p->w*ss] = bad;
+
+		}
+		}
+
+	}
+
+	return 0;
+}
+
+
+static void cbf_fill_in_beam_parameters(struct beam_params *beam,
+                                        struct imagefile *f,
+                                        struct image *image)
+{
+	double eV;
+
+	if ( beam->photon_energy_from == NULL ) {
+
+		/* Explicit value given */
+		eV = beam->photon_energy;
+
+	} else {
+
+		ERROR("Can't get photon energy from CBF yet.\n");
+		eV = 0.0;
+
+	}
+
+	image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale;
+}
+
+
+static void cbf_fill_in_clen(struct detector *det, struct imagefile *f)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		struct panel *p = &det->panels[i];
+
+		if ( p->clen_from != NULL ) {
+
+			ERROR("Can't get clen from CBF yet.\n");
+
+		}
+
+		adjust_centering_for_rail(p);
+
+	}
+}
+
+
+static int read_cbf(struct imagefile *f, struct image *image)
+{
+	cbf_handle cbfh;
+	FILE *fh;
+	int r;
+	unsigned int compression;
+	int binary_id, minelement, maxelement, elsigned, elunsigned;
+	size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding;
+	const char *byteorder;
+	signed int *data;
+
+	if ( image->det == NULL ) {
+		ERROR("read_cbf() needs a geometry\n");
+		return 1;
+	}
+
+	if ( cbf_make_handle(&cbfh) ) {
+		ERROR("Failed to allocate CBF handle\n");
+		return 1;
+	}
+
+	fh = fopen(f->filename, "rb");
+	if ( fh == NULL ) {
+		ERROR("Failed to open '%s'\n", f->filename);
+		return 1;
+	}
+	/* CBFlib calls fclose(fh) when it's ready */
+
+	if ( cbf_read_widefile(cbfh, fh, 0) ) {
+		ERROR("Failed to read CBF file '%s'\n", f->filename);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	/* Select row 0 in data column inside array_data */
+	cbf_find_category(cbfh, "array_data");
+	cbf_find_column(cbfh, "data");
+	cbf_select_row(cbfh, 0);
+
+	/* Get parameters for array read */
+	r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id,
+	                                         &elsize, &elsigned, &elunsigned,
+	                                         &elements,
+	                                         &minelement, &maxelement,
+	                                         &byteorder,
+	                                         &dimfast, &dimmid, &dimslow,
+	                                         &padding);
+	if ( r ) {
+		char *err = cbf_strerr(r);
+		ERROR("Failed to read CBF array parameters: %s\n", err);
+		free(err);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( dimslow != 0 ) {
+		ERROR("CBF data array is 3D - don't know what to do with it\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( dimfast*dimmid*elsize > 10e9 ) {
+		ERROR("CBF data is far too big (%i x %i x %i bytes).\n",
+		      (int)dimfast, (int)dimmid, (int)elsize);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( elsize != 4 ) {
+		STATUS("Don't know what to do with element size %i\n",
+		       (int)elsize);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( strcmp(byteorder, "little_endian") != 0 ) {
+		STATUS("Don't know what to do with non-little-endian datan\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	data = malloc(elsize*dimfast*dimmid);
+	if ( data == NULL ) {
+		ERROR("Failed to allocate memory for CBF data\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1,
+	                         elsize*dimfast*dimmid, &elread);
+	if ( r ) {
+		char *err = cbf_strerr(r);
+		ERROR("Failed to read CBF array: %s\n", err);
+		free(err);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	unpack_panels(image, data, dimfast, dimmid);
+	free(data);
+
+	if ( image->beam != NULL ) {
+		cbf_fill_in_beam_parameters(image->beam, f, image);
+		if ( image->lambda > 1000 ) {
+			ERROR("WARNING: Missing or nonsensical wavelength "
+			      "(%e m) for %s.\n",
+			      image->lambda, image->filename);
+		}
+	}
+	cbf_fill_in_clen(image->det, f);
+	fill_in_adu(image);
+
+	cbf_free_handle(cbfh);
+	return 0;
+}
+
+
+static float *convert_float(signed int *data, int w, int h)
+{
+	float *df;
+	long int i;
+
+	df = malloc(sizeof(float)*w*h);
+	if ( df == NULL ) return NULL;
+
+	for ( i=0; i<w*h; i++ ) {
+		df[i] = data[i];
+	}
+
+	return df;
+}
+
+
+static int read_cbf_simple(struct imagefile *f, struct image *image)
+{
+	cbf_handle cbfh;
+	FILE *fh;
+	int r;
+	unsigned int compression;
+	int binary_id, minelement, maxelement, elsigned, elunsigned;
+	size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding;
+	const char *byteorder;
+	signed int *data;
+
+	if ( cbf_make_handle(&cbfh) ) {
+		ERROR("Failed to allocate CBF handle\n");
+		return 1;
+	}
+
+	fh = fopen(f->filename, "rb");
+	if ( fh == NULL ) {
+		ERROR("Failed to open '%s'\n", f->filename);
+		return 1;
+	}
+	/* CBFlib calls fclose(fh) when it's ready */
+
+	if ( cbf_read_widefile(cbfh, fh, 0) ) {
+		ERROR("Failed to read CBF file '%s'\n", f->filename);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	/* Select row 0 in data column inside array_data */
+	cbf_find_category(cbfh, "array_data");
+	cbf_find_column(cbfh, "data");
+	cbf_select_row(cbfh, 0);
+
+	/* Get parameters for array read */
+	r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id,
+	                                         &elsize, &elsigned, &elunsigned,
+	                                         &elements,
+	                                         &minelement, &maxelement,
+	                                         &byteorder,
+	                                         &dimfast, &dimmid, &dimslow,
+	                                         &padding);
+	if ( r ) {
+		char *err = cbf_strerr(r);
+		ERROR("Failed to read CBF array parameters: %s\n", err);
+		free(err);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( dimslow != 0 ) {
+		ERROR("CBF data array is 3D - don't know what to do with it\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( dimfast*dimmid*elsize > 10e9 ) {
+		ERROR("CBF data is far too big (%i x %i x %i bytes).\n",
+		      (int)dimfast, (int)dimmid, (int)elsize);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( elsize != 4 ) {
+		STATUS("Don't know what to do with element size %i\n",
+		       (int)elsize);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	if ( strcmp(byteorder, "little_endian") != 0 ) {
+		STATUS("Don't know what to do with non-little-endian datan\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	data = malloc(elsize*dimfast*dimmid);
+	if ( data == NULL ) {
+		ERROR("Failed to allocate memory for CBF data\n");
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1,
+	                         elsize*dimfast*dimmid, &elread);
+	if ( r ) {
+		char *err = cbf_strerr(r);
+		ERROR("Failed to read CBF array: %s\n", err);
+		free(err);
+		cbf_free_handle(cbfh);
+		return 1;
+	}
+
+	image->det = simple_geometry(image, dimfast, dimmid);
+	image->dp = malloc(sizeof(float *));
+	if ( image->dp == NULL ) {
+		ERROR("Failed to allocate dp array\n");
+		return 1;
+	}
+	image->dp[0] = convert_float(data, dimfast, dimmid);
+	if ( image->dp[0] == NULL ) {
+		ERROR("Failed to allocate dp array\n");
+		return 1;
+	}
+
+	if ( image->beam != NULL ) {
+		cbf_fill_in_beam_parameters(image->beam, f, image);
+		if ( image->lambda > 1000 ) {
+			ERROR("WARNING: Missing or nonsensical wavelength "
+			      "(%e m) for %s.\n",
+			      image->lambda, image->filename);
+		}
+	}
+	cbf_fill_in_clen(image->det, f);
+	fill_in_adu(image);
+
+	cbf_free_handle(cbfh);
+	return 0;
+}
+
+#else /* HAVE_CBFLIB */
+
+static int read_cbf_simple(struct imagefile *f, struct image *image)
+{
+	ERROR("This version of CrystFEL was compiled without CBF support.\n");
+	return 1;
+}
+
+
+static int read_cbf(struct imagefile *f, struct image *image)
+{
+	ERROR("This version of CrystFEL was compiled without CBF support.\n");
+	return 1;
+}
+
+
+#endif /* HAVE_CBFLIB */
+
+
+/****************************** Image files ***********************************/
+
+
+static signed int is_cbf_file(const char *filename)
+{
+	FILE *fh;
+	char line[1024];
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) return -1;
+
+	if ( fgets(line, 1024, fh) == NULL ) return -1;
+	fclose(fh);
+
+	if ( strstr(line, "CBF") == NULL ) {
+		return 0;
+	}
+
+	return 1;
+}
+
+
+struct imagefile *imagefile_open(const char *filename)
+{
+	struct imagefile *f;
+
+	f = malloc(sizeof(struct imagefile));
+	if ( f == NULL ) return NULL;
+
+	if ( H5Fis_hdf5(filename) > 0 ) {
+
+		/* This is an HDF5, pass through to HDF5 layer */
+		f->type = IMAGEFILE_HDF5;
+		f->hdfile = hdfile_open(filename);
+
+		if ( f->hdfile == NULL ) {
+			free(f);
+			return NULL;
+		}
+
+	} else if ( is_cbf_file(filename) > 0 ) {
+
+		f->type = IMAGEFILE_CBF;
+
+	}
+
+	f->filename = strdup(filename);
+	return f;
+}
+
+
+int imagefile_read(struct imagefile *f, struct image *image,
+                   struct event *event)
+{
+	if ( f->type == IMAGEFILE_HDF5 ) {
+		return hdf5_read2(f->hdfile, image, event, 0);
+	} else if ( f->type == IMAGEFILE_CBF ) {
+		return read_cbf(f, image);
+	} else {
+		ERROR("Unknown file type %i\n", f->type);
+		return 1;
+	}
+}
+
+
+/* Read a simple file, no multi-event, no prior geometry etc, and
+ * generate a geometry for it */
+int imagefile_read_simple(struct imagefile *f, struct image *image)
+{
+	if ( f->type == IMAGEFILE_HDF5 ) {
+		return hdf5_read(f->hdfile, image, NULL, 0);
+	} else {
+		return read_cbf_simple(f, image);
+	}
+}
+
+
+enum imagefile_type imagefile_get_type(struct imagefile *f)
+{
+	assert(f != NULL);
+	return f->type;
+}
+
+
+struct hdfile *imagefile_get_hdfile(struct imagefile *f)
+{
+	if ( f->type != IMAGEFILE_HDF5 ) {
+		ERROR("Not an HDF5 file!\n");
+		return NULL;
+	}
+
+	return f->hdfile;
+}
+
+
+void imagefile_copy_fields(struct imagefile *f,
+                           const struct imagefile_field_list *copyme,
+                           FILE *fh, struct event *ev)
+{
+	int i;
+
+	if ( copyme == NULL ) return;
+
+	for ( i=0; i<copyme->n_fields; i++ ) {
+
+		char *val;
+		char *field;
+
+		field = copyme->fields[i];
+
+		if ( f->type == IMAGEFILE_HDF5 ) {
+			val = hdfile_get_string_value(f->hdfile, field, ev);
+			if ( field[0] == '/' ) {
+				fprintf(fh, "hdf5%s = %s\n", field, val);
+			} else {
+				fprintf(fh, "hdf5/%s = %s\n", field, val);
+			}
+			free(val);
+
+		} else {
+			STATUS("Mock CBF variable\n");
+			fprintf(fh, "cbf/%s = %s\n", field, "(FIXME)");
+		}
+
+
+	}
+}
+
+
+void imagefile_close(struct imagefile *f)
+{
+	if ( f->type == IMAGEFILE_HDF5 ) {
+		hdfile_close(f->hdfile);
+	}
+	free(f->filename);
+	free(f);
+}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 9fd9b495..9719bb59 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -3,11 +3,11 @@
  *
  * Handle images and image features
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2016 Thomas White <taw@physics.org>
+ *   2009-2017 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
  *
  *
@@ -44,6 +44,8 @@ struct detector;
 struct imagefeature;
 struct sample;
 struct image;
+struct imagefile;
+struct imagefile_field_list;
 
 #include "utils.h"
 #include "cell.h"
@@ -51,6 +53,7 @@ struct image;
 #include "reflist.h"
 #include "crystal.h"
 #include "index.h"
+#include "events.h"
 
 /**
  * SpectrumType:
@@ -67,7 +70,26 @@ typedef enum {
 	SPECTRUM_TWOCOLOUR
 } SpectrumType;
 
-/* Structure describing a feature in an image */
+
+/**
+ * imagefeature:
+ *  @parent: Image this feature belongs to
+ *  @fs: Fast scan coordinate
+ *  @ss: Slow scan coordinate
+ *  @p: Pointer to panel
+ *  @intensity: Intensity of peak
+ *  @rx: Reciprocal x coordinate in m^-1
+ *  @ry: Reciprocal y coordinate in m^-1
+ *  @rz: Reciprocal z coordinate in m^-1
+ *  @name: Text name for feature
+ *
+ *  Represents a peak in an image.
+ *
+ *  Note carefully that the @fs and @ss coordinates are the distances, measured
+ *  in pixels, from the corner of the panel.  They are NOT pixel indices.
+ *  If the peak is in the middle of the first pixel, its coordinates would be
+ *  0.5,0.5.
+ */
 struct imagefeature {
 
 	struct image                    *parent;
@@ -81,12 +103,21 @@ struct imagefeature {
 	double                          ry;
 	double                          rz;
 
-	/* Internal use only */
+	const char                      *name;
+
+	/*< private >*/
 	int                             valid;
+};
 
-	const char                      *name;
+
+/* An enum representing the image file formats we can handle */
+enum imagefile_type
+{
+	IMAGEFILE_HDF5,
+	IMAGEFILE_CBF
 };
 
+
 /* An opaque type representing a list of image features */
 typedef struct _imagefeaturelist ImageFeatureList;
 
@@ -98,44 +129,46 @@ struct sample
 };
 
 
+/**
+ * beam_params:
+ *  @photon_energy: eV per photon
+ *  @photon_energy_from: HDF5 dataset name
+ *  @photon_energy_scale: Scale factor for photon energy, if it comes from HDF5
+ */
 struct beam_params
 {
-	double photon_energy;  /* eV per photon */
-	char  *photon_energy_from; /* HDF5 dataset name */
-	double photon_energy_scale;  /* Scale factor for photon energy, if the
-	                              * energy is to be from the HDF5 file */
+	double photon_energy;
+	char  *photon_energy_from;
+	double photon_energy_scale;
 };
 
 
 /**
  * image:
- *
- * <programlisting>
- * struct image
- * {
- *    Crystal                 **crystals;
- *    int                     n_crystals;
- *    IndexingMethod          indexed_by;
- *
- *    struct detector         *det;
- *    struct beam_params      *beam;
- *    char                    *filename;
- *    const struct copy_hdf5_field *copyme;
- *
- *    int                     id;
- *
- *    double                  lambda;
- *    double                  div;
- *    double                  bw;
- *
- *    int                     width;
- *    int                     height;
- *
- *    long long int           num_peaks;
- *    long long int           num_saturated_peaks;
- *    ImageFeatureList        *features;
- * };
- * </programlisting>
+ *   @crystals: Array of crystals in the image
+ *   @n_crystals: The number of crystals in the image
+ *   @indexed_by: Indexing method which indexed this pattern
+ *   @det: Detector structure
+ *   @beam: Beam parameters structure
+ *   @filename: Filename for the image file
+ *   @copyme: Fields to copy from the image file to the stream
+ *   @id: ID number of the thread handling this image
+ *   @serial: Serial number for this image
+ *   @lambda: Wavelength
+ *   @div: Divergence
+ *   @bw: Bandwidth
+ *   @num_peaks: The number of peaks
+ *   @num_saturated_peaks: The number of saturated peaks
+ *   @features: The peaks found in the image
+ *   @dp: The image data, by panel
+ *   @bad: The bad pixel mask, array by panel
+ *   @sat: The per-pixel saturation mask, array by panel
+ *   @event: Event ID for the image
+ *   @stuff_from_stream: Items read back from the stream
+ *   @avg_clen: Mean of camera length values for all panels
+ *   @spectrum: Spectrum information
+ *   @nsamples: Number of spectrum samples
+ *   @spectrum_size: SIze of spectrum array
  *
  * The field <structfield>data</structfield> contains the raw image data, if it
  * is currently available.  The data might be available throughout the
@@ -152,8 +185,8 @@ struct beam_params
  * returned by the low-level indexing system. <structfield>n_crystals</structfield>
  * is the number of crystals which were found in the image.
  *
- * <structfield>copyme</structfield> represents a list of HDF5 fields to copy
- * to the output stream.
+ * <structfield>copyme</structfield> represents a list of fields in the image
+ * file (e.g. HDF5 fields or CBF headers) to copy to the output stream.
  **/
 struct image;
 
@@ -171,7 +204,7 @@ struct image {
 	struct beam_params      *beam;  /* The nominal beam parameters */
 	char                    *filename;
 	struct event            *event;
-	const struct copy_hdf5_field *copyme;
+	const struct imagefile_field_list *copyme;
 	struct stuff_from_stream *stuff_from_stream;
 
 	double                  avg_clen;  /* Average camera length extracted
@@ -192,8 +225,8 @@ struct image {
 	double                  bw;            /* Bandwidth as a fraction */
 
 	/* Detected peaks */
-	long long int           num_peaks;
-	long long int           num_saturated_peaks;
+	long long               num_peaks;
+	long long               num_saturated_peaks;
 	ImageFeatureList        *features;
 
 };
@@ -233,6 +266,25 @@ 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);
 
+/* Image files */
+extern struct imagefile *imagefile_open(const char *filename);
+extern int imagefile_read(struct imagefile *f, struct image *image,
+                          struct event *event);
+extern int imagefile_read_simple(struct imagefile *f, struct image *image);
+extern struct hdfile *imagefile_get_hdfile(struct imagefile *f);
+extern enum imagefile_type imagefile_get_type(struct imagefile *f);
+extern void imagefile_copy_fields(struct imagefile *f,
+                                  const struct imagefile_field_list *copyme,
+                                  FILE *fh, struct event *ev);
+extern void imagefile_close(struct imagefile *f);
+
+/* Field lists */
+extern struct imagefile_field_list *new_imagefile_field_list(void);
+extern void free_imagefile_field_list(struct imagefile_field_list *f);
+
+extern void add_imagefile_field(struct imagefile_field_list *copyme,
+                                const char *name);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index e7fd879d..bdabb062 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -59,6 +59,17 @@
 #include "taketwo.h"
 
 
+struct _indexingprivate
+{
+	UnitCell *target_cell;
+	float tolerance[4];
+
+	int n_methods;
+	IndexingMethod *methods;
+	void **engine_private;
+};
+
+
 static int debug_index(struct image *image)
 {
 	Crystal *cr = crystal_new();
@@ -72,78 +83,119 @@ static int debug_index(struct image *image)
 }
 
 
-IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
-                                   struct detector *det, float *ltl,
-                                   const char *options)
+static void *prepare_method(IndexingMethod *m, UnitCell *cell,
+                            struct detector *det, float *ltl,
+                            const char *options)
 {
-	int n;
-	int nm = 0;
-	IndexingPrivate **iprivs;
+	char *str;
+	IndexingMethod in = *m;
+	void *priv = NULL;
 
-	while ( indm[nm] != INDEXING_NONE ) nm++;
-	iprivs = malloc((nm+1) * sizeof(IndexingPrivate *));
+	switch ( *m & INDEXING_METHOD_MASK ) {
 
-	for ( n=0; n<nm; n++ ) {
+		case INDEXING_NONE :
+		priv = "none";
+		break;
 
-		int i;
-		IndexingMethod in;
-		char *str;
+		case INDEXING_DIRAX :
+		priv = dirax_prepare(m, cell, det, ltl);
+		break;
 
-		in = indm[n];
+		case INDEXING_ASDF :
+		priv = asdf_prepare(m, cell, det, ltl);
+		break;
 
-		switch ( indm[n] & INDEXING_METHOD_MASK ) {
+		case INDEXING_MOSFLM :
+		priv = mosflm_prepare(m, cell, det, ltl);
+		break;
 
-			case INDEXING_DIRAX :
-			iprivs[n] = dirax_prepare(&indm[n], cell, det, ltl);
-			break;
+		case INDEXING_XDS :
+		priv = xds_prepare(m, cell, det, ltl);
+		break;
 
-			case INDEXING_ASDF :
-			iprivs[n] = asdf_prepare(&indm[n], cell, det, ltl);
-			break;
+		case INDEXING_DEBUG :
+		priv = (IndexingPrivate *)strdup("Hello!");
+		break;
 
-			case INDEXING_MOSFLM :
-			iprivs[n] = mosflm_prepare(&indm[n], cell, det, ltl);
-			break;
+		case INDEXING_FELIX :
+		priv = felix_prepare(m, cell, det, ltl, options);
+		break;
 
-			case INDEXING_XDS :
-			iprivs[n] = xds_prepare(&indm[n], cell, det, ltl);
-			break;
+		case INDEXING_TAKETWO :
+		priv = taketwo_prepare(m, cell, det, ltl);
+		break;
 
-			case INDEXING_DEBUG :
-			iprivs[n] = (IndexingPrivate *)strdup("Hello!");
-			break;
+		default :
+		ERROR("Don't know how to prepare indexing method %i\n", *m);
+		break;
 
-			case INDEXING_FELIX :
-			iprivs[n] = felix_prepare(&indm[n], cell, det, ltl,
-			                          options);
-			break;
+	}
 
-			case INDEXING_TAKETWO :
-			iprivs[n] = taketwo_prepare(&indm[n], cell, det, ltl);
-			break;
+	str = indexer_str(*m);
 
-			default :
-			ERROR("Don't know how to prepare indexing method %i\n",
-			      indm[n]);
-			break;
+	if ( priv == NULL ) {
+		ERROR("Failed to prepare indexing method %s\n", str);
+		free(str);
+		return NULL;
+	}
 
-		}
+	if ( *m != INDEXING_NONE ) {
+		STATUS("Prepared indexing method %s\n", str);
+	}
+	free(str);
 
-		if ( iprivs[n] == NULL ) return NULL;
+	if ( in != *m ) {
+		ERROR("Note: flags were altered to take into account "
+		      "the information provided and/or the limitations "
+		      "of the indexing method.\nPlease check the "
+		      "methods listed above carefully.\n");
+	}
 
-		str = indexer_str(indm[n]);
-		STATUS("Prepared indexing method %i: %s\n", n, str);
-		free(str);
+	return priv;
+}
+
+
+IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
+                                struct detector *det, float *ltl,
+                                int no_refine, const char *options)
+{
+	int i, n;
+	char **method_strings;
+	IndexingPrivate *ipriv;
+
+	/* Parse indexing methods */
+	n = assplode(method_list, ",", &method_strings, ASSPLODE_NONE);
+
+	IndexingMethod *methods = malloc(n * sizeof(IndexingMethod));
+	if ( methods == NULL ) {
+		ERROR("Failed to allocate indexing method list\n");
+		return NULL;
+	}
 
-		if ( in != indm[n] ) {
-			ERROR("Note: flags were altered to take into account "
-			      "the information provided and/or the limitations "
-			      "of the indexing method.\nPlease check the "
-			      "methods listed above carefully.\n");
+	for ( i=0; i<n; i++ ) {
+		methods[i] = get_indm_from_string(method_strings[i]);
+		if ( methods[i] == INDEXING_ERROR ) {
+			free(methods);
+			return NULL;
 		}
+	}
+
+	ipriv = malloc(sizeof(struct _indexingprivate));
+	if ( ipriv == NULL ) {
+		ERROR("Failed to allocate indexing data\n");
+		return NULL;
+	}
+
+	ipriv->engine_private = malloc((n+1) * sizeof(void *));
 
-		for ( i=0; i<n; i++ ) {
-			if ( indm[i] == indm[n] ) {
+	for ( i=0; i<n; i++ ) {
+
+		int j;
+
+		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
+		                                          det, ltl, options);
+		for ( j=0; j<i; j++ ) {
+			if ( methods[i] == methods[j] ) {
 				ERROR("Duplicate indexing method.\n");
 				ERROR("Have you specified some flags which "
 				      "aren't accepted by one of your "
@@ -152,59 +204,66 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
 			}
 		}
 
+	}
+
+	ipriv->methods = methods;
+	ipriv->n_methods = n;
 
+	if ( cell != NULL ) {
+		ipriv->target_cell = cell_new_from_cell(cell);
+	} else {
+		ipriv->target_cell = NULL;
 	}
-	iprivs[n] = NULL;
+	for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i];
 
-	return iprivs;
+	return ipriv;
 }
 
 
-void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
+void cleanup_indexing(IndexingPrivate *ipriv)
 {
-	int n = 0;
+	int n;
 
-	if ( indms == NULL ) return;  /* Nothing to do */
-	if ( privs == NULL ) return;  /* Nothing to do */
+	if ( ipriv == NULL ) return;  /* Nothing to do */
 
-	while ( indms[n] != INDEXING_NONE ) {
+	for ( n=0; n<ipriv->n_methods; n++ ) {
 
-		switch ( indms[n] & INDEXING_METHOD_MASK ) {
+		switch ( ipriv->methods[n] & INDEXING_METHOD_MASK ) {
 
 			case INDEXING_NONE :
 			break;
 
 			case INDEXING_DIRAX :
-			dirax_cleanup(privs[n]);
+			dirax_cleanup(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_ASDF :
-			asdf_cleanup(privs[n]);
+			asdf_cleanup(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_MOSFLM :
-			mosflm_cleanup(privs[n]);
+			mosflm_cleanup(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_XDS :
-			xds_cleanup(privs[n]);
+			xds_cleanup(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_FELIX :
-			felix_cleanup(privs[n]);
+			felix_cleanup(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_DEBUG :
-			free(privs[n]);
+			free(ipriv->engine_private[n]);
 			break;
 
 			case INDEXING_TAKETWO :
-			taketwo_cleanup(privs[n]);
+			taketwo_cleanup(ipriv->engine_private[n]);
 			break;
 
 			default :
 			ERROR("Don't know how to clean up indexing method %i\n",
-			      indms[n]);
+			      ipriv->methods[n]);
 			break;
 
 		}
@@ -213,8 +272,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
 
 	}
 
-	free(indms);
-	free(privs);
+	free(ipriv->methods);
+	free(ipriv->engine_private);
+	cell_free(ipriv->target_cell);
+	free(ipriv);
 }
 
 
@@ -243,7 +304,7 @@ void map_all_peaks(struct image *image)
 
 /* Return non-zero for "success" */
 static int try_indexer(struct image *image, IndexingMethod indm,
-                       IndexingPrivate *ipriv)
+                       IndexingPrivate *ipriv, void *mpriv)
 {
 	int i, r;
 	int n_bad = 0;
@@ -254,19 +315,19 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		return 0;
 
 		case INDEXING_DIRAX :
-		r = run_dirax(image, ipriv);
+		r = run_dirax(image, mpriv);
 		break;
 
 		case INDEXING_ASDF :
-		r = run_asdf(image, ipriv);
+		r = run_asdf(image, mpriv);
 		break;
 
 		case INDEXING_MOSFLM :
-		r = run_mosflm(image, ipriv);
+		r = run_mosflm(image, mpriv);
 		break;
 
 		case INDEXING_XDS :
-		r = run_xds(image, ipriv);
+		r = run_xds(image, mpriv);
 		break;
 
 		case INDEXING_DEBUG :
@@ -274,11 +335,11 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		break;
 
 		case INDEXING_FELIX :
-		r = felix_index(image, ipriv);
+		r = felix_index(image, mpriv);
 		break;
 
 		case INDEXING_TAKETWO :
-		r = taketwo_index(image, ipriv);
+		r = taketwo_index(image, mpriv);
 		break;
 
 		default :
@@ -296,11 +357,35 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		crystal_set_mosaicity(cr, 0.0);
 
 		/* Prediction refinement if requested */
-		if ( (indm & INDEXING_REFINE)
-		   && (refine_prediction(image, cr) != 0) )
-		{
-			crystal_set_user_flag(cr, 1);
-			n_bad++;
+		if ( indm & INDEXING_REFINE ) {
+
+			UnitCell *out;
+
+			if ( refine_prediction(image, cr) ) {
+				crystal_set_user_flag(cr, 1);
+				n_bad++;
+				continue;
+			}
+
+			if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS)
+			  || (indm & INDEXING_CHECK_CELL_AXES) )
+			{
+
+				/* Check that the cell parameters are still
+				 * within the tolerance */
+				out = match_cell(crystal_get_cell(cr),
+				                 ipriv->target_cell, 0,
+				                 ipriv->tolerance, 0);
+
+				if ( out == NULL ) {
+					crystal_set_user_flag(cr, 1);
+					n_bad++;
+				}
+
+				cell_free(out);
+
+			}
+
 		}
 
 	}
@@ -434,20 +519,18 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image)
 	}
 }
 
-void index_pattern(struct image *image,
-                   IndexingMethod *indms, IndexingPrivate **iprivs)
+void index_pattern(struct image *image, IndexingPrivate *ipriv)
 {
-	index_pattern_2(image, indms, iprivs, NULL);
+	index_pattern_2(image, ipriv, NULL);
 }
 
 
-void index_pattern_2(struct image *image, IndexingMethod *indms,
-                     IndexingPrivate **iprivs, int *ping)
+void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping)
 {
 	int n = 0;
 	ImageFeatureList *orig;
 
-	if ( indms == NULL ) return;
+	if ( ipriv == NULL ) return;
 
 	map_all_peaks(image);
 	image->crystals = NULL;
@@ -455,7 +538,7 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
 
 	orig = image->features;
 
-	while ( indms[n] != INDEXING_NONE ) {
+	for ( n=0; n<ipriv->n_methods; n++ ) {
 
 		int done = 0;
 		int r;
@@ -466,10 +549,11 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
 
 		do {
 
-			r = try_indexer(image, indms[n], iprivs[n]);
+			r = try_indexer(image, ipriv->methods[n],
+			                ipriv, ipriv->engine_private[n]);
 			success += r;
 			ntry++;
-			done = finished_retry(indms[n], r, image);
+			done = finished_retry(ipriv->methods[n], r, image);
 			if ( ntry > 5 ) done = 1;
 			if ( ping != NULL ) (*ping)++;
 
@@ -481,11 +565,14 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
 		 * crystals with a different indexing method) */
 		if ( success ) break;
 
-		n++;
+	}
 
+	if ( n < ipriv->n_methods ) {
+		image->indexed_by = ipriv->methods[n];
+	} else {
+		image->indexed_by = INDEXING_NONE;
 	}
 
-	image->indexed_by = indms[n];
 	image->features = orig;
 }
 
@@ -656,100 +743,123 @@ char *indexer_str(IndexingMethod indm)
 }
 
 
-IndexingMethod *build_indexer_list(const char *str)
+static IndexingMethod warn_method(const char *str)
+{
+	ERROR("Indexing method must contain exactly one engine name: '%s'\n",
+	      str);
+	return INDEXING_ERROR;
+}
+
+
+IndexingMethod get_indm_from_string(const char *str)
 {
 	int n, i;
-	char **methods;
-	IndexingMethod *list;
-	int nmeth = 0;
+	char **bits;
+	IndexingMethod method = INDEXING_NONE;
+	int have_method = 0;
 
-	n = assplode(str, ",-", &methods, ASSPLODE_NONE);
-	list = malloc((n+1)*sizeof(IndexingMethod));
+	n = assplode(str, "-", &bits, ASSPLODE_NONE);
 
-	nmeth = -1;  /* So that the first method is #0 */
 	for ( i=0; i<n; i++ ) {
 
-		if ( strcmp(methods[i], "dirax") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_DIRAX;
+		if ( strcmp(bits[i], "dirax") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_DIRAX;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "asdf") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_ASDF;
+		} else if ( strcmp(bits[i], "asdf") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_ASDF;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "mosflm") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_MOSFLM;
+		} else if ( strcmp(bits[i], "mosflm") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_MOSFLM;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "xds") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_XDS;
+		} else if ( strcmp(bits[i], "xds") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_XDS;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "felix") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_FELIX;
+		} else if ( strcmp(bits[i], "felix") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_FELIX;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "taketwo") == 0) {
-			list[++nmeth] = INDEXING_DEFAULTS_TAKETWO;
+		} else if ( strcmp(bits[i], "taketwo") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_TAKETWO;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "none") == 0) {
-			list[++nmeth] = INDEXING_NONE;
+		} else if ( strcmp(bits[i], "none") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_NONE;
+			have_method = 1;
 
-		} else if ( strcmp(methods[i], "simulation") == 0) {
-			list[++nmeth] = INDEXING_SIMULATION;
-			return list;
+		} else if ( strcmp(bits[i], "simulation") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_SIMULATION;
+			return method;
 
-		} else if ( strcmp(methods[i], "debug") == 0) {
-			list[++nmeth] = INDEXING_DEBUG;
-			return list;
+		} else if ( strcmp(bits[i], "debug") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEBUG;
+			return method;
 
-		} else if ( strcmp(methods[i], "raw") == 0) {
-			list[nmeth] = set_raw(list[nmeth]);
+		} else if ( strcmp(bits[i], "raw") == 0) {
+			method = set_raw(method);
 
-		} else if ( strcmp(methods[i], "bad") == 0) {
-			list[nmeth] = set_bad(list[nmeth]);
+		} else if ( strcmp(bits[i], "bad") == 0) {
+			method = set_bad(method);
 
-		} else if ( strcmp(methods[i], "comb") == 0) {
-			list[nmeth] = set_comb(list[nmeth]);  /* Default */
+		} else if ( strcmp(bits[i], "comb") == 0) {
+			method = set_comb(method);  /* Default */
 
-		} else if ( strcmp(methods[i], "axes") == 0) {
-			list[nmeth] = set_axes(list[nmeth]);
+		} else if ( strcmp(bits[i], "axes") == 0) {
+			method = set_axes(method);
 
-		} else if ( strcmp(methods[i], "latt") == 0) {
-			list[nmeth] = set_lattice(list[nmeth]);
+		} else if ( strcmp(bits[i], "latt") == 0) {
+			method = set_lattice(method);
 
-		} else if ( strcmp(methods[i], "nolatt") == 0) {
-			list[nmeth] = set_nolattice(list[nmeth]);
+		} else if ( strcmp(bits[i], "nolatt") == 0) {
+			method = set_nolattice(method);
 
-		} else if ( strcmp(methods[i], "cell") == 0) {
-			list[nmeth] = set_cellparams(list[nmeth]);
+		} else if ( strcmp(bits[i], "cell") == 0) {
+			method = set_cellparams(method);
 
-		} else if ( strcmp(methods[i], "nocell") == 0) {
-			list[nmeth] = set_nocellparams(list[nmeth]);
+		} else if ( strcmp(bits[i], "nocell") == 0) {
+			method = set_nocellparams(method);
 
-		} else if ( strcmp(methods[i], "retry") == 0) {
-			list[nmeth] |= INDEXING_RETRY;
+		} else if ( strcmp(bits[i], "retry") == 0) {
+			method |= INDEXING_RETRY;
 
-		} else if ( strcmp(methods[i], "noretry") == 0) {
-			list[nmeth] &= ~INDEXING_RETRY;
+		} else if ( strcmp(bits[i], "noretry") == 0) {
+			method &= ~INDEXING_RETRY;
 
-		} else if ( strcmp(methods[i], "multi") == 0) {
-			list[nmeth] |= INDEXING_MULTI;
+		} else if ( strcmp(bits[i], "multi") == 0) {
+			method |= INDEXING_MULTI;
 
-		} else if ( strcmp(methods[i], "nomulti") == 0) {
-			list[nmeth] &= ~INDEXING_MULTI;
+		} else if ( strcmp(bits[i], "nomulti") == 0) {
+			method &= ~INDEXING_MULTI;
 
-		} else if ( strcmp(methods[i], "refine") == 0) {
-			list[nmeth] |= INDEXING_REFINE;
+		} else if ( strcmp(bits[i], "refine") == 0) {
+			method |= INDEXING_REFINE;
 
-		} else if ( strcmp(methods[i], "norefine") == 0) {
-			list[nmeth] &= ~INDEXING_REFINE;
+		} else if ( strcmp(bits[i], "norefine") == 0) {
+			method &= ~INDEXING_REFINE;
 
 		} else {
 			ERROR("Bad list of indexing methods: '%s'\n", str);
-			return NULL;
+			return INDEXING_ERROR;
 		}
 
-		free(methods[i]);
+		free(bits[i]);
 
 	}
-	free(methods);
-	list[++nmeth] = INDEXING_NONE;
+	free(bits);
+
+	if ( !have_method ) return warn_method(str);
 
-	return list;
+	return method;
 }
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 33d568f6..2107ed80 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -81,6 +81,7 @@
  * @INDEXING_DEBUG: Results injector for debugging
  * @INDEXING_ASDF: Use in-built "asdf" indexer
  * @INDEXING_TAKETWO: Use in-built "taketwo" indexer
+ * @INDEXING_ERROR: Special value for unrecognised indexing engine name
  * @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell
  *   axes for agreement with given cell.
  * @INDEXING_CHECK_CELL_AXES: Check unit cell axes for agreement with given
@@ -90,7 +91,7 @@
  * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to
  *   guide the indexing process.
  * @INDEXING_USE_CELL_PARAMETERS: Use the unit cell parameters to guide the
- *   indexingprocess.
+ *   indexing process.
  * @INDEXING_RETRY: If the indexer doesn't succeed, delete the weakest peaks
  *   and try again.
  * @INDEXING_MULTI: If the indexer succeeds, delete the peaks explained by the
@@ -115,6 +116,8 @@ typedef enum {
 	INDEXING_ASDF = 8,
 	INDEXING_TAKETWO = 9,
 
+	INDEXING_ERROR = 255,  /* Unrecognised indexing engine */
+
 	/* Bits at the top of the IndexingMethod are flags which modify the
 	 * behaviour of the indexer. */
 	INDEXING_CHECK_CELL_COMBINATIONS = 256,
@@ -143,28 +146,28 @@ extern "C" {
  * IndexingPrivate:
  *
  * This is an opaque data structure containing information needed by the
- * indexing method.
+ * indexing system.
  **/
-typedef void *IndexingPrivate;
+typedef struct _indexingprivate IndexingPrivate;
 
-extern IndexingMethod *build_indexer_list(const char *str);
+/* Convert indexing methods to/from text */
 extern char *indexer_str(IndexingMethod indm);
+extern IndexingMethod get_indm_from_string(const char *method);
 
 #include "detector.h"
 #include "cell.h"
 #include "image.h"
 
-extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
-                                          struct detector *det, float *ltl,
-                                          const char *options);
+extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell,
+                                       struct detector *det, float *ltl,
+                                       int no_refine, const char *options);
 
-extern void index_pattern(struct image *image,
-                          IndexingMethod *indms, IndexingPrivate **iprivs);
+extern void index_pattern(struct image *image, IndexingPrivate *ipriv);
 
-extern void index_pattern_2(struct image *image, IndexingMethod *indms,
-                            IndexingPrivate **iprivs, int *ping);
+extern void index_pattern_2(struct image *image, IndexingPrivate *ipriv,
+                            int *ping);
 
-extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs);
+extern void cleanup_indexing(IndexingPrivate *ipriv);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 05f853eb..a8e6e119 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -725,7 +725,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
 }
 
 
-int run_mosflm(struct image *image, IndexingPrivate *ipriv)
+int run_mosflm(struct image *image, void *ipriv)
 {
 	struct mosflm_data *mosflm;
 	unsigned int opts;
@@ -843,8 +843,8 @@ int run_mosflm(struct image *image, IndexingPrivate *ipriv)
 }
 
 
-IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
-                                struct detector *det, float *ltl)
+void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
+                     struct detector *det, float *ltl)
 {
 	struct mosflm_private *mp;
 	int need_cell = 0;
@@ -911,7 +911,7 @@ IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
 }
 
 
-void mosflm_cleanup(IndexingPrivate *pp)
+void mosflm_cleanup(void *pp)
 {
 	struct mosflm_private *p;
 	p = (struct mosflm_private *)pp;
diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h
index 572741dd..39ec5390 100644
--- a/libcrystfel/src/mosflm.h
+++ b/libcrystfel/src/mosflm.h
@@ -3,13 +3,13 @@
  *
  * Invoke the DPS auto-indexing algorithm through MOSFLM
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
  *   2010 Richard Kirian <rkirian@asu.edu>
- *   2012-2014 Thomas White <taw@physics.org>
+ *   2012-2014,2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -41,12 +41,12 @@
 extern "C" {
 #endif
 
-extern int run_mosflm(struct image *image, IndexingPrivate *ipriv);
+extern int run_mosflm(struct image *image, void *ipriv);
 
-extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
-                                       struct detector *det, float *ltl);
+extern void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
+                            struct detector *det, float *ltl);
 
-extern void mosflm_cleanup(IndexingPrivate *pp);
+extern void mosflm_cleanup(void *pp);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c
new file mode 100644
index 00000000..417465df
--- /dev/null
+++ b/libcrystfel/src/peakfinder8.c
@@ -0,0 +1,1195 @@
+/*
+ * peakfinder8.c
+ *
+ * The peakfinder8 algorithm
+ *
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017  Valerio Mariani <valerio.mariani@desy.de>
+ *   2017  Anton Barty <anton.barty@desy.de>
+ *   2017  Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+
+#include "peakfinder8.h"
+
+
+// CrystFEL-only block 1
+struct radius_maps
+{
+	float **r_maps;
+	int n_rmaps;
+};
+
+
+struct peakfinder_mask
+{
+	char **masks;
+	int n_masks;
+};
+
+
+struct peakfinder_panel_data
+{
+	float **panel_data;
+	int *panel_h;
+	int *panel_w;
+	int num_panels;
+};
+// End of CrystFEL-only block 1
+
+
+struct radial_stats
+{
+	float *roffset;
+	float *rthreshold;
+	float *lthreshold;
+	float *rsigma;
+	int *rcount;
+	int n_rad_bins;
+};
+
+
+struct peakfinder_intern_data
+{
+	char *pix_in_peak_map;
+	int *infs;
+	int *inss;
+	int *peak_pixels;
+};
+
+
+struct peakfinder_peak_data
+{
+	int num_found_peaks;
+	int *npix;
+	float *com_fs;
+	float *com_ss;
+	int *com_index;
+	float *tot_i;
+	float *max_i;
+	float *sigma;
+	float *snr;
+};
+
+
+// CrystFEL-only block 2
+static struct radius_maps *compute_radius_maps(struct detector *det)
+{
+	int i, u, iss, ifs;
+	struct panel p;
+	struct radius_maps *rm = NULL;
+
+	rm = (struct radius_maps *)malloc(sizeof(struct radius_maps));
+	if ( rm == NULL ) {
+		return NULL;
+	}
+
+	rm->r_maps = (float **)malloc(det->n_panels*sizeof(float*));
+	if ( rm->r_maps == NULL ) {
+		free(rm);
+		return NULL;
+	}
+
+	rm->n_rmaps = det->n_panels;
+
+	for( i=0 ; i<det->n_panels ; i++ ) {
+
+		p = det->panels[i];
+		rm->r_maps[i] = (float *)malloc(p.h*p.w*sizeof(float));
+
+		if ( rm->r_maps[i] == NULL ) {
+			for ( u = 0; u<i; u++ ) {
+				free(rm->r_maps[u]);
+			}
+			free(rm);
+			return NULL;
+		}
+
+		for ( iss=0 ; iss<p.h ; iss++ ) {
+			for ( ifs=0; ifs<p.w; ifs++ ) {
+
+				int rmi;
+				int x,y;
+
+				rmi = ifs + p.w * iss;
+
+				x = (p.cnx  + ifs * p.fsx + iss * p.ssx);
+				y = (p.cny  + ifs * p.fsy + iss * p.ssy);
+
+				rm->r_maps[i][rmi] = sqrt(x * x + y * y);
+			}
+		}
+	}
+	return rm;
+}
+
+
+static void free_radius_maps(struct radius_maps *r_maps)
+{
+	int i;
+
+	for ( i=0 ; i<r_maps->n_rmaps ; i++ ) {
+		free(r_maps->r_maps[i]);
+	}
+	free(r_maps->r_maps);
+	free(r_maps);
+}
+
+
+static struct peakfinder_mask *create_peakfinder_mask(struct image *img,
+                                                      struct radius_maps *rmps,
+                                                      int min_res,
+                                                      int max_res)
+{
+	int i;
+	struct peakfinder_mask *msk;
+
+	msk = (struct peakfinder_mask *)malloc(sizeof(struct peakfinder_mask));
+	msk->masks =(char **) malloc(img->det->n_panels*sizeof(char*));
+	msk->n_masks = img->det->n_panels;
+	for ( i=0; i<img->det->n_panels; i++) {
+
+		struct panel p;
+		int iss, ifs;
+
+		p = img->det->panels[i];
+
+		msk->masks[i] = (char *)calloc(p.w*p.h,sizeof(char));
+
+		for ( iss=0 ; iss<p.h ; iss++ ) {
+			for ( ifs=0 ; ifs<p.w ; ifs++ ) {
+
+				int idx;
+
+				idx = ifs + iss*p.w;
+
+				if ( rmps->r_maps[i][idx] < max_res
+				  && rmps->r_maps[i][idx] > min_res ) {
+
+					if  (! ( ( img->bad != NULL )
+					      && ( img->bad[i] != NULL )
+					      && ( img->bad[i][idx] != 0 ) ) ) {
+						msk->masks[i][idx] = 1;
+					}
+
+				}
+			}
+		}
+	}
+	return msk;
+}
+
+
+static void free_peakfinder_mask(struct peakfinder_mask * pfmask)
+{
+	int i;
+
+	for ( i=0 ; i<pfmask->n_masks ; i++ ) {
+		free(pfmask->masks[i]);
+	}
+	free(pfmask->masks);
+	free(pfmask);
+}
+
+
+static struct peakfinder_panel_data *allocate_panel_data(int num_panels)
+{
+
+	struct peakfinder_panel_data *pfdata;
+
+	pfdata = (struct peakfinder_panel_data *)malloc(sizeof(struct peakfinder_panel_data));
+	if ( pfdata == NULL ) {
+		return NULL;
+	}
+
+	pfdata->panel_h = (int *)malloc(num_panels*sizeof(int));
+	if ( pfdata->panel_h == NULL ) {
+		free(pfdata);
+		return NULL;
+	}
+
+	pfdata->panel_w = (int *)malloc(num_panels*sizeof(int));
+	if ( pfdata->panel_w == NULL ) {
+		free(pfdata->panel_h);
+		free(pfdata);
+		return NULL;
+	}
+
+	pfdata->panel_data = (float **)malloc(num_panels*sizeof(float*));
+	if ( pfdata->panel_data == NULL ) {
+		free(pfdata->panel_w);
+		free(pfdata->panel_h);
+		free(pfdata);
+		return NULL;
+	}
+
+	pfdata->num_panels = num_panels;
+
+	return pfdata;
+}
+
+
+static void free_panel_data(struct peakfinder_panel_data *pfdata)
+{
+	free(pfdata->panel_data);
+	free(pfdata->panel_w);
+	free(pfdata->panel_h);
+	free(pfdata);
+}
+
+
+static void compute_num_radial_bins(int w, int h, float *r_map, float *max_r)
+{
+	int ifs, iss;
+	int pidx;
+
+	for ( iss=0 ; iss<h ; iss++ ) {
+		for ( ifs=0 ; ifs<w ; ifs++ ) {
+			pidx = iss * w + ifs;
+			if ( r_map[pidx] > *max_r ) {
+				*max_r = r_map[pidx];
+			}
+		}
+	}
+}
+// End of CrystFEL-only block 2
+
+
+static struct radial_stats* allocate_radial_stats(int num_rad_bins)
+{
+	struct radial_stats* rstats;
+
+	rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats));
+	if ( rstats == NULL ) {
+		return NULL;
+	}
+
+	rstats->roffset = (float *)malloc(num_rad_bins*sizeof(float));
+	if ( rstats->roffset == NULL ) {
+		free(rstats);
+		return NULL;
+	}
+
+	rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float));
+	if ( rstats->rthreshold == NULL ) {
+		free(rstats);
+		free(rstats->roffset);
+		return NULL;
+	}
+
+	rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float));
+	if ( rstats->lthreshold == NULL ) {
+		free(rstats);
+		free(rstats->rthreshold);
+		free(rstats->roffset);
+		return NULL;
+	}
+
+	rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float));
+	if ( rstats->rsigma == NULL ) {
+		free(rstats);
+		free(rstats->roffset);
+		free(rstats->rthreshold);
+		free(rstats->lthreshold);
+		return NULL;
+	}
+
+	rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int));
+	if ( rstats->rcount == NULL ) {
+		free(rstats);
+		free(rstats->roffset);
+		free(rstats->rthreshold);
+		free(rstats->lthreshold);
+		free(rstats->rsigma);
+		return NULL;
+	}
+
+	rstats->n_rad_bins = num_rad_bins;
+
+	return rstats;
+}
+
+
+static void free_radial_stats(struct radial_stats *rstats)
+{
+	free(rstats->roffset);
+	free(rstats->rthreshold);
+	free(rstats->lthreshold);
+	free(rstats->rsigma);
+	free(rstats->rcount);
+	free(rstats);
+}
+
+
+static void fill_radial_bins(float *data,
+                             int w,
+                             int h,
+                             float *r_map,
+                             char *mask,
+                             float *rthreshold,
+                             float *lthreshold,
+                             float *roffset,
+                             float *rsigma,
+                             int *rcount)
+{
+	int iss, ifs;
+	int pidx;
+
+	int curr_r;
+	float value;
+
+	for ( iss = 0; iss<h ; iss++ ) {
+		for ( ifs = 0; ifs<w ; ifs++ ) {
+			pidx = iss * w + ifs;
+			if ( mask[pidx] != 0 ) {
+				curr_r = (int)rint(r_map[pidx]);
+				value = data[pidx];
+				if ( value < rthreshold[curr_r ] && value>lthreshold[curr_r]) {
+					roffset[curr_r] += value;
+					rsigma[curr_r] += (value * value);
+					rcount[curr_r] += 1;
+				}
+			}
+		}
+	}
+}
+
+
+static void compute_radial_stats(float *rthreshold,
+                                 float *lthreshold,
+                                 float *roffset,
+                                 float *rsigma,
+                                 int *rcount,
+                                 int num_rad_bins,
+                                 float min_snr,
+                                 float acd_threshold)
+{
+	int ri;
+	float this_offset, this_sigma;
+
+	for ( ri=0 ; ri<num_rad_bins ; ri++ ) {
+
+		if ( rcount[ri] == 0 ) {
+			roffset[ri] = 0;
+			rsigma[ri] = 0;
+			rthreshold[ri] = 1e9;
+			lthreshold[ri] = -1e9;
+		} else {
+			this_offset = roffset[ri] / rcount[ri];
+			this_sigma = rsigma[ri] / rcount[ri] - (this_offset * this_offset);
+			if ( this_sigma >= 0 ) {
+				this_sigma = sqrt(this_sigma);
+			}
+
+			roffset[ri] = this_offset;
+			rsigma[ri] = this_sigma;
+			rthreshold[ri] = roffset[ri] + min_snr*rsigma[ri];
+			lthreshold[ri] = roffset[ri] - min_snr*rsigma[ri];
+
+			if ( rthreshold[ri] < acd_threshold ) {
+				rthreshold[ri] = acd_threshold;
+			}
+		}
+	}
+
+}
+
+
+struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks)
+{
+	struct peakfinder_peak_data *pkdata;
+
+	pkdata = (struct peakfinder_peak_data*)malloc(sizeof(struct peakfinder_peak_data));
+	if ( pkdata == NULL ) {
+		return NULL;
+	}
+
+	pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int));
+	if ( pkdata->npix == NULL ) {
+		free(pkdata);
+		free(pkdata->npix);
+		return NULL;
+	}
+
+	pkdata->com_fs = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->com_fs == NULL ) {
+		free(pkdata->npix);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->com_ss = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->com_ss == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->com_index = (int *)malloc(max_num_peaks*sizeof(int));
+	if ( pkdata->com_ss == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata->com_ss);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->tot_i == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata->com_ss);
+		free(pkdata->com_index);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->max_i = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->max_i == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata->com_ss);
+		free(pkdata->com_index);
+		free(pkdata->tot_i);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->sigma = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->sigma == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata->com_ss);
+		free(pkdata->com_index);
+		free(pkdata->tot_i);
+		free(pkdata->max_i);
+		free(pkdata);
+		return NULL;
+	}
+
+	pkdata->snr = (float *)malloc(max_num_peaks*sizeof(float));
+	if ( pkdata->snr == NULL ) {
+		free(pkdata->npix);
+		free(pkdata->com_fs);
+		free(pkdata->com_ss);
+		free(pkdata->com_index);
+		free(pkdata->tot_i);
+		free(pkdata->max_i);
+		free(pkdata->sigma);
+		free(pkdata);
+		return NULL;
+	}
+
+	return pkdata;
+}
+
+
+static void free_peak_data(struct peakfinder_peak_data *pkdata) {
+	free(pkdata->npix);
+	free(pkdata->com_fs);
+	free(pkdata->com_ss);
+	free(pkdata->com_index);
+	free(pkdata->tot_i);
+	free(pkdata->max_i);
+	free(pkdata->sigma);
+	free(pkdata->snr);
+	free(pkdata);
+}
+
+
+static struct peakfinder_intern_data *allocate_peakfinder_intern_data(int data_size,
+                                                                      int max_pix_count)
+{
+
+	struct peakfinder_intern_data *intern_data;
+
+	intern_data = (struct peakfinder_intern_data *)malloc(sizeof(struct peakfinder_intern_data));
+	if ( intern_data == NULL ) {
+		return NULL;
+	}
+
+	intern_data->pix_in_peak_map =(char *)calloc(data_size, sizeof(char));
+	if ( intern_data->pix_in_peak_map == NULL ) {
+		free(intern_data);
+		return NULL;
+	}
+
+	intern_data->infs =(int *)calloc(data_size, sizeof(int));
+	if ( intern_data->infs == NULL ) {
+		free(intern_data->pix_in_peak_map);
+		free(intern_data);
+		return NULL;
+	}
+
+	intern_data->inss =(int *)calloc(data_size, sizeof(int));
+	if ( intern_data->inss == NULL ) {
+		free(intern_data->pix_in_peak_map);
+		free(intern_data->infs);
+		free(intern_data);
+		return NULL;
+	}
+
+	intern_data->peak_pixels =(int *)calloc(max_pix_count, sizeof(int));
+	if ( intern_data->peak_pixels == NULL ) {
+		free(intern_data->pix_in_peak_map);
+		free(intern_data->infs);
+		free(intern_data->inss);
+		free(intern_data);
+		return NULL;
+	}
+
+	return intern_data;
+}
+
+
+static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid)
+{
+	free(pfid->peak_pixels);
+	free(pfid->pix_in_peak_map);
+	free(pfid->infs);
+	free(pfid->inss);
+	free(pfid);
+}
+
+
+
+static void peak_search(int p,
+                        struct peakfinder_intern_data *pfinter,
+                        float *copy, char *mask, float *r_map,
+                        float *rthreshold, float *roffset,
+                        int *num_pix_in_peak, int asic_size_fs,
+                        int asic_size_ss, int aifs, int aiss,
+                        int num_pix_fs, float *sum_com_fs,
+                        float *sum_com_ss, float *sum_i, int max_pix_count)
+{
+
+	int k, pi;
+	int curr_radius;
+	float curr_threshold;
+	int curr_fs;
+	int curr_ss;
+	float curr_i;
+
+	int search_fs[9] = { 0, -1, 0, 1, -1, 1, -1, 0, 1 };
+	int search_ss[9] = { 0, -1, -1, -1, 0, 0, 1, 1, 1 };
+	int search_n = 9;
+
+	// Loop through search pattern
+	for ( k=0; k<search_n; k++ ) {
+
+		if ( (pfinter->infs[p] + search_fs[k]) < 0 ) continue;
+		if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue;
+		if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue;
+		if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue;
+
+		// Neighbour point in big array
+		curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs;
+		curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss;
+		pi = curr_fs + curr_ss * num_pix_fs;
+
+		curr_radius = (int)rint(r_map[pi]);
+		curr_threshold = rthreshold[curr_radius];
+
+		// Above threshold?
+		if ( copy[pi] > curr_threshold
+		  && pfinter->pix_in_peak_map[pi] == 0
+		  && mask[pi] != 0 ) {
+
+			curr_i = copy[pi] - roffset[curr_radius];
+			*sum_i += curr_i;
+			*sum_com_fs += curr_i * ((float)curr_fs);  // for center of mass x
+			*sum_com_ss += curr_i * ((float)curr_ss);  // for center of mass y
+
+			pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + search_ss[k];
+			pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + search_fs[k];
+			pfinter->pix_in_peak_map[pi] = 1;
+			if ( *num_pix_in_peak < max_pix_count ) {
+				  pfinter->peak_pixels[*num_pix_in_peak] = pi;
+			}
+			*num_pix_in_peak = *num_pix_in_peak + 1;
+		}
+	}
+}
+
+
+static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int,
+                           float *copy, float *r_map,
+                           float *rthreshold, float *roffset,
+                           char *pix_in_peak_map, char *mask, int asic_size_fs,
+                           int asic_size_ss, int aifs, int aiss,
+                           int num_pix_fs,float *local_sigma, float *local_offset,
+                           float *background_max_i, int com_idx,
+                           int local_bg_radius)
+{
+	int ssj, fsi;
+	float pix_radius;
+	int curr_fs, curr_ss;
+	int pi;
+	int curr_radius;
+	float curr_threshold;
+	float curr_i;
+
+	int np_sigma;
+	int np_counted;
+	int local_radius;
+
+	float sum_i;
+	float sum_i_squared;
+
+	ring_width = 2 * local_bg_radius;
+
+	sum_i = 0;
+	sum_i_squared = 0;
+	np_sigma = 0;
+	np_counted = 0;
+	local_radius = 0;
+
+	for ( ssj = -ring_width ; ssj<ring_width ; ssj++ ) {
+		for ( fsi = -ring_width ; fsi<ring_width ; fsi++ ) {
+
+			// Within-ASIC check
+			if ( (com_fs_int + fsi) < 0 ) continue;
+			if ( (com_fs_int + fsi) >= asic_size_fs ) continue;
+			if ( (com_ss_int + ssj) < 0 ) continue;
+			if ( (com_ss_int + ssj) >= asic_size_ss )
+			continue;
+
+			// Within outer ring check
+			pix_radius = sqrt(fsi * fsi + ssj * ssj);
+			if ( pix_radius>ring_width ) continue;
+
+			// Position of this point in data stream
+			curr_fs = com_fs_int + fsi + aifs * asic_size_fs;
+			curr_ss = com_ss_int + ssj + aiss * asic_size_ss;
+			pi = curr_fs + curr_ss * num_pix_fs;
+
+			curr_radius = (int)rint(r_map[pi]);
+			curr_threshold = rthreshold[curr_radius];
+
+			// Intensity above background ??? just intensity?
+			curr_i = copy[pi];
+
+			// Keep track of value and value-squared for offset and sigma calculation
+			if ( curr_i < curr_threshold && pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) {
+
+				np_sigma++;
+				sum_i += curr_i;
+				sum_i_squared += (curr_i * curr_i);
+
+				if ( curr_i > *background_max_i ) {
+					*background_max_i = curr_i;
+				}
+			}
+			np_counted += 1;
+		}
+	}
+
+	// Calculate local background and standard deviation
+	if ( np_sigma != 0 ) {
+		*local_offset = sum_i / np_sigma;
+		*local_sigma = sum_i_squared / np_sigma - (*local_offset * *local_offset);
+		if (*local_sigma >= 0) {
+			*local_sigma = sqrt(*local_sigma);
+		} else {
+			*local_sigma = 0.01;
+		}
+	} else {
+		local_radius = (int)rint(r_map[(int)rint(com_idx)]);
+		*local_offset = roffset[local_radius];
+		*local_sigma = 0.01;
+	}
+}
+
+
+static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs,
+                          int aiss, int aifs, float *rthreshold,
+                          float *roffset, int *peak_count,
+                          float *copy, struct peakfinder_intern_data *pfinter,
+                          float *r_map, char *mask, int *npix, float *com_fs,
+                          float *com_ss, int *com_index, float *tot_i,
+                          float *max_i, float *sigma, float *snr,
+                          int min_pix_count, int max_pix_count,
+                          int local_bg_radius, float min_snr, int max_n_peaks)
+{
+	int pxss, pxfs;
+	int num_pix_in_peak;
+
+	// Loop over pixels within a module
+	for ( pxss=1 ; pxss<asic_size_ss-1 ; pxss++ ) {
+		for ( pxfs=1 ; pxfs<asic_size_fs-1 ; pxfs++ ) {
+
+			float curr_thresh;
+			int pxidx;
+			int curr_rad;
+
+			pxidx = (pxss + aiss * asic_size_ss) * num_pix_fs +
+			pxfs + aifs * asic_size_fs;
+
+			curr_rad = (int)rint(r_map[pxidx]);
+			curr_thresh = rthreshold[curr_rad];
+
+			if ( copy[pxidx] > curr_thresh
+			  && pfinter->pix_in_peak_map[pxidx] == 0
+			  && mask[pxidx] != 0 ) {   //??? not sure if needed
+
+				// This might be the start of a new peak - start searching
+				float sum_com_fs, sum_com_ss;
+				float sum_i;
+				float peak_com_fs, peak_com_ss;
+				float peak_com_fs_int, peak_com_ss_int;
+				float peak_tot_i, pk_tot_i_raw;
+				float peak_max_i, pk_max_i_raw;
+				float peak_snr;
+				float local_sigma, local_offset;
+				float background_max_i;
+				int lt_num_pix_in_pk;
+				int ring_width;
+				int peak_idx;
+				int com_idx;
+				int p;
+
+				pfinter->infs[0] = pxfs;
+				pfinter->inss[0] = pxss;
+				pfinter->peak_pixels[0] = pxidx;
+				num_pix_in_peak = 0; //y 1;
+
+				sum_i = 0;
+				sum_com_fs = 0;
+				sum_com_ss = 0;
+
+				// Keep looping until the pixel count within this peak does not change
+				do {
+					lt_num_pix_in_pk = num_pix_in_peak;
+
+					// Loop through points known to be within this peak
+					for ( p=0; p<=num_pix_in_peak; p++ ) { //changed from 1 to 0 by O.Y.
+						peak_search(p,
+						            pfinter, copy, mask,
+						            r_map,
+						            rthreshold,
+						            roffset,
+						            &num_pix_in_peak,
+						            asic_size_fs,
+						            asic_size_ss,
+						            aifs, aiss,
+						            num_pix_fs,
+						            &sum_com_fs,
+						            &sum_com_ss,
+						            &sum_i,
+						            max_pix_count);
+					}
+
+				} while ( lt_num_pix_in_pk != num_pix_in_peak );
+
+				// Too many or too few pixels means ignore this 'peak'; move on now
+				if ( num_pix_in_peak < min_pix_count || num_pix_in_peak > max_pix_count ) continue;
+
+				// If for some reason sum_i is 0 - it's better to skip
+				if ( fabs(sum_i) < 1e-10 ) continue;
+
+				// Calculate center of mass for this peak from initial peak search
+				peak_com_fs = sum_com_fs / fabs(sum_i);
+				peak_com_ss = sum_com_ss / fabs(sum_i);
+
+				com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * num_pix_fs;
+
+				peak_com_fs_int = (int)rint(peak_com_fs) - aifs * asic_size_fs;
+				peak_com_ss_int = (int)rint(peak_com_ss) - aiss * asic_size_ss;
+
+				// Calculate the local signal-to-noise ratio and local background in an annulus around
+				// this peak (excluding pixels which look like they might be part of another peak)
+				local_sigma = 0.0;
+				local_offset = 0.0;
+				background_max_i = 0.0;
+
+				ring_width = 2 * local_bg_radius;
+
+				search_in_ring(ring_width, peak_com_fs_int,
+				               peak_com_ss_int,
+				               copy, r_map, rthreshold,
+				               roffset,
+				               pfinter->pix_in_peak_map,
+				               mask, asic_size_fs,
+				               asic_size_ss,
+				               aifs, aiss,
+				               num_pix_fs,
+				               &local_sigma,
+				               &local_offset,
+				               &background_max_i,
+				               com_idx, local_bg_radius);
+
+				// Re-integrate (and re-centroid) peak using local background estimates
+				peak_tot_i = 0;
+				pk_tot_i_raw = 0;
+				peak_max_i = 0;
+				pk_max_i_raw = 0;
+				sum_com_fs = 0;
+				sum_com_ss = 0;
+
+				for ( peak_idx = 0 ;
+					peak_idx < num_pix_in_peak && peak_idx < max_pix_count ;
+					peak_idx++ ) {
+
+					int curr_idx;
+					float curr_i;
+					float curr_i_raw;
+					int curr_fs, curr_ss;
+
+					curr_idx = pfinter->peak_pixels[peak_idx];
+					curr_i_raw = copy[curr_idx];
+					curr_i = curr_i_raw - local_offset;
+					peak_tot_i += curr_i;
+					pk_tot_i_raw += curr_i_raw;
+
+					// Remember that curr_idx = curr_fs + curr_ss*num_pix_fs
+					curr_fs = curr_idx % num_pix_fs;
+					curr_ss = curr_idx / num_pix_fs;
+					sum_com_fs += curr_i * ((float)curr_fs);
+					sum_com_ss += curr_i * ((float)curr_ss);
+
+					if ( curr_i_raw > pk_max_i_raw ) pk_max_i_raw = curr_i_raw;
+					if ( curr_i > peak_max_i ) peak_max_i = curr_i;
+				}
+
+
+				// This CAN happen! Better to skip...
+				if ( fabs(peak_tot_i) < 1e-10 ) continue;
+
+				peak_com_fs = sum_com_fs / fabs(peak_tot_i);
+				peak_com_ss = sum_com_ss / fabs(peak_tot_i);
+
+				// Calculate signal-to-noise and apply SNR criteria
+				if ( fabs(local_sigma) > 1e-10 ) {
+					peak_snr = peak_tot_i / local_sigma;
+				} else {
+					peak_snr = 0;
+				}
+
+				if (peak_snr < min_snr) continue;
+
+				// Is the maximum intensity in the peak enough above intensity in background region to
+				// be a peak and not noise? The more pixels there are in the peak, the more relaxed we
+				// are about this criterion
+				//f_background_thresh = background_max_i - local_offset; //!!! Ofiget'!  If I uncomment
+				// if (peak_max_i < f_background_thresh) {               // these lines the result is
+				// different!
+				if (peak_max_i < background_max_i - local_offset) continue;
+
+				// This is a peak? If so, add info to peak list
+				if ( num_pix_in_peak >= min_pix_count
+				  && num_pix_in_peak <= max_pix_count ) {
+
+					// Bragg peaks in the mask
+					for ( peak_idx = 0 ;
+					      peak_idx < num_pix_in_peak &&
+					      peak_idx < max_pix_count ;
+					      peak_idx++ ) {
+						pfinter->pix_in_peak_map[pfinter->peak_pixels[peak_idx]] = 2;
+					}
+
+					int peak_com_idx;
+					peak_com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) *
+						                num_pix_fs;
+					// Remember peak information
+					if ( *peak_count < max_n_peaks ) {
+
+						int pidx;
+						pidx = *peak_count;
+
+						npix[pidx] = num_pix_in_peak;
+						com_fs[pidx] = peak_com_fs;
+						com_ss[pidx] = peak_com_ss;
+						com_index[pidx] = peak_com_idx;
+						tot_i[pidx] = peak_tot_i;
+						max_i[pidx] = peak_max_i;
+						sigma[pidx] = local_sigma;
+						snr[pidx] = peak_snr;
+					}
+					*peak_count += 1;
+				}
+			}
+		}
+	}
+}
+
+
+static int peakfinder8_base(float *roffset, float *rthreshold,
+                            float *data, char *mask, float *r_map,
+                            int asic_size_fs, int num_asics_fs,
+                            int asic_size_ss, int num_asics_ss,
+                            int max_n_peaks, int *num_found_peaks,
+                            int *npix, float *com_fs,
+                            float *com_ss, int *com_index, float *tot_i,
+                            float *max_i, float *sigma, float *snr,
+                            int min_pix_count, int max_pix_count,
+                            int local_bg_radius, float min_snr,
+                            char* outliersMask)
+{
+
+	int num_pix_fs, num_pix_ss, num_pix_tot;
+	int aifs, aiss;
+	int peak_count;
+	struct peakfinder_intern_data *pfinter;
+
+	num_pix_fs = asic_size_fs * num_asics_fs;
+	num_pix_ss = asic_size_ss * num_asics_ss;
+	num_pix_tot = num_pix_fs * num_pix_ss;
+
+	pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count);
+	if ( pfinter == NULL ) {
+		return 1;
+	}
+
+	peak_count = 0;
+
+	// Loop over modules (nxn array)
+	for ( aiss=0 ; aiss<num_asics_ss ; aiss++ ) {
+		for ( aifs=0 ; aifs<num_asics_fs ; aifs++ ) {                 // ??? to change to proper panels need
+			process_panel(asic_size_fs, asic_size_ss, num_pix_fs, // change copy, mask, r_map
+			              aiss, aifs, rthreshold, roffset,
+			              &peak_count, data, pfinter, r_map, mask,
+			              npix, com_fs, com_ss, com_index, tot_i,
+			              max_i, sigma, snr, min_pix_count,
+			              max_pix_count, local_bg_radius, min_snr,
+			              max_n_peaks);
+		}
+	}
+	*num_found_peaks = peak_count;
+
+	if (outliersMask != NULL) {
+		memcpy(outliersMask, pfinter->pix_in_peak_map, num_pix_tot*sizeof(char));
+	}
+
+	free_peakfinder_intern_data(pfinter);
+
+	return 0;
+}
+
+
+int peakfinder8(struct image *img, int max_n_peaks,
+                float threshold, float min_snr,
+                int min_pix_count, int max_pix_count,
+                int local_bg_radius, int min_res,
+                int max_res, int use_saturated)
+{
+	struct radius_maps *rmaps;
+	struct peakfinder_mask *pfmask;
+	struct peakfinder_panel_data *pfdata;
+	struct radial_stats *rstats;
+	struct peakfinder_peak_data *pkdata;
+	int num_rad_bins;
+	int pi;
+	int i, it_counter;
+	int num_found_peaks;
+	int remaining_max_num_peaks;
+	int iterations;
+	float max_r;
+
+	iterations = 5;
+
+	if ( img-> det == NULL) {
+		return 1;
+	}
+
+	rmaps = compute_radius_maps(img->det);
+	if ( rmaps == NULL ) {
+		return 1;
+	}
+
+	pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res);
+	if ( pfmask == NULL ) {
+		free_radius_maps(rmaps);
+		return 1;
+	}
+
+	pfdata = allocate_panel_data(img->det->n_panels);
+	if ( pfdata == NULL) {
+		free_radius_maps(rmaps);
+		free_peakfinder_mask(pfmask);
+		return 1;
+	}
+
+	for ( pi=0 ; pi<img->det->n_panels ; pi++ ) {
+		pfdata->panel_h[pi] = img->det->panels[pi].h;
+		pfdata->panel_w[pi] = img->det->panels[pi].w;
+		pfdata->panel_data[pi] = img->dp[pi];
+		pfdata->num_panels = img->det->n_panels;
+	}
+
+	max_r = -1e9;
+
+	for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
+
+		compute_num_radial_bins(pfdata->panel_w[pi],
+		                        pfdata->panel_h[pi],
+		                        rmaps->r_maps[pi],
+		                        &max_r);
+	}
+
+	num_rad_bins = (int)ceil(max_r) + 1;
+
+	rstats = allocate_radial_stats(num_rad_bins);
+	if ( rstats == NULL ) {
+		free_radius_maps(rmaps);
+		free_peakfinder_mask(pfmask);
+		free_panel_data(pfdata);
+		return 1;
+	}
+
+	for ( i=0 ; i<rstats->n_rad_bins ; i++) {
+		rstats->rthreshold[i] = 1e9;
+	}
+
+	for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) {
+
+		for ( i=0; i<num_rad_bins; i++ ) {
+			rstats->roffset[i] = 0;
+			rstats->rsigma[i] = 0;
+			rstats->rcount[i] = 0;
+		}
+
+		for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
+
+			fill_radial_bins(pfdata->panel_data[pi],
+			                 pfdata->panel_w[pi],
+			                 pfdata->panel_h[pi],
+			                 rmaps->r_maps[pi],
+			                 pfmask->masks[pi],
+			                 rstats->rthreshold,
+			                 rstats->lthreshold,
+			                 rstats->roffset,
+			                 rstats->rsigma,
+			                 rstats->rcount);
+
+		}
+
+		compute_radial_stats(rstats->rthreshold,
+		                     rstats->lthreshold,
+		                     rstats->roffset,
+		                     rstats->rsigma,
+		                     rstats->rcount,
+		                     num_rad_bins,
+		                     min_snr,
+		                     threshold);
+
+	}
+
+	pkdata = allocate_peak_data(max_n_peaks);
+	if ( pkdata == NULL ) {
+		free_radius_maps(rmaps);
+		free_peakfinder_mask(pfmask);
+		free_panel_data(pfdata);
+		free_radial_stats(rstats);
+		return 1;
+	}
+
+	remaining_max_num_peaks = max_n_peaks;
+
+	for ( pi=0 ; pi<img->det->n_panels ; pi++) {
+
+		int peaks_to_add;
+		int pki;
+		int ret;
+
+		num_found_peaks = 0;
+
+		if ( img->det->panels[pi].no_index ) {
+			continue;
+		}
+
+		ret = peakfinder8_base(rstats->roffset,
+		                       rstats->rthreshold,
+		                       pfdata->panel_data[pi],
+		                       pfmask->masks[pi],
+		                       rmaps->r_maps[pi],
+		                       pfdata->panel_w[pi], 1,
+		                       pfdata->panel_h[pi], 1,
+		                       max_n_peaks,
+		                       &num_found_peaks,
+		                       pkdata->npix,
+		                       pkdata->com_fs,
+		                       pkdata->com_ss,
+		                       pkdata->com_index,
+		                       pkdata->tot_i,
+		                       pkdata->max_i,
+		                       pkdata->sigma,
+		                       pkdata->snr,
+		                       min_pix_count,
+		                       max_pix_count,
+		                       local_bg_radius,
+		                       min_snr,
+		                       NULL);
+
+		if ( ret != 0 ) {
+			free_radius_maps(rmaps);
+			free_peakfinder_mask(pfmask);
+			free_panel_data(pfdata);
+			free_radial_stats(rstats);
+			return 1;
+		}
+
+		peaks_to_add = num_found_peaks;
+
+		if ( num_found_peaks > remaining_max_num_peaks ) {
+			peaks_to_add = remaining_max_num_peaks;
+		}
+
+		remaining_max_num_peaks -= peaks_to_add;
+
+		for ( pki=0 ; pki<peaks_to_add ; pki++ ) {
+
+			struct panel *p;
+
+			p = &img->det->panels[pi];
+
+			img->num_peaks += 1;
+			if ( pkdata->max_i[pki] > p->max_adu ) {
+				img->num_saturated_peaks++;
+				if ( !use_saturated ) {
+					continue;
+				}
+			}
+
+			image_add_feature(img->features,
+			                  pkdata->com_fs[pki]+0.5,
+			                  pkdata->com_ss[pki]+0.5,
+			                  p,
+			                  img,
+			                  pkdata->tot_i[pki],
+			                  NULL);
+		}
+	}
+
+	free_radius_maps(rmaps);
+	free_peakfinder_mask(pfmask);
+	free_panel_data(pfdata);
+	free_radial_stats(rstats);
+	free_peak_data(pkdata);
+	return 0;
+}
diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h
new file mode 100644
index 00000000..483eebdf
--- /dev/null
+++ b/libcrystfel/src/peakfinder8.h
@@ -0,0 +1,50 @@
+/*
+ * peakfinder8.h
+ *
+ * The peakfinder8 algorithm
+ *
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017  Valerio Mariani <valerio.mariani@desy.de>
+ *   2017  Anton Barty <anton.barty@desy.de>
+ *   2017  Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef PEAKFINDER8_H
+#define PEAKFINDER8_H
+
+#include "image.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int peakfinder8(struct image *img, int max_n_peaks,
+                       float threshold, float min_snr,
+                       int mix_pix_count, int max_pix_count,
+                       int local_bg_radius, int min_res,
+                       int max_res, int use_saturated);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // PEAKFINDER8_H
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index cb2a9f5a..28c79538 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -3,7 +3,7 @@
  *
  * Peak search and other image analysis
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
@@ -12,6 +12,7 @@
  *   2012      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *   2011      Andrew Martin <andrew.martin@desy.de>
  *   2011      Richard Kirian
+ *   2017      Valerio Mariani <valerio.mariani@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -52,6 +53,7 @@
 #include "reflist-utils.h"
 #include "cell-utils.h"
 #include "geometry.h"
+#include "peakfinder8.h"
 
 
 static int cull_peaks_in_panel(struct image *image, struct panel *p)
@@ -517,8 +519,8 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 
 
 void search_peaks(struct image *image, float threshold, float min_gradient,
-                  float min_snr, double ir_inn, double ir_mid, double ir_out,
-                  int use_saturated)
+                  float min_snr, double ir_inn, double ir_mid,
+                  double ir_out, int use_saturated)
 {
 	int i;
 
@@ -541,6 +543,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 }
 
 
+int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
+                              float threshold, float min_snr,
+                              int min_pix_count, int max_pix_count,
+                              int local_bg_radius, int min_res,
+                              int max_res, int use_saturated)
+{
+	if ( image->features != NULL ) {
+		image_feature_list_free(image->features);
+	}
+	image->features = image_feature_list_new();
+	image->num_peaks = 0;
+	image->num_saturated_peaks = 0;
+
+	return peakfinder8(image, max_n_peaks, threshold, min_snr,
+	                   min_pix_count, max_pix_count,
+	                   local_bg_radius, min_res,
+	                   max_res, use_saturated);
+}
+
+
 int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst)
 {
 	int n_feat = 0;
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index ba724419..a5095127 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -3,11 +3,12 @@
  *
  * Peak search and other image analysis
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2010-2015 Thomas White <taw@physics.org>
+ *   2017      Valerio Mariani <valerio.mariani@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -45,9 +46,14 @@ extern "C" {
 extern int *make_BgMask(struct image *image, struct panel *p, double ir_inn);
 
 extern void search_peaks(struct image *image, float threshold,
-                         float min_gradient, float min_snr,
-                         double ir_inn, double ir_mid, double ir_out,
-                         int use_saturated);
+                         float min_gradient, float min_snr, double ir_inn,
+                         double ir_mid, double ir_out, int use_saturated);
+
+extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
+                                    float threshold, float min_snr,
+                                    int mix_pix_count, int max_pix_count,
+                                    int local_bg_radius, int min_res,
+                                    int max_res, int use_saturated);
 
 extern int peak_sanity_check(struct image *image, Crystal **crystals,
                              int n_cryst);
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 4621f4f4..380a1b8a 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -27,8 +27,6 @@
  *
  */
 
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
 #include <math.h>
 #include <stdio.h>
 #include <assert.h>
diff --git a/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c
index 56273fdb..ccf35194 100644
--- a/libcrystfel/src/statistics.c
+++ b/libcrystfel/src/statistics.c
@@ -30,7 +30,6 @@
 #include <config.h>
 #endif
 
-#define _ISOC99_SOURCE
 #include <math.h>
 #include <stdlib.h>
 #include <gsl/gsl_errno.h>
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index e858172e..fb4b0c70 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -806,8 +806,8 @@ static int write_crystal(Stream *st, Crystal *cr, int include_reflections)
 }
 
 
-int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
-                int include_peaks, int include_reflections, struct event* ev)
+int write_chunk(Stream *st, struct image *i, struct imagefile *imfile,
+                int include_peaks, int include_reflections, struct event *ev)
 {
 	int j;
 	char *indexer;
@@ -832,7 +832,7 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
 	fprintf(st->fh, "beam_divergence = %.2e rad\n", i->div);
 	fprintf(st->fh, "beam_bandwidth = %.2e (fraction)\n", i->bw);
 
-	copy_hdf5_fields(hdfile, i->copyme, st->fh, ev);
+	imagefile_copy_fields(imfile, i->copyme, st->fh, ev);
 
 	if ( i->det != NULL ) {
 
@@ -1196,13 +1196,9 @@ int read_chunk_2(Stream *st, struct image *image,  StreamReadFlags srf)
 		}
 
 		if ( strncmp(line, "indexed_by = ", 13) == 0 ) {
-			IndexingMethod *list;
-			list = build_indexer_list(line+13);
-			if ( list == NULL ) {
+			image->indexed_by = get_indm_from_string(line+13);
+			if ( image->indexed_by == INDEXING_ERROR ) {
 				ERROR("Failed to read indexer list\n");
-			} else {
-				image->indexed_by = list[0];
-				free(list);
 			}
 		}
 
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index ff8628c0..764e3e36 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -3,11 +3,11 @@
  *
  * Stream tools
  *
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
  *   2011      Andrew Aquila
  *
@@ -39,6 +39,7 @@
 struct image;
 struct hdfile;
 struct event;
+struct imagefile;
 #include "cell.h"
 
 #define GEOM_START_MARKER "----- Begin geometry file -----"
@@ -106,10 +107,15 @@ extern int read_chunk(Stream *st, struct image *image);
 extern int read_chunk_2(Stream *st, struct image *image,
                            StreamReadFlags srf);
 
-extern int write_chunk(Stream *st, struct image *image, struct hdfile *hdfile,
+extern int write_chunk(Stream *st, struct image *image, struct imagefile *imfile,
                        int include_peaks, int include_reflections,
                        struct event *ev);
 
+extern int write_chunk_2(Stream *st, struct image *image,
+                         struct imagefile *imfile,
+                         int include_peaks, int include_reflections,
+                         struct event *ev);
+
 extern void write_command(Stream *st, int argc, char *argv[]);
 extern void write_geometry_file(Stream *st, const char *geom_filename);
 extern int rewind_stream(Stream *st);
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 461da0cd..90d81a71 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -387,144 +387,6 @@ int poisson_noise(gsl_rng *rng, double expected)
 }
 
 
-/**
- * SECTION:quaternion
- * @short_description: Simple quaternion handling
- * @title: Quaternion
- * @section_id:
- * @see_also:
- * @include: "utils.h"
- * @Image:
- *
- * There is a simple quaternion structure in CrystFEL.  At the moment, it is
- * only used when simulating patterns, as an argument to cell_rotate() to
- * orient the unit cell.
- */
-
-/**
- * quaternion_modulus:
- * @q: A %quaternion
- *
- * If a quaternion represents a pure rotation, its modulus should be unity.
- *
- * Returns: the modulus of the given quaternion.
- **/
-double quaternion_modulus(struct quaternion q)
-{
-	return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
-}
-
-
-/**
- * normalise_quaternion:
- * @q: A %quaternion
- *
- * Rescales the quaternion such that its modulus is unity.
- *
- * Returns: the normalised version of @q
- **/
-struct quaternion normalise_quaternion(struct quaternion q)
-{
-	double mod;
-	struct quaternion r;
-
-	mod = quaternion_modulus(q);
-
-	r.w = q.w / mod;
-	r.x = q.x / mod;
-	r.y = q.y / mod;
-	r.z = q.z / mod;
-
-	return r;
-}
-
-
-/**
- * random_quaternion:
- * @rng: A GSL random number generator to use
- *
- * Returns: a randomly generated, normalised, quaternion.
- **/
-struct quaternion random_quaternion(gsl_rng *rng)
-{
-	struct quaternion q;
-
-	q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
-	q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
-	q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
-	q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
-	q = normalise_quaternion(q);
-
-	return q;
-}
-
-
-/**
- * quaternion_valid:
- * @q: A %quaternion
- *
- * Checks if the given quaternion is normalised.
- *
- * This function performs a nasty floating point comparison of the form
- * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
- * relied upon to spot anything other than the most obvious input error.
- *
- * Returns: 1 if the quaternion is normalised, 0 if not.
- **/
-int quaternion_valid(struct quaternion q)
-{
-	double qmod;
-
-	qmod = quaternion_modulus(q);
-
-	/* Modulus = 1 to within some tolerance?
-	 * Nasty allowance for floating-point accuracy follows... */
-	if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
-
-	return 0;
-}
-
-
-/**
- * quat_rot
- * @q: A vector (in the form of a "struct rvec")
- * @z: A %quaternion
- *
- * Rotates a vector according to a quaternion.
- *
- * Returns: A rotated version of @p.
- **/
-struct rvec quat_rot(struct rvec q, struct quaternion z)
-{
-	struct rvec res;
-	double t01, t02, t03, t11, t12, t13, t22, t23, t33;
-
-	t01 = z.w*z.x;
-	t02 = z.w*z.y;
-	t03 = z.w*z.z;
-	t11 = z.x*z.x;
-	t12 = z.x*z.y;
-	t13 = z.x*z.z;
-	t22 = z.y*z.y;
-	t23 = z.y*z.z;
-	t33 = z.z*z.z;
-
-	res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
-	            + (2.0 * (t12 + t03)) * q.v
-	            + (2.0 * (t13 - t02)) * q.w;
-
-	res.v =       (2.0 * (t12 - t03)) * q.u
-	      + (1.0 - 2.0 * (t11 + t33)) * q.v
-	            + (2.0 * (t01 + t23)) * q.w;
-
-	res.w =       (2.0 * (t02 + t13)) * q.u
-	            + (2.0 * (t23 - t01)) * q.v
-	      + (1.0 - 2.0 * (t11 + t22)) * q.w;
-
-	return res;
-}
-
-
 /* Return non-zero if c is in delims */
 static int assplode_isdelim(const char c, const char *delims)
 {
@@ -691,3 +553,141 @@ void utils_fudge_gslcblas()
 {
         STATUS("%p\n", cblas_sgemm);
 }
+
+
+/**
+ * SECTION:quaternion
+ * @short_description: Simple quaternion handling
+ * @title: Quaternion
+ * @section_id:
+ * @see_also:
+ * @include: "utils.h"
+ * @Image:
+ *
+ * There is a simple quaternion structure in CrystFEL.  At the moment, it is
+ * only used when simulating patterns, as an argument to cell_rotate() to
+ * orient the unit cell.
+ */
+
+/**
+ * quaternion_modulus:
+ * @q: A %quaternion
+ *
+ * If a quaternion represents a pure rotation, its modulus should be unity.
+ *
+ * Returns: the modulus of the given quaternion.
+ **/
+double quaternion_modulus(struct quaternion q)
+{
+	return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
+}
+
+
+/**
+ * normalise_quaternion:
+ * @q: A %quaternion
+ *
+ * Rescales the quaternion such that its modulus is unity.
+ *
+ * Returns: the normalised version of @q
+ **/
+struct quaternion normalise_quaternion(struct quaternion q)
+{
+	double mod;
+	struct quaternion r;
+
+	mod = quaternion_modulus(q);
+
+	r.w = q.w / mod;
+	r.x = q.x / mod;
+	r.y = q.y / mod;
+	r.z = q.z / mod;
+
+	return r;
+}
+
+
+/**
+ * random_quaternion:
+ * @rng: A GSL random number generator to use
+ *
+ * Returns: a randomly generated, normalised, quaternion.
+ **/
+struct quaternion random_quaternion(gsl_rng *rng)
+{
+	struct quaternion q;
+
+	q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q = normalise_quaternion(q);
+
+	return q;
+}
+
+
+/**
+ * quaternion_valid:
+ * @q: A %quaternion
+ *
+ * Checks if the given quaternion is normalised.
+ *
+ * This function performs a nasty floating point comparison of the form
+ * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
+ * relied upon to spot anything other than the most obvious input error.
+ *
+ * Returns: 1 if the quaternion is normalised, 0 if not.
+ **/
+int quaternion_valid(struct quaternion q)
+{
+	double qmod;
+
+	qmod = quaternion_modulus(q);
+
+	/* Modulus = 1 to within some tolerance?
+	 * Nasty allowance for floating-point accuracy follows... */
+	if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
+
+	return 0;
+}
+
+
+/**
+ * quat_rot
+ * @q: A vector (in the form of a "struct rvec")
+ * @z: A %quaternion
+ *
+ * Rotates a vector according to a quaternion.
+ *
+ * Returns: A rotated version of @p.
+ **/
+struct rvec quat_rot(struct rvec q, struct quaternion z)
+{
+	struct rvec res;
+	double t01, t02, t03, t11, t12, t13, t22, t23, t33;
+
+	t01 = z.w*z.x;
+	t02 = z.w*z.y;
+	t03 = z.w*z.z;
+	t11 = z.x*z.x;
+	t12 = z.x*z.y;
+	t13 = z.x*z.z;
+	t22 = z.y*z.y;
+	t23 = z.y*z.z;
+	t33 = z.z*z.z;
+
+	res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
+	            + (2.0 * (t12 + t03)) * q.v
+	            + (2.0 * (t13 - t02)) * q.w;
+
+	res.v =       (2.0 * (t12 - t03)) * q.u
+	      + (1.0 - 2.0 * (t11 + t33)) * q.v
+	            + (2.0 * (t01 + t23)) * q.w;
+
+	res.w =       (2.0 * (t02 + t13)) * q.u
+	            + (2.0 * (t23 - t01)) * q.v
+	      + (1.0 - 2.0 * (t11 + t22)) * q.w;
+
+	return res;
+}
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index 4955f875..a759ff15 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -61,44 +61,10 @@
 #define THOMSON_LENGTH (2.81794e-15)
 
 
-/* ------------------------------ Quaternions ------------------------------- */
-
-/**
- * quaternion:
- *
- * <programlisting>
- * struct quaternion
- * {
- *    double w
- *    double x
- *    double y
- *    double z
- * };
- * </programlisting>
- *
- * A structure representing a quaternion.
- *
- **/
-struct quaternion;
-
-struct quaternion {
-	double w;
-	double x;
-	double y;
-	double z;
-};
-
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-extern struct quaternion normalise_quaternion(struct quaternion q);
-extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(gsl_rng *rng);
-extern int quaternion_valid(struct quaternion q);
-extern struct rvec quat_rot(struct rvec q, struct quaternion z);
-
-
 /* --------------------------- Useful functions ----------------------------- */
 
 extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v);
@@ -126,17 +92,6 @@ extern double flat_noise(gsl_rng *rng, double expected, double width);
 extern double gaussian_noise(gsl_rng *rng, double expected, double stddev);
 extern int poisson_noise(gsl_rng *rng, double expected);
 
-/* Keep these ones inline, to avoid function call overhead */
-static inline struct quaternion invalid_quaternion(void)
-{
-	struct quaternion quat;
-	quat.w = 0.0;
-	quat.x = 0.0;
-	quat.y = 0.0;
-	quat.z = 0.0;
-	return quat;
-}
-
 static inline double distance(double x1, double y1, double x2, double y2)
 {
 	return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
@@ -275,6 +230,47 @@ extern char *safe_basename(const char *in);
 #define unlikely(x) (x)
 #endif
 
+
+/* ------------------------------ Quaternions ------------------------------- */
+
+/**
+ * quaternion:
+ *  @w: component
+ *  @x: component
+ *  @y: component
+ *  @z: component
+ *
+ * A structure representing a quaternion.
+ *
+ **/
+struct quaternion;
+
+struct quaternion {
+	double w;
+	double x;
+	double y;
+	double z;
+};
+
+extern struct quaternion normalise_quaternion(struct quaternion q);
+extern double quaternion_modulus(struct quaternion q);
+extern struct quaternion random_quaternion(gsl_rng *rng);
+extern int quaternion_valid(struct quaternion q);
+extern struct rvec quat_rot(struct rvec q, struct quaternion z);
+
+/* Keep these ones inline, to avoid function call overhead */
+static inline struct quaternion invalid_quaternion(void)
+{
+	struct quaternion quat;
+	quat.w = 0.0;
+	quat.x = 0.0;
+	quat.y = 0.0;
+	quat.z = 0.0;
+	return quat;
+}
+
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c
index 15826760..8a7daab7 100644
--- a/libcrystfel/src/xds.c
+++ b/libcrystfel/src/xds.c
@@ -3,12 +3,12 @@
  *
  * Invoke xds for crystal autoindexing
  *
- * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2013 Cornelius Gati
  *
  * Authors:
- *   2010-2016 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *        2013 Cornelius Gati <cornelius.gati@cfel.de>
  *
  * This file is part of CrystFEL.
@@ -498,7 +498,7 @@ static int write_inp(struct image *image, struct xds_private *xp)
 }
 
 
-int run_xds(struct image *image, IndexingPrivate *priv)
+int run_xds(struct image *image, void *priv)
 {
 	unsigned int opts;
 	int status;
@@ -619,8 +619,8 @@ int run_xds(struct image *image, IndexingPrivate *priv)
 }
 
 
-IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
-                             struct detector *det, float *ltl)
+void *xds_prepare(IndexingMethod *indm, UnitCell *cell,
+                  struct detector *det, float *ltl)
 {
 	struct xds_private *xp;
 	int need_cell = 0;
@@ -687,11 +687,11 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
 	xp->cell = cell;
 	xp->indm = *indm;
 
-	return (IndexingPrivate *)xp;
+	return xp;
 }
 
 
-void xds_cleanup(IndexingPrivate *pp)
+void xds_cleanup(void *pp)
 {
 	struct xds_private *xp;
 
diff --git a/libcrystfel/src/xds.h b/libcrystfel/src/xds.h
index a0db2054..094d6d2f 100644
--- a/libcrystfel/src/xds.h
+++ b/libcrystfel/src/xds.h
@@ -3,13 +3,13 @@
  *
  * Invoke xds for crystal autoindexing
  *
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  * Copyright © 2013 Cornelius Gati
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
- *        2013 Cornelius Gati <cornelius.gati@cfel.de>
+ *   2010-2013,2017 Thomas White <taw@physics.org>
+ *        2013      Cornelius Gati <cornelius.gati@cfel.de>
  *
  * This file is part of CrystFEL.
  *
@@ -42,12 +42,12 @@
 extern "C" {
 #endif
 
-extern int run_xds(struct image *image, IndexingPrivate *ipriv);
+extern int run_xds(struct image *image, void *ipriv);
 
-extern IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
-                                    struct detector *det, float *ltl);
+extern void *xds_prepare(IndexingMethod *indm, UnitCell *cell,
+                         struct detector *det, float *ltl);
 
-extern void xds_cleanup(IndexingPrivate *pp);
+extern void xds_cleanup(void *pp);
 
 #ifdef __cplusplus
 }