Randomise the initial basis

Other preening


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

2 files changed, 45 insertions, 19 deletions

diff --git a/src/ipr.c b/src/ipr.c
index 79368e1..2d78ca1 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -4,6 +4,7 @@
  * Iterative Prediction-Refinement Reconstruction
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * 	    Gordon Ball <gfb21@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -38,6 +39,30 @@ static int ipr_choose_max(Basis *basis) {
 	
 }
 
+static int ipr_random_image(ControlContext *ctx) {
+	double n;
+	n = (double)random()/RAND_MAX;
+	n *= ctx->n_images;
+	return (int)n;
+}
+
+static Reflection *ipr_random_reflection(ReflectionContext *rctx) {
+	
+	double n;
+	int i;
+	Reflection *cur;
+	
+	n = (double)random()/RAND_MAX;
+	n *= rctx->n_reflections/2;
+	cur = rctx->reflections;	/* Jump a maximum of half way in */
+	for ( i=1; i<n; i++ ) {		/* Skip initial placeholder */
+		cur = cur->next;
+	}
+	
+	return cur;
+	
+}
+
 /* Try to turn the basis into a more sensible one */
 static void ipr_try_compact(Basis *basis) {
 
@@ -114,9 +139,10 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 	basis->a.x = 0;	basis->a.y = 0;	basis->a.z = 0;	basis->a.modulus = 1e30;
 	basis->b.x = 0;	basis->b.y = 0;	basis->b.z = 0;	basis->b.modulus = 1e30;
 	basis->c.x = 0;	basis->c.y = 0;	basis->c.z = 0;	basis->c.modulus = 1e30;
-	reflection = ctx->reflectionctx->reflections->next;	/* First reflection, skip initial placeholder */
 	cur_max = 1;
 	
+	reflection = ipr_random_reflection(ctx->reflectionctx);
+	
 	do {
 		double mod;
 		unsigned int angle_works;
@@ -173,13 +199,17 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 		
 	} while ( reflection || (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) );
 	
+	ipr_try_compact(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;
+
 	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);
 	
-	ipr_try_compact(basis);
-	
 	return basis;
 
 }
@@ -190,6 +220,7 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 	ReflectionContext *ordered;
 	double max_res;
 	signed int h, k, l;
+	int max_order_a, max_order_b, max_order_c;
 	
 	/* NB This assumes the diffraction pattern is "vaguely" centered... */
 	if ( ctx->inputfiletype != INPUT_CACHE) {
@@ -207,13 +238,14 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 	}
 	
 	ordered = reflection_init();
-//	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;
-	for ( h=-30; h<=30; h++ ) {
-		for ( k=-30; k<=30; k++ ) {
-			for ( l=-30; l<=30; l++ ) {
-				
+	
+	max_order_a = max_res/basis->a.modulus;
+	max_order_b = max_res/basis->b.modulus;
+	max_order_c = max_res/basis->c.modulus;
+	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++ ) {
+
 				double x, y, z;
 				
 				x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
@@ -236,13 +268,6 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
 
 }
 
-static int ipr_random_image(ControlContext *ctx) {
-	double n;
-	n = (double)random()/RAND_MAX;
-	n *= ctx->n_images;
-	return (int)n;
-}
-
 static gint ipr_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {
 
 	ImageReflection *refl;
@@ -271,13 +296,14 @@ int ipr_refine(ControlContext *ctx) {
 	Basis *basis;
 	int finished;
 	
+	srand(time(NULL));
+	
 	basis = ipr_choose_initial_basis(ctx);
 	if ( !basis ) return -1;
 	
 	ctx->ipr_basis = basis;
 	ctx->ipr_lat = ipr_generate(ctx, basis);
 	
-	srand(time(NULL));
 	finished = 0;
 	do {
 	
diff --git a/src/reproject.c b/src/reproject.c
index 3106998..7c3ec70 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -23,7 +23,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	ImageReflection *refl;
 	Reflection *reflection;
 	int i;
-	double smax = 0.25e9;
+	double smax = 0.5e8;
 	double tilt, omega;
 	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
 	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */