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

10 files changed, 104 insertions, 48 deletions

diff --git a/src/basis.c b/src/basis.c
index 841ccc6..08e3edb 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -18,7 +18,7 @@
 #include "reflections.h"
 #include "basis.h"
 
-static double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
+double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
 
 	int n_indexed, n_counted;
 	Reflection *cur;
diff --git a/src/basis.h b/src/basis.h
index 9fb1cda..5803cab 100644
--- a/src/basis.h
+++ b/src/basis.h
@@ -32,5 +32,7 @@ typedef struct basis_struct {
 	Vector c;
 } Basis;
 
+extern double basis_efom(struct reflectionlist_struct *reflectionlist, Basis *basis);
+
 #endif	/* BASIS_H */
 
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 3877da8..ceb36ff 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -210,11 +210,14 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 		switch ( cur->type ) {
 			case IMAGEDISPLAY_MARK_CIRCLE_1 : gc = imagedisplay->gc_marks_1; break;
 			case IMAGEDISPLAY_MARK_CIRCLE_2 : gc = imagedisplay->gc_marks_2; break;
+			case IMAGEDISPLAY_MARK_CIRCLE_3 : gc = imagedisplay->gc_marks_3; break;
 			case IMAGEDISPLAY_MARK_LINE_1 : gc = imagedisplay->gc_marks_1; break;
+			case IMAGEDISPLAY_MARK_LINE_2 : gc = imagedisplay->gc_marks_2; break;
 			default : gc = imagedisplay->gc_marks_1; break;
 		}
 		
-		if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2) ) {
+		if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2)
+			 || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_3) ) {
 		
 			gdk_draw_arc(drawingarea->window, gc, FALSE,
 						xoffs + cur->x*scale - 5,
@@ -258,6 +261,10 @@ static gint imagedisplay_realize(GtkWidget *widget, ImageDisplay *imagedisplay)
 	gdk_color_parse("#00dd00", &colour);
 	gdk_gc_set_rgb_fg_color(imagedisplay->gc_marks_2, &colour);
 
+	imagedisplay->gc_marks_3 = gdk_gc_new(imagedisplay->drawingarea->window);
+	gdk_color_parse("#00ddff", &colour);
+	gdk_gc_set_rgb_fg_color(imagedisplay->gc_marks_3, &colour);
+
 	imagedisplay->realised = TRUE;
 
 	return 0;
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index 3b7192a..d404bcb 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -31,6 +31,7 @@ typedef enum {
 typedef enum {
 	IMAGEDISPLAY_MARK_CIRCLE_1,
 	IMAGEDISPLAY_MARK_CIRCLE_2,
+	IMAGEDISPLAY_MARK_CIRCLE_3,
 	IMAGEDISPLAY_MARK_LINE_1,
 	IMAGEDISPLAY_MARK_LINE_2
 } ImageDisplayMarkType;
@@ -66,6 +67,7 @@ typedef struct imagedisplay_struct {
 	GdkGC			*gc_tiltaxis;
 	GdkGC			*gc_marks_1;
 	GdkGC			*gc_marks_2;
+	GdkGC			*gc_marks_3;
 	gboolean		realised;
 	
 	unsigned int		drawingarea_width;
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index e08bd26..67ed96b 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -519,7 +519,7 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) {
 	
 		px = gsl_matrix_get(p,0,i);
 		py = gsl_matrix_get(p,1,i);
-		image_add_feature(flist, (px-imagerecord->x_centre), (py-imagerecord->y_centre), imagerecord, 1.0);
+		image_add_feature(flist, px, py, imagerecord, 1.0);
 		
 	}
 	
diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c
index 17a2cbe..87a461f 100644
--- a/src/itrans-threshold.c
+++ b/src/itrans-threshold.c
@@ -45,7 +45,7 @@ ImageFeatureList *itrans_peaksearch_threshold(ImageRecord *imagerecord, ControlC
 				assert(y<height);
 				assert(x>=0);
 				assert(y>=0);
-				image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]);
+				image_add_feature(flist, x, y, imagerecord, image[x + width*y]);
 			}
 		}
 	}
@@ -100,7 +100,7 @@ ImageFeatureList *itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord,
 			assert(max_y<height);
 			assert(max_x>=0);
 			assert(max_y>=0);
-			image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]);
+			image_add_feature(flist, x, y, imagerecord, image[x + width*y]);
 			
 			/* Remove it and its surroundings */
 			for ( y=((max_y-10>0)?max_y-10:0); y<((max_y+10)<height?max_y+10:height); y++ ) {
diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c
index d984774..79e2468 100644
--- a/src/itrans-zaefferer.c
+++ b/src/itrans-zaefferer.c
@@ -84,8 +84,7 @@ ImageFeatureList *itrans_peaksearch_zaefferer(ImageRecord *imagerecord) {
 					assert(mask_y<height);
 					assert(mask_x>=0);
 					assert(mask_y>=0);
-					image_add_feature(flist, (mask_x-imagerecord->x_centre), (mask_y-imagerecord->y_centre),
-									imagerecord, image[mask_x + width*mask_y]);
+					image_add_feature(flist, mask_x, mask_y, imagerecord, image[mask_x + width*mask_y]);
 				}
 				
 			}
diff --git a/src/mapping.c b/src/mapping.c
index 668bfc4..afe36c4 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -38,7 +38,7 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 	double nx, ny, nz;
 	
 	imagerecord = refl->parent;
-	x = refl->x;	y = refl->y;
+	x = refl->x - imagerecord->x_centre;	y = refl->y - imagerecord->y_centre;
 	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
 	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
 	k = 1/imagerecord->lambda;
@@ -95,7 +95,7 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 static void mapping_map_features(ControlContext *ctx) {
 
 	int i;
-
+	
 	/* Create reflection list for measured reflections */
 	if ( ctx->reflectionlist ) reflectionlist_free(ctx->reflectionlist);
 	ctx->reflectionlist = reflectionlist_new();
@@ -114,7 +114,7 @@ static void mapping_map_features(ControlContext *ctx) {
 				reflection_add(ctx->reflectionlist, nx, ny, nz,
 					ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
 			}
-
+			
 		}
 		
 	}
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);
 	}
diff --git a/src/reproject.c b/src/reproject.c
index c0dfc73..2ff7e6c 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -85,7 +85,6 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 			double theta;
 			double x, y;
 			double rx, ry, rz;
-			double xtest, ytest;
 			
 			/* Determine the intersection point */
 			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
@@ -154,10 +153,11 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 					return NULL;
 				}
 				
+				x += image->x_centre;
+				y += image->y_centre;
+				
 				/* Sanity check */
-				xtest = x + image->x_centre;
-				ytest = y + image->y_centre;
-				if ( (xtest>=0) && (xtest<image->width) && (ytest>=0) && (ytest<image->height) ) {
+				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
 				
 					/* Record the reflection */
 					image_add_feature_index(flist, x, y, image, reflection->intensity,
@@ -171,7 +171,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 		}
 		
 		reflection = reflection->next;
-
+	
 	} while ( reflection );
 	
 	return flist;
@@ -229,29 +229,25 @@ static void reproject_mark_peaks(ControlContext *ctx) {
 	ImageFeatureList *rflist;
 	ImageFeatureList *flist;
 	size_t j;
-	double xc, yc;
-	
-	xc = ctx->images->images[ctx->reproject_cur_image].x_centre;
-	yc = ctx->images->images[ctx->reproject_cur_image].y_centre;
 	
 	/* Draw the reprojected peaks */
 	rflist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice);
 	for ( j=0; j<rflist->n_features; j++ ) {
-		imagedisplay_add_mark(ctx->reproject_id, xc+rflist->features[j].x, yc+rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
+		imagedisplay_add_mark(ctx->reproject_id, rflist->features[j].x, rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
 	}
 	
 	/* Now draw the original measured peaks */
 	flist = ctx->images->images[ctx->reproject_cur_image].features;
 	for ( j=0; j<flist->n_features; j++ ) {
-		imagedisplay_add_mark(ctx->reproject_id, xc+flist->features[j].x, yc+flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2);
+		imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2);
 	}
 	
 	/* Now connect partners */
 	reproject_partner_features(rflist, &ctx->images->images[ctx->reproject_cur_image]);
 	for ( j=0; j<rflist->n_features; j++ ) {
 		if ( rflist->features[j].partner ) {
-			imagedisplay_add_line(ctx->reproject_id, xc+rflist->features[j].x, yc+rflist->features[j].y,
-							 xc+rflist->features[j].partner->x, yc+rflist->features[j].partner->y,
+			imagedisplay_add_line(ctx->reproject_id, rflist->features[j].x, rflist->features[j].y,
+							 rflist->features[j].partner->x, rflist->features[j].partner->y,
 							 IMAGEDISPLAY_MARK_LINE_1);
 		}
 	}