Initial refinement maths

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

9 files changed, 359 insertions, 8 deletions

diff --git a/Makefile.am b/Makefile.am
index 218699f..a0214f9 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,5 +1,5 @@
 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/cache.h src/itrans-threshold.h src/basis.h \
 		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/dirax.h src/image.h src/refine.h src/gtk-valuegraph.h
 SUBDIRS = src data
diff --git a/src/Makefile.am b/src/Makefile.am
index 78f8943..c9d4dd9 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,7 +1,7 @@
 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
+		dirax.c image.c refine.c gtk-valuegraph.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
diff --git a/src/gtk-valuegraph.c b/src/gtk-valuegraph.c
new file mode 100644
index 0000000..a847946
--- /dev/null
+++ b/src/gtk-valuegraph.c
@@ -0,0 +1,235 @@
+/*
+ * gtk-valuegraph.c
+ *
+ * A widget to display a graph of a sequence of values
+ *
+ * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  synth2d - two-dimensional Fourier synthesis
+ *
+ */
+
+#include <gtk/gtk.h>
+#include <math.h>
+
+#include "gtk-valuegraph.h"
+
+static GtkObjectClass *parent_class = NULL;
+
+static void gtk_value_graph_destroy(GtkObject *gtk_value_graph) {
+	parent_class->destroy(gtk_value_graph);
+}
+
+GtkWidget *gtk_value_graph_new() {
+
+	GtkValueGraph *gtk_value_graph;
+	
+	gtk_value_graph = GTK_VALUE_GRAPH(gtk_type_new(gtk_value_graph_get_type()));
+	gtk_value_graph->data = NULL;
+	gtk_value_graph->n = 0;
+	
+	return GTK_WIDGET(gtk_value_graph);
+
+}
+
+static GObject *gtk_value_graph_constructor(GType type, guint n_construct_properties, GObjectConstructParam *construct_properties) {
+
+	GtkValueGraphClass *class;
+	GObjectClass *p_class;
+	GObject *obj;
+	
+	class = GTK_VALUE_GRAPH_CLASS(g_type_class_peek(gtk_value_graph_get_type()));
+	p_class = G_OBJECT_CLASS(g_type_class_peek_parent(class));
+	
+	obj = p_class->constructor(type, n_construct_properties, construct_properties);
+	
+	return obj;
+
+}
+
+static void gtk_value_graph_class_init(GtkValueGraphClass *class) {
+
+	GtkObjectClass *object_class;
+	GObjectClass *g_object_class;
+
+	object_class = (GtkObjectClass *) class;
+	g_object_class = G_OBJECT_CLASS(class);
+	
+	object_class->destroy = gtk_value_graph_destroy;
+	g_object_class->constructor = gtk_value_graph_constructor;
+
+	parent_class = gtk_type_class(gtk_drawing_area_get_type());
+	
+}
+
+
+static gint gtk_value_graph_draw(GtkWidget *graph, GdkEventExpose *event, gpointer data) {
+
+	GtkValueGraph *vg;
+	unsigned int bw_left, bw_right, bw_top, bw_bottom;
+	PangoLayout *y0_layout;
+	PangoLayout *y1_layout;
+	PangoLayout *x0_layout;
+	PangoLayout *x1_layout;
+	PangoRectangle y0_extent, y1_extent, x0_extent, x1_extent;
+	unsigned int width, height;
+	char tmp[32];
+	unsigned int i;
+	
+	vg = GTK_VALUE_GRAPH(graph);
+	
+	/* Blank white background */
+	gdk_draw_rectangle(graph->window, graph->style->white_gc, TRUE, 0, 0, graph->allocation.width, graph->allocation.height);
+	
+	/* Create PangoLayouts for labels */
+	y0_layout = gtk_widget_create_pango_layout(graph, "0");
+	pango_layout_get_pixel_extents(y0_layout, NULL, &y0_extent);
+
+	if ( fabs(log(vg->ymax)/log(10)) < 3 ) {
+		snprintf(tmp, 31, "%.4f", vg->ymax);
+	} else {
+		snprintf(tmp, 31, "%1.1e", vg->ymax);
+	}
+	y1_layout = gtk_widget_create_pango_layout(graph, tmp);
+	pango_layout_get_pixel_extents(y1_layout, NULL, &y1_extent);
+	
+	x0_layout = gtk_widget_create_pango_layout(graph, "0");
+	pango_layout_get_pixel_extents(x0_layout, NULL, &x0_extent);
+
+	if ( vg->xmax < 1000 ) {
+		snprintf(tmp, 31, "%i", vg->xmax);
+	} else {
+		snprintf(tmp, 31, "%1.1e", (double)vg->xmax);
+	}
+	x1_layout = gtk_widget_create_pango_layout(graph, tmp);
+	pango_layout_get_pixel_extents(x1_layout, NULL, &x1_extent);
+	
+	/* Determine border widths */
+	bw_left = 1+((y1_extent.width > y0_extent.width) ? y1_extent.width : y0_extent.width);
+	bw_right = 1+x1_extent.width/2;
+	bw_top = 1+y1_extent.height/2;
+	bw_bottom = 1+((x1_extent.height > x0_extent.height) ? x1_extent.height : x0_extent.height);
+	width = graph->allocation.width;
+	height = graph->allocation.height;
+	
+	/* Draw axis lines */
+	gdk_draw_line(graph->window, graph->style->black_gc, bw_left, height-1-bw_bottom, bw_left, bw_top);
+	gdk_draw_line(graph->window, graph->style->black_gc, bw_left, height-1-bw_bottom, width-1-bw_right, height-1-bw_bottom);
+
+	/* Label axes */
+	gdk_draw_layout(graph->window, graph->style->black_gc, 1+bw_left-x0_extent.width/2, height-1-bw_bottom, x0_layout);
+	gdk_draw_layout(graph->window, graph->style->black_gc, width-bw_right-x1_extent.width/2, height-1-bw_bottom, x1_layout);
+	gdk_draw_layout(graph->window, graph->style->black_gc, bw_left-y0_extent.width-1, height-1-bw_bottom-y0_extent.height/2, y0_layout);
+	gdk_draw_layout(graph->window, graph->style->black_gc, bw_left-y1_extent.width-1, 1, y1_layout);
+	
+	/* Plot data */
+	for ( i=0; i<vg->n; i++ ) {
+		
+		unsigned int x, y;
+		double xd, yd;
+		
+		xd = (((double)width-bw_left-bw_right)/(double)vg->xmax)*(double)(i+1);	/* Graph axes go from 1 */
+		x = bw_left + xd;
+		yd = (((double)height-bw_top-bw_bottom)/(double)vg->ymax)*(double)vg->data[i];
+		y = height-bw_bottom - yd;
+		
+		gdk_draw_point(graph->window, graph->style->black_gc, x, y);
+		
+	}
+	
+	return 0;
+
+}
+
+static void gtk_value_graph_init(GtkValueGraph *gtk_value_graph) {
+	gtk_widget_set_size_request(GTK_WIDGET(gtk_value_graph), 100, 200);
+	g_signal_connect(G_OBJECT(gtk_value_graph), "expose_event", G_CALLBACK(gtk_value_graph_draw), NULL);
+}
+
+guint gtk_value_graph_get_type(void) {
+
+  	static guint gtk_value_graph_type = 0;
+
+	if ( !gtk_value_graph_type ) {
+
+		GtkTypeInfo gtk_value_graph_info = {
+			"GtkValueGraph",
+			sizeof(GtkValueGraph),
+			sizeof(GtkValueGraphClass),
+			(GtkClassInitFunc) gtk_value_graph_class_init,
+			(GtkObjectInitFunc) gtk_value_graph_init,
+			NULL,
+			NULL,
+			(GtkClassInitFunc) NULL,
+		};
+		gtk_value_graph_type = gtk_type_unique(gtk_drawing_area_get_type(), &gtk_value_graph_info);
+
+	}
+
+	return gtk_value_graph_type;
+
+}
+
+static double gtk_value_graph_peak(double *data, unsigned int n) {
+
+	unsigned int i;
+	double max;
+	
+	if ( n == 0 ) return 1;
+	
+	max = 0;
+	for ( i=0; i<n; i++ ) {
+		if ( data[i] > max ) max = data[i];
+	}
+
+	return max;
+
+}
+
+/* Calculate the best range for a axis with maximum value n */
+static double gtk_value_graph_axis_max(double n) {
+
+	double mantissa, exponent, test;
+return n;	
+	if ( n == 0 ) return 1;
+	
+	/* Convert to standard form */
+	exponent = rint(log(n)/log(10));
+	mantissa = n / pow(10, exponent);
+
+	/* Check if the value can be exactly represented */
+	test = mantissa * 10;
+	test = rint(test);
+	test /= 10;
+	if ( fabs(test - mantissa) > 0.001 ) {
+	
+		/* Round the mantissa upwards */
+		mantissa += 0.1;
+		mantissa *= 10;
+		mantissa = rint(mantissa);
+		mantissa /= 10;
+	
+	} /* Else don't touch it */
+	
+	return mantissa*pow(10, exponent);
+	
+}
+
+void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) {
+	
+	double dmax;
+	
+	/* Recalculate axes */
+	dmax = gtk_value_graph_peak(data, n);
+	vg->data = data;
+	vg->n = n;
+	vg->xmax = gtk_value_graph_axis_max(n);
+	vg->ymax = gtk_value_graph_axis_max(dmax);
+	
+	//printf("n=%i, dmax=%f => xmax=%i, ymax=%f\n", n, dmax, vg->xmax, vg->ymax);
+	
+	/* Schedule redraw */
+	gtk_widget_queue_draw_area(GTK_WIDGET(vg), 0, 0, GTK_WIDGET(vg)->allocation.width, GTK_WIDGET(vg)->allocation.height);
+	
+}
+
diff --git a/src/gtk-valuegraph.h b/src/gtk-valuegraph.h
new file mode 100644
index 0000000..b976ef7
--- /dev/null
+++ b/src/gtk-valuegraph.h
@@ -0,0 +1,42 @@
+/*
+ * gtk-valuegraph.c
+ *
+ * A widget to display a graph of a sequence of values
+ *
+ * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  synth2d - two-dimensional Fourier synthesis
+ *
+ */
+
+#ifndef GTKVALUEGRAPH_H
+#define GTKVALUEGRAPH_H
+
+#include <gtk/gtk.h>
+
+typedef struct {
+
+	GtkDrawingArea parent;		/* Parent widget */
+	
+	double *data;			/* Data to be graphed */
+	unsigned int n;			/* Number of data points */
+	unsigned int xmax;		/* Maximum value on x (index) axis */
+	double ymax;			/* Maximum value on y (data) axis */
+	
+} GtkValueGraph;
+
+typedef struct {
+	GtkDrawingAreaClass parent_class;
+	void (* changed) (GtkValueGraph *gtkvaluegraph);
+} GtkValueGraphClass;
+
+extern guint gtk_value_graph_get_type(void);
+extern GtkWidget *gtk_value_graph_new(void);
+extern void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n);
+
+#define GTK_VALUE_GRAPH(obj)		GTK_CHECK_CAST(obj, gtk_value_graph_get_type(), GtkValueGraph)
+#define GTK_VALUE_GRAPH_CLASS(class)	GTK_CHECK_CLASS_CAST(class, gtk_value_graph_get_type(), GtkValueGraphClass)
+#define GTK_IS_VALUE_GRAPH(obj)		GTK_CHECK_TYPE(obj, gtk_value_graph_get_type())
+
+#endif /* GTKVALUEGRAPH_H */
+
diff --git a/src/image.c b/src/image.c
index b2b091d..00848cf 100644
--- a/src/image.c
+++ b/src/image.c
@@ -79,6 +79,7 @@ void image_add_feature(ImageFeatureList *flist, double x, double y, ImageRecord
 	flist->features[flist->n_features].intensity = intensity;
 	flist->features[flist->n_features].parent = parent;
 	flist->features[flist->n_features].partner = NULL;
+	flist->features[flist->n_features].partner_d = 0.0;
 	
 	flist->n_features++;
 
diff --git a/src/image.h b/src/image.h
index 82c6443..c87498c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -25,7 +25,8 @@ typedef struct imagefeature_struct {
 	double				y;
 	double				intensity;
 	
-	struct imagefeature_struct	*partner;
+	struct imagefeature_struct	*partner;	/* Partner for this feature (in another feature list) or NULL */
+	double				partner_d;	/* Distance between this feature and its partner, if any. */
 	
 } ImageFeature;
 
diff --git a/src/refine.c b/src/refine.c
index fba19fb..a0b00e0 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -13,8 +13,74 @@
 #include <config.h>
 #endif
 
+#include <gtk/gtk.h>
+#include <math.h>
+
 #include "displaywindow.h"
+#include "gtk-valuegraph.h"
+#include "basis.h"
+#include "reflections.h"
+#include "image.h"
+#include "reproject.h"
+
+static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_lattice) {
+
+	ImageFeatureList *rflist;
+	ImageFeatureList *flist;
+	double xc, yc;
+	int i;
+	double total;
+	
+	xc = image->x_centre;
+	yc = image->y_centre;
+	rflist = reproject_get_reflections(image, cell_lattice);
+	flist = image->features;
+	
+	reproject_partner_features(rflist, image);
+	
+	total = 0.0;
+	for ( i=0; i<rflist->n_features; i++ ) {
+		
+		double d;
+		
+		d = rflist->features[i].partner_d;
+		total += d*d;
+	
+	}
+	
+	image_feature_list_free(rflist);
+	
+	return sqrt(total);
+
+}
 
 void refine_open(DisplayWindow *dw) {
+
+	GtkWidget *window;
+	GtkWidget *graph;
+	double old_tilt;
+	int n;
+	double values[401];
+	size_t idx;
+	double tilt;
+	
+	window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
+	graph = gtk_value_graph_new();
+	
+	dw->ctx->cell_lattice = reflection_list_from_cell(dw->ctx->cell);
+	n = 0;
+	idx = 0;
+	old_tilt = dw->ctx->images->images[n].tilt;
+	for ( tilt=old_tilt-0.2; tilt<=old_tilt+0.2; tilt+=0.001 ) {
+		dw->ctx->images->images[n].tilt = tilt;
+		values[idx++] = refine_image_deviation(&dw->ctx->images->images[n], dw->ctx->cell_lattice);
+		printf("%f %f\n", tilt, values[idx-1]); fflush(stdout);
+	}
+	dw->ctx->images->images[n].tilt = old_tilt;
+	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph), values, 41);
+	
+	gtk_container_add(GTK_CONTAINER(window), graph);
+	gtk_widget_show_all(window);
+
 }
 
diff --git a/src/reproject.c b/src/reproject.c
index 1726639..e919523 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -19,7 +19,7 @@
 #include "displaywindow.h"
 #include "image.h"
 
-ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) {
+ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
 
 	ImageFeatureList *flist;
 	Reflection *reflection;
@@ -178,7 +178,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 }
 
 /* Attempt to find the partner for "feature" from the feature list of "image" */
-static ImageFeature *reproject_find_partner(ImageFeature *feature, ImageRecord *image) {
+static ImageFeature *reproject_find_partner(ImageFeature *feature, ImageRecord *image, double *d) {
 
 	int i;
 	double dmin = +HUGE_VAL;
@@ -199,6 +199,7 @@ static ImageFeature *reproject_find_partner(ImageFeature *feature, ImageRecord *
 	}
 	
 	if ( dmin <= 70.0 ) {
+		*d = dmin;
 		return &image->features->features[closest];
 	}
 	
@@ -212,7 +213,12 @@ void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image) {
 	int i;
 	
 	for ( i=0; i<flist->n_features; i++ ) {
-		flist->features[i].partner = reproject_find_partner(&flist->features[i], image);
+	
+		double d = 0.0;
+		
+		flist->features[i].partner = reproject_find_partner(&flist->features[i], image, &d);
+		flist->features[i].partner_d = d;
+	
 	}
 
 }
@@ -228,7 +234,7 @@ static void reproject_mark_peaks(ControlContext *ctx) {
 	yc = ctx->images->images[ctx->reproject_cur_image].y_centre;
 	
 	/* Draw the reprojected peaks */
-	rflist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	rflist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
 	for ( j=0; j<rflist->n_features; j++ ) {
 		imagedisplay_add_mark(ctx->reproject_id, xc+rflist->features[j].x, yc+rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
 	}
diff --git a/src/reproject.h b/src/reproject.h
index 39cd64d..a0dfd5f 100644
--- a/src/reproject.h
+++ b/src/reproject.h
@@ -19,7 +19,7 @@
 #include "control.h"
 #include "image.h"
 
-extern ImageFeatureList *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionList *reflectionlist, ControlContext *ctx);
+extern ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist);
 extern void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image);
 extern void reproject_open(ControlContext *ctx);