diff options
-rw-r--r-- | src/image.c | 41 | ||||
-rw-r--r-- | src/image.h | 3 | ||||
-rw-r--r-- | src/prealign.c | 29 | ||||
-rw-r--r-- | src/reproject.c | 10 |
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; + } } |