diff options
| -rw-r--r-- | libcrystfel/src/image.c | 39 | ||||
| -rw-r--r-- | libcrystfel/src/image.h | 1 | ||||
| -rw-r--r-- | libcrystfel/src/index.c | 114 |
3 files changed, 153 insertions, 1 deletions
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; i<flist->n_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 <cornelius.gati@cfel.de> * 2015 Kenneth Beyerlein <kenneth.beyerlein@desy.de> + * 2014 Takanori Nakane <nakane.t@gmail.com> * * This file is part of CrystFEL. * @@ -40,6 +41,7 @@ #include <math.h> #include <string.h> #include <assert.h> +#include <fenv.h> #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; i<np; i++ ) { + if ( image_get_feature(peaks, i) != NULL ) n++; + } + + if ( n < 6 ) return 1; + + /* Delete 10% weakest peaks */ + ndel = n/10; + if ( ndel == 0 ) return 1; + + for ( i=n-1; i>n-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; i<image_feature_count(image->features); 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); } |