Refactor image handling code

Remove itrans-lsq


git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@158 bf6ca9ba-c028-0410-8290-897cf20841d1

30 files changed, 425 insertions, 1004 deletions

diff --git a/Makefile.am b/Makefile.am
index 1e84891..205c21a 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,5 +1,5 @@
 EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src/control.h src/readpng.h src/mrc.h src/imagedisplay.h src/main.h \
-		data/displaywindow.ui src/utils.h src/itrans.h src/qdrp.h src/cache.h src/itrans-threshold.h \
-		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h src/basis.h \
-		src/dirax.h
+		data/displaywindow.ui src/utils.h src/itrans.h src/qdrp.h src/cache.h src/itrans-threshold.h src/basis.h \
+		src/itrans-zaefferer.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h \
+		src/dirax.h src/image.h
 SUBDIRS = src data
diff --git a/src/Makefile.am b/src/Makefile.am
index eef1fa6..0a3fbb2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,7 +1,7 @@
 bin_PROGRAMS = dtr
 dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c cache.c \
-		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \
-		dirax.c
+		itrans-threshold.c itrans-zaefferer.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \
+		dirax.c image.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
diff --git a/src/basis.c b/src/basis.c
index 8a282e3..841ccc6 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -1,7 +1,7 @@
 /*
  * basis.c
  *
- * Find approximate lattices to feed various procedures
+ * Handle basis structures
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
  *
@@ -14,10 +14,7 @@
 #endif
 
 #include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
 
-#include "utils.h"
 #include "reflections.h"
 #include "basis.h"
 
@@ -71,232 +68,3 @@ static double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
 
 }
 
-static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) {
-	
-	Reflection	*tcentre;
-	int		lfom;
-	double		tol;
-	int		j;
-	
-	lfom = 0;
-	tol = modulus(vx, vy, vz)/10.0;
-	tcentre = ctx->reflectionlist->reflections;
-	do {
-	
-		for ( j=-20; j<=20; j++ ) {
-		
-			Reflection *check;
-			
-			check = reflectionlist_find_nearest(ctx->reflectionlist, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j);
-			if ( check && (distance3d(check->x, check->y, check->z, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j) < tol) ) {
-				lfom++;
-			}
-			
-		}
-		
-		tcentre = tcentre->next;
-		
-	} while ( tcentre );
-
-	return lfom;
-
-}
-
-/* Select a suitable number of sensible seed vectors */
-static ReflectionList *basis_find_seeds(ControlContext *ctx) {
-
-	double		tilt_min;
-	double		tilt_max;
-	double		tilt_mid;
-	ImageRecord	*imagerecord;
-	double		x_temp, y_temp, z_temp;
-	double		scale;
-	double		x, y, z;
-	Reflection	*centre;
-	int		i;
-	ReflectionList *seeds;
-	
-	seeds = reflectionlist_new();
-	
-	/* Locate the 'plane' in the middle of the "wedge".
-	 *	This whole procedure assumes there is just one tilt axis. */
-	tilt_min = control_min_tilt(ctx);
-	tilt_max = control_max_tilt(ctx);
-	tilt_mid = tilt_min + (tilt_max-tilt_min)/2;
-	imagerecord = control_image_nearest_tilt(ctx, tilt_mid);
-	
-	/* Apply the last two steps of the mapping transform to get the direction from the origin
-	 *	towards the middle of the wedge */
-	x_temp = 0.0;
-	y_temp = cos(deg2rad(imagerecord->tilt));
-	z_temp = -sin(deg2rad(imagerecord->tilt));
-	x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega));
-	y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega));
-	z = z_temp;
-	
-	/* Find the point in the middle of the "wedge" */
-	scale = reflectionlist_largest_g(ctx->reflectionlist)/6;
-	x *= scale;
-	y *= scale;
-	z *= scale;
-	reflection_add(ctx->reflectionlist, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
-	
-	/* Find an "origin" reflection */
-	centre = reflectionlist_find_nearest(ctx->reflectionlist, x, y, z);
-	if ( !centre ) return NULL;
-	centre->found = 1;
-	reflection_add(ctx->reflectionlist, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED);
-	
-	for ( i=1; i<=10; i++ ) {
-	
-		Reflection *vector;
-		Reflection *new_seed;
-		int accept, lfom;
-		double vx, vy, vz;
-		
-		do {
-			
-			Reflection	*check;
-			
-			accept = 1;
-		
-			/* Find a "candidate vector" reflection */
-			vector = reflectionlist_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9);
-			if ( !vector ) {
-				printf("BS: Only found %i seeds\n", i);
-				return seeds;
-			}
-			vector->found = 1;
-			
-			/* Get vector components (not the coordinates the vector was calculated from!) */
-			vx = vector->x - centre->x;
-			vy = vector->y - centre->y;
-			vz = vector->z - centre->z;
-			
-			/* Proximity test: don't duplicate seeds */
-			check = reflectionlist_find_nearest(seeds, vx, vy, vz);
-			if ( check ) {
-				if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) {
-					/* Too close to another seed */
-					accept = 0;
-					continue;
-				}
-			}
-			
-			/* Record lFOM for later analysis */
-			lfom = basis_lfom(ctx, vx, vy, vz);
-			
-		} while ( !accept );
-
-		/* Add to the list of seeds */
-		new_seed = reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL);
-		new_seed->lfom = lfom;
-		
-		/* Create a marker in the default list for visualisation */
-		reflection_add(ctx->reflectionlist, vx, vy, vz, 1.0, REFLECTION_MARKER);
-
-	}
-
-	return seeds;
-
-}
-
-/* Assemble the most sensible basis from seeds */
-static Basis *basis_assemble_seeds(ReflectionList *seeds) {
-
-	Basis *basis;
-	Reflection *ref;
-	int i, j, b;
-	ReflectionList *seeds_sorted;
-
-	seeds_sorted = reflectionlist_sort_lfom(seeds);
-	
-	basis = malloc(10*sizeof(Basis));
-	
-	for ( b=0; b<10; b++ ) {
-	
-		ref = seeds_sorted->reflections;
-		j = 0;	/* Number of basis components already found */
-		for ( i=1; i<=seeds_sorted->n_reflections; i++ ) {
-		
-			double vx, vy, vz;
-		
-			vx = ref->x;
-			vy = ref->y;
-			vz = ref->z;
-			
-			printf("Seed %2i: lFOM=%6i", i, ref->lfom);
-			
-			switch ( j ) {
-				case 0 : {
-					basis[b].a.x = vx;
-					basis[b].a.y = vy;
-					basis[b].a.z = vz;
-					j++;
-					printf(" *");
-					break;
-				}
-				case 1 : {
-					if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6)
-					  && (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) < M_PI-M_PI/6) ) {
-						basis[b].b.x = vx;
-						basis[b].b.y = vy;
-						basis[b].b.z = vz;
-						j++;
-						printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z)));
-						break;
-					}
-				}
-				case 2 : {
-					double cx, cy, cz;
-					cx = basis[b].a.y*basis[b].b.z - basis[b].a.z*basis[b].b.y;
-					cy =  - basis[b].a.x*basis[b].b.z + basis[b].a.z*basis[b].b.x;
-					cz = basis[b].a.x*basis[b].b.y - basis[b].a.y*basis[b].b.x;
-					if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6)
-					  && (angle_between(vx, vy, vz, basis[b].b.x, basis[b].b.y, basis[b].b.z) > M_PI/6)
-					  && (angle_between(vx, vy, vz, cx, cy, cz) < M_PI/2-M_PI/6) ) {
-						basis[b].c.x = vx;
-						basis[b].c.y = vy;
-						basis[b].c.z = vz;
-						j++;
-						printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, cx, cy, cz)));
-						break;
-					}
-				}
-			}
-			
-			printf("\n");
-			if ( j >= 3 ) break;
-			ref = ref->next;
-			
-		}
-		if ( j < 3 ) {
-			break;
-		}
-	}
-	
-	return basis;
-
-}
-
-Basis *basis_find(ControlContext *ctx) {
-
-	Basis *basis;
-	ReflectionList *seeds;
-	
-	/* Get the shortlist of seeds */
-	seeds = basis_find_seeds(ctx);
-	if ( seeds->n_reflections < 3 ) {
-		printf("BS: Not enough seeds to form a basis\n");
-		return NULL;
-	}
-	printf("BS: Found %i seeds\n", seeds->n_reflections);
-
-	basis = basis_assemble_seeds(seeds);
-	
-	printf("BS: eFOM = %7.3f %%\n", basis_efom(ctx->reflectionlist, basis)*100);
-	
-	return basis;
-
-}
-
diff --git a/src/basis.h b/src/basis.h
index 1d12be1..9fb1cda 100644
--- a/src/basis.h
+++ b/src/basis.h
@@ -1,7 +1,7 @@
 /*
  * basis.h
  *
- * Find approximate lattices to feed various procedures
+ * Handle basis structures
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
  *
@@ -32,7 +32,5 @@ typedef struct basis_struct {
 	Vector c;
 } Basis;
 
-extern Basis *basis_find(ControlContext *ctx);
-
 #endif	/* BASIS_H */
 
diff --git a/src/control.c b/src/control.c
index 4e49de4..e581b06 100644
--- a/src/control.c
+++ b/src/control.c
@@ -14,41 +14,7 @@
 #include <math.h>
 
 #include "control.h"
-
-int control_add_image(ControlContext *ctx, uint16_t *image, int width, int height, double tilt) {
-
-	if ( ctx->images ) {
-		ctx->images = realloc(ctx->images, (ctx->n_images+1)*sizeof(ImageRecord));
-	} else {
-		ctx->images = malloc(sizeof(ImageRecord));
-	}
-	
-	ctx->images[ctx->n_images].tilt = tilt;
-	ctx->images[ctx->n_images].omega = ctx->omega;
-	ctx->images[ctx->n_images].image = image;
-	ctx->images[ctx->n_images].width = width;
-	ctx->images[ctx->n_images].height = height;
-	ctx->images[ctx->n_images].lambda = ctx->lambda;
-	ctx->images[ctx->n_images].fmode = ctx->fmode;
-	ctx->images[ctx->n_images].x_centre = ctx->x_centre;
-	ctx->images[ctx->n_images].y_centre = ctx->y_centre;
-	ctx->images[ctx->n_images].slop = 0.0;
-	
-	if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-		ctx->images[ctx->n_images].pixel_size = ctx->pixel_size;
-		ctx->images[ctx->n_images].camera_len = 0;
-		ctx->images[ctx->n_images].resolution = 0;
-	} else if ( ctx->fmode == FORMULATION_CLEN ) {
-		ctx->images[ctx->n_images].pixel_size = 0;
-		ctx->images[ctx->n_images].camera_len = ctx->camera_length;
-		ctx->images[ctx->n_images].resolution = ctx->resolution;
-	}
-	
-	ctx->n_images++;
-	
-	return ctx->n_images - 1;
-
-}
+#include "image.h"
 
 ControlContext *control_ctx_new() {
 
@@ -61,56 +27,9 @@ ControlContext *control_ctx_new() {
 	ctx->have_centres = 0;
 	ctx->cell = NULL;
 	ctx->dirax = NULL;
-	ctx->n_images = 0;
-	ctx->images = NULL;
+	ctx->images = image_list_new();
 	
 	return ctx;
 
 }
 
-/* Return the minimum (most negative) tilt angle used */
-double control_min_tilt(ControlContext *ctx) {
-
-	int i;
-	double min = 360;
-	
-	for ( i=0; i<ctx->n_images; i++ ) {
-		if ( ctx->images[i].tilt < min ) min = ctx->images[i].tilt;
-	}
-	
-	return min;
-
-}
-
-/* Return the maximum (most positive) tilt angle used */
-double control_max_tilt(ControlContext *ctx) {
-
-	int i;
-	double max = -360;
-	
-	for ( i=0; i<ctx->n_images; i++ ) {
-		if ( ctx->images[i].tilt > max ) max = ctx->images[i].tilt;
-	}
-	
-	return max;
-	
-}
-
-/* Return a reference to the image record with tilt closest to the given value */
-ImageRecord *control_image_nearest_tilt(ControlContext *ctx, double tilt) {
-
-	int i;
-	double dev = 360;
-	ImageRecord *im = NULL;
-	
-	for ( i=0; i<ctx->n_images; i++ ) {
-		if ( fabs(ctx->images[i].tilt - tilt) < dev ) {
-			dev = fabs(ctx->images[i].tilt - tilt);
-			im = &ctx->images[i];
-		}
-	}
-	
-	return im;
-
-}
-
diff --git a/src/control.h b/src/control.h
index 6d318e9..a311b09 100644
--- a/src/control.h
+++ b/src/control.h
@@ -37,39 +37,11 @@ typedef enum {
 	PEAKSEARCH_NONE,
 	PEAKSEARCH_THRESHOLD,
 	PEAKSEARCH_ADAPTIVE_THRESHOLD,
-	PEAKSEARCH_LSQ,
 	PEAKSEARCH_ZAEFFERER,
 	PEAKSEARCH_STAT,
 	PEAKSEARCH_CACHED
 } PeakSearchMode;
 
-typedef struct imagerecord_struct {
-
-	uint16_t	*image;
-	double		tilt;		/* Degrees.  Defines where the pattern lies in reciprocal space */
-	double		omega;		/* Degrees.  Defines where the pattern lies in reciprocal space */
-	double		slop;		/* Degrees.  Defines where the pattern lies "on the negative" */
-	
-	FormulationMode fmode;
-	double		pixel_size;
-	double		camera_len;
-	double		lambda;
-	double		resolution;
-	
-	int		width;
-	int		height;
-	int		x_centre;
-	int		y_centre;
-
-} ImageRecord;
-
-typedef struct {
-	ImageRecord	*parent;
-	int		x;
-	int		y;
-	double		intensity;
-} ImageReflection;
-
 typedef struct cctx_struct {
 	
 	/* Modes */
@@ -101,11 +73,10 @@ typedef struct cctx_struct {
 	unsigned int			lambda_set;
 	
 	/* The input images */
-	int				n_images;
-	ImageRecord			*images;
+	struct imagelist_struct		*images;
 	
-	/* Output */
-	struct reflectionlist_struct	*reflectionlist;
+	/* "Output" */
+	struct reflectionlist_struct	*reflectionlist;	/* Measured reflections (and stuff added to get displayed) */
 	struct dw_struct		*dw;
 	struct basis_struct		*cell;			/* Current estimate of the reciprocal unit cell */
 	struct reflectionlist_struct	*cell_lattice;		/* Reflections calculated from 'cell' */
@@ -141,11 +112,7 @@ typedef struct cctx_struct {
 	
 } ControlContext;
 
-extern int control_add_image(ControlContext *ctx, uint16_t *image, int width, int height, double tilt);
 extern ControlContext *control_ctx_new(void);
-extern double control_min_tilt(ControlContext *ctx);
-extern double control_max_tilt(ControlContext *ctx);
-extern ImageRecord *control_image_nearest_tilt(ControlContext *ctx, double tilt);
 
 #endif	/* CONTROL_H */
 
diff --git a/src/displaywindow.c b/src/displaywindow.c
index d4cb3f7..e9aaa7d 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -521,7 +521,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 		 *	so rotate tilt axis anticlockwise by omega.
 		 *	Since the rotation is about +z, this is already anticlockwise
 		 *	when looking down z. */
-		glRotatef(ctx->images[0].omega, 0.0, 0.0, 1.0);
+		glRotatef(ctx->images->images[0].omega, 0.0, 0.0, 1.0);
 		glScalef(50.0, 1.0, 1.0);
 		glTranslatef(-0.5, 0.0, 0.0);
 		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, yellow);
diff --git a/src/image.c b/src/image.c
new file mode 100644
index 0000000..7205e0b
--- /dev/null
+++ b/src/image.c
@@ -0,0 +1,106 @@
+/*
+ * image.c
+ *
+ * Handle images and image features
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#include <stdlib.h>
+#include <assert.h>
+
+#include "control.h"
+#include "image.h"
+
+int image_add(ImageList *list, uint16_t *image, int width, int height, double tilt, ControlContext *ctx) {
+
+	if ( list->images ) {
+		list->images = realloc(list->images, (list->n_images+1)*sizeof(ImageRecord));
+	} else {
+		assert(list->n_images == 0);
+		list->images = malloc(sizeof(ImageRecord));
+	}
+	
+	list->images[list->n_images].tilt = tilt;
+	list->images[list->n_images].omega = ctx->omega;
+	list->images[list->n_images].image = image;
+	list->images[list->n_images].width = width;
+	list->images[list->n_images].height = height;
+	list->images[list->n_images].lambda = ctx->lambda;
+	list->images[list->n_images].fmode = ctx->fmode;
+	list->images[list->n_images].x_centre = ctx->x_centre;
+	list->images[list->n_images].y_centre = ctx->y_centre;
+	list->images[list->n_images].slop = 0.0;
+	list->images[list->n_images].features = NULL;
+	
+	if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
+		list->images[list->n_images].pixel_size = ctx->pixel_size;
+		list->images[list->n_images].camera_len = 0;
+		list->images[list->n_images].resolution = 0;
+	} else if ( ctx->fmode == FORMULATION_CLEN ) {
+		list->images[list->n_images].pixel_size = 0;
+		list->images[list->n_images].camera_len = ctx->camera_length;
+		list->images[list->n_images].resolution = ctx->resolution;
+	}
+	
+	list->n_images++;
+	
+	return list->n_images - 1;
+
+}
+
+ImageList *image_list_new() {
+
+	ImageList *list;
+	
+	list = malloc(sizeof(ImageList));
+	
+	list->n_images = 0;
+	list->images = NULL;
+	
+	return list;
+
+}
+
+void image_add_feature(ImageFeatureList *flist, double x, double y, ImageRecord *parent, double intensity) {
+
+	if ( flist->features ) {
+		flist->features = realloc(flist->features, (flist->n_features+1)*sizeof(ImageFeature));
+	} else {
+		assert(flist->n_features == 0);
+		flist->features = malloc(sizeof(ImageFeature));
+	}
+	
+	flist->features[flist->n_features].x = x;
+	flist->features[flist->n_features].y = y;
+	flist->features[flist->n_features].intensity = intensity;
+	flist->features[flist->n_features].x = x;
+	flist->features[flist->n_features].parent = parent;
+	
+	flist->n_features++;
+
+}
+
+ImageFeatureList *image_feature_list_new() {
+
+	ImageFeatureList *flist;
+	
+	flist = malloc(sizeof(ImageFeatureList));
+	
+	flist->n_features = 0;
+	flist->features = NULL;
+	
+	return flist;
+}
+
+void image_feature_list_free(ImageFeatureList *flist) {
+
+	if ( !flist ) return;
+	
+	if ( flist->features ) free(flist->features);
+	free(flist);
+
+}
diff --git a/src/image.h b/src/image.h
new file mode 100644
index 0000000..642cc3a
--- /dev/null
+++ b/src/image.h
@@ -0,0 +1,73 @@
+/*
+ * image.h
+ *
+ * Handle images and image features
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+ 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef IMAGE_H
+#define IMAGE_H
+
+#include "control.h"
+
+typedef struct imagefeature_struct {
+
+	struct imagerecord_struct	*parent;
+	int				x;
+	int				y;
+	double				intensity;
+	
+} ImageFeature;
+
+typedef struct {
+
+	ImageFeature		*features;
+	int			n_features;
+
+} ImageFeatureList;
+
+typedef struct imagerecord_struct {
+
+	uint16_t		*image;
+	double			tilt;		/* Degrees.  Defines where the pattern lies in reciprocal space */
+	double			omega;		/* Degrees.  Defines where the pattern lies in reciprocal space */
+	double			slop;		/* Degrees.  Defines where the pattern lies "on the negative" */
+	
+	FormulationMode		fmode;
+	double			pixel_size;
+	double			camera_len;
+	double			lambda;
+	double			resolution;
+	
+	int			width;
+	int			height;
+	int			x_centre;
+	int			y_centre;
+	
+	ImageFeatureList	*features;
+
+} ImageRecord;
+
+typedef struct imagelist_struct {
+
+	int		n_images;
+	ImageRecord	*images;
+
+} ImageList;
+
+extern ImageList *image_list_new(void);
+extern int image_add(ImageList *list, uint16_t *image, int width, int height, double tilt, ControlContext *ctx);
+extern ImageFeatureList *image_feature_list_new(void);
+extern void image_feature_list_free(ImageFeatureList *flist);
+extern void image_add_feature(ImageFeatureList *flist, double x, double y, ImageRecord *parent, double intensity);
+
+#endif	/* IMAGE_H */
+
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index f69c720..981cf06 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -20,6 +20,7 @@
 #include <gtk/gtk.h>
 
 #include "control.h"
+#include "image.h"
 
 typedef enum {
 	IMAGEDISPLAY_NONE		= 0,
diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c
deleted file mode 100644
index cdc2a5c..0000000
--- a/src/itrans-lsq.c
+++ /dev/null
@@ -1,334 +0,0 @@
-/*
- * itrans-lsq.c
- *
- * Peak detection by least-squares fitting
- *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
- *
- *  dtr - Diffraction Tomography Reconstruction
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <gsl/gsl_multifit_nlin.h>
-#include <stdint.h>
-
-#include "control.h"
-#include "imagedisplay.h"
-#include "reflections.h"
-
-#define MAX_LSQ_ITERATIONS 100
-
-typedef struct struct_sampledpoints {
-	gsl_matrix *coordinates;
-	gsl_vector *values;
-	unsigned int n_samples;
-} SampledPoints;
-
-static int itrans_peaksearch_lsq_f(const gsl_vector *gaussian, void *params, gsl_vector *residuals) {
-
-	double cal;
-	double obs;
-	SampledPoints *samples = params;
-	double a, b, c, d, e, f;
-	unsigned int sample;
-	
-	a = gsl_vector_get(gaussian, 0);  b = gsl_vector_get(gaussian, 1);  c = gsl_vector_get(gaussian, 2);
-	d = gsl_vector_get(gaussian, 3);  e = gsl_vector_get(gaussian, 4);  f = gsl_vector_get(gaussian, 5);
-	
-	//printf("Trial parameters: a=%f b=%f c=%f d=%f e=%f f=%f\n", a, b, c, d, e, f);
-	
-	for ( sample=0; sample<samples->n_samples; sample++ ) {
-		int dx = gsl_matrix_get(samples->coordinates, sample, 0);
-		int dy = gsl_matrix_get(samples->coordinates, sample, 1);
-		cal = exp(a + b*dx + c*dy + d*dx*dx + e*dy*dy + f*dx*dy);
-		obs = gsl_vector_get(samples->values, sample);
-		gsl_vector_set(residuals, sample, (cal-obs));
-		//printf("At %3i,%3i: residual=%f - %f = %f\n", dx, dy, cal, obs, cal-obs);
-	}
-	
-	return GSL_SUCCESS;
-
-}
-
-static int itrans_peaksearch_lsq_df(const gsl_vector *gaussian, void *params, gsl_matrix *J) {
-
-	double a, b, c, d, e, f;
-	unsigned int sample;
-	SampledPoints *samples = params;
-	
-	a = gsl_vector_get(gaussian, 0);  b = gsl_vector_get(gaussian, 1);  c = gsl_vector_get(gaussian, 2);
-	d = gsl_vector_get(gaussian, 3);  e = gsl_vector_get(gaussian, 4);  f = gsl_vector_get(gaussian, 5);
-
-	//printf("------a------------b------------c------------d------------e------------f-------\n");
-	for ( sample=0; sample<samples->n_samples; sample++ ) {
-	
-		double ex;
-		int dx, dy;
-		double derivative;
-		
-		dx = gsl_matrix_get(samples->coordinates, sample, 0);
-		dy = gsl_matrix_get(samples->coordinates, sample, 1);
-		ex = exp(a + b*dx + c*dy + d*dx*dx + e*dy*dy + f*dx*dy);
-		
-		derivative = ex;		/* wrt a */
-		gsl_matrix_set(J, sample, 0, derivative);
-		derivative = dx * ex;		/* wrt b */
-		gsl_matrix_set(J, sample, 1, derivative);
-		derivative = dy * ex;		/* wrt c */
-		gsl_matrix_set(J, sample, 2, derivative);
-		derivative = dx * dx * ex;	/* wrt d */
-		gsl_matrix_set(J, sample, 3, derivative);
-		derivative = dy * dy * ex;	/* wrt e */
-		gsl_matrix_set(J, sample, 4, derivative);
-		derivative = dx * dy * ex;	/* wrt f */
-		gsl_matrix_set(J, sample, 5, derivative);
-
-		/*if ( (dx == 0) && (dy == 0) ) {
-			printf("> ");
-		} else {
-			printf("| ");
-		}
-		printf("%10.2e | %10.2e | %10.2e | %10.2e | %10.2e | %10.2e", gsl_matrix_get(J, sample, 0), gsl_matrix_get(J, sample, 1),
-				gsl_matrix_get(J, sample, 2), gsl_matrix_get(J, sample, 3), gsl_matrix_get(J, sample, 4), gsl_matrix_get(J, sample, 5));
-		if ( (dx == 0) && (dy == 0) ) {
-			printf(" <\n");
-		} else {
-			printf(" |\n");
-		}*/
-
-	}
-	//printf("-------------------------------------------------------------------------------\n");
-	
-	return GSL_SUCCESS;
-
-}
-
-static int itrans_peaksearch_lsq_fdf(const gsl_vector *gaussian, void *params, gsl_vector *f, gsl_matrix *J) {
-
-	itrans_peaksearch_lsq_f(gaussian, params, f);
-	itrans_peaksearch_lsq_df(gaussian, params, J);
-	
-	return GSL_SUCCESS;
-
-}
-
-static void itrans_interpolate(uint16_t *image, int width, int x, int y) {
-
-	int a, b, c, d;
-	double av_horiz, av_vert, av;
-	printf("Interpolating...\n");
-	a = image[(x-1) + width*y];  b = image[(x+1) + width*y];
-	c = image[x + width*(y-1)];  d = image[x + width*(y+1)];
-	av_horiz = (a+b)/2;	av_vert = (c+d)/2;
-	av = (av_horiz+av_vert) / 2;
-	image[x + width*y] = av;
-
-}
-
-unsigned int itrans_peaksearch_lsq(ImageRecord *imagerecord, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
-
-	uint16_t max_val = 0;
-	int width, height;
-	uint16_t *image;
-	unsigned int n_reflections = 0;
-	
-	width = imagerecord->width;
-	height = imagerecord->height;
-	image = imagerecord->image;
-
-	/* Least-Squares Craziness. NB Doesn't quite work... */
-	do {
-	
-		int max_x = 0;
-		int max_y = 0;;
-		int x, y;
-		gsl_multifit_fdfsolver *s;
-		gsl_multifit_function_fdf f;
-		unsigned int iter;
-		gsl_vector *gaussian;
-		int iter_status, conv_status;
-		gsl_matrix *covar;
-		SampledPoints samples;
-		int sample;
-		double ga, gb, gc, gd, ge, gf;
-		gboolean sanity;
-		int dx, dy;
-		int bb_lh, bb_rh, bb_top, bb_bot;
-		uint16_t ival;
-		double sx, sy;
-	
-		/* Locate the highest point */
-		max_val = 0;
-		for ( y=10; y<height-10; y++ ) {
-			for ( x=10; x<width-10; x++ ) {
-
-				if ( image[x + width*y] > max_val ) {
-					max_val = image[x + width*y];
-					max_x = x;  max_y = y;
-				}
-			
-			}
-		}
-		x = max_x; y = max_y;
-		
-		/* Try fitting a Gaussian to this region and see what happens... */
-		f.p = 6;
-		
-		/* Set initial estimate of the fit parameters.  I = exp(a + bx + cy + dx^2 + ey^2 + fxy) */
-		gaussian = gsl_vector_alloc(f.p);
-		gsl_vector_set(gaussian, 0, log(max_val));	/* a */
-		gsl_vector_set(gaussian, 1, -0.01);		/* b */
-		gsl_vector_set(gaussian, 2, -0.01);		/* c */
-		gsl_vector_set(gaussian, 3, -0.01);		/* d */
-		gsl_vector_set(gaussian, 4, -0.01);		/* e */
-		gsl_vector_set(gaussian, 5, -0.01);		/* f */
-		
-		/* Determine the bounding box which contains most of the peak */
-		/* Left */	dx = 0;	do { ival = image[(x+dx) + width*(y+0)]; dx--; } while ( (ival>max_val/300) && (x+dx>0) );	bb_lh = dx;
-		/* Right */	dx = 0;	do { ival = image[(x+dx) + width*(y+0)]; dx++; } while ( (ival>max_val/300) && (x+dx<width) );	bb_rh = dx;
-		/* Top */	dy = 0;	do { ival = image[(x+0) + width*(y+dy)]; dy++; } while ( (ival>max_val/300) && (y+dy<height) );	bb_top = dy;
-		/* Bottom */	dy = 0;	do { ival = image[(x+0) + width*(y+dy)]; dy--; } while ( (ival>max_val/300) && (y+dy>0) );	bb_bot = dy;
-		printf("Peak %i,%i: bounding box %i<dx<%i, %i<dy<%i\n", x, y, bb_lh, bb_rh, bb_bot, bb_top);
-		sanity = TRUE;
-		if ( (bb_rh-bb_lh) > 300 ) sanity = FALSE;
-		if ( (bb_top-bb_bot) > 300 ) sanity = FALSE;
-		if ( sanity ) {
-		
-			/* Choose which points inside this bounding box to use for curve fitting */
-			samples.n_samples = ((bb_top-bb_bot)+1)*((bb_rh-bb_lh)+1);
-			samples.coordinates = gsl_matrix_alloc(samples.n_samples, 2);
-			sample = 0;
-			for ( sx=bb_lh; sx<=bb_rh; sx++ ) {
-				for ( sy=bb_bot; sy<=bb_top; sy++ ) {
-					gsl_matrix_set(samples.coordinates, sample, 0, sx);	gsl_matrix_set(samples.coordinates, sample, 1, sy);
-					sample++;
-				}
-			}
-			
-			samples.values = gsl_vector_alloc(samples.n_samples);
-			f.n = samples.n_samples;
-			for ( sample=0; sample<samples.n_samples; sample++ ) {
-				int dx = gsl_matrix_get(samples.coordinates, sample, 0);
-				int dy = gsl_matrix_get(samples.coordinates, sample, 1);
-				gsl_vector_set(samples.values, sample, image[(x+dx) + width*(y+dy)]);
-				//printf("At %3i,%3i: value=%i\n", dx, dy, image[(x+dx) + width*(y+dy)]);
-			}
-			f.params = &samples;
-			
-			/* Execute the LSQ fitting procedure */
-			s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, f.n, f.p);
-			f.f = &itrans_peaksearch_lsq_f;	f.df = &itrans_peaksearch_lsq_df;	f.fdf = &itrans_peaksearch_lsq_fdf;
-			gsl_multifit_fdfsolver_set(s, &f, gaussian);
-			iter = 0; conv_status = GSL_CONTINUE;
-			do {
-			
-				/* Iterate */
-				iter_status = gsl_multifit_fdfsolver_iterate(s);
-				iter++;		
-				
-				/* Check for error */
-				if ( iter_status ) {
-					break;
-				}
-				
-				/* Test for convergence */
-				conv_status = gsl_multifit_test_delta(s->dx, s->x, 1e-6, 1e-6);
-										
-			} while ( (conv_status == GSL_CONTINUE) && (iter < MAX_LSQ_ITERATIONS) );
-			
-			/* See how good the fit is */
-			covar = gsl_matrix_alloc(f.p, f.p);
-			gsl_multifit_covar(s->J, 0.0, covar);
-			ga = gsl_vector_get(s->x, 0);  gb = gsl_vector_get(s->x, 1);  gc = gsl_vector_get(s->x, 2);
-			gd = gsl_vector_get(s->x, 3);  ge = gsl_vector_get(s->x, 4);  gf = gsl_vector_get(s->x, 5);
-			if ( fabs(exp(ga + gb*100 + gc*100 + gd*100*100 + ge*100*100 + gf*100*100)) > 10 ) {
-				printf("Failed sanity check: %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", x, y, ga, gb, gc, gd, ge, gf);
-				sanity = FALSE;
-			} else {
-				sanity = TRUE;
-			}
-			if ( (conv_status == GSL_SUCCESS) && (!iter_status) && (sanity) ) {
-				
-				/* Good fit! */
-				int dx, dy;
-				int bb_top, bb_bot, bb_lh, bb_rh;
-				double frac;
-				double brightness;
-			
-				printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n",
-						iter, x, y, ga, gb, gc, gd, ge, gf);
-				brightness = image[x + width*y];
-				reflection_add_from_dp(ctx, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, brightness);
-				n_reflections++;
-				
-				/* Remove this peak from the image */
-				/* Find right-hand edge of bounding box */
-				dx = 0;	do {
-					double pval, ival;
-					pval = exp(ga + gb*dx + gc*0 + gd*dx*dx + ge*0*0 + gf*dx*-0);
-					ival = image[(x+dx) + width*(y+0)];
-					frac = pval / ival;
-					dx++;
-				} while ( (frac > 0.1) && (x+dx<width) );
-				bb_rh = dx;
-				/* Find left-hand edge of bounding box */
-				dx = 0;	do {
-					double pval, ival;
-					pval = exp(ga + gb*dx + gc*0 + gd*dx*dx + ge*0*0 + gf*dx*-0);
-					ival = image[(x+dx) + width*(y+0)];
-					frac = pval / ival;
-					dx--;
-				} while ( (frac > 0.1) && (x+dx>0) );
-				bb_lh = dx;
-				/* Find top edge of bounding box */
-				dy = 0;	do {
-					double pval, ival;
-					pval = exp(ga + gb*0 + gc*dy + gd*0*0 + ge*dy*dy + gf*0*dy);
-					ival = image[(x+0) + width*(y+dy)];
-					frac = pval / ival;
-					dy++;
-				} while ( (frac > 0.1) && (y+dy<height) );
-				bb_top = dy;
-				/* Find bottom edge of bounding box */
-				dy = 0;	do {
-					double pval, ival;
-					pval = exp(ga + gb*0 + gc*dy + gd*0*0 + ge*dy*dy + gf*0*dy);
-					ival = image[(x+0) + width*(y+dy)];
-					frac = pval / ival;
-					dy--;
-				} while ( (frac > 0.1) && (y+dy>0) );
-				bb_bot = dy;
-				printf("Fitted peak bounding box %i<dx<%i, %i<dy<%i\n", bb_lh, bb_rh, bb_bot, bb_top);
-				for ( dx=bb_lh; dx<bb_rh; dx++ ) {
-					for ( dy=bb_bot; dy<bb_top; dy++ ) {
-						double pval;
-						pval = exp(ga + gb*dx + gc*dy + gd*dx*dx + ge*dy*dy + gf*dx*dy);
-						image[(x+dx) + width*(y+dy)] -= pval;
-					}
-				}
-				
-			} else {
-				itrans_interpolate(image, width, x, y);
-			}
-			
-			gsl_matrix_free(covar);
-			gsl_multifit_fdfsolver_free(s);
-			gsl_vector_free(gaussian);
-			gsl_matrix_free(samples.coordinates);
-			gsl_vector_free(samples.values);
-		
-		} else {
-			printf("Failed bounding box sanity check\n");
-			itrans_interpolate(image, width, x, y);
-		}
-		
-	
-	} while ( max_val > 50 );
-	
-	return n_reflections;
-}
-
diff --git a/src/itrans-lsq.h b/src/itrans-lsq.h
deleted file mode 100644
index 24f4e98..0000000
--- a/src/itrans-lsq.h
+++ /dev/null
@@ -1,24 +0,0 @@
-/*
- * itrans-lsq.h
- *
- * Peak detection by least-squares fitting
- *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
- *
- *  dtr - Diffraction Tomography Reconstruction
- *
- */
- 
-#ifndef ITRANS_LSQ_H
-#define ITRANS_LSQ_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "control.h"
-
-extern unsigned int itrans_peaksearch_lsq(ImageRecord *imagerecord, ControlContext *ctx);
-
-#endif	/* ITRANS_LSQ_H */
-
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index 21cf415..e08bd26 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -500,7 +500,7 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c
 
 }
 
-unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ctx) {
+ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) {
 
 	unsigned int n_reflections;	
 	gsl_matrix *m;
@@ -508,6 +508,9 @@ unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ct
 	int i;
 	double px,py;
 	uint16_t *image = imagerecord->image;
+	ImageFeatureList *flist;
+	
+	flist = image_feature_list_new();
 
 	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord->width, imagerecord->height);
 	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
@@ -516,13 +519,14 @@ unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ct
 	
 		px = gsl_matrix_get(p,0,i);
 		py = gsl_matrix_get(p,1,i);
-		reflection_add_from_dp(ctx, (px-imagerecord->x_centre), (py-imagerecord->y_centre), imagerecord, 1.0);
+		image_add_feature(flist, (px-imagerecord->x_centre), (py-imagerecord->y_centre), imagerecord, 1.0);
 		
 	}
+	
 	gsl_matrix_free(m);
 	gsl_matrix_free(p);
 	
-	return n_reflections;
+	return flist;
 	
 }
 
diff --git a/src/itrans-stat.h b/src/itrans-stat.h
index f475a2c..75daa98 100644
--- a/src/itrans-stat.h
+++ b/src/itrans-stat.h
@@ -17,8 +17,9 @@
 #include <config.h>
 #endif
 
-#include "control.h"
+#include "image.h"
 
-unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ctx);
+extern ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord);
 
 #endif /* ITRANS_STAT_H */
+
diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c
index 30e476c..17a2cbe 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -16,17 +16,17 @@
 #include <stdint.h>
 #include <assert.h>
 
-#include "control.h"
-#include "imagedisplay.h"
-#include "reflections.h"
+#include "image.h"
 
-unsigned int itrans_peaksearch_threshold(ImageRecord *imagerecord, ControlContext *ctx) {
+ImageFeatureList *itrans_peaksearch_threshold(ImageRecord *imagerecord, ControlContext *ctx) {
 
 	int width, height;
 	int x, y;
-	unsigned int n_reflections = 0;
 	uint16_t *image = imagerecord->image;
 	uint16_t max = 0;
+	ImageFeatureList *flist;
+	
+	flist = image_feature_list_new();
 	
 	width = imagerecord->width;
 	height = imagerecord->height;
@@ -45,25 +45,25 @@ unsigned int itrans_peaksearch_threshold(ImageRecord *imagerecord, ControlContex
 				assert(y<height);
 				assert(x>=0);
 				assert(y>=0);
-				reflection_add_from_dp(ctx, (x-imagerecord->x_centre), (y-imagerecord->y_centre),
-							imagerecord, image[x + width*y]);
-				n_reflections++;
+				image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]);
 			}
 		}
 	}
 	
-	return n_reflections;
+	return flist;
 	
 }
 
-unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, ControlContext *ctx) {
+ImageFeatureList *itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, ControlContext *ctx) {
 
 	uint16_t max_val = 0;
 	int width, height;
-	unsigned int n_reflections = 0;
 	uint16_t *image;
 	uint16_t max;
 	int x, y;
+	ImageFeatureList *flist;
+	
+	flist = image_feature_list_new();
 	
 	image = imagerecord->image;
 	width = imagerecord->width;
@@ -100,9 +100,7 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, Cont
 			assert(max_y<height);
 			assert(max_x>=0);
 			assert(max_y>=0);
-			reflection_add_from_dp(ctx, (max_x-imagerecord->x_centre), (max_y-imagerecord->y_centre),
-							imagerecord, image[x + width*y]);
-			n_reflections++;
+			image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]);
 			
 			/* Remove it and its surroundings */
 			for ( y=((max_y-10>0)?max_y-10:0); y<((max_y+10)<height?max_y+10:height); y++ ) {
@@ -114,7 +112,7 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, Cont
 		
 	} while ( max_val > 50 );
 	
-	return n_reflections;
+	return flist;
 
 }
 
diff --git a/src/itrans-threshold.h b/src/itrans-threshold.h
index f19ea57..496ded9 100644
--- a/src/itrans-threshold.h
+++ b/src/itrans-threshold.h
@@ -16,10 +16,10 @@
 #include <config.h>
 #endif
 
-#include "control.h"
+#include "image.h"
 
-extern unsigned int itrans_peaksearch_threshold(ImageRecord *image, ControlContext *ctx);
-extern unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *image, ControlContext *ctx);
+extern ImageFeatureList *itrans_peaksearch_threshold(ImageRecord *image);
+extern ImageFeatureList *itrans_peaksearch_adaptive_threshold(ImageRecord *image);
 
 #endif	/* ITRANS_THRESHOLD_H */
 
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index 889aac4..d984774 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -16,25 +16,24 @@
 #include <stdint.h>
 #include <assert.h>
 
-#include "control.h"
-#include "imagedisplay.h"
-#include "reflections.h"
 #include "utils.h"
+#include "image.h"
 
 #define PEAK_WINDOW_SIZE 20
 
-unsigned int itrans_peaksearch_zaefferer(ImageRecord *imagerecord, ControlContext *ctx) {
+ImageFeatureList *itrans_peaksearch_zaefferer(ImageRecord *imagerecord) {
 
 	int x, y;
-	unsigned int n_reflections;
 	int width, height;
 	uint16_t *image;
+	ImageFeatureList *flist;
+	
+	flist = image_feature_list_new();
 	
 	image = imagerecord->image;
 	width = imagerecord->width;
 	height = imagerecord->height;
 	
-	n_reflections = 0;
 	for ( x=1; x<width-1; x++ ) {
 		for ( y=1; y<height-1; y++ ) {
 
@@ -85,16 +84,15 @@ unsigned int itrans_peaksearch_zaefferer(ImageRecord *imagerecord, ControlContex
 					assert(mask_y<height);
 					assert(mask_x>=0);
 					assert(mask_y>=0);
-					reflection_add_from_dp(ctx, (mask_x-imagerecord->x_centre), (mask_y-imagerecord->y_centre),
-							imagerecord, image[mask_x + width*mask_y]);
-					n_reflections++;
+					image_add_feature(flist, (mask_x-imagerecord->x_centre), (mask_y-imagerecord->y_centre),
+									imagerecord, image[mask_x + width*mask_y]);
 				}
 				
 			}
 		}
 	}
 
-	return n_reflections;
+	return flist;
 
 }
 
diff --git a/src/itrans-zaefferer.h b/src/itrans-zaefferer.h
index 7a3e7ed..dbb819b 100644
--- a/src/itrans-zaefferer.h
+++ b/src/itrans-zaefferer.h
@@ -16,9 +16,9 @@
 #include <config.h>
 #endif
 
-#include "control.h"
+#include "image.h"
 
-extern unsigned int itrans_peaksearch_zaefferer(ImageRecord *image, ControlContext *ctx);
+extern ImageFeatureList *itrans_peaksearch_zaefferer(ImageRecord *image);
 
 #endif	/* ITRANS_ZAEFFERER_H */
 
diff --git a/src/itrans.c b/src/itrans.c
index 05d2793..2d99bbb 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -14,24 +14,19 @@
 #include <config.h>
 #endif
 
-#include "control.h"
-#include "reflections.h"
+#include "image.h"
 #include "itrans-threshold.h"
 #include "itrans-zaefferer.h"
-#include "itrans-lsq.h"
 #include "itrans-stat.h"
 
-void itrans_process_image(ImageRecord *image, ControlContext *ctx) {
+ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode) {
 	
-	unsigned int n_reflections;
-	
-	switch ( ctx->psmode ) {
-		case PEAKSEARCH_THRESHOLD : itrans_peaksearch_threshold(image, ctx); break;
-		case PEAKSEARCH_ADAPTIVE_THRESHOLD : itrans_peaksearch_adaptive_threshold(image, ctx); break;
-		case PEAKSEARCH_LSQ : itrans_peaksearch_lsq(image, ctx); break;
-		case PEAKSEARCH_ZAEFFERER : itrans_peaksearch_zaefferer(image, ctx); break;
-		case PEAKSEARCH_STAT : itrans_peaksearch_stat(image, ctx); break;
-		default: n_reflections = 0;
+	switch ( psmode ) {
+		case PEAKSEARCH_THRESHOLD : return itrans_peaksearch_threshold(image);
+		case PEAKSEARCH_ADAPTIVE_THRESHOLD : return itrans_peaksearch_adaptive_threshold(image);
+		case PEAKSEARCH_ZAEFFERER : return itrans_peaksearch_zaefferer(image);
+		case PEAKSEARCH_STAT : return itrans_peaksearch_stat(image);
+		default: return NULL;
 	}
 	
 }
diff --git a/src/itrans.h b/src/itrans.h
index a293a98..ec2e852 100644
--- a/src/itrans.h
+++ b/src/itrans.h
@@ -18,7 +18,7 @@
 
 #include "control.h"
 
-extern void itrans_process_image(ImageRecord *image, ControlContext *ctx);
+extern ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode);
 
 #endif	/* ITRANS_H */
 
diff --git a/src/main.c b/src/main.c
index 5365601..0a8362c 100644
--- a/src/main.c
+++ b/src/main.c
@@ -37,11 +37,11 @@
 void main_do_reconstruction(ControlContext *ctx) {
 
 	if ( ctx->cache_filename ) {
-		prealign_sum_stack(ctx);
+		prealign_sum_stack(ctx->images, ctx->have_centres);
 		ctx->reflectionlist = cache_load(ctx->cache_filename);
 		printf("Loading cached reflections from '%s'\n", ctx->cache_filename);
 	} else if ( ctx->inputfiletype != INPUT_DRX ) {
-		prealign_sum_stack(ctx);
+		prealign_sum_stack(ctx->images, ctx->have_centres);
 		mapping_create(ctx);
 	}
 	
@@ -63,10 +63,9 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 		switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) ) {
 			case 0 : ctx->psmode = PEAKSEARCH_THRESHOLD; break;
 			case 1 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD; break;
-			case 2 : ctx->psmode = PEAKSEARCH_LSQ; break;
-			case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
-			case 4 : ctx->psmode = PEAKSEARCH_STAT; break;
-			case 5 : ctx->psmode = PEAKSEARCH_CACHED; break;
+			case 2 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
+			case 3 : ctx->psmode = PEAKSEARCH_STAT; break;
+			case 4 : ctx->psmode = PEAKSEARCH_CACHED; break;
 			default: ctx->psmode = PEAKSEARCH_NONE; break;	/* This happens when reading from a cache file */
 		}
 		
@@ -112,7 +111,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 
 static gint main_peaksearch_changed(GtkWidget *method_window, ControlContext *ctx) {
 
-	if ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) == 5 ) {
+	if ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) == 4 ) {
 		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector), TRUE);
 	} else {
 		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector), FALSE);
@@ -149,7 +148,6 @@ void main_method_dialog_open(ControlContext *ctx) {
 	ctx->combo_peaksearch = gtk_combo_box_new_text();
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Simple Thresholding");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Adaptive Thresholding");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Least-Squares Fit");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical Analysis");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get From Cache File");
diff --git a/src/mapping.c b/src/mapping.c
index c2fed73..a4ecd96 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -12,29 +12,117 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 
 #include "reflections.h"
 #include "control.h"
 #include "itrans.h"
+#include "image.h"
+
+static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
+
+	/* "Input" space */
+	double d, x, y;
+	ImageRecord *imagerecord;
+	
+	/* Angular description of reflection */
+	double theta, psi, k;
+	
+	/* Reciprocal space */
+	double tilt;
+	double omega;
+	
+	double x_temp, y_temp, z_temp;
+	double nx, ny, nz;
+	
+	imagerecord = refl->parent;
+	x = refl->x;	y = refl->y;
+	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
+	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
+	k = 1/imagerecord->lambda;
+	
+	/* Calculate an angular description of the reflection */
+	if ( imagerecord->fmode == FORMULATION_CLEN ) {
+		x /= imagerecord->resolution;
+		y /= imagerecord->resolution;	/* Convert pixels to metres */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, imagerecord->camera_len);
+	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
+		x *= imagerecord->pixel_size;
+		y *= imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, k);
+	} else {
+		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
+		return -1;
+	}
+	psi = atan2(y, x);
+	
+	x_temp = k*sin(theta)*cos(psi);
+	y_temp = k*sin(theta)*sin(psi);
+	z_temp = k - k*cos(theta);
+	
+	/* Apply the rotations...
+		First: rotate image clockwise until tilt axis is aligned horizontally. */
+	nx = x_temp*cos(omega) + y_temp*sin(omega);
+	ny = -x_temp*sin(omega) + y_temp*cos(omega);
+	nz = z_temp;
+	
+	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
+		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
+	x_temp = nx; y_temp = ny; z_temp = nz;	
+	nx = x_temp;
+	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
+	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
+	
+	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
+	x_temp = nx; y_temp = ny; z_temp = nz;	
+	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
+	nz = z_temp;
+	
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
+	
+	return 0;
+
+}
 
 ReflectionList *mapping_create(ControlContext *ctx) {
 
-	ReflectionList *reflectionlist;
 	int i;
 	
-	/* Create reflection context */
-	reflectionlist = reflectionlist_new();
+	/* Pass all images through itrans */
+	printf("MP: Analysing images..."); fflush(stdout);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
+	}
+	printf("done.\n");
+	
+	/* Create reflection list for measured reflections */
+	ctx->reflectionlist = reflectionlist_new();
+	printf("MP: Mapping to 3D..."); fflush(stdout);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
 	
-	/* Pass all images through itrans
-	 *   (let itrans add the reflections to reflectionlist for now) */
-	printf("MP: Processing images..."); fflush(stdout);
-	ctx->reflectionlist = reflectionlist;
-	for ( i=0; i<ctx->n_images; i++ ) {
-		itrans_process_image(&ctx->images[i], ctx);
+		int j;
+		
+		/* Iterate over the features in this image */
+		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
+		
+			double nx, ny, nz, twotheta;
+			
+			if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
+				reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
+			}
+
+		}
+		
 	}
-	printf("done: %i 'measured' reflections\n", ctx->reflectionlist->n_reflections);
+	printf("done.\n");
 	
-	return reflectionlist;
+	return ctx->reflectionlist;
 
 }
 
diff --git a/src/mrc.c b/src/mrc.c
index 31c0a5b..880c159 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -24,6 +24,7 @@
 #include "itrans.h"
 #include "reflections.h"
 #include "utils.h"
+#include "image.h"
 
 int mrc_read(ControlContext *ctx) {
 
@@ -88,12 +89,11 @@ int mrc_read(ControlContext *ctx) {
 			for ( y=0; y<mrc.ny; y++ ) {
 				uimage[x + mrc.nx*y] = image[x + mrc.nx*y] + 32767;
 				if ( uimage[x + mrc.nx*y] < 10 ) printf("%i\n", uimage[x + mrc.nx*y]);
-				
 			}
 		}
 		free(image);
 		
-		control_add_image(ctx, uimage, mrc.nx, mrc.ny, ext[i].a_tilt);
+		image_add(ctx->images, uimage, mrc.nx, mrc.ny, ext[i].a_tilt, ctx);
 		
 	}
 
diff --git a/src/prealign.c b/src/prealign.c
index b72943c..ef23c6e 100644
--- a/src/prealign.c
+++ b/src/prealign.c
@@ -17,6 +17,7 @@
 #include "control.h"
 #include "imagedisplay.h"
 #include "main.h"
+#include "image.h"
 
 typedef struct {
 	int		n;
@@ -36,18 +37,18 @@ static gint prealign_clicked(GtkWidget *widget, GdkEventButton *event, PreAlignB
 	x -= xoffs;	y -= yoffs;
 	x /= scale;	y /= scale;
 	y = pb->id->imagerecord.height - y;
-	pb->ctx->images[pb->n].x_centre = x;
-	pb->ctx->images[pb->n].y_centre = y;
+	pb->ctx->images->images[pb->n].x_centre = x;
+	pb->ctx->images->images[pb->n].y_centre = y;
 	
 	pb->n++;
-	if ( pb->n >= pb->ctx->n_images ) {
+	if ( pb->n >= pb->ctx->images->n_images ) {
 		/* Finished */
 		imagedisplay_close(pb->id);
 		main_do_reconstruction(pb->ctx);
 		free(pb);
 	} else {
 		/* Display the next pattern */
-		imagedisplay_put_data(pb->id, pb->ctx->images[pb->n]);
+		imagedisplay_put_data(pb->id, pb->ctx->images->images[pb->n]);
 	}
 	
 	return 0;
@@ -65,13 +66,13 @@ void prealign_do_series(ControlContext *ctx) {
 	pb = malloc(sizeof(PreAlignBlock));
 	pb->n = 0;
 	pb->ctx = ctx;
-	pb->id = imagedisplay_open_with_message(ctx->images[pb->n], "Image Pre-alignment",
+	pb->id = imagedisplay_open_with_message(ctx->images->images[pb->n], "Image Pre-alignment",
 				"Click the centre of the zero-order beam as accurately as you can.", IMAGEDISPLAY_QUIT_IF_CLOSED,
 				G_CALLBACK(prealign_clicked), pb);
 
 }
 
-void prealign_sum_stack(ControlContext *ctx) {
+void prealign_sum_stack(ImageList *list, int have_centres) {
 
 	int twidth, theight;
 	int mnorth, msouth, mwest, meast;
@@ -82,11 +83,11 @@ void prealign_sum_stack(ControlContext *ctx) {
 	
 	/* Determine maximum size of image to accommodate, and allocate memory */
 	mnorth = 0;  msouth = 0;  mwest = 0;  meast = 0;
-	for ( i=0; i<ctx->n_images; i++ ) {
-		if ( ctx->images[i].width-ctx->images[i].x_centre > meast ) meast = ctx->images[i].width-ctx->images[i].x_centre;
-		if ( ctx->images[i].x_centre > mwest ) mwest = ctx->images[i].x_centre;
-		if ( ctx->images[i].height-ctx->images[i].y_centre > mnorth ) mnorth = ctx->images[i].height-ctx->images[i].y_centre;
-		if ( ctx->images[i].y_centre > msouth ) msouth = ctx->images[i].y_centre;
+	for ( i=0; i<list->n_images; i++ ) {
+		if ( list->images[i].width-list->images[i].x_centre > meast ) meast = list->images[i].width-list->images[i].x_centre;
+		if ( list->images[i].x_centre > mwest ) mwest = list->images[i].x_centre;
+		if ( list->images[i].height-list->images[i].y_centre > mnorth ) mnorth = list->images[i].height-list->images[i].y_centre;
+		if ( list->images[i].y_centre > msouth ) msouth = list->images[i].y_centre;
 	}
 	twidth = mwest + meast;
 	theight = mnorth + msouth;
@@ -95,23 +96,23 @@ void prealign_sum_stack(ControlContext *ctx) {
 	memset(image_total, 0, twidth * theight * sizeof(uint16_t));
 	
 	/* Add the image stack together */
-	if ( !ctx->have_centres ) {
+	if ( !have_centres ) {
 		
 		int max_x, max_y;
 		uint16_t max_val;
 		
-		for ( i=0; i<ctx->n_images; i++ ) {
+		for ( i=0; i<list->n_images; i++ ) {
 			int xoffs, yoffs;
-			xoffs = (twidth - ctx->images[i].width)/2;
-			yoffs = (theight - ctx->images[i].height)/2;
-			for ( y=0; y<ctx->images[i].height; y++ ) {
-				for ( x=0; x<ctx->images[i].width; x++ ) {
+			xoffs = (twidth - list->images[i].width)/2;
+			yoffs = (theight - list->images[i].height)/2;
+			for ( y=0; y<list->images[i].height; y++ ) {
+				for ( x=0; x<list->images[i].width; x++ ) {
 					assert(x+xoffs < twidth);
 					assert(y+yoffs < theight);
 					assert(x+xoffs >= 0);
 					assert(y+yoffs >= 0);
 					image_total[(x+xoffs) + twidth*(y+yoffs)] +=
-						ctx->images[i].image[x + ctx->images[i].width*y]/ctx->n_images;
+						list->images[i].image[x + list->images[i].width*y]/list->n_images;
 				}
 			}
 		}
@@ -128,30 +129,31 @@ void prealign_sum_stack(ControlContext *ctx) {
 		}
 		
 		/* Record this measurement on all images */
-		for ( i=0; i<ctx->n_images; i++ ) {
-			ctx->images[i].x_centre = max_x;
-			ctx->images[i].y_centre = max_y;
+		for ( i=0; i<list->n_images; i++ ) {
+			list->images[i].x_centre = max_x;
+			list->images[i].y_centre = max_y;
 		}
 		total_record.x_centre = max_x;
 		total_record.y_centre = max_y;
 	
 	} else {
 		
-		for ( i=0; i<ctx->n_images; i++ ) {
+		for ( i=0; i<list->n_images; i++ ) {
 		
 			int xoffs, yoffs;
 			
-			xoffs = mwest - ctx->images[i].x_centre;
-			yoffs = msouth - ctx->images[i].y_centre;
+			xoffs = mwest - list->images[i].x_centre;
+			yoffs = msouth - list->images[i].y_centre;
 			
-			for ( y=0; y<ctx->images[i].height; y++ ) {
-				for ( x=0; x<ctx->images[i].width; x++ ) {
+			for ( y=0; y<list->images[i].height; y++ ) {
+				for ( x=0; x<list->images[i].width; x++ ) {
 					assert(x+xoffs < twidth);
 					assert(y+yoffs < theight);
 					assert(x+xoffs >= 0);
 					assert(y+yoffs >= 0);
 					image_total[(x+xoffs) + twidth*(y+yoffs)] +=
-						ctx->images[i].image[x + ctx->images[i].width*y]/ctx->n_images;
+						list->images[i].image[x + list->images[i].width*y]/list->n_images;
+					total_record.omega = list->images[i].omega;
 				}
 			}
 			
@@ -165,7 +167,6 @@ void prealign_sum_stack(ControlContext *ctx) {
 	total_record.image = image_total;
 	total_record.width = twidth;
 	total_record.height = theight;
-	total_record.omega = ctx->omega;
 	sum_id = imagedisplay_open(total_record, "Sum of All Images", IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS);
 
 }
diff --git a/src/prealign.h b/src/prealign.h
index cc4c4b5..c21390b 100644
--- a/src/prealign.h
+++ b/src/prealign.h
@@ -16,10 +16,10 @@
 #ifndef PREALIGN_H
 #define PREALIGN_H
 
-#include "control.h"
+#include "image.h"
 
 extern void prealign_do_series(ControlContext *ctx);
-extern void prealign_sum_stack(ControlContext *ctx);
+extern void prealign_sum_stack(ImageList *list, int have_centres);
 
 #endif	/* PREALIGN_H */
 
diff --git a/src/readpng.c b/src/readpng.c
index 41d589d..d310abb 100644
--- a/src/readpng.c
+++ b/src/readpng.c
@@ -144,7 +144,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
 	fclose(fh);
 
-	control_add_image(ctx, image, width, height, tilt);
+	image_add(ctx->images, image, width, height, tilt, ctx);
 	
 	return 0;
 	
diff --git a/src/reflections.c b/src/reflections.c
index 0005332..7fcec2f 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -126,94 +126,6 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 	
 }
 
-int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
-
-	/* "Input" space */
-	double d, x, y;
-	ImageRecord *imagerecord;
-	
-	/* Angular description of reflection */
-	double theta, psi, k;
-	
-	/* Reciprocal space */
-	double tilt;
-	double omega;
-	
-	double x_temp, y_temp, z_temp;
-	double nx, ny, nz;
-	
-	imagerecord = refl->parent;
-	x = refl->x;	y = refl->y;
-	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
-	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
-	k = 1/imagerecord->lambda;
-	
-	/* Calculate an angular description of the reflection */
-	if ( imagerecord->fmode == FORMULATION_CLEN ) {
-		x /= imagerecord->resolution;
-		y /= imagerecord->resolution;	/* Convert pixels to metres */
-		d = sqrt((x*x) + (y*y));
-		theta = atan2(d, imagerecord->camera_len);
-	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
-		x *= imagerecord->pixel_size;
-		y *= imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
-		d = sqrt((x*x) + (y*y));
-		theta = atan2(d, k);
-	} else {
-		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
-		return -1;
-	}
-	psi = atan2(y, x);
-	
-	x_temp = k*sin(theta)*cos(psi);
-	y_temp = k*sin(theta)*sin(psi);
-	z_temp = k - k*cos(theta);
-	
-	/* Apply the rotations...
-		First: rotate image clockwise until tilt axis is aligned horizontally. */
-	nx = x_temp*cos(omega) + y_temp*sin(omega);
-	ny = -x_temp*sin(omega) + y_temp*cos(omega);
-	nz = z_temp;
-	
-	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
-		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
-	x_temp = nx; y_temp = ny; z_temp = nz;	
-	nx = x_temp;
-	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
-	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
-	
-	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
-	x_temp = nx; y_temp = ny; z_temp = nz;	
-	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
-	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
-	nz = z_temp;
-	
-	*ddx = nx;
-	*ddy = ny;
-	*ddz = nz;
-	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
-	
-	return 0;
-
-}
-
-/* x and y in pixels, measured from centre of image */
-void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord *imagerecord, double intensity) {
-
-	double nx, ny, nz, twotheta;
-	ImageReflection refl;
-	
-	refl.parent = imagerecord;
-	refl.x = x;
-	refl.y = y;
-	refl.intensity = intensity;
-	
-	if ( !reflection_map_to_space(&refl, &nx, &ny, &nz, &twotheta) ) {
-		reflection_add(ctx->reflectionlist, nx, ny, nz, intensity, REFLECTION_NORMAL);
-	}
-
-}
-
 void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r) {
 	if ( reflectionlist->last_reflection ) {
 		reflectionlist->last_reflection->next = r;
@@ -315,42 +227,6 @@ Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, dou
 
 }
 
-/* This destroys the lfom values in the input list */
-ReflectionList *reflectionlist_sort_lfom(ReflectionList *in) {
-
-	ReflectionList *out;
-	Reflection *best;
-	
-	out = reflectionlist_new();
-	
-	do {
-	
-		Reflection *reflection;
-		int lfom = 0;
-		
-		reflection = in->reflections;
-		best = NULL;
-		while ( reflection ) {
-			if ( reflection->lfom > lfom ) {
-				best = reflection;
-				lfom = reflection->lfom;
-			}
-			reflection = reflection->next;
-		};
-		
-		if ( best ) {
-			Reflection *new;
-			new = reflection_add(out, best->x, best->y, best->z, best->intensity, best->type);
-			new->lfom = best->lfom;
-			best->lfom = 0;
-		}
-		
-	} while ( best );
-	
-	return out;
-
-}
-
 /* Generate a list of reflections from a unit cell */
 ReflectionList *reflection_list_from_cell(Basis *basis) {
 
diff --git a/src/reflections.h b/src/reflections.h
index 8347574..2cc7f10 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -62,16 +62,14 @@ extern void reflectionlist_clear_markers(ReflectionList *reflectionlist);
 extern void reflectionlist_free(ReflectionList *reflectionlist);
 
 #include "control.h"
+#include "image.h"
 extern Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, double z, double intensity, ReflectionType type);
-extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord *imagerecord, double intensity);
 extern void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r);
-extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta);
 
 extern double reflectionlist_largest_g(ReflectionList *reflectionlist);
 extern Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x, double y, double z);
 extern Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance);
 extern Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type);
-extern ReflectionList *reflectionlist_sort_lfom(ReflectionList *in);
 
 #include "basis.h"
 extern ReflectionList *reflection_list_from_cell(Basis *basis);
diff --git a/src/reproject.c b/src/reproject.c
index 15a3ca9..3604481 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -17,14 +17,12 @@
 #include "utils.h"
 #include "imagedisplay.h"
 #include "displaywindow.h"
+#include "image.h"
 
-#define MAX_IMAGE_REFLECTIONS 8*1024
+ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) {
 
-ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionList *reflectionlist, ControlContext *ctx) {
-
-	ImageReflection *refl;
+	ImageFeatureList *flist;
 	Reflection *reflection;
-	int i;
 	double smax = 0.1e9;
 	double tilt, omega;
 	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
@@ -32,8 +30,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
 	//int done_debug = 0;
 			
-	tilt = deg2rad(image.tilt);
-	omega = deg2rad(image.omega);
+	flist = image_feature_list_new();
+	
+	tilt = deg2rad(image->tilt);
+	omega = deg2rad(image->omega);
 	
 	/* Calculate the (normalised) incident electron wavevector
 	 *	n is rotated "into" the reconstruction, so only one omega step. */
@@ -41,9 +41,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
 	nxt = nx;  nyt = ny;  nzt = nz;
 	nx = nxt*cos(-omega) + nyt*sin(-omega);  ny = -nxt*sin(-omega) + nyt*cos(-omega);  nz = nzt;
-	kx = nx / image.lambda;
-	ky = ny / image.lambda;
-	kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
+	kx = nx / image->lambda;
+	ky = ny / image->lambda;
+	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
 	//reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
 	
 	/* Determine where "up" is
@@ -54,9 +54,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	ux = uxt*cos(-omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
 	//reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
 	
-	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
 	reflection = reflectionlist->reflections;
-	i = 0;
 	do {
 	
 		double xl, yl, zl;
@@ -71,11 +69,11 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 		/* Next, solve the relrod equation to calculate the excitation error */
 		a = 1.0;
 		b = -2.0*(xl*nx + yl*ny + zl*nz);
-		c = xl*xl + yl*yl + zl*zl - 1.0/(image.lambda*image.lambda);
+		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
 		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
 		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
-		s1 = 1.0/image.lambda - A1;
-		s2 = 1.0/image.lambda - A2;
+		s1 = 1.0/image->lambda - A1;
+		s2 = 1.0/image->lambda - A2;
 		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
 		
 		/* Skip this reflection if s is large */
@@ -140,39 +138,31 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			//	if ( (i==20) && !done_debug ) done_debug = 1;
 				
 				/* Calculate image coordinates from polar representation */
-				if ( image.fmode == FORMULATION_CLEN ) {
-					x = image.camera_len*tan(theta)*cos(psi);
-					y = image.camera_len*tan(theta)*sin(psi);
-					x *= image.resolution;
-					y *= image.resolution;
-				} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
-					x = tan(theta)*cos(psi) / image.lambda;
-					y = tan(theta)*sin(psi) / image.lambda;
-					x /= image.pixel_size;
-					y /= image.pixel_size;
+				if ( image->fmode == FORMULATION_CLEN ) {
+					x = image->camera_len*tan(theta)*cos(psi);
+					y = image->camera_len*tan(theta)*sin(psi);
+					x *= image->resolution;
+					y *= image->resolution;
+				} else if ( image->fmode == FORMULATION_PIXELSIZE ) {
+					x = tan(theta)*cos(psi) / image->lambda;
+					y = tan(theta)*sin(psi) / image->lambda;
+					x /= image->pixel_size;
+					y /= image->pixel_size;
 				} else {
 					fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
 					return NULL;
 				}
 				
 				/* Adjust centre */
-				x += image.x_centre;
-				y += image.y_centre;
+				x += image->x_centre;
+				y += image->y_centre;
 				
 				/* Sanity check */
-				if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) {
+				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
 				
 					/* Record the reflection */
-					refl[i].x = x;
-					refl[i].y = y;
-					i++;
-					
-					if ( i > MAX_IMAGE_REFLECTIONS ) {
-						fprintf(stderr, "Too many reflections\n");
-						break;
-					}
-					
-					//printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
+					image_add_feature(flist, x, y, image, reflection->intensity);
+					//printf("Reflection at %f, %f\n", x, y);
 				
 				} /* else it's outside the picture somewhere */
 				
@@ -184,30 +174,28 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 
 	} while ( reflection );
 	
-	//printf("Found %i reflections in image\n", i);
-	*n = i;
-	return refl;
+	return flist;
 
 }
 
 static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {
 
-	ImageReflection *refl;
-	size_t n, j;
+	ImageFeatureList *flist;
+	size_t j;
 	
 	ctx->reproject_cur_image++;
-	if ( ctx->reproject_cur_image == ctx->n_images ) ctx->reproject_cur_image = 0;
+	if ( ctx->reproject_cur_image == ctx->images->n_images ) ctx->reproject_cur_image = 0;
 	
 	imagedisplay_clear_circles(ctx->reproject_id);
 	reflectionlist_clear_markers(ctx->reflectionlist);
 	
-	refl = reproject_get_reflections(ctx->images[ctx->reproject_cur_image], &n, ctx->cell_lattice, ctx);
-	for ( j=0; j<n; j++ ) {
-		imagedisplay_mark_circle(ctx->reproject_id, refl[j].x, refl[j].y);
+	flist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	for ( j=0; j<flist->n_features; j++ ) {
+		imagedisplay_mark_circle(ctx->reproject_id, flist->features[j].x, flist->features[j].y);
 	}
+	image_feature_list_free(flist);
 	
-	imagedisplay_put_data(ctx->reproject_id, ctx->images[ctx->reproject_cur_image]);
-	displaywindow_update(ctx->dw);
+	imagedisplay_put_data(ctx->reproject_id, ctx->images->images[ctx->reproject_cur_image]);
 	
 	return 0;
 
@@ -215,8 +203,8 @@ static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlC
 
 void reproject_open(ControlContext *ctx) {
 	
-	size_t n, j;
-	ImageReflection *refl;
+	size_t j;
+	ImageFeatureList *flist;
 	
 	if ( !ctx->cell ) {
 		printf("RP: No current cell\n");
@@ -233,13 +221,15 @@ void reproject_open(ControlContext *ctx) {
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
 	ctx->reproject_cur_image = 0;
 	
-	ctx->reproject_id = imagedisplay_open_with_message(ctx->images[ctx->reproject_cur_image], "Current Image", "Click to change image",
+	ctx->reproject_id = imagedisplay_open_with_message(ctx->images->images[ctx->reproject_cur_image], "Current Image", "Click to change image",
 				IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS, G_CALLBACK(reproject_clicked), ctx);
 	
-	refl = reproject_get_reflections(ctx->images[ctx->reproject_cur_image], &n, ctx->cell_lattice, ctx);
-	for ( j=0; j<n; j++ ) {
-		imagedisplay_mark_circle(ctx->reproject_id, refl[j].x, refl[j].y);
+	flist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	for ( j=0; j<flist->n_features; j++ ) {
+		imagedisplay_mark_circle(ctx->reproject_id, flist->features[j].x, flist->features[j].y);
 	}
+	printf("%i reflections\n", flist->n_features);
+	image_feature_list_free(flist);
 	
 }
 
diff --git a/src/reproject.h b/src/reproject.h
index 9ff3579..92e7183 100644
--- a/src/reproject.h
+++ b/src/reproject.h
@@ -17,9 +17,9 @@
 #define REPROJECT_H
 
 #include "control.h"
-#include "reflections.h"
+#include "image.h"
 
-extern ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionList *reflectionlist, ControlContext *ctx);
+extern ImageFeatureList *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionList *reflectionlist, ControlContext *ctx);
 extern void reproject_open(ControlContext *ctx);
 
 #endif	/* REPROJECT_H */