Don't recalculate reprojections unnecessarily

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

4 files changed, 55 insertions, 46 deletions

diff --git a/src/image.c b/src/image.c
index 8eba5f9..3bb6ba1 100644
--- a/src/image.c
+++ b/src/image.c
@@ -37,6 +37,7 @@ int image_add(ImageList *list, uint16_t *image, int width, int height, double ti
 	list->images[list->n_images].y_centre = ctx->y_centre;
 	list->images[list->n_images].slop = 0.0;
 	list->images[list->n_images].features = NULL;
+	list->images[list->n_images].rflist = NULL;
 	
 	if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
 		list->images[list->n_images].pixel_size = ctx->pixel_size;
diff --git a/src/image.h b/src/image.h
index ff18580..178f363 100644
--- a/src/image.h
+++ b/src/image.h
@@ -59,7 +59,8 @@ typedef struct imagerecord_struct {
 	int			x_centre;
 	int			y_centre;
 	
-	ImageFeatureList	*features;
+	ImageFeatureList	*features;	/* "Experimental" features */
+	ImageFeatureList	*rflist;	/* "Predicted" features */
 
 } ImageRecord;
 
diff --git a/src/refine.c b/src/refine.c
index da578b4..9b901d0 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -88,10 +88,8 @@ static const char *refine_decode(RefinementIndex i) {
 
 /* Use the IPR algorithm to make "cell" fit the given image */
 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;
@@ -103,11 +101,12 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 	cd.b.x = 0.0;	cd.b.y = 0.0;	cd.b.z = 0.0;
 	cd.c.x = 0.0;	cd.c.y = 0.0;	cd.c.z = 0.0;
 	
-	rflist = reproject_get_reflections(image, cell_lattice);
+	if ( !image->rflist ) {
+		image->rflist = reproject_get_reflections(image, cell_lattice);
+	}
 	flist = image->features;
-	reproject_partner_features(rflist, image);
 	
-	for ( i=0; i<rflist->n_features; i++ ) {
+	for ( i=0; i<image->rflist->n_features; i++ ) {
 		
 		double dix, diy;
 		double dx, dy, dz, twotheta;
@@ -118,25 +117,25 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 		double dh, dk, dl;
 		
 		/* Skip if no partner */
-		if ( !rflist->features[i].partner ) continue;
+		if ( !image->rflist->features[i].partner ) continue;
 		
-		h = rflist->features[i].h;
-		k = rflist->features[i].k;
-		l = rflist->features[i].l;
+		h = image->rflist->features[i].h;
+		k = image->rflist->features[i].k;
+		l = image->rflist->features[i].l;
 		
 		/* Determine the difference vector */
-		dix = rflist->features[i].partner->x - rflist->features[i].x;
-		diy = rflist->features[i].partner->y - rflist->features[i].y;
+		dix = image->rflist->features[i].partner->x - image->rflist->features[i].x;
+		diy = image->rflist->features[i].partner->y - image->rflist->features[i].y;
 	//	printf("RF: Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy);
 		
 		/* Map the difference vector to the relevant tilted plane */
-		old_x = rflist->features[i].partner->x;
-		old_y = rflist->features[i].partner->y;
-		rflist->features[i].partner->x = dix + rflist->features[i].partner->parent->x_centre;
-		rflist->features[i].partner->y = diy + rflist->features[i].partner->parent->y_centre;
-		mapping_map_to_space(rflist->features[i].partner, &dx, &dy, &dz, &twotheta);
-		rflist->features[i].partner->x = old_x;
-		rflist->features[i].partner->y = old_y;
+		old_x = image->rflist->features[i].partner->x;
+		old_y = image->rflist->features[i].partner->y;
+		image->rflist->features[i].partner->x = dix + image->rflist->features[i].partner->parent->x_centre;
+		image->rflist->features[i].partner->y = diy + image->rflist->features[i].partner->parent->y_centre;
+		mapping_map_to_space(image->rflist->features[i].partner, &dx, &dy, &dz, &twotheta);
+		image->rflist->features[i].partner->x = old_x;
+		image->rflist->features[i].partner->y = old_y;
 	//	printf("RF: dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3);
 		
 		/* Select the basis vectors which are allowed to be altered */
@@ -230,8 +229,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 	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 image->rflist->features[i].partner;
 	
 	return NULL;
 	
@@ -307,32 +305,33 @@ static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp)
 
 		/* Ensure lattice is up to date */
 		reproject_lattice_changed(ctx);
+		ctx->images->images[i].rflist = NULL;	/* Invalidate reprojection for this image - it's wrong */
 	
 		if ( is_odd(i) ) {
 			
 			/* Odd-numbered images: measure */
-			ImageFeatureList *rflist;
 			ImageFeatureList *flist;
 			ImageRecord *image;
 			int j, n;
 			double image_dev_mean = 0;
 			
 			image = &ctx->images->images[i];
-			rflist = reproject_get_reflections(image, ctx->cell_lattice);
+			if ( !image->rflist ) {
+				image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+			}
 			flist = image->features;
-			reproject_partner_features(rflist, image);
 			
 			n = 0;
-			for ( j=0; j<rflist->n_features; j++ ) {
+			for ( j=0; j<image->rflist->n_features; j++ ) {
 				
 				double dix, diy;
 				
 				/* Skip if no partner */
-				if ( !rflist->features[j].partner ) continue;
+				if ( !image->rflist->features[j].partner ) continue;
 				
 				/* Determine the difference vector */
-				dix = rflist->features[j].partner->x - rflist->features[j].x;
-				diy = rflist->features[j].partner->y - rflist->features[j].y;
+				dix = image->rflist->features[j].partner->x - image->rflist->features[j].x;
+				diy = image->rflist->features[j].partner->y - image->rflist->features[j].y;
 				
 				image_dev_mean += sqrt(dix*dix + diy*diy);
 				n++;
@@ -345,14 +344,11 @@ static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp)
 			series_dev_mean += image_dev_mean;
 			series_dev_n++;
 			
-			image_feature_list_free(rflist);
-			
 		} else {
 			
 			/* Even-numbered images: fit */
 			refine_fit_image(ctx->cell, &ctx->images->images[i], ctx->cell_lattice);
-			//displaywindow_update(ctx->dw);
-
+			
 		}
 
 	}
@@ -375,6 +371,7 @@ static gint refine_sequence(GtkWidget *step_button, ControlContext *ctx) {
 	GtkWidget *graph_warp;
 	double *warp_vals;
 	size_t idx;
+	ImageDisplay *id;
 	
 	fit_vals = malloc(401*sizeof(double));
 	warp_vals = malloc(401*sizeof(double));
@@ -385,6 +382,10 @@ static gint refine_sequence(GtkWidget *step_button, ControlContext *ctx) {
 		return 0;
 	}
 	
+	/* Temporarily disable ImageDisplay stuff */
+	id = ctx->reproject_id;
+	ctx->reproject_id = NULL;
+	
 	for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) {
 	
 		double fit, warp;
@@ -411,6 +412,9 @@ static gint refine_sequence(GtkWidget *step_button, ControlContext *ctx) {
 		
 	}
 	
+	ctx->reproject_id = id;
+	reproject_lattice_changed(ctx);
+	
 	window_fit = gtk_window_new(GTK_WINDOW_TOPLEVEL);
 	gtk_window_set_default_size(GTK_WINDOW(window_fit), 640, 256);
 	gtk_window_set_title(GTK_WINDOW(window_fit), "Omega-Search Graph: Fit");
diff --git a/src/reproject.c b/src/reproject.c
index db471fa..9efe965 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -18,6 +18,7 @@
 #include "imagedisplay.h"
 #include "displaywindow.h"
 #include "image.h"
+#include "reproject.h"
 
 ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
 
@@ -174,6 +175,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	
 	} while ( reflection );
 	
+	reproject_partner_features(flist, image);
+	
 	return flist;
 
 }
@@ -202,40 +205,40 @@ void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image) {
 
 static void reproject_mark_peaks(ControlContext *ctx) {
 
-	ImageFeatureList *rflist;
 	ImageFeatureList *flist;
 	size_t j;
+	ImageRecord *image;
+	
+	image = &ctx->images->images[ctx->reproject_cur_image];
 	
 	/* 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, rflist->features[j].x, rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
+	if ( !image->rflist ) {
+		image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+	}
+	for ( j=0; j<image->rflist->n_features; j++ ) {
+		imagedisplay_add_mark(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
 	}
 	
 	/* Now draw the original measured peaks */
-	flist = ctx->images->images[ctx->reproject_cur_image].features;
+	flist = image->features;
 	for ( j=0; j<flist->n_features; j++ ) {
 		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, rflist->features[j].x, rflist->features[j].y,
-							 rflist->features[j].partner->x, rflist->features[j].partner->y,
+	for ( j=0; j<image->rflist->n_features; j++ ) {
+		if ( image->rflist->features[j].partner ) {
+			imagedisplay_add_line(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y,
+							 image->rflist->features[j].partner->x, image->rflist->features[j].partner->y,
 							 IMAGEDISPLAY_MARK_LINE_1);
 		}
 	}	
 	
-	image_feature_list_free(rflist);
-	
 }
 	
 static void reproject_update(ControlContext *ctx) {
 
 	imagedisplay_clear_marks(ctx->reproject_id);
-	//reflectionlist_clear_markers(ctx->reflectionlist);
 	
 	reproject_mark_peaks(ctx);