Merge branch 'master' into taketwo

author: Thomas White <taw@physics.org> 2017-07-06 17:22:11 +0200
committer: Thomas White <taw@physics.org> 2017-07-06 17:24:59 +0200
commit: 5292f57d4434c7267e8d945871513d742ff42427 (patch)
tree: d460aa5cef5a501516876850ef243cfc27313d5a
parent: 48d4a6b8e82cce81222ec58fdfb488ed79ce0bcf (diff)
parent: dc3395900fc3ce0d3961757628ff83ad6456be19 (diff)
76 files changed, 3837 insertions, 902 deletions
diff --git a/AUTHORS b/AUTHORS
index 45af8ddf..0299e408 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,3 +1,23 @@
+Funding acknowledgements
+------------------------
+
+Development of CrystFEL is primarily funded by the Helmholtz Association via
+programme-oriented funds.
+
+Additional funding for CrystFEL is provided by "X-Probe", a project of the
+European Union's 2020 Research and Innovation Program Under the Marie
+Skłodowska-Curie grant agreement 637295 (2015-2018).
+
+Additional funding for CrystFEL is provided by the BMBF German-Russian
+Cooperation "SyncFELMed", grant 05K14CHA (2014-2017).
+
+Past funding for CrystFEL has been received from BioStruct-X, a project funded
+by the Seventh Framework Programme (FP7) of the European Commission.
+
+
+Contributors
+------------
+
 * Thomas White <taw@physics.org>
   Lead author, architecture, "vision" and so on.
 
@@ -44,6 +64,7 @@
 
 * Wolfgang Brehm <wolfgang.brehm@gmail.com>
   "Detwinning" algorithm (ambigator)
+  Fixes for d1sig and d2sig
 
 * Valerio Mariani <<alerio.mariani@desy.de>
   Packaging and general bug fixing
@@ -72,3 +93,13 @@
 
 * Thomas Grant <tgrant@hwi.buffalo.edu>
   Original version of scripts/peakogram-stream
+
+* Steve Aplin <steve.aplin@desy.de>
+  Original version of scripts/turbo-index-slurm
+
+* Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+  geoptimiser
+  Peakfinder8 improvements
+
+* Helen Ginn  <helen@strubi.ox.ac.uk>
+  TakeTwo indexing algorithm
diff --git a/ChangeLog b/ChangeLog
index 86788368..f45c1fad 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,48 @@
+Changes in this development version
+-----------------------------------
+
+- indexamajig --peaks=peakfinder8 was added (Valerio Mariani, Oleksandr Yefanov)
+- CBF files are now supported as input.
+- TakeTwo indexing was incorporated (Helen Ginn)
+- indexamajig: Hung worker processes will now be detected and killed.
+- Peak locations from HDF5 or CXI files are now offset by 0.5,0.5.
+- Detector panels no longer need to be perpendicular to the X-ray beam.
+- Detector "rail" direction ("camera length" axis) no longer needs to be along
+  the X-ray beam direction.
+- Lattice parameters are now checked after prediction refinement as well as before.
+- Multi-event HDF5 files can now contain peak lists.
+- The peak list location can now be given in geometry file (instead of --hdf5-peaks)
+- The number of detector units per photon (rather than per eV) can now be
+  specified in the geometry file.
+- indexamajig: Small changes to how peaks from HDF5 files are checked.
+- compare_hkl --fom=d1sig and d2sig were fixed (Wolfgang Brehm).
+- compare_hkl --min-measurements was added.
+- New polarisation correction.
+- Reflection data files now contain audit information (CrystFEL version number
+  and command line parameters).
+- Improvements to enumeration of events in multi-event files
+- scripts/detector-shift: Show a heat map, handle different panel resolutions
+  correctly
+- Add scripts/move-entire-detector, scripts/split-by-mask,
+  scripts/peakogram-stream and scripts/sum-peaks.
+- scripts/turbo-index was renamed to turbo-index-lsf, and turbo-index-slurm
+  was added (Steve Aplin).
+- An example geometry file for Eiger was added.
+- Multiple bug fixes for asdf indexing
+- cell_explorer: Fixed invisible Greek letters on some systems
+- Mask paths with fewer placeholders than data paths (e.g. some SACLA files)
+  are now handled properly.
+- cell_explorer: Added "Save cell" function.
+- pattern_sim: Multiple bug fixes.
+- geoptimiser: Multiple bug fixes.
+- partialator: Fixes for some edge case bugs.
+- indexamajig: Fix files left open when using XDS indexing.
+- {check,compare}_hkl: Add warnings for potentially silly option choices.
+- Remove "data slab" from detector data handling (simplifies internal data
+  structures a lot).
+- Python scripts were made compatible with Python 3.
+
+
 CrystFEL version 0.6.2, 21st March 2016
 ---------------------------------------
 
diff --git a/Makefile.am b/Makefile.am
index eb10db9c..fda28fbc 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -159,6 +159,7 @@ crystfeldoc_DATA = doc/twin-calculator.pdf doc/examples/lcls-dec.geom \
                    doc/examples/cspad-cxiformat.geom \
                    doc/examples/pilatus.geom \
                    doc/examples/cell-example.cell \
+                   doc/examples/Eiger16M-binning2-nativefiles.geom \
                    doc/hitrate.html doc/hitrate.png
 
 EXTRA_DIST += $(crystfeldoc_DATA) doc/twin-calculator.odt \
@@ -180,9 +181,10 @@ script_DATA = scripts/alternate-stream scripts/cell-please \
               scripts/find-pairs scripts/plot-cc-and-scale.R \
               scripts/ave-resolution scripts/crystal-frame-number \
               scripts/plot-radius-resolution \
-              scripts/detector-shift scripts/turbo-index \
+              scripts/detector-shift scripts/turbo-index-lsf \
               scripts/gaincal-to-saturation-map scripts/move-entire-detector \
-              scripts/split-by-mask
+              scripts/split-by-mask scripts/turbo-index-slurm \
+              scripts/sum-peaks
 
 EXTRA_DIST += $(script_DATA)
 
diff --git a/README b/README
index cd911d04..219cca15 100644
--- a/README
+++ b/README
@@ -1,7 +1,7 @@
 CrystFEL - Crystallography with a FEL
 -------------------------------------
 
-Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
                       a research centre of the Helmholtz Association.
 
 Authors:
@@ -26,6 +26,9 @@ Authors:
   Keitaro Yamashita <k.yamashita@spring8.or.jp>
   Marmoru Suzuki <mamoru.suzuki@protein.osaka-u.ac.jp>
   Thomas Grant <tgrant@hwi.buffalo.edu>
+  Steve Aplin <steve.aplin@desy.de>
+  Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+  Helen Ginn  <helen@strubi.ox.ac.uk>
 
 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
@@ -150,11 +153,6 @@ then your version of FFTW3 is not compiled in a suitable way.  You'll need to
 install it again (from source) adding "--enable-shared" to its ./configure
 command line.
 
-If you get an error about "cairo_surface_show_page" in src/scaling-report.c,
-simply edit the file and comment that line out.  Partialator will then not be
-able to produce a useful scaling report, but the rest of CrystFEL will work
-properly.
-
 If you are installing from Git, the following extra things apply:
 
 - You must run "./autogen.sh" to generate "configure" and other files before
@@ -174,3 +172,19 @@ If you are installing from Git, the following extra things apply:
 - You will not be able to use the "--enable-gtk-doc" option to configure unless
   you have at least version 1.9 of gtk-doc installed.
 
+
+Funding acknowledgements
+------------------------
+
+Development of CrystFEL is primarily funded by the Helmholtz Association via
+programme-oriented funds.
+
+Additional funding for CrystFEL is provided by "X-Probe", a project of the
+European Union's 2020 Research and Innovation Program Under the Marie
+Skłodowska-Curie grant agreement 637295 (2015-2018).
+
+Additional funding for CrystFEL is provided by the BMBF German-Russian
+Cooperation "SyncFELMed", grant 05K14CHA (2014-2017).
+
+Past funding for CrystFEL has been received from BioStruct-X, a project funded
+by the Seventh Framework Programme (FP7) of the European Commission.
diff --git a/configure.ac b/configure.ac
index 8044a145..2a5f40cf 100644
--- a/configure.ac
+++ b/configure.ac
@@ -8,6 +8,7 @@ AM_INIT_AUTOMAKE([subdir-objects])
 AC_PROG_CC
 gl_EARLY
 AM_PROG_CC_C_O
+AC_PROG_CC_C99
 AC_PROG_AWK
 AC_PROG_INSTALL
 AC_PROG_LN_S
@@ -36,6 +37,39 @@ AC_ARG_WITH(gsl,
  GSL_LIBS="-L$withval/lib -lgsl -lgslcblas -lm"],
 [PKG_CHECK_MODULES([GSL], [gsl])])
 
+AC_ARG_ENABLE(cbf, AS_HELP_STRING([--disable-cbf], [Disable CBF file support]))
+AC_ARG_WITH(cbflib_dir, AS_HELP_STRING([--with-cbflib],
+            [Specify location of CBFlib headers and libraries]))
+AC_MSG_CHECKING([whether to use CBFlib])
+AS_IF([test "x$enable_cbf" != "xno"],
+[
+   AC_MSG_RESULT([yes])
+   AS_IF([test "x$with_cbflib_dir" != "x"],
+   [
+      dnl NB not ${with_cbflib_dir}/include/cbflib, because cbflib installs
+      dnl its own HDF5 headers which we do not want in the include path
+      CBF_CFLAGS="-I${with_cbflib_dir}/include"
+      dnl Fortunately no libhdf5.so in this folder
+      CBF_LIBS="-lcbf -L${with_cbflib_dir}/lib"
+   ], [
+      AC_CHECK_LIB([cbf], [cbf_make_handle], [
+         CBF_LIBS="-lcbf"
+      ])
+   ])
+   SAVED_CPPFLAGS=$CPPFLAGS
+   CPPFLAGS="$CPPFLAGS $CBF_CFLAGS"
+   AC_CHECK_HEADERS([cbflib/cbf.h cbf/cbf.h],
+   [
+      have_cbflib=true
+   ])
+   CPPFLAGS=$SAVE_CPPFLAGS
+], [
+   AC_MSG_RESULT([no])
+])
+AS_IF([test x$have_cbflib = xtrue],
+[
+   AC_DEFINE([HAVE_CBFLIB], [1], [Define to 1 if CBFlib is available])
+])
 
 AC_ARG_WITH(opencl, AS_HELP_STRING([--with-opencl], [Use OpenCL]))
 AC_ARG_WITH(opencl_dir, AS_HELP_STRING([--with-opencl-dir],
@@ -324,7 +358,7 @@ MAIN_CFLAGS="$MAIN_CFLAGS $GDK_pixbuf_2_CFLAGS $FFTW_CFLAGS $Pango_CFLAGS"
 MAIN_CFLAGS="$MAIN_CFLAGS $PangoCairo_CFLAGS"
 AC_SUBST([MAIN_CFLAGS])
 
-LIBCRYSTFEL_CFLAGS="$CFLAGS $HDF5_CFLAGS $GSL_CFLAGS $FFTW_CFLAGS"
+LIBCRYSTFEL_CFLAGS="$CFLAGS $HDF5_CFLAGS $GSL_CFLAGS $FFTW_CFLAGS $CBF_CFLAGS"
 AC_SUBST([LIBCRYSTFEL_CFLAGS])
 
 MAIN_LIBS="$LIBS $HDF5_LIBS -lm -lz $GSL_LIBS $GTK_LIBS $OPENCL_LIBS -pthread"
@@ -334,7 +368,7 @@ MAIN_LIBS="$MAIN_LIBS $Pango_LIBS $PangoCairo_LIBS $LDFLAGS"
 AC_SUBST([MAIN_LIBS])
 
 LIBCRYSTFEL_LIBS="$LIBS $HDF5_LIBS $GSL_LIBS $FFTW_LIBS $CURSES_LIB $LDFLAGS"
-LIBCRYSTFEL_LIBS="$LIBCRYSTFEL_LIBS -pthread"
+LIBCRYSTFEL_LIBS="$LIBCRYSTFEL_LIBS -pthread $CBF_LIBS"
 AC_SUBST([LIBCRYSTFEL_LIBS])
 
 GTK_DOC_CHECK([1.9],[--flavour no-tmpl])
diff --git a/doc/examples/Eiger16M-binning2-nativefiles.geom b/doc/examples/Eiger16M-binning2-nativefiles.geom
new file mode 100644
index 00000000..0a53a74e
--- /dev/null
+++ b/doc/examples/Eiger16M-binning2-nativefiles.geom
@@ -0,0 +1,55 @@
+; Example geometry file for Eiger 16M detector, using its native file format
+; and binning 2.
+
+; Camera length (in m) and photon energy (eV)
+clen = 0.1
+photon_energy = 22000
+
+; adu_per_photon needs a relatively recent CrystFEL version.  If your version is
+; older, change it to adu_per_eV and set it to one over the photon energy in eV
+adu_per_photon = 1
+res = 13333.3   ; 75 micron pixel size
+
+; These lines describe the data layout for the Eiger native multi-event files
+dim0 = %
+dim1 = ss
+dim2 = fs
+data = /entry/data/data
+
+; Mask out strips between panels
+bad_v0/min_fs = 1030
+bad_v0/min_ss = 0
+bad_v0/max_fs = 1039
+bad_v0/max_ss = 2166
+
+bad_h0/min_fs = 0
+bad_h0/min_ss = 514
+bad_h0/max_fs = 2069
+bad_h0/max_ss = 550
+
+bad_h1/min_fs = 0
+bad_h1/min_ss = 1065
+bad_h1/max_fs = 2069
+bad_h1/max_ss = 1101
+
+bad_h2/min_fs = 0
+bad_h2/min_ss = 1616
+bad_h2/max_fs = 2069
+bad_h2/max_ss = 1652
+
+; Uncomment these lines if you have a separate bad pixel map (recommended!)
+;mask_file = eiger-badmap.h5
+;mask = /data/data
+;mask_good = 0x0
+;mask_bad = 0x1
+
+; corner_{x,y} set the position of the corner of the detector (in pixels)
+; relative to the beam
+panel0/min_fs = 0
+panel0/min_ss = 0
+panel0/max_fs = 2069
+panel0/max_ss = 2166
+panel0/corner_x = -1000.0
+panel0/corner_y = -1000.0
+panel0/fs = x
+panel0/ss = y
diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1
index 8e464d12..30fcd93e 100644
--- a/doc/man/compare_hkl.1
+++ b/doc/man/compare_hkl.1
@@ -72,9 +72,9 @@ Note that this figure of merit compares measurements within one data set, so it
 The ratio of Rano to Rsplit, as defined above.
 .IP "\fBd1sig\fR and \fBd2sig\fR"
 .PD
-The fraction of differences between intensities which are within one times (for \fBd1sig\fR) and two times (for \fBd2sig\fR) the mean of the corresponding sigma(I) values.
+The fraction of differences between intensities which are within one times (for \fBd1sig\fR) and two times (for \fBd2sig\fR) the combination of the corresponding sigma(I) values.
 .PP
-I1 and I2 are the intensities of the same reflection in both reflection lists.  The scale factor, k, is given by sum(I1*i2) / sum(I2^2), unless you use \fB-u\fR.
+I1 and I2 are the intensities of the same reflection in both reflection lists.  The two sets of reflections will be put on a common scale (linear and Debye-Waller terms) unless you use \fB-u\fR.
 .RE
 
 .PD 0
diff --git a/doc/man/crystfel.7 b/doc/man/crystfel.7
index bc545069..b7063b21 100644
--- a/doc/man/crystfel.7
+++ b/doc/man/crystfel.7
@@ -1,7 +1,7 @@
 .\"
 .\" CrystFEL main man page
 .\"
-.\" Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+.\" Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
 .\"                       a research centre of the Helmholtz Association.
 .\"
 .\" Part of CrystFEL - crystallography with a FEL
@@ -125,6 +125,15 @@ to this list if you experience difficulty with "-y" at any time.
 .IP Cubic
 \fB23\fR, \fBm-3\fR, \fB432\fR, \fB-43m\fR, \fBm-3m\fR.
 
+.SH FUNDING ACKNOWLEDGEMENTS
+Development of CrystFEL is primarily funded by the Helmholtz Association via programme-oriented funds.
+
+Additional funding for CrystFEL is provided by "X-Probe", a project of the European Union's 2020 Research and Innovation Program Under the Marie Skłodowska-Curie grant agreement 637295 (2015-2018).
+
+Additional funding for CrystFEL is provided by the BMBF German-Russian Cooperation "SyncFELMed", grant 05K14CHA (2014-2017).
+
+Past funding for CrystFEL has been received from BioStruct-X, a project funded by the Seventh Framework Programme (FP7) of the European Commission.
+
 .SH AUTHOR
 This page was written by Thomas White and Valerio Mariani.
 
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index d7f59f9d..eaf812d6 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -1,7 +1,7 @@
 .\"
 .\" indexamajig man page
 .\"
-.\" Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+.\" Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
 .\"                       a research centre of the Helmholtz Association.
 .\"
 .\" Part of CrystFEL - crystallography with a FEL
@@ -45,8 +45,12 @@ You can control the peak detection on the command line.  Firstly, you can choose
 
 \fB--peaks=cxi\fR works similarly to this, but expects four separate HDF5 datasets beneath \fIpath\fR, \fBnPeaks\fR, \fBpeakXPosRaw\fR, \fBpeakYPosRaw\fR and \fBpeakTotalIntensity\fR.  See the specification for the CXI file format at http://www.cxidb.org/ for more details.
 
+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 HDF5 or CXI files, including Cheetah, may instead consider the peak locations to be pixel indices in the data array.  Therefore, the peak coordinates from \fB--peaks=cxi\fR or \fB--peaks=hdf5\fR will by default have 0.5 added to them.  Use \fB--no-half-pixel-shift\fR if this isn't what you want.
+
 If you use \fB--peaks=zaef\fR, indexamajig will use a simple gradient search after Zaefferer (2000).  You can control the overall threshold and minimum squared gradient for finding a peak using \fB--threshold\fR and \fB--min-gradient\fR.  The threshold has arbitrary units matching the pixel values in the data, and the minimum gradient has the equivalent squared units.  Peaks will be rejected if the 'foot point' is further away from the 'summit' of the peak by more than the inner integration radius (see below).  They will also be rejected if the peak is closer than twice the inner integration radius from another peak.
 
+If you instead use \fB--peaks=peakfinder8\fR, indexamajig will use the "peakfinder8" peak finding algorithm describerd in Barty et al. (2014). Pixels above a radius-dependent intensity threshold are considered as candidate peaks (although the user sets an absolute minimum threshold for candidate peaks). Peaks are then only accepted if their signal to noise level over the local background is sufficiently high. Peaks can include multiple pixels and the user can reject a peak if it includes too many or too few. The distance of a peak from the center of the detector can also be used as a filtering criterion. Note that the peakfinder8 will not report more than 2048 peaks for each panel: any additional peak is ignored.
+
 You can suppress peak detection altogether for a panel in the geometry file by specifying the "no_index" value for the panel as non-zero.
 
 
@@ -229,6 +233,16 @@ these.  The defaults are probably not appropriate for your situation.
 The default is \fB--int-radius=4,5,7\fR.
 
 .PD 0
+.IP \fB--min-peaks=\fIn\fR
+.PD
+Do not try to index frames with fewer than \fIn\fR peaks.  These frames will still be described in the output stream.  To exclude them, use \fB--no-non-hits-in-stream\fR.
+
+.PD 0
+.IP \fB--no-non-hits-in-stream\fR
+.PD
+Completely exclude 'non-hit' frames in the stream.  When this option is given, frames with fewer than the number of peaks given to \fB--min-peaks\fR will not have chunks written to the stream at all.
+
+.PD 0
 .IP \fB--basename\fR
 .PD
 Remove the directory parts of the filenames taken from the input file.  If \fB--prefix\fR or \fB-x\fR is also given, the directory parts of the filename will be removed \fIbefore\fR adding the prefix.
@@ -270,7 +284,7 @@ This option is here for historical purposes only, to disable a correction which
 .PD 0
 .IP \fB--threshold=\fR\fIthres\fR
 .PD
-Set the overall threshold for peak detection using \fB--peaks=zaef\fR to \fIthres\fR, which has the same units as the detector data.  The default is \fB--threshold=800\fR.
+Set the overall threshold for peak detection using \fB--peaks=zaef\fR or \fB--peaks=peakfinder8\fR to \fIthres\fR, which has the same units as the detector data.  The default is \fB--threshold=800\fR.
 
 .PD 0
 .IP \fB--min-gradient=\fR\fIgrad\fR
@@ -280,12 +294,37 @@ Set the square of the gradient threshold for peak detection using \fB--peaks=zae
 .PD 0
 .IP \fB--min-snr=\fR\fIsnr\fR
 .PD
-Set the minimum I/sigma(I) for peak detection when using \fB--peaks=zaef\fR.  The default is \fB--min-snr=5\fR.
+Set the minimum I/sigma(I) for peak detection when using \fB--peaks=zaef\fR or \fB--peaks=peakfinder8\fR.  The default is \fB--min-snr=5\fR.
+
+.PD 0
+.IP \fB--min-pix-count=\fR\fIcnt\fR
+.PD
+Accepts peaks only if they include more than \fR\fIcnt\fR pixels, when using \fB--peaks=peakfinder8\fR.  The default is \fB--min-pix-count=2\fR.
+
+.PD 0
+.IP \fB--max-pix-count=\fR\fIcnt\fR
+.PD
+Accepts peaks only if they include less than \fR\fIcnt\fR pixels, when using \fB--peaks=peakfinder8\fR.  The default is \fB--max-pix-count=200\fR.
+
+.PD 0
+.IP \fB--local-bg-radius=\fR\fIr\fR
+.PD
+Radius (in pixels) used for the estimation of the local background when using \fB--peaks=peakfinder8\fR.  The default is \fB--local-bg-radius=3\fR.
+
+.PD 0
+.IP \fB--min-res=\fR\fIpx\fR
+.PD
+Only accept peaks if they lay at more than \fR\fIpx\fR pixels from the center of the detector when using \fB--peaks=peakfinder8\fR.  The default is \fB--min-res=0\fR.
+
+.PD 0
+.IP \fB--max-res=\fR\fIpx\fR
+.PD
+Only accept peaks if they lay at less than \fR\fIpx\fR pixels from the center of the detector when using \fB--peaks=peakfinder8\fR.  The default is \fB--max-res=1200\fR.
 
 .PD 0
 .IP \fB--copy-hdf5-field=\fR\fIpath\fR
 .PD
-Copy the information from \fIpath\fR in the HDF5 file into the output stream.  The information must be a single scalar value.  This option is sometimes useful to allow data to be separated after indexing according to some condition such the presence of an optical pump pulse.  You can give this option as many times as you need to copy multiple bits of information.
+Copy the information from \fR\fIpath\fR in the HDF5 file into the output stream.  The information must be a single scalar value.  This option is sometimes useful to allow data to be separated after indexing according to some condition such the presence of an optical pump pulse.  You can give this option as many times as you need to copy multiple bits of information.
 
 .PD 0
 .IP "\fB-j\fR \fIn\fR"
@@ -308,6 +347,11 @@ Normally, peaks which contain one or more pixels above max_adu (defined in the d
 When using \fB--peaks=hdf5\fR or \fB--peaks=cxi\fR, the peaks will be put through some of the same checks as if you were using \fB--peaks=zaef\fR.  These checks reject peaks which are too close to panel edges, are saturated (unless you use \fB--use-saturated\fR), have other nearby peaks (closer than twice the inner integration radius, see \fB--int-radius\fR), or have any part in a bad region.  Using this option skips this validation step, and uses the peaks directly.
 
 .PD 0
+.IP \fB--no-half-pixel-shift\fR
+.PD
+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 HDF5 or CXI files, including Cheetah, may instead consider the peak locations to be pixel indices in the data array.  Therefore, the peak coordinates from \fB--peaks=cxi\fR or \fB--peaks=hdf5\fR will by default have 0.5 added to them.  This option \fBdisables\fR this half-pixel offset.
+
+.PD 0
 .IP \fB--check-hdf5-snr\fR
 .PD
 With this option with \fB--peaks=hdf5\fR, the peaks will additionally be checked to see that they satisfy the minimum SNR specified with \fB--min-snr\fR.
diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt
index 4cd922d2..072f8b38 100644
--- a/doc/reference/libcrystfel/CrystFEL-sections.txt
+++ b/doc/reference/libcrystfel/CrystFEL-sections.txt
@@ -8,6 +8,8 @@ reflist_new
 reflist_free
 reflection_new
 reflection_free
+reflist_add_notes
+reflist_get_notes
 <SUBSECTION>
 add_refl
 add_refl_to_list
@@ -18,7 +20,6 @@ next_refl
 find_refl
 next_found_refl
 <SUBSECTION>
-get_excitation_error
 get_detector_pos
 get_panel
 get_partiality
@@ -77,6 +78,7 @@ copy_reflist
 find_equiv_in_list
 resolution_limits
 max_intensity
+reflist_add_command_and_version
 </SECTION>
 
 <SECTION>
@@ -151,32 +153,32 @@ solve_svd
 AssplodeFlag
 C_VACUO
 ELECTRON_CHARGE
-ERROR
-J_to_eV
 PLANCK
-STATUS
 THOMSON_LENGTH
 assplode
-biggest
-check_prefix
 chomp
-deg2rad
-eV_to_J
 gaussian_noise
-safe_basename
 progress_bar
-rad2deg
-is_odd
 poisson_noise
 notrail
+random_flat
+flat_noise
+show_matrix
+STATUS
+ERROR
+J_to_eV
+eV_to_J
+biggest
 smallest
+check_prefix
+safe_basename
+deg2rad
+rad2deg
+is_odd
 ph_en_to_lambda
 ph_lambda_to_en
 ph_eV_to_lambda
 ph_lambda_to_eV
-random_flat
-flat_noise
-show_matrix
 likely
 unlikely
 UNUSED
@@ -218,6 +220,20 @@ ImageFeatureList
 SpectrumType
 sample
 <SUBSECTION>
+imagefile
+imagefile_close
+imagefile_copy_fields
+imagefile_field_list
+imagefile_get_hdfile
+imagefile_get_type
+imagefile_open
+imagefile_read
+imagefile_read_simple
+<SUBSECTION>
+new_imagefile_field_list
+add_imagefile_field
+free_imagefile_field_list
+<SUBSECTION>
 image_add_feature
 image_feature_closest
 image_reflection_closest
@@ -268,10 +284,11 @@ INDEXING_DEFAULTS_ASDF
 INDEXING_DEFAULTS_FELIX
 INDEXING_METHOD_MASK
 INDEXING_CONTROL_FLAGS
-build_indexer_list
+setup_indexing
 cleanup_indexing
-prepare_indexing
+get_indm_from_string
 index_pattern
+index_pattern_2
 indexer_str
 dirax_prepare
 run_dirax
@@ -314,6 +331,7 @@ describe_symmetry
 symmetry_name
 set_symmetry_name
 get_matrix_name
+pointgroup_warning
 is_centrosymmetric
 is_centric
 </SECTION>
@@ -351,19 +369,17 @@ rigid_group
 rg_collection
 <SUBSECTION>
 copy_geom
-fill_in_values
+fill_in_adu
 free_detector_geometry
 get_detector_geometry
 get_detector_geometry_2
 write_detector_geometry
 write_detector_geometry_2
-find_panel
+panel_number
 find_panel_by_name
-find_panel_number
 simple_geometry
 record_image
 get_pixel_extents
-get_q
 get_q_for_panel
 get_tt
 smallest_q
@@ -378,6 +394,7 @@ rigid_group_is_in_collection
 single_panel_data_source
 find_rigid_group_collection_by_name
 detector_has_clen_references
+adjust_centering_for_rail
 </SECTION>
 
 <SECTION>
@@ -398,7 +415,6 @@ find_event
 get_event_string
 get_event_from_event_string
 event_path_placeholder_subst
-partial_event_substitution
 retrieve_full_path
 initialize_filename_plus_event
 free_filename_plus_event
@@ -439,7 +455,9 @@ add_copy_hdf5_field
 new_copy_hdf5_field_list
 free_copy_hdf5_field_list
 get_peaks
+get_peaks_2
 get_peaks_cxi
+get_peaks_cxi_2
 hdfile_is_scalar
 check_path_existence
 fill_event_list
@@ -464,6 +482,7 @@ crystal_get_resolution_limit
 crystal_get_user_flag
 crystal_get_num_implausible_reflections
 crystal_get_notes
+crystal_get_det_shift
 crystal_set_cell
 crystal_set_image
 crystal_set_mosaicity
@@ -477,6 +496,7 @@ crystal_set_user_flag
 crystal_set_num_implausible_reflections
 crystal_set_notes
 crystal_add_notes
+crystal_set_det_shift
 </SECTION>
 
 <SECTION>
@@ -498,7 +518,9 @@ y_gradient
 <FILE>peaks</FILE>
 peak_sanity_check
 search_peaks
+search_peaks_peakfinder8
 make_BgMask
+peakfinder8
 validate_peaks
 sort_peaks
 </SECTION>
@@ -555,6 +577,7 @@ close_stream
 read_chunk
 read_chunk_2
 write_chunk
+write_chunk_2
 rewind_stream
 is_stream
 write_command
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
 }
diff --git a/scripts/Rsplit_surface.py b/scripts/Rsplit_surface.py
index ca68a92e..63a72cd8 100644
--- a/scripts/Rsplit_surface.py
+++ b/scripts/Rsplit_surface.py
@@ -1,13 +1,5 @@
 #!/usr/bin/python
-# coding=utf-8
-
-from __future__ import division
-from array import array
-from numpy import *
-
-import matplotlib
-
-# Rsplit_surface.py
+# -*- coding: utf-8 -*-
 #
 # Plot Rsplit as a contour map
 #
@@ -32,6 +24,12 @@ import matplotlib
 # along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
 
 
+from __future__ import division
+from array import array
+from numpy import *
+
+import matplotlib
+
 """there could be problems with dependencies, in this case, """
 #matplotlib.use('PS')
 
diff --git a/scripts/add-beam-params b/scripts/add-beam-params
index 23c9caaf..e4e38f70 100755
--- a/scripts/add-beam-params
+++ b/scripts/add-beam-params
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Add beam parameters to stream
 #
-# Copyright (c) 2014 Deutsches Elektronen-Synchrotron DESY,
-#                    a research centre of the Helmholtz Association.
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2014 Thomas White <taw@physics.org>
+#    2014,2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -66,7 +66,7 @@ def count_crystals(f, g, start_after, stop_after):
                 in_crystal = 0
                 need_end_chunk = 1
 
-    print "Wrote %i crystals to %s" % (n_crystals_written, opt.ofn)
+    print("Wrote {} crystals to {}".format(n_crystals_written, opt.ofn))
 
 def count_chunks(f, g, start_after, stop_after):
     n_chunks_seen = 0
@@ -95,7 +95,7 @@ def count_chunks(f, g, start_after, stop_after):
                 if n_chunks_written == stop_after:
                     break
 
-    print "Wrote %i chunks to %s" % (n_chunks_written, opt.ofn)
+    print("Wrote {} chunks to {}".format(n_chunks_written, opt.ofn))
 
 op = optparse.OptionParser()
 op.add_option('', '--input', action='store', type='string', dest='ifn',
@@ -111,7 +111,7 @@ op.add_option('', '--chunks', action='store_true', dest='chunks',
 opt,arg = op.parse_args(sys.argv)
 
 if not (opt.ifn and opt.ofn):
-    print "You need at least --input and --output"
+    print("You need at least --input and --output")
     exit(1)
 
 f = open(opt.ifn, 'r')
diff --git a/scripts/ave-resolution b/scripts/ave-resolution
index a69c456e..69c21195 100755
--- a/scripts/ave-resolution
+++ b/scripts/ave-resolution
@@ -1,14 +1,13 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
 #
 # Find mean diffracting resolution
 #
-# Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY,
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
 #                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2014-2015 Thomas White <taw@physics.org>
+#    2014-2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -24,16 +23,16 @@ while True:
         break
     if fline.find("diffraction_resolution_limit") != -1:
         res = float(fline.split('= ')[1].split(' ')[0].rstrip("\r\n"))
-	a.append(res)
+        a.append(res)
         continue
 
 f.close()
 
 b = numpy.array(a)
-print " Mean: %.2f nm^-1 = %.2f A" % (numpy.mean(b),10.0/numpy.mean(b))
-print " Best: %.2f nm^-1 = %.2f A" % (numpy.max(b),10.0/numpy.max(b))
-print "Worst: %.2f nm^-1 = %.2f A" % (numpy.min(b),10.0/numpy.min(b))
-print "Std deviation: %.2f nm^-1" % (numpy.std(b))
+print(" Mean: {:.2} nm^-1 = {:.2} A".format(numpy.mean(b),10.0/numpy.mean(b)))
+print(" Best: {:.2} nm^-1 = {:.2} A".format(numpy.max(b),10.0/numpy.max(b)))
+print("Worst: {:.2} nm^-1 = {:.2} A".format(numpy.min(b),10.0/numpy.min(b)))
+print("Std deviation: {:.2} nm^-1".format(numpy.std(b)))
 
 plt.hist(a,bins=30)
 plt.title('Resolution based on indexing results')
diff --git a/scripts/clean-stream.py b/scripts/clean-stream.py
index 1e906e54..57c0ed93 100644
--- a/scripts/clean-stream.py
+++ b/scripts/clean-stream.py
@@ -1,6 +1,6 @@
 #!/usr/bin/python
 # coding=utf-8
-
+#
 # clean-stream.py
 #
 # Remove non-indexed frames from a stream
@@ -86,5 +86,5 @@ while  (line != ''):
 
 	line = infile_1.readline()
 
-print '%d suited of %d patterns have been extracted and saved as %s' % (num_suited, n_patt, sys.argv[2])
+print('%d suited of %d patterns have been extracted and saved as %s' % (num_suited, n_patt, sys.argv[2]))
 Nfile.write('%d' % num_suited)
diff --git a/scripts/crystal-frame-number b/scripts/crystal-frame-number
index b3ed3f1f..883abb81 100755
--- a/scripts/crystal-frame-number
+++ b/scripts/crystal-frame-number
@@ -1,6 +1,13 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
+#
+# Show sequence numbers of crystals and frames
+#
+# Copyright © 2015-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+#
+# Author:
+#    2015-2017 Thomas White <taw@physics.org>
 
 import sys
 
@@ -16,9 +23,9 @@ while True:
     if fline.find("Image filename") != -1:
         frame_number += 1
         fn = fline.split(': ')[1].split(' ')[0].rstrip("\r\n")
-	print 'Frame %i: %s' % (frame_number, fn)
+        print('Frame {}: {}'.format(frame_number, fn))
     if fline.find("diffraction_resolution_limit") != -1:
         crystal_number += 1
-	print '  Crystal %i: %s' % (crystal_number, fline.rstrip("\r\n"))
+        print('  Crystal {}: {}'.format(crystal_number, fline.rstrip("\r\n")))
 
 f.close()
diff --git a/scripts/detector-shift b/scripts/detector-shift
index a348e3c9..24cbb504 100755
--- a/scripts/detector-shift
+++ b/scripts/detector-shift
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Determine mean detector shift based on prediction refinement results
 #
-# Copyright (c) 2015-2016 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.
 #
 # Author:
-#    2015-2016 Thomas White <taw@physics.org>
+#    2015-2017 Thomas White <taw@physics.org>
 #    2016      Marmoru Suzuki <mamoru.suzuki@protein.osaka-u.ac.jp>
 #
 
@@ -54,13 +54,13 @@ f.close()
 
 mean_x = sum(x_shifts) / len(x_shifts)
 mean_y = sum(y_shifts) / len(y_shifts)
-print 'Mean shifts: dx = %.2f mm,  dy = %.2f mm' % (mean_x,mean_y)
+print('Mean shifts: dx = {:.2} mm,  dy = {:.2} mm'.format(mean_x,mean_y))
 
 # Apply shifts to geometry
 if have_geom:
 
     out = os.path.splitext(geom)[0]+'-predrefine.geom'
-    print 'Applying corrections to %s, output filename %s' % (geom,out)
+    print('Applying corrections to {}, output filename {}'.format(geom,out))
     g = open(geom, 'r')
     h = open(out, 'w')
     panel_resolutions = {}
@@ -99,7 +99,7 @@ if have_geom:
                 res = panel_resolutions[panel]
             else:
                 res = default_res
-                print 'Using default resolution (%f px/m) for panel %s' % (res, panel)
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
             h.write('%s/corner_x = %f\n' % (panel,panel_cnx+(mean_x*res*1e-3)))
             continue
 
@@ -111,7 +111,7 @@ if have_geom:
                 res = panel_resolutions[panel]
             else:
                 res = default_res
-                print 'Using default resolution (%f px/m) for panel %s' % (res, panel)
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
             h.write('%s/corner_y = %f\n' % (panel,panel_cny+(mean_y*res*1e-3)))
             continue
 
diff --git a/scripts/find-pairs b/scripts/find-pairs
index 4629767d..8c3b69d0 100755
--- a/scripts/find-pairs
+++ b/scripts/find-pairs
@@ -1,14 +1,13 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
 #
 # Find pairs of observations from the same pattern
 #
-# Copyright © 2014 Deutsches Elektronen-Synchrotron DESY,
-#                  a research centre of the Helmholtz Association.
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2014 Thomas White <taw@physics.org>
+#    2014-2017 Thomas White <taw@physics.org>
 #
 
 import re as regexp
@@ -41,33 +40,33 @@ def find_405(f):
 
 f = open(fn1, 'r')
 flines = find_405(f)
-print len(flines),"measurements in",fn1
+print("{} measurements in {}".format(len(flines),fn1))
 
 g = open(fn2, 'r')
 glines = find_405(g)
-print len(glines),"measurements in",fn2
+print("{} measurements in {}".format(len(glines),fn2))
 
-print "\nThe common measurements:\n"
+print("\nThe common measurements:\n")
 for fn in flines.keys():
     if fn in glines:
-        print fn
-        print flines[fn].rstrip("\r\n")
-        print glines[fn]
-        print
+        print(fn)
+        print(flines[fn].rstrip("\r\n"))
+        print(glines[fn])
+        print()
         del flines[fn]
         del glines[fn]
 
-print "\nThe measurements only in",fn1,":\n"
+print("\nThe measurements only in {}:\n".format(fn1))
 
 for fn in flines.keys():
-    print fn
-    print flines[fn].rstrip("\r\n")
-    print
+    print(fn)
+    print(flines[fn].rstrip("\r\n"))
+    print()
 
-print "\nThe measurements only in",fn2,":\n"
+print("\nThe measurements only in {}:\n".format(fn2))
 
 for fn in glines.keys():
-    print fn
-    print glines[fn].rstrip("\r\n")
-    print
+    print(fn)
+    print(glines[fn].rstrip("\r\n"))
+    print("")
 
diff --git a/scripts/gaincal-to-saturation-map b/scripts/gaincal-to-saturation-map
index 8a803526..eac8883c 100755
--- a/scripts/gaincal-to-saturation-map
+++ b/scripts/gaincal-to-saturation-map
@@ -1,4 +1,13 @@
 #!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Convert gain map to saturation map
+#
+# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+#
+# Author:
+#    2016-2017 Thomas White <taw@physics.org>
 
 import numpy as np
 import h5py
diff --git a/scripts/move-entire-detector b/scripts/move-entire-detector
index 79ce5a5c..3d06c408 100755
--- a/scripts/move-entire-detector
+++ b/scripts/move-entire-detector
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Shift the detector by a given amount
 #
-# Copyright (c) 2016 Deutsches Elektronen-Synchrotron DESY,
-#                    a research centre of the Helmholtz Association.
+# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2016 Thomas White <taw@physics.org>
+#    2016-2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -58,8 +58,8 @@ if args.beam:
 else:
     bm = 'detector'
 
-print 'Input filename %s\nOutput filename %s' % (args.ifn,out)
-print 'Shifting %s by %f,%f %s' % (bm, args.xshift, args.yshift, units)
+print('Input filename {}\nOutput filename {}'.format(args.ifn,out))
+print('Shifting {} by {},{} {}'.format(bm, args.xshift, args.yshift, units))
 
 if args.beam:
     args.xshift = -args.xshift
@@ -104,7 +104,7 @@ while True:
         else:
             res = default_res
             if not args.px:
-                print 'Using default resolution (%f px/m) for panel %s' % (res, panel)
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
         if args.px:
             h.write('%s/corner_x = %f\n' % (panel,panel_cnx+args.xshift))
         else:
@@ -120,7 +120,7 @@ while True:
         else:
             res = default_res
             if not args.px:
-                print 'Using default resolution (%f px/m) for panel %s' % (res, panel)
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
         if args.px:
             h.write('%s/corner_y = %f\n' % (panel,panel_cny+args.yshift))
         else:
diff --git a/scripts/peak-intensity b/scripts/peak-intensity
index 53610b73..65e066cf 100755
--- a/scripts/peak-intensity
+++ b/scripts/peak-intensity
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Quantify peak intensities
 #
-# Copyright (c) 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.
 #
 # Author:
-#    2015 Thomas White <taw@physics.org>
+#    2015-2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -40,8 +40,8 @@ while True:
 
 f.close()
 
-print '%i patterns, %i peaks' % (n_patt,n_peaks)
-print 'Mean %.2f peaks per pattern' % (n_peaks/n_patt)
-print 'Mean %.2f ADU per peak' % (total_intens/n_peaks)
-print 'Mean %.2f ADU total per pattern' % (total_intens/n_patt)
+print('{} patterns, %i peaks'.format(n_patt,n_peaks))
+print('Mean {:.2f} peaks per pattern'.format(n_peaks/float(n_patt)))
+print('Mean {:.2f} ADU per peak'.format(total_intens/float(n_peaks)))
+print('Mean {:.2f} ADU total per pattern'.format(total_intens/n_patt))
 
diff --git a/scripts/peakogram-stream b/scripts/peakogram-stream
index 209e7e95..1f67b39c 100755
--- a/scripts/peakogram-stream
+++ b/scripts/peakogram-stream
@@ -1,16 +1,14 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
-# peakogram-stream
 #
 # Check a stream for saturation
 #
-# Copyright © 2016      Deutsches Elektronen-Synchrotron DESY,
+# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
 #                       a research centre of the Helmholtz Association.
 # Copyright © 2016      The Research Foundation for SUNY
 #
 # Authors:
-#   2016      Thomas White <taw@physics.org>
+#   2016-2017 Thomas White <taw@physics.org>
 #   2014-2016 Thomas Grant <tgrant@hwi.buffalo.edu>
 #
 # This file is part of CrystFEL.
@@ -120,7 +118,7 @@ data = np.asarray(data,dtype=float)
 sys.stdout.write("\r%i peaks found" % n)
 sys.stdout.flush()
 
-print ""
+print("")
 
 x = data[:,0]
 y = data[:,1]
diff --git a/scripts/split-by-mask b/scripts/split-by-mask
index aede3f43..18b672e4 100755
--- a/scripts/split-by-mask
+++ b/scripts/split-by-mask
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Split stream according to an external mask
 #
-# Copyright (c) 2014-2016 Deutsches Elektronen-Synchrotron DESY,
-#                         a research centre of the Helmholtz Association.
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2014-2016 Thomas White <taw@physics.org>
+#    2014-2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -21,7 +21,7 @@ prog = re.compile("^\s+([\d-]+)\s+([\d-]+)\s+([\d-]+)\s+[\d\.\-]+\s+[\d\.\-]+\s+
 def process_refln(fline, g, h, mask):
     match = prog.match(fline)
     if not match:
-        print 'Line not understood, WTF? %s' % fline
+        print('Line not understood, WTF? {}'.format(fline))
 	return
     fs = int(round(float(match.group(4))))
     ss = int(round(float(match.group(5))))
@@ -48,7 +48,7 @@ op.add_option('', '--mask', action='store', type='string', dest='mask_fn',
 opt,arg = op.parse_args(sys.argv)
 
 if not opt.ifn:
-    print "You need at least --input"
+    print("You need at least --input")
     exit(1)
 
 f = open(opt.ifn, 'r')
diff --git a/scripts/sum-peaks b/scripts/sum-peaks
new file mode 100755
index 00000000..eccb9fc6
--- /dev/null
+++ b/scripts/sum-peaks
@@ -0,0 +1,47 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Generate "peak powder" from CrystFEL stream
+#
+# Copyright © 2017 Deutsches Elektronen-Synchrotron DESY,
+#                  a research centre of the Helmholtz Association.
+#
+# Author:
+#    2017 Thomas White <taw@physics.org>
+#
+
+import numpy as np
+import h5py
+import sys
+import re
+
+f = open(sys.argv[1], 'r')
+powder = np.zeros((512,1024), dtype=float)
+peaks = []
+
+prog1 = re.compile("^\s*([\d\-\.]+)\s+([\d\-\.]+)\s+[\d\-\.]+\s+([\d\-\.]+)\s+\S+$")
+
+while True:
+
+    fline = f.readline()
+    if not fline:
+        break
+
+    if fline == '----- End chunk -----\n':
+        if len(peaks) > 10:
+            for p in peaks:
+                powder[p[1],p[0]] += p[2]
+        peaks = []
+
+    match = prog1.match(fline)
+    if match:
+        fs = int(float(match.group(1)))
+        ss = int(float(match.group(2)))
+        intensity = float(match.group(3))
+        peaks.append((fs,ss,intensity))
+
+f.close()
+
+fh = h5py.File("summed-peaks.h5", "w")
+fh.create_dataset("/data", (512,1024), data=powder)
+fh.close()
diff --git a/scripts/truncate-stream b/scripts/truncate-stream
index 1725082a..60679d30 100755
--- a/scripts/truncate-stream
+++ b/scripts/truncate-stream
@@ -1,13 +1,13 @@
 #!/usr/bin/env python
-
+# -*- coding: utf-8 -*-
 #
 # Chop chunks or crystals from one or both ends of a stream
 #
-# Copyright (c) 2014 Deutsches Elektronen-Synchrotron DESY,
-#                    a research centre of the Helmholtz Association.
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
 #
 # Author:
-#    2014 Thomas White <taw@physics.org>
+#    2014-2017 Thomas White <taw@physics.org>
 #
 
 import sys
@@ -66,7 +66,7 @@ def count_crystals(f, g, start_after, stop_after):
                 in_crystal = 0
                 need_end_chunk = 1
 
-    print "Wrote %i crystals to %s" % (n_crystals_written, opt.ofn)
+    print("Wrote {} crystals to {}".format(n_crystals_written, opt.ofn))
 
 def count_chunks(f, g, start_after, stop_after):
     n_chunks_seen = 0
@@ -95,7 +95,7 @@ def count_chunks(f, g, start_after, stop_after):
                 if n_chunks_written == stop_after:
                     break
 
-    print "Wrote %i chunks to %s" % (n_chunks_written, opt.ofn)
+    print("Wrote {} chunks to {}".format(n_chunks_written, opt.ofn))
 
 op = optparse.OptionParser()
 op.add_option('', '--input', action='store', type='string', dest='ifn',
@@ -111,7 +111,7 @@ op.add_option('', '--chunks', action='store_true', dest='chunks',
 opt,arg = op.parse_args(sys.argv)
 
 if not (opt.ifn and opt.ofn):
-    print "You need at least --input and --output"
+    print("You need at least --input and --output")
     exit(1)
 
 f = open(opt.ifn, 'r')
diff --git a/scripts/turbo-index b/scripts/turbo-index-lsf
index c11f521e..00e4ec15 100644
--- a/scripts/turbo-index
+++ b/scripts/turbo-index-lsf
@@ -3,9 +3,9 @@
 RUN=$1
 NOSAMPLE=`echo $RUN | sed -e 's/\-.*$//'`
 
-GEOM=<name of geometry file>
+GEOM=my.geom  # Name of your geometry file
 
-find <path to CXI files>/$RUN -name '*.cxi' > files-${RUN}.lst
+find /path/to/CXI/files/$RUN -name '*.cxi' > files-${RUN}.lst  # Set location of files
 list_events -i files-${RUN}.lst -g $GEOM -o events-${RUN}.lst
 wc -l events-${RUN}.lst
 rm -f split-events-${RUN}.lst files-${RUN}.lst
@@ -18,6 +18,7 @@ for FILE in split-events-${RUN}.lst*; do
         NAME=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${NOSAMPLE}-/"`
         echo "$NAME: $FILE  --->  $STREAM"
 
+	# Set indexing parameters here
         bsub -q psanaq -o $NAME.log -J $NAME -n 12 -R "span[hosts=1]" \
 	            indexamajig \
                            -i $FILE -o $STREAM -j 32 -g $GEOM --peaks=cxi
diff --git a/scripts/turbo-index-slurm b/scripts/turbo-index-slurm
new file mode 100755
index 00000000..a0f1bec5
--- /dev/null
+++ b/scripts/turbo-index-slurm
@@ -0,0 +1,85 @@
+#!/bin/sh
+
+# Split a large indexing job into many small tasks and submit using SLURM
+
+# ./turbo-index my-files.lst label my.geom /location/for/streams
+
+# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+#
+# Authors:
+#   2016      Steve Aplin <steve.aplin@desy.de>
+#   2016-2017 Thomas White <taw@physics.org>
+
+SPLIT=1000  # Size of job chunks
+MAIL=you@example.org  # Email address for SLURM notifications
+
+INPUT=$1
+RUN=$2
+GEOM=$3
+STREAMDIR=$4
+
+# Set up environment here if necessary
+#source /path/to/crystfel/setup.sh
+
+# Generate event list from file above
+list_events -i $INPUT -g $GEOM -o events-${RUN}.lst
+if [ $? != 0 ]; then
+       echo "list_events failed"
+       exit 1
+fi
+# If you are using single-event files instead of multi-event ("CXI") ones,
+# comment out the above lines and uncomment the following one:
+#cp $INPUT events-${RUN}.lst
+
+# Count total number of events
+wc -l events-${RUN}.lst
+
+# Split the events up, will create files with $SPLIT lines
+split -a 3 -d -l $SPLIT events-${RUN}.lst split-events-${RUN}.lst
+
+# Clean up
+rm -f events-${RUN}.lst
+
+# Loop over the event list files, and submit a batch job for each of them
+for FILE in split-events-${RUN}.lst*; do
+
+    # Stream file is the output of crystfel
+    STREAM=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${RUN}.stream/"`
+
+    # Job name
+    NAME=`echo $FILE | sed -e "s/split-events-${RUN}.lst/${RUN}-/"`
+
+    echo "$NAME: $FILE  --->  $STREAM"
+
+    SLURMFILE="${NAME}.sh"
+
+    echo "#!/bin/sh" > $SLURMFILE
+    echo >> $SLURMFILE
+
+    echo "#SBATCH --partition=mypartition" >> $SLURMFILE  # Set your partition here
+    echo "#SBATCH --time=01:00:00" >> $SLURMFILE
+    echo "#SBATCH --nodes=1" >> $SLURMFILE
+    echo "#SBATCH --nice=100" >> $SLURMFILE # Set priority very low to allow other jobs through
+    echo >> $SLURMFILE
+
+    echo "#SBATCH --workdir   $PWD" >> $SLURMFILE
+    echo "#SBATCH --job-name  $NAME" >> $SLURMFILE
+    echo "#SBATCH --output    $NAME-%N-%j.out" >> $SLURMFILE
+    echo "#SBATCH --error     $NAME-%N-%j.err" >> $SLURMFILE
+    echo "#SBATCH --mail-type END" >> $SLURMFILE
+    echo "#SBATCH --mail-user $MAIL" >> $SLURMFILE
+    echo >> $SLURMFILE
+
+    echo "#source /path/to/crystfel/setup.sh" >> $SLURMFILE  # Set up environment here (again) if necessary
+    echo >> $SLURMFILE
+
+    command="indexamajig -i $FILE -o $STREAMDIR/$STREAM"
+    command="$command -j \`nproc\` -g $GEOM"
+    #command="$command --peaks=zaef"  # Indexing parameters here
+
+    echo $command >> $SLURMFILE
+
+    sbatch $SLURMFILE
+
+done
diff --git a/src/cell_explorer.c b/src/cell_explorer.c
index c8bbb069..f6c2c590 100644
--- a/src/cell_explorer.c
+++ b/src/cell_explorer.c
@@ -658,7 +658,9 @@ static void scan_cells(CellWindow *w)
 
 		if ( ignore ) continue;
 
-		cell_get_parameters(cells[i], &a, &b, &c, &al, &be, &ga);
+		if ( cell_get_parameters(cells[i], &a, &b, &c, &al, &be, &ga) ) {
+			continue;
+		}
 		a *= 1e10;  b *= 1e10;  c *= 1e10;
 		al = rad2deg(al);  be = rad2deg(be);  ga = rad2deg(ga);
 
@@ -719,7 +721,10 @@ static void scan_minmax(CellWindow *w)
 		int j;
 		int found = 0;
 
-		cell_get_parameters(w->cells[i], &a, &b, &c, &al, &be, &ga);
+		if ( cell_get_parameters(w->cells[i], &a, &b, &c, &al, &be, &ga) ) {
+			ERROR("Cell %i is bad\n", i);
+			continue;
+		}
 		a *= 1e10;  b *= 1e10;  c *= 1e10;
 		al = rad2deg(al);  be = rad2deg(be);  ga = rad2deg(ga);
 
@@ -1555,6 +1560,8 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	gsl_set_error_handler_off();
+
 	w.cells = NULL;
 	w.indms = NULL;
 	w.n_cells = 0;
diff --git a/src/check_hkl.c b/src/check_hkl.c
index def986c9..c7afb6cc 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -501,7 +501,8 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
 		signed int hs, ks, ls;
 		int bin;
 
-		d = 2.0 * resolution(cell, h, k, l);
+		get_asymm(sym, h, k, l, &hs, &ks, &ls);
+		d = 2.0 * resolution(cell, hs, ks, ls);
 
 		if ( forbidden_reflection(cell, h, k, l) ) continue;
 
@@ -514,7 +515,6 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
 		}
 		if ( bin == -1 ) continue;
 
-		get_asymm(sym, h, k, l, &hs, &ks, &ls);
 		if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
 		add_refl(counted, hs, ks, ls);
 
@@ -530,7 +530,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		signed int h, k, l;
+		signed int h, k, l, hs, ks, ls;
 		double d, val, esd;
 		int bin;
 		int j;
@@ -538,7 +538,8 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
 		get_indices(refl, &h, &k, &l);
 		if ( forbidden_reflection(cell, h, k, l) ) continue;
 
-		d = resolution(cell, h, k, l) * 2.0;
+		get_asymm(sym, h, k, l, &hs, &ks, &ls);
+		d = resolution(cell, hs, ks, ls) * 2.0;
 		val = get_intensity(refl);
 		esd = get_esd_intensity(refl);
 
diff --git a/src/cl-utils.c b/src/cl-utils.c
index dc0a2713..a7e500cd 100644
--- a/src/cl-utils.c
+++ b/src/cl-utils.c
@@ -50,6 +50,15 @@ const char *clError(cl_int err)
 		case CL_SUCCESS :
 		return "no error";
 
+		case CL_DEVICE_NOT_AVAILABLE :
+		return "device not available";
+
+		case CL_DEVICE_NOT_FOUND :
+		return "device not found";
+
+		case CL_INVALID_DEVICE_TYPE :
+		return "invalid device type";
+
 		case CL_INVALID_PLATFORM :
 		return "invalid platform";
 
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 51e51dc9..81f12c94 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -286,13 +286,13 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
 		break;
 
 		case FOM_D1SIG :
-		if ( fabs(i1-i2) < (sig1+sig2)/2.0 ) {
+		if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) {
 			fctx->n_within[bin]++;
 		}
 		break;
 
 		case FOM_D2SIG :
-		if ( fabs(i1-i2) < sig1+sig2 ) {  /* = 2 * (sig1+sig2)/2 */
+		if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) {
 			fctx->n_within[bin]++;
 		}
 		break;
@@ -1125,7 +1125,7 @@ int main(int argc, char *argv[])
 	char *bfile = NULL;
 	char *sym_str = NULL;
 	SymOpList *sym;
-	int ncom, nrej, nneg, nres, nbij, ncen;
+	int ncom, nrej, nmul, nneg, nres, nbij, ncen;
 	RefList *list1_acc;
 	RefList *list2_acc;
 	RefList *list1;
@@ -1149,6 +1149,7 @@ int main(int argc, char *argv[])
 	double min_I = +INFINITY;
 	double max_I = -INFINITY;
 	float highres, lowres;
+	int mul_cutoff = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -1164,6 +1165,7 @@ int main(int argc, char *argv[])
 		{"shell-file",         1, NULL,                7},
 		{"highres",            1, NULL,                8},
 		{"lowres",             1, NULL,                9},
+		{"min-measurements",   1, NULL,               11},
 		{"ignore-negs",        0, &config_ignorenegs,  1},
 		{"zero-negs",          0, &config_zeronegs,    1},
 		{"intensity-shells",   0, &config_intshells,   1},
@@ -1260,6 +1262,13 @@ int main(int argc, char *argv[])
 			rmin_fix = 1.0 / (lowres/1e10);
 			break;
 
+			case 11 :
+			if ( sscanf(optarg, "%i", &mul_cutoff) != 1 ) {
+				ERROR("Invalid value for --min-measurements\n");
+				return 1;
+			}
+			break;
+
 			case '?' :
 			break;
 
@@ -1413,6 +1422,7 @@ int main(int argc, char *argv[])
 	/* Select reflections to be used */
 	ncom = 0;
 	nrej = 0;
+	nmul = 0;
 	nneg = 0;
 	nres = 0;
 	nbij = 0;
@@ -1426,6 +1436,7 @@ int main(int argc, char *argv[])
 		signed int h, k, l;
 		double val1, val2;
 		double esd1, esd2;
+		int mul1, mul2;
 		Reflection *refl2;
 		Reflection *refl1_acc;
 		Reflection *refl2_acc;
@@ -1441,6 +1452,9 @@ int main(int argc, char *argv[])
 		esd1 = get_esd_intensity(refl1);
 		esd2 = get_esd_intensity(refl2);
 
+		mul1 = get_redundancy(refl1);
+		mul2 = get_redundancy(refl2);
+
 		if ( (val1 < sigma_cutoff * esd1)
 		  || (val2 < sigma_cutoff * esd2) )
 		{
@@ -1453,6 +1467,11 @@ int main(int argc, char *argv[])
 			continue;
 		}
 
+		if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) {
+			nmul++;
+			continue;
+		}
+
 		if ( config_zeronegs ) {
 			int d = 0;
 			if ( val1 < 0.0 ) {
@@ -1544,6 +1563,11 @@ int main(int argc, char *argv[])
 		       " negative intensities which were set to zero.\n", nneg);
 	}
 
+	if ( nmul > 0 ) {
+		STATUS("%i reflection pairs rejected because either or both"
+		       " versions had too few measurements.\n", nmul);
+	}
+
 	if ( nres > 0 ) {
 		STATUS("%i reflection pairs rejected because either or both"
 		       " versions were outside the resolution range.\n", nres);
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 331170ae..22abcfd2 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -379,6 +379,7 @@ struct gpu_context *setup_gpu(int no_sfac,
 	float *flags_ptr;
 	size_t maxwgsize;
 	int i;
+	int have_ctx = 0;
 	char cflags[512] = "";
 	char *insert_stuff = NULL;
 
@@ -393,16 +394,37 @@ struct gpu_context *setup_gpu(int no_sfac,
 		ERROR("Couldn't find at least one platform!\n");
 		return NULL;
 	}
-	prop[0] = CL_CONTEXT_PLATFORM;
-	prop[1] = (cl_context_properties)platforms[0];
-	prop[2] = 0;
 
-	gctx = malloc(sizeof(*gctx));
-	gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
-	                                    NULL, NULL, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create OpenCL context: %i\n", err);
-		free(gctx);
+	/* Find a GPU platform in the list */
+	for ( i=0; i<nplat; i++ ) {
+
+		prop[0] = CL_CONTEXT_PLATFORM;
+		prop[1] = (cl_context_properties)platforms[i];
+		prop[2] = 0;
+
+		gctx = malloc(sizeof(*gctx));
+		gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
+		                                    NULL, NULL, &err);
+
+		if ( err != CL_SUCCESS ) {
+			if ( err == CL_DEVICE_NOT_FOUND ) {
+				/* No GPU device in this platform */
+				continue; /* Try next platform */
+			} else {
+				ERROR("Couldn't create OpenCL context: %s (%i)\n",
+				clError(err), err);
+				free(gctx);
+				return NULL;
+			}
+		} else {
+			STATUS("Using OpenCL platform %i (%i total)\n", i, nplat);
+			have_ctx = 1;
+			break;
+		}
+	}
+
+	if ( !have_ctx ) {
+		ERROR("Couldn't find a GPU device in any platform.\n");
 		return NULL;
 	}
 
diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 87c23f4e..0336f25e 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -74,8 +74,8 @@ static void displaywindow_error(DisplayWindow *dw, const char *message)
 static gint displaywindow_closed(GtkWidget *window, DisplayWindow *dw)
 {
 
-	if ( dw->hdfile != NULL ) {
-		hdfile_close(dw->hdfile);
+	if ( dw->imagefile != NULL ) {
+		imagefile_close(dw->imagefile);
 	}
 
 	if ( dw->surf != NULL ) cairo_surface_destroy(dw->surf);
@@ -725,7 +725,7 @@ static gint displaywindow_set_binning(GtkWidget *widget, DisplayWindow *dw)
 		return 0;
 	}
 
-	if ( dw->hdfile == NULL ) {
+	if ( dw->imagefile == NULL ) {
 		return 0;
 	}
 
@@ -852,8 +852,8 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw)
 		return 0;
 	}
 
-	if ( dw->hdfile == NULL ) {
-	    return 0;
+	if ( dw->imagefile == NULL ) {
+		return 0;
 	}
 
 	bd = malloc(sizeof(BoostIntDialog));
@@ -927,8 +927,8 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event)
 	float **old_dp = dw->image->dp;
 	int **old_bad = dw->image->bad;
 
-	fail = hdf5_read2(dw->hdfile, dw->image,
-                          dw->ev_list->events[new_event], 0);
+	fail = imagefile_read(dw->imagefile, dw->image,
+	                      dw->ev_list->events[new_event]);
 	if ( fail ) {
 		ERROR("Couldn't load image");
 		dw->image->dp = old_dp;
@@ -1046,8 +1046,8 @@ static gint displaywindow_set_newevent(GtkWidget *widget, DisplayWindow *dw)
 		return 0;
 	}
 
-	if ( dw->hdfile == NULL ) {
-	    return 0;
+	if ( dw->imagefile == NULL ) {
+		return 0;
 	}
 
 	ed = malloc(sizeof(EventDialog));
@@ -1162,7 +1162,7 @@ static gint displaywindow_set_ringradius(GtkWidget *widget, DisplayWindow *dw)
 		return 0;
 	}
 
-	if ( dw->hdfile == NULL ) {
+	if ( dw->imagefile == NULL ) {
 		return 0;
 	}
 
@@ -1796,7 +1796,7 @@ static gint displaywindow_show_numbers(GtkWidget *widget, DisplayWindow *dw)
 		return 0;
 	}
 
-	if ( dw->hdfile == NULL ) {
+	if ( dw->imagefile == NULL ) {
 		return 0;
 	}
 
@@ -2306,12 +2306,14 @@ struct newhdf {
 	char name[1024];
 };
 
+/* New HDF5 element selected from menu */
 static gint displaywindow_newhdf(GtkMenuItem *item, struct newhdf *nh)
 {
 	gboolean a;
 	int fail;
 
 	if ( nh->dw->not_ready_yet ) return 0;
+	assert(imagefile_get_type(nh->dw->imagefile) == IMAGEFILE_HDF5);
 
 	a = gtk_check_menu_item_get_active(GTK_CHECK_MENU_ITEM(nh->widget));
 	if ( !a ) return 0;
@@ -2320,7 +2322,8 @@ static gint displaywindow_newhdf(GtkMenuItem *item, struct newhdf *nh)
 	 * one */
 	free_detector_geometry(nh->dw->image->det);
 	nh->dw->image->det = NULL;
-	fail = hdf5_read(nh->dw->hdfile, nh->dw->image, nh->name, 0);
+	fail = hdf5_read(imagefile_get_hdfile(nh->dw->imagefile),
+	                 nh->dw->image, nh->name, 0);
 	if ( fail ) {
 		ERROR("Couldn't load image");
 		return 1;
@@ -2456,7 +2459,12 @@ static int displaywindow_update_menus(DisplayWindow *dw, const char *selectme)
 	GtkWidget *ms;
 	GtkWidget *w;
 
-	ms = displaywindow_createhdfmenus(dw->hdfile, dw, selectme);
+	if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+		return 0;
+	}
+
+	ms = displaywindow_createhdfmenus(imagefile_get_hdfile(dw->imagefile),
+	                                  dw, selectme);
 
 	if ( ms == NULL ) return 1;
 
@@ -2828,8 +2836,8 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
 	dw->image->lambda = 0.0;
 	dw->image->filename = filename;
 
-	dw->hdfile = hdfile_open(filename);
-	if ( dw->hdfile == NULL ) {
+	dw->imagefile = imagefile_open(filename);
+	if ( dw->imagefile == NULL ) {
 		ERROR("Couldn't open file: %s\n", filename);
 		free(dw->geom_filename);
 		free(dw);
@@ -2842,11 +2850,19 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
 	}
 
 	if ( dw->image->det != NULL && ( dw->image->det->path_dim != 0 ||
-	                                 dw->image->det->dim_dim != 0  )) {
+	                                 dw->image->det->dim_dim != 0  ))
+	{
+		struct hdfile *hdfile;
+
+		if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+			ERROR("Multi-event geometry, but not HDF5 file!\n");
+			return NULL;
+		}
+		hdfile = imagefile_get_hdfile(dw->imagefile);
 
 		dw->multi_event = 1;
 
-		dw->ev_list = fill_event_list(dw->hdfile, dw->image->det);
+		dw->ev_list = fill_event_list(hdfile, dw->image->det);
 
 		if ( dw->ev_list == NULL ) {
 			ERROR("Error while parsing file structure\n");
@@ -2883,18 +2899,26 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
 				dw->curr_event = 0;
 				ev = dw->ev_list->events[dw->curr_event];
 			}
-			check = hdf5_read2(dw->hdfile, dw->image, ev, 0);
+			check = imagefile_read(dw->imagefile, dw->image, ev);
 		} else {
-			check = hdf5_read2(dw->hdfile, dw->image, NULL, 0);
+			check = imagefile_read(dw->imagefile, dw->image, NULL);
 		}
 
 	} else {
-		check =	hdf5_read(dw->hdfile, dw->image, element, 0);
+		if ( element != NULL ) {
+			if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+				ERROR("-e/--image requiest an HDF5 file\n");
+				return NULL;
+			}
+			hdfile_set_image(imagefile_get_hdfile(dw->imagefile),
+			                 element);
+		}
+		check =	imagefile_read_simple(dw->imagefile, dw->image);
 		dw->simple = 1;
 	}
 	if ( check ) {
 		ERROR("Couldn't load file\n");
-		hdfile_close(dw->hdfile);
+		imagefile_close(dw->imagefile);
 		free(dw->geom_filename);
 		return NULL;
 	}
diff --git a/src/dw-hdfsee.h b/src/dw-hdfsee.h
index 85c82ac7..fdbea231 100644
--- a/src/dw-hdfsee.h
+++ b/src/dw-hdfsee.h
@@ -3,12 +3,12 @@
  *
  * Quick yet non-crappy HDF viewer
  *
- * 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:
- *   2009-2014 Thomas White <taw@physics.org>
+ *   2009-2017 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
  *   2014      Takanori Nakane <nakane.t@gmail.com>
  *   2012      Richard Kirian
@@ -101,8 +101,8 @@ typedef struct {
 
 	int             simple;
 
-	struct hdfile	*hdfile;
-	struct image	*image;
+	struct imagefile *imagefile;
+	struct image	 *image;
 
 	char            *geom_filename;
 	char            *rg_coll_name;
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
index e0337358..da9b16c4 100644
--- a/src/geoptimiser.c
+++ b/src/geoptimiser.c
@@ -631,10 +631,18 @@ static int add_distance_to_list(struct gpanel *gp,
                                 double *det_shift)
 {
 	int pix_index;
+	int ifs, iss;
 	double rfs, rss;
 	double crx, cry;
 
-	pix_index = ((int)rint(imfe->fs) + gp->p->w*(int)rint(imfe->ss));
+	ifs = imfe->fs;
+	iss = imfe->ss;  /* Explicit rounding towards zero (truncation) */
+	pix_index = ifs + gp->p->w*iss;
+
+	if ( (ifs >= gp->p->w) || (iss >= gp->p->h) ) {
+		ERROR("Peak is outside panel!\n");
+		return 1;
+	}
 
 	if ( gp->num_pix_displ[pix_index] > 0 ) {
 
@@ -782,7 +790,14 @@ static int compute_pixel_displacements(struct image *images, int n_images,
 
 				r = add_distance_to_list(gp, imfe, refl, fx, fy,
 				                         det_shift);
-				if ( r ) return r;
+				if ( r ) {
+					ERROR("Error processing peak %f,%f "
+					      "(panel %s), image %s %s\n",
+					      imfe->fs, imfe->ss, gp->p->name,
+					      images[cp].filename,
+					      get_event_string(images[cp].event));
+					return r;
+				}
 
 			}
 		}
@@ -1070,12 +1085,7 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
 
 			/* Calculate corner adjustment
 			 * NB All panels follow the first one */
-			if ( ip == 0 ) {
-
-				new_cnx = p->cnx * cs;
-				new_cny = p->cny * cs;
-
-			} else {
+			if ( ip > 0 && connected->rigid_groups[di]->n_panels == 2 && !gparams->no_cspad ) {
 
 				struct panel *p0;
 				double delta_x, delta_y, delta;
@@ -1090,6 +1100,10 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
 				new_cnx = p0->cnx + delta*p0->fsx;
 				new_cny = p0->cny + delta*p0->fsy;
 
+			} else {
+
+				new_cnx = p->cnx * cs;
+				new_cny = p->cny * cs;
 			}
 
 			/* The average displacements now need to be updated */
@@ -2488,8 +2502,10 @@ int main(int argc, char *argv[])
 	}
 
 #ifdef HAVE_SAVE_TO_PNG
+#if !GLIB_CHECK_VERSION(2,36,0)
 	g_type_init();
 #endif
+#endif
 
 	ret_val = optimize_geometry(gparams, det, quadrants, connected);
 
diff --git a/src/get_hkl.c b/src/get_hkl.c
index d6efe747..77f34da2 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -98,7 +98,9 @@ static void show_help(const char *s)
 
 static void copy_notes(RefList *out, RefList *in)
 {
-	reflist_add_notes(out, reflist_get_notes(in));
+	if ( reflist_get_notes(in) != NULL ) {
+		reflist_add_notes(out, reflist_get_notes(in));
+	}
 }
 
 
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index f5493453..60400ad8 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -372,10 +372,10 @@ static void run_work(const struct index_args *iargs, Stream *st,
 	}
 
 	time_accounts_set(taccs, TACC_FINALCLEANUP);
-	cleanup_indexing(iargs->indm, iargs->ipriv);
+	cleanup_indexing(iargs->ipriv);
 	free_detector_geometry(iargs->det);
 	free(iargs->hdf5_peak_path);
-	free_copy_hdf5_field_list(iargs->copyme);
+	free_imagefile_field_list(iargs->copyme);
 	cell_free(iargs->cell);
 	if ( iargs->profile ) time_accounts_print(taccs);
 	time_accounts_free(taccs);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index a1997b9e..ce87fde6 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,16 +3,17 @@
  *
  * Index patterns, output hkl+intensity etc.
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *   2011      Richard Kirian
  *   2012      Lorenzo Galli
  *   2012      Chunhong Yoon
+ *   2017      Valerio Mariani <valerio.mariani@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -88,6 +89,7 @@ static void show_help(const char *s)
 "     --peaks=<method>     Use 'method' for finding peaks.  Choose from:\n"
 "                           zaef  : Use Zaefferer (2000) gradient detection.\n"
 "                                    This is the default method.\n"
+"                           peakfinder8: Use Peakfinder8 algorithm.\n"
 "                           hdf5  : Get from a table in HDF5 file.\n"
 "                           cxi   : Get from CXI format HDF5 file.\n"
 "     --hdf5-peaks=<p>     Find peaks table in HDF5 file here.\n"
@@ -95,32 +97,52 @@ static void show_help(const char *s)
 "     --integration=<meth> Perform final pattern integration using <meth>.\n"
 "\n\n"
 "For more control over the process, you might need:\n\n"
-"    --tolerance=<tol>   Set the tolerances for cell comparison.\n"
-"                          Default: 5,5,5,1.5.\n"
-"    --filter-noise      Apply an aggressive noise filter which sets all\n"
-"                         pixels in each 3x3 region to zero if any of them\n"
-"                         have negative values.  Intensity measurement will\n"
-"                         be performed on the image as it was before this.\n"
-"    --median-filter=<n> Apply a median filter to the image data.  Intensity\n"
-"                         measurement will be performed on the image as it\n"
-"                         was before this.  The side length of the median\n"
-"                         filter box will be 2<n>+1.  Default: 0 (no filter).\n"
-"    --no-sat-corr       Don't correct values of saturated peaks using a\n"
-"                         table included in the HDF5 file.\n"
-"    --threshold=<n>     Only accept peaks above <n> ADU.  Default: 800.\n"
-"    --min-gradient=<n>  Minimum squared gradient for Zaefferer peak search.\n"
-"                         Default: 100,000.\n"
-"    --min-snr=<n>       Minimum signal-to-noise ratio for peaks.\n"
-"                         Default: 5.\n"
-"    --check-hdf5-snr    Check SNR for peaks from --peaks=hdf5.\n"
-"    --peak-radius=<r>   Integration radii for peak search.\n"
-"    --int-radius=<r>    Set the integration radii.  Default: 4,5,7.\n"
-"    --push-res=<n>      Integrate higher than apparent resolution cutoff.\n"
-"    --highres=<n>       Absolute resolution cutoff in Angstroms.\n"
-"    --fix-profile-radius Fix the reciprocal space profile radius for spot\n"
-"                         prediction (default: automatically determine.\n"
-"    --fix-bandwidth     Set the bandwidth for spot prediction.\n"
-"    --fix-divergence    Set the divergence (full angle) for spot prediction.\n"
+"    --tolerance=<tol>     Set the tolerances for cell comparison.\n"
+"                            Default: 5,5,5,1.5.\n"
+"    --filter-noise        Apply an aggressive noise filter which sets all\n"
+"                           pixels in each 3x3 region to zero if any of them\n"
+"                           have negative values.  Intensity measurement will\n"
+"                           be performed on the image as it was before this.\n"
+"    --median-filter=<n>   Apply a median filter to the image data. Intensity\n"
+"                           measurement will be performed on the image as it\n"
+"                           was before this.  The side length of the median\n"
+"                           filter box will be 2<n>+1.\n"
+"                           Default: 0 (no filter).\n"
+"    --no-sat-corr         Don't correct values of saturated peaks using a\n"
+"                           table included in the HDF5 file.\n"
+"    --threshold=<n>       Only accept peaks above <n> ADU in both the\n"
+"	                    Zaefferer and Peakfinder8 algorithms.\n"
+"                           Default: 800.\n"
+"    --min-gradient=<n>    Minimum squared gradient for Zaefferer peak\n"
+"                           search. Default: 100,000.\n"
+"    --min-snr=<n>         Minimum signal-to-noise ratio for peaks, with both\n"
+"                           Zaefferer and Peakfinder8 algorithms.\n"
+"                           Default: 5.\n"
+"    --min-pix-count=<n>   Only accept peaks if they include more than <n>\n"
+"                           pixels, in the Peakfinder8 algorithm.\n"
+"                           Default: 2.\n"
+"    --max-pix-count=<n>   Only accept peaks if they include less than <n>\n"
+"                           pixels, in the Peakfinder8 algorithm.\n"
+"                           Default: 200.\n"
+"    --min-peaks=<n>       Minimum number of peaks for indexing.\n"
+"    --local-bg-radius=<n> Radius (in pixels) to use for the estimation of\n"
+"	                    local background in the Peakfinder8 algorithm.\n"
+"                           Default: 3.\n"
+"    --min-res=<n>         Only accept peaks if they lay at more than <n>\n"
+"	                    pixels from the center of the detector, in the\n"
+"                           peakfinder8 algorithm. Default: 0.\n"
+"    --max-res=<n>         Only accept peaks if they lay at less than <n>\n"
+"	                    pixels from the center of the detector, in the\n"
+"                           peakfinder8 algorithm. Default: 1200.\n"
+"    --check-hdf5-snr      Check SNR for peaks from --peaks=hdf5.\n"
+"    --peak-radius=<r>     Integration radii for peak search.\n"
+"    --int-radius=<r>      Set the integration radii.  Default: 4,5,7.\n"
+"    --push-res=<n>        Integrate higher than apparent resolution cutoff.\n"
+"    --highres=<n>         Absolute resolution cutoff in Angstroms.\n"
+"    --fix-profile-radius  Fix the reciprocal space profile radius for spot\n"
+"                           prediction (default: automatically determine.\n"
+"    --fix-bandwidth       Set the bandwidth for spot prediction.\n"
+"    --fix-divergence      Set the divergence (full angle) for spot prediction.\n"
 "\n"
 "\nFor time-resolved stuff, you might want to use:\n\n"
 "     --copy-hdf5-field <f>  Copy the value of field <f> into the stream. You\n"
@@ -131,16 +153,19 @@ static void show_help(const char *s)
 " --temp-dir=<path>        Put the temporary folder under <path>.\n"
 "\n"
 "\nOptions you probably won't need:\n\n"
-"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
-"     --no-use-saturated   During the initial peak search, reject\n"
-"                           peaks which contain pixels above max_adu.\n"
-"     --no-revalidate      Don't re-integrate and check HDF5 peaks for\n"
-"                           validity.\n"
-"     --no-peaks-in-stream Do not record peak search results in the stream.\n"
-"     --no-refls-in-stream Do not record integrated reflections in the stream.\n"
-"     --int-diag=<cond>    Show debugging information about reflections.\n"
-"     --no-refine          Skip the prediction refinement step.\n"
-"     --profile            Show timing data for performance monitoring.\n"
+"     --no-check-prefix     Don't attempt to correct the --prefix.\n"
+"     --no-use-saturated    During the initial peak search, reject\n"
+"                            peaks which contain pixels above max_adu.\n"
+"     --no-revalidate       Don't re-integrate and check HDF5 peaks for\n"
+"                            validity.\n"
+"     --no-peaks-in-stream  Do not record peak search results in the stream.\n"
+"     --no-refls-in-stream  Do not record integrated reflections in the stream.\n"
+"     --no-non-hits-in-stream  Do not include non-hit frames in the stream.\n"
+"                               (see --min-peaks)\n"
+"     --int-diag=<cond>     Show debugging information about reflections.\n"
+"     --no-refine           Skip the prediction refinement step.\n"
+"     --profile             Show timing data for performance monitoring.\n"
+"     --no-half-pixel-shift Don't offset the HDF5 peak locations by 0.5 px.\n"
 "\nLow-level options for the felix indexer:\n\n"
 "     --felix-options      Change the default arguments passed to the indexer.\n"
 "                          Given as a list of comma separated list of \n"
@@ -150,9 +175,9 @@ static void show_help(const char *s)
 }
 
 
-static void add_geom_beam_stuff_to_copy_hdf5(struct copy_hdf5_field *copyme,
-                                             struct detector *det,
-                                             struct beam_params *beam)
+static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copyme,
+                                              struct detector *det,
+                                              struct beam_params *beam)
 {
 	int i;
 
@@ -161,12 +186,12 @@ static void add_geom_beam_stuff_to_copy_hdf5(struct copy_hdf5_field *copyme,
 		struct panel *p = &det->panels[i];
 
 		if ( p->clen_from != NULL ) {
-			add_copy_hdf5_field(copyme, p->clen_from);
+			add_imagefile_field(copyme, p->clen_from);
 		}
 	}
 
 	if ( beam->photon_energy_from != NULL ) {
-		add_copy_hdf5_field(copyme, beam->photon_energy_from);
+		add_imagefile_field(copyme, beam->photon_energy_from);
 	}
 }
 
@@ -181,8 +206,6 @@ int main(int argc, char *argv[])
 	int config_checkprefix = 1;
 	int config_basename = 0;
 	int integrate_saturated = 0;
-	IndexingMethod *indm;
-	IndexingPrivate **ipriv;
 	char *indm_str = NULL;
 	char *cellfile = NULL;
 	char *prefix = NULL;
@@ -214,11 +237,17 @@ int main(int argc, char *argv[])
 	iargs.threshold = 800.0;
 	iargs.min_gradient = 100000.0;
 	iargs.min_snr = 5.0;
+	iargs.min_pix_count = 2;
+	iargs.max_pix_count = 200;
+	iargs.min_res = 0;
+	iargs.max_res = 1200;
+	iargs.local_bg_radius = 3;
 	iargs.check_hdf5_snr = 0;
 	iargs.det = NULL;
 	iargs.peaks = PEAK_ZAEF;
 	iargs.beam = &beam;
 	iargs.hdf5_peak_path = NULL;
+	iargs.half_pixel_shift = 1;
 	iargs.copyme = NULL;
 	iargs.pk_inn = -1.0;
 	iargs.pk_mid = -1.0;
@@ -230,13 +259,14 @@ int main(int argc, char *argv[])
 	iargs.no_revalidate = 0;
 	iargs.stream_peaks = 1;
 	iargs.stream_refls = 1;
+	iargs.stream_nonhits = 1;
 	iargs.int_diag = INTDIAG_NONE;
-	iargs.copyme = new_copy_hdf5_field_list();
+	iargs.copyme = new_imagefile_field_list();
+	iargs.min_peaks = 0;
 	if ( iargs.copyme == NULL ) {
 		ERROR("Couldn't allocate HDF5 field list.\n");
 		return 1;
 	}
-	iargs.indm = NULL;  /* No default */
 	iargs.ipriv = NULL;  /* No default */
 	iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL);
 	iargs.push_res = 0.0;
@@ -268,12 +298,14 @@ int main(int argc, char *argv[])
 		{"basename",           0, &config_basename,          1},
 		{"no-peaks-in-stream", 0, &iargs.stream_peaks,       0},
 		{"no-refls-in-stream", 0, &iargs.stream_refls,       0},
+		{"no-non-hits-in-stream", 0, &iargs.stream_nonhits,  0},
 		{"integrate-saturated",0, &integrate_saturated,      1},
 		{"no-use-saturated",   0, &iargs.use_saturated,      0},
 		{"no-revalidate",      0, &iargs.no_revalidate,      1},
 		{"check-hdf5-snr",     0, &iargs.check_hdf5_snr,     1},
 		{"no-refine",          0, &no_refine,                1},
 		{"profile",            0, &iargs.profile,            1},
+		{"no-half-pixel-shift",0, &iargs.half_pixel_shift,   0},
 
 		/* Long-only options which don't actually do anything */
 		{"no-sat-corr",        0, &iargs.satcorr,            0},
@@ -306,6 +338,12 @@ int main(int argc, char *argv[])
 		{"fix-bandwidth",      1, NULL,               23},
 		{"fix-divergence",     1, NULL,               24},
 		{"felix-options",      1, NULL,               25},
+	        {"min-pix-count",      1, NULL,               26},
+	        {"max-pix-count",      1, NULL,               27},
+	        {"local-bg-radius",    1, NULL,               28},
+	        {"min-res",            1, NULL,               29},
+	        {"max-res",            1, NULL,               30},
+	        {"min-peaks",          1, NULL,               31},
 
 		{0, 0, NULL, 0}
 	};
@@ -400,7 +438,7 @@ int main(int argc, char *argv[])
 			break;
 
 			case 10 :
-			add_copy_hdf5_field(iargs.copyme, optarg);
+			add_imagefile_field(iargs.copyme, optarg);
 			break;
 
 			case 11 :
@@ -489,6 +527,30 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+		        case 26:
+			iargs.min_pix_count = atoi(optarg);
+			break;
+
+		        case 27:
+			iargs.max_pix_count = atoi(optarg);
+			break;
+
+		        case 28:
+			iargs.local_bg_radius = atoi(optarg);
+			break;
+
+		        case 29:
+			iargs.min_res = atoi(optarg);
+			break;
+
+		        case 30:
+			iargs.max_res = atoi(optarg);
+			break;
+
+			case 31:
+			iargs.min_peaks = atoi(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -541,6 +603,8 @@ int main(int argc, char *argv[])
 	}
 	if ( strcmp(speaks, "zaef") == 0 ) {
 		iargs.peaks = PEAK_ZAEF;
+	} else if ( strcmp(speaks, "peakfinder8") == 0 ) {
+		iargs.peaks = PEAK_PEAKFINDER8;
 	} else if ( strcmp(speaks, "hdf5") == 0 ) {
 		iargs.peaks = PEAK_HDF5;
 	} else if ( strcmp(speaks, "cxi") == 0 ) {
@@ -574,7 +638,7 @@ int main(int argc, char *argv[])
 		      geom_filename);
 		return 1;
 	}
-	add_geom_beam_stuff_to_copy_hdf5(iargs.copyme, iargs.det, iargs.beam);
+	add_geom_beam_stuff_to_field_list(iargs.copyme, iargs.det, iargs.beam);
 
 	/* If no peak path from geometry file, use these (but see later) */
 	if ( iargs.hdf5_peak_path == NULL ) {
@@ -591,35 +655,6 @@ int main(int argc, char *argv[])
 		iargs.hdf5_peak_path = command_line_peak_path;
 	}
 
-	/* Parse indexing methods */
-	if ( indm_str == NULL ) {
-
-		STATUS("You didn't specify an indexing method, so I  won't try "
-		       " to index anything.\n"
-		       "If that isn't what you wanted, re-run with"
-		       " --indexing=<methods>.\n");
-		indm = NULL;
-
-	} else {
-
-		int i = 0;
-
-		indm = build_indexer_list(indm_str);
-		if ( indm == NULL ) {
-			ERROR("Invalid indexer list '%s'\n", indm_str);
-			return 1;
-		}
-		free(indm_str);
-
-		/* If --no-refine, unset the refinement flag on all methods */
-		if ( no_refine ) {
-			while ( indm[i] != INDEXING_NONE ) {
-				indm[i] &= ~INDEXING_REFINE;
-				i++;
-			}
-		}
-	}
-
 	/* Parse integration method */
 	if ( int_str != NULL ) {
 
@@ -755,27 +790,34 @@ int main(int argc, char *argv[])
 	}
 	free(outfile);
 
-	/* Prepare the indexer */
-	if ( indm != NULL ) {
-		ipriv = prepare_indexing(indm, iargs.cell, iargs.det,
-		                         iargs.tols, iargs.felix_options);
-		if ( ipriv == NULL ) {
-			ERROR("Failed to prepare indexing.\n");
+	/* Prepare the indexing system */
+	if ( indm_str == NULL ) {
+
+		STATUS("You didn't specify an indexing method, so I  won't try "
+		       " to index anything.\n"
+		       "If that isn't what you wanted, re-run with"
+		       " --indexing=<methods>.\n");
+		iargs.ipriv = NULL;
+
+	} else {
+
+		iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det,
+		                             iargs.tols, no_refine,
+		                             iargs.felix_options);
+		if ( iargs.ipriv == NULL ) {
+			ERROR("Failed to set up indexing system\n");
 			return 1;
 		}
-	} else {
-		ipriv = NULL;
+		free(indm_str);
+
 	}
 
 	gsl_set_error_handler_off();
 
-	iargs.indm = indm;
-	iargs.ipriv = ipriv;
-
 	create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
 	               st, tempdir);
 
-	free_copy_hdf5_field_list(iargs.copyme);
+	free_imagefile_field_list(iargs.copyme);
 	cell_free(iargs.cell);
 	free(iargs.beam->photon_energy_from);
 	free(prefix);
@@ -783,7 +825,7 @@ int main(int argc, char *argv[])
 	free_detector_geometry(iargs.det);
 	free(iargs.hdf5_peak_path);
 	close_stream(st);
-	cleanup_indexing(indm, ipriv);
+	cleanup_indexing(iargs.ipriv);
 
 	return 0;
 }
diff --git a/src/merge.c b/src/merge.c
index 8d1fae0f..9734c469 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -237,10 +237,7 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 	Reflection *refl;
 	RefListIterator *iter;
 
-	if ( n == 0 ) {
-		ERROR("No crystals!\n");
-		return NULL;
-	}
+	if ( n == 0 ) return NULL;
 
 	full = reflist_new();
 
diff --git a/src/partialator.c b/src/partialator.c
index 09feebb4..569145e8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -191,6 +191,12 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 	snprintf(tmp, 1024, "%s1", outfile);
 	split = merge_intensities(crystals1, n_crystals1, nthreads,
 		                  pmodel, min_measurements, push_res, 1);
+
+	if ( split == NULL ) {
+		ERROR("Not enough crystals for two way split!\n");
+		return;
+	}
+
 	STATUS("Writing two-way split to %s ", tmp);
 	write_reflist_2(tmp, split, sym);
 	reflist_free(split);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index e1e04789..e9fc9916 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -745,7 +745,7 @@ int main(int argc, char *argv[])
 		       "each of a, b and c\n");
 	}
 	if ( intfile == NULL ) {
-		STATUS("               Full intensities: all equal");
+		STATUS("               Full intensities: all equal\n");
 	} else {
 		STATUS("               Full intensities: from %s\n", intfile);
 	}
diff --git a/src/process_hkl.c b/src/process_hkl.c
index d2cf640b..301bc6e4 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -287,7 +287,8 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
 		if ( isnan(scale) ) return 1;
 		if ( scale <= 0.0 ) return 1;
 		if ( stat != NULL ) {
-			fprintf(stat, "%s %f %f\n", image->filename, scale, cc);
+			fprintf(stat, "%s %s %f %f\n", image->filename,
+			        get_event_string(image->event), scale, cc);
 		}
 
 	} else {
diff --git a/src/process_image.c b/src/process_image.c
index bcaee543..498b3398 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -8,7 +8,7 @@
  *
  * Authors:
  *   2010-2017 Thomas White <taw@physics.org>
- *   2014      Valerio Mariani
+ *   2014-2017 Valerio Mariani <valerio.mariani@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -102,7 +102,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
                    int serial, struct sb_shm *sb_shared, TimeAccounts *taccs)
 {
 	int check;
-	struct hdfile *hdfile;
+	struct imagefile *imfile;
 	struct image image;
 	int i;
 	int r;
@@ -125,15 +125,15 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	time_accounts_set(taccs, TACC_HDF5OPEN);
 	sb_shared->pings[cookie]++;
-	hdfile = hdfile_open(image.filename);
-	if ( hdfile == NULL ) {
+	imfile = imagefile_open(image.filename);
+	if ( imfile == NULL ) {
 		ERROR("Couldn't open file: %s\n", image.filename);
 		return;
 	}
 
 	time_accounts_set(taccs, TACC_HDF5READ);
 	sb_shared->pings[cookie]++;
-	check = hdf5_read2(hdfile, &image, image.event, 0);
+	check = imagefile_read(imfile, &image, image.event);
 	if ( check ) {
 		return;
 	}
@@ -159,8 +159,14 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	sb_shared->pings[cookie]++;
 	switch ( iargs->peaks ) {
 
+		struct hdfile *hdfile;
+
 		case PEAK_HDF5:
-		if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) {
+		hdfile = imagefile_get_hdfile(imfile);
+		if ( (hdfile == NULL)
+		  || (get_peaks_2(&image, hdfile, iargs->hdf5_peak_path,
+		                  iargs->half_pixel_shift)) )
+		{
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
@@ -172,8 +178,12 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		break;
 
 		case PEAK_CXI:
-		if ( get_peaks_cxi(&image, hdfile, iargs->hdf5_peak_path,
-		                   pargs->filename_p_e) ) {
+		hdfile = imagefile_get_hdfile(imfile);
+		if ( (hdfile == NULL)
+		 ||  (get_peaks_cxi_2(&image, hdfile, iargs->hdf5_peak_path,
+		                      pargs->filename_p_e,
+		                      iargs->half_pixel_shift)) )
+		{
 			ERROR("Failed to get peaks from CXI file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
@@ -191,6 +201,28 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		             iargs->use_saturated);
 		break;
 
+	        case PEAK_PEAKFINDER8:
+		if ( search_peaks_peakfinder8(&image, 2048,
+		                               iargs->threshold,
+		                               iargs->min_snr,
+		                               iargs->min_pix_count,
+		                               iargs->max_pix_count,
+		                               iargs->local_bg_radius,
+		                               iargs->min_res,
+		                               iargs->max_res,
+		                               iargs->use_saturated) ) {
+			if ( image.event != NULL ) {
+				ERROR("Failed to find peaks in image %s"
+				      "(event %s).\n", image.filename,
+				      get_event_string(image.event));
+			} else {
+				ERROR("Failed to find peaks in image %s.",
+				      image.filename);
+			}
+
+		}
+		break;
+
 	}
 
 	restore_image_data(image.dp, image.det, prefilter);
@@ -201,7 +233,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	if ( r ) {
 		ERROR("Failed to chdir to temporary folder: %s\n",
 		      strerror(errno));
-		hdfile_close(hdfile);
+		imagefile_close(imfile);
 		return;
 	}
 
@@ -217,15 +249,30 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		image.bw = 0.00000001;
 	}
 
+	if ( image_feature_count(image.features) < iargs->min_peaks ) {
+		r = chdir(rn);
+		if ( r ) {
+			ERROR("Failed to chdir: %s\n", strerror(errno));
+			imagefile_close(imfile);
+			return;
+		}
+		free(rn);
+
+		if ( iargs->stream_nonhits ) {
+			goto streamwrite;
+		} else {
+			goto out;
+		}
+	}
+
 	/* Index the pattern */
 	time_accounts_set(taccs, TACC_INDEXING);
-	index_pattern_2(&image, iargs->indm, iargs->ipriv,
-	                &sb_shared->pings[cookie]);
+	index_pattern_2(&image, iargs->ipriv, &sb_shared->pings[cookie]);
 
 	r = chdir(rn);
 	if ( r ) {
 		ERROR("Failed to chdir: %s\n", strerror(errno));
-		hdfile_close(hdfile);
+		imagefile_close(imfile);
 		return;
 	}
 	free(rn);
@@ -261,9 +308,10 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	                iargs->int_diag_k, iargs->int_diag_l,
 	                &sb_shared->term_lock);
 
+streamwrite:
 	time_accounts_set(taccs, TACC_WRITESTREAM);
 	sb_shared->pings[cookie]++;
-	ret = write_chunk(st, &image, hdfile,
+	ret = write_chunk(st, &image, imfile,
 	                  iargs->stream_peaks, iargs->stream_refls,
 	                  pargs->filename_p_e->ev);
 	if ( ret != 0 ) {
@@ -280,6 +328,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		       get_event_string(image.event));
 	}
 
+out:
 	/* Count crystals which are still good */
 	time_accounts_set(taccs, TACC_TOTALS);
 	sb_shared->pings[cookie]++;
@@ -313,5 +362,5 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	image_feature_list_free(image.features);
 	free_detector_geometry(image.det);
-	hdfile_close(hdfile);
+	imagefile_close(imfile);
 }
diff --git a/src/process_image.h b/src/process_image.h
index d41c23f5..ec51c188 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -3,12 +3,12 @@
  *
  * The processing pipeline for one image
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2010-2016 Thomas White <taw@physics.org>
- *   2014      Valerio Mariani
+ *   2014-2017 Valerio Mariani <valerio.mariani@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -42,6 +42,7 @@ struct index_args;
 
 
 enum {
+        PEAK_PEAKFINDER8,
 	PEAK_ZAEF,
 	PEAK_HDF5,
 	PEAK_CXI,
@@ -61,24 +62,32 @@ struct index_args
 	float min_snr;
 	int check_hdf5_snr;
 	struct detector *det;
-	IndexingMethod *indm;
-	IndexingPrivate **ipriv;
+	IndexingPrivate *ipriv;
 	int peaks;                /* Peak detection method */
 	float tols[4];
 	struct beam_params *beam;
 	char *hdf5_peak_path;
+	int half_pixel_shift;
 	float pk_inn;
 	float pk_mid;
 	float pk_out;
 	float ir_inn;
 	float ir_mid;
 	float ir_out;
-	struct copy_hdf5_field *copyme;
+	int min_res;
+	int max_res;
+	int max_n_peaks;
+	int min_pix_count;
+	int max_pix_count;
+	int local_bg_radius;
+	int min_peaks;
+	struct imagefile_field_list *copyme;
 	int integrate_saturated;
 	int use_saturated;
 	int no_revalidate;
 	int stream_peaks;
 	int stream_refls;
+	int stream_nonhits;
 	IntegrationMethod int_meth;
 	IntDiag int_diag;
 	signed int int_diag_h;
author	Thomas White <taw@physics.org>	2017-07-06 17:22:11 +0200
committer	Thomas White <taw@physics.org>	2017-07-06 17:24:59 +0200
commit	5292f57d4434c7267e8d945871513d742ff42427 (patch)
tree	d460aa5cef5a501516876850ef243cfc27313d5a
parent	48d4a6b8e82cce81222ec58fdfb488ed79ce0bcf (diff)
parent	dc3395900fc3ce0d3961757628ff83ad6456be19 (diff)