From d2a2f1928752e518ac7b8798175fd9a010674dd0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 18 Dec 2015 14:00:56 +0100 Subject: "Retry" mechanism for indexing This increases the indexing rate a bit in situations where there are lots of weak peaks which, although they may be real, don't help indexing. The weakest 10% of peaks get cut out and the indexing re-run. This also allows multiple hits to be indexed, using the "inelegant peak subtraction method", by retrying indexing in the same way after deleting peaks which are accounted for by the lattice just found. --- libcrystfel/src/image.c | 39 +++++++++++++++++ libcrystfel/src/image.h | 1 + libcrystfel/src/index.c | 114 +++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 153 insertions(+), 1 deletion(-) diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 3dfccea3..cd7047dd 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -95,6 +95,45 @@ ImageFeatureList *image_feature_list_new() } +static int comp(const void *a, const void *b) +{ + const struct imagefeature *ap = a; + const struct imagefeature *bp = b; + + return ap->intensity < bp->intensity; +} + + +/* Strongest first. Returned list is guaranteed not to have any holes + * (feature->valid = 0) */ +ImageFeatureList *sort_peaks(ImageFeatureList *flist) +{ + ImageFeatureList *n = image_feature_list_new(); + int nf, i; + + if ( n == NULL ) return NULL; + + n->features = malloc(flist->n_features*sizeof(struct imagefeature)); + if ( n->features == NULL ) { + free(n); + return NULL; + } + + nf = 0; + for ( i=0; in_features; i++ ) { + struct imagefeature *f; + f = image_get_feature(flist, i); + if ( f == NULL ) continue; + n->features[nf++] = flist->features[i]; + } + n->n_features = nf; + + qsort(n->features, nf, sizeof(struct imagefeature), comp); + + return n; +} + + void image_feature_list_free(ImageFeatureList *flist) { if ( !flist ) return; diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 37904308..e9cd76fc 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -224,6 +224,7 @@ extern Reflection *image_reflection_closest(RefList *rlist, extern int image_feature_count(ImageFeatureList *flist); extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); +extern ImageFeatureList *sort_peaks(ImageFeatureList *flist); extern void image_add_crystal(struct image *image, Crystal *cryst); extern void free_all_crystals(struct image *image); diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 99e40053..2e54859d 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -13,6 +13,7 @@ * 2012 Lorenzo Galli * 2013 Cornelius Gati * 2015 Kenneth Beyerlein + * 2014 Takanori Nakane * * This file is part of CrystFEL. * @@ -40,6 +41,7 @@ #include #include #include +#include #include "image.h" #include "utils.h" @@ -265,10 +267,91 @@ static int try_indexer(struct image *image, IndexingMethod indm, } +static int delete_weakest_peaks(ImageFeatureList *peaks) +{ + int i; + int np, ndel, n; + + np = image_feature_count(peaks); + n = 0; + for ( i=0; in-ndel-1; i-- ) { + image_remove_feature(peaks, i); + } + + return 0; +} + + +static int delete_explained_peaks(struct image *image, Crystal *cr) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + const double min_dist = 0.25; + int i, nspots = 0, nindexed = 0; + + /* Round towards nearest */ + fesetround(1); + + /* Cell basis vectors for this image */ + cell_get_cartesian(crystal_get_cell(cr), &ax, &ay, &az, + &bx, &by, &bz, &cx, &cy, &cz); + + /* Loop over peaks, checking proximity to nearest reflection */ + for ( i=0; ifeatures); i++ ) { + + struct imagefeature *f; + struct rvec q; + double h, k, l, hd, kd, ld; + double dsq; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + nspots++; + + /* Reciprocal space position of found peak */ + q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); + + /* Decimal and fractional Miller indices of nearest + * reciprocal lattice point */ + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + /* Check distance */ + dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); + + if ( sqrt(dsq) < min_dist ) { + image_remove_feature(image->features, i); + nindexed++; + } + } + + /* Return TRUE if not enough peaks to continue */ + return (nspots - nindexed) < 5; + +} + + void index_pattern(struct image *image, IndexingMethod *indms, IndexingPrivate **iprivs) { int n = 0; + ImageFeatureList *peaks; + ImageFeatureList *orig; if ( indms == NULL ) return; @@ -276,14 +359,43 @@ void index_pattern(struct image *image, image->crystals = NULL; image->n_crystals = 0; + peaks = sort_peaks(image->features); + orig = image->features; + image->features = peaks; + while ( indms[n] != INDEXING_NONE ) { - if ( try_indexer(image, indms[n], iprivs[n]) ) break; + int done = 0; + int r; + int ntry = 0; + + do { + + r = try_indexer(image, indms[n], iprivs[n]); + ntry++; + + if ( r == 0 ) { + /* Retry with fewer peaks */ + done = delete_weakest_peaks(image->features); + } else { + /* Remove "used" spots and try for + * another lattice */ + Crystal *cr; + cr = image->crystals[image->n_crystals-1]; + done = delete_explained_peaks(image, cr); + } + + if ( ntry > 5 ) done = 1; + + } while ( !done ); + n++; } image->indexed_by = indms[n]; + image->features = orig; + image_feature_list_free(peaks); } -- cgit v1.2.3