From 25cbb53d692a19e5087f34bcb018ced840421639 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 30 May 2013 10:21:52 -0700 Subject: process_hkl: Add --min-snr --- src/process_hkl.c | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) (limited to 'src') 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= Require at least measurements before a\n" " reflection appears in the output. Default: 2\n" +" --min-snr= 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 ) { -- cgit v1.2.3