From b5afea064fba22f893fd32a5834bdb4017d96d30 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 22 Oct 2012 16:27:03 -0700
Subject: indexamajig: Validate (re-integrate, check bad regions etc) peaks
 from HDF5

---
 libcrystfel/src/peaks.c | 78 +++++++++++++++++++++++++++++++++++++++++++++++++
 libcrystfel/src/peaks.h |  3 ++
 src/im-sandbox.c        |  2 ++
 3 files changed, 83 insertions(+)

diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 503d534d..bf6a92f0 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -831,3 +831,81 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 
 	free(il);
 }
+
+
+void validate_peaks(struct image *image, double min_snr,
+                    int ir_inn, int ir_mid, int ir_out)
+{
+	int i, n;
+	ImageFeatureList *flist;
+	int n_wtf, n_int, n_dft, n_snr, n_prx;
+
+	flist = image_feature_list_new();
+	if ( flist == NULL ) return;
+
+	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;
+	for ( i=0; i<n; i++ ) {
+
+		struct imagefeature *f;
+		int r;
+		double d;
+		int idx;
+		double f_fs, f_ss;
+		double intensity, sigma;
+		struct panel *p;
+
+		f = image_get_feature(image->features, i);
+		if ( f == NULL ) {
+			n_wtf++;
+			continue;
+		}
+
+		p = find_panel(image->det, f->fs, f->ss);
+		if ( p == NULL ) {
+			n_wtf++;
+			continue;
+		}
+
+		r = integrate_peak(image, f->fs, f->ss,
+		                   &f_fs, &f_ss, &intensity, &sigma,
+		                   ir_inn, ir_mid, ir_out, 0, NULL, NULL);
+		if ( r ) {
+			n_int++;
+			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) )
+		{
+			n_dft++;
+			continue;
+		}
+
+		if ( fabs(intensity)/sigma < min_snr ) {
+			n_snr++;
+			continue;
+		}
+
+		/* Check for a nearby feature */
+		image_feature_closest(flist, f_fs, f_ss, &d, &idx);
+		if ( d < 2.0*ir_inn ) {
+			n_prx++;
+			continue;
+		}
+
+		/* Add using "better" coordinates */
+		image_add_feature(flist, f_fs, f_ss, image, intensity, NULL);
+
+	}
+
+	//STATUS("HDF5: %i peaks, validated: %i.  WTF: %i, integration: %i,"
+	//       " drifted: %i, SNR: %i, proximity: %i\n",
+	//       n, image_feature_count(flist),
+	//       n_wtf, n_int, n_dft, n_snr, n_prx);
+	image_feature_list_free(image->features);
+	image->features = flist;
+}
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index afc7b868..39fdf54c 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -53,4 +53,7 @@ extern int peak_sanity_check(struct image *image);
 extern void estimate_resolution(RefList *list, UnitCell *cell,
                                 double *min, double *max);
 
+extern void validate_peaks(struct image *image, double min_snr,
+                           int ir_inn, int ir_mid, int ir_out);
+
 #endif	/* PEAKS_H */
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index c313bda9..31dff291 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -291,6 +291,8 @@ static void process_image(const struct index_args *iargs,
 			iargs->hdf5_peak_path)) {
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
+		validate_peaks(&image, iargs->min_int_snr,
+		               iargs->ir_inn, iargs->ir_mid, iargs->ir_out);
 		break;
 
 		case PEAK_ZAEF:
-- 
cgit v1.2.3