Initial simplex refinement stuff

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

7 files changed, 98 insertions, 529 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index d51f2bb..e9aae1c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,13 +1,10 @@
-bin_PROGRAMS = dtr iprtest
+bin_PROGRAMS = dtr
 
 dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c cache.c \
 		itrans-threshold.c itrans-zaefferer.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \
 		dirax.c image.c refine.c gtk-valuegraph.c intensities.c glbits.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil
 
-iprtest_SOURCES = iprtest.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c
-iprtest_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil
-
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
 
diff --git a/src/displaywindow.c b/src/displaywindow.c
index af17934..5bd00fc 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -339,12 +339,12 @@ static gint displaywindow_image_last(GtkWidget *widget, DisplayWindow *dw) {
 }
 
 static gint displaywindow_refinestep(GtkWidget *widget, DisplayWindow *dw) {
-	refine_do_image(dw->ctx);
+	refine_do_cell(dw->ctx);
 	return 0;
 }
 
 static gint displaywindow_refinestack(GtkWidget *widget, DisplayWindow *dw) {
-	refine_do_stack(dw->ctx);
+	refine_do_sequence(dw->ctx);
 	return 0;
 }
 
@@ -392,8 +392,8 @@ static void displaywindow_addmenubar(DisplayWindow *dw) {
 		{ "DirAxAction", "dtr-dirax", "Start _DirAx", "<Ctrl>D", NULL, G_CALLBACK(displaywindow_dirax) },
 		{ "DirAxReRunAction", NULL, "Run another DirAx cycle", NULL, NULL, G_CALLBACK(displaywindow_dirax_rerun) },
 		{ "StopDirAxAction", NULL, "Stop DirAx", NULL, NULL, G_CALLBACK(displaywindow_dirax_stop) },
-		{ "RefineStepAction", "dtr-refine", "Refine Current Image", NULL, NULL, G_CALLBACK(displaywindow_refinestep) },
-		{ "RefineSeqAction", "dtr-refineall", "Refine Entire Stack", NULL, NULL, G_CALLBACK(displaywindow_refinestack) },
+		{ "RefineStepAction", "dtr-refine", "Refine Unit Cell", NULL, NULL, G_CALLBACK(displaywindow_refinestep) },
+		{ "RefineSeqAction", "dtr-refineall", "Run Refinement Sequence", NULL, NULL, G_CALLBACK(displaywindow_refinestack) },
 		{ "SetAxisAction", "dtr-tiltaxis", "Set Tilt Axis Position...", NULL, NULL, G_CALLBACK(displaywindow_setaxis) },
 		{ "IncrAxisAction", NULL, "Increase Tilt Axis Position", "<Ctrl>Up", NULL, G_CALLBACK(displaywindow_incraxis) },
 		{ "DecrAxisAction", NULL, "Decrease Tilt Axis Position", "<Ctrl>Down", NULL, G_CALLBACK(displaywindow_decraxis) },
@@ -404,8 +404,8 @@ static void displaywindow_addmenubar(DisplayWindow *dw) {
 		
 		{ "ButtonDirAxAction", "dtr-dirax", "Run _DirAx", NULL, NULL, G_CALLBACK(displaywindow_dirax) },
 		{ "ButtonTiltAxisAction", "dtr-tiltaxis", "Tilt Axis", NULL, NULL, G_CALLBACK(displaywindow_setaxis) },
-		{ "ButtonRefineStepAction", "dtr-refine", "Refine Image", NULL, NULL, G_CALLBACK(displaywindow_refinestep) },
-		{ "ButtonRefineSeqAction", "dtr-refineall", "Refine Stack", NULL, NULL, G_CALLBACK(displaywindow_refinestack) },
+		{ "ButtonRefineStepAction", "dtr-refine", "Refine", NULL, NULL, G_CALLBACK(displaywindow_refinestep) },
+		{ "ButtonRefineSeqAction", "dtr-refineall", "Sequence", NULL, NULL, G_CALLBACK(displaywindow_refinestack) },
 		{ "ButtonFirstImageAction", GTK_STOCK_GOTO_FIRST, "First Image", NULL, NULL, G_CALLBACK(displaywindow_image_first) },
 		{ "ButtonPrevImageAction", GTK_STOCK_GO_BACK, "Previous Image", NULL, NULL, G_CALLBACK(displaywindow_image_prev) },
 		{ "ButtonNextImageAction", GTK_STOCK_GO_FORWARD, "Next Image", NULL, NULL, G_CALLBACK(displaywindow_image_next) },
diff --git a/src/image.c b/src/image.c
index ebeff15..d8bba0c 100644
--- a/src/image.c
+++ b/src/image.c
@@ -3,7 +3,7 @@
  *
  * Handle images and image features
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
diff --git a/src/image.h b/src/image.h
index 6b2f4e7..62a4d2c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -3,7 +3,7 @@
  *
  * Handle images and image features
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
diff --git a/src/iprtest.c b/src/iprtest.c
deleted file mode 100644
index 369fd51..0000000
--- a/src/iprtest.c
+++ /dev/null
@@ -1,180 +0,0 @@
-/*
- * iprtest.c
- *
- * Unit test for IPR
- *
- * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
- *
- *  dtr - Diffraction Tomography Reconstruction
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <gtk/gtk.h>
-#include <stdlib.h>
-#include <math.h>
-
-#include "basis.h"
-#include "reflections.h"
-#include "image.h"
-#include "reproject.h"
-#include "refine.h"
-#include "utils.h"
-#include "mapping.h"
-
-typedef struct {
-	
-	ImageRecord *image;
-	
-	Basis *cell;
-	ReflectionList *reflections;
-	ReflectionList *gen;
-	
-} IPRTestContext;
-
-static gint main_window_closed(GtkWidget *window, IPRTestContext *ctx) {
-	gtk_exit(0);
-	return 0;
-}
-
-static gboolean main_window_expose(GtkWidget *drawing, GdkEventExpose *event, IPRTestContext *ctx) {
-
-	cairo_t *cr;
-	ImageFeature *f;
-	size_t i;
-	
-	cr = gdk_cairo_create(drawing->window);
-
-	cairo_new_path(cr);
-	cairo_rectangle(cr, 0.0, 0.0, 640.0, 640.0);
-	cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);
-	cairo_fill(cr);
-	
-	f = ctx->image->features->features;
-	cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
-	for ( i=0; i<ctx->image->features->n_features; i++ ) {
-		cairo_new_path(cr);
-		cairo_arc(cr, f[i].x, f[i].y, 4.0, 0, 2*M_PI);
-		cairo_fill(cr);
-	}
-	
-	f = ctx->image->rflist->features;
-	cairo_set_line_width(cr, 1.0);
-	for ( i=0; i<ctx->image->rflist->n_features; i++ ) {
-	
-		cairo_new_path(cr);
-		cairo_arc(cr, f[i].x, f[i].y, 8.0, 0, 2*M_PI);
-		
-		if ( f[i].partner != NULL ) {
-		
-			cairo_set_source_rgb(cr, 0.0, 0.8, 0.0);
-			cairo_stroke(cr);
-		
-			cairo_new_path(cr);
-			cairo_move_to(cr, f[i].x, f[i].y);
-			cairo_line_to(cr, f[i].partner->x, f[i].partner->y);
-			cairo_stroke(cr);
-			
-		} else {
-
-			cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
-			cairo_stroke(cr);
-		
-		}
-		
-	}
-	
-	cairo_destroy(cr);
-
-	return FALSE;
-
-}
-
-static gboolean main_button_press(GtkWidget *drawing, GdkEventButton *event, IPRTestContext *ctx) {
-
-	if ( (event->type == GDK_BUTTON_PRESS) && (event->button == 1) ) {
-		refine_fit_image(ctx->cell, ctx->image);
-		image_feature_list_free(ctx->image->rflist);
-		reflection_list_from_new_cell(ctx->gen, ctx->cell);
-		ctx->image->rflist = reproject_get_reflections(ctx->image, ctx->gen);
-		gdk_window_invalidate_rect(drawing->window, &drawing->allocation, FALSE);
-	}
-	
-	return FALSE;
-}
-
-int main(int argc, char *argv[]) {
-
-	GtkWidget *window;
-	GtkWidget *drawing;
-	IPRTestContext *ctx;
-	
-	gtk_init(&argc, &argv);
-	
-	ctx = malloc(sizeof(IPRTestContext));
-	
-	ctx->image = malloc(sizeof(ImageRecord));
-	ctx->image->image = NULL;
-	ctx->image->tilt = deg2rad(0.0);
-	ctx->image->omega = deg2rad(0.0);
-	ctx->image->slop = deg2rad(0.0);
-	ctx->image->fmode = FORMULATION_PIXELSIZE;
-	ctx->image->pixel_size = 0.6e8;
-	ctx->image->camera_len = 0.0;
-	ctx->image->lambda = lambda(1000000);
-	ctx->image->resolution = 1.0;
-	ctx->image->width = 640;
-	ctx->image->height = 640;
-	ctx->image->x_centre = 320;
-	ctx->image->y_centre = 320;
-	ctx->image->features = NULL;
-	ctx->image->rflist = NULL;
-	
-	ctx->cell = malloc(sizeof(Basis));
-	/* The "true" cell */
-	ctx->cell->a.x = 5.0e9;  ctx->cell->a.y = 0.0e9;  ctx->cell->a.z = 0.0e9;
-	ctx->cell->b.x = 0.0e9;  ctx->cell->b.y = 5.0e9;  ctx->cell->b.z = 0.0e9;
-	ctx->cell->c.x = 0.0e9;  ctx->cell->c.y = 0.0e9;  ctx->cell->c.z = 5.0e9;
-	ctx->reflections = reflection_list_from_cell(ctx->cell);
-	ctx->image->features = reproject_get_reflections(ctx->image, ctx->reflections);
-	/* The "model" cell */
-	ctx->cell->a.x = 5.2e9;  ctx->cell->a.y = 0.1e9;  ctx->cell->a.z = 0.0e9;
-	ctx->cell->b.x = 0.2e9;  ctx->cell->b.y = 4.8e9;  ctx->cell->b.z = 0.0e9;
-	ctx->cell->c.x = 0.0e9;  ctx->cell->c.y = 0.0e9;  ctx->cell->c.z = 5.0e9;
-	ctx->gen = reflection_list_from_cell(ctx->cell);
-	ctx->image->rflist = reproject_get_reflections(ctx->image, ctx->gen);
-	
-	window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-	gtk_window_set_title(GTK_WINDOW(window), "IPR Test");
-	gtk_widget_set_size_request(GTK_WIDGET(window), 640, 640);
-	g_signal_connect(G_OBJECT(window), "destroy", G_CALLBACK(main_window_closed), ctx);
-	gtk_window_set_resizable(GTK_WINDOW(window), FALSE);
-	
-	drawing = gtk_drawing_area_new();
-	gtk_container_add(GTK_CONTAINER(window), drawing);
-	g_signal_connect(G_OBJECT(drawing), "expose_event", G_CALLBACK(main_window_expose), ctx);
-	gtk_widget_add_events(drawing, GDK_BUTTON_PRESS_MASK);
-	g_signal_connect(GTK_OBJECT(drawing), "button_press_event", G_CALLBACK(main_button_press), ctx);
-	
-	gtk_widget_show_all(window);
-	
-	gtk_main();
-	
-	return 0;
-	
-}
-
-/* Dummy function stubs */
-#include "displaywindow.h"
-#include "gtk-valuegraph.h"
-void displaywindow_update_imagestack(DisplayWindow *dw) { };
-void displaywindow_enable_cell_functions(DisplayWindow *dw, gboolean g) { };
-void displaywindow_update(DisplayWindow *dw) { };
-void displaywindow_error(const char *msg, DisplayWindow *dw) { };
-guint gtk_value_graph_get_type() { return 0; };
-GtkWidget *gtk_value_graph_new() { return NULL; };
-void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) { };
-
diff --git a/src/refine.c b/src/refine.c
index b18e5fb..45aec8e 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -19,394 +19,149 @@
 #include <stdio.h>
 #include <assert.h>
 #include <string.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
 
+#include "control.h"
 #include "displaywindow.h"
-#include "gtk-valuegraph.h"
-#include "basis.h"
-#include "reflections.h"
 #include "image.h"
 #include "reproject.h"
-#include "control.h"
 #include "mapping.h"
-#include "imagedisplay.h"
-#include "utils.h"
 
-typedef enum {
-	AXIS_X = 1,
-	AXIS_Y = 2
-} Axis;
+/* A simplex is an array of ten of these */
+typedef struct {
+	double dax;	double dbx;	double dcx;
+	double day;	double dby;	double dcy;
+	double daz;	double dbz;	double dcz;
+} SimplexVertex;
+
+typedef struct {
+	signed int h;	signed int k;	signed int l;
+	double dx;	double dy;	double dz;
+} Deviation;
+
+void refine_do_sequence(ControlContext *ctx) {
+
+	
+}
 
-static void refine_fix_unconstrained(gsl_matrix *M) {
+void refine_do_cell(ControlContext *ctx) {
 
-	int row, j;
-	int modified = 0;
+	SimplexVertex s[10];
+	Deviation *d;
+	const double delta = 0.1e9;
+	int i, nf, f, v_worst;
+	double fom_worst;
 	
-	/* Find a row which is all zero */
-	for ( row=0; row<M->size1; row++ ) {
+	if ( !ctx->cell_lattice ) {
+		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
+		return;
+	}
 	
-		int zero = 1;
+	if ( ctx->images->n_images == 0 ) {
+		displaywindow_error("There are no images to refine against.", ctx->dw);
+		return;
+	}
+	
+	/* Initialise the simplex */
+	s[0].dax = 0.0;	s[0].dbx = 0.0;	s[0].dcx = 0.0;
+	s[0].day = 0.0;	s[0].dby = 0.0;	s[0].dcy = 0.0;
+	s[0].daz = 0.0;	s[0].dbz = 0.0;	s[0].dcz = 0.0;
+	memcpy(&s[1], &s[0], sizeof(SimplexVertex));	s[1].dax = delta;
+	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].dbx = delta;
+	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].dcx = delta;
+	memcpy(&s[4], &s[0], sizeof(SimplexVertex));	s[4].day = delta;
+	memcpy(&s[5], &s[0], sizeof(SimplexVertex));	s[5].dby = delta;
+	memcpy(&s[6], &s[0], sizeof(SimplexVertex));	s[6].dcy = delta;
+	memcpy(&s[7], &s[0], sizeof(SimplexVertex));	s[7].daz = delta;
+	memcpy(&s[8], &s[0], sizeof(SimplexVertex));	s[8].dbz = delta;
+	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
+	
+	/* Create the table of indicies and deviations */
+	nf = 0;
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		int j;
 		
-		for ( j=0; j<M->size2; j++ ) {
-			if ( gsl_matrix_get(M, row, j) != 0.0 ) {
-				zero = 0;
-				break;
-			}
+		if ( !ctx->images->images[i].rflist ) {
+			ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice);
 		}
 		
-		if ( zero ) {
-			
-			/* Find a column which is all zero */
-			int i, col;
-			int row_altered = 0;
-			
-			for ( col=0; col<M->size2; col++ ) {
-			
-				int zero2 = 1;
-				
-				if ( row_altered ) break; /* This row is now fixed, so move on to the next */
-				
-				for ( i=0; i<M->size1; i++ ) {
-					if ( gsl_matrix_get(M, i, col) != 0.0 ) {
-						zero2 = 0;
-						break;
-					}
-				}
-				
-				if ( zero2 ) {
-					gsl_matrix_set(M, row, col, 1.0);
-					modified = 1;
-					row_altered = 1;
-				}
-			
-			}
-			
+		for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) {
+			if ( ctx->images->images[i].rflist->features[j].partner != NULL ) nf++;
 		}
 		
 	}
+	printf("RF: There are %i partnered features in total\n", nf);
 	
-	if ( modified) printf("RF: M was modified to fix unconstrained variables.\n");
-
-}
-
-/* Use the IPR algorithm to make "cell" fit the given image */
-ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
-	
-	ImageFeatureList *flist;
-	int j, ns;
-	Axis axis;
-	double dax = 0.0, dbx = 0.0, dcx = 0.0;
-	double day = 0.0, dby = 0.0, dcy = 0.0;
-	double dlx = 0.0, dly = 0.0, dlz = 0.0;
-	
-	flist = image->features;
-		
-	for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) {
-	
-		gsl_permutation *perm;
-		int s;
-		double det;
-		gsl_matrix *M;
-		gsl_vector *q;
-		gsl_vector *p;
-		gsl_vector *scratch;
-		gsl_matrix *V;
-		gsl_vector *S;
+	d = malloc(nf*sizeof(Deviation));
+	f = 0;
+	for ( i=0; i<ctx->images->n_images; i++ ) {
 		
-		if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n");
-		if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n");
+		ImageRecord *image;
+		int j;
 		
-		M = gsl_matrix_alloc(3, 3);
-		p = gsl_vector_alloc(3);
-		gsl_matrix_set_zero(M);
-		gsl_vector_set_zero(p);
+		image = &ctx->images->images[i];
 		
-		ns = 0;
-		for ( j=0; j<image->rflist->n_features; j++ ) {
+		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
 		
-			double val;
 			ImageFeature *rf;
-			signed int h, k, l;
-			double xy;
 			double dix, diy, dx, dy;
+			double dlx, dly, dlz;
 			double old_x, old_y;
 			
 			rf = &image->rflist->features[j];
-			
 			if ( !rf->partner ) continue;
 			
-			h = rf->reflection->h;
-			k = rf->reflection->k;
-			l = rf->reflection->l;
+			d[f].h = rf->reflection->h;
+			d[f].k = rf->reflection->k;
+			d[f].l = rf->reflection->l;
 			
 			/* Determine the difference vector */
 			dix = rf->partner->x - rf->x;
 			diy = rf->partner->y - rf->y;
-			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy);
+			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy);
 			
 			old_x = rf->partner->x;
 			old_y = rf->partner->y;
 			rf->partner->x = dix + rf->partner->parent->x_centre;
 			rf->partner->y = diy + rf->partner->parent->y_centre;
 			mapping_scale(rf->partner, &dx, &dy);
+			mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
 			rf->partner->x = old_x;
 			rf->partner->y = old_y;
-			printf("=> %+10.5f %+10.5f nm^-1\n", dx/1e9, dy/1e9);
-			
-			xy = 0;
-			switch ( axis ) {
-				case AXIS_X : xy = dx; break;
-				case AXIS_Y : xy = dy; break;
-			}
+			printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9);
 			
-			/* Elements of "p" */
-			val = gsl_vector_get(p, 0);	val += xy * h;	gsl_vector_set(p, 0, val);
-			val = gsl_vector_get(p, 1);	val += xy * k;	gsl_vector_set(p, 1, val);
-			val = gsl_vector_get(p, 2);	val += xy * l;	gsl_vector_set(p, 2, val);
-			gsl_matrix_get(M, 2, 2);
-			/* Elements of "M" */
-			val = gsl_matrix_get(M, 0, 0);	val += h * h;	gsl_matrix_set(M, 0, 0, val);
-			val = gsl_matrix_get(M, 0, 1);	val += k * h;	gsl_matrix_set(M, 0, 1, val);
-			val = gsl_matrix_get(M, 0, 2);	val += l * h;	gsl_matrix_set(M, 0, 2, val);
+			d[f].dx = dlx;
+			d[f].dy = dly;
+			d[f].dz = dlz;
 			
-			val = gsl_matrix_get(M, 1, 0);	val += h * k;	gsl_matrix_set(M, 1, 0, val);
-			val = gsl_matrix_get(M, 1, 1);	val += k * k;	gsl_matrix_set(M, 1, 1, val);
-			val = gsl_matrix_get(M, 1, 2);	val += l * k;	gsl_matrix_set(M, 1, 2, val);
-			
-			val = gsl_matrix_get(M, 2, 0);	val += h * l;	gsl_matrix_set(M, 2, 0, val);
-			val = gsl_matrix_get(M, 2, 1);	val += k * l;	gsl_matrix_set(M, 2, 1, val);
-			val = gsl_matrix_get(M, 2, 2);	val += l * l;	gsl_matrix_set(M, 2, 2, val);
-			
-			ns++;
+			f++;
 		
 		}
 		
-		if ( ns == 0 ) {
-			printf("RF: No partners found\n");
-			gsl_matrix_free(M);
-			gsl_vector_free(p);
-			return NULL;
-		}	
-		
-		/* Prepare for solving */
-		refine_fix_unconstrained(M);
-	//	gsl_matrix_set(M, 2, 0, 0.0);
-	//	gsl_matrix_set(M, 2, 1, 0.0);
-	//	gsl_matrix_set(M, 2, 2, 1.0);
-	//	gsl_matrix_set(M, 0, 2, 0.0);
-	//	gsl_matrix_set(M, 1, 2, 0.0);
-	//	gsl_vector_set(p, 2, 0.0);
-		matrix_vector_show(M, p, "RF: ");
-		
-		/* Calculate and display the determinant */
-		perm = gsl_permutation_alloc(M->size1);
-		gsl_linalg_LU_decomp(M, perm, &s);
-		det = gsl_linalg_LU_det(M, s);
-		printf("RF: Determinant of M = %f\n", det);
-		gsl_permutation_free(perm);		
-		
-		/* Solve the equations */
-		V = gsl_matrix_alloc(M->size2, M->size2);
-		S = gsl_vector_alloc(M->size2);
-		scratch = gsl_vector_alloc(M->size2);
-		gsl_linalg_SV_decomp(M, V, S, scratch);
-		q = gsl_vector_alloc(3);	/* This is the "answer" */
-		gsl_vector_set_zero(q);
-		gsl_linalg_SV_solve(M, V, S, p, q);
-		gsl_vector_free(scratch);
-		gsl_matrix_free(V);
-		gsl_vector_free(S);
-		gsl_matrix_free(M);
-		gsl_vector_free(p);
-		
-		switch ( axis ) {
-			case AXIS_X : {
-				dax = gsl_vector_get(q, 0);
-				dbx = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "x" of the image coordinate  */
-				dcx = gsl_vector_get(q, 2);	/* system, of a,b and c */
-				break;
-			}
-			case AXIS_Y : {
-				day = gsl_vector_get(q, 0);
-				dby = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "y" of the image coordinate  */
-				dcy = gsl_vector_get(q, 2);	/* system, of a,b and c */
-				break;
-			}
-		}
-		
-		gsl_vector_free(q);
-
 	}
+	assert( f == nf );
 	
-	printf("RF: ------------------------------------------------------------------ Refinement results ---------\n");
-	printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9);
-	printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9);
-	printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9);
-	
-	/* Update the cell */
-	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->a.x += dlx;  cell->a.y += dly; cell->a.z += dlz;
-	
-	mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->b.x += dlx;  cell->b.y += dly;  cell->b.z += dlz;
+	/* Find the least favourable vertex of the simplex */
+	v_worst = 0;
+	fom_worst = 0;
+	for ( i=0; i<10; i++ ) {
 	
-	mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
-	
-	return NULL;
-	
-}
-
-static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) {
-
-	int i;
-	double series_dev_max = 0;
-	double series_dev_min = +HUGE_VAL;
-	double series_dev_mean = 0;
-	int series_dev_n = 0;
-	
-	/* Ensure that ctx->cell_lattice is set up */
-	reproject_cell_to_lattice(ctx);
-	
-	for ( i=0; i<ctx->images->n_images; i++ ) {
-
-		ImageRecord *image;
-		int j, n;
-		double image_dev_mean = 0;
-		
-		image = &ctx->images->images[i];
+		double fom = 0;
+		int j;
 		
-		/* Fit this image and update ctx->cell_lattice, index the selected pattern */
-		if ( !image->rflist ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
-		refine_fit_image(ctx->cell, image);
-		reproject_cell_to_lattice(ctx);
-		image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+		for ( j=0; j<nf; j++ ) {
 		
-		n = 0;
-		for ( j=0; j<image->rflist->n_features; j++ ) {
+			double xdev, ydev, zdev;
 			
-			double dix, diy;
+			xdev = 0;
+			ydev = 0;
+			zdev = 0;
 			
-			/* Skip if no partner */
-			if ( !image->rflist->features[j].partner ) continue;
+			fom += sqrt(xdev*xdev + ydev*ydev + zdev*zdev);
 			
-			/* Determine the difference vector */
-			dix = image->rflist->features[j].partner->x - image->rflist->features[j].x;
-			diy = image->rflist->features[j].partner->y - image->rflist->features[j].y;
-			
-			image_dev_mean += sqrt(dix*dix + diy*diy);
-			n++;
-			
-		}
-		image_dev_mean /= n;
-		
-		if ( image_dev_mean > series_dev_max ) series_dev_max = image_dev_mean;
-		if ( image_dev_mean < series_dev_min ) series_dev_min = image_dev_mean;
-		series_dev_mean += image_dev_mean;
-		series_dev_n++;
-		
-	}
-	
-	series_dev_mean /= series_dev_n;
-	*fit = series_dev_mean;
-	*warp = (series_dev_max - series_dev_min)/series_dev_min;
-	
-	return 0;
-
-}
-
-void refine_do_stack(ControlContext *ctx) {
-
-	double omega_offs;
-	GtkWidget *window_fit;
-	GtkWidget *graph_fit;
-	double *fit_vals;
-	GtkWidget *window_warp;
-	GtkWidget *graph_warp;
-	double *warp_vals;
-	size_t idx;
-	
-	fit_vals = malloc(401*sizeof(double));
-	warp_vals = malloc(401*sizeof(double));
-	idx = 0;
-
-	if ( !ctx->cell_lattice ) {
-		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
-		return;
-	}
-	
-	if ( ctx->images->n_images == 0 ) {
-		displaywindow_error("There are no images to refine against.", ctx->dw);
-		return;
-	}
-	
-	for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) {
-	
-		double fit, warp;
-		int i;
-		Basis cell_copy;
-		
-		memcpy(&cell_copy, ctx->cell, sizeof(Basis));
-		for ( i=0; i<ctx->images->n_images; i++ ) {
-			ctx->images->images[i].omega += omega_offs;
-		}
-		
-		if ( refine_sequence_sweep(ctx, &fit, &warp) ) {
-			printf("RF: Sequencer sweep failed\n");
-		}
-		printf("RF: omega_offs=%f, fit=%f, warp=%f\n", omega_offs, fit, warp);
-		fit_vals[idx] = fit;
-		warp_vals[idx++] = warp;
-		
-		for ( i=0; i<ctx->images->n_images; i++ ) {
-			ctx->images->images[i].omega -= omega_offs;
 		}
-		memcpy(ctx->cell, &cell_copy, sizeof(Basis));
 		
 	}
 	
-	displaywindow_update(ctx->dw);
-	reproject_lattice_changed(ctx);
-	
-	window_fit = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-	gtk_window_set_default_size(GTK_WINDOW(window_fit), 640, 256);
-	gtk_window_set_title(GTK_WINDOW(window_fit), "Omega-Search Graph: Fit");
-	graph_fit = gtk_value_graph_new();
-	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_fit), fit_vals, idx);
-	gtk_container_add(GTK_CONTAINER(window_fit), graph_fit);
-	gtk_widget_show_all(window_fit);
-	
-	window_warp = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-	gtk_window_set_default_size(GTK_WINDOW(window_warp), 640, 256);
-	gtk_window_set_title(GTK_WINDOW(window_warp), "Omega-Search Graph: Warp");
-	graph_warp = gtk_value_graph_new();
-	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_warp), warp_vals, idx);
-	gtk_container_add(GTK_CONTAINER(window_warp), graph_warp);
-	gtk_widget_show_all(window_warp);
-	
-}
-
-void refine_do_image(ControlContext *ctx) {
-
-	if ( !ctx->cell_lattice ) {
-		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
-		return;
-	}
-	
-	if ( ctx->images->n_images == 0 ) {
-		displaywindow_error("There are no images to refine against.", ctx->dw);
-		return;
-	}
-		
-	ImageFeature *fitted;
-	
-	fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->dw->cur_image]);
-	
 	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
 	reproject_lattice_changed(ctx);
 	displaywindow_update(ctx->dw);
diff --git a/src/refine.h b/src/refine.h
index b15e56e..d3ca6d9 100644
--- a/src/refine.h
+++ b/src/refine.h
@@ -3,7 +3,7 @@
  *
  * Refine the reconstruction
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -17,12 +17,9 @@
 #endif
 
 #include "control.h"
-#include "image.h"
-#include "reflections.h"
 
-extern ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image);
-extern void refine_do_stack(ControlContext *ctx);
-extern void refine_do_image(ControlContext *ctx);
+extern void refine_do_sequence(ControlContext *ctx);
+extern void refine_do_cell(ControlContext *ctx);
 
 #endif	/* REFINE_H */