Calculate CC for each frame. Reject bad frames (--min-cc). Output stats (--stat=)

But these should be in partialator?

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;
 }