path: root/src/reproject.c

/*
 * reproject.c
 *
 * Synthesize diffraction patterns
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

#include <stdlib.h>
#include <math.h>

#include "control.h"
#include "reflections.h"
#include "utils.h"
#include "imagedisplay.h"
#include "displaywindow.h"
#include "image.h"

ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {

	ImageFeatureList *flist;
	Reflection *reflection;
	double smax = 0.1e9;
	double tilt, omega;
	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */
	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
	//int done_debug = 0;
			
	flist = image_feature_list_new();
	
	tilt = deg2rad(image->tilt);
	omega = deg2rad(image->omega);
	
	/* Calculate the (normalised) incident electron wavevector
	 *	n is rotated "into" the reconstruction, so only one omega step. */
	nxt = 0.0;  nyt = 0.0;  nzt = 1.0;
	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
	nxt = nx;  nyt = ny;  nzt = nz;
	nx = nxt*cos(-omega) + nyt*sin(-omega);  ny = -nxt*sin(-omega) + nyt*cos(-omega);  nz = nzt;
	kx = nx / image->lambda;
	ky = ny / image->lambda;
	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
	//reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
	
	/* Determine where "up" is
	 *	See above. */
	uxt = 0.0;  uyt = 1.0;  uzt = 0.0;
	ux = uxt;  uy = cos(tilt)*uyt + sin(tilt)*uzt;  uz = -sin(tilt)*uyt + cos(tilt)*uzt;
	uxt = ux;  uyt = uy;  uzt = uz;
	ux = uxt*cos(-omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
	//reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
	
	reflection = reflectionlist->reflections;
	do {
	
		double xl, yl, zl;
		double a, b, c;
		double A1, A2, s1, s2, s;
		
		/* Get the coordinates of the reciprocal lattice point */
		xl = reflection->x;
		yl = reflection->y;
		zl = reflection->z;
		
		/* Next, solve the relrod equation to calculate the excitation error */
		a = 1.0;
		b = -2.0*(xl*nx + yl*ny + zl*nz);
		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
		s1 = 1.0/image->lambda - A1;
		s2 = 1.0/image->lambda - A2;
		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
		
		/* Skip this reflection if s is large */
		if ( fabs(s) <= smax ) {
		
			double xi, yi, zi;
			double gx, gy, gz;
			double cx, cy, cz;
			double theta;
			double x, y;
			double rx, ry, rz;
			
			/* Determine the intersection point */
			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
			
			/* Calculate Bragg angle */
			gx = xi - kx;
			gy = yi - ky;
			gz = zi - kz;	/* This is the vector from the centre of the sphere to the intersection */
			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
			
			if ( theta > 0.0 ) {
			
				double psi, disc;
				
				//reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED);
				//reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER);
				
				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
				cx = yi*nz-zi*ny;  cy = nx*zi-nz*xi;  cz = ny*xi-nx*yi;  /* c = i x n */
				psi = angle_between(cx, cy, cz, ux, uy, uz);
				
				/* Finally, resolve the +/- Pi ambiguity from the previous step */
				rx = cy*nz-cz*ny;  ry = nx*cz-nz*cx;  rz = ny*cx-nx*cy;  /* r = [i x n] x n */
				disc = angle_between(rx, ry, rz, ux, uy, uz);
			//	if ( (i==20) && !done_debug ) {
			//		reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3);
			//		reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4);
			//		reflection_add(ctx->reflectionlist, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4);
			//		printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc));
			//	}
				if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) {
					psi -= M_PI_2;		/* Case 1 */
			//		if ( (i==20) && !done_debug ) printf("case 1\n");
				} else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) {
					psi = 3*M_PI_2 - psi;	/* Case 2 */
			//		if ( (i==20) && !done_debug ) printf("case 2\n");
				} else if ( (psi < M_PI_2) && (disc < M_PI_2) ) {
					psi = 3*M_PI_2 - psi;	/* Case 3 */
			//		if ( (i==20) && !done_debug ) printf("case 3\n");
				} else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) {
					psi = 3*M_PI_2 + psi;	/* Case 4 */
			//		if ( (i==20) && !done_debug ) printf("case 4\n");
				}
				
			//	if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi));
				psi = M_PI_2 - psi;	/* Anticlockwise from "+x" instead of clockwise from "up" */
			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi));
				
				psi += omega;
			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi));
			//	if ( (i==20) && !done_debug ) done_debug = 1;
				
				/* Calculate image coordinates from polar representation */
				if ( image->fmode == FORMULATION_CLEN ) {
					x = image->camera_len*tan(theta)*cos(psi);
					y = image->camera_len*tan(theta)*sin(psi);
					x *= image->resolution;
					y *= image->resolution;
				} else if ( image->fmode == FORMULATION_PIXELSIZE ) {
					x = tan(theta)*cos(psi) / image->lambda;
					y = tan(theta)*sin(psi) / image->lambda;
					x /= image->pixel_size;
					y /= image->pixel_size;
				} else {
					fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
					return NULL;
				}
				
				x += image->x_centre;
				y += image->y_centre;
				
				/* Sanity check */
				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
				
					/* Record the reflection */
					image_add_feature_index(flist, x, y, image, reflection->intensity,
									reflection->h, reflection->k, reflection->l);
					//printf("Reflection at %f, %f\n", x, y);
				
				} /* else it's outside the picture somewhere */
				
			} /* else this is the central beam so don't worry about it */
		
		}
		
		reflection = reflection->next;
	
	} while ( reflection );
	
	return flist;

}

/* Attempt to find partners from the feature list of "image" for each feature in "flist". */
void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image) {

	int i;
	
	for ( i=0; i<flist->n_features; i++ ) {
	
		double d = 0.0;
		
		flist->features[i].partner = image_feature_closest(image->features, flist->features[i].x, flist->features[i].y, &d);
		flist->features[i].partner_d = d;
	
	}

}

static void reproject_mark_peaks(ControlContext *ctx) {

	ImageFeatureList *rflist;
	ImageFeatureList *flist;
	size_t j;
	
	/* Draw the reprojected peaks */
	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, rflist->features[j].x, rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
	}
	
	/* Now draw the original measured peaks */
	flist = ctx->images->images[ctx->reproject_cur_image].features;
	for ( j=0; j<flist->n_features; j++ ) {
		imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2);
	}
	
	/* Now connect partners */
	reproject_partner_features(rflist, &ctx->images->images[ctx->reproject_cur_image]);
	for ( j=0; j<rflist->n_features; j++ ) {
		if ( rflist->features[j].partner ) {
			imagedisplay_add_line(ctx->reproject_id, rflist->features[j].x, rflist->features[j].y,
							 rflist->features[j].partner->x, rflist->features[j].partner->y,
							 IMAGEDISPLAY_MARK_LINE_1);
		}
	}	
	
	image_feature_list_free(rflist);
	
}
	
static void reproject_update(ControlContext *ctx) {

	imagedisplay_clear_marks(ctx->reproject_id);
	//reflectionlist_clear_markers(ctx->reflectionlist);
	
	reproject_mark_peaks(ctx);
	
	imagedisplay_put_data(ctx->reproject_id, ctx->images->images[ctx->reproject_cur_image]);
	
}

static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {

	ctx->reproject_cur_image++;
	if ( ctx->reproject_cur_image == ctx->images->n_images ) ctx->reproject_cur_image = 0;
	
	reproject_update(ctx);
	
	return 0;

}

void reproject_lattice_changed(ControlContext *ctx) {

	if ( ctx->cell_lattice ) {
		reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell);
	} else {
		ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
	}
	
	if ( ctx->reproject_id ) reproject_update(ctx);

}

static gint reproject_closed(GtkWidget *widget, ControlContext *ctx) {
	ctx->reproject_id = NULL;
	return 0;
}

void reproject_open(ControlContext *ctx) {
	
	if ( ctx->reproject_id ) {
		displaywindow_error("Reprojection window is already open.", ctx->dw);
		return;
	}
	
	if ( !ctx->cell ) {
		printf("RP: No current cell\n");
		displaywindow_error("No reciprocal unit cell has been specified.", ctx->dw);
		return;
	}
	
	if ( !ctx->images ) {
		printf("RP: No images!\n");
		displaywindow_error("No images to reproject!", ctx->dw);
		return;
	}
	
	ctx->reproject_cur_image = 0;
	reproject_lattice_changed(ctx);
	
	ctx->reproject_id = imagedisplay_open_with_message(ctx->images->images[ctx->reproject_cur_image],
				"Reprojected Diffraction Pattern", "Click to change image",
				IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS, G_CALLBACK(reproject_clicked), ctx);
	
	reproject_mark_peaks(ctx);
	
	g_signal_connect(GTK_OBJECT(ctx->reproject_id->drawingarea), "destroy", G_CALLBACK(reproject_closed), ctx);
	
}