Nicer pattern centering

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

4 files changed, 69 insertions, 14 deletions

diff --git a/src/image.c b/src/image.c
index 5290cb4..8eba5f9 100644
--- a/src/image.c
+++ b/src/image.c
@@ -116,12 +116,11 @@ void image_feature_list_free(ImageFeatureList *flist) {
 
 }
 
-ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y, double *d) {
+ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y, double *d, int *idx) {
 
 	int i;
 	double dmin = +HUGE_VAL;
-	int closest;
-	closest = 0;
+	int closest = 0;
 	
 	for ( i=0; i<flist->n_features; i++ ) {
 	
@@ -136,11 +135,45 @@ ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y,
 	
 	}
 	
-	if ( dmin <= 20.0 ) {
+	if ( dmin < +HUGE_VAL ) {
 		*d = dmin;
+		*idx = closest;
 		return &flist->features[closest];
 	}
+	
+	return NULL;
+
+}
 
+ImageFeature *image_feature_second_closest(ImageFeatureList *flist, double x, double y, double *d, int *idx) {
+
+	int i;
+	double dmin = +HUGE_VAL;
+	int closest = 0;
+	double dfirst;
+	int idxfirst;
+	
+	image_feature_closest(flist, x, y, &dfirst, &idxfirst);
+	
+	for ( i=0; i<flist->n_features; i++ ) {
+	
+		double d;
+		
+		d = distance(flist->features[i].x, flist->features[i].y, x, y);
+	
+		if ( (d < dmin) && (i != idxfirst) ) {
+			dmin = d;
+			closest = i;
+		}
+	
+	}
+	
+	if ( dmin < +HUGE_VAL ) {
+		*d = dmin;
+		*idx = closest;
+		return &flist->features[closest];
+	}
+	
 	return NULL;
 
 }
diff --git a/src/image.h b/src/image.h
index ebabee6..ff18580 100644
--- a/src/image.h
+++ b/src/image.h
@@ -77,7 +77,8 @@ extern void image_feature_list_free(ImageFeatureList *flist);
 extern void image_add_feature(ImageFeatureList *flist, double x, double y, ImageRecord *parent, double intensity);
 extern void image_add_feature_index(ImageFeatureList *flist, double x, double y, ImageRecord *parent, double intensity,
 					signed int h, signed int k, signed int l);
-extern ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y, double *d);
+extern ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y, double *d, int *idx);
+extern ImageFeature *image_feature_second_closest(ImageFeatureList *flist, double x, double y, double *d, int *idx);
 
 #endif	/* IMAGE_H */
 
diff --git a/src/prealign.c b/src/prealign.c
index 955ef4d..09785d2 100644
--- a/src/prealign.c
+++ b/src/prealign.c
@@ -234,7 +234,7 @@ void prealign_fine_centering(ImageList *list) {
 			assert(mask_x>=0);
 			assert(mask_y>=0);
 			
-			printf("Image %3i: centre offset by %i,%i\n", i, mask_x-list->images[i].x_centre, mask_y-list->images[i].y_centre);
+			printf("AL: Image %3i: centre offset by %i,%i\n", i, mask_x-list->images[i].x_centre, mask_y-list->images[i].y_centre);
 			
 			list->images[i].x_centre = mask_x;
 			list->images[i].y_centre = mask_y;
@@ -253,13 +253,28 @@ void prealign_feature_centering(ImageList *list) {
 	
 	for ( i=0; i<list->n_images; i++ ) {
 	
-		double d;
-		ImageFeature *feature;
+		double d1, d2;
+		ImageFeature *feature1;
+		ImageFeature *feature2;
+		int idx;
+		
+		feature1 = image_feature_closest(list->images[i].features, list->images[i].x_centre, list->images[i].y_centre, &d1, &idx);
+		feature2 = image_feature_second_closest(list->images[i].features, list->images[i].x_centre, list->images[i].y_centre, &d2, &idx);
+		
+		printf("AL: Image %i, d1=%f, d2=%f\n", i, d1, d2);
+		
+		if ( (fabs(d2-d1) <= 19.0) && feature1 && feature2 ) {
+			list->images[i].x_centre = (feature1->x + feature2->x)/2;
+			list->images[i].y_centre = (feature1->y + feature2->y)/2;
+		} else {
+			if ( feature1 ) {
+				list->images[i].x_centre = feature1->x;
+				list->images[i].y_centre = feature1->y;
+			} else {
+				printf("AL: Couldn't find centre feature for image %i\n", i);
+			}
+		}
 		
-		feature = image_feature_closest(list->images[i].features, list->images[i].x_centre, list->images[i].y_centre, &d);
-		list->images[i].x_centre = feature->x;
-		list->images[i].y_centre = feature->y;
-	
 	}
 
 }
diff --git a/src/reproject.c b/src/reproject.c
index 095b39a..db471fa 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -186,9 +186,15 @@ void reproject_partner_features(ImageFeatureList *flist, ImageRecord *image) {
 	for ( i=0; i<flist->n_features; i++ ) {
 	
 		double d = 0.0;
+		ImageFeature *partner;
+		int idx;
 		
-		flist->features[i].partner = image_feature_closest(image->features, flist->features[i].x, flist->features[i].y, &d);
-		flist->features[i].partner_d = d;
+		partner = image_feature_closest(image->features, flist->features[i].x, flist->features[i].y, &d, &idx);
+		
+		if ( (d <= 20.0) && partner ) {
+			flist->features[i].partner = partner;
+			flist->features[i].partner_d = d;
+		}
 	
 	}