Start work on binary tree

15 files changed, 419 insertions, 306 deletions

diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 2252842f..9aeae803 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -168,8 +168,7 @@ static void add_image(void *args, int cookie)
 	image.flags = NULL;
 	image.indexed_cell = NULL;
 	image.filename = pargs->filename;
-	image.cpeaks = NULL;
-	image.n_cpeaks = 0;
+	image.reflections = NULL;
 	image.det = NULL;
 
 	STATUS("%3i: Processing '%s'\n", cookie, pargs->filename);
diff --git a/src/cubeit.c b/src/cubeit.c
index 808e2e23..7b40f8e5 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -218,8 +218,7 @@ static void sum_image(void *pg, int cookie)
 	image.flags = NULL;
 	image.indexed_cell = NULL;
 	image.filename = apargs->filename;
-	image.cpeaks = NULL;
-	image.n_cpeaks = 0;
+	image.reflections = NULL;
 	image.det = pargs->det;
 
 	STATUS("Processing '%s'\n", apargs->filename);
diff --git a/src/geometry.c b/src/geometry.c
index 9bb41d12..9f9b953c 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -22,6 +22,7 @@
 #include "image.h"
 #include "peaks.h"
 #include "beam-parameters.h"
+#include "reflist.h"
 
 
 #define MAX_CPEAKS (256 * 256)
@@ -117,7 +118,7 @@ static double partiality(double r1, double r2, double r)
 
 
 static int check_reflection(struct image *image, double mres, int output,
-                            struct cpeak *cpeaks, int np,
+                            RefList *reflections,
                             signed int h, signed int k, signed int l,
                             double asx, double asy, double asz,
                             double bsx, double bsy, double bsz,
@@ -136,6 +137,7 @@ static int check_reflection(struct image *image, double mres, int output,
 	double divergence = image->div;
 	double lambda = image->lambda;
 	double klow, kcen, khigh;    /* Wavenumber */
+	Reflection *refl;
 
 	/* "low" gives the largest Ewald sphere,
 	 * "high" gives the smallest Ewald sphere. */
@@ -207,17 +209,9 @@ static int check_reflection(struct image *image, double mres, int output,
 	if ( p == -1 ) return 0;
 
 	/* Add peak to list */
-	cpeaks[np].h = h;
-	cpeaks[np].k = k;
-	cpeaks[np].l = l;
-	cpeaks[np].x = xda;
-	cpeaks[np].y = yda;
-	cpeaks[np].r1 = rlow;
-	cpeaks[np].r2 = rhigh;
-	cpeaks[np].p = part;
-	cpeaks[np].clamp1 = clamp_low;
-	cpeaks[np].clamp2 = clamp_high;
-	np++;
+	refl = add_refl(reflections, h, k, l);
+	set_detector_pos(refl, 0.0, xda, yda);
+	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
 
 	if ( output ) {
 		printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
@@ -228,23 +222,18 @@ static int check_reflection(struct image *image, double mres, int output,
 }
 
 
-struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                 int *n, int output)
+RefList *find_intersections(struct image *image, UnitCell *cell,
+                            int output)
 {
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
-	struct cpeak *cpeaks;
-	int np = 0;
+	RefList *reflections;
 	int hmax, kmax, lmax;
 	double mres;
 	signed int h, k, l;
 
-	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
-	if ( cpeaks == NULL ) {
-		*n = 0;
-		return NULL;
-	}
+	reflections = reflist_new();
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
 	                          &bsx, &bsy, &bsz,
@@ -260,29 +249,30 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	for ( l=-lmax; l<lmax; l++ ) {
 		/* Ignore central beam */
 		if ( (h==0) && (k==0) && (l==0) ) continue;
-		np += check_reflection(image, mres, output, cpeaks, np, h, k, l,
-		                       asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
-		if ( np == MAX_CPEAKS ) goto out;
+		check_reflection(image, mres, output, reflections, h, k, l,
+		                 asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
 	}
 	}
 	}
 
-out:
-	*n = np;
-	return cpeaks;
+	return reflections;
 }
 
 
-double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
+double integrate_all(struct image *image, RefList *reflections)
 {
 	double itot = 0.0;
-	int i;
+	Reflection *refl;
 
-	for ( i=0; i<n; i++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		float x, y, intensity;
+		double xp, yp;
+		get_detector_pos(refl, &xp, &yp);
 
-		if ( integrate_peak(image, cpeaks[i].x, cpeaks[i].y, &x, &y,
+		if ( integrate_peak(image, xp, yp, &x, &y,
                                     &intensity, NULL, NULL, 0, 0) ) continue;
 
 		itot += intensity;
diff --git a/src/geometry.h b/src/geometry.h
index 0782a3e1..a0761456 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -17,10 +17,12 @@
 #include <config.h>
 #endif
 
-extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                        int *n, int output);
+#include "reflist.h"
 
-extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
+extern RefList *find_intersections(struct image *image, UnitCell *cell,
+                                   int output);
+
+extern double integrate_all(struct image *image, RefList *reflections);
 
 
 #endif	/* GEOMETRY_H */
diff --git a/src/image.h b/src/image.h
index 10fb2406..ce0e7343 100644
--- a/src/image.h
+++ b/src/image.h
@@ -24,6 +24,7 @@
 #include "utils.h"
 #include "cell.h"
 #include "detector.h"
+#include "reflist.h"
 
 
 #define MAX_CELL_CANDIDATES (32)
@@ -52,35 +53,6 @@ struct imagefeature {
 typedef struct _imagefeaturelist ImageFeatureList;
 
 
-/* This structure represents a predicted peak in an image */
-struct cpeak
-{
-	/* Indices */
-	signed int h;
-	signed int k;
-	signed int l;
-
-	double min_distance;
-	int valid;
-
-	/* Partiality */
-	double r1;  /* First excitation error */
-	double r2;  /* Second excitation error */
-	double p;   /* Partiality */
-	int clamp1; /* Clamp status for r1 */
-	int clamp2; /* Clamp status for r2 */
-
-	/* Location in image */
-	int x;
-	int y;
-
-	int scalable;
-
-	/* Intensity */
-	double intensity;
-};
-
-
 /* Structure describing an image */
 struct image {
 
@@ -93,8 +65,7 @@ struct image {
 	struct detector         *det;
 	struct beam_params      *beam;  /* The nominal beam parameters */
 	char                    *filename;
-	struct cpeak            *cpeaks;
-	int                     n_cpeaks;
+	RefList                 *reflections;
 
 	int                     id;   /* ID number of the thread
 	                               * handling this image */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 6c01f95e..40d7cbd8 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -263,8 +263,7 @@ static struct image *get_simage(struct image *template, int alternate)
 	image->f0 = template->f0;
 
 	/* Prevent muppetry */
-	image->cpeaks = NULL;
-	image->n_cpeaks = 0;
+	image->reflections = NULL;
 
 	return image;
 }
@@ -326,8 +325,7 @@ static void process_image(void *pp, int cookie)
 	image.indexed_cell = NULL;
 	image.id = cookie;
 	image.filename = filename;
-	image.cpeaks = NULL;
-	image.n_cpeaks = 0;
+	image.reflections = NULL;
 	image.det = pargs->static_args.det;
 
 	STATUS("Processing '%s'\n", image.filename);
@@ -434,7 +432,7 @@ done:
 	free(image.data);
 	free(image.flags);
 	image_feature_list_free(image.features);
-	free(image.cpeaks);
+	reflist_free(image.reflections);
 	hdfile_close(hdfile);
 }
 
diff --git a/src/partialator.c b/src/partialator.c
index d1ad3768..a4f3381c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -34,6 +34,7 @@
 #include "beam-parameters.h"
 #include "post-refinement.h"
 #include "hrs-scaling.h"
+#include "reflist.h"
 
 
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
@@ -84,9 +85,9 @@ static void refine_image(int mytask, void *tasks)
 	struct image *image = pargs->image;
 	double nominal_photon_energy = pargs->image->beam->photon_energy;
 	struct hdfile *hdfile;
-	struct cpeak *spots;
-	int n, i;
+	int i;
 	double dev, last_dev;
+	RefList *reflections;
 
 	hdfile = hdfile_open(image->filename);
 	if ( hdfile == NULL ) {
@@ -110,16 +111,17 @@ static void refine_image(int mytask, void *tasks)
 	       a/1.0e-9, b/1.0e-9, c/1.0e-9,
 	       rad2deg(al), rad2deg(be), rad2deg(ga));
 
-	spots = find_intersections(image, image->indexed_cell, &n, 0);
+	/* FIXME: Don't do this each time */
+	reflections = find_intersections(image, image->indexed_cell, 0);
 	dev = +INFINITY;
 	i = 0;
 	do {
 		last_dev = dev;
-		dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n);
+		dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections);
 		STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
 		i++;
 	} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
-	mean_partial_dev(image, spots, n, pargs->sym,
+	mean_partial_dev(image, reflections, pargs->sym,
 	                 pargs->i_full, pargs->graph);
 	if ( pargs->pgraph ) {
 		fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
@@ -128,7 +130,7 @@ static void refine_image(int mytask, void *tasks)
 	free(image->data);
 	if ( image->flags != NULL ) free(image->flags);
 	hdfile_close(hdfile);
-	free(spots);
+	reflist_free(reflections);
 
 	/* Muppet proofing */
 	image->data = NULL;
@@ -163,22 +165,23 @@ static void refine_all(struct image *images, int n_total_patterns,
 }
 
 
-static void uniquify(struct cpeak *spot, const char *sym)
+static void uniquify(Reflection *refl, const char *sym)
 {
+	signed int h, k, l;
 	signed int ha, ka, la;
 
-	get_asymm(spot->h, spot->k, spot->l, &ha, &ka, &la, sym);
-	spot->h = ha;
-	spot->k = ka;
-	spot->l = la;
+	get_indices(refl, &h, &k, &l);
+	get_asymm(h, k, l, &ha, &ka, &la, sym);
+	set_indices(refl, h, k, l);
 }
 
 
+/* FIXME: Get rid of this */
 static void integrate_image(struct image *image, ReflItemList *obs,
                             const char *sym)
 {
-	struct cpeak *spots;
-	int j, n;
+	RefList *reflections;
+	Reflection *refl;
 	struct hdfile *hdfile;
 	double nominal_photon_energy = image->beam->photon_energy;
 
@@ -199,42 +202,39 @@ static void integrate_image(struct image *image, ReflItemList *obs,
 	}
 
 	/* Figure out which spots should appear in this pattern */
-	spots = find_intersections(image, image->indexed_cell, &n, 0);
+	reflections = find_intersections(image, image->indexed_cell, 0);
 
 	/* For each reflection, estimate the partiality */
-	for ( j=0; j<n; j++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		signed int h, k, l;
 		float i_partial;
 		float xc, yc;
+		double x, y;
 
-		uniquify(&spots[j], sym);
-
-		h = spots[j].h;
-		k = spots[j].k;
-		l = spots[j].l;
+		uniquify(refl, sym);
+		get_indices(refl, &h, &k, &l);
 
 		/* Don't attempt to use spots with very small
 		 * partialities, since it won't be accurate. */
-		if ( spots[j].p < 0.1 ) continue;
+		if ( get_partiality(refl) < 0.1 ) continue;
 
-		/* Actual measurement of this reflection from this
-		 * pattern? */
-		/* FIXME: Coordinates aren't whole numbers */
-		if ( integrate_peak(image, spots[j].x, spots[j].y,
+		/* Actual measurement of this reflection from this pattern? */
+		get_detector_pos(refl, &x, &y);
+		if ( integrate_peak(image, x, y,
 		                    &xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
-			spots[j].valid = 0;
+			delete_refl(refl);
 			continue;
 		}
 
-		spots[j].intensity = i_partial;
-		spots[j].valid = 1;
+		set_int(refl, i_partial);
 
 		if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
 
 	}
-	image->cpeaks = spots;
-	image->n_cpeaks = n;
+	image->reflections = reflections;
 
 	free(image->data);
 	if ( image->flags != NULL ) free(image->flags);
@@ -253,21 +253,20 @@ static void select_scalable_reflections(struct image *images, int n)
 
 	for ( m=0; m<n; m++ ) {
 
-		int j;
+		Reflection *refl;
 
-		for ( j=0; j<images[m].n_cpeaks; j++ ) {
+		for ( refl = first_refl(images[m].reflections);
+		      refl != NULL;
+		      refl = next_refl(refl) ) {
 
 			int scalable = 1;
+			double v;
 
-			if ( images[m].cpeaks[j].p < 0.1 ) scalable = 0;
-			if ( !images[m].cpeaks[j].valid ) {
-				scalable = 0;
-			} else {
-				double v = fabs(images[m].cpeaks[j].intensity);
-				if ( v < 0.1 ) scalable = 0;
-			}
+			if ( get_partiality(refl) < 0.1 ) scalable = 0;
+			v = fabs(get_intensity(refl));
+			if ( v < 0.1 ) scalable = 0;
 
-			images[m].cpeaks[j].scalable = scalable;
+			set_scalable(refl, scalable);
 
 		}
 
@@ -525,7 +524,7 @@ int main(int argc, char *argv[])
 
 	/* Clean up */
 	for ( i=0; i<n_total_patterns; i++ ) {
-		free(images[i].cpeaks);
+		reflist_free(images[i].reflections);
 	}
 	free(I_full);
 	delete_items(obs);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 4584bb3c..90425328 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -567,7 +567,7 @@ int main(int argc, char *argv[])
 			find_projected_peaks(&image, cell, 0, 0.1);
 			output_intensities(&image, cell, NULL, 0, 0, stdout,
 			                   0, 0.1);
-			free(image.cpeaks);
+			reflist_free(image.reflections);
 		}
 
 		if ( powder_fn != NULL ) {
diff --git a/src/peaks.c b/src/peaks.c
index 88ac50bf..d1132a3a 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -32,9 +32,6 @@
 #include "diffraction.h"
 
 
-/* Maximum number of peaks which may be predicted by find_projected_peaks() */
-#define MAX_CPEAKS (8192)
-
 /* How close a peak must be to an indexed position to be considered "close"
  * for the purposes of double hit detection and sanity checking. */
 #define PEAK_CLOSE (30.0)
@@ -483,12 +480,11 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
-	struct cpeak *cpeaks;
-	int n_cpeaks = 0;
+	RefList *reflections;
 	double alen, blen, clen;
+	int n_reflections = 0;
 
-	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
-	if ( cpeaks == NULL ) return 0;
+	reflections = reflist_new();
 
 	/* "Borrow" direction values to get reciprocal lengths */
 	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -507,8 +503,8 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
 		signed int h, k, l;
 		struct rvec q;
 		double dist;
-		int found = 0;
-		int j;
+		Reflection *refl;
+		double cur_dist;
 
 		q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
 
@@ -535,44 +531,26 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
 			if ( dist > domain_r ) continue;
 		}
 
-		for ( j=0; j<n_cpeaks; j++ ) {
-			if ( (cpeaks[j].h == h) && (cpeaks[j].k == k)
-			                      && (cpeaks[j].l == l) ) {
-
-				if ( dist < cpeaks[j].min_distance ) {
-					cpeaks[j].min_distance = dist;
-					cpeaks[j].x = x;
-					cpeaks[j].y = y;
-				}
-
-				found = 1;
-
-			}
-		}
-
-		if ( !found ) {
-			cpeaks[n_cpeaks].min_distance = dist;
-			cpeaks[n_cpeaks].x = x;
-			cpeaks[n_cpeaks].y = y;
-			cpeaks[n_cpeaks].h = h;
-			cpeaks[n_cpeaks].k = k;
-			cpeaks[n_cpeaks].l = l;
-			n_cpeaks++;
-			if ( n_cpeaks == MAX_CPEAKS ) {
-				ERROR("Unit cell is insanely large!\n");
-				goto out;
+		refl = find_refl(reflections, h, k, l);
+		if ( refl != NULL ) {
+			cur_dist = get_excitation_error(refl);
+			if ( dist < cur_dist ) {
+				set_detector_pos(refl, dist, x, y);
 			}
+		} else {
+			add_refl_with_det_pos(reflections, h, k, l, x, y, dist);
+			n_reflections++;
 		}
 
 	}
 	}
 
-out:
-	STATUS("Found %i reflections\n", n_cpeaks);
-	image->cpeaks = cpeaks;
-	image->n_cpeaks = n_cpeaks;
+	optimise_reflist(reflections);
+
+	STATUS("Found %i reflections\n", n_reflections);
+	image->reflections = reflections;
 
-	return n_cpeaks;
+	return n_reflections;
 }
 
 
@@ -685,20 +663,14 @@ void output_intensities(struct image *image, UnitCell *cell,
                         int use_closer, FILE *ofh,
                         int circular_domain, double domain_r)
 {
-	int i;
-	int n_found;
-	int n_indclose = 0;
-	int n_foundclose = 0;
-	int n_veto = 0;
-	int n_veto_second = 0;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
+	Reflection *refl;
 
-	if ( image->n_cpeaks == 0 ) {
+	if ( image->reflections == NULL ) {
 		find_projected_peaks(image, cell, circular_domain, domain_r);
 	}
-	if ( image->n_cpeaks == 0 ) return;
 
 	/* Get exclusive access to the output stream if necessary */
 	if ( mutex != NULL ) pthread_mutex_lock(mutex);
@@ -709,16 +681,20 @@ void output_intensities(struct image *image, UnitCell *cell,
 	                          &bsx, &bsy, &bsz,
 	                          &csx, &csy, &csz);
 
-	for ( i=0; i<image->n_cpeaks; i++ ) {
+	for ( refl = first_refl(image->reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		float x, y, intensity;
 		double d;
 		int idx;
 		double bg, max;
 		struct panel *p;
+		double px, py;
+		signed int h, k, l;
 
-		p = find_panel(image->det, image->cpeaks[i].x,
-		                           image->cpeaks[i].y);
+		get_detector_pos(refl, &px, &py);
+		p = find_panel(image->det, px, py);
 		if ( p == NULL ) continue;
 		if ( p->no_index ) continue;
 
@@ -729,9 +705,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 			if ( image->features != NULL ) {
 				f = image_feature_closest(image->features,
-					                  image->cpeaks[i].x,
-					                  image->cpeaks[i].y,
-					                  &d, &idx);
+					                  px, py, &d, &idx);
 			} else {
 				f = NULL;
 			}
@@ -754,7 +728,6 @@ void output_intensities(struct image *image, UnitCell *cell,
 					 * region.  Integration of the same
 					 * peak included a bad region this time.
 					 */
-					n_veto_second++;
 					continue;
 				}
 				intensity = f->intensity;
@@ -763,81 +736,36 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 				int r;
 
-				r = integrate_peak(image,
-					           image->cpeaks[i].x,
-					           image->cpeaks[i].y,
-					           &x, &y, &intensity, &bg, &max,
+				r = integrate_peak(image, px, py, &x, &y,
+				                   &intensity, &bg, &max,
 					           polar, 1);
 				if ( r ) {
 					/* Plain old ordinary peak veto */
-					n_veto++;
 					continue;
 				}
 
 			}
 
-			if ( (f != NULL) && (d < PEAK_CLOSE) ) {
-				n_indclose++;
-			}
-
 		} else {
 
 			int r;
 
-			r = integrate_peak(image,
-				           image->cpeaks[i].x,
-				           image->cpeaks[i].y,
-				           &x, &y, &intensity, &bg, &max, polar,
-				           0);
+			r = integrate_peak(image, px, py, &x, &y,
+			                   &intensity, &bg, &max, polar, 0);
 			if ( r ) {
 				/* Plain old ordinary peak veto */
-				n_veto++;
 				continue;
 			}
 
 		}
 
 		/* Write h,k,l, integrated intensity and centroid coordinates */
+		get_indices(refl, &h, &k, &l);
 		fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n",
-		       image->cpeaks[i].h, image->cpeaks[i].k, image->cpeaks[i].l,
-		       intensity, x, y, max, bg);
-
-		image->cpeaks[i].x = x;
-		image->cpeaks[i].y = y;
-
-	}
-
-	n_found = image_feature_count(image->features);
-	for ( i=0; i<n_found; i++ ) {
-
-		struct imagefeature *f;
-		int j;
-
-		f = image_get_feature(image->features, i);
-
-		if ( f == NULL ) continue;
-
-		for ( j=0; j<image->n_cpeaks; j++ ) {
-
-			double d;
-
-			d = pow(image->cpeaks[j].x-f->x, 2.0)
-			  + pow(image->cpeaks[j].y-f->y, 2.0);
-
-			if ( d < PEAK_CLOSE ) n_foundclose++;
-
-		}
+		        h, l, l, intensity, x, y, max, bg);
 
 	}
 
-	fprintf(ofh, "Peak statistics: %i peaks found by the peak search out of "
-	       "%i were close to indexed positions. "
-	       "%i indexed positions out of %i were close to detected peaks.\n",
-	       n_foundclose, n_found, n_indclose, image->n_cpeaks);
-	fprintf(ofh, "%i integrations using indexed locations were aborted because "
-	       "they hit one or more bad pixels.\n", n_veto);
-	fprintf(ofh, "%i integrations using peak search locations were aborted "
-	       "because they hit one or more bad pixels.\n", n_veto_second);
 	/* Blank line at end */
 	fprintf(ofh, "\n");
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index b860bb42..25664e69 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -65,7 +65,7 @@ static double partiality_rgradient(double r, double profile_radius)
 
 
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, struct cpeak spot, double r)
+double gradient(struct image *image, int k, Reflection *refl, double r)
 {
 	double ds, tt, azi;
 	double nom, den;
@@ -74,24 +74,32 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	double xl, yl, zl;
+	signed int hi, ki, li;
+	double r1, r2, p;
+	int clamp_low, clamp_high;
+
+	get_indices(refl, &hi, &ki, &li);
 
 	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
 	                                         &bsx, &bsy, &bsz,
 	                                         &csx, &csy, &csz);
-	xl = spot.h*asx + spot.k*bsx + spot.l*csx;
-	yl = spot.h*asy + spot.k*bsy + spot.l*csy;
-	zl = spot.h*asz + spot.k*bsz + spot.l*csz;
+	xl = hi*asx + ki*bsx + li*csx;
+	yl = hi*asy + ki*bsy + li*csy;
+	zl = hi*asz + ki*bsz + li*csz;
+
 
-	ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
+	ds = 2.0 * resolution(image->indexed_cell, hi, ki, li);
 	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+k);
 	azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
 
+	get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
+
 	/* Calculate the gradient of partiality wrt excitation error. */
-	if ( spot.clamp1 == 0 ) {
-		g += partiality_gradient(spot.r1, r);
+	if ( clamp_low == 0 ) {
+		g += partiality_gradient(r1, r);
 	}
-	if ( spot.clamp2 == 0 ) {
-		g += partiality_gradient(spot.r2, r);
+	if ( clamp_high == 0 ) {
+		g += partiality_gradient(r2, r);
 	}
 
 	/* For many gradients, just multiply the above number by the gradient
@@ -99,7 +107,7 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
 	switch ( k ) {
 
 	case REF_SCALE :
-		return -spot.p*pow(image->osf, -2.0);
+		return -p*pow(image->osf, -2.0);
 
 	case REF_DIV :
 		nom = sqrt(2.0) * ds * sin(image->div/2.0);
@@ -107,33 +115,33 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
 		return (nom/den) * g;
 
 	case REF_R :
-		if ( spot.clamp1 == 0 ) {
-			g += partiality_rgradient(spot.r1, r);
+		if ( clamp_low == 0 ) {
+			g += partiality_rgradient(r1, r);
 		}
-		if ( spot.clamp2 == 0 ) {
-			g += partiality_rgradient(spot.r2, r);
+		if ( clamp_high == 0 ) {
+			g += partiality_rgradient(r2, r);
 		}
 		return g;
 
 	/* Cell parameters and orientation */
 	case REF_ASX :
-		return spot.h * pow(sin(tt), -1.0) * cos(azi) * g;
+		return hi * pow(sin(tt), -1.0) * cos(azi) * g;
 	case REF_BSX :
-		return spot.k * pow(sin(tt), -1.0) * cos(azi) * g;
+		return ki * pow(sin(tt), -1.0) * cos(azi) * g;
 	case REF_CSX :
-		return spot.l * pow(sin(tt), -1.0) * cos(azi) * g;
+		return li * pow(sin(tt), -1.0) * cos(azi) * g;
 	case REF_ASY :
-		return spot.h * pow(sin(tt), -1.0) * sin(azi) * g;
+		return hi * pow(sin(tt), -1.0) * sin(azi) * g;
 	case REF_BSY :
-		return spot.k * pow(sin(tt), -1.0) * sin(azi) * g;
+		return ki * pow(sin(tt), -1.0) * sin(azi) * g;
 	case REF_CSY :
-		return spot.l * pow(sin(tt), -1.0) * sin(azi) * g;
+		return li * pow(sin(tt), -1.0) * sin(azi) * g;
 	case REF_ASZ :
-		return spot.h * pow(cos(tt), -1.0) * g;
+		return hi * pow(cos(tt), -1.0) * g;
 	case REF_BSZ :
-		return spot.k * pow(cos(tt), -1.0) * g;
+		return ki * pow(cos(tt), -1.0) * g;
 	case REF_CSZ :
-		return spot.l * pow(cos(tt), -1.0) * g;
+		return li * pow(cos(tt), -1.0) * g;
 
 	}
 
@@ -218,31 +226,34 @@ void apply_shift(struct image *image, int k, double shift)
 }
 
 
-double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
+double mean_partial_dev(struct image *image, RefList *reflections,
                         const char *sym, double *i_full, FILE *graph)
 {
-	int h, n_used;
+	int n_used;
 	double delta_I = 0.0;
+	Reflection *refl;
 
 	n_used = 0;
-	for ( h=0; h<n; h++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		signed int hind, kind, lind;
 		signed int ha, ka, la;
 		double I_full;
 		float I_partial;
 		float xc, yc;
+		double x, y;
+		double p;
 
-		hind = spots[h].h;
-		kind = spots[h].k;
-		lind = spots[h].l;
+		get_indices(refl, &hind, &kind, &lind);
 
-		if ( !spots[h].scalable ) continue;
+		if ( !get_scalable(refl) ) continue;
 
-		/* Actual measurement of this reflection from this
-		 * pattern? */
-		/* FIXME: Coordinates aren't whole numbers */
-		if ( integrate_peak(image, spots[h].x, spots[h].y,
+		/* Actual measurement of this reflection from this pattern? */
+		get_detector_pos(refl, &x, &y);
+		/* FIXME: Get rid of this */
+		if ( integrate_peak(image, x, y,
 		                    &xc, &yc, &I_partial, NULL, NULL,
 		                    1, 0) ) {
 			continue;
@@ -250,14 +261,15 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
 
 		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
 		I_full = lookup_intensity(i_full, ha, ka, la);
-		delta_I += fabs(I_partial - (spots[h].p * I_full / image->osf));
+		p = get_partiality(refl);
+		delta_I += fabs(I_partial - (p * I_full / image->osf));
 		n_used++;
 
 		if ( graph != NULL ) {
 			fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
 			               " %5.2f %5.2f\n",
-			       hind, kind, lind, I_partial/spots[h].p, xc, yc,
-			       spots[h].p, I_partial / I_full);
+			       hind, kind, lind, I_partial/p, xc, yc,
+			       p, I_partial / I_full);
 		}
 
 	}
@@ -268,19 +280,21 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
 
 /* Perform one cycle of post refinement on 'image' against 'i_full' */
 double pr_iterate(struct image *image, double *i_full, const char *sym,
-                  struct cpeak **pspots, int *n)
+                  RefList *reflections)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
 	gsl_vector *shifts;
-	int h, param;
-	struct cpeak *spots = *pspots;
+	int param;
+	Reflection *refl;
 
 	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
 	v = gsl_vector_calloc(NUM_PARAMS);
 
 	/* Construct the equations, one per reflection in this image */
-	for ( h=0; h<*n; h++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		signed int hind, kind, lind;
 		signed int ha, ka, la;
@@ -288,17 +302,17 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 		float I_partial;
 		float xc, yc;
 		int k;
+		double x, y;
+		double p;
 
-		hind = spots[h].h;
-		kind = spots[h].k;
-		lind = spots[h].l;
+		get_indices(refl, &hind, &kind, &lind);
 
-		if ( !spots[h].scalable ) continue;
+		if ( !get_scalable(refl) ) continue;
 
-		/* Actual measurement of this reflection from this
-		 * pattern? */
-		/* FIXME: Coordinates aren't whole numbers */
-		if ( integrate_peak(image, spots[h].x, spots[h].y,
+		/* Actual measurement of this reflection from this pattern? */
+		get_detector_pos(refl, &x, &y);
+		/* FIXME: Get rid of this */
+		if ( integrate_peak(image, x, y,
 		                    &xc, &yc, &I_partial, NULL, NULL,
 		                    1, 0) ) {
 			continue;
@@ -306,7 +320,8 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 
 		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
 		I_full = lookup_intensity(i_full, ha, ka, la);
-		delta_I = I_partial - (spots[h].p * I_full / image->osf);
+		p = get_partiality(refl);
+		delta_I = I_partial - (p * I_full / image->osf);
 
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 
@@ -319,9 +334,9 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 
 				M_curr = gsl_matrix_get(M, g, k);
 
-				M_c = gradient(image, g, spots[h],
+				M_c = gradient(image, g, refl,
 				               image->profile_radius)
-				    * gradient(image, k, spots[h],
+				    * gradient(image, k, refl,
 				               image->profile_radius);
 				M_c *= pow(I_full, 2.0);
 
@@ -329,8 +344,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 
 			}
 
-			gr = gradient(image, k, spots[h],
-			              image->profile_radius);
+			gr = gradient(image, k, refl, image->profile_radius);
 			v_c = delta_I * I_full * gr;
 			gsl_vector_set(v, k, v_c);
 
@@ -350,8 +364,5 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 	gsl_vector_free(v);
 	gsl_vector_free(shifts);
 
-	free(spots);
-	spots = find_intersections(image, image->indexed_cell, n, 0);
-	*pspots = spots;
-	return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
+	return mean_partial_dev(image, reflections, sym, i_full, NULL);
 }
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 983ecc28..b3c6c074 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -21,6 +21,7 @@
 #include <stdio.h>
 
 #include "image.h"
+#include "reflist.h"
 
 
 /* Refineable parameters */
@@ -44,12 +45,12 @@ enum {
 void apply_shift(struct image *image, int k, double shift);
 
 
-double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
+double mean_partial_dev(struct image *image, RefList *reflections,
                         const char *sym, double *i_full, FILE *graph);
 
 
 double pr_iterate(struct image *image, double *i_full, const char *sym,
-                  struct cpeak **pspots, int *n);
+                  RefList *reflections);
 
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/reflist.c b/src/reflist.c
new file mode 100644
index 00000000..990bc535
--- /dev/null
+++ b/src/reflist.c
@@ -0,0 +1,139 @@
+/*
+ * reflist.c
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+
+#include "reflist.h"
+
+
+struct _reflection {
+
+	unsigned int serial;  /* Serial number */
+
+	signed int h;
+	signed int k;
+	signed int l;
+
+	double excitation_error;
+
+	/* Partiality */
+	double r1;  /* First excitation error */
+	double r2;  /* Second excitation error */
+	double p;   /* Partiality */
+	int clamp1; /* Clamp status for r1 */
+	int clamp2; /* Clamp status for r2 */
+
+	/* Location in image */
+	int x;
+	int y;
+
+	int scalable;
+
+	/* Intensity */
+	double intensity;
+
+
+};
+
+
+struct _reflist {
+
+	struct ref *head;
+
+};
+
+
+/**************************** Creation / deletion *****************************/
+
+/* Create a reflection list */
+RefList *reflist_new()
+{
+	RefList *new;
+
+	new = malloc(sizeof(struct _reflist));
+
+	return new;
+}
+
+
+void reflist_free(RefList *list)
+{
+}
+
+
+/********************************** Search ************************************/
+
+Reflection *find_refl(RefList *list, INDICES)
+{
+}
+
+
+/********************************** Getters ***********************************/
+
+double get_excitation_error(Reflection *refl)
+{
+}
+
+
+void get_detector_pos(Reflection *refl, double *x, double *y)
+{
+}
+
+
+void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
+{
+}
+
+
+/********************************** Setters ***********************************/
+
+void set_detector_pos(Reflection *refl, double exerr, double x, double y)
+{
+}
+
+
+void set_partial(Reflection *refl, double r1, double r2, double p,
+                 double clamp_low, double clamp_high)
+{
+}
+
+
+/********************************* Insertion **********************************/
+
+Reflection *add_refl(RefList *list, INDICES)
+{
+}
+
+
+Reflection *add_refl_with_det_pos(RefList *refl, INDICES, double exerr,
+                                  double x, double y)
+{
+}
+
+
+/********************************* Iteration **********************************/
+
+Reflection *first_refl(RefList *list)
+{
+}
+
+
+Reflection *next_refl(Reflection *refl)
+{
+}
+
+
+/*********************************** Voodoo ***********************************/
+
+void optimise_reflist(RefList *list)
+{
+}
diff --git a/src/reflist.h b/src/reflist.h
new file mode 100644
index 00000000..0ba09a3b
--- /dev/null
+++ b/src/reflist.h
@@ -0,0 +1,69 @@
+/*
+ * reflist.h
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef REFLIST_H
+#define REFLIST_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+typedef struct _reflist RefList;
+typedef struct _reflection Reflection;
+
+#define INDICES signed int h, signed int k, signed int l
+
+/* Creation/deletion */
+extern RefList *reflist_new(void);
+extern void reflist_free(RefList *list);
+
+/* Search */
+extern Reflection *find_refl(RefList *list, INDICES);
+
+/* Get */
+extern double get_excitation_error(Reflection *refl);
+extern void get_detector_pos(Reflection *refl, double *x, double *y);
+extern void get_indices(Reflection *refl,
+                        signed int *h, signed int *k, signed int *l);
+extern double get_partiality(Reflection *refl);
+extern double get_intensity(Reflection *refl);
+extern void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+                        int *clamp_low, int *clamp_high);
+extern int get_scalable(Reflection *refl);
+
+/* Set */
+extern void set_detector_pos(Reflection *refl, double exerr,
+                             double x, double y);
+extern void set_partial(Reflection *refl, double r1, double r2, double p,
+                        double clamp_low, double clamp_high);
+extern void set_indices(Reflection *refl,
+                        signed int h, signed int k, signed int l);
+extern void set_int(Reflection *refl, double intensity);
+extern void set_scalable(Reflection *refl, int scalable);
+
+/* Insertion */
+extern Reflection *add_refl(RefList *list, INDICES);
+extern Reflection *add_refl_with_det_pos(RefList *list, INDICES, double exerr,
+                                        double x, double y);
+
+
+/* Deletion */
+extern void delete_refl(Reflection *refl);
+
+/* Iteration */
+extern Reflection *first_refl(RefList *list);
+extern Reflection *next_refl(Reflection *refl);
+
+/* Voodoo */
+extern void optimise_reflist(RefList *list);
+
+#endif	/* REFLIST_H */
diff --git a/src/reintegrate.c b/src/reintegrate.c
index 9d882d4b..09634cd7 100644
--- a/src/reintegrate.c
+++ b/src/reintegrate.c
@@ -110,8 +110,7 @@ static void process_image(void *pg, int cookie)
 	image.flags = NULL;
 	image.indexed_cell = NULL;
 	image.filename = apargs->filename;
-	image.cpeaks = NULL;
-	image.n_cpeaks = 0;
+	image.reflections = NULL;
 	image.det = pargs->det;
 
 	STATUS("Processing '%s'\n", apargs->filename);
diff --git a/src/templates.c b/src/templates.c
index 4bd85a76..9abbcb17 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -21,6 +21,7 @@
 #include "geometry.h"
 #include "hdf5-file.h"
 #include "peaks.h"
+#include "reflist.h"
 
 #include <assert.h>
 
@@ -42,8 +43,7 @@ struct _indexingprivate_template
 struct template {
 	double omega;
 	double phi;
-	int n;
-	struct cpeak *spots;
+	RefList *spots;
 };
 
 
@@ -164,24 +164,22 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
 	for ( omega = 0.0; omega < omega_max-omega_step; omega += omega_step ) {
 	for ( phi = 0.0; phi < phi_max-phi_step; phi += phi_step ) {
 
-		int n;
-		struct cpeak *cpeaks;
+		RefList *reflections;
 		UnitCell *cell_rot;
 
 		assert(i < n_templates);
 
 		cell_rot = rotate_cell(cell, omega, phi, 0.0);
 
-		cpeaks = find_intersections(&image, cell_rot, &n, 0);
-		if ( cpeaks == NULL ) {
+		reflections = find_intersections(&image, cell_rot, 0);
+		if ( reflections == NULL ) {
 			ERROR("Template calculation failed.\n");
 			return NULL;
 		}
 
 		priv->templates[i].omega = omega;
 		priv->templates[i].phi = phi;
-		priv->templates[i].n = n;
-		priv->templates[i].spots = cpeaks;
+		priv->templates[i].spots = reflections;
 		i++;
 
 		free(cell_rot);
@@ -227,23 +225,28 @@ static int fast_integrate_peak(struct image *image, int xp, int yp)
 }
 
 
-static double integrate_all_rot(struct image *image, struct cpeak *cpeaks,
-                                int n, double rot)
+static double integrate_all_rot(struct image *image, RefList *reflections,
+                                double rot)
 {
 	double itot = 0.0;
-	int i;
 	double cosr, sinr;
 	int num_int = 0;
+	Reflection *refl;
 
 	cosr = cos(rot);
 	sinr = sin(rot);
 
-	for ( i=0; i<n; i++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		float xp, yp;
+		double x, y;
+
+		get_detector_pos(refl, &x, &y);
 
-		xp =  cosr*cpeaks[i].x + sinr*cpeaks[i].y;
-		yp = -sinr*cpeaks[i].x + cosr*cpeaks[i].y;
+		xp =  cosr*x + sinr*y;
+		yp = -sinr*x + cosr*y;
 
 		itot += fast_integrate_peak(image, rint(xp), rint(yp));
 		num_int++;
@@ -256,26 +259,31 @@ static double integrate_all_rot(struct image *image, struct cpeak *cpeaks,
 
 /* Return the mean of the distances between peaks in the image and peaks from
  * the given template. */
-static double mean_distance(struct image *image, struct cpeak *cpeaks,
-                            int n, double rot)
+static double mean_distance(struct image *image, RefList *reflections,
+                            double rot)
 {
 	double dtot = 0.0;
-	int i;
 	double cosr, sinr;
 	int num_dist = 0;
+	Reflection *refl;
 
 	cosr = cos(rot);
 	sinr = sin(rot);
 
 	/* For each template peak */
-	for ( i=0; i<n; i++ ) {
+	for ( refl = first_refl(reflections);
+	      refl != NULL;
+	      refl = next_refl(refl) ) {
 
 		float xp, yp;
 		int j;
 		double min_dsq;
+		double x, y;
+
+		get_detector_pos(refl, &x, &y);
 
-		xp =  cosr*cpeaks[i].x + sinr*cpeaks[i].y;
-		yp = -sinr*cpeaks[i].x + cosr*cpeaks[i].y;
+		xp =  cosr*x + sinr*y;
+		yp = -sinr*x + cosr*y;
 
 		/* Compare to every real peak */
 		min_dsq = +INFINITY;
@@ -328,11 +336,11 @@ void match_templates(struct image *image, IndexingPrivate *ipriv)
 
 		if ( !peaks ) {
 			val = integrate_all_rot(image, priv->templates[i].spots,
-			                        priv->templates[i].n, rot);
+			                        rot);
 			best = val > max;
 		} else {
 			val = mean_distance(image, priv->templates[i].spots,
-			                    priv->templates[i].n, rot);
+			                    rot);
 			best = val < max;
 		}