Implement fast reprojection

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

11 files changed, 128 insertions, 59 deletions

diff --git a/src/control.c b/src/control.c
index b626d7c..e0318c0 100644
--- a/src/control.c
+++ b/src/control.c
@@ -13,31 +13,25 @@
 
 #include "control.h"
 
-int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size) {
+int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt) {
 
 	ctx->images[ctx->n_images].tilt = tilt;
-	ctx->images[ctx->n_images].omega = omega;
-	ctx->images[ctx->n_images].pixel_size = pixel_size;
-	ctx->images[ctx->n_images].camera_len = 0;
+	ctx->images[ctx->n_images].omega = ctx->omega;
 	ctx->images[ctx->n_images].image = image;
 	ctx->images[ctx->n_images].width = width;
 	ctx->images[ctx->n_images].height = height;
+	ctx->images[ctx->n_images].lambda = ctx->lambda;
+	ctx->images[ctx->n_images].fmode = ctx->fmode;
 	
-	ctx->n_images++;
-	
-	return ctx->n_images - 1;
-
-}
-
-int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len) {
-
-	ctx->images[ctx->n_images].tilt = tilt;
-	ctx->images[ctx->n_images].omega = omega;
-	ctx->images[ctx->n_images].pixel_size = 0;
-	ctx->images[ctx->n_images].camera_len = camera_len;
-	ctx->images[ctx->n_images].image = image;
-	ctx->images[ctx->n_images].width = width;
-	ctx->images[ctx->n_images].height = height;
+	if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
+		ctx->images[ctx->n_images].pixel_size = ctx->pixel_size;
+		ctx->images[ctx->n_images].camera_len = 0;
+		ctx->images[ctx->n_images].resolution = ctx->resolution;
+	} else if ( ctx->fmode == FORMULATION_CLEN ) {
+		ctx->images[ctx->n_images].pixel_size = 0;
+		ctx->images[ctx->n_images].camera_len = ctx->camera_length;
+		ctx->images[ctx->n_images].resolution = 0;
+	}
 	
 	ctx->n_images++;
 	
diff --git a/src/control.h b/src/control.h
index 1952dab..83ba0fb 100644
--- a/src/control.h
+++ b/src/control.h
@@ -52,10 +52,17 @@ typedef struct imagerecord_struct {
 	int16_t		*image;
 	double		tilt;
 	double		omega;
+	
+	FormulationMode fmode;
 	double		pixel_size;
 	double		camera_len;
+	double		lambda;
+	double		resolution;
+	
 	int		width;
 	int		height;
+	int		x_centre;
+	int		y_centre;
 
 } ImageRecord;
 
@@ -87,12 +94,10 @@ typedef struct cctx_struct {
 	unsigned int		started;
 	unsigned int		camera_length_set;
 	unsigned int		omega_set;
-	unsigned int		max_d_set;
 	unsigned int		resolution_set;
 	unsigned int		lambda_set;
 
 	/* Miscellaneous modes and operations */
-	unsigned int		max_d;
 	unsigned int		first_image;
 	
 	/* Size of input images - assumed the same throughout. */
@@ -114,8 +119,7 @@ typedef struct cctx_struct {
 
 } ControlContext;
 
-extern int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size);
-extern int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len);
+extern int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt);
 
 #endif	/* CONTROL_H */
 
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index 8ac845b..aaa8110 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -34,12 +34,6 @@ unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, do
 			double dx1, dx2, dy1, dy2;
 			double dxs, dys;
 			double grad;
-
-			if ( ctx->max_d ) {
-				if ( (x-ctx->x_centre)*(x-ctx->x_centre) + (y-ctx->y_centre)*(y-ctx->y_centre) > ctx->max_d*ctx->max_d ) {
-					continue;
-				}
-			}
 			
 			/* Get gradients */
 			dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y];
diff --git a/src/main.c b/src/main.c
index 4b42809..4e5723a 100644
--- a/src/main.c
+++ b/src/main.c
@@ -208,7 +208,7 @@ int main(int argc, char *argv[]) {
 	#endif
 	
 	ctx->filename = strdup(argv[1]);
-	ctx->max_d = 0;
+	ctx->n_images = 0;
 	
 	if ( stat(argv[1], &stat_buffer) == -1 ) {
 		fprintf(stderr, "File '%s' not found\n", argv[1]);
diff --git a/src/mapping.c b/src/mapping.c
index 0a0dcdf..a8d50b3 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -75,6 +75,8 @@ ReflectionContext *mapping_create(ControlContext *ctx) {
 	 *   (let itrans add the reflections to rctx for now) */
 	ctx->reflectionctx = rctx;
 	for ( i=0; i<ctx->n_images; i++ ) {
+		ctx->images[i].x_centre = ctx->x_centre;
+		ctx->images[i].y_centre = ctx->y_centre;
 		itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt);
 	}
 	
diff --git a/src/mrc.c b/src/mrc.c
index 26a31de..abc096c 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -85,7 +85,7 @@ int mrc_read(ControlContext *ctx) {
 		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
 		fread(image, mrc.nx*mrc.ny*2, 1, fh);
 		
-		control_add_image_ps(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt, ctx->omega, ctx->pixel_size);
+		control_add_image(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt);
 		
 	}
 
diff --git a/src/qdrp.c b/src/qdrp.c
index 6cc9296..96adfda 100644
--- a/src/qdrp.c
+++ b/src/qdrp.c
@@ -166,26 +166,6 @@ static int qdrp_parseline(ControlContext *ctx, const char *line) {
 			
 		}
 		
-		if ( strncasecmp(line+pos, "max_d", 5) == 0 ) {
-			
-			char *max_d_s;
-			
-			skip_chars;
-			skip_whitespace;
-			mark = pos;
-			skip_chars;
-			max_d_s = strndup(line+mark, pos-mark);
-			ctx->max_d = strtod(max_d_s, NULL);
-			free(max_d_s);
-			ctx->max_d_set = 1;
-			if ( ctx->max_d == 0 ) {
-				printf("QD: max_d=0 is Not Allowed!\n");
-				return 1;
-			}
-			printf("QD: max_d = %i pixels\n", ctx->max_d);	
-			
-		}
-		
 		if ( strncasecmp(line+pos, "omega", 5) == 0 ) {
 			
 			char *omega_s;
@@ -216,10 +196,8 @@ int qdrp_read(ControlContext *ctx) {
 
 	ctx->camera_length_set = 0;
 	ctx->omega_set = 0;
-	ctx->max_d_set = 0;
 	ctx->resolution_set = 0;
 	ctx->lambda_set = 0;
-	ctx->max_d = 0;
 	ctx->started = 0;
 	ctx->first_image = 1;
 	
diff --git a/src/readpng.c b/src/readpng.c
index 68a92c3..20fc182 100644
--- a/src/readpng.c
+++ b/src/readpng.c
@@ -141,7 +141,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) {
 	png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
 	fclose(fh);
 
-	control_add_image_cl(ctx, image, width, height, tilt, ctx->omega, ctx->camera_length);
+	control_add_image(ctx, image, width, height, tilt);
 	
 	return 0;
 	
diff --git a/src/reflections.c b/src/reflections.c
index 03d901b..e16bcff 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -111,10 +111,10 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt
 	double tilt;
 	double omega;
 	double x_temp, y_temp, z_temp;
-
+	
 	tilt = 2*M_PI*(tilt_degrees/360);	/* Convert to Radians */
 	omega = 2*M_PI*(ctx->omega/360);	/* Likewise */
-
+	
 	d = sqrt((x*x) + (y*y));
 	L = ctx->camera_length;
 	theta = atan2(d, L);
@@ -124,7 +124,7 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt
 	x_temp = r*sin(theta)*cos(psi);
 	y_temp = -r*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
 	z_temp = r- r*cos(theta);
-
+	
 	/* Apply the rotations...
 		First: rotate image clockwise until tilt axis is aligned horizontally. */
 	nx = x_temp*cos(omega) + y_temp*-sin(omega);
diff --git a/src/reproject.c b/src/reproject.c
index f62e8f3..c9023af 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -10,21 +10,118 @@
  */
 
 #include <stdlib.h>
+#include <math.h>
 
 #include "control.h"
 #include "reflections.h"
 
 #define MAX_IMAGE_REFLECTIONS 1024
 
-ImageReflection *reproject_get_reflections(ImageRecord *image, size_t *n, ReflectionContext *ref) {
+ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx) {
 
 	ImageReflection *refl;
+	Reflection *reflection;
 	int i;
+	double smax = 1e9;
+	double tilt, omega;
+	
+	tilt = 2*M_PI*(image.tilt/360);
+	omega = 2*M_PI*(image.omega/360);	/* Convert to Radians */
 	
 	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
 	
 	i = 0;
 	
+	reflection = rctx->reflections;
+	do {
+	
+		double xl, yl, zl;
+		double nx, ny, nz;
+		double xt, yt, zt;
+		double a, b, c;
+		double s1, s2, s;
+		
+		/* Get the coordinates of the reciprocal lattice point */
+		xl = reflection->x;
+		yl = reflection->y;
+		zl = reflection->z;
+		
+		/* Now calculate the (normalised) incident electron wavevector */
+		xt = 0;
+		yt = - sin(tilt);
+		zt = - cos(tilt);
+		nx = xt*cos(omega) + yt*-sin(omega);
+		ny = xt*sin(omega) + yt*cos(omega);
+		nz = zt;
+		
+		/* Next, solve the reprojection equation to calculate the excitation error */
+		a = nx*nx + ny*ny + nz*nz;
+		b = 2*(xl*nx + yl*ny + zl*nz - nz/image.lambda);
+		c = xl*xl + yl*yl + zl*zl - 2*zl/image.lambda;
+		s1 = (-b + sqrt(b*b-4*a*c))/(2*a);
+		s2 = (-b - sqrt(b*b-4*a*c))/(2*a);
+		if ( s1 < s2 ) s = s1; else s = s2;
+		
+		/* Skip this reflection if s is large */
+		if ( s <= smax ) {
+		
+			double xddd, yddd, zddd;
+			double xdd, ydd, zdd;
+			double xd, yd, zd;
+			double theta, psi;
+			double x, y;
+			
+			/* Determine the intersection point */
+			xddd = xl + s*nx;  yddd = yl + s*ny;  zddd = zl + s*nz;
+			
+			/* Invert the image->3D mapping to get the image coordinates */
+			xdd = xddd;
+			ydd = (yddd/cos(tilt) - zddd*tan(tilt)/cos(tilt))/(1+tan(tilt)*tan(tilt));
+			zdd = (-zddd-ydd*sin(tilt))/cos(tilt);
+			
+			yd = (ydd-xdd*tan(omega))/(sin(omega)*tan(omega)+cos(omega));
+			xd = (xdd+yd*sin(omega))/cos(omega);
+			zd = zdd;
+			
+			if ( image.fmode == FORMULATION_CLEN ) {
+				psi = atan2(-yd, xd);
+				theta = acos(1-zd*image.lambda);
+				x = image.camera_len*sin(theta)*cos(psi);
+				y = image.camera_len*sin(theta)*sin(psi);
+				x *= image.resolution;
+			} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
+				x = xd / image.pixel_size;
+				y = yd / image.pixel_size;
+			} else {
+				fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
+				return NULL;
+			}
+			
+			/* Adjust centre */
+			x += image.x_centre;
+			y += image.y_centre;
+			
+			/* Sanity check */
+			if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) {
+			
+				/* Record the reflection */
+				refl[i].x = x;
+				refl[i].y = y;
+				i++;
+				
+				if ( i > MAX_IMAGE_REFLECTIONS ) break;
+				
+				printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
+			
+			} else {
+				fprintf(stderr, "Reflection failed sanity test\n");
+			}
+		
+		}
+		
+		reflection = reflection->next;
+
+	} while ( reflection );
 	
 	*n = i;
 	return refl;
diff --git a/src/reproject.h b/src/reproject.h
index c65549c..ac7b240 100644
--- a/src/reproject.h
+++ b/src/reproject.h
@@ -19,7 +19,7 @@
 #include "control.h"
 #include "reflections.h"
 
-extern ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *ref);
+extern ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx);
 
 #endif	/* REPROJECT_H */