Lots of stuff

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

1 files changed, 97 insertions, 47 deletions

diff --git a/src/basis.c b/src/basis.c
index 939dc93..19004ee 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -102,6 +102,7 @@ static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) {
 
 }
 
+/* Select a suitable number of sensible seed vectors */
 static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 
 	double		tilt_min;
@@ -115,7 +116,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 	int		i;
 	ReflectionList *seeds;
 	
-	seeds = reflection_init();
+	seeds = reflectionlist_new();
 	
 	/* Locate the 'plane' in the middle of the "wedge".
 	 *	This whole procedure assumes there is just one tilt axis. */
@@ -134,7 +135,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 	z = z_temp;
 	
 	/* Find the point in the middle of the "wedge" */
-	scale = reflection_largest_g(ctx->reflectionlist)/4;
+	scale = reflection_largest_g(ctx->reflectionlist)/6;
 	x *= scale;
 	y *= scale;
 	z *= scale;
@@ -149,21 +150,21 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 	for ( i=1; i<=10; i++ ) {
 	
 		Reflection *vector;
-		int accept;
+		Reflection *new_seed;
+		int accept, lfom;
 		double vx, vy, vz;
 		
 		do {
 			
 			Reflection	*check;
-			int		lfom;
 			
 			accept = 1;
 		
 			/* Find a "candidate vector" reflection */
-			vector = reflection_find_nearest_longer(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */
+			vector = reflection_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */
 			if ( !vector ) {
-				printf("BS: Couldn't find enough seeds\n");
-				return NULL;
+				printf("BS: Only found %i seeds\n", i);
+				return seeds;
 			}
 			vector->found = 1;
 			
@@ -172,8 +173,8 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 			vy = vector->y - centre->y;
 			vz = vector->z - centre->z;
 			
-			/* Proximity test */
-			check = reflection_find_nearest_type(ctx->reflectionlist, vx, vy, vz, REFLECTION_NORMAL);
+			/* Proximity test: don't duplicate seeds */
+			check = reflection_find_nearest(seeds, vx, vy, vz);
 			if ( check ) {
 				if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) {
 					/* Too close to another seed */
@@ -182,17 +183,16 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 				}
 			}
 			
-			/* lFOM test */
+			/* Record lFOM for later analysis */
 			lfom = basis_lfom(ctx, vx, vy, vz);
-			if ( lfom < 1 ) {
-				accept = 0;
-				continue;
-			}
-			printf("lfom=%i\n", lfom);
 			
 		} while ( !accept );
 
-		reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL);
+		/* Add to the list of seeds */
+		new_seed = reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL);
+		new_seed->lfom = lfom;
+		
+		/* Create a marker in the default list for visualisation */
 		reflection_add(ctx->reflectionlist, vx, vy, vz, 1.0, REFLECTION_MARKER);
 
 	}
@@ -201,48 +201,98 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) {
 
 }
 
-Basis *basis_find(ControlContext *ctx) {
+/* Assemble the most sensible basis from seeds */
+static Basis *basis_assemble_seeds(ReflectionList *seeds) {
 
 	Basis *basis;
-	ReflectionList *seeds;
 	Reflection *ref;
-	int i;
-	
-	/* Get the shortlist of seeds */
-	seeds = basis_find_seeds(ctx);
+	int i, j, b;
+	ReflectionList *seeds_sorted;
 
-	/* Assemble the seeds into a basis */
-	basis = malloc(sizeof(Basis));
-	ref = seeds->reflections;
-	for ( i=1; i<=3; i++ ) {
+	seeds_sorted = reflection_sort_lfom(seeds);
 	
-		double vx, vy, vz;
+	basis = malloc(10*sizeof(Basis));
 	
-		vx = ref->x;
-		vy = ref->y;
-		vz = ref->z;
+	for ( b=0; b<10; b++ ) {
 	
-		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;
+		ref = seeds_sorted->reflections;
+		j = 0;	/* Number of basis components already found */
+		for ( i=1; i<=seeds_sorted->n_reflections; i++ ) {
+		
+			double vx, vy, vz;
+		
+			vx = ref->x;
+			vy = ref->y;
+			vz = ref->z;
+			
+			printf("Seed %2i: lFOM=%6i", i, ref->lfom);
+			
+			switch ( j ) {
+				case 0 : {
+					basis[b].a.x = vx;
+					basis[b].a.y = vy;
+					basis[b].a.z = vz;
+					j++;
+					printf(" *");
+					break;
+				}
+				case 1 : {
+					if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6)
+					  && (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) < M_PI-M_PI/6) ) {
+						basis[b].b.x = vx;
+						basis[b].b.y = vy;
+						basis[b].b.z = vz;
+						j++;
+						printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z)));
+						break;
+					}
+				}
+				case 2 : {
+					double cx, cy, cz;
+					cx = basis[b].a.y*basis[b].b.z - basis[b].a.z*basis[b].b.y;
+					cy =  - basis[b].a.x*basis[b].b.z + basis[b].a.z*basis[b].b.x;
+					cz = basis[b].a.x*basis[b].b.y - basis[b].a.y*basis[b].b.x;
+					if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6)
+					  && (angle_between(vx, vy, vz, basis[b].b.x, basis[b].b.y, basis[b].b.z) > M_PI/6)
+					  && (angle_between(vx, vy, vz, cx, cy, cz) < M_PI/2-M_PI/6) ) {
+						basis[b].c.x = vx;
+						basis[b].c.y = vy;
+						basis[b].c.z = vz;
+						j++;
+						printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, cx, cy, cz)));
+						break;
+					}
+				}
 			}
+			
+			printf("\n");
+			if ( j >= 3 ) break;
+			ref = ref->next;
+			
 		}
-		
-		ref = ref->next;
+		if ( j < 3 ) {
+			break;
+		}
+	}
 	
+	return basis;
+
+}
+
+Basis *basis_find(ControlContext *ctx) {
+
+	Basis *basis;
+	ReflectionList *seeds;
+	
+	/* Get the shortlist of seeds */
+	seeds = basis_find_seeds(ctx);
+	if ( seeds->n_reflections < 3 ) {
+		printf("BS: Not enough seeds to form a basis\n");
+		return NULL;
 	}
+	printf("BS: Found %i seeds\n", seeds->n_reflections);
+
+	basis = basis_assemble_seeds(seeds);
 	
 	printf("BS: eFOM = %7.3f %%\n", basis_efom(ctx->reflectionlist, basis)*100);