path: root/src/refine.c

/*
 * refine.c
 *
 * Refine the reconstruction
 *
 * (c) 2007-2008 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 <string.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.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"
#include "imagedisplay.h"
#include "utils.h"

typedef enum {
	AXIS_X = 1,
	AXIS_Y = 2
} Axis;

static void refine_fix_unconstrained(gsl_matrix *M) {

	int row, j;
	int modified = 0;
	
	/* Find a row which is all zero */
	for ( row=0; row<M->size1; row++ ) {
	
		int zero = 1;
		
		for ( j=0; j<M->size2; j++ ) {
			if ( gsl_matrix_get(M, row, j) != 0.0 ) {
				zero = 0;
				break;
			}
		}
		
		if ( zero ) {
			
			/* Find a column which is all zero */
			int i, col;
			int row_altered = 0;
			
			for ( col=0; col<M->size2; col++ ) {
			
				int zero2 = 1;
				
				if ( row_altered ) break; /* This row is now fixed, so move on to the next */
				
				for ( i=0; i<M->size1; i++ ) {
					if ( gsl_matrix_get(M, i, col) != 0.0 ) {
						zero2 = 0;
						break;
					}
				}
				
				if ( zero2 ) {
					gsl_matrix_set(M, row, col, 1.0);
					modified = 1;
					row_altered = 1;
				}
			
			}
			
		}
		
	}
	
	if ( modified) printf("RF: M was modified to fix unconstrained variables.\n");

}

/* Use the IPR algorithm to make "cell" fit the given image */
ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
	
	ImageFeatureList *flist;
	int j, ns;
	Axis axis;
	double dax = 0.0, dbx = 0.0, dcx = 0.0;
	double day = 0.0, dby = 0.0, dcy = 0.0;
	double dlx = 0.0, dly = 0.0, dlz = 0.0;
	
	flist = image->features;
		
	for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) {
	
		gsl_permutation *perm;
		int s;
		double det;
		gsl_matrix *M;
		gsl_vector *q;
		gsl_vector *p;
		gsl_vector *scratch;
		gsl_matrix *V;
		gsl_vector *S;
		
		if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n");
		if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n");
		
		M = gsl_matrix_alloc(3, 3);
		p = gsl_vector_alloc(3);
		gsl_matrix_set_zero(M);
		gsl_vector_set_zero(p);
		
		ns = 0;
		for ( j=0; j<image->rflist->n_features; j++ ) {
		
			double val;
			ImageFeature *rf;
			signed int h, k, l;
			double xy;
			double dix, diy, dx, dy;
			double old_x, old_y;
			
			rf = &image->rflist->features[j];
			
			if ( !rf->partner ) continue;
			
			h = rf->reflection->h;
			k = rf->reflection->k;
			l = rf->reflection->l;
			
			/* Determine the difference vector */
			dix = rf->partner->x - rf->x;
			diy = rf->partner->y - rf->y;
			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy);
			
			old_x = rf->partner->x;
			old_y = rf->partner->y;
			rf->partner->x = dix + rf->partner->parent->x_centre;
			rf->partner->y = diy + rf->partner->parent->y_centre;
			mapping_scale(rf->partner, &dx, &dy);
			rf->partner->x = old_x;
			rf->partner->y = old_y;
			printf("=> %+10.5f %+10.5f nm^-1\n", dx/1e9, dy/1e9);
			
			xy = 0;
			switch ( axis ) {
				case AXIS_X : xy = dx; break;
				case AXIS_Y : xy = dy; break;
			}
			
			/* Elements of "p" */
			val = gsl_vector_get(p, 0);	val += xy * h;	gsl_vector_set(p, 0, val);
			val = gsl_vector_get(p, 1);	val += xy * k;	gsl_vector_set(p, 1, val);
			val = gsl_vector_get(p, 2);	val += xy * l;	gsl_vector_set(p, 2, val);
			gsl_matrix_get(M, 2, 2);
			/* Elements of "M" */
			val = gsl_matrix_get(M, 0, 0);	val += h * h;	gsl_matrix_set(M, 0, 0, val);
			val = gsl_matrix_get(M, 0, 1);	val += k * h;	gsl_matrix_set(M, 0, 1, val);
			val = gsl_matrix_get(M, 0, 2);	val += l * h;	gsl_matrix_set(M, 0, 2, val);
			
			val = gsl_matrix_get(M, 1, 0);	val += h * k;	gsl_matrix_set(M, 1, 0, val);
			val = gsl_matrix_get(M, 1, 1);	val += k * k;	gsl_matrix_set(M, 1, 1, val);
			val = gsl_matrix_get(M, 1, 2);	val += l * k;	gsl_matrix_set(M, 1, 2, val);
			
			val = gsl_matrix_get(M, 2, 0);	val += h * l;	gsl_matrix_set(M, 2, 0, val);
			val = gsl_matrix_get(M, 2, 1);	val += k * l;	gsl_matrix_set(M, 2, 1, val);
			val = gsl_matrix_get(M, 2, 2);	val += l * l;	gsl_matrix_set(M, 2, 2, val);
			
			ns++;
		
		}
		
		if ( ns == 0 ) {
			printf("RF: No partners found\n");
			gsl_matrix_free(M);
			gsl_vector_free(p);
			return NULL;
		}	
		
		/* Prepare for solving */
		refine_fix_unconstrained(M);
		matrix_vector_show(M, p, "RF: ");
		
		/* Calculate and display the determinant */
		perm = gsl_permutation_alloc(M->size1);
		gsl_linalg_LU_decomp(M, perm, &s);
		det = gsl_linalg_LU_det(M, s);
		printf("RF: Determinant of M = %f\n", det);
		gsl_permutation_free(perm);		
		
		/* Solve the equations */
		V = gsl_matrix_alloc(M->size2, M->size2);
		S = gsl_vector_alloc(M->size2);
		scratch = gsl_vector_alloc(M->size2);
		gsl_linalg_SV_decomp(M, V, S, scratch);
		q = gsl_vector_alloc(3);	/* This is the "answer" */
		gsl_vector_set_zero(q);
		gsl_linalg_SV_solve(M, V, S, p, q);
		gsl_vector_free(scratch);
		gsl_matrix_free(V);
		gsl_vector_free(S);
		gsl_matrix_free(M);
		gsl_vector_free(p);
		
		switch ( axis ) {
			case AXIS_X : {
				dax = gsl_vector_get(q, 0);
				dbx = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "x" of the image coordinate  */
				dcx = gsl_vector_get(q, 2);	/* system, of a,b and c */
				break;
			}
			case AXIS_Y : {
				day = gsl_vector_get(q, 0);
				dby = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "y" of the image coordinate  */
				dcy = gsl_vector_get(q, 2);	/* system, of a,b and c */
				break;
			}
		}
		
		gsl_vector_free(q);

	}
	
	printf("RF: ------------------------------------------------------------------ Refinement results ---------\n");
	printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9);
	printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9);
	printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9);
	
	/* Update the cell */
	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
	printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
	cell->a.x += dlx;  cell->a.y += dly; cell->a.z += dlz;
	
	mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
	printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
	cell->b.x += dlx;  cell->b.y += dly;  cell->b.z += dlz;
	
	mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
	printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
	
	return NULL;
	
}

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;
	
	/* Ensure that ctx->cell_lattice is set up */
	reproject_cell_to_lattice(ctx);
	
	for ( i=0; i<ctx->images->n_images; i++ ) {

		ImageRecord *image;
		int j, n;
		double image_dev_mean = 0;
		
		image = &ctx->images->images[i];
		
		/* Fit this image and update ctx->cell_lattice, index the selected pattern */
		if ( !image->rflist ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
		refine_fit_image(ctx->cell, image);
		reproject_cell_to_lattice(ctx);
		image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
		
		n = 0;
		for ( j=0; j<image->rflist->n_features; j++ ) {
			
			double dix, diy;
			
			/* Skip if no partner */
			if ( !image->rflist->features[j].partner ) continue;
			
			/* Determine the difference vector */
			dix = image->rflist->features[j].partner->x - image->rflist->features[j].x;
			diy = image->rflist->features[j].partner->y - image->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++;
		
	}
	
	series_dev_mean /= series_dev_n;
	*fit = series_dev_mean;
	*warp = (series_dev_max - series_dev_min)/series_dev_min;
	
	return 0;

}

void refine_do_stack(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_lattice ) {
		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
		return;
	}
	
	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");
		}
		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));
		
	}
	
	displaywindow_update(ctx->dw);
	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);
	
}

void refine_do_image(ControlContext *ctx) {

	if ( !ctx->cell_lattice ) {
		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
		return;
	}
	
	ImageFeature *fitted;
	
	fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->dw->cur_image]);
	
	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
	reproject_lattice_changed(ctx);
	displaywindow_update(ctx->dw);
		
}