Beginnings of a new initial basis choice (needs fixing)

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

1 files changed, 86 insertions, 16 deletions

diff --git a/src/reflections.c b/src/reflections.c
index 1a80f71..63eaa54 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -13,11 +13,14 @@
 #include <config.h>
 #endif
 
+#define _GNU_SOURCE
 #include <stdlib.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 
 #include "reflections.h"
+#include "utils.h"
 
 static void reflection_addfirst(ReflectionContext *reflectionctx) {
 	/* Create first items on lists - saves faffing later.  Corresponds to a central marker.
@@ -124,11 +127,11 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 	
 }
 
-/* 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) {
+int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
 
 	/* "Input" space */
-	double d;
+	double d, x, y;
+	ImageRecord *imagerecord;
 	
 	/* Angular description of reflection */
 	double theta, psi, k;
@@ -140,24 +143,26 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 	double x_temp, y_temp, z_temp;
 	double nx, ny, nz;
 	
-	tilt = 2*M_PI*(imagerecord.tilt/360);	/* Convert to Radians */
-	omega = 2*M_PI*(imagerecord.omega/360);	/* Likewise */
-	k = 1/imagerecord.lambda;
+	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 ( ctx->fmode == FORMULATION_CLEN ) {
-		x /= imagerecord.resolution;
-		y /= imagerecord.resolution;	/* Convert pixels to metres */
+	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 ( ctx->fmode == FORMULATION_PIXELSIZE ) {
-		x *= imagerecord.pixel_size;
-		y *= imagerecord.pixel_size;	/* Convert pixels to metres^-1 */
+		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;
+		return -1;
 	}
 	psi = atan2(y, x);
 	
@@ -177,14 +182,36 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 	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;
 	
-	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
+	
+	return 0;
+
+}
+
+/* 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, twotheta;
+	ImageReflection refl;
+	
+	refl.parent = imagerecord;
+	refl.x = x;
+	refl.y = y;
+	refl.intensity = intensity;
+	
+	if ( !reflection_map_to_space(&refl, &nx, &ny, &nz, &twotheta) ) {
+		reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
+	}
 
 }
 
@@ -194,3 +221,46 @@ void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) {
 	rctx->last_reflection = r;
 }
 
+double reflection_largest_g(ReflectionContext *reflectionctx) {
+
+	double max = 0.0;
+	Reflection *reflection;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( reflection->type == REFLECTION_NORMAL ) {
+			double mod;
+			mod = modulus(reflection->x, reflection->y, reflection->z);
+			if ( mod > max ) max = mod;
+		}
+		reflection = reflection->next;
+	};
+	
+	return max;
+
+}
+
+Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z) {
+
+	double max = +INFINITY;
+	Reflection *reflection;
+	Reflection *best = NULL;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( reflection->type == REFLECTION_NORMAL ) {
+			double mod;
+			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+			if ( mod < max ) {
+				max = mod;
+				best = reflection;
+			}
+		}
+		reflection = reflection->next;
+	};
+	
+	return best;
+
+
+}
+