Add a new refinement test

4 files changed, 193 insertions, 2 deletions

diff --git a/Makefile.am b/Makefile.am
index 2ff295e..7ecf791 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -3,5 +3,5 @@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src
 		src/itrans-zaefferer.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h \
 		src/dirax.h src/image.h src/refine.h src/gtk-valuegraph.h src/intensities.h src/glbits.h src/gtk-valuegraph.h
 SUBDIRS = src data
-TESTS = src/refinetest2d src/refinetest3d1 src/refinetest3d2
+TESTS = src/refinetest2d src/refinetest3d1 src/refinetest3d2 src/refinetest3d3
 
diff --git a/src/.gitignore b/src/.gitignore
index 14a1da8..1fbde8f 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -3,4 +3,5 @@
 refinetest2d
 refinetest3d1
 refinetest3d2
+refinetest3d3
 dtr
diff --git a/src/Makefile.am b/src/Makefile.am
index 52986fc..ef32492 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = dtr refinetest2d refinetest3d1 refinetest3d2
+bin_PROGRAMS = dtr refinetest2d refinetest3d1 refinetest3d2 refinetest3d3
 
 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 \
@@ -14,6 +14,9 @@ refinetest3d1_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil
 refinetest3d2_SOURCES = refinetest3d2.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c control.c
 refinetest3d2_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil
 
+refinetest3d3_SOURCES = refinetest3d3.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c control.c
+refinetest3d3_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil
+
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I/usr/include/GL
 
diff --git a/src/refinetest3d3.c b/src/refinetest3d3.c
new file mode 100644
index 0000000..ffb72a5
--- /dev/null
+++ b/src/refinetest3d3.c
@@ -0,0 +1,187 @@
+/*
+ * refinetest3d3.c
+ *
+ * Unit test for refinement procedure in 3D
+ *
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#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"
+#include "control.h"
+#include "displaywindow.h"
+#include "mrc.h"
+
+static int check_cell(Basis *cell, Basis *cell_real) {
+	
+	int fail;
+	
+	printf("     Refinement    Actual\n");
+	printf("---------------------------\n");
+	printf("ax  %+8f     %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9);
+	printf("ay  %+8f     %+8f\n", cell->a.y/1e9, cell_real->a.y/1e9);
+	printf("az  %+8f     %+8f\n", cell->a.z/1e9, cell_real->a.z/1e9);
+	printf("bx  %+8f     %+8f\n", cell->b.x/1e9, cell_real->b.x/1e9);
+	printf("by  %+8f     %+8f\n", cell->b.y/1e9, cell_real->b.y/1e9);
+	printf("bz  %+8f     %+8f\n", cell->b.z/1e9, cell_real->b.z/1e9);
+	printf("cx  %+8f     %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9);
+	printf("cy  %+8f     %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9);
+	printf("cz  %+8f     %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9);
+	
+	fail = 0;
+	if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: ax not determined correctly (got %8f, should be %8f)\n",
+										cell->a.x/1e9, cell_real->a.x/1e9);
+		fail = 1;
+	} else {
+		printf("ax is OK.\n");
+	}
+	if ( fabs(cell->a.y - cell_real->a.y) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: ay not determined correctly (got %8f, should be %8f)\n",
+										cell->a.y/1e9, cell_real->a.y/1e9);
+		fail = 1;
+	} else {
+		printf("ay is OK.\n");
+	}
+	if ( fabs(cell->a.z - cell_real->a.z) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: az not determined correctly (got %8f, should be %8f)\n",
+										cell->a.z/1e9, cell_real->a.z/1e9);
+		fail = 1;
+	} else {
+		printf("az is OK.\n");
+	}
+	if ( fabs(cell->b.x - cell_real->b.x) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: bx not determined correctly (got %8f, should be %8f)\n",
+										cell->b.x/1e9, cell_real->b.x/1e9);
+		fail = 1;
+	} else {
+		printf("bx is OK.\n");
+	}
+	if ( fabs(cell->b.y - cell_real->b.y) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: by not determined correctly (got %8f, should be %8f)\n",
+										cell->b.y/1e9, cell_real->b.y/1e9);
+		fail = 1;
+	} else {
+		printf("by is OK.\n");
+	}
+	if ( fabs(cell->b.z - cell_real->b.z) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: bz not determined correctly (got %8f, should be %8f)\n",
+										cell->b.z/1e9, cell_real->b.z/1e9);
+		fail = 1;
+	} else {
+		printf("bz is OK.\n");
+	}
+	if ( fabs(cell->c.x - cell_real->c.x) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: cx not determined correctly (got %8f, should be %8f)\n",
+										cell->c.x/1e9, cell_real->c.x/1e9);
+		fail = 1;
+	} else {
+		printf("cx is OK.\n");
+	}
+	if ( fabs(cell->c.y - cell_real->c.y) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: cy not determined correctly (got %8f, should be %8f)\n",
+										cell->c.y/1e9, cell_real->c.y/1e9);
+		fail = 1;
+	} else {
+		printf("cx is OK.\n");
+	}
+	if ( fabs(cell->c.z - cell_real->c.z) > 0.01e9 ) {
+		fprintf(stderr, "refinetest3d: cz not determined correctly (got %8f, should be %8f)\n",
+										cell->c.z/1e9, cell_real->c.z/1e9);
+		fail = 1;
+	} else {
+		printf("cz is OK.\n");
+	}
+	
+	return fail;
+
+}	
+
+int main(int argc, char *argv[]) {
+
+	ControlContext *ctx;
+	ReflectionList *reflections_real;
+	Basis *cell_real;
+	int fail;
+	int i;
+	
+	ctx = control_ctx_new();
+	ctx->omega = deg2rad(45.0);
+	ctx->lambda = lambda(300.0e3);	/* 300 keV */
+	ctx->fmode = FORMULATION_PIXELSIZE;
+	ctx->x_centre = 256;
+	ctx->y_centre = 256;
+	ctx->pixel_size = 5e7;
+	for ( i=0; i<90; i++ ) {
+		image_add(ctx->images, NULL, 512, 512, deg2rad(i), ctx);
+		ctx->images->images[i].features = image_feature_list_new();
+	}
+	
+	/* Fudge to avoid horrifying pointer-related death */
+	ctx->dw = malloc(sizeof(DisplayWindow));
+	ctx->dw->cur_image = 0;
+	
+	/* The "true" cell */
+	cell_real = malloc(sizeof(Basis));
+	cell_real->a.x = 5.0e9;  cell_real->a.y = 0.0e9;  cell_real->a.z = 0.0e9;
+	cell_real->b.x = 0.0e9;  cell_real->b.y = 5.0e9;  cell_real->b.z = 0.0e9;
+	cell_real->c.x = 0.0e9;  cell_real->c.y = 0.0e9;  cell_real->c.z = 5.0e9;
+	/* The "real" reflections */
+	reflections_real = reflection_list_from_cell(cell_real);
+	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
+	ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real);
+	ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real);
+	
+	/* The "model" cell to be refined */
+	ctx->cell = malloc(sizeof(Basis));
+	ctx->cell->a.x = 5.0e9;  ctx->cell->a.y = 0.1e9;  ctx->cell->a.z = 0.1e9;
+	ctx->cell->b.x = 0.1e9;  ctx->cell->b.y = 5.0e9;  ctx->cell->b.z = 0.1e9;
+	ctx->cell->c.x = 0.1e9;  ctx->cell->c.y = 0.1e9;  ctx->cell->c.z = 5.0e9;
+	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice);
+		reproject_partner_features(ctx->images->images[i].rflist, &ctx->images->images[i]);
+	}
+	
+	refine_do_cell(ctx);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		image_feature_list_free(ctx->images->images[i].rflist);
+	}
+	
+	fail = check_cell(ctx->cell, cell_real);
+	
+	free(ctx);
+	
+	if ( fail ) return 1;
+	
+	printf("\n3D refinement test OK.\n");
+	
+	return 0;
+	
+}
+
+/* Dummy function stubs */
+#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) { };
+