aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2019-08-23 12:21:14 +0200
committerThomas White <taw@physics.org>2019-08-26 14:13:57 +0200
commita403d930065382247f83fd14d5f8eb88504ad35b (patch)
tree601ff3d64b10821ebfbb648bbde257a725c59667
parentb6af09d78d49d552a05b81bd20a7854a180373f5 (diff)
Fix resolution estimation
Didn't work when the number of peaks was very low.
-rw-r--r--libcrystfel/src/peaks.c24
1 files changed, 13 insertions, 11 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 6da79e4f..f7433383 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -812,34 +812,36 @@ static int compare_double(const void *av, const void *bv)
double estimate_peak_resolution(ImageFeatureList *peaks, double lambda)
{
- int i, n, j;
+ int i, npk, ncut;
double *rns;
double max_res;
- n = image_feature_count(peaks);
- rns = malloc(n*sizeof(double));
+ npk = image_feature_count(peaks);
+
+ /* No peaks -> no resolution! */
+ if ( npk == 0 ) return 0.0;
+
+ rns = malloc(npk*sizeof(double));
if ( rns == NULL ) return -1.0;
/* Get resolution values for all peaks */
- j = 0;
- for ( i=0; i<n; i++ ) {
+ for ( i=0; i<npk; 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);
+ rns[i] = 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];
+ qsort(rns, npk, sizeof(double), compare_double);
+ ncut = npk/50;
+ if ( ncut < 2 ) ncut = 0;
+ max_res = rns[(npk-1)-ncut];
free(rns);
return max_res;