Beginnings of a new initial basis choice (needs fixing)

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

18 files changed, 250 insertions, 248 deletions

diff --git a/src/control.c b/src/control.c
index ac3de0b..97ad34c 100644
--- a/src/control.c
+++ b/src/control.c
@@ -25,6 +25,7 @@ int control_add_image(ControlContext *ctx, uint16_t *image, int width, int heigh
 	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;
@@ -56,3 +57,49 @@ ControlContext *control_ctx_new() {
 
 }
 
+/* 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 negative) 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 ( ctx->images[i].tilt - tilt < dev ) {
+			dev = ctx->images[i].tilt - tilt;
+			im = &ctx->images[i];
+		}
+	}
+	
+	return im;
+
+}
+
diff --git a/src/control.h b/src/control.h
index cf3b8c4..50c7230 100644
--- a/src/control.h
+++ b/src/control.h
@@ -51,8 +51,9 @@ typedef enum {
 typedef struct imagerecord_struct {
 
 	uint16_t	*image;
-	double		tilt;		/* Degrees */
-	double		omega;		/* Degrees */
+	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;
@@ -68,55 +69,56 @@ typedef struct imagerecord_struct {
 } ImageRecord;
 
 typedef struct {
-	int	x;
-	int	y;
-	double	intensity;
+	ImageRecord	*parent;
+	int		x;
+	int		y;
+	double		intensity;
 } ImageReflection;
 
 typedef struct cctx_struct {
 	
 	/* Modes */
-	InputFileType		inputfiletype;
-	FormulationMode		fmode;
-	ReconstructionMode 	rmode;
-	PeakSearchMode 		psmode;
-	unsigned int		prealign;
-	unsigned int		savecache;
-	unsigned int		have_centres;
+	InputFileType			inputfiletype;
+	FormulationMode			fmode;
+	ReconstructionMode 		rmode;
+	PeakSearchMode 			psmode;
+	unsigned int			prealign;
+	unsigned int			savecache;
+	unsigned int			have_centres;
 	
 	/* Input filename */
-	char 			*filename;
+	char 				*filename;
 	
 	/* Basic parameters, stored here solely so they can be copied
 	 * into the ImageRecord(s) more easily */
-	double			camera_length;
-	double			omega;		/* Degrees */
-	double			resolution;
-	double			lambda;
-	double			pixel_size;
-	double			x_centre;
-	double			y_centre;
+	double				camera_length;
+	double				omega;		/* Degrees */
+	double				resolution;
+	double				lambda;
+	double				pixel_size;
+	double				x_centre;
+	double				y_centre;
 	
 	/* QDRP Parser flags */
-	unsigned int		started;
-	unsigned int		camera_length_set;
-	unsigned int		omega_set;
-	unsigned int		resolution_set;
-	unsigned int		lambda_set;
+	unsigned int			started;
+	unsigned int			camera_length_set;
+	unsigned int			omega_set;
+	unsigned int			resolution_set;
+	unsigned int			lambda_set;
 	
 	/* The input images */
-	int			n_images;
-	ImageRecord		images[MAX_IMAGES];
+	int				n_images;
+	ImageRecord			images[MAX_IMAGES];
 	
 	/* Output */
-	struct rctx_struct	*reflectionctx;
-	struct dw_struct	*dw;
+	struct rctx_struct		*reflectionctx;
+	struct dw_struct		*dw;
 	
 	/* GTK bits */
-	GtkWidget		*combo_algorithm;
-	GtkWidget		*combo_peaksearch;
-	GtkWidget		*checkbox_prealign;
-	GtkWidget		*checkbox_savecache;
+	GtkWidget			*combo_algorithm;
+	GtkWidget			*combo_peaksearch;
+	GtkWidget			*checkbox_prealign;
+	GtkWidget			*checkbox_savecache;
 	
 	/* IPR stuff */
 	int				ipr_cur_image;
@@ -128,6 +130,9 @@ typedef struct cctx_struct {
 
 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 4ab7d23..6a3adfe 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -265,7 +265,7 @@ static gint displaywindow_gl_motion_notify(GtkWidget *widget, GdkEventMotion *ev
 	
 #define DRAW_BLOB							\
 	double step = M_PI/(double)BLOB_BITS;				\
-	double size = 0.2;						\
+	double size = 0.15;						\
 	int is, js;							\
 	for ( is=0; is<BLOB_BITS; is++ ) {				\
 		for ( js=0; js<BLOB_BITS*2; js++ ) {			\
diff --git a/src/ipr.c b/src/ipr.c
index 71444af..03824de 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -27,18 +27,6 @@
 #include "ipr.h"
 #include "displaywindow.h"
 
-static int ipr_choose_max(Basis *basis) {
-
-	if ( (basis->a.modulus >= basis->b.modulus) && (basis->a.modulus >= basis->c.modulus) ) return 1;
-	if ( (basis->b.modulus >= basis->a.modulus) && (basis->b.modulus >= basis->c.modulus) ) return 2;
-	if ( (basis->c.modulus >= basis->a.modulus) && (basis->c.modulus >= basis->b.modulus) ) return 3;
-	
-	/* Never reaches here... */
-	fprintf(stderr, "IP: Couldn't determine maximum basis vector modulus\n");
-	abort();
-	
-}
-
 #if 0
 static int ipr_random_image(ControlContext *ctx) {
 	double n;
@@ -48,160 +36,47 @@ static int ipr_random_image(ControlContext *ctx) {
 }
 #endif
 
-static Reflection *ipr_random_reflection(ReflectionContext *rctx) {
-	
-	double n;
-	int i;
-	Reflection *cur;
-	
-	n = (double)random()/RAND_MAX;
-	n *= rctx->n_reflections/2;
-	cur = rctx->reflections;	/* Jump a maximum of half way in */
-	for ( i=1; i<n; i++ ) {		/* Skip initial placeholder */
-		cur = cur->next;
-	}
-	
-	return cur;
-	
-}
-
-/* Try to turn the basis into a more sensible one */
-static void ipr_try_compact(Basis *basis) {
-
-	int i;
-	
-	for ( i=1; i<5; i++ ) {
-	
-		if ( modulus(basis->a.x - basis->b.x, basis->a.y - basis->b.y, basis->a.z - basis->b.z) < basis->a.modulus ) {
-			basis->a.x = basis->a.x - basis->b.x; basis->a.y = basis->a.y - basis->b.y; basis->a.y = basis->a.y - basis->b.y;
-			basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
-		}
-		if ( modulus(basis->a.x - basis->c.x, basis->a.z - basis->c.z, basis->a.z - basis->c.z) < basis->a.modulus ) {
-			basis->a.x = basis->a.x - basis->c.x; basis->a.y = basis->a.y - basis->c.y; basis->a.y = basis->a.y - basis->c.y;
-			basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
-		}
-		if ( modulus(basis->a.x + basis->b.x, basis->a.y + basis->b.y, basis->a.z + basis->b.z) < basis->a.modulus ) {
-			basis->a.x = basis->a.x + basis->b.x; basis->a.y = basis->a.y + basis->b.y; basis->a.y = basis->a.y + basis->b.y;
-			basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
-		}
-		if ( modulus(basis->a.x + basis->c.x, basis->a.z + basis->c.z, basis->a.z + basis->c.z) < basis->a.modulus ) {
-			basis->a.x = basis->a.x + basis->c.x; basis->a.y = basis->a.y + basis->c.y; basis->a.y = basis->a.y + basis->c.y;
-			basis->a.modulus = modulus(basis->a.x, basis->a.y, basis->a.z);
-		}
-
-		if ( modulus(basis->b.x - basis->a.x, basis->b.y - basis->a.y, basis->b.z - basis->a.z) < basis->b.modulus ) {
-			basis->b.x = basis->b.x - basis->a.x; basis->b.y = basis->b.y - basis->a.y; basis->b.y = basis->b.y - basis->a.y;
-			basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
-		}
-		if ( modulus(basis->b.x - basis->c.x, basis->b.y - basis->c.y, basis->b.z - basis->c.z) < basis->b.modulus ) {
-			basis->b.x = basis->b.x - basis->c.x; basis->b.y = basis->b.y - basis->c.y; basis->b.y = basis->b.y - basis->c.y;
-			basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
-		}
-		if ( modulus(basis->b.x + basis->a.x, basis->b.y + basis->a.y, basis->b.z + basis->a.z) < basis->b.modulus ) {
-			basis->b.x = basis->b.x + basis->a.x; basis->b.y = basis->b.y + basis->a.y; basis->b.y = basis->b.y + basis->a.y;
-			basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
-		}
-		if ( modulus(basis->b.x + basis->c.x, basis->b.y + basis->c.y, basis->b.z + basis->c.z) < basis->b.modulus ) {
-			basis->b.x = basis->b.x + basis->c.x; basis->b.y = basis->b.y + basis->c.y; basis->b.y = basis->b.y + basis->c.y;
-			basis->b.modulus = modulus(basis->b.x, basis->b.y, basis->b.z);
-		}
-
-		if ( modulus(basis->c.x - basis->a.x, basis->c.y - basis->a.z, basis->c.z - basis->a.z) < basis->c.modulus ) {
-			basis->c.x = basis->c.x - basis->a.x; basis->c.y = basis->c.y - basis->a.y; basis->c.y = basis->c.y - basis->a.y;
-			basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
-		}
-		if ( modulus(basis->c.x - basis->b.x, basis->c.y - basis->b.y, basis->c.z - basis->b.z) < basis->c.modulus ) {
-			basis->c.x = basis->c.x - basis->b.x; basis->c.y = basis->c.y - basis->b.y; basis->c.y = basis->c.y - basis->b.y;
-			basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
-		}
-		if ( modulus(basis->c.x + basis->a.x, basis->c.y + basis->a.z, basis->c.z + basis->a.z) < basis->c.modulus ) {
-			basis->c.x = basis->c.x + basis->a.x; basis->c.y = basis->c.y + basis->a.y; basis->c.y = basis->c.y + basis->a.y;
-			basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
-		}
-		if ( modulus(basis->c.x + basis->b.x, basis->c.y + basis->b.y, basis->c.z + basis->b.z) < basis->c.modulus ) {
-			basis->c.x = basis->c.x + basis->b.x; basis->c.y = basis->c.y + basis->b.y; basis->c.y = basis->c.y + basis->b.y;
-			basis->c.modulus = modulus(basis->c.x, basis->c.y, basis->c.z);
-		}
-
-	}
-
-}
-
 static Basis *ipr_choose_random_basis(ControlContext *ctx) {
-	
-	Basis *basis;
-	Reflection *reflection;
-	int cur_max;
+
+	Basis		*basis;
+	double		tilt_min;
+	double		tilt_max;
+	double		tilt_mid;
+	ImageRecord	*imagerecord;
+	double		nx, ny;
+	double		scale;
+	ImageReflection refl;
+	double		x, y, z, twotheta;
+	Reflection	*centre;
+	
+	/* Locate the 'plane' in the middle of the "wedge" */
+	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_min);
+	
+	/* Find the point in the middle of the "wedge" */
+	scale = reflection_largest_g(ctx->reflectionctx);
+	nx = imagerecord->x_centre - scale*sin(imagerecord->omega);
+	ny = imagerecord->y_centre + scale*cos(imagerecord->omega);
+	
+	/* Find where this point is in 3D space */
+	refl.x = nx;
+	refl.y = ny;
+	refl.parent = imagerecord;
+	if ( reflection_map_to_space(&refl, &x, &y, &z, &twotheta) ) return NULL;
+	
+	/* Find the nearest reflection */
+	centre = reflection_find_nearest(ctx->reflectionctx, x, y, z);
+	if ( !centre ) return NULL;
+	reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_MARKER);
+	printf("Nearest %f,%f,%f is %f,%f,%f\n", x, y, z, centre->x, centre->y, centre->z);
 	
 	basis = malloc(sizeof(Basis));
 	
-restart:
-	basis->a.x = 0;	basis->a.y = 0;	basis->a.z = 0;	basis->a.modulus = 1e30;
-	basis->b.x = 0;	basis->b.y = 0;	basis->b.z = 0;	basis->b.modulus = 1e30;
-	basis->c.x = 0;	basis->c.y = 0;	basis->c.z = 0;	basis->c.modulus = 1e30;
-	cur_max = 1;
-	
-	reflection = ipr_random_reflection(ctx->reflectionctx);
-	
-	do {
-		double mod;
-		unsigned int angle_works;
-		
-		mod = modulus(reflection->x, reflection->y, reflection->z);
-		
-		angle_works = 1;
-		if ( basis->a.modulus != 1e30 ) {
-			double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->a.x, basis->a.y, basis->a.z);
-			if ( angle < 20 ) angle_works = 0;
-			if ( angle > 160 ) angle_works = 0;
-		}
-		if ( basis->b.modulus != 1e30 ) {
-			double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->b.x, basis->b.y, basis->b.z);
-			if ( angle < 20 ) angle_works = 0;
-			if ( angle > 160 ) angle_works = 0;
-		}
-		if ( basis->c.modulus != 1e30 ) {
-			double angle = angle_between_d(reflection->x, reflection->y, reflection->z, basis->c.x, basis->c.y, basis->c.z);
-			if ( angle < 20 ) angle_works = 0;
-			if ( angle > 160 ) angle_works = 0;
-		}
-		if ( (mod > 1e9) && angle_works ) {
-			cur_max = ipr_choose_max(basis);
-			if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 1) ) {
-				basis->a.x = reflection->x;
-				basis->a.y = reflection->y;
-				basis->a.z = reflection->z;
-				basis->a.modulus = mod;
-			}
-			if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 2) ) {
-				basis->b.x = reflection->x;
-				basis->b.y = reflection->y;
-				basis->b.z = reflection->z;
-				basis->b.modulus = mod;
-			}
-			if ( (mod < basis->a.modulus) && (mod < basis->b.modulus) && (mod < basis->c.modulus) && (cur_max == 3) ) {
-				basis->c.x = reflection->x;
-				basis->c.y = reflection->y;
-				basis->c.z = reflection->z;
-				basis->c.modulus = mod;
-			}
-		}
-		reflection = reflection->next;
-		
-		/* End of list? */
-		if ( !reflection ) {
-			if ( (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) ) {
-				goto restart;	/* Probably the one circumstance under which this is allowed... */
-			}
-		}
-		
-	} while ( reflection || (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) );
-	
-	ipr_try_compact(basis);
-	
-//	basis->a.x = 1.842e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
-//	basis->b.x = 0.0;  basis->b.y = 1.842e9;  basis->b.z = 0.0;
-//	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.842e9;
+	basis->a.x = 1.842e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
+	basis->b.x = 0.0;  basis->b.y = 1.842e9;  basis->b.z = 0.0;
+	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.842e9;
 	
 	return basis;
 
@@ -335,6 +210,9 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 					if ( ref ) {	/* Check it's not NULL */
 						ref->h = h;	ref->k = k;	ref->l = l;
 					}
+					if ( (h!=0) || (k!=0) || (l!=0) ) {
+					//	reflection_add(ctx->reflectionctx, x, y, z, 1, REFLECTION_GENERATED);
+					}
 				}
 				
 			}
diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c
index baac6b6..cdc2a5c 100644
--- a/src/itrans-lsq.c
+++ b/src/itrans-lsq.c
@@ -129,16 +129,16 @@ static void itrans_interpolate(uint16_t *image, int width, int x, int y) {
 
 }
 
-unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+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;
+	width = imagerecord->width;
+	height = imagerecord->height;
+	image = imagerecord->image;
 
 	/* Least-Squares Craziness. NB Doesn't quite work... */
 	do {
@@ -262,7 +262,7 @@ unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx,
 				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);
+				reflection_add_from_dp(ctx, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, brightness);
 				n_reflections++;
 				
 				/* Remove this peak from the image */
diff --git a/src/itrans-lsq.h b/src/itrans-lsq.h
index 44f80d6..24f4e98 100644
--- a/src/itrans-lsq.h
+++ b/src/itrans-lsq.h
@@ -18,7 +18,7 @@
 
 #include "control.h"
 
-extern unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx);
+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 7099c8d..319e565 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -488,23 +488,23 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c
 
 }
 
-unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx) {
+unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ctx) {
 
 	unsigned int n_reflections;	
 	gsl_matrix *m;
 	gsl_matrix *p;
 	int i;
 	double px,py;
-	uint16_t *image = imagerecord.image;
+	uint16_t *image = imagerecord->image;
 
-	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord.width, imagerecord.height);
+	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord->width, imagerecord->height);
 	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
 	
 	for ( i=0; i<n_reflections; i++ ) {
 	
 		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);
+		reflection_add_from_dp(ctx, (px-imagerecord->x_centre), (py-imagerecord->y_centre), imagerecord, 1.0);
 		
 	}
 	gsl_matrix_free(m);
diff --git a/src/itrans-stat.h b/src/itrans-stat.h
index 3e6b5a5..f475a2c 100644
--- a/src/itrans-stat.h
+++ b/src/itrans-stat.h
@@ -19,6 +19,6 @@
 
 #include "control.h"
 
-unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx);
+unsigned int itrans_peaksearch_stat(ImageRecord *imagerecord, ControlContext *ctx);
 
 #endif /* ITRANS_STAT_H */
diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c
index 519acc8..30e476c 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -20,16 +20,16 @@
 #include "imagedisplay.h"
 #include "reflections.h"
 
-unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext *ctx) {
+unsigned int 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 *image = imagerecord->image;
 	uint16_t max = 0;
 	
-	width = imagerecord.width;
-	height = imagerecord.height;
+	width = imagerecord->width;
+	height = imagerecord->height;
 	
 	/* Simple Thresholding */
 	for ( y=0; y<height; y++ ) {
@@ -45,7 +45,7 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 				assert(y<height);
 				assert(x>=0);
 				assert(y>=0);
-				reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre),
+				reflection_add_from_dp(ctx, (x-imagerecord->x_centre), (y-imagerecord->y_centre),
 							imagerecord, image[x + width*y]);
 				n_reflections++;
 			}
@@ -56,7 +56,7 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 	
 }
 
-unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, ControlContext *ctx) {
+unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, ControlContext *ctx) {
 
 	uint16_t max_val = 0;
 	int width, height;
@@ -65,9 +65,9 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr
 	uint16_t max;
 	int x, y;
 	
-	image = imagerecord.image;
-	width = imagerecord.width;
-	height = imagerecord.height;
+	image = imagerecord->image;
+	width = imagerecord->width;
+	height = imagerecord->height;
 	
 	max = 0;
 	for ( y=0; y<height; y++ ) {
@@ -100,7 +100,7 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr
 			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),
+			reflection_add_from_dp(ctx, (max_x-imagerecord->x_centre), (max_y-imagerecord->y_centre),
 							imagerecord, image[x + width*y]);
 			n_reflections++;
 			
diff --git a/src/itrans-threshold.h b/src/itrans-threshold.h
index 24b3cc1..f19ea57 100644
--- a/src/itrans-threshold.h
+++ b/src/itrans-threshold.h
@@ -18,8 +18,8 @@
 
 #include "control.h"
 
-extern unsigned int itrans_peaksearch_threshold(ImageRecord image, ControlContext *ctx);
-extern unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord image, ControlContext *ctx);
+extern unsigned int itrans_peaksearch_threshold(ImageRecord *image, ControlContext *ctx);
+extern unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord *image, ControlContext *ctx);
 
 #endif	/* ITRANS_THRESHOLD_H */
 
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index 490f644..889aac4 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -23,16 +23,16 @@
 
 #define PEAK_WINDOW_SIZE 20
 
-unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext *ctx) {
+unsigned int itrans_peaksearch_zaefferer(ImageRecord *imagerecord, ControlContext *ctx) {
 
 	int x, y;
 	unsigned int n_reflections;
 	int width, height;
 	uint16_t *image;
 	
-	image = imagerecord.image;
-	width = imagerecord.width;
-	height = imagerecord.height;
+	image = imagerecord->image;
+	width = imagerecord->width;
+	height = imagerecord->height;
 	
 	n_reflections = 0;
 	for ( x=1; x<width-1; x++ ) {
@@ -85,7 +85,7 @@ unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext
 					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),
+					reflection_add_from_dp(ctx, (mask_x-imagerecord->x_centre), (mask_y-imagerecord->y_centre),
 							imagerecord, image[mask_x + width*mask_y]);
 					n_reflections++;
 				}
diff --git a/src/itrans-zaefferer.h b/src/itrans-zaefferer.h
index 77ae704..7a3e7ed 100644
--- a/src/itrans-zaefferer.h
+++ b/src/itrans-zaefferer.h
@@ -18,7 +18,7 @@
 
 #include "control.h"
 
-extern unsigned int itrans_peaksearch_zaefferer(ImageRecord image, ControlContext *ctx);
+extern unsigned int itrans_peaksearch_zaefferer(ImageRecord *image, ControlContext *ctx);
 
 #endif	/* ITRANS_ZAEFFERER_H */
 
diff --git a/src/itrans.c b/src/itrans.c
index db89736..05d2793 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -21,7 +21,7 @@
 #include "itrans-lsq.h"
 #include "itrans-stat.h"
 
-void itrans_process_image(ImageRecord image, ControlContext *ctx) {
+void itrans_process_image(ImageRecord *image, ControlContext *ctx) {
 	
 	unsigned int n_reflections;
 	
diff --git a/src/itrans.h b/src/itrans.h
index 3db2b5c..a293a98 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 void itrans_process_image(ImageRecord *image, ControlContext *ctx);
 
 #endif	/* ITRANS_H */
 
diff --git a/src/mapping.c b/src/mapping.c
index 251c804..cac1ce5 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -30,7 +30,7 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	printf("MP: Processing images..."); fflush(stdout);
 	ctx->reflectionctx = rctx;
 	for ( i=0; i<ctx->n_images; i++ ) {
-		itrans_process_image(ctx->images[i], ctx);
+		itrans_process_image(&ctx->images[i], ctx);
 	}
 	printf("done\n");
 	
diff --git a/src/reflections.c b/src/reflections.c
index 1a80f71..63eaa54 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -13,11 +13,14 @@
 #include <config.h>
 #endif
 
+#define _GNU_SOURCE
 #include <stdlib.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 
 #include "reflections.h"
+#include "utils.h"
 
 static void reflection_addfirst(ReflectionContext *reflectionctx) {
 	/* Create first items on lists - saves faffing later.  Corresponds to a central marker.
@@ -124,11 +127,11 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 	
 }
 
-/* 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) {
+int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
 
 	/* "Input" space */
-	double d;
+	double d, x, y;
+	ImageRecord *imagerecord;
 	
 	/* Angular description of reflection */
 	double theta, psi, k;
@@ -140,24 +143,26 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 	double x_temp, y_temp, z_temp;
 	double nx, ny, nz;
 	
-	tilt = 2*M_PI*(imagerecord.tilt/360);	/* Convert to Radians */
-	omega = 2*M_PI*(imagerecord.omega/360);	/* Likewise */
-	k = 1/imagerecord.lambda;
+	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 ( ctx->fmode == FORMULATION_CLEN ) {
-		x /= imagerecord.resolution;
-		y /= imagerecord.resolution;	/* Convert pixels to metres */
+	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 ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-		x *= imagerecord.pixel_size;
-		y *= imagerecord.pixel_size;	/* Convert pixels to metres^-1 */
+		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;
+		return -1;
 	}
 	psi = atan2(y, x);
 	
@@ -177,14 +182,36 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 	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;
 	
-	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
+	*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->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
+	}
 
 }
 
@@ -194,3 +221,46 @@ void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) {
 	rctx->last_reflection = r;
 }
 
+double reflection_largest_g(ReflectionContext *reflectionctx) {
+
+	double max = 0.0;
+	Reflection *reflection;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( reflection->type == REFLECTION_NORMAL ) {
+			double mod;
+			mod = modulus(reflection->x, reflection->y, reflection->z);
+			if ( mod > max ) max = mod;
+		}
+		reflection = reflection->next;
+	};
+	
+	return max;
+
+}
+
+Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z) {
+
+	double max = +INFINITY;
+	Reflection *reflection;
+	Reflection *best = NULL;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( reflection->type == REFLECTION_NORMAL ) {
+			double mod;
+			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+			if ( mod < max ) {
+				max = mod;
+				best = reflection;
+			}
+		}
+		reflection = reflection->next;
+	};
+	
+	return best;
+
+
+}
+
diff --git a/src/reflections.h b/src/reflections.h
index 2dd363c..824fcba 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -63,8 +63,11 @@ extern Reflection *reflection_add(ReflectionContext *reflectionctx, double x, do
 
 #include "control.h"
 
-extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity);
+extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord *imagerecord, double intensity);
 extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r);
+extern double reflection_largest_g(ReflectionContext *reflectionctx);
+extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta);
+extern Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z);
 
 #endif	/* REFLECTION_H */
 
diff --git a/src/reproject.c b/src/reproject.c
index 7c3ec70..0d79900 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -77,7 +77,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
 		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 */
@@ -103,8 +102,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			
 				double psi, disc;
 				
-				reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_GENERATED);
-				reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
+				//reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_GENERATED);
+				//reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
 				
 				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
 				cx = yi*nz-zi*ny;  cy = nx*zi-nz*xi;  cz = ny*xi-nx*yi;  /* c = i x n */