Break off integrate_rings_once()

5 files changed, 584 insertions, 272 deletions

diff --git a/Makefile.am b/Makefile.am
index 69baeaee..270ddf21 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -10,7 +10,7 @@ noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/pr_p_gradient_check tests/symmetry_check \
                   tests/centering_check tests/transformation_check \
                   tests/cell_check tests/pr_l_gradient_check \
-                  tests/pr_pl_gradient_check
+                  tests/pr_pl_gradient_check tests/ring_check
 
 MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \
                tests/third_merge_check tests/fourth_merge_check
@@ -25,7 +25,7 @@ PARTIAL_CHECKS = tests/partialator_merge_check_1 \
 TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \
         tests/integration_check \
         tests/symmetry_check tests/centering_check tests/transformation_check \
-        tests/cell_check
+        tests/cell_check tests/ring_check
 
 EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS)
 
@@ -103,6 +103,8 @@ tests_centering_check_SOURCES = tests/centering_check.c
 
 tests_transformation_check_SOURCES = tests/transformation_check.c
 
+tests_ring_check_SOURCES = tests/ring_check.c
+
 tests_cell_check_SOURCES = tests/cell_check.c
 
 INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 1b68cf0b..5bd520e4 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -43,6 +43,7 @@
 #endif
 
 #include "reflist.h"
+#include "reflist-utils.h"
 #include "cell.h"
 #include "crystal.h"
 #include "cell-utils.h"
@@ -155,6 +156,8 @@ enum boxmask_val
 
 struct intcontext
 {
+	IntegrationMethod meth;
+
 	int halfw;
 	int w;
 	enum boxmask_val *bm;  /* Box mask */
@@ -175,6 +178,10 @@ struct intcontext
 	int ir_inn;
 	int ir_mid;
 	int ir_out;
+
+	int n_saturated;
+	int n_implausible;
+	double limit;
 };
 
 
@@ -1099,8 +1106,8 @@ static double calc_sigma(struct intcontext *ic, struct peak_box *bx)
 }
 
 
-static void mean_var_background(struct intcontext *ic, struct peak_box *bx,
-                                double *pmean, double *pvar)
+static void mean_var_area(struct intcontext *ic, struct peak_box *bx,
+                          enum boxmask_val v, double *pmean, double *pvar)
 {
 	int p, q;
 	double sum = 0.0;
@@ -1110,7 +1117,7 @@ static void mean_var_background(struct intcontext *ic, struct peak_box *bx,
 
 	for ( p=0; p<ic->w; p++ ) {
 	for ( q=0; q<ic->w; q++ ) {
-		if ( bx->bm[p + ic->w*q] != BM_BG ) continue;
+		if ( bx->bm[p + ic->w*q] != v ) continue;
 		sum += bx->a*p + bx->b*q + bx->c;
 		n++;
 	}
@@ -1120,7 +1127,7 @@ static void mean_var_background(struct intcontext *ic, struct peak_box *bx,
 	n = 0;
 	for ( p=0; p<ic->w; p++ ) {
 	for ( q=0; q<ic->w; q++ ) {
-		if ( bx->bm[p + ic->w*q] != BM_BG ) continue;
+		if ( bx->bm[p + ic->w*q] != v ) continue;
 		var += pow(bx->a*p + bx->b*q + bx->c - mean, 2.0);
 		n++;
 	}
@@ -1336,6 +1343,9 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 	ic.halfw = ir_out;
 	ic.image = image;
 	ic.k = 1.0/image->lambda;
+	ic.meth = meth;
+	ic.n_saturated = 0;
+	ic.n_implausible = 0;
 	ic.cell = cell;
 	if ( init_intcontext(&ic) ) {
 		ERROR("Failed to initialise integration.\n");
@@ -1483,6 +1493,132 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 }
 
 
+static void integrate_rings_once(Reflection *refl, struct image *image,
+                                 struct intcontext *ic, UnitCell *cell)
+{
+	double pfs, pss;
+	signed int h, k, l;
+	struct peak_box *bx;
+	int pn;
+	struct panel *p;
+	int fid_fs, fid_ss;  /* Center coordinates, rounded,
+				* in overall data block */
+	int cfs, css;  /* Corner coordinates */
+	double intensity;
+	double sigma;
+	int saturated;
+	double one_over_d;
+	int r;
+	double bgmean, sig2_bg, sig2_poisson, aduph;
+	double pkmean, sig2_pk;
+
+	set_redundancy(refl, 0);
+
+	get_detector_pos(refl, &pfs, &pss);
+
+	/* Explicit truncation of digits after the decimal point.
+		* This is actually the correct thing to do here, not
+		* e.g. lrint().  pfs/pss is the position of the spot, measured
+		* in numbers of pixels, from the panel corner (not the center
+		* of the first pixel).  So any coordinate from 2.0 to 2.9999
+		* belongs to pixel index 2. */
+	fid_fs = pfs;
+	fid_ss = pss;
+	pn = find_panel_number(image->det, fid_fs, fid_ss);
+	p = &image->det->panels[pn];
+
+	cfs = (fid_fs-p->min_fs) - ic->halfw;
+	css = (fid_ss-p->min_ss) - ic->halfw;
+
+	bx = add_box(ic);
+	bx->refl = refl;
+	bx->cfs = cfs;
+	bx->css = css;
+	bx->p = p;
+	bx->pn = pn;
+	get_indices(refl, &h, &k, &l);
+	bx->verbose = VERBOSITY;
+
+	if ( ic->meth & INTEGRATION_CENTER ) {
+		r = center_and_check_box(ic, bx, &saturated);
+	} else {
+		r = check_box(ic, bx, &saturated);
+		if ( !r ) {
+			fit_bg(ic, bx);
+			if ( bx->verbose ) show_peak_box(ic, bx);
+		}
+		bx->offs_fs = 0.0;
+		bx->offs_ss = 0.0;
+	}
+	if ( r ) {
+		delete_box(ic, bx);
+		return;
+	}
+
+	if ( saturated ) {
+		ic->n_saturated++;
+		if ( !(ic->meth & INTEGRATION_SATURATED) ) {
+			delete_box(ic, bx);
+			return;
+		}
+	}
+
+	intensity = tentative_intensity(ic, bx);
+	mean_var_area(ic, bx, BM_BG, &bgmean, &sig2_bg);
+
+	mean_var_area(ic, bx, BM_PK, &pkmean, &sig2_pk);
+	if ( bx->verbose ) {
+		STATUS("bg mean, var = %.2f, %.2f\n", bgmean, sig2_bg);
+		STATUS("pk mean, var = %.2f, %.2f\n", pkmean, sig2_pk);
+	}
+
+	aduph = bx->p->adu_per_eV * ph_lambda_to_eV(ic->image->lambda);
+	sig2_poisson = aduph * intensity;
+
+	/* If intensity is within one photon of nothing, set the Poisson
+		* error to be one photon */
+	if ( fabs(intensity / aduph) < 1.0 ) {
+		sig2_poisson = aduph;
+		if ( bx->verbose ) {
+			STATUS("number of photons less than 1\n");
+		}
+	}
+
+	/* If intensity is negative by more than one photon, assume that
+		* the peak is in the background and add a Poisson error of
+		* appropriate size */
+	if ( intensity < -aduph ) {
+		sig2_poisson = -aduph*intensity;
+		if ( bx->verbose ) {
+			STATUS("very negative (%.2f photons)\n",
+				intensity/aduph);
+		}
+	}
+
+	sigma = sqrt(sig2_poisson + bx->m*sig2_bg);
+
+	if ( intensity < -5.0*sigma ) {
+		delete_box(ic, bx);
+		ic->n_implausible++;
+		return;
+	}
+
+	/* Record intensity and set redundancy to 1 */
+	bx->intensity = intensity;
+	set_intensity(refl, intensity);
+	set_esd_intensity(refl, sigma);
+	set_redundancy(refl, 1);
+
+	one_over_d = resolution(cell, h, k, l);
+	if ( one_over_d > ic->limit ) ic->limit = one_over_d;
+
+	/* Update position */
+	pfs += bx->offs_fs;
+	pss += bx->offs_ss;
+	set_detector_pos(refl, 0.0, pfs, pss);
+}
+
+
 static void integrate_rings(IntegrationMethod meth, Crystal *cr,
                             struct image *image,
                             double ir_inn, double ir_mid, double ir_out)
@@ -1492,8 +1628,6 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
 	RefListIterator *iter;
 	UnitCell *cell;
 	struct intcontext ic;
-	int n_saturated = 0;
-	double limit = 0.0;
 
 	list = find_intersections(image, cr);
 	if ( list == NULL ) return;
@@ -1505,10 +1639,13 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
 	ic.halfw = ir_out;
 	ic.image = image;
 	ic.k = 1.0/image->lambda;
+	ic.n_saturated = 0;
+	ic.n_implausible = 0;
 	ic.cell = cell;
 	ic.ir_inn = ir_inn;
 	ic.ir_mid = ir_mid;
 	ic.ir_out = ir_out;
+	ic.limit = 0.0;
 	if ( init_intcontext(&ic) ) {
 		ERROR("Failed to initialise integration.\n");
 		return;
@@ -1519,101 +1656,131 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double pfs, pss;
-		signed int h, k, l;
-		struct peak_box *bx;
-		int pn;
-		struct panel *p;
-		int fid_fs, fid_ss;  /* Center coordinates, rounded,
-		                      * in overall data block */
-		int cfs, css;  /* Corner coordinates */
-		double intensity;
-		double sigma;
-		int saturated;
-		double one_over_d;
-		int r;
-		double bgmean, sig2_bg, sig2_poisson, aduph;
+		integrate_rings_once(refl, image, &ic, cell);
+	}
 
-		set_redundancy(refl, 0);
+	refine_rigid_groups(&ic);
 
-		get_detector_pos(refl, &pfs, &pss);
+	free_intcontext(&ic);
 
-		/* Explicit truncation of digits after the decimal point.
-		 * This is actually the correct thing to do here, not
-		 * e.g. lrint().  pfs/pss is the position of the spot, measured
-		 * in numbers of pixels, from the panel corner (not the center
-		 * of the first pixel).  So any coordinate from 2.0 to 2.9999
-		 * belongs to pixel index 2. */
-		fid_fs = pfs;
-		fid_ss = pss;
-		pn = find_panel_number(image->det, fid_fs, fid_ss);
-		p = &image->det->panels[pn];
+	crystal_set_num_saturated_reflections(cr, ic.n_saturated);
+	crystal_set_resolution_limit(cr, ic.limit);
+	crystal_set_reflections(cr, list);
 
-		cfs = (fid_fs-p->min_fs) - ic.halfw;
-		css = (fid_ss-p->min_ss) - ic.halfw;
+	if ( ic.n_implausible ) {
+		STATUS("Warning: %i implausibly negative reflections.\n",
+		       ic.n_implausible);
+	}
+}
 
-		bx = add_box(&ic);
-		bx->refl = refl;
-		bx->cfs = cfs;
-		bx->css = css;
-		bx->p = p;
-		bx->pn = pn;
-		get_indices(refl, &h, &k, &l);
-		bx->verbose = VERBOSITY;
 
-		if ( meth & INTEGRATION_CENTER ) {
-			r = center_and_check_box(&ic, bx, &saturated);
-		} else {
-			r = check_box(&ic, bx, &saturated);
-			if ( !r ) {
-				fit_bg(&ic, bx);
-				if ( bx->verbose ) show_peak_box(&ic, bx);
-			}
-			bx->offs_fs = 0.0;
-			bx->offs_ss = 0.0;
-		}
-		if ( r ) {
-			delete_box(&ic, bx);
-			continue;
-		}
+static void resolution_cutoff(RefList *list, UnitCell *cell)
+{
+	double *rmins;
+	double *rmaxs;
+	double rmin, rmax;
+	double total_vol, vol_per_shell;
+	int i;
+	int *sh_nref;
+	double *sh_snr;
+	Reflection **highest_snr;
+	Reflection *refl;
+	RefListIterator *iter;
 
-		if ( saturated ) {
-			n_saturated++;
-			if ( !(meth & INTEGRATION_SATURATED) ) {
-				delete_box(&ic, bx);
-				continue;
-			}
-		}
+	const int nshells = 100;
+	const double min_snr = 10.0;
+
+	rmins = malloc(nshells*sizeof(double));
+	rmaxs = malloc(nshells*sizeof(double));
+	sh_nref = malloc(nshells*sizeof(int));
+	sh_snr = malloc(nshells*sizeof(double));
+	highest_snr = malloc(nshells*sizeof(Reflection *));
+
+	if ( (rmins==NULL) || (rmaxs==NULL) || (sh_nref==NULL)
+	  || (sh_snr==NULL) || (highest_snr==NULL) )
+	{
+		ERROR("Couldn't allocate memory for resolution shells.\n");
+		return;
+	}
 
-		intensity = tentative_intensity(&ic, bx);
-		mean_var_background(&ic, bx, &bgmean, &sig2_bg);
-		aduph = bx->p->adu_per_eV * ph_lambda_to_eV(ic.image->lambda);
-		sig2_poisson = aduph * intensity;
-		sigma = sqrt(sig2_poisson + bx->m*sig2_bg);
+	resolution_limits(list, cell, &rmin, &rmax);
 
-		/* Record intensity and set redundancy to 1 */
-		bx->intensity = intensity;
-		set_intensity(refl, intensity);
-		set_esd_intensity(refl, sigma);
-		set_redundancy(refl, 1);
+	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
+	vol_per_shell = total_vol / nshells;
+	rmins[0] = rmin;
+	for ( i=1; i<nshells; i++ ) {
 
-		one_over_d = resolution(cell, h, k, l);
-		if ( one_over_d > limit ) limit = one_over_d;
+		double r;
 
-		/* Update position */
-		pfs += bx->offs_fs;
-		pss += bx->offs_ss;
-		set_detector_pos(refl, 0.0, pfs, pss);
+		r = vol_per_shell + pow(rmins[i-1], 3.0);
+		r = pow(r, 1.0/3.0);
+
+		/* Shells of constant volume */
+		rmaxs[i-1] = r;
+		rmins[i] = r;
 
 	}
+	rmaxs[nshells-1] = rmax;
 
-	refine_rigid_groups(&ic);
+	for ( i=0; i<nshells; i++ ) {
+		sh_nref[i] = 0;
+		sh_snr[i] = 0.0;
+	}
 
-	free_intcontext(&ic);
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+		double res, intensity, sigma, snr;
+		int bin, j;
 
-	crystal_set_num_saturated_reflections(cr, n_saturated);
-	crystal_set_resolution_limit(cr, limit);
-	crystal_set_reflections(cr, list);
+		intensity = get_intensity(refl);
+		sigma = get_esd_intensity(refl);
+		snr = intensity / sigma;
+
+		if ( snr < min_snr ) continue;
+
+		get_indices(refl, &h, &k, &l);
+		res = 2.0*resolution(cell, h, k, l);
+
+		bin = -1;
+		for ( j=0; j<nshells; j++ ) {
+			if ( (res>rmins[j]) && (res<=rmaxs[j]) ) {
+				bin = j;
+				break;
+			}
+		}
+
+		/* Allow for slight rounding errors */
+		if ( (bin == -1) && (res <= rmins[0]) ) bin = 0;
+		if ( (bin == -1) && (res >= rmaxs[nshells-1]) ) bin = 0;
+		assert(bin != -1);
+
+		sh_nref[bin]++;
+		sh_snr[bin] += snr;
+	}
+
+	for ( i=nshells-1; i>=0; i-- ) {
+		if ( sh_nref[i] > 2 ) break;
+	}
+	double cen = (rmins[i]+rmaxs[i])/2.0;
+	STATUS("The highest shell is %i\n", i);
+	printf("%10.3f nm^-1 or %.2f A\n", cen/1e9, 10e9/cen);
+
+	FILE *fh;
+	fh = fopen("cutoff.dat", "w");
+	for ( i=0; i<nshells; i++ ) {
+		double cen = (rmins[i]+rmaxs[i])/2.0;
+		double val = sh_nref[i];
+		fprintf(fh, "%10.3f %.2f %.4f\n", cen/1e9, 10e9/cen, val);
+	}
+	fclose(fh);
+
+	free(rmins);
+	free(rmaxs);
+	free(sh_nref);
+	free(sh_snr);
 }
 
 
@@ -1627,24 +1794,29 @@ void integrate_all(struct image *image, IntegrationMethod meth,
 		switch ( meth & INTEGRATION_METHOD_MASK ) {
 
 			case INTEGRATION_NONE :
-			return;
+			break;
 
 			case INTEGRATION_RINGS :
 			integrate_rings(meth, image->crystals[i], image,
 			                ir_inn, ir_mid, ir_out);
-			return;
+			break;
 
 			case INTEGRATION_PROF2D :
 			integrate_prof2d(meth, image->crystals[i], image,
 			                 ir_inn, ir_mid, ir_out);
-			return;
+			break;
 
 			default :
 			ERROR("Unrecognised integration method %i\n", meth);
-			return;
+			break;
 
 		}
 
+		/* Set resolution limit */
+		resolution_cutoff(crystal_get_reflections(image->crystals[i]),
+		                  crystal_get_cell(image->crystals[i]));
+
+
 	}
 }
 
diff --git a/tests/.gitignore b/tests/.gitignore
index 31ce158f..fabb48f3 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -10,4 +10,5 @@ symmetry_check
 centering_check
 transformation_check
 cell_check
+ring_check
 .dirstamp
diff --git a/tests/integration_check.c b/tests/integration_check.c
index ac871581..4b302955 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -1,10 +1,9 @@
 /*
  * integration_check.c
  *
- * Check peak integration
+ * Check reflection integration
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
- * Copyright © 2012 Andrew Martin <andrew.martin@desy.de>
+ * Copyright © 2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -35,126 +34,9 @@
 #include <utils.h>
 #include <beam-parameters.h>
 
-#include "../libcrystfel/src/peaks.c"
+#include "../libcrystfel/src/integration.c"
 
 
-/* The third integration check draws a Poisson background and checks that, on
- * average, it gets subtracted by the background subtraction. */
-static void third_integration_check(struct image *image, int n_trials,
-                                    int *fail)
-{
-	double mean_intensity = 0.0;
-	double mean_bg = 0.0;
-	double mean_max = 0.0;
-	double mean_sigma = 0.0;
-	int i;
-	int fs, ss;
-	int nfail = 0;
-
-	for ( i=0; i<n_trials; i++ ) {
-
-		double intensity, sigma;
-		double fsp, ssp;
-		int r;
-
-		for ( fs=0; fs<image->width; fs++ ) {
-		for ( ss=0; ss<image->height; ss++ ) {
-			image->data[fs+image->width*ss] = poisson_noise(1000.0);
-		}
-		}
-
-		r = integrate_peak(image, 64, 64, &fsp, &ssp,
-		                   &intensity, &sigma, 10.0, 15.0, 17.0,
-		                   NULL, NULL);
-
-		if ( r == 0 ) {
-			mean_intensity += intensity;
-			mean_sigma += sigma;
-		} else {
-			nfail++;
-		}
-
-	}
-	mean_intensity /= n_trials;
-	mean_bg /= n_trials;
-	mean_max /= n_trials;
-	mean_sigma /= n_trials;
-
-	STATUS("  Third check (mean values): intensity = %.2f, sigma = %.2f,"
-	       " integration failed %i/%i times\n",
-	       mean_intensity, mean_sigma, nfail, n_trials);
-
-/* These values are always wrong, because the integration sucks */
-//	if ( fabs(mean_intensity) > 5.0 ) {
-//		ERROR("Mean intensity should be close to zero.\n");
-//		*fail = 1;
-//	}
-//	if ( fabs(mean_intensity) > mean_sigma/10.0 ) {
-//		ERROR("Mean intensity should be much less than mean sigma.\n");
-//		*fail = 1;
-//	}
-}
-
-
-/* The fourth integration check draws a Poisson background and draws a peak on
- * top of it, then checks that the intensity of the peak is correctly recovered
- * accounting for the background. */
-static void fourth_integration_check(struct image *image, int n_trials,
-                                     int *fail)
-{
-	double mean_intensity = 0.0;
-	double mean_sigma = 0.0;
-	int i;
-	int fs, ss;
-	int pcount = 0;
-	int nfail = 0;
-
-	for ( i=0; i<n_trials; i++ ) {
-
-		double intensity, sigma;
-		double fsp, ssp;
-		int r;
-
-		for ( fs=0; fs<image->width; fs++ ) {
-		for ( ss=0; ss<image->height; ss++ ) {
-			int idx = fs+image->width*ss;
-			image->data[idx] = poisson_noise(1000.0);
-			if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
-			image->data[idx] += 1000.0;
-			pcount++;
-		}
-		}
-
-		r = integrate_peak(image, 64, 64, &fsp, &ssp,
-		                   &intensity, &sigma, 10.0, 15.0, 17.0,
-		                   NULL, NULL);
-
-		if ( r == 0 ) {
-			mean_intensity += intensity;
-			mean_sigma += sigma;
-		} else {
-			nfail++;
-		}
-
-	}
-	mean_intensity /= n_trials;
-	mean_sigma /= n_trials;
-	pcount /= n_trials;
-
-	STATUS(" Fourth check (mean values): intensity = %.2f, sigma = %.2f,"
-	       " integration failed %i/%i times\n",
-	       mean_intensity, mean_sigma, nfail, n_trials);
-
-	if ( fabs(mean_intensity - pcount*1000.0) > 4000.0 ) {
-		ERROR("Mean intensity should be close to %f\n", pcount*1000.0);
-		*fail = 1;
-	}
-	if ( fabs(mean_intensity) < mean_sigma ) {
-		ERROR("Mean intensity should be greater than mean sigma.\n");
-		*fail = 1;
-	}
-}
-
 
 int main(int argc, char *argv[])
 {
@@ -164,7 +46,6 @@ int main(int argc, char *argv[])
 	FILE *fh;
 	unsigned int seed;
 	int fail = 0;
-	const int n_trials = 100;
 	int r, npx;
 	double ex;
 
@@ -208,65 +89,7 @@ int main(int argc, char *argv[])
 	image.n_crystals = 0;
 	image.crystals = NULL;
 
-	/* First check: no intensity -> no peak, or very low intensity */
-	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0, NULL, NULL);
-	STATUS("  First check: integrate_peak() returned %i", r);
-	if ( r == 0 ) {
-
-		STATUS(", intensity = %.2f, sigma = %.2f\n", intensity, sigma);
-
-		if ( fabs(intensity) > 0.01 ) {
-			ERROR("Intensity should be very close to zero.\n");
-			fail = 1;
-		}
-
-	} else {
-		STATUS(" (correct)\n");
-	}
-
-	/* Second check: uniform peak gives correct I and low sigma(I) */
-	npx = 0;
-	for ( fs=0; fs<image.width; fs++ ) {
-	for ( ss=0; ss<image.height; ss++ ) {
-		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
-		image.data[fs+image.width*ss] = 1000.0;
-		npx++;
-	}
-	}
-
-	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0, NULL, NULL);
-	if ( r ) {
-		ERROR(" Second check: integrate_peak() returned %i (wrong).\n",
-		      r);
-		fail = 1;
-	} else {
-
-		STATUS(" Second check: intensity = %.2f, sigma = %.2f\n",
-		       intensity, sigma);
-
-		ex = npx*1000.0;
-		if ( within_tolerance(ex, intensity, 1.0) == 0 ) {
-			ERROR("Intensity should be close to %f\n", ex);
-			fail = 1;
-		}
-
-		ex = sqrt(npx*1000.0);
-		if ( within_tolerance(ex, sigma, 1.0) == 0 ) {
-			ERROR("Sigma should be roughly %f.\n", ex);
-			fail = 1;
-		}
-
-	}
-
-	/* Third check: Poisson background should get mostly subtracted */
-	third_integration_check(&image, n_trials, &fail);
-
-	/* Fourth check: peak on Poisson background */
-	fourth_integration_check(&image, n_trials, &fail);
-
-	/* Fifth check: uniform peak on uniform background */
+	/* Uniform peak on uniform background */
 	npx = 0;
 	for ( fs=0; fs<image.width; fs++ ) {
 	for ( ss=0; ss<image.height; ss++ ) {
diff --git a/tests/ring_check.c b/tests/ring_check.c
new file mode 100644
index 00000000..ac871581
--- /dev/null
+++ b/tests/ring_check.c
@@ -0,0 +1,314 @@
+/*
+ * integration_check.c
+ *
+ * Check peak integration
+ *
+ * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012 Andrew Martin <andrew.martin@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 <stdlib.h>
+#include <stdio.h>
+
+#include <image.h>
+#include <utils.h>
+#include <beam-parameters.h>
+
+#include "../libcrystfel/src/peaks.c"
+
+
+/* The third integration check draws a Poisson background and checks that, on
+ * average, it gets subtracted by the background subtraction. */
+static void third_integration_check(struct image *image, int n_trials,
+                                    int *fail)
+{
+	double mean_intensity = 0.0;
+	double mean_bg = 0.0;
+	double mean_max = 0.0;
+	double mean_sigma = 0.0;
+	int i;
+	int fs, ss;
+	int nfail = 0;
+
+	for ( i=0; i<n_trials; i++ ) {
+
+		double intensity, sigma;
+		double fsp, ssp;
+		int r;
+
+		for ( fs=0; fs<image->width; fs++ ) {
+		for ( ss=0; ss<image->height; ss++ ) {
+			image->data[fs+image->width*ss] = poisson_noise(1000.0);
+		}
+		}
+
+		r = integrate_peak(image, 64, 64, &fsp, &ssp,
+		                   &intensity, &sigma, 10.0, 15.0, 17.0,
+		                   NULL, NULL);
+
+		if ( r == 0 ) {
+			mean_intensity += intensity;
+			mean_sigma += sigma;
+		} else {
+			nfail++;
+		}
+
+	}
+	mean_intensity /= n_trials;
+	mean_bg /= n_trials;
+	mean_max /= n_trials;
+	mean_sigma /= n_trials;
+
+	STATUS("  Third check (mean values): intensity = %.2f, sigma = %.2f,"
+	       " integration failed %i/%i times\n",
+	       mean_intensity, mean_sigma, nfail, n_trials);
+
+/* These values are always wrong, because the integration sucks */
+//	if ( fabs(mean_intensity) > 5.0 ) {
+//		ERROR("Mean intensity should be close to zero.\n");
+//		*fail = 1;
+//	}
+//	if ( fabs(mean_intensity) > mean_sigma/10.0 ) {
+//		ERROR("Mean intensity should be much less than mean sigma.\n");
+//		*fail = 1;
+//	}
+}
+
+
+/* The fourth integration check draws a Poisson background and draws a peak on
+ * top of it, then checks that the intensity of the peak is correctly recovered
+ * accounting for the background. */
+static void fourth_integration_check(struct image *image, int n_trials,
+                                     int *fail)
+{
+	double mean_intensity = 0.0;
+	double mean_sigma = 0.0;
+	int i;
+	int fs, ss;
+	int pcount = 0;
+	int nfail = 0;
+
+	for ( i=0; i<n_trials; i++ ) {
+
+		double intensity, sigma;
+		double fsp, ssp;
+		int r;
+
+		for ( fs=0; fs<image->width; fs++ ) {
+		for ( ss=0; ss<image->height; ss++ ) {
+			int idx = fs+image->width*ss;
+			image->data[idx] = poisson_noise(1000.0);
+			if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
+			image->data[idx] += 1000.0;
+			pcount++;
+		}
+		}
+
+		r = integrate_peak(image, 64, 64, &fsp, &ssp,
+		                   &intensity, &sigma, 10.0, 15.0, 17.0,
+		                   NULL, NULL);
+
+		if ( r == 0 ) {
+			mean_intensity += intensity;
+			mean_sigma += sigma;
+		} else {
+			nfail++;
+		}
+
+	}
+	mean_intensity /= n_trials;
+	mean_sigma /= n_trials;
+	pcount /= n_trials;
+
+	STATUS(" Fourth check (mean values): intensity = %.2f, sigma = %.2f,"
+	       " integration failed %i/%i times\n",
+	       mean_intensity, mean_sigma, nfail, n_trials);
+
+	if ( fabs(mean_intensity - pcount*1000.0) > 4000.0 ) {
+		ERROR("Mean intensity should be close to %f\n", pcount*1000.0);
+		*fail = 1;
+	}
+	if ( fabs(mean_intensity) < mean_sigma ) {
+		ERROR("Mean intensity should be greater than mean sigma.\n");
+		*fail = 1;
+	}
+}
+
+
+int main(int argc, char *argv[])
+{
+	struct image image;
+	double fsp, ssp, intensity, sigma;
+	int fs, ss;
+	FILE *fh;
+	unsigned int seed;
+	int fail = 0;
+	const int n_trials = 100;
+	int r, npx;
+	double ex;
+
+	fh = fopen("/dev/urandom", "r");
+	fread(&seed, sizeof(seed), 1, fh);
+	fclose(fh);
+	srand(seed);
+
+	image.data = malloc(128*128*sizeof(float));
+	image.flags = NULL;
+	image.beam = NULL;
+	image.lambda = ph_eV_to_lambda(1000.0);
+
+	image.det = calloc(1, sizeof(struct detector));
+	image.det->n_panels = 1;
+	image.det->panels = calloc(1, sizeof(struct panel));
+
+	image.det->panels[0].min_fs = 0;
+	image.det->panels[0].max_fs = 128;
+	image.det->panels[0].min_ss = 0;
+	image.det->panels[0].max_ss = 128;
+	image.det->panels[0].fsx = 1.0;
+	image.det->panels[0].fsy = 0.0;
+	image.det->panels[0].ssx = 0.0;
+	image.det->panels[0].ssy = 1.0;
+	image.det->panels[0].xfs = 1.0;
+	image.det->panels[0].yfs = 0.0;
+	image.det->panels[0].xss = 0.0;
+	image.det->panels[0].yss = 1.0;
+	image.det->panels[0].cnx = -64.0;
+	image.det->panels[0].cny = -64.0;
+	image.det->panels[0].clen = 1.0;
+	image.det->panels[0].res = 1.0;
+	image.det->panels[0].adu_per_eV = 1.0/1000.0;  /* -> 1 adu per photon */
+	image.det->panels[0].max_adu = +INFINITY;  /* No cutoff */
+
+	image.width = 128;
+	image.height = 128;
+	memset(image.data, 0, 128*128*sizeof(float));
+
+	image.n_crystals = 0;
+	image.crystals = NULL;
+
+	/* First check: no intensity -> no peak, or very low intensity */
+	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
+	                   10.0, 15.0, 17.0, NULL, NULL);
+	STATUS("  First check: integrate_peak() returned %i", r);
+	if ( r == 0 ) {
+
+		STATUS(", intensity = %.2f, sigma = %.2f\n", intensity, sigma);
+
+		if ( fabs(intensity) > 0.01 ) {
+			ERROR("Intensity should be very close to zero.\n");
+			fail = 1;
+		}
+
+	} else {
+		STATUS(" (correct)\n");
+	}
+
+	/* Second check: uniform peak gives correct I and low sigma(I) */
+	npx = 0;
+	for ( fs=0; fs<image.width; fs++ ) {
+	for ( ss=0; ss<image.height; ss++ ) {
+		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
+		image.data[fs+image.width*ss] = 1000.0;
+		npx++;
+	}
+	}
+
+	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
+	                   10.0, 15.0, 17.0, NULL, NULL);
+	if ( r ) {
+		ERROR(" Second check: integrate_peak() returned %i (wrong).\n",
+		      r);
+		fail = 1;
+	} else {
+
+		STATUS(" Second check: intensity = %.2f, sigma = %.2f\n",
+		       intensity, sigma);
+
+		ex = npx*1000.0;
+		if ( within_tolerance(ex, intensity, 1.0) == 0 ) {
+			ERROR("Intensity should be close to %f\n", ex);
+			fail = 1;
+		}
+
+		ex = sqrt(npx*1000.0);
+		if ( within_tolerance(ex, sigma, 1.0) == 0 ) {
+			ERROR("Sigma should be roughly %f.\n", ex);
+			fail = 1;
+		}
+
+	}
+
+	/* Third check: Poisson background should get mostly subtracted */
+	third_integration_check(&image, n_trials, &fail);
+
+	/* Fourth check: peak on Poisson background */
+	fourth_integration_check(&image, n_trials, &fail);
+
+	/* Fifth check: uniform peak on uniform background */
+	npx = 0;
+	for ( fs=0; fs<image.width; fs++ ) {
+	for ( ss=0; ss<image.height; ss++ ) {
+		image.data[fs+image.width*ss] = 1000.0;
+		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
+		image.data[fs+image.width*ss] += 1000.0;
+		npx++;
+	}
+	}
+
+	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
+	                   10.0, 15.0, 17.0, NULL, NULL);
+	if ( r ) {
+		ERROR("   Fifth check: integrate_peak() returned %i (wrong).\n",
+		      r);
+		fail = 1;
+	} else {
+
+		STATUS("  Fifth check: intensity = %.2f, sigma = %.2f\n",
+		       intensity, sigma);
+
+		ex = npx*1000.0;
+		if ( within_tolerance(ex, intensity, 1.0) == 0 ) {
+			ERROR("Intensity should be close to %f\n", ex);
+			fail = 1;
+		}
+
+		ex = sqrt(npx*1000.0);
+		if ( within_tolerance(ex, sigma, 1.0) == 0 ) {
+			ERROR("Sigma should be roughly %f.\n", ex);
+			fail = 1;
+		}
+
+	}
+
+
+	free(image.beam);
+	free(image.det->panels);
+	free(image.det);
+	free(image.data);
+
+	if ( fail ) return 1;
+
+	return 0;
+}