indexamajig: Add --use-saturated

5 files changed, 45 insertions, 18 deletions

diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 7ab194c7..6bf956e2 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -385,7 +385,8 @@ static int integrate_peak(struct image *image, int cfs, int css,
 static void search_peaks_in_panel(struct image *image, float threshold,
                                   float min_gradient, float min_snr,
                                   struct panel *p,
-                                  double ir_inn, double ir_mid, double ir_out)
+                                  double ir_inn, double ir_mid, double ir_out,
+                                  int use_saturated)
 {
 	int fs, ss, stride;
 	float *data;
@@ -400,6 +401,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 	int nrej_fra = 0;
 	int nrej_bad = 0;
 	int nrej_snr = 0;
+	int nrej_sat = 0;
 	int nacc = 0;
 	int ncull;
 
@@ -417,6 +419,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		double max;
 		unsigned int did_something;
 		int r;
+		int saturated;
 
 		/* Overall threshold */
 		if ( data[fs+stride*ss] < threshold ) continue;
@@ -489,7 +492,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		/* Centroid peak and get better coordinates. */
 		r = integrate_peak(image, mask_fs, mask_ss,
 		                   &f_fs, &f_ss, &intensity, &sigma,
-		                   ir_inn, ir_mid, ir_out, 0, NULL, NULL);
+		                   ir_inn, ir_mid, ir_out, 0, NULL, &saturated);
 
 		if ( r ) {
 			/* Bad region - don't detect peak */
@@ -516,6 +519,11 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 			continue;
 		}
 
+		if ( saturated && !use_saturated ) {
+			nrej_sat++;
+			continue;
+		}
+
 		/* Add using "better" coordinates */
 		image_add_feature(image->features, f_fs, f_ss, image, intensity,
 		                  NULL);
@@ -533,10 +541,11 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		ncull = 0;
 	}
 
-//	STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
-//	       "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
-//	       nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
-//	       nrej_snr, min_snr, ncull);
+	//STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
+	//       "%i in bad regions, %i with SNR < %g, %i badrow culled, "
+	//        "%i saturated.\n",
+	//       nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
+	//       nrej_snr, min_snr, ncull, nrej_sat);
 
 	if ( ncull != 0 ) {
 		STATUS("WARNING: %i peaks were badrow culled.  This feature"
@@ -547,7 +556,8 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 
 
 void search_peaks(struct image *image, float threshold, float min_gradient,
-                  float min_snr, double ir_inn, double ir_mid, double ir_out)
+                  float min_snr, double ir_inn, double ir_mid, double ir_out,
+                  int use_saturated)
 {
 	int i;
 
@@ -562,7 +572,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 
 		if ( p->no_index ) continue;
 		search_peaks_in_panel(image, threshold, min_gradient,
-		                      min_snr, p, ir_inn, ir_mid, ir_out);
+		                      min_snr, p, ir_inn, ir_mid, ir_out,
+		                      use_saturated);
 
 	}
 }
@@ -838,11 +849,11 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 
 
 void validate_peaks(struct image *image, double min_snr,
-                    int ir_inn, int ir_mid, int ir_out)
+                    int ir_inn, int ir_mid, int ir_out, int use_saturated)
 {
 	int i, n;
 	ImageFeatureList *flist;
-	int n_wtf, n_int, n_dft, n_snr, n_prx;
+	int n_wtf, n_int, n_dft, n_snr, n_prx, n_sat;
 
 	flist = image_feature_list_new();
 	if ( flist == NULL ) return;
@@ -850,7 +861,7 @@ void validate_peaks(struct image *image, double min_snr,
 	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;
+	n_wtf = 0;  n_int = 0;  n_dft = 0;  n_snr = 0;  n_prx = 0;  n_sat = 0;
 	for ( i=0; i<n; i++ ) {
 
 		struct imagefeature *f;
@@ -860,6 +871,7 @@ void validate_peaks(struct image *image, double min_snr,
 		double f_fs, f_ss;
 		double intensity, sigma;
 		struct panel *p;
+		int saturated;
 
 		f = image_get_feature(image->features, i);
 		if ( f == NULL ) {
@@ -875,7 +887,7 @@ void validate_peaks(struct image *image, double min_snr,
 
 		r = integrate_peak(image, f->fs, f->ss,
 		                   &f_fs, &f_ss, &intensity, &sigma,
-		                   ir_inn, ir_mid, ir_out, 0, NULL, NULL);
+		                   ir_inn, ir_mid, ir_out, 0, NULL, &saturated);
 		if ( r ) {
 			n_int++;
 			continue;
@@ -901,15 +913,20 @@ void validate_peaks(struct image *image, double min_snr,
 			continue;
 		}
 
+		if ( saturated && !use_saturated ) {
+			n_sat++;
+			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",
+	//       " drifted: %i, SNR: %i, proximity: %i, saturated: %i\n",
 	//       n, image_feature_count(flist),
-	//       n_wtf, n_int, n_dft, n_snr, n_prx);
+	//       n_wtf, n_int, n_dft, n_snr, n_prx, n_sat);
 	image_feature_list_free(image->features);
 	image->features = flist;
 }
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index ec76fd3b..3fae13fb 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -39,7 +39,8 @@
 
 extern void search_peaks(struct image *image, float threshold,
                          float min_gradient, float min_snr,
-                         double ir_inn, double ir_mid, double ir_out);
+                         double ir_inn, double ir_mid, double ir_out,
+                         int use_saturated);
 
 extern void integrate_reflections(struct image *image,
                                   int use_closer, int bgsub, double min_snr,
@@ -55,6 +56,7 @@ 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);
+                           int ir_inn, int ir_mid, int ir_out,
+                           int use_saturated);
 
 #endif	/* PEAKS_H */
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 80bb78e0..aed0b04b 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -292,13 +292,15 @@ static void process_image(const struct index_args *iargs,
 			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);
+		               iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+		               iargs->use_saturated);
 		break;
 
 		case PEAK_ZAEF:
 		search_peaks(&image, iargs->threshold,
 		             iargs->min_gradient, iargs->min_snr,
-		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out);
+		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+		             iargs->use_saturated);
 		break;
 
 	}
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 9f23f88b..50dc599b 100644
--- a/src/im-sandbox.h
+++ b/src/im-sandbox.h
@@ -67,6 +67,7 @@ struct index_args
 	double ir_out;
 	struct copy_hdf5_field *copyme;
 	int integrate_saturated;
+	int use_saturated;
 };
 
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 2e5c8904..57c16932 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -161,6 +161,8 @@ static void show_help(const char *s)
 "                           least 10%% of the located peaks.\n"
 "     --no-bg-sub          Don't subtract local background estimates from\n"
 "                           integrated intensities.\n"
+"     --use-saturated       During the initial peak search, don't reject\n"
+"                            peaks which contain pixels above max_adu.\n"
 "     --integrate-saturated  During the final integration stage, don't reject\n"
 "                            peaks which contain pixels above max_adu.\n"
 );
@@ -245,6 +247,7 @@ int main(int argc, char *argv[])
 	float ir_mid = 5.0;
 	float ir_out = 7.0;
 	int integrate_saturated = 0;
+	int use_saturated = 0;
 
 	copyme = new_copy_hdf5_field_list();
 	if ( copyme == NULL ) {
@@ -292,6 +295,7 @@ int main(int argc, char *argv[])
 		{"tolerance",          1, NULL,               13},
 		{"int-radius",         1, NULL,               14},
 		{"integrate-saturated",0, &integrate_saturated,1},
+		{"use-saturated",0, &use_saturated,            1},
 		{0, 0, NULL, 0}
 	};
 
@@ -638,6 +642,7 @@ int main(int argc, char *argv[])
 	iargs.ir_inn = ir_inn;
 	iargs.ir_mid = ir_mid;
 	iargs.ir_out = ir_out;
+	iargs.use_saturated = use_saturated;
 	iargs.integrate_saturated = integrate_saturated;
 
 	create_sandbox(&iargs, n_proc, prefix, config_basename, fh,