Add skeleton prealignment procedure

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

11 files changed, 122 insertions, 49 deletions

diff --git a/Makefile.am b/Makefile.am
index 6fa95bf..bd1c637 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/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h src/itrans-threshold.h \
-		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h
+		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h
 SUBDIRS = src data
diff --git a/Makefile.in b/Makefile.in
index 2fe4168..c54ac8d 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/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h src/itrans-threshold.h \
-		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h
+		src/itrans-zaefferer.h src/itrans-lsq.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h
 
 SUBDIRS = src data
 all: config.h
diff --git a/src/Makefile.am b/src/Makefile.am
index 3705cfc..59c995c 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 \
-		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c
+		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.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 e8f6064..0ccaf65 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -53,7 +53,7 @@ am_dtr_OBJECTS = main.$(OBJEXT) displaywindow.$(OBJEXT) \
 	itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT) cache.$(OBJEXT) \
 	itrans-threshold.$(OBJEXT) itrans-zaefferer.$(OBJEXT) \
 	itrans-lsq.$(OBJEXT) itrans-stat.$(OBJEXT) control.$(OBJEXT) \
-	mapping.$(OBJEXT) reproject.$(OBJEXT)
+	mapping.$(OBJEXT) reproject.$(OBJEXT) prealign.$(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 \
-		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c
+		itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c
 
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
@@ -237,6 +237,7 @@ distclean-compile:
 @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)/prealign.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@
diff --git a/src/control.h b/src/control.h
index 83ba0fb..68eebed 100644
--- a/src/control.h
+++ b/src/control.h
@@ -79,6 +79,7 @@ typedef struct cctx_struct {
 	FormulationMode		fmode;
 	ReconstructionMode 	rmode;
 	PeakSearchMode 		psmode;
+	unsigned int		prealign;
 	
 	/* Input filename */
 	char 			*filename;
@@ -116,6 +117,7 @@ typedef struct cctx_struct {
 	/* GTK bits */
 	GtkWidget		*combo_algorithm;
 	GtkWidget		*combo_peaksearch;
+	GtkWidget		*checkbox_prealign;
 
 } ControlContext;
 
diff --git a/src/ipr.c b/src/ipr.c
index 586a468..0296765 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -202,9 +202,9 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 	}
 	
 	ordered = reflection_init();
-	basis->a.x = 5e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
-	basis->b.x = 0.0;  basis->b.y = 5e9;  basis->b.z = 0.0;
-	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 3e9;
+	//basis->a.x = 5e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
+	//basis->b.x = 0.0;  basis->b.y = 5e9;  basis->b.z = 0.0;
+	//basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 3e9;
 	for ( h=-30; h<=30; h++ ) {
 		for ( k=-30; k<=30; k++ ) {
 			for ( l=-30; l<=30; l++ ) {
@@ -258,7 +258,7 @@ int ipr_refine(ControlContext *ctx) {
 		
 		/* Select an image */
 		i = ipr_random_image(ctx);
-		i = 1;
+		i = 25;
 		cur = imagedisplay_open(ctx->images[i].image, ctx->images[i].width, ctx->images[i].height, "Current Image");
 		
 		refl = reproject_get_reflections(ctx->images[i], &n, lat, ctx);
diff --git a/src/main.c b/src/main.c
index 4e5723a..5df4bcf 100644
--- a/src/main.c
+++ b/src/main.c
@@ -31,6 +31,7 @@
 #include "ipr.h"
 #include "cache.h"
 #include "mapping.h"
+#include "prealign.h"
 
 static gint main_method_window_response(GtkWidget *method_window, gint response, ControlContext *ctx) {
 
@@ -53,6 +54,12 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			default: ctx->psmode = PEAKSEARCH_NONE; break;	/* This happens when reading from a cache file */
 		}
 		
+		if ( gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(ctx->checkbox_prealign)) ) {
+			ctx->prealign = TRUE;
+		} else {
+			ctx->prealign = FALSE;
+		}
+		
 		gtk_widget_destroy(method_window);
 		while ( gtk_events_pending() ) gtk_main_iteration();
 		
@@ -66,6 +73,10 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			val=0;
 		}
 		
+		if ( ctx->prealign ) {
+			prealign_do_series(ctx);
+		}
+		
 		if ( ctx->inputfiletype != INPUT_CACHE ) {
 			mapping_create(ctx);
 		}
@@ -108,9 +119,10 @@ void main_method_dialog_open(ControlContext *ctx) {
 	vbox = gtk_vbox_new(FALSE, 0);
 	hbox = gtk_hbox_new(TRUE, 0);
 	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(method_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
-	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5);
+	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 10);
 	
-	table = gtk_table_new(3, 2, FALSE);
+	table = gtk_table_new(4, 2, FALSE);
+	gtk_table_set_row_spacings(GTK_TABLE(table), 5);
 
 	method_label = gtk_label_new("Reconstruction Algorithm: ");
 	gtk_table_attach_defaults(GTK_TABLE(table), method_label, 1, 2, 1, 2);
@@ -134,6 +146,9 @@ void main_method_dialog_open(ControlContext *ctx) {
 	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3);
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 2, 3);
 	
+	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);
+	
 	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 a8d50b3..d9ca94d 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -31,52 +31,57 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	/* 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;
+	if ( !ctx->prealign ) {
+		/* 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++ ) {
-				image_total[x + twidth*y] += ctx->images[i].image[x + twidth*y]/ctx->n_images;
+				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;
 	
-	/* 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;
-			}
-		}
+		/* 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);
+		
 	}
-	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++ ) {
-		ctx->images[i].x_centre = ctx->x_centre;
-		ctx->images[i].y_centre = ctx->y_centre;
+		if ( !ctx->prealign ) {
+			ctx->images[i].x_centre = ctx->x_centre;
+			ctx->images[i].y_centre = ctx->y_centre;
+		}
 		itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt);
 	}
 	
diff --git a/src/prealign.c b/src/prealign.c
new file mode 100644
index 0000000..86c0089
--- /dev/null
+++ b/src/prealign.c
@@ -0,0 +1,26 @@
+/*
+ * prealign.c
+ *
+ * Rough initial alignment of images
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "control.h"
+
+/* No peak-detection nor 3D mapping has been done yet.
+	Ask the user to give a rough idea (i.e. as accurately as possible...)
+	of the centre of each image. */
+void prealign_do_series(ControlContext *ctx) {
+
+	
+
+}
+
diff --git a/src/prealign.h b/src/prealign.h
new file mode 100644
index 0000000..53ac000
--- /dev/null
+++ b/src/prealign.h
@@ -0,0 +1,24 @@
+/*
+ * prealign.h
+ *
+ * Rough initial alignment of images
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+ 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef PREALIGN_H
+#define PREALIGN_H
+
+#include "control.h"
+
+extern void prealign_do_series(ControlContext *ctx);
+
+#endif	/* PREALIGN_H */
+
diff --git a/src/reproject.c b/src/reproject.c
index f3c3480..bec9a31 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -81,7 +81,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			ydd = (yddd/cos(tilt) - zddd*tan(tilt)/cos(tilt))/(1+tan(tilt)*tan(tilt));
 			zdd = (-zddd-ydd*sin(tilt))/cos(tilt);
 			
-			yd = (ydd-xdd*tan(omega))/(sin(omega)*tan(omega)+cos(omega));
+			yd = (ydd-xdd*tan(omega))/(sin(omega)*tan(omega)+cos(omega));	/* Needs to be done first */
 			xd = (xdd+yd*sin(omega))/cos(omega);
 			zd = zdd;
 			
@@ -117,7 +117,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 					break;
 				}
 				
-				printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
+				//printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
 			
 			} else {
 				//fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y);
@@ -129,7 +129,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 
 	} while ( reflection );
 	
-	printf("Found %i reflections in image\n", i);
+	//printf("Found %i reflections in image\n", i);
 	*n = i;
 	return refl;