Reduce the scope of "count"

Lists of counts had pervaded every corner of CrystFEL, being used as markers
for the presence of reflections.  Now we have a better way of doing this,
the ReflItemList, and few parts of the suite apart from process_hkl have any
business knowing how many observations were made of a particular reflection.

20 files changed, 345 insertions, 433 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index e6caee6f..e69f3753 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -49,13 +49,9 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	char *afile = NULL;
 	char *bfile = NULL;
-	signed int h, k, l;
 	double scale, R2, Rmerge, pearson;
-	unsigned int *c1;
-	unsigned int *c2;
-	int i;
-	int nc1, nc2, ncom;
-	unsigned int *outcounts;
+	int i, ncom;
+	ReflItemList *i1, *i2, *icommon;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -94,72 +90,53 @@ int main(int argc, char *argv[])
 	bfile = strdup(argv[optind]);
 
 	cell = load_cell_from_pdb("molecule.pdb");
-	c1 = new_list_count();
-	ref1 = read_reflections(afile, c1, NULL);
+	ref1 = new_list_intensity();
+	i1 = read_reflections(afile, ref1, NULL, NULL);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
 		return 1;
 	}
-	c2 = new_list_count();
-	ref2 = read_reflections(bfile, c2, NULL);
+	ref2 = new_list_intensity();
+	i2 = read_reflections(bfile, ref2, NULL, NULL);
 	if ( ref2 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", bfile);
 		return 1;
 	}
-	out = new_list_intensity();
-	outcounts = new_list_count();
 
-	/* Knock out the zero beam, in case it's present */
-	set_count(c1, 0, 0, 0, 0);
-	set_count(c2, 0, 0, 0, 0);
+	/* Find common reflections */
+	icommon = intersection_items(i1, i2);
+	ncom = num_items(icommon);
 
-	/* Divide by number of counts, since we're not interested in them */
-	divide_down(ref1, c1);
-	divide_down(ref2, c2);
+	/* List for output scale factor map */
+	out = new_list_intensity();
 
-	for ( h=-INDMAX; h<INDMAX; h++ ) {
-	for ( k=-INDMAX; k<INDMAX; k++ ) {
-	for ( l=-INDMAX; l<INDMAX; l++ ) {
+	for ( i=0; i<ncom; i++ ) {
 
 		double i1, i2;
-		unsigned int c1s, c2s;
+		struct refl_item *it;
+		signed int h, k, l;
 
-		c1s = lookup_count(c1, h, k, l);
-		c2s = lookup_count(c2, h, k, l);
+		it = get_item(icommon, i);
+		h = it->h;  k = it->k;  l = it->l;
 
 		i1 = lookup_intensity(ref1, h, k, l);
 		i2 = lookup_intensity(ref2, h, k, l);
 
-		if ( c1s && c2s ) {
-			if ( (i1 != 0.0) && (i2 != 0.0) ) {
-				set_intensity(out, h, k, l,
-				              (i1/(double)c1s)/i2/(double)c2s);
-				set_count(outcounts, h, k, l, 1);
-			}
-		}
+		set_intensity(out, h, k, l, i1/i2);
 
 	}
-	}
-	}
 
-	nc1 = 0;
-	nc2 = 0;
-	ncom = 0;
-	for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-		nc1 += c1[i];
-		nc2 += c2[i];
-		ncom += c1[i] && c2[i];
-	}
-	STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
-	R2 = stat_r2(ref1, c1, ref2, c2, &scale);
+	STATUS("%i,%i reflections: %i in common\n",
+	       num_items(i1), num_items(i2), ncom);
+	R2 = stat_r2(ref1, ref2, icommon, &scale);
 	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
-	Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
+	Rmerge = stat_rmerge(ref1, ref2, icommon, &scale);
 	STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
-	pearson = stat_pearson(ref1, c1, ref2, c2);
+	pearson = stat_pearson(ref1, ref2, icommon);
 	STATUS("Pearson r = %5.4f\n", pearson);
 
 	if ( outfile != NULL ) {
-		write_reflections(outfile, outcounts, out, NULL, 0, cell, 1);
+		write_reflections(outfile, icommon, out, NULL, NULL, cell);
 	}
 
 	return 0;
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 86594a47..7a5bef85 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -342,8 +342,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 
 /* Setup the OpenCL stuff, create buffers, load the structure factor table */
 struct gpu_context *setup_gpu(int no_sfac, struct image *image,
-                              const double *intensities,
-                              const unsigned int *counts)
+                              const double *intensities)
 {
 	struct gpu_context *gctx;
 	cl_uint nplat;
@@ -404,11 +403,7 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
 	intensities_ptr = malloc(intensities_size);
 	if ( intensities != NULL ) {
 		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			if ( counts[i] == 1 ) {
-				intensities_ptr[i] = intensities[i];
-			} else {
-				intensities_ptr[i] = 1.0e20;
-			}
+			intensities_ptr[i] = intensities[i];
 		}
 	} else {
 		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
index 003c4c3a..2438e8f1 100644
--- a/src/diffraction-gpu.h
+++ b/src/diffraction-gpu.h
@@ -26,8 +26,7 @@ struct gpu_context;
 extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
                                 int na, int nb, int nc, UnitCell *ucell);
 extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
-                                     const double *intensities,
-                                     const unsigned int *counts);
+                                     const double *intensities);
 extern void cleanup_gpu(struct gpu_context *gctx);
 
 #else
@@ -40,8 +39,7 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 }
 
 static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
-                                     const double *intensities,
-                                     const unsigned int *counts)
+                                     const double *intensities)
 {
 	return NULL;
 }
diff --git a/src/diffraction.c b/src/diffraction.c
index fb512e62..efc56231 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -71,13 +71,11 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
 
 
 static double interpolate_linear(const double *ref,
-                                 const unsigned int *counts,
                                  float hd, signed int k, signed int l)
 {
 	signed int h;
 	double val1, val2;
 	float f;
-	unsigned int c1, c2;
 
 	h = (signed int)hd;
 	if ( hd < 0.0 ) h -= 1;
@@ -86,30 +84,15 @@ static double interpolate_linear(const double *ref,
 
 	val1 = lookup_intensity(ref, h, k, l);
 	val2 = lookup_intensity(ref, h+1, k, l);
-	c1 = lookup_count(counts, h, k, l);
-	c2 = lookup_count(counts, h+1, k, l);
 
-	if ( c1 == 0 ) {
-		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-		      h, k, l);
-		return 1.0e20;
-	}
-
-	if ( c2 == 0 ) {
-		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-		      h+1, k, l);
-		return 1.0e20;
-	}
-
-	val1 = val1 / (double)c1;
-	val2 = val2 / (double)c2;
+	val1 = val1;
+	val2 = val2;
 
 	return (1.0-f)*val1 + f*val2;
 }
 
 
 static double interpolate_bilinear(const double *ref,
-                                   const unsigned int *counts,
                                    float hd, float kd, signed int l)
 {
 	signed int k;
@@ -121,15 +104,14 @@ static double interpolate_bilinear(const double *ref,
 	f = kd - (float)k;
 	assert(f >= 0.0);
 
-	val1 = interpolate_linear(ref, counts, hd, k, l);
-	val2 = interpolate_linear(ref, counts, hd, k+1, l);
+	val1 = interpolate_linear(ref, hd, k, l);
+	val2 = interpolate_linear(ref, hd, k+1, l);
 
 	return (1.0-f)*val1 + f*val2;
 }
 
 
 static double interpolate_intensity(const double *ref,
-                                    const unsigned int *counts,
                                     float hd, float kd, float ld)
 {
 	signed int l;
@@ -141,15 +123,14 @@ static double interpolate_intensity(const double *ref,
 	f = ld - (float)l;
 	assert(f >= 0.0);
 
-	val1 = interpolate_bilinear(ref, counts, hd, kd, l);
-	val2 = interpolate_bilinear(ref, counts, hd, kd, l+1);
+	val1 = interpolate_bilinear(ref, hd, kd, l);
+	val2 = interpolate_bilinear(ref, hd, kd, l+1);
 
 	return (1.0-f)*val1 + f*val2;
 }
 
 
 static double complex interpolate_phased_linear(const double *ref,
-                                                const unsigned int *counts,
                                                 const double *phases,
                                                 float hd,
                                                 signed int k, signed int l)
@@ -157,7 +138,6 @@ static double complex interpolate_phased_linear(const double *ref,
 	signed int h;
 	double val1, val2;
 	float f;
-	unsigned int c1, c2;
 	double ph1, ph2;
 	double re1, re2, im1, im2;
 	double re, im;
@@ -169,25 +149,11 @@ static double complex interpolate_phased_linear(const double *ref,
 
 	val1 = lookup_intensity(ref, h, k, l);
 	val2 = lookup_intensity(ref, h+1, k, l);
-	c1 = lookup_count(counts, h, k, l);
-	c2 = lookup_count(counts, h+1, k, l);
 	ph1 = lookup_phase(phases, h, k, l);
 	ph2 = lookup_phase(phases, h+1, k, l);
 
-	if ( c1 == 0 ) {
-		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-		      h, k, l);
-		return 1.0e20;
-	}
-
-	if ( c2 == 0 ) {
-		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-		      h+1, k, l);
-		return 1.0e20;
-	}
-
-	val1 = val1 / (double)c1;
-	val2 = val2 / (double)c2;
+	val1 = val1;
+	val2 = val2;
 
 	/* Calculate real and imaginary parts */
 	re1 = val1 * cos(ph1);
@@ -203,7 +169,6 @@ static double complex interpolate_phased_linear(const double *ref,
 
 
 static double complex interpolate_phased_bilinear(const double *ref,
-                                                  const unsigned int *counts,
                                                   const double *phases,
                                                   float hd, float kd,
                                                   signed int l)
@@ -217,15 +182,14 @@ static double complex interpolate_phased_bilinear(const double *ref,
 	f = kd - (float)k;
 	assert(f >= 0.0);
 
-	val1 = interpolate_phased_linear(ref, counts, phases, hd, k, l);
-	val2 = interpolate_phased_linear(ref, counts, phases, hd, k+1, l);
+	val1 = interpolate_phased_linear(ref, phases, hd, k, l);
+	val2 = interpolate_phased_linear(ref, phases, hd, k+1, l);
 
 	return (1.0-f)*val1 + f*val2;
 }
 
 
 static double interpolate_phased_intensity(const double *ref,
-                                           const unsigned int *counts,
                                            const double *phases,
                                            float hd, float kd, float ld)
 {
@@ -238,16 +202,15 @@ static double interpolate_phased_intensity(const double *ref,
 	f = ld - (float)l;
 	assert(f >= 0.0);
 
-	val1 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l);
-	val2 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l+1);
+	val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l);
+	val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1);
 
 	return cabs((1.0-f)*val1 + f*val2);
 }
 
 
 /* Look up the structure factor for the nearest Bragg condition */
-static double molecule_factor(const double *intensities,
-                              const unsigned int *counts, const double *phases,
+static double molecule_factor(const double *intensities,const double *phases,
                               struct rvec q,
                               double ax, double ay, double az,
                               double bx, double by, double bz,
@@ -267,18 +230,13 @@ static double molecule_factor(const double *intensities,
 		h = (signed int)rint(hd);
 		k = (signed int)rint(kd);
 		l = (signed int)rint(ld);
-		if ( lookup_count(counts, h, k, l) == 0 ) {
-			ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-			      h, k, l);
-			return 1.0e20;
-		}
 		r = lookup_intensity(intensities, h, k, l);
 		break;
 	case GRADIENT_INTERPOLATE :
-		r = interpolate_intensity(intensities, counts, hd, kd, ld);
+		r = interpolate_intensity(intensities, hd, kd, ld);
 		break;
 	case GRADIENT_PHASED :
-		r = interpolate_phased_intensity(intensities, counts, phases,
+		r = interpolate_phased_intensity(intensities, phases,
 		                                 hd, kd, ld);
 		break;
 	default:
@@ -291,7 +249,7 @@ static double molecule_factor(const double *intensities,
 
 
 double water_diffraction(struct rvec q, double en,
-                                 double beam_r, double water_r)
+                         double beam_r, double water_r)
 {
 	double complex fH, fO;
 	double s, modq;
@@ -390,9 +348,8 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
 
 
 void get_diffraction(struct image *image, int na, int nb, int nc,
-                     const double *intensities, const unsigned int *counts,
-                     const double *phases, UnitCell *cell, int do_water,
-                     GradientMethod m)
+                     const double *intensities, const double *phases,
+                     UnitCell *cell, int do_water, GradientMethod m)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -447,7 +404,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 				I_molecule = 1.0e10;
 			} else {
 				I_molecule = molecule_factor(intensities,
-				                             counts, phases, q,
+				                             phases, q,
 			                                     ax,ay,az,
 			                                     bx,by,bz,cx,cy,cz,
 				                             m);
diff --git a/src/diffraction.h b/src/diffraction.h
index c7afaa3e..ffe3a19d 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -26,8 +26,7 @@ typedef enum {
 } GradientMethod;
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            const double *intensities,
-                            const unsigned int *counts, const double *phases,
+                            const double *intensities,const double *phases,
                             UnitCell *cell, int do_water, GradientMethod m);
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index c98508ec..78d79f76 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -38,8 +38,6 @@ static void show_help(const char *s)
 "      --poisson              Simulate Poisson samples.\n"
 "      --twin                 Generate twinned data.\n"
 "  -o, --output=<filename>    Output filename (default: stdout).\n"
-"      --zone-axis            Generate hk0 intensities only (and add\n"
-"                              Synth2D-style header.\n"
 "  -i, --intensities=<file>   Read intensities from file instead of\n"
 "                              calculating them from scratch.  You might use\n"
 "                              this if you need to apply noise or twinning.\n"
@@ -48,37 +46,6 @@ static void show_help(const char *s)
 }
 
 
-static int template_reflections(const char *filename, unsigned int *counts)
-{
-	char *rval;
-	FILE *fh;
-
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		return 1;
-	}
-
-	do {
-
-		char line[1024];
-		int r;
-		signed int h, k, l;
-
-		rval = fgets(line, 1023, fh);
-
-		r = sscanf(line, "%i %i %i", &h, &k, &l);
-		if ( r != 3 ) continue;
-
-		set_count(counts, h, k, l, 1);
-
-	} while ( rval != NULL );
-
-	fclose(fh);
-
-	return 0;
-}
-
-
 /* Apply Poisson noise to all reflections */
 static void noisify_reflections(double *ref)
 {
@@ -111,13 +78,12 @@ int main(int argc, char *argv[])
 	char *template = NULL;
 	int config_noisify = 0;
 	int config_twin = 0;
-	int config_za = 0;
 	char *output = NULL;
-	unsigned int *counts;
-	unsigned int *cts;
 	char *input = NULL;
 	signed int h, k, l;
 	char *filename = NULL;
+	ReflItemList *input_items;
+	ReflItemList *write_items;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -126,7 +92,6 @@ int main(int argc, char *argv[])
 		{"poisson",            0, &config_noisify,     1},
 		{"output",             1, NULL,               'o'},
 		{"twin",               0, &config_twin,        1},
-		{"zone-axis",          0, &config_za,          1},
 		{"intensities",        1, NULL,               'i'},
 		{"pdb",                1, NULL,               'p'},
 		{0, 0, NULL, 0}
@@ -170,41 +135,18 @@ int main(int argc, char *argv[])
 	}
 
 	mol = load_molecule(filename);
-	cts = new_list_count();
-	phases = new_list_intensity(); /* "intensity" type used for phases */
+	phases = new_list_phase();
 	if ( input == NULL ) {
+		input_items = new_items();
 		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
-		                            cts, phases);
+		                            phases, input_items);
 	} else {
-		ideal_ref = read_reflections(input, cts, phases);
+		ideal_ref = new_list_intensity();
+		phases = new_list_phase();
+		input_items = read_reflections(input, ideal_ref, phases, NULL);
 		free(input);
 	}
 
-	counts = new_list_count();
-
-	if ( template != NULL ) {
-
-		if ( template_reflections(template, counts) != 0 ) {
-			ERROR("Failed to template reflections.\n");
-			return 1;
-		}
-
-	} else {
-
-		/* No template? Then only mark reflections which were
-		 * calculated. */
-		for ( h=-INDMAX; h<=INDMAX; h++ ) {
-		for ( k=-INDMAX; k<=INDMAX; k++ ) {
-		for ( l=-INDMAX; l<=INDMAX; l++ ) {
-			unsigned int c;
-			c = lookup_count(cts, h, k, l);
-			set_count(counts, h, k, l, c);
-		}
-		}
-		}
-
-	}
-
 	if ( config_noisify ) noisify_reflections(ideal_ref);
 
 	if ( config_twin ) {
@@ -239,8 +181,24 @@ int main(int argc, char *argv[])
 
 	}
 
-	write_reflections(output, counts, ideal_ref, phases,
-	                  config_za, mol->cell, 1);
+	if ( template ) {
+		/* Write out only reflections which are in the template
+		 * (and which we have in the input) */
+		ReflItemList *template_items;
+		template_items = read_reflections(template, NULL, NULL, NULL);
+		write_items = intersection_items(input_items, template_items);
+		delete_items(template_items);
+	} else {
+		/* Write out all reflections */
+		write_items = new_items();
+		union_items(write_items, input_items);
+	}
+
+	write_reflections(output, write_items, ideal_ref, phases, NULL,
+	                  mol->cell);
+
+	delete_items(input_items);
+	delete_items(write_items);
 
 	return 0;
 }
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 3a3dfb56..5df08ad3 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -65,7 +65,6 @@ struct process_args
 	struct detector *det;
 	IndexingMethod indm;
 	const double *intensities;
-	const unsigned int *counts;
 	struct gpu_context *gctx;
 
 	/* Thread control and output */
@@ -225,19 +224,18 @@ static struct image *get_simage(struct image *template, int alternate)
 
 
 static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
-                               const double *intensities,
-                               const unsigned int *counts, UnitCell *cell)
+                               const double *intensities, UnitCell *cell)
 {
 	/* Set up GPU if necessary */
 	if ( (gctx != NULL) && (*gctx == NULL) ) {
-		*gctx = setup_gpu(0, simage, intensities, counts);
+		*gctx = setup_gpu(0, simage, intensities);
 	}
 
 	if ( (gctx != NULL) && (*gctx != NULL) ) {
 		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
 	} else {
 		get_diffraction(simage, 24, 24, 40,
-		                intensities, counts, NULL, cell, 0,
+		                intensities, NULL, cell, 0,
 		                GRADIENT_MOSAIC);
 	}
 	record_image(simage, 0);
@@ -270,7 +268,6 @@ static struct process_result process_image(struct process_args *pargs)
 	int config_polar = pargs->config_polar;
 	IndexingMethod indm = pargs->indm;
 	const double *intensities = pargs->intensities;
-	const unsigned int *counts = pargs->counts;
 	struct gpu_context *gctx = pargs->gctx;
 
 	image.features = NULL;
@@ -368,11 +365,11 @@ static struct process_result process_image(struct process_args *pargs)
 		if ( config_gpu ) {
 			pthread_mutex_lock(pargs->gpu_mutex);
 			simulate_and_write(simage, &gctx, intensities,
-			                   counts, image.indexed_cell);
+			                   image.indexed_cell);
 			pthread_mutex_unlock(pargs->gpu_mutex);
 		} else {
 			simulate_and_write(simage, NULL, intensities,
-			                   counts, image.indexed_cell);
+			                   image.indexed_cell);
 		}
 	}
 
@@ -469,7 +466,6 @@ int main(int argc, char *argv[])
 	UnitCell *cell;
 	double *intensities = NULL;
 	char *intfile = NULL;
-	unsigned int *counts;
 	char *pdb = NULL;
 	char *prefix = NULL;
 	int nthreads = 1;
@@ -568,11 +564,11 @@ int main(int argc, char *argv[])
 	free(filename);
 
 	if ( intfile != NULL ) {
-		counts = new_list_count();
-		intensities = read_reflections(intfile, counts, NULL);
+		ReflItemList *items;
+		items = read_reflections(intfile, intensities, NULL, NULL);
+		delete_items(items);
 	} else {
 		intensities = NULL;
-		counts = NULL;
 	}
 
 	if ( pdb == NULL ) {
@@ -671,7 +667,6 @@ int main(int argc, char *argv[])
 		pargs->det = det;
 		pargs->indm = indm;
 		pargs->intensities = intensities;
-		pargs->counts = counts;
 		pargs->gctx = gctx;
 		pargs->id = i;
 		pthread_mutex_lock(&pargs->control_mutex);
diff --git a/src/likelihood.c b/src/likelihood.c
index b6a37994..cb5595f8 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -18,17 +18,17 @@
 #include "utils.h"
 
 
-void scale_intensities(const double *model, double *new_pattern,
-                       const unsigned int *model_counts,
-                       ReflItemList *items, double f0, int f0_valid)
+void scale_intensities(const double *model, ReflItemList *model_items,
+                       double *new_pattern, ReflItemList *new_items,
+                       double f0, int f0_valid)
 {
 	double s;
 	unsigned int i;
-	unsigned int *new_counts;
+	ReflItemList *items;
 
-	new_counts = items_to_counts(items);
-
-	s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
+	items = intersection_items(model_items, new_items);
+	s = stat_scale_intensity(model, new_pattern, items);
+	delete_items(items);
 	if ( f0_valid ) printf("%f %f\n", s, f0);
 
 	/* NaN -> abort */
@@ -38,6 +38,4 @@ void scale_intensities(const double *model, double *new_pattern,
 	for ( i=0; i<LIST_SIZE; i++ ) {
 		new_pattern[i] *= s;
 	}
-
-	free(new_counts);
 }
diff --git a/src/likelihood.h b/src/likelihood.h
index 08956fe3..70aef098 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -20,10 +20,9 @@
 #include "utils.h"
 
 
-extern void scale_intensities(const double *model, double *new_pattern,
-	                      const unsigned int *model_counts,
-                              ReflItemList *items, double f0,
-                              int f0_valid);
+extern void scale_intensities(const double *model, ReflItemList *model_items,
+                              double *new_pattern, ReflItemList *new_items,
+                              double f0, int f0_valid);
 
 
 #endif	/* LIKELIHOOD_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 5d59d71b..616eec2c 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -195,7 +195,6 @@ int main(int argc, char *argv[])
 	int n_images = 1; /* Generate one image by default */
 	int done = 0;
 	UnitCell *cell;
-	unsigned int *counts;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -323,17 +322,19 @@ int main(int argc, char *argv[])
 		STATUS("reflection intensities (with --intensities).\n");
 		STATUS("I'll simulate a flat intensity distribution.\n");
 		intensities = NULL;
-		counts = NULL;
 		phases = NULL;
 	} else {
-		counts = new_list_count();
+		ReflItemList *items;
 		if ( grad == GRADIENT_PHASED ) {
 			phases = new_list_phase();
 		} else {
 			phases = NULL;
 		}
-		intensities = read_reflections(intfile, counts, phases);
+		intensities = new_list_intensity();
+		phases = new_list_phase();
+		items = read_reflections(intfile, intensities, phases, NULL);
 		free(intfile);
+		delete_items(items);
 	}
 
 	/* Define image parameters */
@@ -401,12 +402,12 @@ int main(int argc, char *argv[])
 		if ( config_gpu ) {
 			if ( gctx == NULL ) {
 				gctx = setup_gpu(config_nosfac, &image,
-				                 intensities, counts);
+				                 intensities);
 			}
 			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
 		} else {
-			get_diffraction(&image, na, nb, nc, intensities, counts,
-			                phases, cell, !config_nowater, grad);
+			get_diffraction(&image, na, nb, nc, intensities, phases,
+			                cell, !config_nowater, grad);
 		}
 		if ( image.data == NULL ) {
 			ERROR("Diffraction calculation failed.\n");
diff --git a/src/process_hkl.c b/src/process_hkl.c
index b5b6df9a..986477fb 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -89,6 +89,8 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
 	twins = get_twins(test_items, holo, mero);
 	delete_items(test_items);
 
+	/* Idiot check.  Wouldn't be necessary if I could prove that the above
+	 * set of arbitrarily chosen reflections were always general. */
 	if ( num_items(twins) != np ) {
 		ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
 		abort();
@@ -98,10 +100,9 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
 }
 
 
-static int resolve_twin(ReflItemList *items, ReflItemList *twins,
-                        const double *patt, const double *model,
-                        const unsigned int *model_counts,
-                        const char *holo, const char *mero)
+static int resolve_twin(const double *model, ReflItemList *observed,
+                        const double *patt, ReflItemList *items,
+                        ReflItemList *twins, const char *holo, const char *mero)
 {
 	int n, i;
 	double best_fom = 0.0;
@@ -116,6 +117,7 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
 		double *trial_ints = new_list_intensity();
 		unsigned int *trial_counts = new_list_count();
 		double fom;
+		ReflItemList *intersection;
 
 		op = get_item(twins, i)->op;
 
@@ -134,8 +136,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
 
 		}
 
-		fom = stat_pearson(trial_ints, trial_counts,
-		                   model, model_counts);
+		intersection = intersection_items(observed, items);
+		fom = stat_pearson(trial_ints, model, intersection);
+		delete_items(intersection);
 
 		free(trial_ints);
 		free(trial_counts);
@@ -153,9 +156,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
 }
 
 
-static void merge_pattern(double *model, const double *new,
-                          unsigned int *model_counts,
-                          ReflItemList *items, int mo, int sum,
+static void merge_pattern(double *model, ReflItemList *observed,
+                          const double *new,  ReflItemList *items,
+                          unsigned int *model_counts,  int mo,
                           ReflItemList *twins,
                           const char *holo, const char *mero)
 {
@@ -163,8 +166,8 @@ static void merge_pattern(double *model, const double *new,
 	int twin;
 
 	if ( twins != NULL ) {
-		twin = resolve_twin(items, twins, new, model,
-		                     model_counts, holo, mero);
+		twin = resolve_twin(model, observed, new, items,
+		                    twins, holo, mero);
 	} else {
 		twin = 0;
 	}
@@ -187,25 +190,23 @@ static void merge_pattern(double *model, const double *new,
 
 		intensity = lookup_intensity(new, h, k, l);
 
+		/* User asked for max only? */
 		if ( !mo ) {
 			integrate_intensity(model, h, k, l, intensity);
-			if ( sum ) {
-				set_count(model_counts, h, k, l, 1);
-			} else {
-				integrate_count(model_counts, h, k, l, 1);
-			}
 		} else {
 			if ( intensity > lookup_intensity(model, h, k, l) ) {
 				set_intensity(model, h, k, l, intensity);
 			}
-			set_count(model_counts, h, k, l, 1);
 		}
 
+		integrate_count(model_counts, h, k, l, 1);
+
 	}
 }
 
 
-static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
+static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
+                      unsigned int **pcounts,
                       int config_maxonly, int config_scale, int config_sum,
                       int config_stopafter,
                       ReflItemList *twins, const char *holo, const char *mero,
@@ -218,6 +219,10 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
 	int n_patterns = 0;
 	double *new_pattern = new_list_intensity();
 	ReflItemList *items = new_items();
+	ReflItemList *observed = new_items();
+	double *model = new_list_intensity();
+	unsigned int *counts = new_list_count();
+	int i;
 
 	do {
 
@@ -244,19 +249,20 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
 
 			/* Scale if requested */
 			if ( config_scale ) {
-				scale_intensities(model, new_pattern,
-				                  model_counts, items, f0,
-				                  f0_valid);
+				scale_intensities(model, observed,
+				                  new_pattern, items,
+				                  f0, f0_valid);
 			}
 
 			/* Start of second or later pattern */
-			merge_pattern(model, new_pattern, model_counts,
-			              items, config_maxonly, config_sum, twins,
-			              holo, mero);
+			merge_pattern(model, observed, new_pattern, items,
+			              counts, config_maxonly,
+			              twins, holo, mero);
 
 			if ( n_patterns == config_stopafter ) break;
 
 			n_patterns++;
+			union_items(observed, items);
 			clear_items(items);
 
 			progress_bar(n_patterns, n_total_patterns, "Merging");
@@ -292,6 +298,19 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
 	delete_items(items);
 	free(new_pattern);
 
+	if ( !config_sum ) {
+		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+			if ( counts[i] > 0 ) {
+				model[i] /= (double)counts[i];
+				counts[i] = 1;
+			}
+		}
+	}
+
+	*pmodel = model;
+	*pcounts = counts;
+	*pobserved = observed;
+
 	STATUS("%i patterns had no f0 valid value.\n", n_nof0);
 }
 
@@ -322,7 +341,7 @@ int main(int argc, char *argv[])
 	char *output = NULL;
 	FILE *fh;
 	double *model;
-	unsigned int *model_counts;
+	unsigned int *counts;
 	UnitCell *cell;
 	int config_maxonly = 0;
 	int config_stopafter = 0;
@@ -332,6 +351,7 @@ int main(int argc, char *argv[])
 	char *sym = NULL;
 	char *pdb = NULL;
 	ReflItemList *twins;
+	ReflItemList *observed;
 	int i;
 
 	/* Long options */
@@ -400,9 +420,6 @@ int main(int argc, char *argv[])
 		sym = strdup("1");
 	}
 
-	model = new_list_intensity();
-	model_counts = new_list_count();
-
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
@@ -443,7 +460,7 @@ int main(int argc, char *argv[])
 	STATUS("There are %i patterns to process\n", n_total_patterns);
 	rewind(fh);
 
-	merge_all(fh, model, model_counts,
+	merge_all(fh, &model, &observed, &counts,
 	          config_maxonly, config_scale, config_sum, config_stopafter,
                   twins, holo, sym, n_total_patterns);
 	rewind(fh);
@@ -451,13 +468,12 @@ int main(int argc, char *argv[])
 	fclose(fh);
 
 	if ( output != NULL ) {
-		write_reflections(output, model_counts, model, NULL,
-		                  0, cell, 1);
+		write_reflections(output, observed, model, NULL, counts, cell);
 	}
 
 	free(sym);
 	free(model);
-	free(model_counts);
+	free(counts);
 	free(output);
 	free(cell);
 
diff --git a/src/reflections.c b/src/reflections.c
index bd9b446d..0ab7a91e 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -21,12 +21,12 @@
 #include "reflections.h"
 
 
-void write_reflections(const char *filename, unsigned int *counts,
-                       double *ref, double *phases, int zone_axis,
-                       UnitCell *cell, unsigned int min_counts)
+void write_reflections(const char *filename, ReflItemList *items,
+                       double *intensities, double *phases,
+                       unsigned int *counts, UnitCell *cell)
 {
 	FILE *fh;
-	signed int h, k, l;
+	int i;
 
 	if ( filename == NULL ) {
 		fh = stdout;
@@ -40,23 +40,16 @@ void write_reflections(const char *filename, unsigned int *counts,
 	}
 
 	/* Write spacings and angle if zone axis pattern */
-	if ( zone_axis ) {
-		double a, b, c, alpha, beta, gamma;
-		cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-		fprintf(fh, "a %5.3f nm\n",
-		            (0.5/resolution(cell, 1, 0, 0))*1e9);
-		fprintf(fh, "b %5.3f nm\n",
-		            (0.5/resolution(cell, 0, 1, 0))*1e9);
-		fprintf(fh, "angle %5.3f deg\n", rad2deg(alpha));
-		fprintf(fh, "scale 10\n");
-	} else {
-		fprintf(fh, "  h   k   l          I    phase   sigma(I) "
-		            " 1/d(nm^-1)  counts\n");
-	}
+	fprintf(fh, "  h   k   l          I    phase   sigma(I) "
+		    " 1/d(nm^-1)  counts\n");
+
+	for ( i=0; i<num_items(items); i++ ) {
+
+		struct refl_item *it;
+		signed int h, k, l;
 
-	for ( h=-INDMAX; h<INDMAX; h++ ) {
-	for ( k=-INDMAX; k<INDMAX; k++ ) {
-	for ( l=-INDMAX; l<INDMAX; l++ ) {
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
 
 		int N;
 		double intensity, s;
@@ -64,14 +57,12 @@ void write_reflections(const char *filename, unsigned int *counts,
 
 		if ( counts ) {
 			N = lookup_count(counts, h, k, l);
-			if ( N < min_counts ) continue;
 		} else {
 			N = 1;
 		}
-		if ( zone_axis && (l != 0) ) continue;
 
-		/* Divide measured intensity by the number of counts */
-		intensity = lookup_intensity(ref, h, k, l) / N;
+		intensity = lookup_intensity(intensities, h, k, l);
+
 		if ( phases != NULL ) {
 			double p;
 			p = lookup_phase(phases, h, k, l);
@@ -91,18 +82,25 @@ void write_reflections(const char *filename, unsigned int *counts,
 		        h, k, l, intensity, ph, 0.0, s/1.0e9, N);
 
 	}
-	}
-	}
 	fclose(fh);
 }
 
 
-double *read_reflections(const char *filename, unsigned int *counts,
-                         double *phases)
+/* Read reflections from file.  Returns the list of reflections successfully
+ * read in.  "intensities", "phases" and "counts" are lists which will be
+ * populated with the values read from the file.  Existing values in either list
+ * will be overwritten if the reflection is read from the file, but other values
+ * will be left intact.
+ *
+ * "intensities", "phases" or "counts" can be NULL, if you don't need them.
+ */
+ReflItemList *read_reflections(const char *filename,
+                               double *intensities, double *phases,
+                               unsigned int *counts)
 {
-	double *ref;
 	FILE *fh;
 	char *rval;
+	ReflItemList *items;
 
 	fh = fopen(filename, "r");
 	if ( fh == NULL ) {
@@ -110,7 +108,7 @@ double *read_reflections(const char *filename, unsigned int *counts,
 		return NULL;
 	}
 
-	ref = new_list_intensity();
+	items = new_items();
 
 	do {
 
@@ -147,24 +145,24 @@ double *read_reflections(const char *filename, unsigned int *counts,
 			continue;
 		}
 
-		set_intensity(ref, h, k, l, intensity);
+		add_item(items, h, k, l);
+
+		if ( intensities != NULL ) {
+			set_intensity(intensities, h, k, l, intensity);
+		}
 		if ( phases != NULL ) {
 			ph = atof(phs);
 			set_phase(phases, h, k, l, ph);
 		}
 		if ( counts != NULL ) {
 			set_count(counts, h, k, l, cts);
-			/* In this case, the intensity must be multiplied up
-			 * because other parts of the program will try to
-			 * divide it down. */
-			set_intensity(ref, h, k, l, intensity*(double)cts);
 		}
 
 	} while ( rval != NULL );
 
 	fclose(fh);
 
-	return ref;
+	return items;
 }
 
 
@@ -188,16 +186,3 @@ double *ideal_intensities(double complex *sfac)
 
 	return ref;
 }
-
-
-void divide_down(double *intensities, unsigned int *counts)
-{
-	int i;
-
-	for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-		if ( counts[i] > 0 ) {
-			intensities[i] /= (double)counts[i];
-			counts[i] = 1;
-		}
-	}
-}
diff --git a/src/reflections.h b/src/reflections.h
index 84946f1a..de5feb72 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -18,18 +18,18 @@
 
 
 #include "cell.h"
+#include "utils.h"
 
 
-extern void write_reflections(const char *filename, unsigned int *counts,
-                              double *ref, double *phases, int zone_axis,
-                              UnitCell *cell, unsigned int min_counts);
+extern void write_reflections(const char *filename, ReflItemList *items,
+                              double *intensities, double *phases,
+                              unsigned int *counts, UnitCell *cell);
 
-extern double *read_reflections(const char *filename, unsigned int *counts,
-                                double *phases);
+extern ReflItemList *read_reflections(const char *filename,
+                                      double *intensities, double *phases,
+                                      unsigned int *counts);
 
 extern double *ideal_intensities(double complex *sfac);
 
-extern void divide_down(double *intensities, unsigned int *counts);
-
 
 #endif	/* REFLECTIONS_H */
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 60153652..829ceb4e 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -128,10 +128,10 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
 
 		switch ( wght ) {
 		case WGHT_I :
-			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			val = lookup_intensity(ref, h, k, 0);
 			break;
 		case WGHT_SQRTI :
-			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			val = lookup_intensity(ref, h, k, 0);
 			val = (val>0.0) ? sqrt(val) : 0.0;
 			break;
 		case WGHT_COUNTS :
@@ -290,7 +290,6 @@ int main(int argc, char *argv[])
 	UnitCell *cell;
 	char *infile;
 	double *ref;
-	unsigned int *cts;
 	int config_povray = 0;
 	int config_zoneaxis = 0;
 	int config_sqrt = 0;
@@ -303,6 +302,7 @@ int main(int argc, char *argv[])
 	int wght;
 	int colscale;
 	char *cscale = NULL;
+	unsigned int *cts;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -417,8 +417,10 @@ int main(int argc, char *argv[])
 		ERROR("Couldn't load unit cell from %s\n", pdb);
 		return 1;
 	}
+	ref = new_list_intensity();
 	cts = new_list_count();
-	ref = read_reflections(infile, cts, NULL);
+	ReflItemList *items = read_reflections(infile, ref, NULL, cts);
+	delete_items(items);
 	if ( ref == NULL ) {
 		ERROR("Couldn't open file '%s'\n", infile);
 		return 1;
diff --git a/src/sfac.c b/src/sfac.c
index a8b94051..b5ae7ace 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -491,7 +491,7 @@ void free_molecule(struct molecule *mol)
 
 
 double *get_reflections(struct molecule *mol, double en, double res,
-                        unsigned int *counts, double *phases)
+                        double *phases, ReflItemList *items)
 {
 	double *reflections;
 	double asx, asy, asz;
@@ -556,7 +556,7 @@ double *get_reflections(struct molecule *mol, double en, double res,
 
 		set_intensity(reflections, h, k, l, pow(cabs(F), 2.0));
 		if ( phases != NULL ) set_phase(phases, h, k, l, carg(F));
-		set_count(counts, h, k, l, 1);
+		if ( items != NULL ) add_item(items, h, k, l);
 
 	}
 	progress_bar((k+INDMAX)+IDIM*(h+INDMAX), IDIM*IDIM-1,
diff --git a/src/sfac.h b/src/sfac.h
index ea30c83e..115dcb5f 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -19,6 +19,7 @@
 #include <complex.h>
 
 #include "cell.h"
+#include "utils.h"
 
 
 #define MAX_ATOMS (128*1024)
@@ -61,7 +62,7 @@ extern double complex get_sfac(const char *n, double s, double en);
 extern struct molecule *load_molecule(const char *filename);
 extern void free_molecule(struct molecule *mol);
 extern double *get_reflections(struct molecule *mol, double en, double res,
-                               unsigned int *counts, double *phases);
+                               double *phases, ReflItemList *items);
 
 
 #endif	/* SFAC_H */
diff --git a/src/statistics.c b/src/statistics.c
index 889ade51..1acbb718 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -26,10 +26,9 @@
 
 struct r_params {
 	const double *ref1;
-	const unsigned int *c1;
 	const double *ref2;
-	const unsigned int *c2;
-	int fom;
+	ReflItemList *items;  /* Which reflections to use */
+	int fom;              /* Which FoM to use (see the enum just below) */
 };
 
 enum {
@@ -38,23 +37,31 @@ enum {
 };
 
 
-double stat_scale_intensity(const double *ref1, const unsigned int *c1,
-                            const double *ref2, const unsigned int *c2)
+/* Return the least squares optimal scaling factor when comparing intensities.
+ * ref1,ref2 are the two intensity lists to compare.  "items" is a ReflItemList
+ * containing the reflections which should be taken into account.
+ */
+double stat_scale_intensity(const double *ref1, const double *ref2,
+                            ReflItemList *items)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	int i;
 
-	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<LIST_SIZE; i++ ) {
+	for ( i=0; i<num_items(items); i++ ) {
+
+		double i1, i2;
+		struct refl_item *it;
+		signed int h, k, l;
+
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
+
+		i1 = lookup_intensity(ref1, h, k, l);
+		i2 = lookup_intensity(ref2, h, k, l);
 
-		if ( c1[i] && c2[i] ) {
-			double i1, i2;
-			i1 = ref1[i] / (double)c1[i];
-			i2 = ref2[i] / (double)c2[i];
-			top += i1 * i2;
-			bot += i2 * i2;
-		} /* else reflection not common so don't worry about it */
+		top += i1 * i2;
+		bot += i2 * i2;
 
 	}
 
@@ -62,29 +69,36 @@ double stat_scale_intensity(const double *ref1, const unsigned int *c1,
 }
 
 
-double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
-                            const double *ref2, const unsigned int *c2)
+/* Return the least squares optimal scaling factor when comparing the square
+ * roots of the intensities (i.e. one approximation to the structure factor
+ * moduli).
+ * ref1,ref2 are the two intensity lists to compare (they contain intensities,
+ * not square rooted intensities).  "items" is a ReflItemList containing the
+ * reflections which should be taken into account.
+ */
+static double stat_scale_sqrti(const double *ref1, const double *ref2,
+                               ReflItemList *items)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	int i;
 
-	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<LIST_SIZE; i++ ) {
-
-		if ( c1[i] && c2[i] ) {
+	for ( i=0; i<num_items(items); i++ ) {
 
-			double f1, f2;
+		double i1, i2, f1, f2;
+		struct refl_item *it;
+		signed int h, k, l;
 
-			if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
 
-			f1 = sqrt(ref1[i]) / (double)c1[i];
-			f2 = sqrt(ref2[i]) / (double)c2[i];
-
-			top += f1 * f2;
-			bot += f2 * f2;
+		i1 = lookup_intensity(ref1, h, k, l);
+		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+		i2 = lookup_intensity(ref2, h, k, l);
+		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 
-		} /* else reflection not common so don't worry about it */
+		top += i1 * i2;
+		bot += i2 * i2;
 
 	}
 
@@ -92,27 +106,27 @@ double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
 }
 
 
-static double internal_r2(const double *ref1, const unsigned int *c1,
-                          const double *ref2, const unsigned int *c2,
-                          double scale)
+static double internal_r2(const double *ref1, const double *ref2,
+                          ReflItemList *items, double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	int i;
 
-	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<LIST_SIZE; i++ ) {
+	for ( i=0; i<num_items(items); i++ ) {
 
-		if ( c1[i] && c2[i] ) {
+		double i1, i2;
+		struct refl_item *it;
+		signed int h, k, l;
 
-			double i1, i2;
-			i1 = ref1[i] / (scale*(double)c1[i]);
-			i2 = ref2[i] / (double)c2[i];
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
 
-			top += pow(i1 - i2, 2.0);
-			bot += pow(i1, 2.0);
+		i1 = scale * lookup_intensity(ref1, h, k, l);
+		i2 = lookup_intensity(ref2, h, k, l);
 
-		} /* else reflection not measured so don't worry about it */
+		top += pow(i1 - i2, 2.0);
+		bot += pow(i1, 2.0);
 
 	}
 
@@ -120,30 +134,29 @@ static double internal_r2(const double *ref1, const unsigned int *c1,
 }
 
 
-static double internal_rmerge(const double *ref1, const unsigned int *c1,
-                              const double *ref2, const unsigned int *c2,
-                              double scale)
+static double internal_rmerge(const double *ref1, const double *ref2,
+                              ReflItemList *items, double scale)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	int i;
 
-	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<LIST_SIZE; i++ ) {
-
-		if ( c1[i] && c2[i] ) {
+	for ( i=0; i<num_items(items); i++ ) {
 
-			double f1, f2;
+		double i1, i2, f1, f2;
+		struct refl_item *it;
+		signed int h, k, l;
 
-			if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
 
-			f1 = sqrt(ref1[i]) / (scale*(double)c1[i]);
-			f2 = sqrt(ref2[i]) / (double)c2[i];
-
-			top += fabs(f1 - f2);
-			bot += f1 + f2;
+		i1 = lookup_intensity(ref1, h, k, l);
+		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+		i2 = lookup_intensity(ref2, h, k, l);
+		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 
-		} /* else reflection not measured so don't worry about it */
+		top += fabs(f1 - f2);
+		bot += f1 + f2;
 
 	}
 
@@ -157,11 +170,9 @@ static double calc_r(double scale, void *params)
 
 	switch ( rp->fom) {
 	case R_MERGE :
-		return internal_rmerge(rp->ref1, rp->c1,
-		                       rp->ref2, rp->c2, scale);
+		return internal_rmerge(rp->ref1, rp->ref2, rp->items, scale);
 	case R_2 :
-		return internal_r2(rp->ref1, rp->c1,
-		                   rp->ref2, rp->c2, scale);
+		return internal_r2(rp->ref1, rp->ref2, rp->items, scale);
 	}
 
 	ERROR("No such FoM!\n");
@@ -169,9 +180,8 @@ static double calc_r(double scale, void *params)
 }
 
 
-static double r_minimised(const double *ref1, const unsigned int *c1,
-                          const double *ref2, const unsigned int *c2,
-                          double *scalep, int fom)
+static double r_minimised(const double *ref1, const double *ref2,
+                          ReflItemList *items, double *scalep, int fom)
 {
 	gsl_function F;
 	gsl_min_fminimizer *s;
@@ -182,8 +192,7 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
 
 	rp.ref1 = ref1;
 	rp.ref2 = ref2;
-	rp.c1 = c1;
-	rp.c2 = c2;
+	rp.items = items;
 	rp.fom = fom;
 
 	F.function = &calc_r;
@@ -194,10 +203,10 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
 	/* Initial guess */
 	switch ( fom ) {
 	case R_MERGE :
-		scale = stat_scale_sqrti(ref1, c1, ref2, c2);
+		scale = stat_scale_sqrti(ref1, ref2, items);
 		break;
 	case R_2 :
-		scale = stat_scale_intensity(ref1, c1, ref2, c2);
+		scale = stat_scale_intensity(ref1, ref2, items);
 		break;
 	}
 	//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -236,52 +245,52 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
 }
 
 
-double stat_rmerge(const double *ref1, const unsigned int *c1,
-                   const double *ref2, const unsigned int *c2,
-                   double *scalep)
+double stat_rmerge(const double *ref1, const double *ref2,
+                   ReflItemList *items, double *scalep)
 {
-	return r_minimised(ref1, c1, ref2, c2, scalep, R_MERGE);
+	return r_minimised(ref1, ref2, items, scalep, R_MERGE);
 }
 
 
-double stat_r2(const double *ref1, const unsigned int *c1,
-               const double *ref2, const unsigned int *c2,
-               double *scalep)
+double stat_r2(const double *ref1, const double *ref2,
+               ReflItemList *items, double *scalep)
 {
-	return r_minimised(ref1, c1, ref2, c2, scalep, R_2);
+	return r_minimised(ref1, ref2, items, scalep, R_2);
 }
 
 
-double stat_pearson(const double *ref1, const unsigned int *c1,
-                    const double *ref2, const unsigned int *c2)
+double stat_pearson(const double *ref1, const double *ref2, ReflItemList *items)
 {
-	double vec1[4096];
-	double vec2[4096];
-	signed int h, k, l;
+	double *vec1, *vec2;
 	int i = 0;
+	int ni = num_items(items);
+	double val;
 
-	for ( l=-INDMAX; l<INDMAX; l++ ) {
-	for ( k=-INDMAX; k<INDMAX; k++ ) {
-	for ( h=-INDMAX; h<INDMAX; h++ ) {
+	vec1 = malloc(ni*sizeof(double));
+	vec2 = malloc(ni*sizeof(double));
 
-		double i1, i2;
-		unsigned int c1s, c2s;
+	for ( i=0; i<ni; i++ ) {
+
+		double i1, i2, f1, f2;
+		struct refl_item *it;
+		signed int h, k, l;
 
-		c1s = lookup_count(c1, h, k, l);
-		c2s = lookup_count(c2, h, k, l);
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
 
 		i1 = lookup_intensity(ref1, h, k, l);
+		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
 		i2 = lookup_intensity(ref2, h, k, l);
+		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
 
-		if ( c1s && c2s && (i1>0.0) && (i2>0.0) ) {
-			vec1[i] = sqrt(i1 / (double)c1s);
-			vec2[i] = sqrt(i2 / (double)c2s);
-			i++;
-		}
+		vec1[i] = f1;
+		vec2[i] = f2;
 
 	}
-	}
-	}
 
-	return gsl_stats_correlation(vec1, 1, vec2, 1, i);
+	val = gsl_stats_correlation(vec1, 1, vec2, 1, i);
+	free(vec1);
+	free(vec2);
+
+	return val;
 }
diff --git a/src/statistics.h b/src/statistics.h
index 5e3ca52c..c6c2f09c 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -16,18 +16,20 @@
 #ifndef STATISTICS_H
 #define STATISTICS_H
 
-extern double stat_scale_intensity(const double *ref1, const unsigned int *c1,
-                                   const double *ref2, const unsigned int *c2);
 
-extern double stat_r2(const double *ref1, const unsigned int *c1,
-                      const double *ref2, const unsigned int *c2,
-                      double *scalep);
+#include "utils.h"
 
-extern double stat_rmerge(const double *ref1, const unsigned int *c1,
-                          const double *ref2, const unsigned int *c2,
-                          double *scalep);
+extern double stat_scale_intensity(const double *ref1, const double *ref2,
+                                   ReflItemList *items);
+
+extern double stat_rmerge(const double *ref1, const double *ref2,
+                          ReflItemList *items, double *scalep);
+
+extern double stat_r2(const double *ref1, const double *ref2,
+                      ReflItemList *items, double *scalep);
+
+extern double stat_pearson(const double *ref1, const double *ref2,
+                           ReflItemList *items);
 
-extern double stat_pearson(const double *ref1, const unsigned int *c1,
-                           const double *ref2, const unsigned int *c2);
 
 #endif	/* STATISTICS_H */
diff --git a/src/utils.c b/src/utils.c
index 21ba882f..0f20435e 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -387,24 +387,23 @@ int num_items(const ReflItemList *items)
 }
 
 
-unsigned int *items_to_counts(ReflItemList *items)
+void union_op_items(ReflItemList *items, ReflItemList *newi)
 {
-	unsigned int *c;
-	int i;
+	int n, i;
 
-	c = new_list_count();
+	n = num_items(newi);
+	for ( i=0; i<n; i++ ) {
 
-	for ( i=0; i<num_items(items); i++ ) {
-		struct refl_item *r;
-		r = get_item(items, i);
-		set_count(c, r->h, r->k, r->l, 1);
-	}
+		struct refl_item *r = get_item(newi, i);
+		if ( find_op(items, r->op) ) continue;
+
+		add_item_with_op(items, r->h, r->k, r->l, r->op);
 
-	return c;
+	}
 }
 
 
-void union_op_items(ReflItemList *items, ReflItemList *newi)
+void union_items(ReflItemList *items, ReflItemList *newi)
 {
 	int n, i;
 
@@ -412,9 +411,28 @@ void union_op_items(ReflItemList *items, ReflItemList *newi)
 	for ( i=0; i<n; i++ ) {
 
 		struct refl_item *r = get_item(newi, i);
-		if ( find_op(items, r->op) ) continue;
+		if ( find_item(items, r->h, r->k, r->l) ) continue;
 
 		add_item_with_op(items, r->h, r->k, r->l, r->op);
 
 	}
 }
+
+
+ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2)
+{
+	int n, i;
+	ReflItemList *res = new_items();
+
+	n = num_items(i1);
+	for ( i=0; i<n; i++ ) {
+
+		struct refl_item *r = get_item(i1, i);
+		if ( find_item(i2, r->h, r->k, r->l) ) {
+			add_item_with_op(res, r->h, r->k, r->l, r->op);
+		}
+
+	}
+
+	return res;
+}
diff --git a/src/utils.h b/src/utils.h
index a7433d28..f7e98469 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -200,6 +200,8 @@ extern struct refl_item *get_item(ReflItemList *items, int i);
 extern int num_items(const ReflItemList *items);
 extern unsigned int *items_to_counts(ReflItemList *items);
 extern void union_op_items(ReflItemList *items, ReflItemList *newi);
+extern void union_items(ReflItemList *items, ReflItemList *newi);
+extern ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2);
 
 /* ------------------------------ Message macros ---------------------------- */