'Satisfactory-ish' basis finding

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

8 files changed, 195 insertions, 109 deletions

diff --git a/src/displaywindow.c b/src/displaywindow.c
index eedad54..6a3adfe 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -484,7 +484,6 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	glEnd();
 	
 	/* x, y, z pointers */
-#if 0
 	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
 	DRAW_SHORT_POINTER
 	
@@ -499,8 +498,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	glRotatef(-90.0, 0.0, 1.0, 0.0);
 	DRAW_SHORT_POINTER
 	glPopMatrix();
-#endif
-
+	
 	/* Tilt axis */
 	glPushMatrix();
 	/* Images rotate clockwise by omega to put tilt axis at +x,
diff --git a/src/ipr.c b/src/ipr.c
index 47b57e2..dd6ce59 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -27,64 +27,6 @@
 #include "ipr.h"
 #include "displaywindow.h"
 
-#if 0
-static int ipr_random_image(ControlContext *ctx) {
-	double n;
-	n = (double)random()/RAND_MAX;
-	n *= ctx->n_images;
-	return (int)n;
-}
-#endif
-
-static Basis *ipr_choose_random_basis(ControlContext *ctx) {
-
-	Basis		*basis;
-	double		tilt_min;
-	double		tilt_max;
-	double		tilt_mid;
-	ImageRecord	*imagerecord;
-	double		x_temp, y_temp, z_temp;
-	double		scale;
-	double		x, y, z;
-	Reflection	*centre;
-	
-	/* Locate the 'plane' in the middle of the "wedge".
-	 *	This whole procedure assumes there is just one tilt axis. */
-	tilt_min = control_min_tilt(ctx);
-	tilt_max = control_max_tilt(ctx);
-	tilt_mid = tilt_min + (tilt_max-tilt_min)/2;
-	imagerecord = control_image_nearest_tilt(ctx, tilt_mid);
-	
-	/* Apply the last two steps of the mapping transform to get the direction from the origin
-	 *	towards the middle of the wedge */
-	x_temp = 0.0;
-	y_temp = cos(deg2rad(imagerecord->tilt));
-	z_temp = -sin(deg2rad(imagerecord->tilt));
-	x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega));
-	y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega));
-	z = z_temp;
-	
-	/* Find the point in the middle of the "wedge" */
-	scale = reflection_largest_g(ctx->reflectionctx)/4;
-	x *= scale;
-	y *= scale;
-	z *= scale;
-	reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
-	
-	centre = reflection_find_nearest(ctx->reflectionctx, x, y, z);
-	if ( !centre ) return NULL;
-	reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_MARKER);
-	
-	basis = malloc(sizeof(Basis));
-	
-	basis->a.x = 1.842e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
-	basis->b.x = 0.0;  basis->b.y = 1.842e9;  basis->b.z = 0.0;
-	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.842e9;
-	
-	return basis;
-
-}
-
 static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
 
 	int n_indexed, n_counted;
@@ -97,7 +39,7 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
 
 		if ( cur->type == REFLECTION_NORMAL ) {
 			
-			/* Can this basis approximately account for this reflection? */
+			/* Can this basis "approximately" account for this reflection? */
 			double det;
 			double a11, a12, a13, a21, a22, a23, a31, a32, a33;
 			double h, k, l;
@@ -135,40 +77,126 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
 
 }
 
-static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
+static Basis *ipr_find_basis(ControlContext *ctx) {
 
-	Basis *basis;
-	Basis *best_basis;
-	double best_efom;
-	int i;
+	Basis		*basis;
+	double		tilt_min;
+	double		tilt_max;
+	double		tilt_mid;
+	ImageRecord	*imagerecord;
+	double		x_temp, y_temp, z_temp;
+	double		scale;
+	double		x, y, z;
+	Reflection	*centre;
+	int		i;
 	
-	/* Track the best basis throughout the whole procedure */
-	best_basis = NULL;  best_efom = 0;
+	/* Locate the 'plane' in the middle of the "wedge".
+	 *	This whole procedure assumes there is just one tilt axis. */
+	tilt_min = control_min_tilt(ctx);
+	tilt_max = control_max_tilt(ctx);
+	tilt_mid = tilt_min + (tilt_max-tilt_min)/2;
+	imagerecord = control_image_nearest_tilt(ctx, tilt_mid);
 	
-	printf("IP: Finding initial basis using eFOM...\n");
-	for ( i=1; i<=30; i++ ) {
-		
-		double efom;
+	/* Apply the last two steps of the mapping transform to get the direction from the origin
+	 *	towards the middle of the wedge */
+	x_temp = 0.0;
+	y_temp = cos(deg2rad(imagerecord->tilt));
+	z_temp = -sin(deg2rad(imagerecord->tilt));
+	x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega));
+	y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega));
+	z = z_temp;
+	
+	/* Find the point in the middle of the "wedge" */
+	scale = reflection_largest_g(ctx->reflectionctx)/4;
+	x *= scale;
+	y *= scale;
+	z *= scale;
+	reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
+	
+	/* Find an "origin" reflection */
+	centre = reflection_find_nearest(ctx->reflectionctx, x, y, z);
+	if ( !centre ) return NULL;
+	reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED);
+	
+	basis = malloc(sizeof(Basis));
+	for ( i=1; i<=3; i++ ) {
+	
+		Reflection *vector;
+		int accept;
+		double vx, vy, vz;
 		
-		basis = ipr_choose_random_basis(ctx);
-		efom = ipr_efom(ctx->reflectionctx, basis);
+		do {
+			
+			Reflection	*check;
+			int		orders, j, lfom;
+			double		tol;
+			
+			accept = 1;
 		
-		if ( (efom > best_efom) || !best_basis ) {
-			if ( best_basis ) free(best_basis);
-			best_efom = efom;
-			best_basis = basis;
-			printf("IP:    %3i: eFOM = %7.3f %%\n", i, efom*100);
+			/* Find a "candidate vector" reflection */
+			vector = reflection_find_nearest_longer(ctx->reflectionctx, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */
+			if ( !vector ) {
+				printf("IP: Couldn't find enough seeds\n");
+				return NULL;
+			}
+			
+			/* Get vector components (not the coordinates the vector was calculated from!) */
+			vx = vector->x - centre->x;
+			vy = vector->y - centre->y;
+			vz = vector->z - centre->z;
+			
+			/* Proximity test */
+			check = reflection_find_nearest_type(ctx->reflectionctx, vx, vy, vz, REFLECTION_MARKER);
+			if ( check ) {
+				if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) {
+					/* Too close to another seed */
+					accept = 0;
+					continue;
+				}
+			}
+			
+			/* Line FoM test */
+			orders = 100e9 / modulus(vx, vy, vz);	/* Go out by +/- this number of iterations.  Rounding error ignored */
+			lfom = 0;
+			tol = modulus(vx, vy, vz)/10.0;
+			for ( j=-orders; j<=orders; j++ ) {
+				check = reflection_find_nearest(ctx->reflectionctx, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j);
+				if ( check && (distance3d(check->x, check->y, check->z, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j) < tol) ) {
+					lfom++;
+				}
+			}
+			if ( lfom < 1 ) {
+				accept = 0;
+				continue;
+			}
+			
+		} while ( !accept );
+		reflection_add(ctx->reflectionctx, vx, vy, vz, 1.0, REFLECTION_MARKER);
 		
+		switch ( i ) {
+			case 1 : {
+				basis->a.x = vx;
+				basis->a.y = vy;
+				basis->a.z = vz;
+			}
+			case 2 : {
+				basis->b.x = vx;
+				basis->b.y = vy;
+				basis->b.z = vz;
+			}
+			case 3 : {
+				basis->c.x = vx;
+				basis->c.y = vy;
+				basis->c.z = vz;
+			}
 		}
 		
 	}
 	
-	printf("IP: Initial basis:\n");
-	printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus);
-	printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus);
-	printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus);
+	printf("IP: eFOM = %7.3f %%\n", ipr_efom(ctx->reflectionctx, basis)*100);
 	
 	return basis;
+
 }
 
 /* Generate list of reflections to analyse */
@@ -191,14 +219,14 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 		}
 		max_res = max_res / 2;
 	} else {
-		max_res = 2e10; //works until I put some more in the reflect.cache header
+		max_res = 1e10; //works until I put some more in the reflect.cache header
 	}
 	
 	ordered = reflection_init();
 	
-	max_order_a = max_res/basis->a.modulus;
-	max_order_b = max_res/basis->b.modulus;
-	max_order_c = max_res/basis->c.modulus;
+	max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z);
+	max_order_b = max_res/modulus(basis->b.x, basis->b.y, basis->b.z);
+	max_order_c = max_res/modulus(basis->c.x, basis->c.y, basis->c.z);
 	for ( h=-max_order_a; h<=max_order_a; h++ ) {
 		for ( k=-max_order_b; k<=max_order_b; k++ ) {
 			for ( l=-max_order_c; l<=max_order_c; l++ ) {
@@ -217,6 +245,7 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 					}
 					if ( (h!=0) || (k!=0) || (l!=0) ) {
 					//	reflection_add(ctx->reflectionctx, x, y, z, 1, REFLECTION_GENERATED);
+						reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
 					}
 				}
 				
@@ -256,10 +285,11 @@ int ipr_refine(ControlContext *ctx) {
 	Basis *basis;
 	int finished;
 	
-	srand(time(NULL));
-	
-	basis = ipr_choose_initial_basis(ctx);
-	if ( !basis ) return -1;
+	basis = ipr_find_basis(ctx);
+	if ( !basis ) {
+		printf("IP: Unable to find basis\n");
+		return -1;
+	}
 	
 	ctx->ipr_basis = basis;
 	ctx->ipr_lat = ipr_generate(ctx, basis);
diff --git a/src/ipr.h b/src/ipr.h
index 3de3724..e511756 100644
--- a/src/ipr.h
+++ b/src/ipr.h
@@ -23,8 +23,6 @@ typedef struct {
 	double y;
 	double z;
 	
-	double modulus;	/* Convenience */
-	
 } Vector;
 
 typedef struct basis_struct {
diff --git a/src/reflections.c b/src/reflections.c
index 63eaa54..a1f6a2c 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -102,6 +102,7 @@ void reflection_free(ReflectionContext *reflectionctx) {
 Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) {
 
 	Reflection *new_reflection;
+	Reflection *nearest;
 	
 	if ( reflectionctx->list_capped ) return NULL;
 	
@@ -110,6 +111,10 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 		fprintf(stderr, "No further reflections will be stored.  Go and fix the peak detection.\n");
 		reflectionctx->list_capped = 1;
 	}
+	
+	nearest = reflection_find_nearest_type(reflectionctx, x, y, z, type);
+	if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL;
+	
 	reflectionctx->n_reflections++;
 	
 	new_reflection = malloc(sizeof(Reflection));
@@ -119,6 +124,7 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 	new_reflection->z = z;
 	new_reflection->intensity = intensity;
 	new_reflection->type = type;
+	new_reflection->found = 0;
 	
 	reflectionctx->last_reflection->next = new_reflection;
 	reflectionctx->last_reflection = new_reflection;
@@ -248,7 +254,56 @@ Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x,
 	
 	reflection = reflectionctx->reflections;
 	while ( reflection ) {
-		if ( reflection->type == REFLECTION_NORMAL ) {
+		if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
+			double mod;
+			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+			if ( mod < max ) {
+				max = mod;
+				best = reflection;
+			}
+		}
+		reflection = reflection->next;
+	};
+	
+	if ( best ) best->found = 1;
+	return best;
+
+
+}
+
+Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance) {
+
+	double max = +INFINITY;
+	Reflection *reflection;
+	Reflection *best = NULL;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
+			double mod;
+			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
+			if ( (mod < max) && (mod >= min_distance) ) {
+				max = mod;
+				best = reflection;
+			}
+		}
+		reflection = reflection->next;
+	};
+	
+	if ( best ) best->found = 1;
+	return best;
+
+}
+
+Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type) {
+
+	double max = +INFINITY;
+	Reflection *reflection;
+	Reflection *best = NULL;
+	
+	reflection = reflectionctx->reflections;
+	while ( reflection ) {
+		if ( reflection->type == type ) {
 			double mod;
 			mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z);
 			if ( mod < max ) {
diff --git a/src/reflections.h b/src/reflections.h
index 824fcba..9dd281a 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -29,16 +29,17 @@ typedef enum {
 
 typedef struct reflection_struct {
 
-	double x;
-	double y;
-	double z;
-	double intensity;
+	double		x;
+	double		y;
+	double		z;
+	double 		intensity;
 
-	signed int h;
-	signed int k;
-	signed int l;
+	signed int	h;
+	signed int	k;
+	signed int	l;
 
-	ReflectionType type;
+	ReflectionType	type;
+	int		found;		/* This reflection has been used in the seed-finding process */
 
 	struct reflection_struct *next;	/* MUST BE LAST in order for caching to work */
 
@@ -68,6 +69,8 @@ extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *
 extern double reflection_largest_g(ReflectionContext *reflectionctx);
 extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta);
 extern Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z);
+extern Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance);
+extern Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type);
 
 #endif	/* REFLECTION_H */
 
diff --git a/src/reproject.c b/src/reproject.c
index 0971b25..ed359fc 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -33,12 +33,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	tilt = deg2rad(image.tilt);
 	omega = deg2rad(image.omega);
 	
-	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
-	
-	i = 0;
-	
-	reflection = rctx->reflections;
-	
 	/* 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;
@@ -58,6 +52,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	ux = uxt*cos(-omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
 	//reflection_add(ctx->reflectionctx, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
 	
+	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
+	reflection = rctx->reflections;
+	i = 0;
 	do {
 	
 		double xl, yl, zl;
diff --git a/src/utils.c b/src/utils.c
index eaaa036..4618b44 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -39,6 +39,10 @@ double modulus(double x, double y, double z) {
 	return sqrt(x*x + y*y + z*z);
 }
 
+double distance3d(double x1, double y1, double z1, double x2, double y2, double z2) {
+	return modulus(x1-x2, y1-y2, z1-z2);
+}
+
 /* Angle between two vectors.  Answer in radians */
 double angle_between(double x1, double y1, double z1, double x2, double y2, double z2) {
 
diff --git a/src/utils.h b/src/utils.h
index 99fc748..1d268e6 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -27,6 +27,7 @@ extern double angle_between(double x1, double y1, double z1, double x2, double y
 extern double angle_between_d(double x1, double y1, double z1, double x2, double y2, double z2);
 extern double lambda(double voltage);
 extern void matrix_renormalise(gsl_matrix *m);
+extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2);
 
 #define rad2deg(a) ((a)*180/M_PI)
 #define deg2rad(a) ((a)*M_PI/180)