Merge the two reflection_add routines for different formulations

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

8 files changed, 53 insertions, 128 deletions

diff --git a/src/control.h b/src/control.h
index bae9f4b..ac3a357 100644
--- a/src/control.h
+++ b/src/control.h
@@ -50,8 +50,8 @@ typedef enum {
 typedef struct imagerecord_struct {
 
 	uint16_t	*image;
-	double		tilt;
-	double		omega;
+	double		tilt;		/* Degrees */
+	double		omega;		/* Degrees */
 	
 	FormulationMode fmode;
 	double		pixel_size;
@@ -89,7 +89,7 @@ typedef struct cctx_struct {
 	/* Basic parameters, stored here solely so they can be copied
 	 * into the ImageRecord(s) more easily */
 	double			camera_length;
-	double			omega;
+	double			omega;		/* Degrees */
 	double			resolution;
 	double			lambda;
 	double			pixel_size;
diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c
index 168a5b3..baac6b6 100644
--- a/src/itrans-lsq.c
+++ b/src/itrans-lsq.c
@@ -257,14 +257,12 @@ unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx,
 				int dx, dy;
 				int bb_top, bb_bot, bb_lh, bb_rh;
 				double frac;
+				double brightness;
 			
-				printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", iter, x, y, ga, gb, gc, gd, ge, gf);
-				if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
-					reflection_add_from_reciprocal(ctx, (x-imagerecord.x_centre)*imagerecord.pixel_size, (y-imagerecord.y_centre)*imagerecord.pixel_size,
-										tilt_degrees, image[x + width*y]);
-				} else {
-					reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), tilt_degrees, image[x + width*y]);
-				}
+				printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n",
+						iter, x, y, ga, gb, gc, gd, ge, gf);
+				brightness = image[x + width*y];
+				reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), imagerecord, brightness);
 				n_reflections++;
 				
 				/* Remove this peak from the image */
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index f0677d4..7099c8d 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -495,7 +495,6 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx
 	gsl_matrix *p;
 	int i;
 	double px,py;
-	double tilt_degrees = imagerecord.tilt;
 	uint16_t *image = imagerecord.image;
 
 	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord.width, imagerecord.height);
@@ -505,12 +504,7 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx
 	
 		px = gsl_matrix_get(p,0,i);
 		py = gsl_matrix_get(p,1,i);
-		if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
-			reflection_add_from_reciprocal(ctx, (px-imagerecord.x_centre)*imagerecord.pixel_size, (py-imagerecord.y_centre)*imagerecord.pixel_size,
-								tilt_degrees, 1.0);
-		} else {
-			reflection_add_from_dp(ctx, (px-imagerecord.x_centre)/imagerecord.resolution, (py-imagerecord.y_centre)/imagerecord.resolution, tilt_degrees, 1.0);
-		}
+		reflection_add_from_dp(ctx, (px-imagerecord.x_centre), (py-imagerecord.y_centre), imagerecord, 1.0);
 		
 	}
 	gsl_matrix_free(m);
@@ -519,3 +513,4 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx
 	return n_reflections;
 	
 }
+
diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c
index 71c39a1..519acc8 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -26,7 +26,6 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 	int x, y;
 	unsigned int n_reflections = 0;
 	uint16_t *image = imagerecord.image;
-	double tilt_degrees = imagerecord.tilt;
 	uint16_t max = 0;
 	
 	width = imagerecord.width;
@@ -46,17 +45,8 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext
 				assert(y<height);
 				assert(x>=0);
 				assert(y>=0);
-				if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
-					reflection_add_from_reciprocal(ctx,
-								(signed)(x-imagerecord.x_centre)*imagerecord.pixel_size,
-								(signed)(y-imagerecord.y_centre)*imagerecord.pixel_size,
-								tilt_degrees, image[x + width*y]);
-				} else {
-					reflection_add_from_dp(ctx,
-								(signed)(x-imagerecord.x_centre)/imagerecord.resolution,
-								(signed)(y-imagerecord.y_centre)/imagerecord.resolution,
-								tilt_degrees, image[x + width*y]);
-				}
+				reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre),
+							imagerecord, image[x + width*y]);
 				n_reflections++;
 			}
 		}
@@ -72,7 +62,6 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr
 	int width, height;
 	unsigned int n_reflections = 0;
 	uint16_t *image;
-	double tilt_degrees = imagerecord.tilt;
 	uint16_t max;
 	int x, y;
 	
@@ -111,16 +100,8 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr
 			assert(max_y<height);
 			assert(max_x>=0);
 			assert(max_y>=0);
-			if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
-				reflection_add_from_reciprocal(ctx,
-								(signed)(max_x-imagerecord.x_centre)*imagerecord.pixel_size,
-								(signed)(max_y-imagerecord.y_centre)*imagerecord.pixel_size,
-								tilt_degrees, image[max_x + width*max_y]);
-			} else {
-				reflection_add_from_dp(ctx, (signed)(max_x-imagerecord.x_centre)/imagerecord.resolution,
-							(signed)(max_y-imagerecord.y_centre)/imagerecord.resolution,
-							tilt_degrees, image[max_x + width*max_y]);
-			}
+			reflection_add_from_dp(ctx, (max_x-imagerecord.x_centre), (max_y-imagerecord.y_centre),
+							imagerecord, image[x + width*y]);
 			n_reflections++;
 			
 			/* Remove it and its surroundings */
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index 1454fd4..c6e7b39 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -82,14 +82,8 @@ unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext
 				}
 				
 				if ( !did_something ) {
-					if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) {
-						reflection_add_from_reciprocal(ctx, (signed)(mask_x-imagerecord.x_centre)*imagerecord.pixel_size,
-										(signed)(mask_y-imagerecord.y_centre)*imagerecord.pixel_size,
-										tilt_degrees, image[mask_x + width*mask_y]);
-					} else {
-						reflection_add_from_dp(ctx, (signed)(mask_x-imagerecord.x_centre)/imagerecord.resolution, (signed)(mask_y-imagerecord.y_centre)/imagerecord.resolution,
-										tilt_degrees, image[mask_x + width*mask_y]);
-					}
+					reflection_add_from_dp(ctx, ((double)mask_x-imagerecord.x_centre), ((double)mask_y-imagerecord.y_centre),
+							imagerecord, image[mask_x + width*mask_y]);
 					n_reflections++;
 				}
 				
diff --git a/src/reflections.c b/src/reflections.c
index d53af10..a1a8267 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -93,85 +93,46 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 	
 }
 
-/* x and y in metres (in real space!) */
-void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {
+/* x and y in pixels, measured from centre of image */
+void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity) {
 
-	double nx, ny, nz;
-
-	/* Real space */
-	double L;
+	/* "Input" space */
 	double d;
-
-	/* Shared */
-	double theta;
-	double psi;
-	
-	/* Reciprocal space */
-	double r;
-	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);
-	psi = atan2(y, x);
-	
-	r = 1/ctx->lambda;
-	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);
-	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;
-	
-	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
-
-}
-
-/* Alternative formulation for when the input data is already in reciprocal space units 
-	x and y in metres^-1 (in reciprocal space) */
-void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {
-
-	double nx, ny, nz;
-	
-	/* Input (reciprocal space) */
-	double dr;
-	
-	/* Shared */
-	double theta;
-	double psi;
-	double r;
+	/* Angular description of reflection */
+	double theta, psi, k;
 	
 	/* Reciprocal space */
 	double tilt;
 	double omega;
-	double x_temp, y_temp, z_temp;
 	
-	tilt = M_PI*(tilt_degrees/180);	/* Convert to Radians */
-	omega = M_PI*(ctx->omega/180);	/* Likewise */
+	double x_temp, y_temp, z_temp;
+	double nx, ny, nz;
 	
-	dr = sqrt((x*x) + (y*y));
-	r = 1/ctx->lambda;
-	theta = atan2(dr, r);
+	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 ( ctx->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 ( ctx->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;
+	}
 	psi = atan2(y, x);
 	
-	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);
+	x_temp = k*sin(theta)*cos(psi);
+	y_temp = -k*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
+	z_temp = k- k*cos(theta);
 	
 	/* Apply the rotations...
 		First: rotate image clockwise until tilt axis is aligned horizontally. */
diff --git a/src/reflections.h b/src/reflections.h
index 571f754..ae19362 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -76,8 +76,7 @@ extern Reflection *reflection_add(ReflectionContext *reflectionctx, double x, do
 
 #include "control.h"
 
-extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity);
-extern void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity);
+extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity);
 extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r);
 
 #endif	/* REFLECTION_H */
diff --git a/src/reproject.c b/src/reproject.c
index 1deb044..faa9433 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -26,8 +26,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	double smax = 0.5e9;
 	double tilt, omega;
 	
-	tilt = 2*M_PI*(image.tilt/360);
-	omega = 2*M_PI*(image.omega/360);	/* Convert to Radians */
+	tilt = deg2rad(image.tilt);
+	omega = deg2rad(image.omega);
 	
 	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
 	
@@ -88,8 +88,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			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);
-			printf("theta=%f   ", theta);
-			theta = deg2rad(theta);
 			cx = nz*gy - ny*gz;
 			cy = nz*gx - nx*gz;
 			cz = ny*gx - nx*gy;
@@ -99,10 +97,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			ux = uyt*-sin(omega);
 			uy = uyt*cos(omega);
 			uz = uzt;
-			psi = angle_between(cx, cy, cz, ux, uy, uz);
-			psi -= 90;
-			printf("psi=%f\n", psi);
-			psi = deg2rad(psi);
+			psi = angle_between(cx, cy, cz, ux, uy, uz) - M_PI_2;
 			
 			if ( image.fmode == FORMULATION_CLEN ) {
 				x = image.camera_len*sin(theta)*cos(psi);
@@ -110,8 +105,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 				x *= image.resolution;
 				y *= image.resolution;
 			} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
-				x = 0;
-				y = 0;
+				x = sin(theta)*cos(psi) / image.lambda;
+				y = sin(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;