Feature partnering algorithm (of sorts)

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

1 files changed, 55 insertions, 5 deletions

diff --git a/src/reproject.c b/src/reproject.c
index b13dd66..535dc0e 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -177,8 +177,49 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 }
 
+/* Attempt to find the partner for "feature" from the feature list of "image" */
+static ImageFeature *reproject_find_partner(ImageFeature *feature, ImageRecord *image) {
+
+	int i;
+	double dmin = +HUGE_VAL;
+	int closest;
+	closest = 0;
+	
+	for ( i=0; i<image->features->n_features; i++ ) {
+	
+		double d;
+		
+		d = distance(image->features->features[i].x, image->features->features[i].x, feature->x, feature->y);
+	
+		if ( d < dmin ) {
+			dmin = d;
+			closest = i;
+		}
+	
+	}
+	
+	if ( dmin < 20 ) {
+		return &image->features->features[closest];
+	}
+	
+	return NULL; /* No suitable partner found */
+
+}
+
+/* Attempt to find partners from the feature list of "image" for each feature in "flist". */
+void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image) {
+
+	int i;
+	
+	for ( i=0; i<flist->n_features; i++ ) {
+		flist->features[i].partner = reproject_find_partner(&flist->features[i], image);
+	}
+
+}
+
 static void reproject_mark_peaks(ControlContext *ctx) {
 
+	ImageFeatureList *rflist;
 	ImageFeatureList *flist;
 	size_t j;
 	double xc, yc;
@@ -187,11 +228,10 @@ static void reproject_mark_peaks(ControlContext *ctx) {
 	yc = ctx->images->images[ctx->reproject_cur_image].y_centre;
 	
 	/* Draw the reprojected peaks */
-	flist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
-	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_1);
+	rflist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	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);
 	}
-	image_feature_list_free(flist);
 	
 	/* Now draw the original measured peaks */
 	flist = ctx->images->images[ctx->reproject_cur_image].features;
@@ -199,6 +239,16 @@ static void reproject_mark_peaks(ControlContext *ctx) {
 		imagedisplay_add_mark(ctx->reproject_id, xc+flist->features[j].x, yc+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++ ) {
+		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_MARK_LINE_1);
+	}	
+	
+	image_feature_list_free(rflist);
+	
 }
 
 static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {
@@ -207,7 +257,7 @@ static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlC
 	ctx->reproject_cur_image++;
 	if ( ctx->reproject_cur_image == ctx->images->n_images ) ctx->reproject_cur_image = 0;
 	
-	imagedisplay_clear_circles(ctx->reproject_id);
+	imagedisplay_clear_marks(ctx->reproject_id);
 	reflectionlist_clear_markers(ctx->reflectionlist);
 	
 	reproject_mark_peaks(ctx);