Break off cell delta routine

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

1 files changed, 223 insertions, 91 deletions

diff --git a/src/refine.c b/src/refine.c
index 12f0851..4ae7ec8 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -57,12 +57,113 @@ static const char *refine_decode(RefinementIndex i) {
 }
 #endif
 
+/* Return the "cell delta" to add to "cell" in order to make "feature" fit its partner */
+static Basis refine_cell_delta(Basis *cell, ImageFeature *feature, RefinementIndex *sharedr) {
+
+	double dix, diy;
+	double dx, dy, dz, twotheta;
+	double old_x, old_y;
+	RefinementIndex index, delta, shared;
+	signed int h, k, l;
+	double a11, a12, a13, a21, a22, a23, a31, a32, a33, det;
+	double dh, dk, dl;
+	Basis cd;
+	
+	h = feature->reflection->h;
+	k = feature->reflection->k;
+	l = feature->reflection->l;
+	
+	cd.a.x = 0.0;	cd.a.y = 0.0;	cd.a.z = 0.0;
+	cd.b.x = 0.0;	cd.b.y = 0.0;	cd.b.z = 0.0;
+	cd.c.x = 0.0;	cd.c.y = 0.0;	cd.c.z = 0.0;
+	
+	/* Determine the difference vector */
+	dix = feature->partner->x - feature->x;
+	diy = feature->partner->y - feature->y;
+//	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 = feature->partner->x;
+	old_y = feature->partner->y;
+	feature->partner->x = dix + feature->partner->parent->x_centre;
+	feature->partner->y = diy + feature->partner->parent->y_centre;
+	mapping_map_to_space(feature->partner, &dx, &dy, &dz, &twotheta);
+	feature->partner->x = old_x;
+	feature->partner->y = old_y;
+//	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;
+	if ( h ) index |= INDEX_A;
+	if ( k ) index |= INDEX_B;
+	if ( l ) index |= INDEX_C;
+	assert(index != 0);	/* Can't refine using the central beam! */
+	
+	/* Set up the coordinate transform from hkl to xyz */
+	a11 = cell->a.x;  a12 = cell->a.y;  a13 = cell->a.z;
+	a21 = cell->b.x;  a22 = cell->b.y;  a23 = cell->b.z;
+	a31 = cell->c.x;  a32 = cell->c.y;  a33 = cell->c.z;
+	
+	/* Invert the matrix to get dh,dk,dl from dx,dy,dz */
+	det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+	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("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;
+	if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B;
+	if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C;
+	
+	shared = index & delta;
+	
+//	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("RF: Pure shear.\n");
+	}
+	
+	if ( shared & INDEX_A ) {
+		double w = (double)abs(h) / (abs(h)+abs(k)+abs(l));
+	//	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;
+	}
+	if ( shared & INDEX_B ) {
+		double w = (double)abs(k) / (abs(h)+abs(k)+abs(l));
+	//	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;
+	}
+	if ( shared & INDEX_C ) {
+		double w = (double)abs(l) / (abs(h)+abs(k)+abs(l));
+	//	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;
+	}
+	
+//	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);
+	
+	*sharedr = shared;
+	
+	return cd;
+
+}
+
 /* Use the IPR algorithm to make "cell" fit the given image */
 static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) {
 	
 	ImageFeatureList *flist;
 	int i;
-	Basis cd; /* Cell delta */
+	Basis cd; /* Overall cell delta */
 	int n_a = 0;
 	int n_b = 0;
 	int n_c = 0;
@@ -79,98 +180,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 	
 	for ( i=0; i<image->rflist->n_features; i++ ) {
 		
-		double dix, diy;
-		double dx, dy, dz, twotheta;
-		double old_x, old_y;
-		RefinementIndex index, delta, shared;
-		signed int h, k, l;
-		double a11, a12, a13, a21, a22, a23, a31, a32, a33, det;
-		double dh, dk, dl;
+		Basis this_cd;
+		RefinementIndex shared;
 		
 		/* Skip if no partner */
 		if ( !image->rflist->features[i].partner ) continue;
 		
-		h = image->rflist->features[i].reflection->h;
-		k = image->rflist->features[i].reflection->k;
-		l = image->rflist->features[i].reflection->l;
-		
-		/* Determine the difference vector */
-		dix = image->rflist->features[i].partner->x - image->rflist->features[i].x;
-		diy = image->rflist->features[i].partner->y - image->rflist->features[i].y;
-	//	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 = image->rflist->features[i].partner->x;
-		old_y = image->rflist->features[i].partner->y;
-		image->rflist->features[i].partner->x = dix + image->rflist->features[i].partner->parent->x_centre;
-		image->rflist->features[i].partner->y = diy + image->rflist->features[i].partner->parent->y_centre;
-		mapping_map_to_space(image->rflist->features[i].partner, &dx, &dy, &dz, &twotheta);
-		image->rflist->features[i].partner->x = old_x;
-		image->rflist->features[i].partner->y = old_y;
-	//	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;
-		if ( h ) index |= INDEX_A;
-		if ( k ) index |= INDEX_B;
-		if ( l ) index |= INDEX_C;
-		assert(index != 0);	/* Can't refine using the central beam! */
-		
-		/* Set up the coordinate transform from hkl to xyz */
-		a11 = cell->a.x;  a12 = cell->a.y;  a13 = cell->a.z;
-		a21 = cell->b.x;  a22 = cell->b.y;  a23 = cell->b.z;
-		a31 = cell->c.x;  a32 = cell->c.y;  a33 = cell->c.z;
-		
-		/* Invert the matrix to get dh,dk,dl from dx,dy,dz */
-		det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
-		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("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;
-		if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B;
-		if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C;
-		
-		shared = index & delta;
-		
-	//	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("RF: Pure shear.\n");
-		}
-		
-		if ( shared & INDEX_A ) {
-			double w = (double)abs(h) / (abs(h)+abs(k)+abs(l));
-		//	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;
-			n_a++;
-		}
-		if ( shared & INDEX_B ) {
-			double w = (double)abs(k) / (abs(h)+abs(k)+abs(l));
-		//	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;
-			n_b++;
-		}
-		if ( shared & INDEX_C ) {
-			double w = (double)abs(l) / (abs(h)+abs(k)+abs(l));
-		//	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("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);
+		this_cd = refine_cell_delta(cell, &image->rflist->features[i], &shared);
+		cd = basis_add(cd, this_cd);
+		if ( shared & INDEX_A ) n_a++;
+		if ( shared & INDEX_B ) n_b++;
+		if ( shared & INDEX_C ) n_c++;
 		
 		done = TRUE;
 		
@@ -232,6 +252,112 @@ static gint refine_step(GtkWidget *step_button, ControlContext *ctx) {
 
 }
 
+static int refine_sequence_random(ControlContext *ctx, double *fit, double *warp) {
+
+	int i;
+	
+	for ( i=0; i<1000; i++ ) {
+	
+		int ni;
+		int nf;
+		ImageRecord *image;
+		ImageFeature *feature;
+		
+		/* Ensure lattice is up to date */
+		reproject_lattice_changed(ctx);
+		ctx->images->images[i].rflist = NULL;	/* Invalidate reprojection for this image - it's wrong */
+		
+		/* Select a random feature from a random image */
+		ni = ((float)ctx->images->n_images * random()) / RAND_MAX;
+		image = &ctx->images->images[ni];
+		nf = ((float)image->features->n_features * random()) / RAND_MAX;
+		feature = &image->features->features[nf];
+		
+		
+		
+	}
+	
+	*fit = 0.0;
+	*warp = 0.0;
+	
+	return 0;
+
+}
+
+static gint refine_random_sequence(GtkWidget *widget, 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;
+	ImageDisplay *id;
+	
+	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;
+	}
+	
+	/* Temporarily disable ImageDisplay stuff */
+	id = ctx->reproject_id;
+	ctx->reproject_id = NULL;
+	
+	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_random(ctx, &fit, &warp) ) {
+			printf("RF: Sequencer run 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));
+		
+	}
+	
+	ctx->reproject_id = id;
+	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);
+	
+	return 0;
+
+}
+
 static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) {
 
 	int i;
@@ -298,7 +424,7 @@ static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp)
 
 }
 
-static gint refine_sequence(GtkWidget *step_button, ControlContext *ctx) {
+static gint refine_sequence(GtkWidget *widget, ControlContext *ctx) {
 
 	double omega_offs;
 	GtkWidget *window_fit;
@@ -389,6 +515,7 @@ void refine_open(ControlContext *ctx) {
 	GtkWidget *label;
 	GtkWidget *step_button;
 	GtkWidget *sequence_button;
+	GtkWidget *sequence_random_button;
 	
 	if ( ctx->refine_window ) return;
 	
@@ -418,10 +545,15 @@ void refine_open(ControlContext *ctx) {
 	gtk_misc_set_alignment(GTK_MISC(label), 0.0, 0.5);
 	gtk_label_set_markup(GTK_LABEL(label), "<span weight=\"bold\">Sequencing</span>");
 	gtk_table_attach_defaults(GTK_TABLE(table), label, 1, 2, 4, 5);
-	sequence_button = gtk_button_new_with_label("Run Sequencer");
+
+	sequence_button = gtk_button_new_with_label("Run Sweeping Sequencer");
 	gtk_table_attach_defaults(GTK_TABLE(table), sequence_button, 1, 2, 5, 6);
 	g_signal_connect(G_OBJECT(sequence_button), "clicked", G_CALLBACK(refine_sequence), ctx);
 		
+	sequence_random_button = gtk_button_new_with_label("Run Random Sequencer");
+	gtk_table_attach_defaults(GTK_TABLE(table), sequence_random_button, 1, 2, 6, 7);
+	g_signal_connect(G_OBJECT(sequence_random_button), "clicked", G_CALLBACK(refine_random_sequence), ctx);
+	
 	g_signal_connect(G_OBJECT(ctx->refine_window), "response", G_CALLBACK(refine_response), ctx);
 	gtk_widget_show_all(ctx->refine_window);