Add option to limit peak SNR

4 files changed, 37 insertions, 12 deletions

diff --git a/src/indexamajig.c b/src/indexamajig.c
index 79bf69e6..26ea158d 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -68,6 +68,7 @@ struct static_index_args
 	int config_bgsub;
 	float threshold;
 	float min_gradient;
+	float min_snr;
 	struct detector *det;
 	IndexingMethod *indm;
 	IndexingPrivate **ipriv;
@@ -182,6 +183,8 @@ static void show_help(const char *s)
 "     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
 "     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
 "                           Default: 100,000.\n"
+"     --min-snr=<n>        Minimum signal-to-noise ratio for peaks.\n"
+"                           Default: 5.\n"
 " -e, --image=<element>    Use this image from the HDF5 file.\n"
 "                           Example: /data/data0.\n"
 "                           Default: The first one found.\n"
@@ -310,7 +313,8 @@ static void process_image(void *pp, int cookie)
 		break;
 	case PEAK_ZAEF :
 		search_peaks(&image, pargs->static_args.threshold,
-		             pargs->static_args.min_gradient);
+		             pargs->static_args.min_gradient,
+		             pargs->static_args.min_snr);
 		break;
 	}
 
@@ -506,6 +510,7 @@ int main(int argc, char *argv[])
 	int config_basename = 0;
 	float threshold = 800.0;
 	float min_gradient = 100000.0;
+	float min_snr = 5;
 	struct detector *det;
 	char *geometry = NULL;
 	IndexingMethod *indm;
@@ -563,6 +568,7 @@ int main(int argc, char *argv[])
 		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
 		{"threshold",          1, NULL,               't'},
 		{"min-gradient",       1, NULL,                4},
+		{"min-snr",            1, NULL,               11},
 		{"no-check-prefix",    0, &config_checkprefix, 0},
 		{"no-closer-peak",     0, &config_closer,      0},
 		{"insane",             0, &config_insane,      1},
@@ -641,6 +647,10 @@ int main(int argc, char *argv[])
 			min_gradient = strtof(optarg, NULL);
 			break;
 
+		case 11 :
+			min_snr = strtof(optarg, NULL);
+			break;
+
 		case 'e' :
 			element = strdup(optarg);
 			break;
@@ -899,6 +909,7 @@ int main(int argc, char *argv[])
 	qargs.static_args.cellr = cellr;
 	qargs.static_args.threshold = threshold;
 	qargs.static_args.min_gradient = min_gradient;
+	qargs.static_args.min_snr = min_snr;
 	qargs.static_args.det = det;
 	qargs.static_args.indm = indm;
 	qargs.static_args.ipriv = ipriv;
diff --git a/src/peaks.c b/src/peaks.c
index fd1b0983..92556ea6 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -281,7 +281,8 @@ int integrate_peak(struct image *image, int cfs, int css,
 
 
 static void search_peaks_in_panel(struct image *image, float threshold,
-                                  float min_gradient, struct panel *p)
+                                  float min_gradient, float min_snr,
+                                  struct panel *p)
 {
 	int fs, ss, stride;
 	float *data;
@@ -297,6 +298,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 	int nrej_pro = 0;
 	int nrej_fra = 0;
 	int nrej_bad = 0;
+	int nrej_snr = 0;
 	int nacc = 0;
 	int ncull;
 	const int pws = p->peak_sep/2;
@@ -404,8 +406,8 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 			continue;
 		}
 
-		if (intensity/sigma < 5) {
-			//printf("SNR: %g\n",intensity/sigma);
+		if (intensity/sigma < min_snr) {
+			nrej_snr++;
 			continue;
 		}
 
@@ -434,12 +436,13 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 	}
 
 //	STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
-//	       "%i in bad regions, %i badrow culled.\n",
-//	       nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, ncull);
+//	       "%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);
 }
 
 
-void search_peaks(struct image *image, float threshold, float min_gradient)
+void search_peaks(struct image *image, float threshold, float min_gradient,
+                  float min_snr)
 {
 	int i;
 
@@ -453,7 +456,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
 		struct panel *p = &image->det->panels[i];
 
 		if ( p->no_index ) continue;
-		search_peaks_in_panel(image, threshold, min_gradient, p);
+		search_peaks_in_panel(image, threshold, min_gradient, min_snr, p);
 
 	}
 }
diff --git a/src/peaks.h b/src/peaks.h
index 07c257e6..7b66b9ac 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -22,7 +22,7 @@
 #include "reflist.h"
 
 extern void search_peaks(struct image *image, float threshold,
-                         float min_gradient);
+                         float min_gradient, float min_snr);
 
 extern void integrate_reflections(struct image *image,
                                   int polar, int use_closer, int bgsub);
diff --git a/src/sum_stack.c b/src/sum_stack.c
index 2e5b138e..a6e7599a 100644
--- a/src/sum_stack.c
+++ b/src/sum_stack.c
@@ -49,6 +49,7 @@ struct sum_args
 	SumMethod sum_method;
 	double threshold;
 	double min_gradient;
+	double min_snr;
 	char *element;
 };
 
@@ -64,6 +65,7 @@ struct queue_args
 	SumMethod sum_method;
 	double threshold;
 	double min_gradient;
+	double min_snr;
 	char *element;
 
 	char *use_this_one_instead;
@@ -100,6 +102,8 @@ static void show_help(const char *s)
 "                           method.  Default: 400 adu\n"
 "      --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
 "                            Default: 100,000.\n"
+"      --min-snr=<n>        Minimum signal-to-noise ration for peaks.\n"
+"                            Default: 5.\n"
 " -e, --image=<element>    Use this image from the HDF5 file.\n"
 "                           Example: /data/data0.\n"
 "                           Default: The first one found.\n"
@@ -116,14 +120,14 @@ static void show_help(const char *s)
 
 
 static void sum_peaks(struct image *image, double *sum, double threshold,
-                      double min_gradient)
+                      double min_gradient, double min_snr)
 {
 	int i;
 	int w = image->width;
 	int h = image->height;
 	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
 
-	search_peaks(image, threshold, min_gradient);
+	search_peaks(image, threshold, min_gradient, min_snr);
 
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
@@ -231,7 +235,7 @@ static void add_image(void *args, int cookie)
 
 	case SUM_PEAKS :
 		sum_peaks(&image, pargs->sum, pargs->threshold,
-		          pargs->min_gradient);
+		          pargs->min_gradient, pargs->min_snr);
 		break;
 
 	}
@@ -261,6 +265,7 @@ static void *get_image(void *qp)
 	pargs->sum_method = qargs->sum_method;
 	pargs->threshold = qargs->threshold;
 	pargs->min_gradient = qargs->min_gradient;
+	pargs->min_snr = qargs->min_snr;
 	pargs->config_cmfilter = qargs->config_cmfilter;
 	pargs->config_noisefilter = qargs->config_noisefilter;
 	pargs->sum = qargs->sum;
@@ -312,6 +317,7 @@ int main(int argc, char *argv[])
 	char *intermediate = NULL;
 	double threshold = 400.0;
 	float min_gradient = 100000.0;
+	float min_snr = 5;
 	SumMethod sum;
 	int nthreads = 1;
 	struct queue_args qargs;
@@ -337,6 +343,7 @@ int main(int argc, char *argv[])
 		{"intermediate",       1, NULL,               'p'},
 		{"threshold",          1, NULL,               't'},
 		{"min-gradient",       1, NULL,                4},
+		{"min-snr",            1, NULL,               11},
 		{"basename",           0, &config_basename,    1},
 		{"image",              1, NULL,               'e'},
 		{0, 0, NULL, 0}
@@ -387,6 +394,10 @@ int main(int argc, char *argv[])
 			min_gradient = strtof(optarg, NULL);
 			break;
 
+		case 11 :
+			min_snr = strtof(optarg, NULL);
+			break;
+
 		case 0 :
 			break;