Be a lot more clear about which things are properties of images and which

belong to the control context.


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

19 files changed, 143 insertions, 148 deletions

diff --git a/src/control.h b/src/control.h
index 68eebed..5de0b86 100644
--- a/src/control.h
+++ b/src/control.h
@@ -80,34 +80,27 @@ typedef struct cctx_struct {
 	ReconstructionMode 	rmode;
 	PeakSearchMode 		psmode;
 	unsigned int		prealign;
+	unsigned int		savecache;
 	
 	/* Input filename */
 	char 			*filename;
 	
-	/* Basic parameters */
+	/* Basic parameters, stored here solely so they can be copied
+	 * into the ImageRecord(s) more easily */
 	double			camera_length;
 	double			omega;
 	double			resolution;
 	double			lambda;
 	double			pixel_size;
-
-	/* (QDRP) Parser flags */
+	
+	/* QDRP Parser flags */
 	unsigned int		started;
 	unsigned int		camera_length_set;
 	unsigned int		omega_set;
 	unsigned int		resolution_set;
 	unsigned int		lambda_set;
-
-	/* Miscellaneous modes and operations */
-	unsigned int		first_image;
-	
-	/* Size of input images - assumed the same throughout. */
-	int			width;
-	int			height;
-	int			x_centre;
-	int			y_centre;
 	
-	/* Information about the input images */
+	/* The input images */
 	int			n_images;
 	ImageRecord		images[MAX_IMAGES];
 	
@@ -118,6 +111,7 @@ typedef struct cctx_struct {
 	GtkWidget		*combo_algorithm;
 	GtkWidget		*combo_peaksearch;
 	GtkWidget		*checkbox_prealign;
+	GtkWidget		*checkbox_savecache;
 
 } ControlContext;
 
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 9c7aee5..3a87aa5 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -76,29 +76,34 @@ static void imagedisplay_put_data(ImageDisplay *imagedisplay, int16_t *image16)
 
 	unsigned int x, y;
 	unsigned int w, h;
+	int16_t min, max;
 	
 	h = imagedisplay->height;
 	w = imagedisplay->width;
-
+	
+	min = 0;  max = 0;
+	for ( y=0; y<h; y++ ) {
+		for ( x=0; x<w; x++ ) {
+			int16_t val;
+			val = image16[x+w*y];
+			if ( val > max ) max = val;
+			if ( val < min ) min = val;
+		}
+	}
+	
 	/* Turn 16-bit image data into 8-bit display data */
 	imagedisplay->data = malloc(3*w*h);
 	for ( y=0; y<h; y++ ) {
 		for ( x=0; x<w; x++ ) {
 		
-			int16_t val16 = image16[x+w*y];
-			if ( val16 > 255 ) {
-				imagedisplay->data[3*( x+w*(h-1-y) )] = 255;
-				imagedisplay->data[3*( x+w*(h-1-y) )+1] = 255;
-				imagedisplay->data[3*( x+w*(h-1-y) )+2] = 255;
-			} else if ( val16 < 0 ) {
-				imagedisplay->data[3*( x+w*(h-1-y) )] = 0;
-				imagedisplay->data[3*( x+w*(h-1-y) )+1] = 0;
-				imagedisplay->data[3*( x+w*(h-1-y) )+2] = 0;
-			} else {
-				imagedisplay->data[3*( x+w*(h-1-y) )] = val16;
-				imagedisplay->data[3*( x+w*(h-1-y) )+1] = val16;
-				imagedisplay->data[3*( x+w*(h-1-y) )+2] = val16;
-			}
+			int16_t val16, val8;
+			
+			val16 = image16[x+w*y];
+			val8 = (255*(val16-min)) / (max-min);
+			
+			imagedisplay->data[3*( x+w*(h-1-y) )] = val8;
+			imagedisplay->data[3*( x+w*(h-1-y) )+1] = val8;
+			imagedisplay->data[3*( x+w*(h-1-y) )+2] = val8;
 			
 		}
 	}
@@ -149,7 +154,7 @@ ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsigned int
 	
 }
 
-void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx, double omega) {
+void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, double xc, double yc, double omega) {
 
 	guchar *image_eightbit;
 	int w, h;
@@ -170,7 +175,7 @@ void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx,
 		gradient = 1/gradient;
 		/* Start at the centre and draw a line out in each direction until it hits an edge.
 			This makes the whole thing a lot easier. */
-		xs = ctx->x_centre;  y = ctx->y_centre;
+		xs = xc;  y = yc;
 		do {
 			x = xs;
 			image_eightbit[3*(x+w*(h-1-y))+0] = 255;
@@ -179,7 +184,7 @@ void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx,
 			y++;
 			xs += gradient;
 		} while ( (xs<w) && (y<h) && (xs>=0) && (y>=0) );
-		xs = ctx->x_centre;  y = ctx->y_centre;
+		xs = xc;  y = yc;
 		do {
 			x = xs;
 			image_eightbit[3*(x+w*(h-1-y))+0] = 255;
@@ -191,7 +196,7 @@ void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx,
 	} else {
 		double ys;
 		signed int x, y;
-		x = ctx->x_centre;  ys = ctx->y_centre;
+		x = xc;  ys = yc;
 		do {
 			y = ys;
 			image_eightbit[3*(x+w*(h-1-y))+0] = 255;
@@ -201,7 +206,7 @@ void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx,
 			ys += gradient;
 		} while ( (x<w) && (ys<h) && (x>=0) && (ys>=0) );
 		
-		x = ctx->x_centre;  ys = ctx->y_centre;
+		x = xc;  ys = yc;
 		do {
 			y = ys;
 			image_eightbit[3*(x+w*(h-1-y))+0] = 255;
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index 84ffb03..61aa157 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -41,7 +41,7 @@ extern ImageDisplay *imagedisplay_open(int16_t *image, unsigned int width, unsig
 extern void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x, unsigned int y);
 extern void imagedisplay_mark_circle(ImageDisplay *imagedisplay, unsigned int x, unsigned int y);
 #include "control.h"
-extern void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx, double omega);
+extern void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, double xc, double yc, double omega);
 
 #endif	/* IMAGEDISPLAY_H */
 
diff --git a/src/ipr.c b/src/ipr.c
index 0296765..c8f9d8a 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -191,10 +191,13 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 	
 	/* NB This assumes the diffraction pattern is "vaguely" centered... */
 	if ( ctx->inputfiletype != INPUT_CACHE) {
+		int max_width, max_height;
+		max_width = ctx->images[0].width;
+		max_height = ctx->images[0].height;
 		if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-			max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) * ctx->pixel_size;
+			max_res = sqrt(max_width*max_width + max_height*max_height) * ctx->images[0].pixel_size;
 		} else { 
-			max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) / (ctx->lambda * ctx->camera_length);
+			max_res = sqrt(max_width*max_width + max_height*max_height) / (ctx->images[0].lambda * ctx->images[0].camera_len);
 		}
 		max_res = max_res / 2;
 	} else {
diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c
index 1bf3d80..841a8b9 100644
--- a/src/itrans-lsq.c
+++ b/src/itrans-lsq.c
@@ -129,13 +129,16 @@ static void itrans_interpolate(int16_t *image, int width, int x, int y) {
 
 }
 
-unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
 
 	int16_t max_val = 0;
 	int width, height;
+	int16_t *image;
 	unsigned int n_reflections = 0;
 	
-	width = ctx->width;	height = ctx->height;
+	width = imagerecord.width;
+	height = imagerecord.height;
+	image = imagerecord.image;
 
 	/* Least-Squares Craziness. NB Doesn't quite work... */
 	do {
@@ -256,14 +259,11 @@ unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, double t
 				double frac;
 			
 				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);
-				if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-					reflection_add_from_reciprocal(ctx, (x-ctx->x_centre)*ctx->pixel_size, (y-ctx->y_centre)*ctx->pixel_size,
+				if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
+					reflection_add_from_reciprocal(ctx, (x-imagerecord.x_centre)*imagerecord.pixel_size, (y-imagerecord.y_centre)*imagerecord.pixel_size,
 										tilt_degrees, image[x + width*y]);
 				} else {
-					reflection_add_from_dp(ctx, (x-ctx->x_centre), (y-ctx->y_centre), tilt_degrees, image[x + width*y]);
-				}
-				if ( ctx->first_image ) {
-					imagedisplay_mark_point(imagedisplay, x, y);
+					reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), tilt_degrees, image[x + width*y]);
 				}
 				n_reflections++;
 				
@@ -333,3 +333,4 @@ unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, double t
 	
 	return n_reflections;
 }
+
diff --git a/src/itrans-lsq.h b/src/itrans-lsq.h
index 949a58b..44f80d6 100644
--- a/src/itrans-lsq.h
+++ b/src/itrans-lsq.h
@@ -16,11 +16,9 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-
 #include "control.h"
-#include "imagedisplay.h"
 
-extern unsigned int itrans_peaksearch_lsq(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
+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 c2369ec..6d64a08 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -29,15 +29,15 @@
  * can be done as matrices to save effort
  * Renormalises matrix to 0->1
  */
-static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(ControlContext *ctx, int16_t *image) {
+static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(int16_t *image, int width, int height) {
 
 	gsl_matrix *raw;
 	int i, j;
 
-	raw = gsl_matrix_alloc(ctx->width,ctx->height);
+	raw = gsl_matrix_alloc(width, height);
 	for ( i=0; i<raw->size1; i++ ) {
 		for ( j=0; j<raw->size2; j++ ) {
-			gsl_matrix_set(raw, i, j, (double)image[i+j*ctx->width]);
+			gsl_matrix_set(raw, i, j, (double)image[i+j*width]);
 		}
 	}
 	matrix_renormalise(raw);
@@ -488,28 +488,29 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c
 
 }
 
-unsigned int itrans_peaksearch_stat(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx) {
 
 	unsigned int n_reflections;	
 	gsl_matrix *m;
 	gsl_matrix *p;
 	int i;
 	double px,py;
+	double tilt_degrees = imagerecord.tilt;
+	int16_t *image = imagerecord.image;
 
-	m = itrans_peaksearch_stat_createImageMatrix(ctx, image);
+	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);
-		if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-			reflection_add_from_reciprocal(ctx, (px-ctx->x_centre)*ctx->pixel_size, (py-ctx->y_centre)*ctx->pixel_size,
+		if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
+			reflection_add_from_reciprocal(ctx, (px-imagerecord.x_centre)*imagerecord.pixel_size, (py-imagerecord.y_centre)*imagerecord.pixel_size,
 								tilt_degrees, 1.0);
 		} else {
-			reflection_add_from_dp(ctx, (px-ctx->x_centre), (py-ctx->y_centre), tilt_degrees, 1.0);
+			reflection_add_from_dp(ctx, (px-imagerecord.x_centre)/imagerecord.resolution, (py-imagerecord.y_centre)/imagerecord.resolution, tilt_degrees, 1.0);
 		}
-		if (ctx->first_image) imagedisplay_mark_point(imagedisplay, (unsigned int)px, (unsigned int)py);
 		
 	}
 	gsl_matrix_free(m);
diff --git a/src/itrans-stat.h b/src/itrans-stat.h
index ad2f8e3..3e6b5a5 100644
--- a/src/itrans-stat.h
+++ b/src/itrans-stat.h
@@ -17,11 +17,8 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-
 #include "control.h"
-#include "imagedisplay.h"
 
-unsigned int itrans_peaksearch_stat(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
+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 3a1d788..6aadea2 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -19,27 +19,27 @@
 #include "imagedisplay.h"
 #include "reflections.h"
 
-unsigned int itrans_peaksearch_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext *ctx) {
 
 	int width, height;
 	int x, y;
 	unsigned int n_reflections = 0;
+	int16_t *image = imagerecord.image;
+	double tilt_degrees = imagerecord.tilt;
 	
-	width = ctx->width;	height = ctx->height;
+	width = imagerecord.width;
+	height = imagerecord.height;
 
 	/* Simple Thresholding */
 	for ( y=0; y<height; y++ ) {
 		for ( x=0; x<width; x++ ) {
 
 			if ( image[x + width*y] > 10 ) {
-				if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-					reflection_add_from_reciprocal(ctx, (signed)(x-ctx->x_centre)*ctx->pixel_size, (signed)(y-ctx->y_centre)*ctx->pixel_size,
+				if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
+					reflection_add_from_reciprocal(ctx, (signed)(x-imagerecord.x_centre)*imagerecord.pixel_size, (signed)(y-imagerecord.y_centre)*imagerecord.pixel_size,
 									tilt_degrees, image[x + width*y]);
 				} else {
-					reflection_add_from_dp(ctx, (signed)(x-ctx->x_centre), (signed)(y-ctx->y_centre), tilt_degrees, image[x + width*y]);
-				}
-				if ( ctx->first_image ) {
-					imagedisplay_mark_point(imagedisplay, x, y);
+					reflection_add_from_dp(ctx, (signed)(x-imagerecord.x_centre)/imagerecord.resolution, (signed)(y-imagerecord.y_centre)/imagerecord.resolution, tilt_degrees, image[x + width*y]);
 				}
 				n_reflections++;
 			}
@@ -51,14 +51,17 @@ unsigned int itrans_peaksearch_threshold(int16_t *image, ControlContext *ctx, do
 	
 }
 
-unsigned int itrans_peaksearch_adaptive_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, ControlContext *ctx) {
 
 	int16_t max_val = 0;
 	int width, height;
 	unsigned int n_reflections = 0;
+	int16_t *image = imagerecord.image;
+	double tilt_degrees = imagerecord.tilt;
+	
+	width = imagerecord.width;
+	height = imagerecord.height;
 	
-	width = ctx->width;	height = ctx->height;
-
 	/* Adaptive Thresholding */
 	do {
 	
@@ -80,14 +83,11 @@ unsigned int itrans_peaksearch_adaptive_threshold(int16_t *image, ControlContext
 		}
 		
 		if ( max_val > 50 ) {
-			if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-				reflection_add_from_reciprocal(ctx, (signed)(max_x-ctx->x_centre)*ctx->pixel_size, (signed)(max_y-ctx->y_centre)*ctx->pixel_size,
-								tilt_degrees, max_val);
+			if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
+				reflection_add_from_reciprocal(ctx, (signed)(x-imagerecord.x_centre)*imagerecord.pixel_size, (signed)(y-imagerecord.y_centre)*imagerecord.pixel_size,
+									tilt_degrees, image[x + width*y]);
 			} else {
-				reflection_add_from_dp(ctx, (signed)(max_x-ctx->x_centre), (signed)(max_y-ctx->y_centre), tilt_degrees, max_val);
-			}
-			if ( ctx->first_image ) {
-				imagedisplay_mark_point(imagedisplay, max_x, max_y);
+				reflection_add_from_dp(ctx, (signed)(x-imagerecord.x_centre)/imagerecord.resolution, (signed)(y-imagerecord.y_centre)/imagerecord.resolution, tilt_degrees, image[x + width*y]);
 			}
 			n_reflections++;
 			
@@ -104,3 +104,4 @@ unsigned int itrans_peaksearch_adaptive_threshold(int16_t *image, ControlContext
 	return n_reflections;
 
 }
+
diff --git a/src/itrans-threshold.h b/src/itrans-threshold.h
index 32c3faf..24b3cc1 100644
--- a/src/itrans-threshold.h
+++ b/src/itrans-threshold.h
@@ -16,12 +16,10 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-
 #include "control.h"
-#include "imagedisplay.h"
 
-extern unsigned int itrans_peaksearch_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
-extern unsigned int itrans_peaksearch_adaptive_threshold(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
+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 aaa8110..a096e74 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -22,24 +22,32 @@
 
 #define PEAK_WINDOW_SIZE 20
 
-unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext *ctx) {
 
 	int x, y;
 	unsigned int n_reflections;
+	int width, height;
+	int16_t *image;
+	double tilt_degrees;
+	
+	tilt_degrees = imagerecord.tilt;
+	image = imagerecord.image;
+	width = imagerecord.width;
+	height = imagerecord.height;
 	
 	n_reflections = 0;
-	for ( x=1; x<ctx->width-1; x++ ) {
-		for ( y=1; y<ctx->height-1; y++ ) {
+	for ( x=1; x<width-1; x++ ) {
+		for ( y=1; y<height-1; y++ ) {
 
 			double dx1, dx2, dy1, dy2;
 			double dxs, dys;
 			double grad;
 			
 			/* Get gradients */
-			dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y];
-			dx2 = image[(x-1)+ctx->width*y] - image[x+ctx->width*y];
-			dy1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*(y+1)];
-			dy2 = image[x+ctx->width*(y-1)] - image[x+ctx->width*y];
+			dx1 = image[x+width*y] - image[(x+1)+width*y];
+			dx2 = image[(x-1)+width*y] - image[x+width*y];
+			dy1 = image[x+width*y] - image[(x+1)+width*(y+1)];
+			dy2 = image[x+width*(y-1)] - image[x+width*y];
 			
 			/* Average gradient measurements from both sides */
 			dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
@@ -59,12 +67,12 @@ unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, do
 				mask_y = y;
 				
 				while ( (did_something) && (distance(mask_x, mask_y, x, y)<50) ) {
-					max = image[mask_x+ctx->width*mask_y];
+					max = image[mask_x+width*mask_y];
 					did_something = 0;
-					for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, ctx->height); sy++ ) {
-						for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, ctx->width); sx++ ) {
-							if ( image[sx+ctx->width*sy] > max ) {
-								max = image[sx+ctx->width*sy];
+					for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, height); sy++ ) {
+						for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, width); sx++ ) {
+							if ( image[sx+width*sy] > max ) {
+								max = image[sx+width*sy];
 								mask_x = sx;
 								mask_y = sy;
 								did_something = 1;
@@ -74,16 +82,13 @@ unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, do
 				}
 				
 				if ( !did_something ) {
-					if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-						reflection_add_from_reciprocal(ctx, (signed)(mask_x-ctx->x_centre)*ctx->pixel_size,
-										(signed)(mask_y-ctx->y_centre)*ctx->pixel_size,
-										tilt_degrees, image[mask_x + ctx->width*mask_y]);
+					if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
+						reflection_add_from_reciprocal(ctx, (signed)(mask_x-imagerecord.x_centre)*imagerecord.pixel_size,
+										(signed)(mask_y-imagerecord.y_centre)*imagerecord.pixel_size,
+										tilt_degrees, image[mask_x + width*mask_y]);
 					} else {
-						reflection_add_from_dp(ctx, (signed)(mask_x-ctx->x_centre), (signed)(mask_y-ctx->y_centre),
-										tilt_degrees, image[mask_x + ctx->width*mask_y]);
-					}
-					if ( ctx->first_image ) {
-						imagedisplay_mark_point(imagedisplay, mask_x, mask_y);
+						reflection_add_from_dp(ctx, (signed)(mask_x-imagerecord.x_centre)/imagerecord.resolution, (signed)(mask_y-imagerecord.y_centre)/imagerecord.resolution,
+										tilt_degrees, image[mask_x + width*mask_y]);
 					}
 					n_reflections++;
 				}
@@ -95,3 +100,4 @@ unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, do
 	return n_reflections;
 
 }
+
diff --git a/src/itrans-zaefferer.h b/src/itrans-zaefferer.h
index 2726b1b..77ae704 100644
--- a/src/itrans-zaefferer.h
+++ b/src/itrans-zaefferer.h
@@ -16,11 +16,9 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-
 #include "control.h"
-#include "imagedisplay.h"
 
-extern unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
+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 dbb164f..db89736 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -14,38 +14,25 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-#include <gsl/gsl_matrix.h>
-
 #include "control.h"
-#include "imagedisplay.h"
 #include "reflections.h"
 #include "itrans-threshold.h"
 #include "itrans-zaefferer.h"
 #include "itrans-lsq.h"
 #include "itrans-stat.h"
 
-void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees) {
+void itrans_process_image(ImageRecord image, ControlContext *ctx) {
 	
 	unsigned int n_reflections;
-	ImageDisplay *imagedisplay = NULL;
-	
-	ctx->first_image = 0;
-	if ( ctx->first_image ) {
-		imagedisplay = imagedisplay_open(image, ctx->width, ctx->height, "Image Display");
-		imagedisplay_add_tilt_axis(imagedisplay, ctx, ctx->omega);
-	}
 	
 	switch ( ctx->psmode ) {
-		case PEAKSEARCH_THRESHOLD : n_reflections = itrans_peaksearch_threshold(image, ctx, tilt_degrees, imagedisplay); break;
-		case PEAKSEARCH_ADAPTIVE_THRESHOLD : n_reflections = itrans_peaksearch_adaptive_threshold(image, ctx, tilt_degrees, imagedisplay); break;
-		case PEAKSEARCH_LSQ : n_reflections = itrans_peaksearch_lsq(image, ctx, tilt_degrees, imagedisplay); break;
-		case PEAKSEARCH_ZAEFFERER : n_reflections = itrans_peaksearch_zaefferer(image, ctx, tilt_degrees, imagedisplay); break;
-		case PEAKSEARCH_STAT : n_reflections = itrans_peaksearch_stat(image, ctx, tilt_degrees, imagedisplay); break;
+		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;
 	}
 	
-	ctx->first_image = 0;
-	
 }
 
diff --git a/src/itrans.h b/src/itrans.h
index 558ba8a..3db2b5c 100644
--- a/src/itrans.h
+++ b/src/itrans.h
@@ -16,11 +16,9 @@
 #include <config.h>
 #endif
 
-#include <stdint.h>
-
 #include "control.h"
 
-extern void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees);
+extern void itrans_process_image(ImageRecord image, ControlContext *ctx);
 
 #endif	/* ITRANS_H */
 
diff --git a/src/main.c b/src/main.c
index 5df4bcf..125909b 100644
--- a/src/main.c
+++ b/src/main.c
@@ -60,6 +60,12 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			ctx->prealign = FALSE;
 		}
 		
+		if ( gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(ctx->checkbox_savecache)) ) {
+			ctx->savecache = TRUE;
+		} else {
+			ctx->savecache = FALSE;
+		}
+		
 		gtk_widget_destroy(method_window);
 		while ( gtk_events_pending() ) gtk_main_iteration();
 		
@@ -81,7 +87,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			mapping_create(ctx);
 		}
 		
-		if ( (ctx->inputfiletype != INPUT_CACHE) && !val && (ctx->reflectionctx) ) {
+		if ( (ctx->inputfiletype != INPUT_CACHE) && !val && ctx->reflectionctx && ctx->savecache ) {
 			cache_save(ctx->filename, ctx->reflectionctx);			
 		}
 		
@@ -149,6 +155,9 @@ void main_method_dialog_open(ControlContext *ctx) {
 	ctx->checkbox_prealign = gtk_check_button_new_with_label("Pre-align image stack");
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_prealign, 1, 3, 3, 4);
 	
+	ctx->checkbox_savecache = gtk_check_button_new_with_label("Save 3D mapping cache file");
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_savecache, 1, 3, 4, 5);
+	
 	if ( ctx->inputfiletype == INPUT_CACHE ) {
 		gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get from cache file");
 		gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch), FALSE);
diff --git a/src/mapping.c b/src/mapping.c
index d9ca94d..24a6470 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -27,11 +27,14 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	int twidth, theight;
 	int x, y, i;
 	int16_t *image_total;
+	double x_centre = 0;
+	double y_centre = 0;
 	
 	/* Create reflection context */
 	rctx = reflection_init();
 	
 	if ( !ctx->prealign ) {
+	
 		/* Find centre of image stack
 		 * Determine maximum size of image to accommodate, and allocate memory */
 		twidth = 0;  theight = 0;
@@ -64,13 +67,13 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 				}
 			}
 		}
-		ctx->x_centre = max_x;
-		ctx->y_centre = max_y;
-	
+		x_centre = max_x;
+		y_centre = max_y;
+		
 		/* Display */
 		sum_id = imagedisplay_open(image_total, twidth, theight, "Sum of All Images");
-		imagedisplay_mark_point(sum_id, ctx->x_centre, ctx->y_centre);
-		imagedisplay_add_tilt_axis(sum_id, ctx, ctx->omega);
+		imagedisplay_mark_point(sum_id, x_centre, y_centre);
+		imagedisplay_add_tilt_axis(sum_id, x_centre, y_centre, ctx->omega);
 		
 	}
 	
@@ -79,10 +82,10 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	ctx->reflectionctx = rctx;
 	for ( i=0; i<ctx->n_images; i++ ) {
 		if ( !ctx->prealign ) {
-			ctx->images[i].x_centre = ctx->x_centre;
-			ctx->images[i].y_centre = ctx->y_centre;
+			ctx->images[i].x_centre = x_centre;
+			ctx->images[i].y_centre = y_centre;
 		}
-		itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt);
+		itrans_process_image(ctx->images[i], ctx);
 	}
 	
 	return rctx;
diff --git a/src/mrc.c b/src/mrc.c
index abc096c..50596a5 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -64,9 +64,6 @@ int mrc_read(ControlContext *ctx) {
 	pixel_size = ext[0].pixel_size;
 	printf("pixel_size = %f m^-1\n", pixel_size);
 	ctx->fmode = FORMULATION_PIXELSIZE;
-	ctx->first_image = 1;
-	ctx->width = mrc.nx;
-	ctx->height = mrc.ny;
 	
 	for ( i=0; i<mrc.nz; i++ ) {
 	
diff --git a/src/qdrp.c b/src/qdrp.c
index fc484ee..936f2b6 100644
--- a/src/qdrp.c
+++ b/src/qdrp.c
@@ -199,7 +199,6 @@ int qdrp_read(ControlContext *ctx) {
 	ctx->resolution_set = 0;
 	ctx->lambda_set = 0;
 	ctx->started = 0;
-	ctx->first_image = 1;
 	
 	FILE *fh;
 	
diff --git a/src/readpng.c b/src/readpng.c
index 14e5fe2..dc72e6e 100644
--- a/src/readpng.c
+++ b/src/readpng.c
@@ -97,7 +97,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	height = png_get_image_height(png_ptr, info_ptr);
 	bit_depth = png_get_bit_depth(png_ptr, info_ptr);
 	channels = png_get_channels(png_ptr, info_ptr);
-//	printf("RI: width=%i, height=%i, depth=%i, channels=%i\n", width, height, bit_depth, channels);
+	printf("RI: width=%i, height=%i, depth=%i, channels=%i\n", width, height, bit_depth, channels);
 	if ( (bit_depth != 16) && (bit_depth != 8) ) {
 		printf("RI: Whoops! Can't handle images with other than 8 or 16 bpp yet...\n");
 		png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
@@ -109,7 +109,6 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	row_pointers = png_get_rows(png_ptr, info_ptr);
 	
 	image = malloc(height * width * sizeof(int16_t));
-	ctx->width = width; ctx->height = height;
 	ctx->fmode = FORMULATION_CLEN;
 	
 	for ( y=0; y<height; y++ ) {
@@ -121,8 +120,9 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 				int i;
 				val = 0;
 				for ( i=0; i<channels; i++ ) {
-					val += row_pointers[y][(channels*x)+i] << 8;
-					val += row_pointers[y][(channels*x)+i];
+					/* PNG files are big-endian: */
+					val += row_pointers[y][(channels*x*2)+(2*i)] << 8;
+					val += row_pointers[y][(channels*x*2)+(2*i)+1];
 				}
 			}
 			if ( bit_depth == 8 ) {