diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/peaks.c | 199 | ||||
-rw-r--r-- | libcrystfel/src/peaks.h | 3 |
2 files changed, 113 insertions, 89 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index dabc3b1a..c350388b 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -3,12 +3,12 @@ * * Peak search and other image analysis * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de> * 2011 Andrew Martin <andrew.martin@desy.de> * 2011 Richard Kirian @@ -147,24 +147,15 @@ static int cull_peaks(struct image *image) } -/* cfs, css relative to panel origin */ -static int *make_BgMask(struct image *image, struct panel *p, - double ir_out, double ir_inn) +static void add_crystal_to_mask(struct image *image, struct panel *p, + double ir_out, double ir_inn, int w, int h, + int *mask, Crystal *cr) { Reflection *refl; RefListIterator *iter; - int *mask; - int w, h; - - w = p->max_fs - p->min_fs + 1; - h = p->max_ss - p->min_ss + 1; - mask = calloc(w*h, sizeof(int)); - if ( mask == NULL ) return NULL; - - if ( image->reflections == NULL ) return mask; /* Loop over all reflections */ - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -205,6 +196,28 @@ static int *make_BgMask(struct image *image, struct panel *p, } } +} + + +/* cfs, css relative to panel origin */ +static int *make_BgMask(struct image *image, struct panel *p, + double ir_out, double ir_inn) +{ + int *mask; + int w, h; + int i; + + w = p->max_fs - p->min_fs + 1; + h = p->max_ss - p->min_ss + 1; + mask = calloc(w*h, sizeof(int)); + if ( mask == NULL ) return NULL; + + if ( image->crystals == NULL ) return mask; + + for ( i=0; i<image->n_crystals; i++ ) { + add_crystal_to_mask(image, p, ir_inn, ir_out, + w, h, mask, image->crystals[i]); + } return mask; } @@ -581,29 +594,19 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } -double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) +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; - double stot = 0.0; - - /* Round towards nearest */ - fesetround(1); - - /* Cell basis vectors for this image */ - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + int i; + const double min_dist = 0.25; - /* 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; + int j; /* Assume all image "features" are genuine peaks */ f = image_get_feature(image->features, i); @@ -613,38 +616,42 @@ double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) /* 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 */ - if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist) - && (fabs(l - ld) < min_dist) ) - { - double sval; - n_sane++; - sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); - stot += 1.0 - sval; - continue; - } + for ( j=0; j<image->n_crystals; j++ ) { + + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + + cell_get_cartesian(crystal_get_cell(image->crystals[j]), + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + /* 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 */ + if ( (fabs(h - hd) < min_dist) + && (fabs(k - kd) < min_dist) + && (fabs(l - ld) < min_dist) ) + { + n_sane++; + continue; + } - } + } - *pst = stot; - return (double)n_sane / (float)n_feat; -} + } -int peak_sanity_check(struct image *image) -{ - double stot; /* 0 means failed test, 1 means passed test */ - return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5; + return ((double)n_sane / n_feat) >= 0.5; } @@ -705,43 +712,28 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell, } -/* Integrate the list of predicted reflections in "image" */ -void integrate_reflections(struct image *image, int use_closer, int bgsub, - double min_snr, - double ir_inn, double ir_mid, double ir_out, - int integrate_saturated) +static void integrate_crystal(Crystal *cr, struct image *image, int use_closer, + int bgsub, double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int **bgMasks) { + RefList *reflections; struct integr_ind *il; int n, i; double av = 0.0; int first = 1; - int **bgMasks; + int n_saturated = 0; + + reflections = crystal_get_reflections(cr); - if ( num_reflections(image->reflections) == 0 ) return; + if ( num_reflections(reflections) == 0 ) return; - il = sort_reflections(image->reflections, image->indexed_cell, &n); + il = sort_reflections(reflections, crystal_get_cell(cr), &n); if ( il == NULL ) { ERROR("Couldn't sort reflections\n"); return; } - /* Make background masks for all panels */ - bgMasks = calloc(image->det->n_panels, sizeof(int *)); - if ( bgMasks == NULL ) { - ERROR("Couldn't create list of background masks.\n"); - return; - } - for ( i=0; i<image->det->n_panels; i++ ) { - int *mask; - mask = make_BgMask(image, &image->det->panels[i], - ir_out, ir_inn); - if ( mask == NULL ) { - ERROR("Couldn't create background mask.\n"); - return; - } - bgMasks[i] = mask; - } - for ( i=0; i<n; i++ ) { double fs, ss, intensity; @@ -809,7 +801,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, 1, bgMasks[pnum], &saturated); if ( saturated ) { - image->n_saturated++; + n_saturated++; if ( !integrate_saturated ) r = 1; } @@ -844,14 +836,49 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } } + crystal_set_num_saturated_reflections(cr, n_saturated); + crystal_set_resolution_limit(cr, 0.0); + + free(il); +} + + +/* Integrate the list of predicted reflections in "image" */ +void integrate_reflections(struct image *image, int use_closer, int bgsub, + double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated) +{ + int i; + int **bgMasks; + + /* Make background masks for all panels */ + bgMasks = calloc(image->det->n_panels, sizeof(int *)); + if ( bgMasks == NULL ) { + ERROR("Couldn't create list of background masks.\n"); + return; + } + for ( i=0; i<image->det->n_panels; i++ ) { + int *mask; + mask = make_BgMask(image, &image->det->panels[i], + ir_out, ir_inn); + if ( mask == NULL ) { + ERROR("Couldn't create background mask.\n"); + return; + } + bgMasks[i] = mask; + } + + for ( i=0; i<image->n_crystals; i++ ) { + integrate_crystal(image->crystals[i], image, use_closer, + bgsub, min_snr, ir_inn, ir_mid, ir_out, + integrate_saturated, bgMasks); + } + for ( i=0; i<image->det->n_panels; i++ ) { free(bgMasks[i]); } free(bgMasks); - - image->diffracting_resolution = 0.0; - - free(il); } diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 3fae13fb..faee54f6 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -47,9 +47,6 @@ extern void integrate_reflections(struct image *image, double ir_inn, double ir_mid, double ir_out, int integrate_saturated); -extern double peak_lattice_agreement(struct image *image, UnitCell *cell, - double *pst); - extern int peak_sanity_check(struct image *image); extern void estimate_resolution(RefList *list, UnitCell *cell, |