check_hkl: Move "single-list" FoMs into API

Reasons for differences:

1. Resolution shells slightly different
The binning calculation needs to take into account small rounding errors
in the resolution calculation, when not using an explicit resolution
range (--highres).  The old version did this by taking a min/max
resolution range slightly larger than the resolution of the data.  The
new version handles the rounding errors explicitly, so does not need
this.

2. Number of reflections with infinite/invalid I/sigI values halved
The number reported for this count was twice what it should have been,
due to a bug in the old check_hkl.

3. Overall SNR is different
When the above warning message applies, the old version still allowed
the reflections with invalid I/sigI values to contribute to the
denominator of the mean SNR calculation.  The new version does not
include them in the SNR calculation at all.  Note that the reflections
contribute to the other figures of merit unless otherwise stated.

4. Standard deviation of intensity is not calculated
It would've been a lot of work to include this in the new version, and
it's a totally useless number.  If you disagree, please get in touch!

4 files changed, 715 insertions, 488 deletions

diff --git a/libcrystfel/src/fom.c b/libcrystfel/src/fom.c
index 19f35a86..c2de704a 100644
--- a/libcrystfel/src/fom.c
+++ b/libcrystfel/src/fom.c
@@ -43,6 +43,33 @@
 #include "reflist.h"
 #include "reflist-utils.h"
 
+struct fom_context
+{
+	enum fom_type fom;
+	int nshells;
+	int *cts;
+
+	/* For R-factors */
+	double *num;
+	double *den;
+
+	/* For "double" R-factors */
+	double *num2;
+	double *den2;
+
+	/* For CCs */
+	double **vec1;
+	double **vec2;
+	int *n;
+	int nmax;
+
+	/* For "counting" things e.g. d1sig or d2sig */
+	int *n_within;
+
+	long int *n_meas;
+	long int *possible;
+};
+
 enum fom_type fom_type_from_string(const char *s)
 {
 	if ( strcasecmp(s, "r1i") == 0 ) return FOM_R1I;
@@ -78,6 +105,12 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells)
 		fctx->cts[i] = 0;
 	}
 
+	fctx->num = NULL;
+	fctx->den = NULL;
+	fctx->num2 = NULL;
+	fctx->den2 = NULL;
+	fctx->possible = NULL;
+
 	switch ( fctx->fom ) {
 
 		case FOM_RANORSPLIT :
@@ -95,6 +128,10 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells)
 		case FOM_R2 :
 		case FOM_RSPLIT :
 		case FOM_RANO :
+		case FOM_MEAN_INTENSITY :
+		case FOM_STDDEV_INTENSITY :
+		case FOM_SNR :
+		case FOM_REDUNDANCY :
 		fctx->num = malloc(nshells*sizeof(double));
 		fctx->den = malloc(nshells*sizeof(double));
 		if ( (fctx->num == NULL) || (fctx->den == NULL) ) return NULL;
@@ -104,6 +141,15 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells)
 		}
 		break;
 
+		case FOM_COMPLETENESS :
+		/* Uses 'cts' and 'possible' only */
+		break;
+
+		case FOM_NUM_MEASUREMENTS :
+		fctx->n_meas = calloc(nshells, sizeof(long int));
+		if ( fctx->n_meas == NULL ) return NULL;
+		break;
+
 		case FOM_CC :
 		case FOM_CCSTAR :
 		case FOM_CCANO :
@@ -140,37 +186,45 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells)
 }
 
 
-static void add_to_fom(struct fom_context *fctx, double i1, double i2,
-                       double i1bij, double i2bij, double sig1, double sig2,
-                       int bin)
+static int add_to_fom(struct fom_context *fctx,
+                      Reflection *refl1,
+                      Reflection *refl2,
+                      Reflection *refl1bij,
+                      Reflection *refl2bij,
+                      int bin)
 {
-	double f1, f2;
+	double i1, i2, i1bij, i2bij, sig1, sig2;
 	double im, imbij;
+	int bad = 0;
 
 	fctx->cts[bin]++;
 
-	/* Negative intensities have already been weeded out. */
-	f1 = sqrt(i1);
-	f2 = sqrt(i2);
-
 	switch ( fctx->fom ) {
 
 		case FOM_R1I :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		fctx->num[bin] += fabs(i1 - i2);
 		fctx->den[bin] += i1;
 		break;
 
 		case FOM_R1F :
-		fctx->num[bin] += fabs(f1 - f2);
-		fctx->den[bin] += f1;
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
+		fctx->num[bin] += fabs(sqrt(i1) - sqrt(i2));
+		fctx->den[bin] += sqrt(i1);
 		break;
 
 		case FOM_R2 :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		fctx->num[bin] += pow(i1 - i2, 2.0);
 		fctx->den[bin] += pow(i1, 2.0);
 		break;
 
 		case FOM_RSPLIT :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		fctx->num[bin] += fabs(i1 - i2);
 		fctx->den[bin] += i1 + i2;
 		break;
@@ -178,6 +232,8 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
 		case FOM_CC :
 		case FOM_CCSTAR :
 		assert(fctx->n[bin] < fctx->nmax);
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		fctx->vec1[bin][fctx->n[bin]] = i1;
 		fctx->vec2[bin][fctx->n[bin]] = i2;
 		fctx->n[bin]++;
@@ -186,17 +242,27 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
 		case FOM_CCANO :
 		case FOM_CRDANO :
 		assert(fctx->n[bin] < fctx->nmax);
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
+		i1bij = get_intensity(refl1bij);
+		i2bij = get_intensity(refl2bij);
 		fctx->vec1[bin][fctx->n[bin]] = i1 - i1bij;
 		fctx->vec2[bin][fctx->n[bin]] = i2 - i2bij;
 		fctx->n[bin]++;
 		break;
 
 		case FOM_RANORSPLIT :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
 		fctx->num2[bin] += fabs(i1 - i2);
 		fctx->den2[bin] += i1 + i2;
 		/* Intentional fall-through (no break) */
 
 		case FOM_RANO :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
+		i1bij = get_intensity(refl1bij);
+		i2bij = get_intensity(refl2bij);
 		im = (i1 + i2)/2.0;
 		imbij = (i1bij + i2bij)/2.0;
 		fctx->num[bin] += fabs(im - imbij);
@@ -204,22 +270,69 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
 		break;
 
 		case FOM_D1SIG :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
+		sig1 = get_esd_intensity(refl1);
+		sig2 = get_esd_intensity(refl2);
 		if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) {
 			fctx->n_within[bin]++;
 		}
 		break;
 
 		case FOM_D2SIG :
+		i1 = get_intensity(refl1);
+		i2 = get_intensity(refl2);
+		sig1 = get_esd_intensity(refl1);
+		sig2 = get_esd_intensity(refl2);
 		if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) {
 			fctx->n_within[bin]++;
 		}
 		break;
 
+		case FOM_NUM_MEASUREMENTS :
+		fctx->n_meas[bin] += get_redundancy(refl1);
+		break;
+
+		case FOM_REDUNDANCY :
+		fctx->num[bin] += get_redundancy(refl1);
+		fctx->den[bin] += 1.0;
+		break;
+
+		case FOM_SNR :
+		i1 = get_intensity(refl1);
+		sig1 = get_esd_intensity(refl1);
+		if ( isfinite(i1/sig1) ) {
+			fctx->num[bin] += i1/sig1;
+			fctx->den[bin] += 1.0;
+		} else {
+			bad = 1;
+		}
+		break;
+
+		case FOM_MEAN_INTENSITY :
+		i1 = get_intensity(refl1);
+		fctx->num[bin] += i1;
+		fctx->den[bin] += 1.0;
+		break;
+
+		/* FIXME: Delete this FoM? */
+		case FOM_STDDEV_INTENSITY :
+		fctx->num[bin] += 1.0;
+		fctx->den[bin] += 1.0;
+		break;
+
+		case FOM_COMPLETENESS :
+		/* fctx->cts already incremented, as needed.
+		 * Will calculate possible reflections later */
+		break;
+
 	}
+
+	return bad;
 }
 
 
-double fom_overall(struct fom_context *fctx)
+double fom_overall_value(struct fom_context *fctx)
 {
 	double overall_num = INFINITY;
 	double overall_den = 0.0;
@@ -234,6 +347,9 @@ double fom_overall(struct fom_context *fctx)
 	double variance_signal;
 	double variance_error;
 	double cc = INFINITY;
+	long int total_meas = 0;
+	long int overall_cts = 0;
+	long int overall_possible = 0;
 
 	switch ( fctx->fom ) {
 
@@ -242,6 +358,10 @@ double fom_overall(struct fom_context *fctx)
 		case FOM_R2 :
 		case FOM_RSPLIT :
 		case FOM_RANO :
+		case FOM_REDUNDANCY :
+		case FOM_SNR :
+		case FOM_MEAN_INTENSITY :
+		case FOM_STDDEV_INTENSITY :
 		overall_num = 0.0;
 		overall_den = 0.0;
 		for ( i=0; i<fctx->nshells; i++ ) {
@@ -321,14 +441,38 @@ double fom_overall(struct fom_context *fctx)
 		}
 		break;
 
+		case FOM_NUM_MEASUREMENTS :
+		total_meas = 0;
+		for ( i=0; i<fctx->nshells; i++ ) {
+			total_meas += fctx->n_meas[i];
+		}
+		break;
+
+		case FOM_COMPLETENESS :
+		for ( i=0; i<fctx->nshells; i++ ) {
+			overall_cts += fctx->cts[i];
+			overall_possible += fctx->possible[i];
+		}
+		break;
+
 	}
 
 	switch ( fctx->fom ) {
 
 		case FOM_R1I :
 		case FOM_R1F :
+		case FOM_REDUNDANCY :
+		case FOM_SNR :
+		case FOM_MEAN_INTENSITY :
+		case FOM_STDDEV_INTENSITY :
 		return overall_num/overall_den;
 
+		case FOM_COMPLETENESS :
+		return (double)overall_cts / overall_possible;
+
+		case FOM_NUM_MEASUREMENTS :
+		return total_meas;
+
 		case FOM_R2 :
 		return sqrt(overall_num/overall_den);
 
@@ -361,7 +505,7 @@ double fom_overall(struct fom_context *fctx)
 }
 
 
-double fom_shell(struct fom_context *fctx, int i)
+double fom_shell_value(struct fom_context *fctx, int i)
 {
 	double cc;
 	int j;
@@ -374,6 +518,10 @@ double fom_shell(struct fom_context *fctx, int i)
 
 		case FOM_R1I :
 		case FOM_R1F :
+		case FOM_REDUNDANCY :
+		case FOM_SNR :
+		case FOM_MEAN_INTENSITY :
+		case FOM_STDDEV_INTENSITY :
 		return fctx->num[i]/fctx->den[i];
 
 		case FOM_R2 :
@@ -420,6 +568,12 @@ double fom_shell(struct fom_context *fctx, int i)
 		case FOM_D2SIG :
 		return (double)fctx->n_within[i] / fctx->cts[i];
 
+		case FOM_NUM_MEASUREMENTS :
+		return fctx->n_meas[i];
+
+		case FOM_COMPLETENESS :
+		return (double)fctx->cts[i] / fctx->possible[i];
+
 	}
 
 	ERROR("This point is never reached.\n");
@@ -469,7 +623,7 @@ struct fom_shells *fom_make_resolution_shells(double rmin, double rmax,
 }
 
 
-double fom_shell_label(struct fom_shells *s, int i)
+double fom_shell_centre(struct fom_shells *s, int i)
 {
 	return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0;
 }
@@ -605,15 +759,143 @@ static int wilson_scale(RefList *list1, RefList *list2, UnitCell *cell)
 }
 
 
+static int calculate_possible(struct fom_context *fctx,
+                              struct fom_shells *shells,
+                              UnitCell *cell,
+                              const SymOpList *sym)
+{
+	RefList *counted;
+	int hmax, kmax, lmax;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	signed int h, k, l;
+
+	fctx->possible = calloc(fctx->nshells, sizeof(long int));
+	if ( fctx->possible == NULL ) return 1;
+
+	counted = reflist_new();
+	if ( counted == NULL ) {
+		free(fctx->possible);
+		return 1;
+	}
+
+	cell_get_cartesian(cell, &ax, &ay, &az,
+	                         &bx, &by, &bz,
+	                         &cx, &cy, &cz);
+	hmax = shells->rmaxs[fctx->nshells-1] * modulus(ax, ay, az);
+	kmax = shells->rmaxs[fctx->nshells-1] * modulus(bx, by, bz);
+	lmax = shells->rmaxs[fctx->nshells-1] * modulus(cx, cy, cz);
+	for ( h=-hmax; h<=hmax; h++ ) {
+	for ( k=-kmax; k<=kmax; k++ ) {
+	for ( l=-lmax; l<=lmax; l++ ) {
+
+		double d;
+		signed int hs, ks, ls;
+		int bin;
+		int i;
+
+		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;
+
+		bin = -1;
+		for ( i=0; i<fctx->nshells; i++ ) {
+			if ( (d>shells->rmins[i]) && (d<=shells->rmaxs[i]) ) {
+				bin = i;
+				break;
+			}
+		}
+		if ( bin == -1 ) continue;
+
+		if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
+		add_refl(counted, hs, ks, ls);
+
+		fctx->possible[bin]++;
+
+	}
+	}
+	}
+	reflist_free(counted);
+
+	return 0;
+}
+
+
+static int is_anomalous(enum fom_type fom)
+{
+	switch ( fom ) {
+
+		case FOM_CCANO:
+		case FOM_RANO:
+		case FOM_CRDANO:
+		case FOM_RANORSPLIT:
+		return 1;
+
+		case FOM_R1I:
+		case FOM_R1F:
+		case FOM_R2:
+		case FOM_RSPLIT:
+		case FOM_CC:
+		case FOM_CCSTAR:
+		case FOM_D1SIG:
+		case FOM_D2SIG:
+		case FOM_NUM_MEASUREMENTS:
+		case FOM_REDUNDANCY:
+		case FOM_SNR:
+		case FOM_MEAN_INTENSITY:
+		case FOM_STDDEV_INTENSITY:
+		case FOM_COMPLETENESS:
+		return 0;
+	}
+
+	ERROR("This point never reached\n");
+	abort();
+}
+
+
+static int is_single_list(enum fom_type fom)
+{
+	switch ( fom ) {
+
+		case FOM_CCANO:
+		case FOM_RANO:
+		case FOM_CRDANO:
+		case FOM_RANORSPLIT:
+		case FOM_R1I:
+		case FOM_R1F:
+		case FOM_R2:
+		case FOM_RSPLIT:
+		case FOM_CC:
+		case FOM_CCSTAR:
+		case FOM_D1SIG:
+		case FOM_D2SIG:
+		return 0;
+
+		case FOM_NUM_MEASUREMENTS:
+		case FOM_REDUNDANCY:
+		case FOM_SNR:
+		case FOM_MEAN_INTENSITY:
+		case FOM_STDDEV_INTENSITY:
+		case FOM_COMPLETENESS:
+		return 1;
+	}
+
+	ERROR("This point never reached\n");
+	abort();
+}
+
 
 struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell,
                                   struct fom_shells *shells, enum fom_type fom,
-                                  int noscale, SymOpList *sym)
+                                  int noscale, const SymOpList *sym)
 {
 	Reflection *refl1;
 	RefListIterator *iter;
 	struct fom_context *fctx;
-	int n_out;
+	long int n_out = 0;
+	long int n_rej = 0;
 
 	fctx = init_fom(fom, num_reflections(list1), shells->nshells);
 
@@ -622,40 +904,44 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell
 		return NULL;
 	}
 
-	if ( !noscale && wilson_scale(list1, list2, cell) ) {
-		ERROR("Error with scaling.\n");
-		return NULL;
-	}
+	if ( !is_single_list(fom) ) {
+		if ( !noscale && wilson_scale(list1, list2, cell) ) {
+			ERROR("Error with scaling.\n");
+			return NULL;
+		}
 
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		Reflection *refl2;
-		signed int h, k, l;
-		set_flag(refl1, 0);
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		assert(refl2 != NULL);
-		set_flag(refl2, 0);
+		for ( refl1 = first_refl(list1, &iter);
+		      refl1 != NULL;
+		      refl1 = next_refl(refl1, iter) )
+		{
+			Reflection *refl2;
+			signed int h, k, l;
+			set_flag(refl1, 0);
+			get_indices(refl1, &h, &k, &l);
+			refl2 = find_refl(list2, h, k, l);
+			assert(refl2 != NULL);
+			set_flag(refl2, 0);
+		}
 	}
 
-	n_out = 0;
 	for ( refl1 = first_refl(list1, &iter);
 	      refl1 != NULL;
 	      refl1 = next_refl(refl1, iter) )
 	{
 		signed int h, k, l;
 		int bin;
-		double i1, i2;
-		double i1bij, i2bij;
-		double sig1, sig2;
 		Reflection *refl2;
+		Reflection *refl1_bij = NULL;
+		Reflection *refl2_bij = NULL;
 
 		get_indices(refl1, &h, &k, &l);
 
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;
+		if ( is_single_list(fom) ) {
+			refl2 = NULL;
+		} else {
+			refl2 = find_refl(list2, h, k, l);
+			if ( refl2 == NULL ) continue;
+		}
 
 		bin = get_bin(shells, refl1, cell);
 		if ( bin == -1 ) {
@@ -663,17 +949,8 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell
 			continue;
 		}
 
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-		sig1 = get_esd_intensity(refl1);
-		sig2 = get_esd_intensity(refl2);
+		if ( is_anomalous(fom) ) {
 
-		if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO)
-		  || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) )
-		{
-
-			Reflection *refl1_bij = NULL;
-			Reflection *refl2_bij = NULL;
 			signed int hb, kb, lb;
 
 			if ( find_equiv_in_list(list1, -h, -k, -l, sym,
@@ -700,46 +977,59 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell
 			assert(refl1_bij != NULL);
 			assert(refl2_bij != NULL);
 
-			i1bij = get_intensity(refl1_bij);
-			i2bij = get_intensity(refl2_bij);
-
-		} else {
-
-			/* Make it obvious if these get used by mistake */
-			i1bij = +INFINITY;
-			i2bij = +INFINITY;
-
 		}
 
-		add_to_fom(fctx, i1, i2, i1bij, i2bij, sig1, sig2, bin);
+		n_rej += add_to_fom(fctx, refl1, refl2, refl1_bij, refl2_bij, bin);
+
 	}
-	if ( n_out)  {
+	if ( n_out )  {
 		ERROR("WARNING: %i reflection pairs outside range.\n", n_out);
 	}
+	if ( n_rej ) {
+		if ( fom == FOM_SNR ) {
+			ERROR("WARNING: %li reflections had infinite or "
+			      "invalid values of I/sigma(I).\n", n_rej);
+		} else {
+			ERROR("WARNING: %li reflections rejected by add_to_fom\n",
+			      n_rej);
+		}
+	}
+
+	if ( fom == FOM_COMPLETENESS ) {
+		calculate_possible(fctx, shells, cell, sym);
+	}
 
 	return fctx;
 }
 
 
-int fom_select_reflections(RefList *list1, RefList *list2,
-                           RefList *list1_acc, RefList *list2_acc,
-                           UnitCell *cell, SymOpList *sym,
-                           int anom, double rmin_fix, double rmax_fix,
-                           double sigma_cutoff, int ignore_negs,
-                           int zero_negs, int mul_cutoff)
+struct fom_rejections fom_select_reflection_pairs(RefList *list1, RefList *list2,
+                                                  RefList **plist1_acc,
+                                                  RefList **plist2_acc,
+                                                  UnitCell *cell, SymOpList *sym,
+                                                  int anom, double rmin_fix, double rmax_fix,
+                                                  double sigma_cutoff, int ignore_negs,
+                                                  int zero_negs, int mul_cutoff)
 {
 	Reflection *refl1;
 	RefListIterator *iter;
-	int ncom, nrej, nmul, nneg, nres, nbij, ncen;
-
-	/* Select reflections to be used */
-	ncom = 0;
-	nrej = 0;
-	nmul = 0;
-	nneg = 0;
-	nres = 0;
-	nbij = 0;
-	ncen = 0;
+	struct fom_rejections rej;
+	RefList *list1_acc;
+	RefList *list2_acc;
+
+	rej.common = 0;
+	rej.low_snr = 0;
+	rej.negative_deleted = 0;
+	rej.negative_zeroed = 0;
+	rej.few_measurements = 0;
+	rej.outside_resolution_range = 0;
+	rej.no_bijvoet = 0;
+	rej.centric = 0;
+	rej.nan_inf_value = 0;
+
+	list1_acc = reflist_new();
+	list2_acc = reflist_new();
+
 	for ( refl1 = first_refl(list1, &iter);
 	      refl1 != NULL;
 	      refl1 = next_refl(refl1, iter) )
@@ -766,20 +1056,27 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 		mul1 = get_redundancy(refl1);
 		mul2 = get_redundancy(refl2);
 
+		if ( !isfinite(val1) || !isfinite(val2)
+		  || !isfinite(esd1) || !isfinite(esd2) )
+		{
+			rej.nan_inf_value++;
+			continue;
+		}
+
 		if ( (val1 < sigma_cutoff * esd1)
 		  || (val2 < sigma_cutoff * esd2) )
 		{
-			nrej++;
+			rej.low_snr++;
 			continue;
 		}
 
 		if ( ignore_negs && ((val1 < 0.0) || (val2 < 0.0)) ) {
-			nneg++;
+			rej.negative_deleted++;
 			continue;
 		}
 
 		if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) {
-			nmul++;
+			rej.few_measurements++;
 			continue;
 		}
 
@@ -793,13 +1090,14 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 				val2 = 0.0;
 				d = 1;
 			}
-			if ( d ) nneg++;
+			if ( d ) rej.negative_zeroed++;
+			continue;
 		}
 
 		if ( rmin_fix > 0.0 ) {
 			double res = 2.0*resolution(cell, h, k, l);
 			if ( res < rmin_fix ) {
-				nres++;
+				rej.outside_resolution_range++;
 				continue;
 			}
 		}
@@ -807,7 +1105,7 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 		if ( rmax_fix > 0.0 ) {
 			double res = 2.0*resolution(cell, h, k, l);
 			if ( res > rmax_fix ) {
-				nres++;
+				rej.outside_resolution_range++;
 				continue;
 			}
 		}
@@ -820,7 +1118,7 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 		copy_data(refl2_acc, refl2);
 		set_intensity(refl2_acc, val2);
 
-		ncom++;
+		rej.common++;
 
 	}
 
@@ -833,7 +1131,7 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 		list1_acc = reflist_new();
 		list2_acc = reflist_new();
 
-		ncom = 0;
+		rej.common = 0;
 
 		for ( refl1 = first_refl(list1, &iter);
 		      refl1 != NULL;
@@ -857,7 +1155,7 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 			val2 = get_intensity(refl2);
 
 			if ( is_centric(h, k, l, sym) ) {
-				ncen++;
+				rej.centric++;
 				continue;
 			}
 
@@ -874,7 +1172,7 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 			}
 
 			if ( (refl1_bij == NULL) || (refl2_bij == NULL) ) {
-				nbij++;
+				rej.no_bijvoet++;
 				continue;
 			}
 
@@ -886,44 +1184,138 @@ int fom_select_reflections(RefList *list1, RefList *list2,
 			copy_data(refl2_acc, refl2);
 			set_intensity(refl2_acc, val2);
 
-			ncom++;
+			rej.common++;
 		}
 	}
 
-	if ( nrej > 0 ) {
-		STATUS("Discarded %i reflection pairs because either or both"
-		       " versions had I/sigma(I) < %f.\n", nrej, sigma_cutoff);
-	}
+	*plist1_acc = list1_acc;
+	*plist2_acc = list2_acc;
+	return rej;
+}
 
-	if ( ignore_negs && (nneg > 0) ) {
-		STATUS("Discarded %i reflection pairs because either or both"
-		       " versions had negative intensities.\n", nneg);
-	}
 
-	if ( zero_negs && (nneg > 0) ) {
-		STATUS("For %i reflection pairs, either or both versions had"
-		       " negative intensities which were set to zero.\n", nneg);
-	}
+struct fom_rejections fom_select_reflections(RefList *raw_list,
+                                             RefList **plist_acc,
+                                             UnitCell *cell, SymOpList *sym,
+                                             double rmin_fix, double rmax_fix,
+                                             double sigma_cutoff, int ignore_negs,
+                                             int zero_negs, int mul_cutoff)
+{
+	RefList *list;
+	Reflection *refl;
+	RefListIterator *iter;
+	struct fom_rejections rej;
 
-	if ( nmul > 0 ) {
-		STATUS("%i reflection pairs rejected because either or both"
-		       " versions had too few measurements.\n", nmul);
-	}
+	*plist_acc = NULL;
+
+	rej.common = 0;
+	rej.low_snr = 0;
+	rej.negative_deleted = 0;
+	rej.negative_zeroed = 0;
+	rej.few_measurements = 0;
+	rej.outside_resolution_range = 0;
+	rej.no_bijvoet = 0;
+	rej.centric = 0;
+	rej.nan_inf_value = 0;
+
+	list = reflist_new();
+	if ( list == NULL ) return rej;
+
+	for ( refl = first_refl(raw_list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		double val, sig;
+		int ig = 0;
+		Reflection *new;
+
+		get_indices(refl, &h, &k, &l);
+
+		val = get_intensity(refl);
+		sig = get_esd_intensity(refl);
+
+		if ( !isfinite(val) || !isfinite(sig) ) {
+			rej.nan_inf_value++;
+			continue;
+		}
+
+		if ( val < sigma_cutoff * sig ) {
+			rej.low_snr++;
+			ig = 1;
+		}
+
+		if ( ignore_negs && (val < 0.0) ) {
+			rej.negative_deleted++;
+			ig = 1;
+		}
+
+		if ( zero_negs && (val < 0.0) ) {
+			set_intensity(refl, 0.0);
+			rej.negative_zeroed++;
+		}
+
+		if ( rmin_fix > 0.0 ) {
+			double res = 2.0*resolution(cell, h, k, l);
+			if ( res < rmin_fix ) {
+				rej.outside_resolution_range++;
+				continue;
+			}
+		}
 
-	if ( nres > 0 ) {
-		STATUS("%i reflection pairs rejected because either or both"
-		       " versions were outside the resolution range.\n", nres);
+		if ( rmax_fix > 0.0 ) {
+			double res = 2.0*resolution(cell, h, k, l);
+			if ( res > rmax_fix ) {
+				rej.outside_resolution_range++;
+				continue;
+			}
+		}
+
+		if ( ig ) continue;
+
+		new = add_refl(list, h, k, l);
+		copy_data(new, refl);
 	}
 
-	if ( nbij > 0 ) {
-		STATUS("%i reflection pairs rejected because either or both"
-		       " versions did not have Bijvoet partners.\n", nres);
+	*plist_acc = list;
+	return rej;
+}
+
+
+int fom_overall_num_reflections(struct fom_context *fctx)
+{
+	int i;
+	long int n = 0;
+
+	for ( i=0; i<fctx->nshells; i++ ) {
+		n += fctx->cts[i];
 	}
+	return n;
+}
+
 
-	if ( ncen > 0 ) {
-		STATUS("%i reflection pairs rejected because they were"
-		       " centric.\n", ncen);
+int fom_shell_num_reflections(struct fom_context *fctx, int i)
+{
+	return fctx->cts[i];
+}
+
+
+int fom_overall_num_possible(struct fom_context *fctx)
+{
+	int i;
+	long int n = 0;
+
+	assert(fctx->fom == FOM_COMPLETENESS);
+
+	for ( i=0; i<fctx->nshells; i++ ) {
+		n += fctx->possible[i];
 	}
+	return n;
+}
+
 
-	return ncom;
+int fom_shell_num_possible(struct fom_context *fctx, int i)
+{
+	assert(fctx->fom == FOM_COMPLETENESS);
+	return fctx->possible[i];
 }
diff --git a/libcrystfel/src/fom.h b/libcrystfel/src/fom.h
index d3373044..c9c8dbfe 100644
--- a/libcrystfel/src/fom.h
+++ b/libcrystfel/src/fom.h
@@ -38,6 +38,19 @@
 #include <reflist.h>
 #include <symmetry.h>
 
+struct fom_rejections
+{
+	int common;
+	int low_snr;
+	int negative_deleted;
+	int negative_zeroed;
+	int few_measurements;
+	int outside_resolution_range;
+	int no_bijvoet;
+	int centric;
+	int nan_inf_value;
+};
+
 enum fom_type
 {
 	FOM_R1I,
@@ -51,7 +64,13 @@ enum fom_type
 	FOM_RANO,
 	FOM_RANORSPLIT,
 	FOM_D1SIG,
-	FOM_D2SIG
+	FOM_D2SIG,
+	FOM_NUM_MEASUREMENTS,
+	FOM_REDUNDANCY,
+	FOM_SNR,
+	FOM_MEAN_INTENSITY,
+	FOM_STDDEV_INTENSITY,
+	FOM_COMPLETENESS,
 };
 
 struct fom_shells
@@ -61,51 +80,53 @@ struct fom_shells
 	double *rmaxs;
 };
 
-struct fom_context
-{
-	enum fom_type fom;
-	int nshells;
-	int *cts;
-
-	/* For R-factors */
-	double *num;
-	double *den;
+struct fom_context;
+
+extern struct fom_rejections fom_select_reflection_pairs(RefList *list1,
+                                                         RefList *list2,
+                                                         RefList **plist1_acc,
+                                                         RefList **plist2_acc,
+                                                         UnitCell *cell,
+                                                         SymOpList *sym,
+                                                         int anom,
+                                                         double rmin_fix,
+                                                         double rmax_fix,
+                                                         double sigma_cutoff,
+                                                         int ignore_negs,
+                                                         int zero_negs,
+                                                         int mul_cutoff);
+
+extern struct fom_rejections fom_select_reflections(RefList *list,
+                                                    RefList **plist_acc,
+                                                    UnitCell *cell,
+                                                    SymOpList *sym,
+                                                    double rmin_fix,
+                                                    double rmax_fix,
+                                                    double sigma_cutoff,
+                                                    int ignore_negs,
+                                                    int zero_negs,
+                                                    int mul_cutoff);
 
-	/* For "double" R-factors */
-	double *num2;
-	double *den2;
-
-	/* For CCs */
-	double **vec1;
-	double **vec2;
-	int *n;
-	int nmax;
-
-	/* For "counting" things e.g. d1sig or d2sig */
-	int *n_within;
-};
-
-extern int fom_select_reflections(RefList *list1, RefList *list2,
-                                  RefList *list1_acc, RefList *list2_acc,
-                                  UnitCell *cell, SymOpList *sym,
-                                  int anom, double rmin_fix, double rmax_fix,
-                                  double sigma_cutoff, int ignore_negs,
-                                  int zero_negs, int mul_cutoff);
 
 extern struct fom_context *fom_calculate(RefList *list1, RefList *list2,
                                          UnitCell *cell,
                                          struct fom_shells *shells,
                                          enum fom_type fom, int noscale,
-                                         SymOpList *sym);
+                                         const SymOpList *sym);
 
 extern struct fom_shells *fom_make_resolution_shells(double rmin, double rmax,
                                                      int nshells);
 
-extern double fom_shell_label(struct fom_shells *s, int i);
+extern double fom_shell_centre(struct fom_shells *s, int i);
+
+extern double fom_overall_value(struct fom_context *fctx);
+extern double fom_shell_value(struct fom_context *fctx, int i);
 
-extern double fom_shell(struct fom_context *fctx, int i);
+extern int fom_overall_num_reflections(struct fom_context *fctx);
+extern int fom_shell_num_reflections(struct fom_context *fctx, int i);
 
-extern double fom_overall(struct fom_context *fctx);
+extern int fom_overall_num_possible(struct fom_context *fctx);
+extern int fom_shell_num_possible(struct fom_context *fctx, int i);
 
 extern enum fom_type fom_type_from_string(const char *s);
 
diff --git a/src/check_hkl.c b/src/check_hkl.c
index f743d575..d95fc51f 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -45,6 +45,7 @@
 #include <reflist.h>
 #include <reflist-utils.h>
 #include <cell-utils.h>
+#include <fom.h>
 
 #include "version.h"
 
@@ -390,67 +391,16 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
                         double rmin_fix, double rmax_fix, int nshells,
 			const char *shell_file)
 {
-	unsigned long *possible;
-	unsigned long *measurements;
-	unsigned long *measured;
-	unsigned long *snr_measured;
-	double total_vol, vol_per_shell;
-	double *rmins;
-	double *rmaxs;
-	double *snr;
-	double *mean;
-	double *var;
 	double rmin, rmax;
-	signed int h, k, l;
 	int i;
 	FILE *fh;
-	double snr_total = 0;
-	unsigned long nrefl = 0;
-	unsigned long nmeastot = 0;
-	unsigned long nout = 0;
-	unsigned long possible_tot = 0;
-	unsigned long nsilly = 0;
-	Reflection *refl;
-	RefListIterator *iter;
-	RefList *counted;
-	int hmax, kmax, lmax;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-
-	possible = malloc(nshells*sizeof(unsigned long));
-	measurements = malloc(nshells*sizeof(unsigned long));
-	measured = malloc(nshells*sizeof(unsigned long));
-	snr_measured = malloc(nshells*sizeof(unsigned long));
-	if ( (possible == NULL) || (measurements == NULL)
-	  || (measured == NULL) || (snr_measured == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		free(possible);
-		free(measurements);
-		free(measured);
-		free(snr_measured);
-		return;
-	}
-
-	rmins = malloc(nshells*sizeof(double));
-	rmaxs = malloc(nshells*sizeof(double));
-	snr = malloc(nshells*sizeof(double));
-	mean = malloc(nshells*sizeof(double));
-	var = malloc(nshells*sizeof(double));
-	if ( (rmins == NULL) || (rmaxs == NULL) || (snr == NULL)
-	  || (mean == NULL) || (var == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		free(possible);
-		free(measurements);
-		free(measured);
-		free(snr_measured);
-		free(rmins);
-		free(rmaxs);
-		free(snr);
-		free(mean);
-		free(var);
-		return;
-	}
+	struct fom_shells *shells;
+	struct fom_context *nmeas_ctx;
+	struct fom_context *red_ctx;
+	struct fom_context *snr_ctx;
+	struct fom_context *mean_ctx;
+	struct fom_context *stddev_ctx;
+	struct fom_context *compl_ctx;
 
 	fh = fopen(shell_file, "w");
 	if ( fh == NULL ) {
@@ -458,224 +408,72 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
 		return;
 	}
 
-	for ( i=0; i<nshells; i++ ) {
-		possible[i] = 0;
-		measured[i] = 0;
-		snr_measured[i] = 0;
-		measurements[i] = 0;
-		snr[i] = 0;
-		var[i] = 0;
-		mean[i] = 0;
-	}
-
 	resolution_limits(list, cell, &rmin, &rmax);
 	STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
 
-	/* Widen the range just a little bit */
-	rmin -= 0.001e9;
-	rmax += 0.001e9;
-
 	/* Fixed resolution shells if needed */
 	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
 	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
 
-	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++ ) {
-
-		double r;
-
-		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;
-
-	/* Count the number of reflections possible in each shell */
-	counted = reflist_new();
-	cell_get_cartesian(cell, &ax, &ay, &az,
-	                         &bx, &by, &bz,
-	                         &cx, &cy, &cz);
-	hmax = rmax * modulus(ax, ay, az);
-	kmax = rmax * modulus(bx, by, bz);
-	lmax = rmax * modulus(cx, cy, cz);
-	for ( h=-hmax; h<=hmax; h++ ) {
-	for ( k=-kmax; k<=kmax; k++ ) {
-	for ( l=-lmax; l<=lmax; l++ ) {
-
-		double d;
-		signed int hs, ks, ls;
-		int bin;
-
-		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;
-
-		bin = -1;
-		for ( i=0; i<nshells; i++ ) {
-			if ( (d>rmins[i]) && (d<=rmaxs[i]) ) {
-				bin = i;
-				break;
-			}
-		}
-		if ( bin == -1 ) continue;
-
-		if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
-		add_refl(counted, hs, ks, ls);
-
-		possible[bin]++;
-		possible_tot++;
-
-	}
-	}
-	}
-	reflist_free(counted);
-
-	/* Calculate means */
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l, hs, ks, ls;
-		double d, val, esd;
-		int bin;
-		int j;
-
-		get_indices(refl, &h, &k, &l);
-		if ( forbidden_reflection(cell, h, k, l) ) continue;
-
-		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);
-
-		bin = -1;
-		for ( j=0; j<nshells; j++ ) {
-			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
-				bin = j;
-				break;
-			}
-		}
-		if ( bin == -1 ) continue;
-
-		measured[bin]++;
-		mean[bin] += get_intensity(refl);
-
-		if ( !isfinite(val/esd) ) nsilly++;
-
-	}
-
-	for ( i=0; i<nshells; i++ ) {
-		mean[i] /= (double)measured[i];
-	}
-
-	/* Characterise the data set */
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l;
-		double d;
-		int bin;
-		int j;
-		double val, esd;
-
-		get_indices(refl, &h, &k, &l);
-		if ( forbidden_reflection(cell, h, k, l) ) continue;
-
-		d = resolution(cell, h, k, l) * 2.0;
-		val = get_intensity(refl);
-		esd = get_esd_intensity(refl);
-
-		bin = -1;
-		for ( j=0; j<nshells; j++ ) {
-			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
-				bin = j;
-				break;
-			}
-		}
-		if ( bin == -1 ) {
-			nout++;
-			continue;
-		}
-
-		/* measured[bin] was done earlier */
-		measurements[bin] += get_redundancy(refl);
-
-		if ( isfinite(val/esd) ) {
-			snr[bin] += val / esd;
-			snr_total += val / esd;
-			snr_measured[bin]++;
-		} else {
-			nsilly++;
-		}
-
-		nrefl++;
-		nmeastot += get_redundancy(refl);
-
-		var[bin] += pow(val-mean[bin], 2.0);
-
-	}
+	shells = fom_make_resolution_shells(rmin, rmax,  nshells);
 
 	STATUS("Overall values within specified resolution range:\n");
-	STATUS("%li measurements in total.\n", nmeastot);
-	STATUS("%li reflections in total.\n", nrefl);
-	STATUS("%li reflections possible.\n", possible_tot);
-	STATUS("Overall <snr> = %f\n", snr_total/(double)nrefl);
-	STATUS("Overall redundancy = %f measurements/unique reflection\n",
-	       nmeastot/(double)nrefl);
-	STATUS("Overall completeness = %f %%\n", 100.0*nrefl/(double)possible_tot);
 
-	if ( nout ) {
-		STATUS("WARNING; %li reflections outside resolution range.\n",
-		       nout);
-	}
-
-	if ( nsilly ) {
-		STATUS("WARNING; %li reflections had infinite or invalid values"
-		       " of I/sigma(I).\n", nsilly);
-	}
+	nmeas_ctx = fom_calculate(list, NULL, cell, shells,
+	                          FOM_NUM_MEASUREMENTS, 0, sym);
+	red_ctx = fom_calculate(list, NULL, cell, shells,
+	                        FOM_REDUNDANCY, 0, sym);
+	snr_ctx = fom_calculate(list, NULL, cell, shells,
+	                        FOM_SNR, 0, sym);
+	mean_ctx = fom_calculate(list, NULL, cell, shells,
+	                         FOM_MEAN_INTENSITY, 0, sym);
+	stddev_ctx = fom_calculate(list, NULL, cell, shells,
+	                           FOM_STDDEV_INTENSITY, 0, sym);
+	compl_ctx = fom_calculate(list, NULL, cell, shells,
+	                          FOM_COMPLETENESS, 0, sym);
+
+	STATUS("%.0f measurements in total.\n",
+	       fom_overall_value(nmeas_ctx));
+	STATUS("%li reflections in total.\n",
+	       fom_overall_num_reflections(compl_ctx));
+	STATUS("%li reflections possible.\n",
+	       fom_overall_num_possible(compl_ctx));
+	STATUS("Overall <snr> = %f\n", fom_overall_value(snr_ctx));
+	STATUS("Overall redundancy = %f measurements/unique reflection\n",
+	       fom_overall_value(red_ctx));
+	STATUS("Overall completeness = %f %%\n",
+	       100.0*fom_overall_value(compl_ctx));
 
 	fprintf(fh, "Center 1/nm  # refs Possible  Compl       "
 		    "Meas   Red   SNR    Std dev       Mean     d(A)    "
 		    "Min 1/nm   Max 1/nm\n");
 	for ( i=0; i<nshells; i++ ) {
 
-		double cen;
-		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
-		fprintf(fh, "%10.3f %8li %8li %6.2f %10li %5.1f"
+		long int measured, possible;
+
+		measured = fom_shell_num_reflections(compl_ctx, i);
+		possible = fom_shell_num_possible(compl_ctx, i);
+
+		fprintf(fh, "%10.3f %8li %8li %6.2f %10.0f %5.1f"
 		            " %5.2f %10.2f %10.2f %8.2f  %10.3f %10.3f\n",
-		        cen*1.0e-9,
-		        measured[i],
-		        possible[i],
-		        100.0*(double)measured[i]/possible[i],
-		        measurements[i],
-		        (double)measurements[i]/measured[i],
-		        snr[i]/(double)snr_measured[i],
-		        sqrt(var[i]/measured[i]),
-		        mean[i], (1.0/cen)*1e10,
-			rmins[i]*1.0e-9, rmaxs[i]*1.0e-9);
+		        fom_shell_centre(shells, i)*1.0e-9,
+		        measured,
+		        possible,
+		        100.0*fom_shell_value(compl_ctx, i),
+		        fom_shell_value(nmeas_ctx, i),
+		        fom_shell_value(red_ctx, i),
+		        fom_shell_value(snr_ctx, i),
+		        fom_shell_value(stddev_ctx, i),
+		        fom_shell_value(mean_ctx, i),
+		        (1.0/fom_shell_centre(shells, i))*1e10,
+			shells->rmins[i]*1.0e-9,
+			shells->rmaxs[i]*1.0e-9);
 
 	}
 
 	fclose(fh);
 
 	STATUS("Resolution shell information written to %s.\n", shell_file);
-
-	free(possible);
-	free(measurements);
-	free(measured);
-	free(snr_measured);
-	free(rmins);
-	free(rmaxs);
-	free(snr);
-	free(mean);
-	free(var);
 }
 
 
@@ -711,10 +509,7 @@ int main(int argc, char *argv[])
 	SymOpList *sym;
 	RefList *raw_list;
 	RefList *list;
-	Reflection *refl;
-	RefListIterator *iter;
 	char *cellfile = NULL;
-	int rej = 0;
 	float rmin_fix = -1.0;
 	float rmax_fix = -1.0;
 	float sigma_cutoff = -INFINITY;
@@ -725,9 +520,8 @@ int main(int argc, char *argv[])
 	int ltest = 0;
 	int ignorenegs = 0;
 	int zeronegs = 0;
-	int nneg = 0;
-	int nres = 0;
 	float highres, lowres;
+	struct fom_rejections rej;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -905,74 +699,32 @@ int main(int argc, char *argv[])
 	}
 
 	/* Reject some reflections */
-	list = reflist_new();
-	for ( refl = first_refl(raw_list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		double val, sig;
-		int ig = 0;
-		Reflection *new;
-
-		get_indices(refl, &h, &k, &l);
+	rej = fom_select_reflections(raw_list, &list, cell, sym,
+	                             rmin_fix, rmax_fix, sigma_cutoff,
+	                             ignorenegs, zeronegs, 0);
 
-		val = get_intensity(refl);
-		sig = get_esd_intensity(refl);
-
-		if ( val < sigma_cutoff * sig ) {
-			rej++;
-			ig = 1;
-		}
-
-		if ( ignorenegs && (val < 0.0) ) {
-			nneg++;
-			ig = 1;
-		}
-
-		if ( zeronegs && (val < 0.0) ) {
-			set_intensity(refl, 0.0);
-			nneg++;
-		}
-
-		if ( rmin_fix > 0.0 ) {
-			double res = 2.0*resolution(cell, h, k, l);
-			if ( res < rmin_fix ) {
-				nres++;
-				continue;
-			}
-		}
-
-		if ( rmax_fix > 0.0 ) {
-			double res = 2.0*resolution(cell, h, k, l);
-			if ( res > rmax_fix ) {
-				nres++;
-				continue;
-			}
-		}
-
-		if ( ig ) continue;
-
-		new = add_refl(list, h, k, l);
-		copy_data(new, refl);
-
-	}
 	STATUS("Discarded %i reflections (out of %i) with I/sigma(I) < %f\n",
-	       rej, num_reflections(raw_list), sigma_cutoff);
+	       rej.low_snr, num_reflections(raw_list), sigma_cutoff);
 	reflist_free(raw_list);
 
-	if ( ignorenegs && (nneg > 0) ) {
+	if ( rej.negative_deleted > 0 ) {
 		STATUS("Discarded %i reflections because they had negative "
-		       "intensities.\n", nneg);
+		       "intensities.\n", rej.negative_deleted);
 	}
 
-	if ( zeronegs && (nneg > 0) ) {
-		STATUS("Set %i negative intensities to zero\n", nneg);
+	if ( rej.negative_zeroed > 0 ) {
+		STATUS("Set %i negative intensities to zero\n",
+		       rej.negative_zeroed);
 	}
 
-	if ( nres > 0 ) {
+	if ( rej.outside_resolution_range > 0 ) {
 		STATUS("%i reflections rejected because they were outside the "
-		       "resolution range.\n", nres);
+		       "resolution range.\n", rej.outside_resolution_range);
+	}
+
+	if ( rej.nan_inf_value ) {
+		STATUS("WARNING: %i reflections had infinite or invalid values"
+		       " of I or sigma(I).\n", rej.nan_inf_value);
 	}
 
 	if ( wilson ) {
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 2cd7207b..cb5bee38 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -106,55 +106,58 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 	switch ( fom ) {
 
 		case FOM_R1I :
-		STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall(fctx));
+		STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_R1F :
-		STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall(fctx));
+		STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_R2 :
-		STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall(fctx));
+		STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_RSPLIT :
-		STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall(fctx));
+		STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_CC :
-		STATUS("Overall CC = %.7f\n", fom_overall(fctx));
+		STATUS("Overall CC = %.7f\n", fom_overall_value(fctx));
 		break;
 
 		case FOM_CCSTAR :
-		STATUS("Overall CC* = %.7f\n", fom_overall(fctx));
+		STATUS("Overall CC* = %.7f\n", fom_overall_value(fctx));
 		break;
 
 		case FOM_CCANO :
-		STATUS("Overall CCano = %.7f\n", fom_overall(fctx));
+		STATUS("Overall CCano = %.7f\n", fom_overall_value(fctx));
 		break;
 
 		case FOM_CRDANO :
-		STATUS("Overall CRDano = %.7f\n", fom_overall(fctx));
+		STATUS("Overall CRDano = %.7f\n", fom_overall_value(fctx));
 		break;
 
 		case FOM_RANO :
-		STATUS("Overall Rano =  %.2f %%\n", 100.0*fom_overall(fctx));
+		STATUS("Overall Rano =  %.2f %%\n", 100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_RANORSPLIT :
-		STATUS("Overall Rano/Rsplit =  %.7f\n", fom_overall(fctx));
+		STATUS("Overall Rano/Rsplit =  %.7f\n", fom_overall_value(fctx));
 		break;
 
 		case FOM_D1SIG :
 		STATUS("Fraction of differences less than 1 sigma = %.7f %%\n",
-		       100.0*fom_overall(fctx));
+		       100.0*fom_overall_value(fctx));
 		break;
 
 		case FOM_D2SIG :
 		STATUS("Fraction of differences less than 2 sigma = %.7f %%\n",
-		       100.0*fom_overall(fctx));
+		       100.0*fom_overall_value(fctx));
 		break;
 
+		default :
+		ERROR("Unhandled figure of merit type\n");
+		break;
 
 	}
 
@@ -217,14 +220,17 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		fprintf(fh, "%s D<2sigma/%%     nref%s\n", t1, t2);
 		break;
 
+		default :
+		break;
+
 	}
 
 	for ( i=0; i<nshells; i++ ) {
 
 		double r, cen;
 
-		cen = fom_shell_label(shells, i);
-		r = fom_shell(fctx, i);
+		cen = fom_shell_centre(shells, i);
+		r = fom_shell_value(fctx, i);
 
 		switch ( fom ) {
 
@@ -235,7 +241,8 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 			case FOM_RANO :
 			fprintf(fh, "%10.3f %10.2f %10i %10.2f "
 			        "%10.3f  %10.3f\n",
-			        cen*1.0e-9, r*100.0, fctx->cts[i],
+			        cen*1.0e-9, r*100.0,
+			        fom_shell_num_reflections(fctx, i),
 			        (1.0/cen)*1e10,
 			        shells->rmins[i]*1.0e-9,
 			        shells->rmaxs[i]*1.0e-9);
@@ -247,7 +254,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 			case FOM_CRDANO :
 			fprintf(fh, "%10.3f %10.7f %10i %10.2f "
 			        "%10.3f  %10.3f\n",
-			        cen*1.0e-9, r, fctx->cts[i], (1.0/cen)*1e10,
+			        cen*1.0e-9, r,
+			        fom_shell_num_reflections(fctx, i),
+			        (1.0/cen)*1e10,
 			        shells->rmins[i]*1.0e-9,
 			        shells->rmaxs[i]*1.0e-9);
 			break;
@@ -255,7 +264,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 			case FOM_RANORSPLIT :
 			fprintf(fh, "%10.3f    %10.7f %10i %10.2f "
 			        "%10.3f  %10.3f\n",
-			        cen*1.0e-9, r, fctx->cts[i], (1.0/cen)*1e10,
+			        cen*1.0e-9, r,
+			        fom_shell_num_reflections(fctx, i),
+			        (1.0/cen)*1e10,
 			        shells->rmins[i]*1.0e-9,
 			        shells->rmaxs[i]*1.0e-9);
 			break;
@@ -264,12 +275,16 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 			case FOM_D2SIG :
 			fprintf(fh, "%10.3f %10.2f %10i %10.2f "
 			        "%10.3f  %10.3f\n",
-			        cen*1.0e-9, r*100.0, fctx->cts[i],
+			        cen*1.0e-9, r*100.0,
+			        fom_shell_num_reflections(fctx, i),
 			        (1.0/cen)*1e10,
 			        shells->rmins[i]*1.0e-9,
 			        shells->rmaxs[i]*1.0e-9);
 			break;
 
+			default :
+			break;
+
 		}
 
 	}
@@ -311,7 +326,6 @@ int main(int argc, char *argv[])
 	char *sym_str_fromfile1 = NULL;
 	char *sym_str_fromfile2 = NULL;
 	SymOpList *sym;
-	int ncom;
 	RefList *list1_acc;
 	RefList *list2_acc;
 	RefList *list1;
@@ -332,6 +346,7 @@ int main(int argc, char *argv[])
 	float highres, lowres;
 	int mul_cutoff = 0;
 	int anom;
+	struct fom_rejections rej;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -478,7 +493,7 @@ int main(int argc, char *argv[])
 			      " intensities.\n");
 			ERROR("Please try again with --ignore-negs or"
 			      " --zero-negs.\n");
-			exit(1);
+			return 1;
 
 			case FOM_R2 :
 			case FOM_R1I :
@@ -492,6 +507,10 @@ int main(int argc, char *argv[])
 			case FOM_D1SIG :
 			case FOM_D2SIG :
 			break;
+
+			default :
+			ERROR("Unhandled figure of merit!\n");
+			return 1;
 		}
 	}
 
@@ -579,6 +598,10 @@ int main(int argc, char *argv[])
 			      " try again using a non-centrosymmetric point"
 			      " group for '-y'.\n");
 			return 1;
+
+			default :
+			ERROR("Unhandled figure of merit type!\n");
+			return 1;
 		}
 	}
 
@@ -622,21 +645,60 @@ int main(int argc, char *argv[])
 	reflist_free(list1_raw);
 	reflist_free(list2_raw);
 
-	list1_acc = reflist_new();
-	list2_acc = reflist_new();
 	anom = ( (fom == FOM_CCANO) || (fom == FOM_CRDANO)
 	        || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) );
-	ncom = fom_select_reflections(list1, list2, list1_acc, list2_acc,
-	                              cell, sym,
-	                              anom, rmin_fix, rmax_fix, sigma_cutoff,
-	                              config_ignorenegs, config_zeronegs,
-	                              mul_cutoff);
+	rej = fom_select_reflection_pairs(list1, list2, &list1_acc, &list2_acc,
+	                                  cell, sym,
+	                                  anom, rmin_fix, rmax_fix, sigma_cutoff,
+	                                  config_ignorenegs, config_zeronegs,
+	                                  mul_cutoff);
 	reflist_free(list1);
 	reflist_free(list2);
 
 	gsl_set_error_handler_off();
 
-	STATUS("%i reflection pairs accepted.\n", ncom);
+	if ( rej.low_snr > 0 ) {
+		STATUS("Discarded %i reflection pairs because either or both"
+		       " versions had I/sigma(I) < %f.\n",
+		       rej.low_snr, sigma_cutoff);
+	}
+
+	if ( rej.negative_deleted > 0 ) {
+		STATUS("Discarded %i reflection pairs because either or both"
+		       " versions had negative intensities.\n",
+		       rej.negative_deleted);
+	}
+
+	if ( rej.negative_zeroed > 0 ) {
+		STATUS("For %i reflection pairs, either or both versions had"
+		       " negative intensities which were set to zero.\n",
+		       rej.negative_zeroed);
+	}
+
+	if ( rej.few_measurements > 0 ) {
+		STATUS("%i reflection pairs rejected because either or both"
+		       " versions had too few measurements.\n",
+		       rej.few_measurements);
+	}
+
+	if ( rej.outside_resolution_range > 0 ) {
+		STATUS("%i reflection pairs rejected because either or both"
+		       " versions were outside the resolution range.\n",
+		       rej.outside_resolution_range);
+	}
+
+	if ( rej.no_bijvoet > 0 ) {
+		STATUS("%i reflection pairs rejected because either or both"
+		       " versions did not have Bijvoet partners.\n",
+		       rej.no_bijvoet);
+	}
+
+	if ( rej.centric > 0 ) {
+		STATUS("%i reflection pairs rejected because they were"
+		       " centric.\n", rej.centric);
+	}
+
+	STATUS("%i reflection pairs accepted.\n", rej.common);
 
 	resolution_limits(list1_acc, cell, &rmin, &rmax);
 	resolution_limits(list2_acc, cell, &rmin, &rmax);