path: root/src/refine.c

/*
 * refine.c
 *
 * Refine the reconstruction
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <gtk/gtk.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.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"

/* 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) {

	ImageFeatureList *rflist;
	ImageFeatureList *flist;
	int i;
	double total;
	
	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);

}

typedef enum {
	INDEX_A = 1<<0,
	INDEX_B = 1<<1,
	INDEX_C = 1<<2
} RefinementIndex;

/* Use the IPR algorithm to make "cell" fit the given image */
static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) {

	ImageFeatureList *rflist;
	ImageFeatureList *flist;
	int i;
	Basis cd; /* Cell delta */
	int n_a = 0;
	int n_b = 0;
	int n_c = 0;
	
	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;
	
	rflist = reproject_get_reflections(image, cell_lattice);
	flist = image->features;
	reproject_partner_features(rflist, image);
	
	for ( i=0; i<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;
		
		/* Skip if no partner */
		if ( !rflist->features[i].partner ) continue;
		
		/* Determine the difference vector */
		dix = rflist->features[i].partner->x - rflist->features[i].x;
		diy = rflist->features[i].partner->y - rflist->features[i].y;
		
		/* Map the difference vector to the relevant tilted plane */
		old_x = rflist->features[i].partner->x;
		old_y = rflist->features[i].partner->y;
		rflist->features[i].partner->x = dix;
		rflist->features[i].partner->y = diy;
		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;
		
		h = rflist->features[i].h;
		k = rflist->features[i].k;
		l = rflist->features[i].l;
		printf("Feature %3i: %3i %3i %3i dev=%8e %8e %8e (%5f mrad)\n", i, h, k, l, 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("dev(hkl) = %f %f %f\n", dh, dk, dl);
	
		delta = 0;
		if ( fabs(dh) < 0.001 ) delta |= INDEX_A;
		if ( fabs(dk) < 0.001 ) delta |= INDEX_B;
		if ( fabs(dl) < 0.001 ) delta |= INDEX_C;
		/* Typically, 'delta' will have all three bits set */
		
		shared = index & delta;
		
		if ( shared == 0 ) {
			/* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */
			shared = index;
		}
		
		if ( shared & INDEX_A ) {
			double w = abs(h) / (abs(h)+abs(k)+abs(l));
			cd.a.x += w*dx / h;
			cd.a.y += w*dy / h;
			cd.a.z += w*dz / h;
			n_a++;
		}
		if ( shared & INDEX_B ) {
			double w = abs(k) / (abs(h)+abs(k)+abs(l));
			cd.b.x += w*dx / k;
			cd.b.y += w*dy / k;
			cd.b.z += w*dz / k;
			n_b++;
		}
		if ( shared & INDEX_C ) {
			double w = abs(l) / (abs(h)+abs(k)+abs(l));
			cd.c.x += w*dx / l;
			cd.c.y += w*dy / l;
			cd.c.z += w*dz / l;
			n_c++;
		}

	}
	
	image_feature_list_free(rflist);
	
	if ( n_a ) {
		cd.a.x /= n_a;	cd.a.y /= n_a;	cd.a.z /= n_a;
	}
	if ( n_b ) {
		cd.b.x /= n_b;	cd.b.y /= n_b;	cd.b.z /= n_b;
	}
	if ( n_c ) {
		cd.c.x /= n_c;	cd.c.y /= n_c;	cd.c.z /= 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);
	
	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;
	cell->c.x += cd.c.x;	cell->c.y += cd.c.y;	cell->c.z += cd.c.z;
	
}

/* Display a graph of root sum squared deviation distance against some other parameter */
static void refine_show_graph(ControlContext *ctx, int n) {

	GtkWidget *window;
	GtkWidget *graph;
	double old_tilt;
	double *values;
	size_t idx;
	double tilt;
	
	window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
	gtk_window_set_default_size(GTK_WINDOW(window), 640, 256);
	gtk_window_set_title(GTK_WINDOW(window), "Refinement Graph");
	graph = gtk_value_graph_new();
	
	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
	old_tilt = ctx->images->images[n].tilt;
	values = malloc(401*sizeof(double));
	tilt = old_tilt-0.2;
	for ( idx=0; idx<401; idx++ ) {
		ctx->images->images[n].tilt = tilt;
		values[idx] = refine_image_deviation(&ctx->images->images[n], ctx->cell_lattice);
		tilt += 0.001;
	}
	ctx->images->images[n].tilt = old_tilt;
	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph), values, idx);
	
	gtk_container_add(GTK_CONTAINER(window), graph);
	gtk_widget_show_all(window);

}

static gint refine_graph(GtkWidget *step_button, ControlContext *ctx) {
	refine_show_graph(ctx, ctx->reproject_cur_image);
	return 0;
}

static gint refine_step(GtkWidget *step_button, ControlContext *ctx) {

	if ( (ctx->reproject_id) && (ctx->cell_lattice) ) {
		refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
		reproject_lattice_changed(ctx);
	} else {
		displaywindow_error("Please first open the reprojection window and select the image to fit", ctx->dw);
	}

	return 0;

}

static gint refine_response(GtkWidget *refine_window, gint response, ControlContext *ctx) {

	ctx->refine_window = NULL;
	gtk_widget_destroy(refine_window);
	
	return 0;

}

void refine_open(ControlContext *ctx) {

	GtkWidget *vbox;
	GtkWidget *hbox;
	GtkWidget *table;
	GtkWidget *label;
	GtkWidget *step_button;
	GtkWidget *graph_button;
	GtkWidget *sequence_button;
	
	if ( ctx->refine_window ) return;
	
	ctx->refine_window = gtk_dialog_new_with_buttons("Refine Reconstruction", GTK_WINDOW(ctx->dw->window),
		GTK_DIALOG_DESTROY_WITH_PARENT,	GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL);
	gtk_window_set_default_size(GTK_WINDOW(ctx->refine_window), 256, -1);
	
	vbox = gtk_vbox_new(FALSE, 0);
	hbox = gtk_hbox_new(TRUE, 0);
	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(ctx->refine_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5);
	
	table = gtk_table_new(5, 1, FALSE);
	gtk_table_set_row_spacings(GTK_TABLE(table), 5);
	gtk_table_set_col_spacings(GTK_TABLE(table), 5);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), FALSE, FALSE, 0);
	
	label = gtk_label_new("Steps");
	gtk_label_set_markup(GTK_LABEL(label), "<span weight=\"bold\">Steps</span>");
	gtk_misc_set_alignment(GTK_MISC(label), 0.0, 0.5);
	gtk_table_attach_defaults(GTK_TABLE(table), label, 1, 2, 1, 2);
	step_button = gtk_button_new_with_label("Refine Lattice to Fit Current Pattern");
	gtk_table_attach_defaults(GTK_TABLE(table), step_button, 1, 2, 2, 3);
	g_signal_connect(G_OBJECT(step_button), "clicked", G_CALLBACK(refine_step), ctx);
	
	graph_button = gtk_button_new_with_label("Show Graph of Deviation Against Parameter");
	gtk_table_attach_defaults(GTK_TABLE(table), graph_button, 1, 2, 3, 4);
	g_signal_connect(G_OBJECT(graph_button), "clicked", G_CALLBACK(refine_graph), ctx);
	
	label = gtk_label_new("Sequencing");
	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");
	gtk_table_attach_defaults(GTK_TABLE(table), sequence_button, 1, 2, 5, 6);
		
	g_signal_connect(G_OBJECT(ctx->refine_window), "response", G_CALLBACK(refine_response), ctx);
	gtk_widget_show_all(ctx->refine_window);
	

}