indexamajig: Estimate resolution based on peaks only

4 files changed, 55 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);