process_hkl: Add --min-snr

1 files changed, 25 insertions, 10 deletions

diff --git a/src/process_hkl.c b/src/process_hkl.c
index fda02137..730fa14d 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -79,6 +79,8 @@ static void show_help(const char *s)
 "      --no-polarisation     Disable polarisation correction.\n"
 "      --min-measurements=<n> Require at least <n> measurements before a\n"
 "                             reflection appears in the output.  Default: 2\n"
+"      --min-snr=<n>         Require individual intensity measurements to\n"
+"                             have I > n * sigma(I).  Default: -infinity.\n"
 );
 }
 
@@ -177,7 +179,7 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
                         RefList *reference, const SymOpList *sym,
                         double *hist_vals, signed int hist_h,
                         signed int hist_k, signed int hist_l, int *hist_n,
-                        int config_nopolar)
+                        int config_nopolar, double min_snr)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -208,6 +210,12 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
 		double w;
 		double temp, delta, R, mean, M2, sumweight;
 
+		refl_intensity = scale * get_intensity(refl);
+		refl_sigma = scale * get_esd_intensity(refl);
+		w = 1.0;//pow(refl_sigma, -2.0);
+
+		if ( refl_intensity < min_snr * refl_sigma ) continue;
+
 		get_indices(refl, &h, &k, &l);
 
 		/* Put into the asymmetric unit for the target group */
@@ -218,10 +226,6 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
 			model_version = add_refl(model, h, k, l);
 		}
 
-		refl_intensity = scale * get_intensity(refl);
-		refl_sigma = scale * get_esd_intensity(refl);
-		w = 1.0;//pow(refl_sigma, -2.0);
-
 		mean = get_intensity(model_version);
 		sumweight = get_temp1(model_version);
 		M2 = get_temp2(model_version);
@@ -271,7 +275,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
                      double *hist_vals, signed int hist_h,
                      signed int hist_k, signed int hist_l,
                      int *hist_i, int config_nopolar, int min_measurements,
-                     int start_after, int stop_after)
+                     double min_snr, int start_after, int stop_after)
 {
 	int rval;
 	int n_images = 0;
@@ -305,7 +309,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
 			n_crystals++;
 			r = merge_crystal(model, &image, cr, reference, sym,
 			                  hist_vals, hist_h, hist_k, hist_l,
-			                  hist_i, config_nopolar);
+			                  hist_i, config_nopolar, min_snr);
 
 			if ( r == 0 ) n_crystals_used++;
 
@@ -375,6 +379,7 @@ int main(int argc, char *argv[])
 	int r;
 	int start_after = 0;
 	int stop_after = 0;
+	double min_snr = -INFINITY;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -393,6 +398,7 @@ int main(int argc, char *argv[])
 		{"histogram",          1, NULL,               'g'},
 		{"hist-parameters",    1, NULL,               'z'},
 		{"min-measurements",   1, NULL,                2},
+		{"min-snr",            1, NULL,                3},
 		{0, 0, NULL, 0}
 	};
 
@@ -446,13 +452,22 @@ int main(int argc, char *argv[])
 
 			case 2 :
 			errno = 0;
-			min_measurements = strtod(optarg, &rval);
+			min_measurements = strtol(optarg, &rval, 10);
 			if ( *rval != '\0' ) {
 				ERROR("Invalid value for --min-measurements.\n");
 				return 1;
 			}
 			break;
 
+			case 3 :
+			errno = 0;
+			min_snr = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --min-snr.\n");
+				return 1;
+			}
+			break;
+
 			case '?' :
 			break;
 
@@ -533,7 +548,7 @@ int main(int argc, char *argv[])
 
 	hist_i = 0;
 	r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l,
-	              &hist_i, config_nopolar, min_measurements,
+	              &hist_i, config_nopolar, min_measurements, min_snr,
 	              start_after, stop_after);
 	fprintf(stderr, "\n");
 	if ( r ) {
@@ -561,7 +576,7 @@ int main(int argc, char *argv[])
 
 			r = merge_all(st, model, reference, sym,
 				     hist_vals, hist_h, hist_k, hist_l, &hist_i,
-				     config_nopolar, min_measurements,
+				     config_nopolar, min_measurements, min_snr,
 				     start_after, stop_after);
 			fprintf(stderr, "\n");
 			if ( r ) {