More work on refinement

Store absolute image coordinates in cache (not relative to centre)


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

1 files changed, 75 insertions, 25 deletions

diff --git a/src/refine.c b/src/refine.c
index a5ae3d5..7838e56 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -27,6 +27,7 @@
 #include "reproject.h"
 #include "control.h"
 #include "mapping.h"
+#include "imagedisplay.h"
 
 /* Return the root sum squared deviation distance for all the "reprojectable" features in an image */
 static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_lattice) {
@@ -63,16 +64,36 @@ typedef enum {
 	INDEX_C = 1<<2
 } RefinementIndex;
 
+static const char *refine_decode(RefinementIndex i) {
+
+	switch ( i ) {
+	
+		case INDEX_A			 : return "a--";
+		case INDEX_B			 : return "-b-";
+		case INDEX_C			 : return "--c";
+		case INDEX_A | INDEX_B		 : return "ab-";
+		case INDEX_A | INDEX_C		 : return "a-c";
+		case INDEX_B | INDEX_C		 : return "-bc";
+		case INDEX_A | INDEX_B | INDEX_C : return "abc";
+	
+		default				 : return "---";
+	
+	}
+
+}
+
 /* Use the IPR algorithm to make "cell" fit the given image */
-static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) {
+static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) {
 
 	ImageFeatureList *rflist;
 	ImageFeatureList *flist;
+	ImageFeature *ret;
 	int i;
 	Basis cd; /* Cell delta */
 	int n_a = 0;
 	int n_b = 0;
 	int n_c = 0;
+	int done = FALSE;
 	
 	cd.a.x = 0.0;	cd.a.y = 0.0;	cd.a.z = 0.0;
 	cd.b.x = 0.0;	cd.b.y = 0.0;	cd.b.z = 0.0;
@@ -131,54 +152,71 @@ static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *ce
 		dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det;
 		dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det;
 		printf("dev(hkl) = %f %f %f\n", dh, dk, dl);
-	
+		
 		delta = 0;
-		if ( fabs(dh) < 0.001 ) delta |= INDEX_A;
-		if ( fabs(dk) < 0.001 ) delta |= INDEX_B;
-		if ( fabs(dl) < 0.001 ) delta |= INDEX_C;
-		/* Typically, 'delta' will have all three bits set */
+		if ( fabs(dh)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_A;
+		if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B;
+		if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C;
 		
 		shared = index & delta;
 		
+		printf(" index: %s\n delta: %s\nshared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared));
+		
 		if ( shared == 0 ) {
 			/* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */
 			shared = index;
+			printf("Pure shear.\n");
 		}
 		
 		if ( shared & INDEX_A ) {
-			double w = abs(h) / (abs(h)+abs(k)+abs(l));
-			cd.a.x += w*dx / h;
-			cd.a.y += w*dy / h;
-			cd.a.z += w*dz / h;
+			double w = (double)abs(h) / (abs(h)+abs(k)+abs(l));
+			printf("w(a) = %f\n", w);
+			cd.a.x += w*dx / (double)h;
+			cd.a.y += w*dy / (double)h;
+			cd.a.z += w*dz / (double)h;
 			n_a++;
 		}
 		if ( shared & INDEX_B ) {
-			double w = abs(k) / (abs(h)+abs(k)+abs(l));
-			cd.b.x += w*dx / k;
-			cd.b.y += w*dy / k;
-			cd.b.z += w*dz / k;
+			double w = (double)abs(k) / (abs(h)+abs(k)+abs(l));
+			printf("w(b) = %f\n", w);
+			cd.b.x += w*dx / (double)k;
+			cd.b.y += w*dy / (double)k;
+			cd.b.z += w*dz / (double)k;
 			n_b++;
 		}
 		if ( shared & INDEX_C ) {
-			double w = abs(l) / (abs(h)+abs(k)+abs(l));
-			cd.c.x += w*dx / l;
-			cd.c.y += w*dy / l;
-			cd.c.z += w*dz / l;
+			double w = (double)abs(l) / (abs(h)+abs(k)+abs(l));
+			printf("w(c) = %f\n", w);
+			cd.c.x += w*dx / (double)l;
+			cd.c.y += w*dy / (double)l;
+			cd.c.z += w*dz / (double)l;
 			n_c++;
 		}
-
+		
+		printf("Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx);
+		printf("Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy);
+		printf("Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz);
+		
+		done = TRUE;
+		
+		break;
+		
 	}
 	
-	image_feature_list_free(rflist);
+	if ( !done ) {
+		printf("RF: No partners found.\n");
+		return NULL;
+	}
 	
+	printf("n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c);
 	if ( n_a ) {
-		cd.a.x /= n_a;	cd.a.y /= n_a;	cd.a.z /= n_a;
+		cd.a.x /= (double)n_a;	cd.a.y /= (double)n_a;	cd.a.z /= (double)n_a;
 	}
 	if ( n_b ) {
-		cd.b.x /= n_b;	cd.b.y /= n_b;	cd.b.z /= n_b;
+		cd.b.x /= (double)n_b;	cd.b.y /= (double)n_b;	cd.b.z /= (double)n_b;
 	}
 	if ( n_c ) {
-		cd.c.x /= n_c;	cd.c.y /= n_c;	cd.c.z /= n_c;
+		cd.c.x /= (double)n_c;	cd.c.y /= (double)n_c;	cd.c.z /= (double)n_c;
 	}
 	
 	printf("Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z);
@@ -189,6 +227,11 @@ static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *ce
 	cell->b.x += cd.b.x;	cell->b.y += cd.b.y;	cell->b.z += cd.b.z;
 	cell->c.x += cd.c.x;	cell->c.y += cd.c.y;	cell->c.z += cd.c.z;
 	
+	ret = rflist->features[i].partner;
+	image_feature_list_free(rflist);
+	
+	return ret;
+	
 }
 
 /* Display a graph of root sum squared deviation distance against some other parameter */
@@ -231,8 +274,15 @@ static gint refine_graph(GtkWidget *step_button, ControlContext *ctx) {
 static gint refine_step(GtkWidget *step_button, ControlContext *ctx) {
 
 	if ( (ctx->reproject_id) && (ctx->cell_lattice) ) {
-		refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
-		reproject_lattice_changed(ctx);
+		
+		ImageFeature *fitted;
+		
+		fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
+		if ( fitted ) {
+			reproject_lattice_changed(ctx);
+			imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3);
+		}
+		
 	} else {
 		displaywindow_error("Please first open the reprojection window and select the image to fit", ctx->dw);
 	}