aboutsummaryrefslogtreecommitdiff
path: root/src/process_hkl.c
diff options
context:
space:
mode:
authorTakanori Nakane <nakane.t@gmail.com>2014-07-04 16:42:14 +0100
committerThomas White <taw@physics.org>2014-09-18 16:50:00 +0200
commit944c01964d87e9f3830a5079a78e0e5d4903a12e (patch)
tree4577d3f1ab329d430a9e6e2a71fc8852f47fbbc7 /src/process_hkl.c
parentd0443d87cd0fb6a85b1bd26d0746853af19c853e (diff)
Calculate CC for each frame. Reject bad frames (--min-cc). Output stats (--stat=)
But these should be in partialator?
Diffstat (limited to 'src/process_hkl.c')
-rw-r--r--src/process_hkl.c109
1 files changed, 102 insertions, 7 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index b6a2d8cb..b72b84f1 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -67,6 +67,7 @@ static void show_help(const char *s)
" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n"
" -o, --output=<filename> Specify output filename for merged intensities\n"
" Default: processed.hkl).\n"
+" --stat=<filename> Specify output filename for merging statistics.\n"
" -y, --symmetry=<sym> Merge according to point group <sym>.\n"
"\n"
" --start-after=<n> Skip <n> crystals at the start of the stream.\n"
@@ -83,6 +84,8 @@ static void show_help(const char *s)
" 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"
+" --min-cc=<n> Reject frames with CC less than n. Default: infinity.\n"
+" have I > n * sigma(I). Default: -infinity.\n"
" --max-adu=<n> Maximum peak value. Default: infinity.\n"
" --min-res=<n> Merge only crystals which diffract above <n> A.\n"
" --push-res=<n> Integrate higher than apparent resolution cutoff.\n"
@@ -179,6 +182,57 @@ static double scale_intensities(RefList *reference, RefList *new,
return s;
}
+static double cc_intensities(RefList *reference, RefList *new,
+ const SymOpList *sym)
+{
+ // x is reference
+ float s_xy = 0.0;
+ float s_x = 0.0;
+ float s_y = 0.0;
+ float s_x2 = 0.0;
+ float s_y2 = 0.0;
+ int n = 0;
+ float t1, t2;
+
+ Reflection *refl;
+ RefListIterator *iter;
+
+ for ( refl = first_refl(new, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+
+ double i1, i2;
+ signed int hu, ku, lu;
+ signed int h, k, l;
+ Reflection *reference_version;
+
+ get_indices(refl, &h, &k, &l);
+ get_asymm(sym, h, k, l, &hu, &ku, &lu);
+
+ reference_version = find_refl(reference, hu, ku, lu);
+ if ( reference_version == NULL ) continue;
+
+ i1 = get_intensity(reference_version);
+ i2 = get_intensity(refl);
+
+
+ s_xy += i1 * i2;
+ s_x += i1;
+ s_y += i2;
+ s_x2 += i1 * i1;
+ s_y2 += i2 * i2;
+ n++;
+
+ }
+
+ t1 = s_x2 - s_x*s_x / n;
+ t2 = s_y2 - s_y*s_y / n;
+
+ if ( (t1 <= 0.0) || (t2 <= 0.0) ) return 0.0;
+
+ return (s_xy - s_x*s_y/n) / sqrt(t1*t2);
+}
static double *check_hist_size(int n, double *hist_vals)
{
@@ -203,11 +257,11 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
double **hist_vals, signed int hist_h,
signed int hist_k, signed int hist_l, int *hist_n,
int config_nopolar, double min_snr, double max_adu,
- double push_res)
+ double push_res, double min_cc, int do_scale, FILE *stat)
{
Reflection *refl;
RefListIterator *iter;
- double scale;
+ double scale, cc;
/* First, correct for polarisation */
if ( !config_nopolar ) {
@@ -218,9 +272,16 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
if ( reference != NULL ) {
scale = scale_intensities(reference,
crystal_get_reflections(cr), sym);
+ cc = cc_intensities(reference, crystal_get_reflections(cr), sym);
+ if (stat != NULL) {
+ fprintf(stat, "%s,%f,%f\n",image->filename, scale, cc);
+ }
+ if (cc < min_cc || scale <= 0) return 1;
} else {
scale = 1.0;
}
+ if (!do_scale) scale = 1.0;
+
if ( isnan(scale) ) return 1;
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
@@ -319,7 +380,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
int *hist_i, int config_nopolar, int min_measurements,
double min_snr, double max_adu,
int start_after, int stop_after, double min_res,
- double push_res)
+ double push_res, double min_cc, int do_scale, char *stat_output)
{
int rval;
int n_images = 0;
@@ -328,6 +389,14 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
Reflection *refl;
RefListIterator *iter;
int n_crystals_seen = 0;
+ FILE *stat = NULL;
+
+ if (stat_output != NULL) {
+ stat = fopen(stat_output, "w");
+ if (stat == NULL) {
+ ERROR("Failed to open statistics output file %s\n", stat_output);
+ }
+ }
do {
@@ -359,7 +428,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
r = merge_crystal(model, &image, cr, reference, sym,
hist_vals, hist_h, hist_k, hist_l,
hist_i, config_nopolar, min_snr,
- max_adu, push_res);
+ max_adu, push_res, min_cc, do_scale, stat);
if ( r == 0 ) n_crystals_used++;
@@ -398,6 +467,9 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
set_esd_intensity(refl, sqrt(var)/sqrt(red));
}
+ if (stat != NULL)
+ fclose(stat);
+
return 0;
}
@@ -407,6 +479,7 @@ int main(int argc, char *argv[])
int c;
char *filename = NULL;
char *output = NULL;
+ char *stat_output = NULL;
Stream *st;
RefList *model;
int config_scale = 0;
@@ -430,6 +503,8 @@ int main(int argc, char *argv[])
double max_adu = +INFINITY;
double min_res = 0.0;
double push_res = +INFINITY;
+ double min_cc = -INFINITY;
+ int do_scale = 0;
/* Long options */
const struct option longopts[] = {
@@ -438,7 +513,7 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"start-after", 1, NULL, 's'},
{"stop-after", 1, NULL, 'f'},
- {"scale", 0, &config_scale, 1},
+ {"scale", 0, &do_scale, 1},
{"no-polarisation", 0, &config_nopolar, 1},
{"no-polarization", 0, &config_nopolar, 1},
{"symmetry", 1, NULL, 'y'},
@@ -451,6 +526,8 @@ int main(int argc, char *argv[])
{"push-res", 1, NULL, 6},
{"res-push", 1, NULL, 6}, /* compat */
{"version", 0, NULL, 7},
+ {"min_cc", 1, NULL, 8},
+ {"stat", 1, NULL, 9},
{0, 0, NULL, 0}
};
@@ -472,6 +549,10 @@ int main(int argc, char *argv[])
output = strdup(optarg);
break;
+ case 9 :
+ stat_output = strdup(optarg);
+ break;
+
case 's' :
errno = 0;
start_after = strtod(optarg, &rval);
@@ -561,6 +642,16 @@ int main(int argc, char *argv[])
case '?' :
break;
+ case 8 :
+ errno = 0;
+ min_cc = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --min-cc.\n");
+ return 1;
+ }
+ config_scale = 1;
+ break;
+
case 0 :
break;
@@ -572,6 +663,8 @@ int main(int argc, char *argv[])
}
+ if (do_scale) config_scale = 1;
+
if ( filename == NULL ) {
ERROR("Please specify filename using the -i option\n");
return 1;
@@ -639,7 +732,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, min_snr,
- max_adu, start_after, stop_after, min_res, push_res);
+ max_adu, start_after, stop_after, min_res, push_res, min_cc, do_scale, stat_output);
fprintf(stderr, "\n");
if ( r ) {
ERROR("Error while reading stream.\n");
@@ -674,7 +767,7 @@ int main(int argc, char *argv[])
&hist_vals, hist_h, hist_k, hist_l, &hist_i,
config_nopolar, min_measurements, min_snr,
max_adu, start_after, stop_after, min_res,
- push_res);
+ push_res, min_cc, do_scale, stat_output);
fprintf(stderr, "\n");
if ( r ) {
ERROR("Error while reading stream.\n");
@@ -706,6 +799,8 @@ int main(int argc, char *argv[])
reflist_free(model);
free(output);
free(filename);
+ if (stat_output != NULL)
+ free(stat_output);
return 0;
}