Preparation for handing precessed data

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

26 files changed, 300 insertions, 684 deletions

diff --git a/Makefile.am b/Makefile.am
index c90b55e..62f85be 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,4 +1,4 @@
 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/peakdetect.h src/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h src/structure.h src/itrans-threshold.h \
-		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h
+		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h
 SUBDIRS = src data
diff --git a/Makefile.in b/Makefile.in
index 7347fa6..eda9bf4 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -155,7 +155,7 @@ sysconfdir = @sysconfdir@
 target_alias = @target_alias@
 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/peakdetect.h src/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h src/structure.h src/itrans-threshold.h \
-		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h
+		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h
 
 SUBDIRS = src data
 all: config.h
diff --git a/src/Makefile.am b/src/Makefile.am
index c7174e4..65c8218 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,6 +1,6 @@
 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 ipr.c cache.c \
-		structure.c itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c
+		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
diff --git a/src/Makefile.in b/src/Makefile.in
index d2c3512..dc8fe2f 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -51,9 +51,9 @@ am_dtr_OBJECTS = main.$(OBJEXT) displaywindow.$(OBJEXT) \
 	trackball.$(OBJEXT) reflections.$(OBJEXT) readpng.$(OBJEXT) \
 	mrc.$(OBJEXT) imagedisplay.$(OBJEXT) utils.$(OBJEXT) \
 	itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT) cache.$(OBJEXT) \
-	structure.$(OBJEXT) itrans-threshold.$(OBJEXT) \
-	itrans-zaefferer.$(OBJEXT) itrans-lsq.$(OBJEXT) \
-	itrans-stat.$(OBJEXT)
+	itrans-threshold.$(OBJEXT) itrans-zaefferer.$(OBJEXT) \
+	itrans-lsq.$(OBJEXT) itrans-stat.$(OBJEXT) control.$(OBJEXT) \
+	mapping.$(OBJEXT)
 dtr_OBJECTS = $(am_dtr_OBJECTS)
 dtr_DEPENDENCIES =
 DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)
@@ -153,7 +153,7 @@ sharedstatedir = @sharedstatedir@
 sysconfdir = @sysconfdir@
 target_alias = @target_alias@
 dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c cache.c \
-		structure.c itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c
+		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c
 
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
@@ -225,6 +225,7 @@ distclean-compile:
 	-rm -f *.tab.c
 
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cache.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/control.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/imagedisplay.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ipr.Po@am__quote@
@@ -234,11 +235,11 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/itrans-zaefferer.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/itrans.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/main.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mapping.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mrc.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/qdrp.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/readpng.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/reflections.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/structure.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/trackball.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/utils.Po@am__quote@
 
diff --git a/src/control.c b/src/control.c
new file mode 100644
index 0000000..b626d7c
--- /dev/null
+++ b/src/control.c
@@ -0,0 +1,47 @@
+/*
+ * control.c
+ *
+ * Common control structure
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#include <inttypes.h>
+
+#include "control.h"
+
+int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size) {
+
+	ctx->images[ctx->n_images].tilt = tilt;
+	ctx->images[ctx->n_images].omega = omega;
+	ctx->images[ctx->n_images].pixel_size = pixel_size;
+	ctx->images[ctx->n_images].camera_len = 0;
+	ctx->images[ctx->n_images].image = image;
+	ctx->images[ctx->n_images].width = width;
+	ctx->images[ctx->n_images].height = height;
+	
+	ctx->n_images++;
+	
+	return ctx->n_images - 1;
+
+}
+
+int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len) {
+
+	ctx->images[ctx->n_images].tilt = tilt;
+	ctx->images[ctx->n_images].omega = omega;
+	ctx->images[ctx->n_images].pixel_size = 0;
+	ctx->images[ctx->n_images].camera_len = camera_len;
+	ctx->images[ctx->n_images].image = image;
+	ctx->images[ctx->n_images].width = width;
+	ctx->images[ctx->n_images].height = height;
+	
+	ctx->n_images++;
+	
+	return ctx->n_images - 1;
+
+}
+
diff --git a/src/control.h b/src/control.h
index ed9e128..4566315 100644
--- a/src/control.h
+++ b/src/control.h
@@ -4,6 +4,7 @@
  * Common control structure
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -16,6 +17,9 @@
 #define CONTROL_H
 
 #include <gtk/gtk.h>
+#include <inttypes.h>
+
+#define MAX_IMAGES 256
 
 typedef enum ift_enum {
 	INPUT_NONE,
@@ -43,55 +47,69 @@ typedef enum {
 	RECONSTRUCTION_PREDICTION
 } ReconstructionMode;
 
+typedef struct imagerecord_struct {
+
+	int16_t		*image;
+	double		tilt;
+	double		omega;
+	double		pixel_size;
+	double		camera_len;
+	int		width;
+	int		height;
+
+} ImageRecord;
+
 typedef struct cctx_struct {
 	
 	/* Modes */
-	InputFileType inputfiletype;
-	FormulationMode fmode;
-	ReconstructionMode rmode;
-	PeakSearchMode psmode;
-
+	InputFileType		inputfiletype;
+	FormulationMode		fmode;
+	ReconstructionMode 	rmode;
+	PeakSearchMode 		psmode;
+	
 	/* Input filename */
-	char *filename;
+	char 			*filename;
 	
 	/* Basic parameters */
-	double camera_length;
-	double omega;
-	double resolution;
-	double lambda;
-	double pixel_size;
+	double			camera_length;
+	double			omega;
+	double			resolution;
+	double			lambda;
+	double			pixel_size;
 
 	/* (QDRP) Parser flags */
-	unsigned int started;
-	unsigned int camera_length_set;
-	unsigned int omega_set;
-	unsigned int max_d_set;
-	unsigned int resolution_set;
-	unsigned int lambda_set;
+	unsigned int		started;
+	unsigned int		camera_length_set;
+	unsigned int		omega_set;
+	unsigned int		max_d_set;
+	unsigned int		resolution_set;
+	unsigned int		lambda_set;
 
 	/* Miscellaneous modes and operations */
-	unsigned int max_d;
-	unsigned int first_image;
+	unsigned int		max_d;
+	unsigned int		first_image;
 	
 	/* Size of input images - assumed the same throughout. */
-	unsigned int width;
-	unsigned int height;
-	unsigned int x_centre;
-	unsigned int y_centre;
+	int			width;
+	int			height;
+	int			x_centre;
+	int			y_centre;
 	
 	/* Information about the input images */
-	unsigned int n_images;
-	
+	int			n_images;
+	ImageRecord		images[MAX_IMAGES];
 	
 	/* Output */
-	struct rctx_struct *reflectionctx;
+	struct rctx_struct	*reflectionctx;
 	
 	/* GTK bits */
-	GtkWidget *combo_algorithm;
-	GtkWidget *combo_peaksearch;
+	GtkWidget		*combo_algorithm;
+	GtkWidget		*combo_peaksearch;
 
 } ControlContext;
 
-extern ControlContext *control_read(const char *filename);
+extern int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size);
+extern int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len);
 
 #endif	/* CONTROL_H */
+
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 8f7d226..cc8fa0e 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -4,6 +4,7 @@
  * The display window
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -25,7 +26,6 @@
 #include "trackball.h"
 #include "reflections.h"
 #include "main.h"
-#include "structure.h"
 
 enum {
 	DW_ORTHO,
@@ -43,28 +43,6 @@ GLfloat		displaywindow_distance = 150;
 GLfloat		displaywindow_x_pos = 0;
 GLfloat		displaywindow_y_pos = 0;
 
-OctTree *o;
-OctTreeLinkedList *ol;
-
-void disp_oct_stack(OctTree *o) {
-	int i;
-	for (i=0;i<8;i++) {
-		if (o->child[i] != NULL) {
-			glVertex3f(o->ox/1e9,o->oy/1e9,o->oz/1e9);
-			glVertex3f(o->child[i]->ox/1e9,o->child[i]->oy/1e9,o->child[i]->oz/1e9);
-			disp_oct_stack(o->child[i]);
-		}	
-	}
-}
-
-void disp_oct_list(OctTreeLinkedList *ol) {
-	OctTree *o;
-	do {
-		o=ol->o;
-		glVertex3f(o->ox/1e9, o->oy/1e9, o->oz/1e9);
-	} while((ol=ol->next)!=NULL);
-}
-
 void error_report(const char *message) {
  
 	GtkWidget *window;
@@ -379,6 +357,7 @@ static gint displaywindow_gl_expose(GtkWidget *widget, GdkEventExpose *event, Co
 		}
 		if ( reflection->type == REFLECTION_GENERATED ) {
 		
+		#if 0
 			/* Generated reflections are displayed by index, not coordinates.
 				xyz field in Reflection will be used later for "measured" location. */
 			double x, y, z;
@@ -394,8 +373,9 @@ static gint displaywindow_gl_expose(GtkWidget *widget, GdkEventExpose *event, Co
 			z = h*ctx->reflectionctx->basis->a.z + k*ctx->reflectionctx->basis->b.z + l*ctx->reflectionctx->basis->c.z;
 			glTranslatef(x/1e9, y/1e9, z/1e9);
 			
-			gluSphere(quadric, 0.2, 32, 32);
+			gluSphere(quadric, 0.2, 8, 8);
 			glPopMatrix();
+		#endif
 			
 		}
 		reflection = reflection->next;
@@ -456,18 +436,6 @@ static gint displaywindow_gl_expose(GtkWidget *widget, GdkEventExpose *event, Co
 		glVertex3f(-50, 50, 0.0);
 	glEnd();
 	
-	//octtree bit - temporary
-	/*glBegin(GL_LINES);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue);
-		
-		disp_oct_stack(o);
-	glEnd();
-	glBegin(GL_POINTS);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
-		disp_oct_list(ol);
-	glEnd();*/
-		
-
 	if ( gdk_gl_drawable_is_double_buffered(gldrawable) ) {
 		gdk_gl_drawable_swap_buffers(gldrawable);
 	} else {
@@ -583,3 +551,4 @@ void displaywindow_open(ControlContext *ctx) {
 	gtk_widget_show_all(displaywindow_window);
 	
 }
+
diff --git a/src/displaywindow.h b/src/displaywindow.h
index f4a7f27..3d20d4b 100644
--- a/src/displaywindow.h
+++ b/src/displaywindow.h
@@ -4,6 +4,7 @@
  * The display window
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -25,3 +26,4 @@ extern void displaywindow_statusbar(const char *message);
 extern void error_report(const char *message);
 
 #endif	/* DISPLAYWINDOW_H */
+
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 0c46123..6ca646c 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -4,6 +4,7 @@
  * Show raw and processed images
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -242,3 +243,4 @@ void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x, unsigne
 	}
 
 }
+
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index 0aa140f..9d75ded 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -4,6 +4,7 @@
  * Show raw and processed images
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -42,3 +43,4 @@ extern void imagedisplay_mark_point(ImageDisplay *imagedisplay, unsigned int x,
 extern void imagedisplay_add_tilt_axis(ImageDisplay *imagedisplay, ControlContext *ctx, double omega);
 
 #endif	/* IMAGEDISPLAY_H */
+
diff --git a/src/ipr.c b/src/ipr.c
index aadc75c..31de6bf 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -4,6 +4,7 @@
  * Iterative Prediction-Refinement Reconstruction
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -20,7 +21,7 @@
 #include "reflections.h"
 #include "itrans.h"
 #include "utils.h"
-#include "structure.h"
+#include "imagedisplay.h"
 
 static int ipr_choose_max(Basis *basis) {
 
@@ -99,8 +100,7 @@ static void ipr_try_compact(Basis *basis) {
 /* Choose an initial set of three vectors for the expansion.
  * This would probably be just as easy for the user to do... */
 static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
-
-	printf("ipr: c_i_b: starting\n");
+	
 	Basis *basis;
 	Reflection *reflection;
 	int cur_max;
@@ -113,11 +113,11 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 	basis->c.x = 0;	basis->c.y = 0;	basis->c.z = 0;	basis->c.modulus = 1e30;
 	reflection = ctx->reflectionctx->reflections->next;	/* First reflection, skip initial placeholder */
 	cur_max = 1;
-	printf("ipr: c_i_b: starting loop\n");
+	
 	do {
 		double mod;
 		unsigned int angle_works;
-		//printf("ipr: c_i_b: iteration - reflection (%f,%f,%f)\n",reflection->x,reflection->y,reflection->z);
+		
 		mod = modulus(reflection->x, reflection->y, reflection->z);
 		
 		angle_works = 1;
@@ -175,11 +175,6 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 	printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus);
 	printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus);
 	
-	reflection_clear(ctx->reflectionctx);
-	reflection_add(ctx->reflectionctx, basis->a.x, basis->a.y, basis->a.z, 1, REFLECTION_GENERATED);
-	reflection_add(ctx->reflectionctx, basis->b.x, basis->b.y, basis->b.z, 1, REFLECTION_GENERATED);
-	reflection_add(ctx->reflectionctx, basis->c.x, basis->c.y, basis->c.z, 1, REFLECTION_GENERATED);
-	
 	ipr_try_compact(basis);
 	
 	return basis;
@@ -226,49 +221,41 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 
 }
 
-/* Sort the list of reflections into order of increasing g */
-static void ipr_sort(ReflectionContext *rctx) {
-
-	Reflection *prev;
-	Reflection *this;
-	Reflection *next;
-	
-	this = rctx->reflections;
-	next = this->next;
-	prev = NULL;
-	
-	while ( next ) {
-	
-		next = next->next;	/* ! */
-		prev = prev->next;
-		this = this->next;
-		
-	}
-
+static int ipr_random_image(ControlContext *ctx) {
+	double n;
+	n = (double)random()/RAND_MAX;
+	n *= ctx->n_images;
+	return (int)n;
 }
 
-int ipr_reduce(ControlContext *ctx) {
-	printf("ipr: starting\n");
-	Basis *basis;
+int ipr_refine(ControlContext *ctx) {
 	
-	printf("ipr: calling choose_initial_basis\n");
+	Basis *basis;
+	int finished;
 	basis = ipr_choose_initial_basis(ctx);
 	ctx->reflectionctx->basis = basis;
 	if ( !basis ) return -1;
 	
-	
 	/* Get rid of the original list and replace it with the prediction list */
-	printf("ipr: calling reflection_clear\n");
-	reflection_clear(ctx->reflectionctx);
-	free(ctx->reflectionctx);
-	printf("ipr: calling ipr_generate\n");
-	ctx->reflectionctx = ipr_generate(ctx, basis);
+	//reflection_clear(ctx->reflectionctx);
+	//free(ctx->reflectionctx);
+	//ctx->reflectionctx = ipr_generate(ctx, basis);
 	ctx->reflectionctx->basis = basis;
 	
-	/* Sort the reflections into order of increasing g */
-	printf("ipr: calling ipr_sort\n");
-	ipr_sort(ctx->reflectionctx);
+	/* Select an image */
+	finished = 0;
+	do {
+	
+		int i;
+		ImageDisplay *cur;
+		
+		i = ipr_random_image(ctx);
+		cur = imagedisplay_open(ctx->images[i].image, ctx->images[i].width, ctx->images[i].height, "Current Image");
+		finished = 1;
+		
+	} while ( !finished );
 	
 	return 0;
 
 }
+
diff --git a/src/ipr.h b/src/ipr.h
index 3fcb40d..8fb4862 100644
--- a/src/ipr.h
+++ b/src/ipr.h
@@ -4,6 +4,7 @@
  * Iterative prediction-refinement reconstruction
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -18,6 +19,7 @@
 
 #include "control.h"
 
-extern int ipr_reduce(ControlContext *ctx);
+extern int ipr_refine(ControlContext *ctx);
 
 #endif	/* IPR_H */
+
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index b51c5cc..8ac845b 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -81,11 +81,11 @@ 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->width/2)*ctx->pixel_size,
-										(signed)(mask_y-ctx->height/2)*ctx->pixel_size,
+						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]);
 					} else {
-						reflection_add_from_dp(ctx, (signed)(mask_x-ctx->width/2), (signed)(mask_y-ctx->height/2),
+						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 ) {
diff --git a/src/itrans.c b/src/itrans.c
index d1e43ab..dbb164f 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -30,6 +30,7 @@ void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degre
 	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);
@@ -45,6 +46,6 @@ void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degre
 	}
 	
 	ctx->first_image = 0;
-	printf(" ..... %i peaks found\n", n_reflections);
 	
 }
+
diff --git a/src/itrans.h b/src/itrans.h
index a11eadf..558ba8a 100644
--- a/src/itrans.h
+++ b/src/itrans.h
@@ -23,3 +23,4 @@
 extern void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees);
 
 #endif	/* ITRANS_H */
+
diff --git a/src/main.c b/src/main.c
index 528d508..4b42809 100644
--- a/src/main.c
+++ b/src/main.c
@@ -30,7 +30,7 @@
 #include "qdrp.h"
 #include "ipr.h"
 #include "cache.h"
-#include "structure.h"
+#include "mapping.h"
 
 static gint main_method_window_response(GtkWidget *method_window, gint response, ControlContext *ctx) {
 
@@ -56,6 +56,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 		gtk_widget_destroy(method_window);
 		while ( gtk_events_pending() ) gtk_main_iteration();
 		
+		/* Load the input */
 		if ( ctx->inputfiletype == INPUT_QDRP ) {
 			val = qdrp_read(ctx);
 		} else if ( ctx->inputfiletype == INPUT_MRC ) {
@@ -65,16 +66,18 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			val=0;
 		}
 		
+		if ( ctx->inputfiletype != INPUT_CACHE ) {
+			mapping_create(ctx);
+		}
+		
 		if ( (ctx->inputfiletype != INPUT_CACHE) && !val && (ctx->reflectionctx) ) {
 			cache_save(ctx->filename, ctx->reflectionctx);			
 		}
 		
 		if ( !val && (ctx->rmode == RECONSTRUCTION_PREDICTION) && (ctx->reflectionctx) ) {
-			val = ipr_reduce(ctx);
+			val = ipr_refine(ctx);
 		}
 		
-		//dump_histogram(ctx->reflectionctx);
-		
 		if ( !val && (ctx->reflectionctx) ) {
 			displaywindow_open(ctx);
 		} else {
@@ -222,3 +225,4 @@ int main(int argc, char *argv[]) {
 	return 0;
 	
 }
+
diff --git a/src/mapping.c b/src/mapping.c
new file mode 100644
index 0000000..adf3205
--- /dev/null
+++ b/src/mapping.c
@@ -0,0 +1,84 @@
+/*
+ * mapping.c
+ *
+ * 3D Mapping
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "reflections.h"
+#include "control.h"
+#include "imagedisplay.h"
+#include "itrans.h"
+
+ReflectionContext *mapping_create(ControlContext *ctx) {
+
+	ReflectionContext *rctx;
+	ImageDisplay *sum_id;
+	int16_t max_val;
+	int max_x, max_y;
+	int twidth, theight;
+	int x, y, i;
+	int16_t *image_total;
+	
+	/* Create reflection context */
+	rctx = reflection_init();
+	
+	/* Find centre of image stack
+	 * Determine maximum size of image to accommodate, and allocate memory */
+	twidth = 0;  theight = 0;
+	for ( i=0; i<ctx->n_images; i++ ) {
+		if ( ctx->images[i].width > twidth ) twidth = ctx->images[i].width;
+		if ( ctx->images[i].height > theight ) theight = ctx->images[i].height;
+	}
+	image_total = malloc(twidth * theight * sizeof(int16_t));
+	memset(image_total, 0, twidth * theight * sizeof(int16_t));
+	
+	/* Add the image stack together */
+	for ( i=0; i<ctx->n_images; i++ ) {
+		int xoffs, yoffs;
+		xoffs = (twidth - ctx->images[i].width)/2;
+		yoffs = (theight - ctx->images[i].height)/2;
+		for ( y=0; y<theight; y++ ) {
+			for ( x=0; x<twidth; x++ ) {
+				image_total[x + twidth*y] += ctx->images[i].image[x + twidth*y]/ctx->n_images;
+			}
+		}
+	}
+	
+	/* Locate the highest point */
+	max_val = 0;	max_x = 0;	max_y = 0;
+	for ( y=0; y<theight; y++ ) {
+		for ( x=0; x<twidth; x++ ) {
+			if ( image_total[x + twidth*y] > max_val ) {
+				max_val = image_total[x + twidth*y];
+				max_x = x;  max_y = y;
+			}
+		}
+	}
+	ctx->x_centre = max_x;
+	ctx->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);
+	
+	/* Pass all images through itrans
+	 *   (let itrans add the reflections to rctx for now) */
+	ctx->reflectionctx = rctx;
+	for ( i=0; i<ctx->n_images; i++ ) {
+		itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt);
+	}
+	
+	return rctx;
+
+}
+
diff --git a/src/mapping.h b/src/mapping.h
new file mode 100644
index 0000000..f40d8f1
--- /dev/null
+++ b/src/mapping.h
@@ -0,0 +1,25 @@
+/*
+ * mapping.h
+ *
+ * 3D Mapping
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+ 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef MAPPING_H
+#define MAPPING_H
+
+#include "reflections.h"
+#include "control.h"
+
+extern ReflectionContext *mapping_create(ControlContext *ctx);
+
+#endif	/* MAPPING_H */
+
diff --git a/src/mrc.c b/src/mrc.c
index b9b6efc..26a31de 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -15,7 +15,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
-#include <stdint.h>
+#include <inttypes.h>
 #include <string.h>
 
 #include "mrc.h"
@@ -33,11 +33,6 @@ int mrc_read(ControlContext *ctx) {
 	unsigned int i;
 	unsigned int extsize;
 	double pixel_size;
-	int16_t *image_total;
-	unsigned int x, y;
-	int max_x, max_y;
-	int16_t max_val;
-	ImageDisplay *sum_id;
 	
 	fh = fopen(ctx->filename, "rb");
 
@@ -68,50 +63,16 @@ int mrc_read(ControlContext *ctx) {
 	
 	pixel_size = ext[0].pixel_size;
 	printf("pixel_size = %f m^-1\n", pixel_size);
-	ctx->reflectionctx = reflection_init();
 	ctx->fmode = FORMULATION_PIXELSIZE;
 	ctx->first_image = 1;
 	ctx->width = mrc.nx;
 	ctx->height = mrc.ny;
 	
-	printf("Judging centre...");
-	image_total = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
-	for ( y=0; y<mrc.ny; y++ ) {
-		for ( x=0; x<mrc.nx; x++ ) {
-			image_total[x + mrc.nx*y] = 0;
-		}
-	}
-	for ( i=0; i<mrc.nz; i++ ) {
-		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
-		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
-		fread(image, mrc.nx*mrc.ny*2, 1, fh);
-		for ( y=0; y<mrc.ny; y++ ) {
-			for ( x=0; x<mrc.nx; x++ ) {
-				image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz;
-			}
-		}
-	}
-	sum_id = imagedisplay_open(image_total, mrc.nx, mrc.ny, "Sum of All Images");
-	/* Locate the highest point */
-	max_val = 0;	max_x = 0;	max_y = 0;
-	for ( y=0; y<ctx->height; y++ ) {
-		for ( x=0; x<ctx->width; x++ ) {
-			if ( image_total[x + ctx->width*y] > max_val ) {
-				max_val = image_total[x + ctx->width*y];
-				max_x = x;  max_y = y;
-			}
-		}
-	}
-	imagedisplay_mark_point(sum_id, max_x, max_y);
-	ctx->x_centre = max_x;
-	ctx->y_centre = max_y;
-	printf("done\n");
-	
 	for ( i=0; i<mrc.nz; i++ ) {
 	
 		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
 		
-		printf("Image #%3i: tilt=%f omega=%f L=%f\t", i, ext[i].a_tilt, ext[i].tilt_axis, ext[i].magnification);
+		printf("Image #%3i: tilt=%f omega=%f L=%f\n", i, ext[i].a_tilt, ext[i].tilt_axis, ext[i].magnification);
 		ctx->camera_length = ext[i].magnification;
 		if ( ext[i].voltage == 0 ) {
 			ctx->lambda = lambda(200000);
@@ -123,15 +84,8 @@ int mrc_read(ControlContext *ctx) {
 		
 		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
 		fread(image, mrc.nx*mrc.ny*2, 1, fh);
-
-		for ( y=0; y<mrc.ny; y++ ) {
-			for ( x=0; x<mrc.nx; x++ ) {
-				image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz;
-			}
-		}
-
-		itrans_process_image(image, ctx, ext[i].a_tilt);
-		free(image);
+		
+		control_add_image_ps(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt, ctx->omega, ctx->pixel_size);
 		
 	}
 
diff --git a/src/mrc.h b/src/mrc.h
index ff5f0f7..16c7d9d 100644
--- a/src/mrc.h
+++ b/src/mrc.h
@@ -121,3 +121,4 @@ extern int mrc_read(ControlContext *ctx);
 extern unsigned int mrc_is_mrcfile(const char *filename);
 
 #endif	/* MRC_H */
+
diff --git a/src/qdrp.c b/src/qdrp.c
index 6579369..6cc9296 100644
--- a/src/qdrp.c
+++ b/src/qdrp.c
@@ -83,7 +83,7 @@ static int qdrp_parseline(ControlContext *ctx, const char *line) {
 	if ( ctx->started ) {
 		return qdrp_parsefileline(ctx, line);						
 	}
-
+	
 	if ( ( line[0] == '-' ) && ( line[1] == '-') ) {
 		
 		if ( !ctx->camera_length_set ) {
@@ -234,8 +234,9 @@ int qdrp_read(ControlContext *ctx) {
 	while ( !feof(fh) ) {			
 		fgets(line, 256, fh);			
 		if ( !feof(fh) ) {				
-			qdrp_chomp(line);				
-			if ( line[0] != '#' ) {					
+			qdrp_chomp(line);
+			if ( strlen(line) == 0 ) continue;			
+			if ( line[0] != '#' ) {
 				if ( qdrp_parseline(ctx, line) ) {
 					fclose(fh);
 					free(line);
@@ -262,3 +263,4 @@ unsigned int qdrp_is_qdrprc(const char *filename) {
 	return 0;
 
 }
+
diff --git a/src/qdrp.h b/src/qdrp.h
index 7987f6d..8d438aa 100644
--- a/src/qdrp.h
+++ b/src/qdrp.h
@@ -21,3 +21,4 @@ extern int qdrp_read(ControlContext *ctx);
 extern unsigned int qdrp_is_qdrprc(const char *filename);
 
 #endif	/* QDRP_H */
+
diff --git a/src/readpng.c b/src/readpng.c
index bb5a60b..68a92c3 100644
--- a/src/readpng.c
+++ b/src/readpng.c
@@ -1,9 +1,10 @@
 /*
  * readpng.c
  *
- * Read images
+ * Read PNG images
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -36,8 +37,6 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	png_bytep *row_pointers;
 	unsigned int x;
 	unsigned int y;
-	unsigned int xorig;
-	unsigned int yorig;
 	int16_t *image;
 	
 	/* Open file */
@@ -109,8 +108,6 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	/* Get image data */
 	row_pointers = png_get_rows(png_ptr, info_ptr);
 	
-	xorig = 320;
-	yorig = 320;	/* For now */
 	image = malloc(height * width * sizeof(int16_t));
 	ctx->width = width; ctx->height = height;
 	ctx->fmode = FORMULATION_CLEN;
@@ -120,27 +117,33 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 		
 			int val = 0;
 
-			if ( (x-xorig)*(x-xorig) + (y-yorig)*(y-yorig) < ctx->max_d*ctx->max_d ) {
-			
-				if ( bit_depth == 16 ) {
-					val = row_pointers[y][channels*x] << 8;
-					val += row_pointers[y][(channels*x)+1];
+			if ( bit_depth == 16 ) {
+				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];
 				}
-				if ( bit_depth == 8 ) {
-					val = row_pointers[y][channels*x];
+			}
+			if ( bit_depth == 8 ) {
+				int i;
+				val = 0;
+				for ( i=0; i<channels; i++ ) {
+					val += row_pointers[y][(channels*x)+i];
 				}
-				
-				image[x + width*y] = val;
-					
 			}
+			
+			image[x + width*y] = val;
+
 		}
 	}
 	
 	png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
 	fclose(fh);
 
-	itrans_process_image(image, ctx, tilt);
+	control_add_image_cl(ctx, image, width, height, tilt, ctx->omega, ctx->camera_length);
 	
 	return 0;
 	
 }
+
diff --git a/src/reflections.c b/src/reflections.c
index 4ecc268..c1b557e 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -189,3 +189,4 @@ void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) {
 	rctx->last_reflection->next = r;
 	rctx->last_reflection = r;
 }
+
diff --git a/src/structure.c b/src/structure.c
deleted file mode 100644
index 8ce0888..0000000
--- a/src/structure.c
+++ /dev/null
@@ -1,444 +0,0 @@
-/*
- * structure.c
- *
- * 3D analysis
- *
- * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
- *  dtr - Diffraction Tomography Reconstruction
- *
- */
- #include <string.h>
- #include <stdlib.h>
- #include <math.h>
- #include "reflections.h"
- #include "utils.h"
- #include "structure.h"
- 
- static int reflect_count(ReflectionContext *rctx) {
- 	Reflection *r;
- 	r = rctx->reflections;
- 	int count=0;
- 	do {
- 		count++;
- 	} while ((r=r->next)!=NULL);	
- 	return count;
- }
- 
- 
- /*
-  * Find the largest single dimension variable in the world
-  * Relevant for octtree
-  */
- static double largest_dimension(ReflectionContext *rctx) {
- 	Reflection *r;
- 	double max=0;
- 	double val;
- 	r = rctx->reflections;
- 	do {
- 		if ((val = r->x) > max) max = val;
- 		if ((val = r->y) > max) max = val;
- 		if ((val = r->z) > max) max = val;
- 	} while ((r = r->next) != NULL);
- 	return max;
- }
- 
- /*
-  * Calculate the volume of a shell from r1->r2
-  * Requires r2 > r1
-  */
- static double get_shell_volume(double r2, double r1) {
- 	return (4. * M_PI * (r2*r2*r2 - r1*r1*r1))/3.; 
- }
- 
- /*
-  * Attempts to calculate the optimum mask radius
-  * We have a problem here with non-isotropicness since the points in question usually are close to planar, and not equally distributed about the origin
-  * It may be more necessary to try and construct a bounding polygon - although that's going to get horribly messy
-  * Will attempt an octtree based method, but for the moment this will return the radius for which 2/3 of the points are contained.
-  */
- 
- double get_mask_radius(ReflectionContext *rctx) {
- 	Reflection *r;
- 	double maxmod = 0.;
- 	double modul;
- 	double density;
- 	double size;
- 	double mask=0.;
- 	int num_intervals=10;
- 	int interval;
- 	int count=0;
- 	r = rctx->reflections;
- 	do {
-		if (r->type == REFLECTION_NORMAL) { 		
- 			if ((modul = modulus(r->x,r->y,r->z)) > maxmod) maxmod = modul;
- 			count++;
-		}
- 	} while ((r = r->next) != NULL);
- 	
- 	printf("g_m_r: count=%d, maxmod=%f\n",count,maxmod);
- 	
- 	int *bucket;
- 	bucket = calloc(num_intervals,sizeof(int));
- 	maxmod *= 1.001;
- 	size = maxmod/num_intervals;
- 	r = rctx->reflections;
- 	do {
- 		if (r->type == REFLECTION_NORMAL) {
- 			modul = modulus(r->x,r->y,r->z);
- 			interval = (int)((modul/maxmod)* (double)num_intervals);
- 			bucket[interval] += 1;
- 		}
- 	} while ((r = r->next) != NULL);
- 	
- 	int i;
- 	int sum=0;
- 	for (i=0;i<num_intervals;i++) {
- 		density = ( bucket[i] * 1e30 )/get_shell_volume(size*(i+1),size*i);
- 		printf("interval %d - count=%d density=%f volume=%f r2=%f r1=%f\n",i,bucket[i],density,get_shell_volume(size*(i+1),size*i), size*(i+1), size*i);
- 		if (( sum += bucket[i]) < (0.66*count)) mask = size*(i+1);
- 	}
- 	printf("g_m_r: returning %f\n",mask);
- 	return mask;
- }
- 
- /*
-  * Returns a new ReflectionContext with all normal reflections outside mask_radius dumped
-  * The new ReflectionContext will occupy a non-contigious block of memory
-  */
- ReflectionContext *apply_mask_radius(double mask_radius, ReflectionContext *rctx) {
- 	Reflection *r, *newr;
- 	ReflectionContext *nrc;
- 	double modul;
- 	//printf("a_m_r: points starting %d\n",reflect_count(nrc));
- 	
- 	nrc = reflection_init();
- 	r = rctx->reflections;
- 	newr = nrc->reflections;
- 	do {
- 		if (r->type == REFLECTION_NORMAL) {
- 			modul = modulus(r->x,r->y,r->z);
- 			printf("a_m_r: modul/mask_radius %f\n",modul/mask_radius);
- 			if (modul < mask_radius) {
- 				newr = malloc(sizeof(Reflection));
- 				memcpy(newr,r,sizeof(Reflection));
- 				reflection_add_from_reflection(nrc,newr);
- 			}
- 		}
- 	} while ((r=r->next) != NULL);
- 	//printf("a_m_r: points remaining %d\n",reflect_count(nrc));
- 	return nrc;
- }
- 
- /*
-  * checks if the supplied reflection falls in octtree volume vol
-  * returns -1 if it falls outside it
-  * returns 0-7 if it falls within depending which child volume it would occupy
-  */
- static int in_octtree_volume(Reflection *r, OctTree *vol) {
- 	int val=0;
- 	double minx,miny,minz;
- 	double maxx,maxy,maxz;
- 	minx = vol->ox-vol->halfedge;
- 	miny = vol->oy-vol->halfedge;
- 	minz = vol->oz-vol->halfedge;
- 	maxx = vol->ox+vol->halfedge;
- 	maxy = vol->oy+vol->halfedge;
- 	maxz = vol->oz+vol->halfedge;
- 	
- 	if ( r->x > maxx || r->x <= minx || r->y > maxy || r->y <= miny || r->z > maxz || r->z <= minz ) {
- 		return -1;
- 	} else {
- 		if (r->x <= vol->ox) val += 1;
- 		if (r->y <= vol->oy) val += 2;
- 		if (r->z <= vol->oz) val += 4;
- 	}
- 	return val;
- }
- /*
-  * set the x,y,z and halfedge params for a octtree child based on the parent and child#
-  */
- static void set_octtree_origin(OctTree *parent, OctTree *child, int childnum) {
- 	int val = childnum;
- 	child->halfedge = parent->halfedge * 0.5;
- 	if (val >= 4) {
- 		child->oz = parent->oz - child->halfedge;
- 		val %= 4;	
- 	} else {
- 		child->oz = parent->oz + child->halfedge;
- 	}
- 	if (val >= 2) {
- 		child->oy = parent->oy - child->halfedge;
- 		val %= 2;	
- 	} else {
- 		child->oy = parent->oy + child->halfedge;
- 	}
- 	if (val >= 1) {
- 		child->ox = parent->ox - child->halfedge;	
- 	} else {
- 		child->ox = parent->ox + child->halfedge;
- 	}
- }
- 
- /*
-  * checks to see if there are any reflections in this volume
-  * if there are, create a new octtree and attach it to the appropriate branch of the parent
-  * else attach null
-  * if we haven't reached maximum depth, spawn 8 new requests until the desired resolution is reached
-  */
- static void stack_octtree(Reflection **rl, int rcount, OctTree *parent, int childnum, int maxdepth) {
- 	//if there are reflections here
- 	//printf("s_o: starting depth=%d child=%d\n",parent->depth,childnum);
- 	if (rcount > 0) {
- 		//printf("s_o: reflections=%d\n",rcount);
- 		OctTree *here = malloc(sizeof(OctTree)); //create a new OctTree node
- 		here->child = calloc(8,sizeof(OctTree *));
- 		here->parent = parent;
- 		here->childnum = childnum;
- 		here->list = NULL;
- 		parent->child[childnum] = here; //attach it to the parent
- 		set_octtree_origin(parent,here,childnum);
- 		//printf("s_o: set origin (%f,%f,%f) halfedge %f\n",here->ox,here->oy,here->oz,here->halfedge);
- 		
- 		if ((here->depth = parent->depth+1) < maxdepth) { //only process children if we're not at maxdepth
- 			//printf("s_o: allocating rcount=%d\n",rcount);
- 			//int *dest = calloc(rcount,sizeof(int)); //list of the reflections and which child to route them to
- 			//int *distrib = calloc(8,sizeof(int)); //count of reflections to route to each child
- 			int dest[rcount];
- 			int distrib[8] = {0,0,0,0,0,0,0,0};
- 			Reflection **list;
-	 		
- 			int i,j;
- 			for (i=0;i<rcount;i++) {
- 				dest[i] = in_octtree_volume(rl[i],here);
- 				//printf("s_o: reflection %d (%f,%f,%f) in volume %d\n",i,rl[i]->x,rl[i]->y,rl[i]->z,dest[i]);
-		 		distrib[dest[i]] += 1;	
- 			}
- 			int n;
- 			for (i=0;i<8;i++) {
- 				if (distrib[i] > 0) {
-	 				//printf("s_o: creating list for child %d, %d members\n",i,distrib[i]);
- 					list = malloc(sizeof(Reflection *)*distrib[i]);
- 					n=0;
- 					for (j=0;j<rcount;j++) {
- 						if (dest[j] == i) {
- 							list[n++] = rl[j];
-	 					}
- 					}
- 					stack_octtree(list,distrib[i],here,i,maxdepth);
-	 				free(list);
- 				} else {
- 					//printf("s_o: no reflections for child %d\n",i);
- 					here->child[i]=NULL;
- 				}
- 			}
- 			//printf("s_o: [%d] ready to free distrib=%d dest=%d\n",here->depth,distrib,dest);
- 			//free(distrib);
- 			//printf("s_o: freed distrib\n");
- 			//free(dest);
- 			//printf("s_o: freed dest\n");
- 		} else { //add a reflection list of the children
- 			ReflectionList *l = malloc(sizeof(ReflectionList));
- 			l->r = malloc(rcount*sizeof(Reflection *));
- 			memcpy(l->r,rl,rcount*sizeof(Reflection *));
- 			here->list = l;
- 		}
- 	} else { //if there are no reflections, just attach NULL and return
- 		printf("s_o: no reflections, attaching null\n");
- 		parent->child[childnum] = NULL;
- 	}
- }
- 
- 
- /*
-  * generate an octtree filling all space out to the largest dimension in the reflectionlist
-  * then eliminate all volumes containing no reflections down to the desired accuracy
-  * TODO: use the same basis as the reflection
-  */
- 
- OctTree *gen_octtree(ReflectionContext *rctx, int depth) {
- 	printf("g_o: starting\n");
- 	double max = largest_dimension(rctx)*1.01;
- 	int count = reflect_count(rctx);
- 	Reflection *r;
- 	
- 	OctTree *top;
- 	top = malloc(sizeof(OctTree));
- 	top->child = calloc(8,sizeof(OctTree *));
- 	top->parent=NULL;
- 	top->halfedge = max;
- 	top->ox = 0;
- 	top->oy = 0;
- 	top->oz = 0;
- 	top->depth=0;
- 	top->childnum=-1;
- 	
- 	Reflection *rl[count];
- 	r = rctx->reflections;
- 	int n=0;
- 	do {
- 		rl[n++] = r;
- 	} while ((r=r->next)!=NULL);
- 	
- 	int i;
- 	for (i=0;i<8;i++) {
- 		//printf("g_o: stack %d\n",i);
- 		stack_octtree(rl,count,top,i,depth);
- 	}
- 	return top;
- }
- 
- static void print_octtree_stack(OctTree *here, int* at_level) {
- 		int i;
- 		for (i=0;i<8;i++) {
- 			if (here->child[i] != NULL) print_octtree_stack(here->child[i],at_level);	
- 		}
- 		at_level[here->depth]++;
- }
- 
- void print_octtree(OctTree *tree) {
- 		int* at_level = calloc(16,sizeof(int));
- 		print_octtree_stack(tree,at_level);
- 		int i;
- 		for (i=0;i<16;i++) {
- 			printf("level %d nodes %d\n",i,at_level[i]);
- 		}
- }
- 
- //return a list of pointers to the 26 surrounding nodes
- OctTreeLinkedList *get_adjacent_nodes(OctTree *o, int *count) {
- 		
- }
- 
- //return a list of all level x nodes
- void stack_get_depth(OctTree *o, int *maxdepth) {
- 	int i;
- 	if (o->depth > *maxdepth) *maxdepth = o->depth;
- 	for (i=0;i<8;i++) {
- 		if (o->child[i] != NULL) stack_get_depth(o->child[i],maxdepth);
- 	}
- }
- 
- OctTreeLinkedList *get_bottom_nodes(OctTree *o, int *count, int level) {
- 		int depth=0;
- 		stack_get_depth(o,&depth);
- 		
- }
- 
- void free_linked_list(OctTreeLinkedList *otll) {
- 	OctTreeLinkedList *o,*next;
- 	o = otll;
- 	do {
- 		next = o;
- 		free(o);
- 		o = next;
- 	} while (o != NULL);	
- }
- 
- void dump_histogram(ReflectionContext *rctx) {
- 	Reflection *r;
- 	int count = reflect_count(rctx);
- 	int n=0;
- 	double dist;
- 	Reflection **rl = malloc(count*sizeof(Reflection *));
- 	r = rctx->reflections;
- 	do {
- 		rl[n++] = r;
- 	} while ((r = r->next) != NULL);
- 	
- 	FILE *f;
- 	
- 	f = fopen("histogram","w");
- 	int i,j;
- 	for (i=0;i<count;i++) {
- 		for (j=i+1;j<count;j++) {
- 			dist = modulus(rl[i]->x-rl[j]->x,rl[i]->y-rl[j]->y,rl[i]->z-rl[j]->z);
- 			fprintf(f,"%f\n",dist); 
- 		}	
- 	}
- 	fclose(f);	
- 	
- }
- 
- /*
-  * look for sections of the tree with gaps of at least req_length between branches
-  * add a node that branches after at least req_length or terminates
-  */  
- OctTreeLinkedList *stack_find_sparse_trees(OctTree *o, OctTreeLinkedList *l, int *count, int req_length, int *cur_length, int allow_end) {
- 	OctTreeLinkedList *nl;
- 	int i;
- 	int children=0;
- 	for (i=0;i<8;i++) {
- 		if (o->child[i] != NULL) children++;
- 	}
- 	if (children==0) { //end of a chain, add if sufficiently long
- 		(*cur_length)++;
- 		if (allow_end==1) {
- 			if (*cur_length >= req_length) {
- 				nl = malloc(sizeof(OctTreeLinkedList));
- 				nl->o = o;
- 				nl->next = NULL;
- 				l->next = nl;
- 				(*count)++;
- 				printf("s_f_s_t: found end-of-chain length=%d depth=%d\n",*cur_length,o->depth);
- 			} else {
- 				nl = l;
- 			}
- 		} else {
- 			nl = l;
- 		}
- 		(*cur_length)=0;
- 		return nl; //return the current list pointer
- 	} else {
- 		if (children==1) { //middle of a singular chain, add to length counter
- 			(*cur_length)++;
- 			nl = l;
- 		} else {
- 			if (*cur_length >= req_length) { //branch point, see if the current chain is long enough to add
-	 			nl = malloc(sizeof(OctTreeLinkedList));
- 				nl->o = o;
- 				nl->next = NULL;
- 				l->next = nl;	
- 				(*count)++;
- 				printf("s_f_s_t: found branch point length=%d depth=%d\n",*cur_length,o->depth);
- 			} else {
- 				nl = l;
- 			}
- 			(*cur_length)=0; //regardless whether this section was long enough, zero the counter
- 		}
- 		
- 		for (i=0;i<8;i++) {
- 			if (o->child[i] != NULL) {
- 				nl = stack_find_sparse_trees(o->child[i],nl,count,req_length,cur_length,allow_end);	
- 			}	
- 		}
- 		return nl;	
- 	}
- 	
- 	
- }
- 
- OctTreeLinkedList *find_sparse_trees(OctTree *o, int req_length, int allow_end, int *count) {
- 	printf("f_s_t: starting\n");
- 	int cur_length=0;
- 	OctTreeLinkedList *ol = malloc(sizeof(OctTreeLinkedList));
- 	OctTreeLinkedList *ol2;
- 	ol->o = NULL;
- 	ol->next = NULL;
- 	stack_find_sparse_trees(o,ol,count,req_length,&cur_length,allow_end);
- 	ol2 = ol->next;
- 	free(ol);
- 	return ol2;
- }
- 
- ReflectionContext *change_reflection_basis(ReflectionContext *rctx, Basis *basis) {
- 	ReflectionContext *new = reflection_init();
- 	Reflection *r;
- 
- 	//calculate a change-of-basis matrix, replace all the reflection coordinates thus
- 	//hence we hopefully have a square-basis representation, which octtree statistics might work on. maybe.
- }
- 
- 
- 
- 
- 
diff --git a/src/structure.h b/src/structure.h
deleted file mode 100644
index feac924..0000000
--- a/src/structure.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/*
- * structure.h
- *
- * 3D analysis
- *
- * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
- *  dtr - Diffraction Tomography Reconstruction
- *
- */
- #ifndef REFLECTION_H
- #include "reflection.h"
- #endif
- 
-#ifndef STRUCTURE_H_
-#define STRUCTURE_H_
-
-typedef struct octtree {
-	struct octtree *parent;
-	struct octtree **child;
-	struct reflist *list;
-	double ox;
-	double oy;
-	double oz;
-	double halfedge;
-	int depth;
-	int childnum;
-} OctTree;
-
-typedef struct reflist {
-	int count;
-	struct Reflection **r;
-} ReflectionList;
-
-typedef struct octtreelist {
-	struct octtree *o;
-	struct octtreelist *next;
-} OctTreeLinkedList;
-
-extern double get_mask_radius(ReflectionContext *rctx);
-extern ReflectionContext *apply_mask_radius(double mask_radius, ReflectionContext *rctx);
-extern OctTree *gen_octtree(ReflectionContext *rctx, int depth);
-extern void print_octtree(OctTree *tree);
-extern void dump_histogram(ReflectionContext *rctx);
-extern OctTreeLinkedList *find_sparse_trees(OctTree *o, int req_length, int allow_end, int *count);
-
-
-#endif /*STRUCTURE_H_*/