diff options
| -rw-r--r-- | libcrystfel/src/image.h | 6 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.c | 46 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.h | 4 | ||||
| -rw-r--r-- | libcrystfel/src/stream.c | 2 | ||||
| -rwxr-xr-x | scripts/ave-peak-resolution | 46 | ||||
| -rw-r--r-- | src/process_image.c | 3 |
6 files changed, 104 insertions, 3 deletions
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 9769c66a..40717ab3 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -3,11 +3,11 @@ * * Handle images and image features * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2017 Thomas White <taw@physics.org> + * 2009-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * @@ -221,6 +221,8 @@ struct image { /* Detected peaks */ long long num_peaks; long long num_saturated_peaks; + double peak_resolution; /* Estimate of resolution + * based on peaks only */ ImageFeatureList *features; }; diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 6c48bd3a..3d7e8189 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -796,3 +796,49 @@ void validate_peaks(struct image *image, double min_snr, image->num_saturated_peaks = n_sat; image->num_peaks = image_feature_count(flist); } + + +static int compare_double(const void *av, const void *bv) +{ + double a = *(double *)av; + double b = *(double *)bv; + if ( a > b ) return 1; + if ( a < b ) return -1; + return 0; +} + + +double estimate_peak_resolution(ImageFeatureList *peaks, double lambda) +{ + int i, n, j; + double *rns; + double max_res; + + n = image_feature_count(peaks); + rns = malloc(n*sizeof(double)); + if ( rns == NULL ) return -1.0; + + /* Get resolution values for all peaks */ + j = 0; + for ( i=0; i<n; i++ ) { + + struct imagefeature *f; + struct rvec r; + + f = image_get_feature(peaks, i); + if ( f == NULL ) continue; + + r = get_q_for_panel(f->p, f->fs, f->ss, NULL, 1.0/lambda); + rns[j++] = modulus(r.u, r.v, r.w); + + } + + /* Slightly horrible outlier removal */ + qsort(rns, j, sizeof(double), compare_double); + n = j/50; + if ( n < 2 ) n = 2; + max_res = rns[(j-1)-n]; + + free(rns); + return max_res; +} diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 3abef729..a82d8d25 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -7,7 +7,7 @@ * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * 2017 Valerio Mariani <valerio.mariani@desy.de> * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de> * @@ -70,6 +70,8 @@ extern void validate_peaks(struct image *image, double min_snr, int ir_inn, int ir_mid, int ir_out, int use_saturated, int check_snr); +extern double estimate_peak_resolution(ImageFeatureList *peaks, double lambda); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index 9d97ff98..c11e38df 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -881,6 +881,8 @@ int write_chunk(Stream *st, struct image *i, struct imagefile *imfile, fprintf(st->fh, "num_peaks = %lli\n", i->num_peaks); fprintf(st->fh, "num_saturated_peaks = %lli\n", i->num_saturated_peaks); + fprintf(st->fh, "peak_resolution = %f nm^-1 or %f A\n", + i->peak_resolution/1e9, 1e10/i->peak_resolution); if ( include_peaks ) { if ( AT_LEAST_VERSION(st, 2, 3) ) { ret = write_peaks_2_3(i, st->fh); diff --git a/scripts/ave-peak-resolution b/scripts/ave-peak-resolution new file mode 100755 index 00000000..468cffd9 --- /dev/null +++ b/scripts/ave-peak-resolution @@ -0,0 +1,46 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Find mean diffracting resolution +# +# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Author: +# 2014-2017 Thomas White <taw@physics.org> +# + +import sys +import numpy +import matplotlib.pyplot as plt + +f = open(sys.argv[1]) +a = [] +pulseId = 0 + +while True: + fline = f.readline() + if not fline: + break + if fline.find("hdf5/XFEL/pulseId") != -1: + pulseId = int(fline.split('= ')[1]) + if fline.find("peak_resolution") != -1: + res = float(fline.split('= ')[1].split(' ')[0].rstrip("\r\n")) + if pulseId != 0 and pulseId != 488 and pulseId != 496: + a.append(res) + continue + +f.close() + +b = numpy.array(a) +print(" Mean: {:.2} nm^-1 = {:.2} A".format(numpy.mean(b),10.0/numpy.mean(b))) +print(" Best: {:.2} nm^-1 = {:.2} A".format(numpy.max(b),10.0/numpy.max(b))) +print("Worst: {:.2} nm^-1 = {:.2} A".format(numpy.min(b),10.0/numpy.min(b))) +print("Std deviation: {:.2} nm^-1".format(numpy.std(b))) + +plt.hist(a,bins=30) +plt.title('Resolution based on peak search') +plt.xlabel('Resolution / nm^-1') +plt.ylabel('Frequency') +plt.grid(True) +plt.show() diff --git a/src/process_image.c b/src/process_image.c index 1d41ae6e..b0a03e46 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -245,6 +245,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } + image.peak_resolution = estimate_peak_resolution(image.features, + image.lambda); + restore_image_data(image.dp, image.det, prefilter); rn = getcwd(NULL, 0); |