diff options
| -rw-r--r-- | src/index.c | 12 | ||||
| -rw-r--r-- | src/peaks.c | 104 | ||||
| -rw-r--r-- | src/peaks.h | 2 |
3 files changed, 89 insertions, 29 deletions
diff --git a/src/index.c b/src/index.c index ced8bbb8..d5e76c50 100644 --- a/src/index.c +++ b/src/index.c @@ -205,17 +205,15 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, if ( new_cell == NULL ) continue; /* Sanity check */ - image->reflections = find_intersections(image, - new_cell); - if ( !config_insane && - !peak_sanity_check(image->reflections, - image->features) ) - { + image->reflections = find_intersections(image, new_cell); + image->indexed_cell = new_cell; + + if ( !config_insane && !peak_sanity_check(image) ) { cell_free(new_cell); + image->indexed_cell = NULL; continue; } - image->indexed_cell = new_cell; goto done; /* Success */ } diff --git a/src/peaks.c b/src/peaks.c index 92556ea6..a5e1983a 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -462,46 +462,108 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } -int peak_sanity_check(RefList *rlist, ImageFeatureList *flist) +int peak_sanity_check(struct image *image) { + int i; int n_feat = 0; int n_sane = 0; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double min_dist = 0.25; + + /* Round towards nearest */ + fesetround(1); + + /* Cell basis vectors for this image */ + cell_get_cartesian(image->indexed_cell, &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++ ) { - for ( i=0; i<image_feature_count(flist); i++ ) { - - double dist; struct imagefeature *f; - Reflection *refl; - RefListIterator *iter; + struct rvec q; + double h,k,l,hd,kd,ld; - f = image_get_feature(flist, i); + /* Assume all image "features" are genuine peaks */ + f = image_get_feature(image->features, i); if ( f == NULL ) continue; n_feat++; - /* Find closest predicted peak */ - - for ( refl = first_refl(rlist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double fs, ss; - get_detector_pos(refl, &fs, &ss); - dist = sqrt(pow(fs-f->fs, 2.0) + pow(ss-f->ss, 2.0)); - if ( dist < PEAK_CLOSE ) { - n_sane++; - continue; - } + /* 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 Bragg reflection */ + 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 */ + if ( fabs(h - hd) < min_dist && + fabs(k - kd) < min_dist && + fabs(l - ld) < min_dist ) { + // printf("%6.1f %6.1f %3g %3g %3g %6.2f %6.2f %6.2f\n",f->fs,f->ss,h,k,l,hd,kd,ld); + n_sane++; + continue; } } - if ( (float)n_sane / (float)n_feat < 0.1 ) return 0; + /* return 0 means fail test, return 1 means pass test */ + // printf("%d out of %d peaks are \"sane\"\n",n_sane,n_feat); + if ( (float)n_sane / (float)n_feat < 0.5 ) return 0; return 1; } + +//int peak_sanity_check(RefList *rlist, ImageFeatureList *flist) +//{ +// int i; +// int n_feat = 0; +// int n_sane = 0; +// +// for ( i=0; i<image_feature_count(flist); i++ ) { +// +// double dist; +// struct imagefeature *f; +// Reflection *refl; +// RefListIterator *iter; +// +// f = image_get_feature(flist, i); +// if ( f == NULL ) continue; +// n_feat++; +// +// /* Find closest predicted peak */ +// +// for ( refl = first_refl(rlist, &iter); +// refl != NULL; +// refl = next_refl(refl, iter) ) +// { +// double fs, ss; +// get_detector_pos(refl, &fs, &ss); +// dist = sqrt(pow(fs-f->fs, 2.0) + pow(ss-f->ss, 2.0)); +// if ( dist < PEAK_CLOSE ) { +// n_sane++; +// continue; +// } +// } +// +// } +// +// if ( (float)n_sane / (float)n_feat < 0.1 ) return 0; +// +// return 1; +//} + + /* Integrate the list of predicted reflections in "image" */ void integrate_reflections(struct image *image, int polar, int use_closer, int bgsub) diff --git a/src/peaks.h b/src/peaks.h index 7b66b9ac..9d475ea9 100644 --- a/src/peaks.h +++ b/src/peaks.h @@ -27,7 +27,7 @@ extern void search_peaks(struct image *image, float threshold, extern void integrate_reflections(struct image *image, int polar, int use_closer, int bgsub); -extern int peak_sanity_check(RefList *rlist, ImageFeatureList *flist); +extern int peak_sanity_check(struct image * image); /* Exported so it can be poked by integration_check */ extern int integrate_peak(struct image *image, int cfs, int css, |