Lots of stuff

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

7 files changed, 168 insertions, 79 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);
 	
diff --git a/src/cache.c b/src/cache.c
index eb164aa..078fa72 100644
--- a/src/cache.c
+++ b/src/cache.c
@@ -38,7 +38,7 @@ ReflectionList *cache_load(const char *filename) {
 
 	cachedreflection_size = sizeof(Reflection) - sizeof(Reflection *);
 
-	reflectionlist = reflection_init();
+	reflectionlist = reflectionlist_new();
 	f = fopen(filename, "rb");
 	if ( !f ) {
 		fprintf(stderr, "Couldn't read reflections from cache\n");
diff --git a/src/ipr.c b/src/ipr.c
index a5845b0..b74ecca 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -51,7 +51,7 @@ static ReflectionList *ipr_generate(ControlContext *ctx, Basis *basis) {
 		max_res = 1e10; //works until I put some more in the reflect.cache header
 	}
 	
-	ordered = reflection_init();
+	ordered = reflectionlist_new();
 	
 	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);
diff --git a/src/mapping.c b/src/mapping.c
index 8570966..c2fed73 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -23,7 +23,7 @@ ReflectionList *mapping_create(ControlContext *ctx) {
 	int i;
 	
 	/* Create reflection context */
-	reflectionlist = reflection_init();
+	reflectionlist = reflectionlist_new();
 	
 	/* Pass all images through itrans
 	 *   (let itrans add the reflections to reflectionlist for now) */
diff --git a/src/qdrp.c b/src/qdrp.c
index 5fe0c50..5c7df98 100644
--- a/src/qdrp.c
+++ b/src/qdrp.c
@@ -107,7 +107,7 @@ static int qdrp_parseline(ControlContext *ctx, const char *line) {
 		}
 				
 		ctx->started = 1;
-		ctx->reflectionlist = reflection_init();
+		ctx->reflectionlist = reflectionlist_new();
 						
 	}
 					
diff --git a/src/reflections.c b/src/reflections.c
index a0ff76d..56871f0 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -22,26 +22,18 @@
 #include "reflections.h"
 #include "utils.h"
 
-static void reflection_addfirst(ReflectionList *reflectionlist) {
-	/* Create first items on lists - saves faffing later.  Corresponds to a central marker.
-		Reflections are only stored if they have non-zero value. */
-	reflectionlist->reflections = malloc(sizeof(Reflection));
-	reflectionlist->reflections->next = NULL;
-	reflectionlist->reflections->x = 0;
-	reflectionlist->reflections->y = 0;
-	reflectionlist->reflections->z = 0;
-	reflectionlist->reflections->type = REFLECTION_CENTRAL;
-	reflectionlist->last_reflection = reflectionlist->reflections;
-	reflectionlist->n_reflections = 1;
+static void reflectionlist_init(ReflectionList *reflectionlist) {
+	reflectionlist->n_reflections = 0;
 	reflectionlist->list_capped = 0;
+	reflectionlist->reflections = NULL;
+	reflectionlist->last_reflection = NULL;
 }
 
-ReflectionList *reflection_init() {
+ReflectionList *reflectionlist_new() {
 
 	ReflectionList *reflectionlist = malloc(sizeof(ReflectionList));
-	reflectionlist->n_reflections = 0;
-	reflection_addfirst(reflectionlist);
-	reflectionlist->list_capped = 0;
+	
+	reflectionlist_init(reflectionlist);
 	
 	return reflectionlist;
 
@@ -87,15 +79,12 @@ void reflection_clear(ReflectionList *reflectionlist) {
 		reflection = next;
 	} while ( reflection );
 	
-	reflectionlist->n_reflections = 0;
-	reflectionlist->list_capped = 0;
-	reflection_addfirst(reflectionlist);
+	reflectionlist_init(reflectionlist);
 
 }
 
 void reflection_free(ReflectionList *reflectionlist) {
 	reflection_clear(reflectionlist);
-	free(reflectionlist->reflections);
 	free(reflectionlist);
 }
 
@@ -115,8 +104,6 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 	nearest = reflection_find_nearest_type(reflectionlist, x, y, z, type);
 	if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL;
 	
-	reflectionlist->n_reflections++;
-	
 	new_reflection = malloc(sizeof(Reflection));
 	new_reflection->next = NULL;
 	new_reflection->x = x;
@@ -126,8 +113,14 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 	new_reflection->type = type;
 	new_reflection->found = 0;
 	
-	reflectionlist->last_reflection->next = new_reflection;
-	reflectionlist->last_reflection = new_reflection;
+	if ( reflectionlist->last_reflection ) {
+		reflectionlist->last_reflection->next = new_reflection;
+		reflectionlist->last_reflection = new_reflection;
+	} else {
+		reflectionlist->reflections = new_reflection;
+		reflectionlist->last_reflection = new_reflection;
+	}
+	reflectionlist->n_reflections++;
 	
 	return new_reflection;
 	
@@ -222,9 +215,15 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 }
 
 void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r) {
+	if ( reflectionlist->last_reflection ) {
+		reflectionlist->last_reflection->next = r;
+		reflectionlist->last_reflection = r;
+	} else {
+		reflectionlist->reflections = r;
+		reflectionlist->last_reflection = r;
+	}
 	r->next = NULL;
-	reflectionlist->last_reflection->next = r;
-	reflectionlist->last_reflection = r;
+	reflectionlist->n_reflections++;
 }
 
 double reflection_largest_g(ReflectionList *reflectionlist) {
@@ -269,7 +268,7 @@ Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, do
 
 }
 
-Reflection *reflection_find_nearest_longer(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) {
+Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) {
 
 	double max = +INFINITY;
 	Reflection *reflection;
@@ -316,3 +315,39 @@ Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double
 
 }
 
+/* This destroys the lfom values in the input list */
+ReflectionList *reflection_sort_lfom(ReflectionList *in) {
+
+	ReflectionList *out;
+	Reflection *best;
+	
+	out = reflectionlist_new();
+	
+	do {
+	
+		Reflection *reflection;
+		int lfom = 0;
+		
+		reflection = in->reflections;
+		best = NULL;
+		while ( reflection ) {
+			if ( reflection->lfom > lfom ) {
+				best = reflection;
+				lfom = reflection->lfom;
+			}
+			reflection = reflection->next;
+		};
+		
+		if ( best ) {
+			Reflection *new;
+			new = reflection_add(out, best->x, best->y, best->z, best->intensity, best->type);
+			new->lfom = best->lfom;
+			best->lfom = 0;
+		}
+		
+	} while ( best );
+	
+	return out;
+
+}
+
diff --git a/src/reflections.h b/src/reflections.h
index 9ac0cec..c753430 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -39,7 +39,10 @@ typedef struct reflection_struct {
 	signed int	l;
 
 	ReflectionType	type;
+	
+	/* Stuff used when finding bases */
 	int		found;		/* This reflection has been used in the seed-finding process */
+	int		lfom;		/* Line FoM for this seed */
 
 	struct reflection_struct *next;	/* MUST BE LAST in order for caching to work */
 
@@ -55,7 +58,7 @@ typedef struct reflectionlist_struct {
 	
 } ReflectionList;
 
-extern ReflectionList *reflection_init(void);
+extern ReflectionList *reflectionlist_new(void);
 extern void reflection_clear(ReflectionList *reflectionlist);
 extern void reflection_clear_markers(ReflectionList *reflectionlist);
 extern void reflection_free(ReflectionList *reflectionlist);
@@ -69,8 +72,9 @@ extern void reflection_add_from_reflection(ReflectionList *reflectionlist, Refle
 extern double reflection_largest_g(ReflectionList *reflectionlist);
 extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta);
 extern Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, double y, double z);
-extern Reflection *reflection_find_nearest_longer(ReflectionList *reflectionlist, double x, double y, double z, double min_distance);
+extern Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance);
 extern Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type);
+extern ReflectionList *reflection_sort_lfom(ReflectionList *in);
 
 #endif	/* REFLECTION_H */