diff options
Diffstat (limited to 'libcrystfel/src/peaks.c')
-rw-r--r-- | libcrystfel/src/peaks.c | 284 |
1 files changed, 170 insertions, 114 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index bf6a92f0..6c09053b 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_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,27 @@ 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_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, + w, h, mask, image->crystals[i]); + } return mask; } @@ -236,6 +248,8 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( p == NULL ) return 1; if ( p->no_index ) return 1; + if ( saturated != NULL ) *saturated = 0; + /* Determine regions where there is expected to be a peak */ p_cfs = cfs - p->min_fs; p_css = css - p->min_ss; /* Panel-relative coordinates */ @@ -293,10 +307,9 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx]; - /* Veto peak if it contains saturation in bg region */ + /* Check if peak contains saturation in bg region */ if ( use_max_adu && (val > p->max_adu) ) { if ( saturated != NULL ) *saturated = 1; - return 1; } bg_tot += val; @@ -353,10 +366,9 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx] - bg_mean; - /* Veto peak if it contains saturation */ + /* Check if peak contains saturation */ if ( use_max_adu && (val > p->max_adu) ) { if ( saturated != NULL ) *saturated = 1; - return 1; } pk_counts++; @@ -387,7 +399,8 @@ static int integrate_peak(struct image *image, int cfs, int css, static void search_peaks_in_panel(struct image *image, float threshold, float min_gradient, float min_snr, struct panel *p, - double ir_inn, double ir_mid, double ir_out) + double ir_inn, double ir_mid, double ir_out, + int use_saturated) { int fs, ss, stride; float *data; @@ -400,7 +413,7 @@ static void search_peaks_in_panel(struct image *image, float threshold, int nrej_dis = 0; int nrej_pro = 0; int nrej_fra = 0; - int nrej_bad = 0; + int nrej_fail = 0; int nrej_snr = 0; int nacc = 0; int ncull; @@ -419,13 +432,11 @@ static void search_peaks_in_panel(struct image *image, float threshold, double max; unsigned int did_something; int r; + int saturated; /* Overall threshold */ if ( data[fs+stride*ss] < threshold ) continue; - /* Immediate rejection of pixels above max_adu */ - if ( data[fs+stride*ss] > p->max_adu ) continue; - /* Get gradients */ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss]; dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss]; @@ -491,12 +502,11 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Centroid peak and get better coordinates. */ r = integrate_peak(image, mask_fs, mask_ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, NULL); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); - if ( r ) { - /* Bad region - don't detect peak */ - nrej_bad++; - continue; + if ( saturated ) { + image->num_saturated_peaks++; + if ( !use_saturated ) continue; } /* It is possible for the centroid to fall outside the image */ @@ -518,6 +528,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, continue; } + if ( r ) { + /* Bad region - don't detect peak */ + nrej_fail++; + continue; + } + /* Add using "better" coordinates */ image_add_feature(image->features, f_fs, f_ss, image, intensity, NULL); @@ -535,10 +551,13 @@ static void search_peaks_in_panel(struct image *image, float threshold, ncull = 0; } -// STATUS("%i accepted, %i box, %i proximity, %i outside panel, " -// "%i in bad regions, %i with SNR < %g, %i badrow culled.\n", -// nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, -// nrej_snr, min_snr, ncull); + image->num_peaks += nacc; + + //STATUS("%i accepted, %i box, %i proximity, %i outside panel, " + // "%i failed integration, %i with SNR < %g, %i badrow culled, " + // "%i saturated.\n", + // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_fail, + // nrej_snr, min_snr, ncull, nrej_sat); if ( ncull != 0 ) { STATUS("WARNING: %i peaks were badrow culled. This feature" @@ -549,7 +568,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, void search_peaks(struct image *image, float threshold, float min_gradient, - float min_snr, double ir_inn, double ir_mid, double ir_out) + float min_snr, double ir_inn, double ir_mid, double ir_out, + int use_saturated) { int i; @@ -557,6 +577,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient, image_feature_list_free(image->features); } image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; for ( i=0; i<image->det->n_panels; i++ ) { @@ -564,35 +586,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient, if ( p->no_index ) continue; search_peaks_in_panel(image, threshold, min_gradient, - min_snr, p, ir_inn, ir_mid, ir_out); + min_snr, p, ir_inn, ir_mid, ir_out, + use_saturated); } } -double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) +int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { - 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); @@ -602,38 +615,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<n_cryst; j++ ) { + + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + + cell_get_cartesian(crystal_get_cell(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; } @@ -694,40 +711,29 @@ 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) +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, + int res_cutoff) { + 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); - il = sort_reflections(image->reflections, image->indexed_cell, &n); + if ( num_reflections(reflections) == 0 ) return; + + 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; @@ -794,10 +800,19 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, &intensity, &sigma, ir_inn, ir_mid, ir_out, 1, bgMasks[pnum], &saturated); - if ( saturated ) image->n_saturated++; + if ( saturated ) { + n_saturated++; + if ( !integrate_saturated ) r = 1; + } - /* I/sigma(I) cutoff */ - if ( intensity/sigma < min_snr ) r = 1; + /* I/sigma(I) cutoff + * Rejects reflections below --min-integration-snr, or if the + * SNR is clearly silly. Silly indicates that the intensity + * was zero. */ + snr = fabs(intensity)/sigma; + if ( isnan(snr) || (snr < min_snr) ) { + r = 1; + } /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { @@ -808,7 +823,6 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, set_redundancy(refl, 0); } - snr = intensity / sigma; if ( snr > 1.0 ) { if ( first ) { av = snr; @@ -818,27 +832,61 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } //STATUS("%5.2f A, %5.2f, av %5.2f\n", // 1e10/il[i].res, snr, av); - //if ( av < 1.0 ) break; + if ( res_cutoff && (av < 1.0) ) break; } } + 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 res_cutoff) +{ + 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++ ) { - free(bgMasks[i]); + int *mask; + mask = make_BgMask(image, &image->det->panels[i], ir_inn); + if ( mask == NULL ) { + ERROR("Couldn't create background mask.\n"); + return; + } + bgMasks[i] = mask; } - free(bgMasks); - image->diffracting_resolution = 0.0; + 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, res_cutoff); + } - free(il); + for ( i=0; i<image->det->n_panels; i++ ) { + free(bgMasks[i]); + } + free(bgMasks); } void validate_peaks(struct image *image, double min_snr, - int ir_inn, int ir_mid, int ir_out) + int ir_inn, int ir_mid, int ir_out, int use_saturated) { int i, n; ImageFeatureList *flist; - int n_wtf, n_int, n_dft, n_snr, n_prx; + int n_wtf, n_int, n_dft, n_snr, n_prx, n_sat; flist = image_feature_list_new(); if ( flist == NULL ) return; @@ -846,7 +894,7 @@ void validate_peaks(struct image *image, double min_snr, n = image_feature_count(image->features); /* Loop over peaks, putting each one through the integrator */ - n_wtf = 0; n_int = 0; n_dft = 0; n_snr = 0; n_prx = 0; + n_wtf = 0; n_int = 0; n_dft = 0; n_snr = 0; n_prx = 0; n_sat = 0; for ( i=0; i<n; i++ ) { struct imagefeature *f; @@ -856,6 +904,7 @@ void validate_peaks(struct image *image, double min_snr, double f_fs, f_ss; double intensity, sigma; struct panel *p; + int saturated; f = image_get_feature(image->features, i); if ( f == NULL ) { @@ -871,12 +920,17 @@ void validate_peaks(struct image *image, double min_snr, r = integrate_peak(image, f->fs, f->ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, NULL); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); if ( r ) { n_int++; continue; } + if ( saturated ) { + n_sat++; + if ( !use_saturated ) continue; + } + /* It is possible for the centroid to fall outside the image */ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs) || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) @@ -903,9 +957,11 @@ void validate_peaks(struct image *image, double min_snr, } //STATUS("HDF5: %i peaks, validated: %i. WTF: %i, integration: %i," - // " drifted: %i, SNR: %i, proximity: %i\n", + // " drifted: %i, SNR: %i, proximity: %i, saturated: %i\n", // n, image_feature_count(flist), - // n_wtf, n_int, n_dft, n_snr, n_prx); + // n_wtf, n_int, n_dft, n_snr, n_prx, n_sat); image_feature_list_free(image->features); image->features = flist; + image->num_saturated_peaks = n_sat; + image->num_peaks = image_feature_count(flist); } |