Hooks for estimation of maximum resolution for each image

8 files changed, 94 insertions, 0 deletions

diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt
index 1d5b4016..576c4a27 100644
--- a/doc/reference/CrystFEL-sections.txt
+++ b/doc/reference/CrystFEL-sections.txt
@@ -59,6 +59,7 @@ write_reflections_to_file
 read_reflections
 read_reflections_from_file
 asymmetric_indices
+res_cutoff
 </SECTION>
 
 <SECTION>
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 4d4ac2d7..5f11b425 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -84,6 +84,7 @@ typedef struct _imagefeaturelist ImageFeatureList;
  *    double                  osf;
  *    double                  profile_radius;
  *    int                     pr_dud;
+ *    double                  diffracting_resolution;
  *
  *    int                     width;
  *    int                     height;
@@ -150,6 +151,7 @@ struct image {
 	double                  osf;           /* Overall scaling factor */
 	double                  profile_radius; /* Radius of reflection */
 	int                     pr_dud;        /* Post refinement failed */
+	double                  diffracting_resolution;  /* Max 1/d in m^-1 */
 
 	int                     width;
 	int                     height;
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 82ce1cf3..0d72e514 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -290,6 +290,34 @@ int integrate_peak(struct image *image, int cfs, int css,
 }
 
 
+void estimate_resolution(RefList *list, UnitCell *cell,
+                         double *min, double *max)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double one_over_d;
+		signed int h, k, l;
+
+		get_indices(refl, &h, &k, &l);
+		one_over_d = 2.0 * resolution(cell, h, k, l);
+
+		if ( one_over_d > *max ) *max = one_over_d;
+		if ( one_over_d < *min ) *min = one_over_d;
+
+		/* FIXME: Implement this */
+
+	}
+
+	*min = 0.0;
+	*max = 0.0;
+}
+
+
 static void search_peaks_in_panel(struct image *image, float threshold,
                                   float min_gradient, float min_snr,
                                   struct panel *p)
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 9d475ea9..88dab26d 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -35,4 +35,7 @@ extern int integrate_peak(struct image *image, int cfs, int css,
                           double *pbg, double *pmax, double *sigma,
                           int do_polar, int centroid, int bgsub);
 
+extern void estimate_resolution(RefList *list, UnitCell *cell,
+                                double *min, double *max);
+
 #endif	/* PEAKS_H */
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index c13b7bcf..66bd377e 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -434,3 +434,40 @@ double max_intensity(RefList *list)
 
 	return max;
 }
+
+
+/**
+ * res_cutoff:
+ * @list: A %RefList
+ *
+ * Returns: A new %RefList with resolution cutoff applied
+ **/
+RefList *res_cutoff(RefList *list, UnitCell *cell, double min, double max)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *new;
+
+	new = reflist_new();
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double one_over_d;
+		signed int h, k, l;
+		Reflection *n;
+
+		get_indices(refl, &h, &k, &l);
+
+		one_over_d = 2.0 * resolution(cell, h, k, l);
+		if ( one_over_d < min ) continue;
+		if ( one_over_d > max ) continue;
+
+		n = add_refl(new, h, k, l);
+		copy_data(n, refl);
+	}
+
+	reflist_free(list);
+	return new;
+}
diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h
index c451954b..2da0dec1 100644
--- a/libcrystfel/src/reflist-utils.h
+++ b/libcrystfel/src/reflist-utils.h
@@ -45,4 +45,7 @@ extern void resolution_limits(RefList *list, UnitCell *cell,
 
 extern double max_intensity(RefList *list);
 
+extern RefList *res_cutoff(RefList *list, UnitCell *cell,
+                           double min, double max);
+
 #endif	/* REFLIST_UTILS_H */
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index ce948109..b713e399 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -241,6 +241,15 @@ void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f)
 	fprintf(ofh, "photon_energy_eV = %f\n",
 	        J_to_eV(ph_lambda_to_en(i->lambda)));
 
+	if ( i->reflections != NULL ) {
+
+		fprintf(ofh, "diffraction_resolution_limit"
+		             " = %.2f nm^-1 or %.2f A\n",
+		             i->diffracting_resolution/1e9,
+		             1e9 / i->diffracting_resolution);
+
+	}
+
 	if ( i->det != NULL ) {
 
 		int j;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index a742eb65..361d0514 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -358,9 +358,20 @@ static void process_image(void *pp, int cookie)
 		pargs->indexable = 1;
 
 		if ( image.reflections != NULL ) {
+
+			double min, max;
+
 			integrate_reflections(&image, config_polar,
 					      pargs->static_args.config_closer,
 					      pargs->static_args.config_bgsub);
+
+			estimate_resolution(image.reflections,
+			                    image.indexed_cell, &min, &max);
+			image.reflections = res_cutoff(image.reflections,
+			                               image.indexed_cell,
+			                               min, max);
+			image.diffracting_resolution = max;
+
 		}
 
 	} else {