Refactor image handling code

Remove itrans-lsq


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

1 files changed, 99 insertions, 11 deletions

diff --git a/src/mapping.c b/src/mapping.c
index c2fed73..a4ecd96 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -12,29 +12,117 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 
 #include "reflections.h"
 #include "control.h"
 #include "itrans.h"
+#include "image.h"
+
+static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
+
+	/* "Input" space */
+	double d, x, y;
+	ImageRecord *imagerecord;
+	
+	/* Angular description of reflection */
+	double theta, psi, k;
+	
+	/* Reciprocal space */
+	double tilt;
+	double omega;
+	
+	double x_temp, y_temp, z_temp;
+	double nx, ny, nz;
+	
+	imagerecord = refl->parent;
+	x = refl->x;	y = refl->y;
+	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
+	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
+	k = 1/imagerecord->lambda;
+	
+	/* Calculate an angular description of the reflection */
+	if ( imagerecord->fmode == FORMULATION_CLEN ) {
+		x /= imagerecord->resolution;
+		y /= imagerecord->resolution;	/* Convert pixels to metres */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, imagerecord->camera_len);
+	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
+		x *= imagerecord->pixel_size;
+		y *= imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, k);
+	} else {
+		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
+		return -1;
+	}
+	psi = atan2(y, x);
+	
+	x_temp = k*sin(theta)*cos(psi);
+	y_temp = k*sin(theta)*sin(psi);
+	z_temp = k - k*cos(theta);
+	
+	/* Apply the rotations...
+		First: rotate image clockwise until tilt axis is aligned horizontally. */
+	nx = x_temp*cos(omega) + y_temp*sin(omega);
+	ny = -x_temp*sin(omega) + y_temp*cos(omega);
+	nz = z_temp;
+	
+	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
+		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
+	x_temp = nx; y_temp = ny; z_temp = nz;	
+	nx = x_temp;
+	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
+	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
+	
+	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
+	x_temp = nx; y_temp = ny; z_temp = nz;	
+	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
+	nz = z_temp;
+	
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
+	
+	return 0;
+
+}
 
 ReflectionList *mapping_create(ControlContext *ctx) {
 
-	ReflectionList *reflectionlist;
 	int i;
 	
-	/* Create reflection context */
-	reflectionlist = reflectionlist_new();
+	/* Pass all images through itrans */
+	printf("MP: Analysing images..."); fflush(stdout);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
+	}
+	printf("done.\n");
+	
+	/* Create reflection list for measured reflections */
+	ctx->reflectionlist = reflectionlist_new();
+	printf("MP: Mapping to 3D..."); fflush(stdout);
+	for ( i=0; i<ctx->images->n_images; i++ ) {
 	
-	/* Pass all images through itrans
-	 *   (let itrans add the reflections to reflectionlist for now) */
-	printf("MP: Processing images..."); fflush(stdout);
-	ctx->reflectionlist = reflectionlist;
-	for ( i=0; i<ctx->n_images; i++ ) {
-		itrans_process_image(&ctx->images[i], ctx);
+		int j;
+		
+		/* Iterate over the features in this image */
+		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
+		
+			double nx, ny, nz, twotheta;
+			
+			if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
+				reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
+			}
+
+		}
+		
 	}
-	printf("done: %i 'measured' reflections\n", ctx->reflectionlist->n_reflections);
+	printf("done.\n");
 	
-	return reflectionlist;
+	return ctx->reflectionlist;
 
 }