Refinement sequencing

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

1 files changed, 142 insertions, 15 deletions

diff --git a/src/refine.c b/src/refine.c
index b2dd5d0..85573b6 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -18,6 +18,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <assert.h>
+#include <string.h>
 
 #include "displaywindow.h"
 #include "gtk-valuegraph.h"
@@ -28,6 +29,7 @@
 #include "control.h"
 #include "mapping.h"
 #include "imagedisplay.h"
+#include "utils.h"
 
 /* Return the root sum squared deviation distance for all the "reprojectable" features in an image */
 static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_lattice) {
@@ -64,6 +66,7 @@ typedef enum {
 	INDEX_C = 1<<2
 } RefinementIndex;
 
+#if 0
 static const char *refine_decode(RefinementIndex i) {
 
 	switch ( i ) {
@@ -81,6 +84,7 @@ static const char *refine_decode(RefinementIndex i) {
 	}
 
 }
+#endif
 
 /* Use the IPR algorithm to make "cell" fit the given image */
 static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) {
@@ -123,7 +127,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		/* Determine the difference vector */
 		dix = rflist->features[i].partner->x - rflist->features[i].x;
 		diy = rflist->features[i].partner->y - rflist->features[i].y;
-		printf("Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy);
+	//	printf("RF: Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy);
 		
 		/* Map the difference vector to the relevant tilted plane */
 		old_x = rflist->features[i].partner->x;
@@ -133,7 +137,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		mapping_map_to_space(rflist->features[i].partner, &dx, &dy, &dz, &twotheta);
 		rflist->features[i].partner->x = old_x;
 		rflist->features[i].partner->y = old_y;
-		printf("dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3);
+	//	printf("RF: dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3);
 		
 		/* Select the basis vectors which are allowed to be altered */
 		index = 0;
@@ -152,7 +156,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		dh = ((a22*a33-a23*a32)*dx + (a23*a31-a21*a33)*dy + (a21*a32-a22*a31)*dz) / det;
 		dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det;
 		dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det;
-		printf("dev(hkl) = %f %f %f\n", dh, dk, dl);
+	//	printf("RF: dev(hkl) = %f %f %f\n", dh, dk, dl);
 		
 		delta = 0;
 		if ( fabs(dh)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_A;
@@ -161,17 +165,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		
 		shared = index & delta;
 		
-		printf(" index: %s\n delta: %s\nshared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared));
+	//	printf("RF: index: %s\nRF: delta: %s\nRF: shared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared));
 		
 		if ( shared == 0 ) {
 			/* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */
 			shared = index;
-			printf("Pure shear.\n");
+		//	printf("RF: Pure shear.\n");
 		}
 		
 		if ( shared & INDEX_A ) {
 			double w = (double)abs(h) / (abs(h)+abs(k)+abs(l));
-			printf("w(a) = %f\n", w);
+		//	printf("RF: w(a) = %f\n", w);
 			cd.a.x += w*dx / (double)h;
 			cd.a.y += w*dy / (double)h;
 			cd.a.z += w*dz / (double)h;
@@ -179,7 +183,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		}
 		if ( shared & INDEX_B ) {
 			double w = (double)abs(k) / (abs(h)+abs(k)+abs(l));
-			printf("w(b) = %f\n", w);
+		//	printf("RF: w(b) = %f\n", w);
 			cd.b.x += w*dx / (double)k;
 			cd.b.y += w*dy / (double)k;
 			cd.b.z += w*dz / (double)k;
@@ -187,16 +191,16 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		}
 		if ( shared & INDEX_C ) {
 			double w = (double)abs(l) / (abs(h)+abs(k)+abs(l));
-			printf("w(c) = %f\n", w);
+		//	printf("RF: w(c) = %f\n", w);
 			cd.c.x += w*dx / (double)l;
 			cd.c.y += w*dy / (double)l;
 			cd.c.z += w*dz / (double)l;
 			n_c++;
 		}
 		
-		printf("Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx);
-		printf("Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy);
-		printf("Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz);
+	//	printf("RF: Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx);
+	//	printf("RF: Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy);
+	//	printf("RF: Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz);
 		
 		done = TRUE;
 		
@@ -207,7 +211,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		return NULL;
 	}
 	
-	printf("n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c);
+	//printf("RF: n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c);
 	if ( n_a ) {
 		cd.a.x /= (double)n_a;	cd.a.y /= (double)n_a;	cd.a.z /= (double)n_a;
 	}
@@ -218,9 +222,9 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		cd.c.x /= (double)n_c;	cd.c.y /= (double)n_c;	cd.c.z /= (double)n_c;
 	}
 	
-	printf("Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z);
-	printf("Total distortion(b) = %+8e %+8e %+8e\n", cd.b.x, cd.b.y, cd.b.z);
-	printf("Total distortion(c) = %+8e %+8e %+8e\n", cd.c.x, cd.c.y, cd.c.z);
+//	printf("RF: Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z);
+//	printf("RF: Total distortion(b) = %+8e %+8e %+8e\n", cd.b.x, cd.b.y, cd.b.z);
+//	printf("RF: Total distortion(c) = %+8e %+8e %+8e\n", cd.c.x, cd.c.y, cd.c.z);
 	
 	cell->a.x += cd.a.x;	cell->a.y += cd.a.y;	cell->a.z += cd.a.z;
 	cell->b.x += cd.b.x;	cell->b.y += cd.b.y;	cell->b.z += cd.b.z;
@@ -291,14 +295,137 @@ static gint refine_step(GtkWidget *step_button, ControlContext *ctx) {
 
 }
 
+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;
+	
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+
+		/* Ensure lattice is up to date */
+		reproject_lattice_changed(ctx);
+	
+		if ( is_odd(i) ) {
+			
+			/* Odd-numbered images: measure */
+			ImageFeatureList *rflist;
+			ImageFeatureList *flist;
+			ImageRecord *image;
+			int j, n;
+			double image_dev_mean = 0;
+			
+			image = &ctx->images->images[i];
+			rflist = reproject_get_reflections(image, ctx->cell_lattice);
+			flist = image->features;
+			reproject_partner_features(rflist, image);
+			
+			n = 0;
+			for ( j=0; j<rflist->n_features; j++ ) {
+				
+				double dix, diy;
+				
+				/* Skip if no partner */
+				if ( !rflist->features[j].partner ) continue;
+				
+				/* Determine the difference vector */
+				dix = rflist->features[j].partner->x - rflist->features[j].x;
+				diy = rflist->features[j].partner->y - 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++;
+			
+			image_feature_list_free(rflist);
+
+		} else {
+
+			/* Even-numbered images: fit */
+			refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
+			displaywindow_update(ctx->dw);
+
+		}
+
+	}
+	
+	series_dev_mean /= series_dev_n;
+	*fit = series_dev_mean;
+	*warp = (series_dev_max - series_dev_min)/series_dev_min;
+	
+	return 0;
+
+}
+
 static gint refine_sequence(GtkWidget *step_button, 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 ) {
 		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
 		return 0;
 	}
 	
+	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");
+			return 0;
+		}
+		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));
+
+	}
 	
+	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);
 	
 	return 0;