Initial integration stuff

10 files changed, 804 insertions, 285 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 738614ef..550d04e8 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -300,11 +300,6 @@ with old scripts because this option selects the behaviour which is now the
 default.
 
 .PD 0
-.IP \fB--no-bg-sub\fR
-.PD
-Don't subtract local background estimates from integrated intensities.  This is almost never useful, but might help to detect problems from troublesome background noise.
-
-.PD 0
 .IP \fB--use-saturated\fR
 .PD
 Normally, peaks which contain one or more pixels above max_adu (defined in the detector geometry file) will not be used for indexing.  Using this option skips this check, possibly improving the indexing rate if there is a large proportion of saturated peaks.
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index a9bda55b..2812fa52 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -9,7 +9,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
                          src/reflist-utils.c src/filters.c \
                          src/render.c src/index.c src/dirax.c src/mosflm.c \
                          src/cell-utils.c src/integer_matrix.c src/crystal.c \
-                         src/grainspotter.c src/xds.c
+                         src/grainspotter.c src/xds.c src/integration.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/reax.c
@@ -26,7 +26,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \
                                  src/filters.h src/dirax.h src/mosflm.h \
                                  src/reax.h src/cell-utils.h \
                                  src/integer_matrix.h src/crystal.h \
-                                 src/grainspotter.h src/xds.h
+                                 src/grainspotter.h src/xds.h src/integration.h
 
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
 AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index efb3f2ce..273ac7f0 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -135,11 +135,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	double cet, cez;
 	double pr;
 	double L;
+	double del;
 
 	/* Don't predict 000 */
 	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
 
 	pr = crystal_get_profile_radius(cryst);
+	del = image->div + crystal_get_mosaicity(cryst);
 
 	/* "low" gives the largest Ewald sphere (wavelength short => k large)
 	 * "high" gives the smallest Ewald sphere (wavelength long => k small)
@@ -152,12 +154,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 	tl = sqrt(xl*xl + yl*yl);
 
-	cet = -sin(image->div/2.0) * khigh;
-	cez = -cos(image->div/2.0) * khigh;
+	cet = -sin(del/2.0) * khigh;
+	cez = -cos(del/2.0) * khigh;
 	rhigh = khigh - distance(cet, cez, tl, zl);  /* Loss of precision */
 
-	cet =  sin(image->div/2.0) * klow;
-	cez = -cos(image->div/2.0) * klow;
+	cet =  sin(del/2.0) * klow;
+	cez = -cos(del/2.0) * klow;
 	rlow = klow - distance(cet, cez, tl, zl);  /* Loss of precision */
 
 	if ( (signbit(rlow) == signbit(rhigh))
@@ -168,7 +170,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		ERROR("Reflection with rlow < rhigh!\n");
 		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
 		      h, k, l, rlow, rhigh);
-		ERROR("div = %e\n", image->div);
+		ERROR("div + m = %e\n", del);
 		return NULL;
 	}
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
new file mode 100644
index 00000000..d6860fec
--- /dev/null
+++ b/libcrystfel/src/integration.c
@@ -0,0 +1,672 @@
+/*
+ * integration.c
+ *
+ * Integration of intensities
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * 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 <stdlib.h>
+#include <assert.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_eigen.h>
+
+#include "reflist.h"
+#include "cell.h"
+#include "crystal.h"
+#include "cell-utils.h"
+#include "geometry.h"
+#include "image.h"
+#include "peaks.h"
+#include "integration.h"
+
+
+struct integr_ind
+{
+	double res;
+	Reflection *refl;
+};
+
+
+static int compare_resolution(const void *av, const void *bv)
+{
+	const struct integr_ind *a = av;
+	const struct integr_ind *b = bv;
+
+	return a->res > b->res;
+}
+
+
+static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell,
+                                           int *np)
+{
+	struct integr_ind *il;
+	Reflection *refl;
+	RefListIterator *iter;
+	int i, n;
+
+	n = num_reflections(list);
+	*np = 0;  /* For now */
+
+	if ( n == 0 ) return NULL;
+
+	il = calloc(n, sizeof(struct integr_ind));
+	if ( il == NULL ) return NULL;
+
+	i = 0;
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+		double res;
+
+		if ( get_redundancy(refl) == 0 ) continue;
+
+		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
+
+		il[i].res = res;
+		il[i].refl = refl;
+
+		i++;
+		assert(i <= n);
+	}
+
+	qsort(il, n, sizeof(struct integr_ind), compare_resolution);
+
+	*np = n;
+	return il;
+}
+
+
+static void check_eigen(gsl_vector *e_val)
+{
+	int i;
+	double vmax, vmin;
+	const int n = e_val->size;
+	const double max_condition = 1e6;
+	const int verbose = 0;
+
+	if ( verbose ) STATUS("Eigenvalues:\n");
+	vmin = +INFINITY;
+	vmax = 0.0;
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( verbose ) STATUS("%i: %e\n", i, val);
+		if ( val > vmax ) vmax = val;
+		if ( val < vmin ) vmin = val;
+	}
+
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( val < vmax/max_condition ) {
+			gsl_vector_set(e_val, i, 0.0);
+		}
+	}
+
+	vmin = +INFINITY;
+	vmax = 0.0;
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( val == 0.0 ) continue;
+		if ( val > vmax ) vmax = val;
+		if ( val < vmin ) vmin = val;
+	}
+	if ( verbose ) {
+		STATUS("Condition number: %e / %e = %5.2f\n",
+		       vmax, vmin, vmax/vmin);
+	}
+}
+
+
+static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
+{
+	gsl_matrix *s_vec;
+	gsl_vector *s_val;
+	int err, n;
+	gsl_vector *shifts;
+
+	n = v->size;
+	if ( v->size != M->size1 ) return NULL;
+	if ( v->size != M->size2 ) return NULL;
+
+	s_val = gsl_vector_calloc(n);
+	s_vec = gsl_matrix_calloc(n, n);
+
+	err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val);
+	if ( err ) {
+		ERROR("SVD failed: %s\n", gsl_strerror(err));
+		gsl_matrix_free(s_vec);
+		gsl_vector_free(s_val);
+		return NULL;
+	}
+	/* "M" is now "U" */
+
+	check_eigen(s_val);
+
+	shifts = gsl_vector_calloc(n);
+	err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
+	if ( err ) {
+		ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
+		gsl_matrix_free(s_vec);
+		gsl_vector_free(s_val);
+		gsl_vector_free(shifts);
+		return NULL;
+	}
+
+	gsl_matrix_free(s_vec);
+	gsl_vector_free(s_val);
+
+	return shifts;
+}
+
+
+struct intcontext
+{
+	int halfw;
+	int w;
+	struct image *image;
+	gsl_matrix *bgm;
+};
+
+
+static void addm(gsl_matrix *m, int i, int j, double val)
+{
+	double v = gsl_matrix_get(m, i, j);
+	gsl_matrix_set(m, i, j, v+val);
+}
+
+
+static void addv(gsl_vector *v, int i, double val)
+{
+	double k = gsl_vector_get(v, i);
+	gsl_vector_set(v, i, k+val);
+}
+
+
+static double boxi(struct intcontext *ic,
+                   int fid_fs, int fid_ss, int p, int q)
+{
+	int fs, ss;
+
+	fs = fid_fs + p;
+	ss = fid_ss + q;
+
+	/* FIXME: Wrong panel? */
+	return ic->image->data[fs + ic->image->width*ss];
+}
+
+
+static void fit_bg(struct intcontext *ic,
+                   int fid_fs, int fid_ss,
+                   double *pa, double *pb, double *pc)
+{
+	int p, q;
+	gsl_vector *v;
+	gsl_vector *ans;
+
+	v = gsl_vector_calloc(3);
+
+	for ( p=0; p<ic->w; p++ ) {
+	for ( q=0; q<ic->w; q++ ) {
+		double bi;
+		bi = boxi(ic, fid_fs, fid_ss, p, q);
+
+		addv(v, 0, bi*p);
+		addv(v, 1, bi*q);
+		addv(v, 2, bi);
+	}
+	}
+
+	ans = solve_svd(v, ic->bgm);
+	gsl_vector_free(v);
+
+	*pa = gsl_vector_get(ans, 0);
+	*pb = gsl_vector_get(ans, 1);
+	*pc = gsl_vector_get(ans, 2);
+
+	gsl_vector_free(ans);
+}
+
+
+static void init_intcontext(struct intcontext *ic)
+{
+	int p, q;
+
+	ic->w = 2*ic->halfw + 1;
+
+	ic->bgm = gsl_matrix_calloc(3, 3);
+	if ( ic->bgm == NULL ) {
+		ERROR("Failed to initialise matrix.\n");
+		return;
+	}
+
+	for ( p=0; p<ic->w; p++ ) {
+	for ( q=0; q<ic->w; q++ ) {
+		addm(ic->bgm, 0, 0, p*p);
+		addm(ic->bgm, 0, 1, p*q);
+		addm(ic->bgm, 0, 2, p);
+		addm(ic->bgm, 1, 0, p*q);
+		addm(ic->bgm, 1, 1, q*q);
+		addm(ic->bgm, 1, 2, q);
+		addm(ic->bgm, 2, 0, p);
+		addm(ic->bgm, 2, 1, q);
+		addm(ic->bgm, 2, 2, 1);
+	}
+	}
+}
+
+
+static void measure_all_intensities(RefList *list, struct image *image)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	struct intcontext ic;
+
+	ic.halfw = 4;
+	ic.image = image;
+	init_intcontext(&ic);
+
+	for ( refl = first_refl(list, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+	{
+		double pfs, pss;
+		int fid_fs, fid_ss;
+		double a, b, c;
+
+		get_detector_pos(refl, &pfs, &pss);
+		fid_fs = lrint(pfs);
+		fid_ss = lrint(pss);
+		fit_bg(&ic, fid_fs, fid_ss, &a, &b, &c);
+
+		//set_intensity(refl, intensity);
+	}
+}
+
+
+static void estimate_mosaicity(Crystal *cr, struct image *image)
+{
+	int msteps = 50;
+	int i;
+	const double mest = crystal_get_mosaicity(cr);
+	const double mmax = 2.0 * mest;
+	RefList *list;
+
+	STATUS("Initial estimate: m = %f\n", mest);
+
+	crystal_set_mosaicity(cr, mmax);
+	list = find_intersections(image, cr);
+	crystal_set_reflections(cr, list);
+	measure_all_intensities(list, image);
+
+	for ( i=1; i<=msteps; i++ ) {
+
+		/* "m" varies from just over zero up to 2x the given estimate */
+		Reflection *refl;
+		RefListIterator *iter;
+		const double m = mmax*((double)i/msteps);
+		int n_gained = 0;
+		int n_lost = 0;
+		double i_gained = 0.0;
+		double i_lost = 0.0;
+
+		crystal_set_mosaicity(cr, m);
+		update_partialities(cr, PMODEL_SPHERE);
+
+		for ( refl = first_refl(list, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			if ( get_redundancy(refl) == 0 ) {
+				if ( get_temp1(refl) > 0.0 ) {
+					i_lost += get_intensity(refl);
+					n_lost++;
+				}
+				set_temp1(refl, -1.0);
+			} else if ( get_temp1(refl) < 0.0 ) {
+				i_gained += get_intensity(refl);
+				n_gained++;
+				set_temp1(refl, 1.0);
+			}
+		}
+
+		if ( i > 1 ) {
+			STATUS("%.2e %10.2f %4i %10.2f %4i %10.2f\n", m,
+			       i_gained, n_gained, i_lost, n_lost,
+		               i_gained - i_lost);
+		}
+
+	}
+}
+
+
+static void estimate_resolution(RefList *reflections, Crystal *cr,
+                                struct image *image)
+{
+	struct integr_ind *il;
+	int n, i;
+	int score = 1000;  /* FIXME */
+	int cutoff = 0;
+	double limit = 0.0;
+
+	if ( num_reflections(reflections) == 0 ) return;
+
+	il = sort_reflections(reflections, crystal_get_cell(cr), &n);
+	if ( il == NULL ) {
+		ERROR("Couldn't sort reflections\n");
+		return;
+	}
+
+	for ( i=0; i<n; i++ ) {
+
+		double intensity, sigma, snr;
+		Reflection *refl;
+
+		refl = il[i].refl;
+		intensity = get_intensity(refl);
+		sigma = get_esd_intensity(refl);
+
+		/* I/sigma(I) cutoff
+		 * Rejects reflections below --min-integration-snr, or if the
+		 * SNR is clearly silly.  Silly indicates that the intensity
+		 * was zero. */
+		snr = fabs(intensity)/sigma;
+
+		/* Record intensity and set redundancy to 1 on success */
+		if ( !cutoff ) {
+			set_redundancy(refl, 1);
+		} else {
+			set_redundancy(refl, 0);
+		}
+
+		if ( snr > 3.0 ) {
+			score++;
+		} else {
+			score--;
+		}
+
+		//STATUS("%5.2f A, %5.2f, %i\n", 1e10/il[i].res, snr, score);
+		if ( score == 0 ) {
+			limit = il[i].res;
+			cutoff = 1;
+		}
+
+	}
+
+	crystal_set_resolution_limit(cr, limit);
+
+	free(il);
+}
+
+
+static void integrate_refine(Crystal *cr, struct image *image, int use_closer,
+                             double min_snr,
+                             double ir_inn, double ir_mid, double ir_out,
+                             int integrate_saturated, int **bgMasks)
+{
+	RefList *reflections;
+
+	/* Create initial list of reflections with nominal parameters */
+	reflections = find_intersections(image, cr);
+	measure_all_intensities(reflections, image);
+
+	/* Find resolution limit of pattern using this list */
+	estimate_resolution(reflections, cr, image);
+
+	reflist_free(reflections);
+
+	STATUS("Initial resolution estimate = %.2f nm^-1 or %.2f A\n",
+	       crystal_get_resolution_limit(cr)/1e9,
+	       1e9 / crystal_get_resolution_limit(cr));
+
+	/* Estimate the mosaicity of the crystal using this resolution limit */
+	estimate_mosaicity(cr, image);
+
+	/* Create new list of reflections with refined mosaicity */
+	reflections = find_intersections(image, cr);
+	measure_all_intensities(reflections, image);
+
+	estimate_resolution(reflections, cr, image);
+}
+
+
+static void integrate_rings(Crystal *cr, struct image *image, int use_closer,
+                            double min_snr,
+                            double ir_inn, double ir_mid, double ir_out,
+                            int integrate_saturated, int **bgMasks)
+{
+	RefList *list;
+	Reflection *refl;
+	RefListIterator *iter;
+	UnitCell *cell;
+	int n_saturated = 0;
+	double limit = 0.0;
+
+	list = find_intersections(image, cr);
+	if ( list == NULL ) return;
+
+	if ( num_reflections(list) == 0 ) return;
+
+	cell = crystal_get_cell(cr);
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double fs, ss, intensity;
+		double d;
+		int idx;
+		double sigma, snr;
+		double pfs, pss;
+		int r;
+		signed int h, k, l;
+		struct panel *p;
+		int pnum, j, found;
+		int saturated;
+		double one_over_d;
+
+		get_detector_pos(refl, &pfs, &pss);
+		get_indices(refl, &h, &k, &l);
+
+		/* Is there a really close feature which was detected? */
+		if ( use_closer ) {
+
+			struct imagefeature *f;
+
+			if ( image->features != NULL ) {
+				f = image_feature_closest(image->features,
+					                  pfs, pss, &d, &idx);
+			} else {
+				f = NULL;
+			}
+
+			/* FIXME: Horrible hardcoded value */
+			if ( (f != NULL) && (d < 10.0) ) {
+
+				double exe;
+
+				exe = get_excitation_error(refl);
+
+				pfs = f->fs;
+				pss = f->ss;
+
+				set_detector_pos(refl, exe, pfs, pss);
+
+			}
+
+		}
+
+		p = find_panel(image->det, pfs, pss);
+		if ( p == NULL ) continue;  /* Next peak */
+		found = 0;
+		for ( j=0; j<image->det->n_panels; j++ ) {
+			if ( &image->det->panels[j] == p ) {
+				pnum = j;
+				found = 1;
+				break;
+			}
+		}
+		if ( !found ) {
+			ERROR("Couldn't find panel %p in list.\n", p);
+			return;
+		}
+
+		r = integrate_peak(image, pfs, pss, &fs, &ss,
+		                   &intensity, &sigma, ir_inn, ir_mid, ir_out,
+		                   bgMasks[pnum], &saturated);
+
+		if ( !r && saturated ) {
+			n_saturated++;
+			if ( !integrate_saturated ) r = 1;
+		}
+
+		/* I/sigma(I) cutoff
+		 * Rejects reflections below --min-integration-snr, or if the
+		 * SNR is clearly silly.  Silly indicates that the intensity
+		 * was zero. */
+		snr = fabs(intensity)/sigma;
+		if ( !r && (isnan(snr) || (snr < min_snr)) ) {
+			r = 1;
+		}
+
+		/* Record intensity and set redundancy to 1 on success */
+		if ( !r ) {
+			set_intensity(refl, intensity);
+			set_esd_intensity(refl, sigma);
+			set_redundancy(refl, 1);
+		} else {
+			set_redundancy(refl, 0);
+		}
+
+		one_over_d = resolution(cell, h, k, l);
+		if ( one_over_d > limit ) limit = one_over_d;
+
+	}
+
+	crystal_set_num_saturated_reflections(cr, n_saturated);
+	crystal_set_resolution_limit(cr, limit);
+	crystal_set_reflections(cr, list);
+}
+
+
+void integrate_all(struct image *image, IntegrationMethod meth,
+                   int use_closer, double min_snr,
+                   double ir_inn, double ir_mid, double ir_out,
+                   int integrate_saturated)
+{
+	int i;
+	int **bgMasks;
+
+	/* Make background masks for all panels */
+	bgMasks = calloc(image->det->n_panels, sizeof(int *));
+	if ( bgMasks == NULL ) {
+		ERROR("Couldn't create list of background masks.\n");
+		return;
+	}
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		int *mask;
+		mask = make_BgMask(image, &image->det->panels[i], ir_inn);
+		if ( mask == NULL ) {
+			ERROR("Couldn't create background mask.\n");
+			return;
+		}
+		bgMasks[i] = mask;
+	}
+
+	for ( i=0; i<image->n_crystals; i++ ) {
+
+		switch ( meth & INTEGRATION_METHOD_MASK ) {
+
+			case INTEGRATION_NONE :
+			return;
+
+			case INTEGRATION_RINGS :
+			integrate_rings(image->crystals[i], image, use_closer,
+			                min_snr, ir_inn, ir_mid, ir_out,
+				        integrate_saturated, bgMasks);
+			return;
+
+			case INTEGRATION_REFINE :
+			integrate_refine(image->crystals[i], image, use_closer,
+			                 min_snr, ir_inn, ir_mid, ir_out,
+				         integrate_saturated, bgMasks);
+			return;
+
+			default :
+			ERROR("Unrecognised integration method %i\n", meth);
+			return;
+
+		}
+
+	}
+
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		free(bgMasks[i]);
+	}
+	free(bgMasks);
+}
+
+
+IntegrationMethod integration_method(const char *str, int *err)
+{
+	int n, i;
+	char **methods;
+	IntegrationMethod meth = INTEGRATION_NONE;
+
+	if ( err != NULL ) *err = 0;
+	n = assplode(str, ",-", &methods, ASSPLODE_NONE);
+
+	for ( i=0; i<n; i++ ) {
+
+		if ( strcmp(methods[i], "rings") == 0) {
+			meth = INTEGRATION_DEFAULTS_RINGS;
+
+		} else if ( strcmp(methods[i], "refine") == 0) {
+			meth = INTEGRATION_DEFAULTS_REFINE;
+
+		} else if ( strcmp(methods[i], "none") == 0) {
+			return INTEGRATION_NONE;
+
+		} else {
+			ERROR("Bad integration method: '%s'\n", str);
+			if ( err != NULL ) *err = 1;
+			return INTEGRATION_NONE;
+		}
+
+		free(methods[i]);
+
+	}
+	free(methods);
+
+	return meth;
+
+}
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
new file mode 100644
index 00000000..12e47ff5
--- /dev/null
+++ b/libcrystfel/src/integration.h
@@ -0,0 +1,75 @@
+/*
+ * integration.h
+ *
+ * Integration of intensities
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * 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 INTEGRATION_H
+#define INTEGRATION_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#define INTEGRATION_DEFAULTS_RINGS (INTEGRATION_RINGS | INTEGRATION_SATURATED)
+#define INTEGRATION_DEFAULTS_REFINE (INTEGRATION_REFINE | INTEGRATION_SATURATED)
+
+/**
+ * IntegrationMethod:
+ * @INTEGRATION_NONE: No integration at all
+ * @INTEGRATION_RINGS: Predict reflections and integrate them all
+ * @INTEGRATION_REFINE: As @INTEGRATION_RINGS, but
+ *
+ * @INTEGRATION_SATURATED: Integrate saturated reflections
+ *
+ * An enumeration of all the available integration methods.
+ **/
+typedef enum {
+
+	INTEGRATION_NONE   = 0,
+
+	/* The core integration methods themselves */
+	INTEGRATION_RINGS  = 1,
+	INTEGRATION_REFINE = 2,
+
+	/* Bits at the top of the IntegrationMethod are flags which modify the
+	 * behaviour of the integration. */
+	INTEGRATION_SATURATED = 256,
+	INTEGRATION_CLOSER = 512,
+
+} IntegrationMethod;
+
+/* This defines the bits in "IntegrationMethod" which are used to represent the
+ * core of the integration method */
+#define INTEGRATION_METHOD_MASK (0xff)
+
+extern IntegrationMethod integration_method(const char *t, int *err);
+
+extern void integrate_all(struct image *image, IntegrationMethod meth,
+	                  int use_closer, double min_snr,
+	                  double ir_inn, double ir_mid, double ir_out,
+	                  int integrate_saturated);
+
+#endif	/* INTEGRATION_H */
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 4b8a8ee8..05c564ba 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -52,6 +52,7 @@
 #include "reflist-utils.h"
 #include "beam-parameters.h"
 #include "cell-utils.h"
+#include "geometry.h"
 
 
 /* Degree of polarisation of X-ray beam */
@@ -200,7 +201,7 @@ static void add_crystal_to_mask(struct image *image, struct panel *p,
 
 
 /* cfs, css relative to panel origin */
-static int *make_BgMask(struct image *image, struct panel *p, double ir_inn)
+int *make_BgMask(struct image *image, struct panel *p, double ir_inn)
 {
 	int *mask;
 	int w, h;
@@ -224,11 +225,11 @@ static int *make_BgMask(struct image *image, struct panel *p, double ir_inn)
 
 /* Returns non-zero if peak has been vetoed.
  * i.e. don't use result if return value is not zero. */
-static int integrate_peak(struct image *image, int cfs, int css,
-                          double *pfs, double *pss,
-                          double *intensity, double *sigma,
-                          double ir_inn, double ir_mid, double ir_out,
-                          int *bgPkMask, int *saturated)
+int integrate_peak(struct image *image, int cfs, int css,
+                   double *pfs, double *pss,
+                   double *intensity, double *sigma,
+                   double ir_inn, double ir_mid, double ir_out,
+                   int *bgPkMask, int *saturated)
 {
 	signed int dfs, dss;
 	double lim_sq, out_lim_sq, mid_lim_sq;
@@ -653,233 +654,6 @@ int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst)
 }
 
 
-struct integr_ind
-{
-	double res;
-	Reflection *refl;
-};
-
-
-static int compare_resolution(const void *av, const void *bv)
-{
-	const struct integr_ind *a = av;
-	const struct integr_ind *b = bv;
-
-	return a->res > b->res;
-}
-
-
-static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell,
-                                           int *np)
-{
-	struct integr_ind *il;
-	Reflection *refl;
-	RefListIterator *iter;
-	int i, n;
-
-	n = num_reflections(list);
-	*np = 0;  /* For now */
-
-	if ( n == 0 ) return NULL;
-
-	il = calloc(n, sizeof(struct integr_ind));
-	if ( il == NULL ) return NULL;
-
-	i = 0;
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l;
-		double res;
-
-		get_indices(refl, &h, &k, &l);
-		res = resolution(cell, h, k, l);
-
-		il[i].res = res;
-		il[i].refl = refl;
-
-		i++;
-		assert(i <= n);
-	}
-
-	qsort(il, n, sizeof(struct integr_ind), compare_resolution);
-
-	*np = n;
-	return il;
-}
-
-
-static void integrate_crystal(Crystal *cr, struct image *image, int use_closer,
-                              int bgsub, double min_snr,
-                              double ir_inn, double ir_mid, double ir_out,
-                              int integrate_saturated, int **bgMasks,
-                              int res_cutoff)
-{
-	RefList *reflections;
-	struct integr_ind *il;
-	int n, i;
-	double av = 0.0;
-	int first = 1;
-	int n_saturated = 0;
-
-	reflections = crystal_get_reflections(cr);
-
-	if ( num_reflections(reflections) == 0 ) return;
-
-	il = sort_reflections(reflections, crystal_get_cell(cr), &n);
-	if ( il == NULL ) {
-		ERROR("Couldn't sort reflections\n");
-		return;
-	}
-
-	for ( i=0; i<n; i++ ) {
-
-		double fs, ss, intensity;
-		double d;
-		int idx;
-		double sigma, snr;
-		double pfs, pss;
-		int r;
-		Reflection *refl;
-		signed int h, k, l;
-		struct panel *p;
-		int pnum, j, found;
-		int saturated;
-
-		refl = il[i].refl;
-
-		get_detector_pos(refl, &pfs, &pss);
-		get_indices(refl, &h, &k, &l);
-
-		/* Is there a really close feature which was detected? */
-		if ( use_closer ) {
-
-			struct imagefeature *f;
-
-			if ( image->features != NULL ) {
-				f = image_feature_closest(image->features,
-					                  pfs, pss, &d, &idx);
-			} else {
-				f = NULL;
-			}
-
-			/* FIXME: Horrible hardcoded value */
-			if ( (f != NULL) && (d < 10.0) ) {
-
-				double exe;
-
-				exe = get_excitation_error(refl);
-
-				pfs = f->fs;
-				pss = f->ss;
-
-				set_detector_pos(refl, exe, pfs, pss);
-
-			}
-
-		}
-
-		p = find_panel(image->det, pfs, pss);
-		if ( p == NULL ) continue;  /* Next peak */
-		found = 0;
-		for ( j=0; j<image->det->n_panels; j++ ) {
-			if ( &image->det->panels[j] == p ) {
-				pnum = j;
-				found = 1;
-				break;
-			}
-		}
-		if ( !found ) {
-			ERROR("Couldn't find panel %p in list.\n", p);
-			return;
-		}
-
-		r = integrate_peak(image, pfs, pss, &fs, &ss,
-		                   &intensity, &sigma, ir_inn, ir_mid, ir_out,
-		                   bgMasks[pnum], &saturated);
-
-		if ( !r && saturated ) {
-			n_saturated++;
-			if ( !integrate_saturated ) r = 1;
-		}
-
-		/* I/sigma(I) cutoff
-		 * Rejects reflections below --min-integration-snr, or if the
-		 * SNR is clearly silly.  Silly indicates that the intensity
-		 * was zero. */
-		snr = fabs(intensity)/sigma;
-		if ( isnan(snr) || (snr < min_snr) ) {
-			r = 1;
-		}
-
-		/* Record intensity and set redundancy to 1 on success */
-		if ( r == 0 ) {
-			set_intensity(refl, intensity);
-			set_esd_intensity(refl, sigma);
-			set_redundancy(refl, 1);
-		} else {
-			set_redundancy(refl, 0);
-		}
-
-		if ( snr > 1.0 ) {
-			if ( first ) {
-				av = snr;
-				first = 0;
-			} else {
-				av = av + 0.1*(snr - av);
-			}
-			//STATUS("%5.2f A, %5.2f, av %5.2f\n",
-			//       1e10/il[i].res, snr, av);
-			if ( res_cutoff && (av < 1.0) ) break;
-		}
-	}
-
-	crystal_set_num_saturated_reflections(cr, n_saturated);
-	crystal_set_resolution_limit(cr, 0.0);
-
-	free(il);
-}
-
-
-/* Integrate the list of predicted reflections in "image" */
-void integrate_reflections(struct image *image, int use_closer, int bgsub,
-                           double min_snr,
-                           double ir_inn, double ir_mid, double ir_out,
-                           int integrate_saturated, int res_cutoff)
-{
-	int i;
-	int **bgMasks;
-
-	/* Make background masks for all panels */
-	bgMasks = calloc(image->det->n_panels, sizeof(int *));
-	if ( bgMasks == NULL ) {
-		ERROR("Couldn't create list of background masks.\n");
-		return;
-	}
-	for ( i=0; i<image->det->n_panels; i++ ) {
-		int *mask;
-		mask = make_BgMask(image, &image->det->panels[i], ir_inn);
-		if ( mask == NULL ) {
-			ERROR("Couldn't create background mask.\n");
-			return;
-		}
-		bgMasks[i] = mask;
-	}
-
-	for ( i=0; i<image->n_crystals; i++ ) {
-		integrate_crystal(image->crystals[i], image, use_closer,
-		                  bgsub, min_snr, ir_inn, ir_mid, ir_out,
-		                  integrate_saturated, bgMasks, res_cutoff);
-	}
-
-	for ( i=0; i<image->det->n_panels; i++ ) {
-		free(bgMasks[i]);
-	}
-	free(bgMasks);
-}
-
-
 void validate_peaks(struct image *image, double min_snr,
                     int ir_inn, int ir_mid, int ir_out, int use_saturated)
 {
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index f788bae5..b94ebee5 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -37,16 +37,13 @@
 
 #include "reflist.h"
 
+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);
 
-extern void integrate_reflections(struct image *image,
-                                  int use_closer, int bgsub, double min_snr,
-                                  double ir_inn, double ir_mid, double ir_out,
-                                  int integrate_saturated, int res_cutoff);
-
 extern int peak_sanity_check(struct image *image, Crystal **crystals,
                              int n_cryst);
 
@@ -54,4 +51,10 @@ extern void validate_peaks(struct image *image, double min_snr,
                            int ir_inn, int ir_mid, int ir_out,
                            int use_saturated);
 
+extern int integrate_peak(struct image *image, int cfs, int css,
+                          double *pfs, double *pss,
+                          double *intensity, double *sigma,
+                          double ir_inn, double ir_mid, double ir_out,
+                          int *bgPkMask, int *saturated);
+
 #endif	/* PEAKS_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 208c0ed0..7d26c15d 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -66,6 +66,7 @@
 #include "stream.h"
 #include "reflist-utils.h"
 #include "cell-utils.h"
+#include "integration.h"
 
 #include "im-sandbox.h"
 
@@ -100,6 +101,7 @@ static void show_help(const char *s)
 "                           hdf5  : Get from a table in HDF5 file.\n"
 "     --hdf5-peaks=<p>     Find peaks table in HDF5 file here.\n"
 "                           Default: /processing/hitfinder/peakinfo\n"
+"     --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"
@@ -139,17 +141,12 @@ static void show_help(const char *s)
 "     --no-check-prefix    Don't attempt to correct the --prefix.\n"
 "     --closer-peak        Don't integrate from the location of a nearby peak\n"
 "                           instead of the predicted spot.  Don't use.\n"
-"     --no-bg-sub          Don't subtract local background estimates from\n"
-"                           integrated intensities.\n"
 "     --use-saturated      During the initial peak search, don't reject\n"
 "                           peaks which contain pixels above max_adu.\n"
 "     --integrate-saturated During the final integration stage, don't reject\n"
 "                           peaks which contain pixels above max_adu.\n"
 "     --no-revalidate      Don't re-integrate and check HDF5 peaks for\n"
 "                           validity.\n"
-"     --integrate-found    Skip the spot prediction step, and just integrate\n"
-"                           the intensities of the spots found by the initial\n"
-"                           peak search.\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"
 );
@@ -179,6 +176,7 @@ int main(int argc, char *argv[])
 	char *use_this_one_instead;
 	struct index_args iargs;
 	char *intrad = NULL;
+	char *int_str = NULL;
 
 	/* Defaults */
 	iargs.cell = NULL;
@@ -186,7 +184,6 @@ int main(int argc, char *argv[])
 	iargs.median_filter = 0;
 	iargs.satcorr = 1;
 	iargs.closer = 0;
-	iargs.bgsub = 1;
 	iargs.tols[0] = 5.0;
 	iargs.tols[1] = 5.0;
 	iargs.tols[2] = 5.0;
@@ -207,10 +204,8 @@ int main(int argc, char *argv[])
 	iargs.use_saturated = 0;
 	iargs.integrate_saturated = 0;
 	iargs.no_revalidate = 0;
-	iargs.integrate_found = 0;
 	iargs.stream_peaks = 1;
 	iargs.stream_refls = 1;
-	iargs.res_cutoff = 0;
 	iargs.copyme = new_copy_hdf5_field_list();
 	if ( iargs.copyme == NULL ) {
 		ERROR("Couldn't allocate HDF5 field list.\n");
@@ -218,6 +213,7 @@ int main(int argc, char *argv[])
 	}
 	iargs.indm = NULL;  /* No default */
 	iargs.ipriv = NULL;  /* No default */
+	iargs.int_meth = integration_method("rings", NULL);
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -242,15 +238,11 @@ int main(int argc, char *argv[])
 		{"no-closer-peak",     0, &iargs.closer,             0},
 		{"closer-peak",        0, &iargs.closer,             1},
 		{"basename",           0, &config_basename,          1},
-		{"bg-sub",             0, &iargs.bgsub,              1},
-		{"no-bg-sub",          0, &iargs.bgsub,              0},
 		{"no-peaks-in-stream", 0, &iargs.stream_peaks,       0},
 		{"no-refls-in-stream", 0, &iargs.stream_refls,       0},
-		{"res-cutoff",         0, &iargs.res_cutoff,         1},
 		{"integrate-saturated",0, &iargs.integrate_saturated,1},
 		{"use-saturated",      0, &iargs.use_saturated,      1},
 		{"no-revalidate",      0, &iargs.no_revalidate,      1},
-		{"integrate-found",    0, &iargs.integrate_found,    1},
 
 		/* Long-only options with arguments */
 		{"peaks",              1, NULL,                2},
@@ -267,6 +259,7 @@ int main(int argc, char *argv[])
 		{"tolerance",          1, NULL,               13},
 		{"int-radius",         1, NULL,               14},
 		{"median-filter",      1, NULL,               15},
+		{"integration",        1, NULL,               16},
 
 		{0, 0, NULL, 0}
 	};
@@ -391,6 +384,10 @@ int main(int argc, char *argv[])
 			iargs.median_filter = atoi(optarg);
 			break;
 
+			case 16 :
+			int_str = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -465,6 +462,18 @@ int main(int argc, char *argv[])
 		free(indm_str);
 	}
 
+	if ( int_str != NULL ) {
+
+		int err;
+
+		iargs.int_meth = integration_method(int_str, &err);
+		if ( err ) {
+			ERROR("Invalid integration method '%s'\n", int_str);
+			return 1;
+		}
+		free(int_str);
+	}
+
 	if ( toler != NULL ) {
 		int ttt;
 		ttt = sscanf(toler, "%f,%f,%f,%f",
diff --git a/src/process_image.c b/src/process_image.c
index ef673c84..bf694623 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -46,6 +46,7 @@
 #include "stream.h"
 #include "reflist-utils.h"
 #include "process_image.h"
+#include "integration.h"
 
 
 void process_image(const struct index_args *iargs, struct pattern_args *pargs,
@@ -188,34 +189,21 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	/* Integrate each crystal's diffraction spots */
 	for ( i=0; i<image.n_crystals; i++ ) {
 
-		RefList *reflections;
-
 		/* Set default crystal parameter(s) */
 		crystal_set_profile_radius(image.crystals[i],
 		                           image.beam->profile_radius);
-
-		if ( iargs->integrate_found ) {
-			reflections = select_intersections(&image,
-			                                   image.crystals[i]);
-		} else {
-			reflections = find_intersections(&image,
-			                                 image.crystals[i]);
-		}
-
-		crystal_set_reflections(image.crystals[i], reflections);
+		crystal_set_mosaicity(image.crystals[i], 2e-3);  /* radians */
+		crystal_set_image(image.crystals[i], &image);
 
 	}
 
 	/* Integrate all the crystals at once - need all the crystals so that
 	 * overlaps can be detected. */
-	integrate_reflections(&image, iargs->closer,
-	                              iargs->bgsub,
-	                              iargs->min_int_snr,
-	                              iargs->ir_inn,
-	                              iargs->ir_mid,
-	                              iargs->ir_out,
-	                              iargs->integrate_saturated,
-	                              iargs->res_cutoff);
+	integrate_all(&image, iargs->int_meth,
+	                      iargs->closer,
+	                      iargs->min_int_snr,
+	                      iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+	                      iargs->integrate_saturated);
 
 	write_chunk(st, &image, hdfile,
 	            iargs->stream_peaks, iargs->stream_refls);
diff --git a/src/process_image.h b/src/process_image.h
index 5fe11e33..03041562 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -34,6 +34,9 @@
 #endif
 
 
+#include "integration.h"
+
+
 enum {
 	PEAK_ZAEF,
 	PEAK_HDF5,
@@ -49,7 +52,6 @@ struct index_args
 	int median_filter;
 	int satcorr;
 	int closer;
-	int bgsub;
 	float threshold;
 	float min_gradient;
 	float min_snr;
@@ -69,10 +71,9 @@ struct index_args
 	int integrate_saturated;
 	int use_saturated;
 	int no_revalidate;
-	int integrate_found;
 	int stream_peaks;
 	int stream_refls;
-	int res_cutoff;
+	IntegrationMethod int_meth;
 };