Merge branch 'tom/multicrystal'

Conflicts:
	libcrystfel/src/mosflm.c

1 files changed, 116 insertions, 92 deletions

diff --git a/src/process_hkl.c b/src/process_hkl.c
index 8bfbea2b..8705ce0f 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -3,12 +3,12 @@
  *
  * Assemble and process FEL Bragg intensities
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2013 Thomas White <taw@physics.org>
  *   2011      Andrew Martin <andrew.martin@desy.de>
  *   2012      Lorenzo Galli <lorenzo.galli@desy.de>
  *
@@ -48,6 +48,8 @@
 #include "stream.h"
 #include "reflist.h"
 #include "image.h"
+#include "crystal.h"
+#include "thread-pool.h"
 
 
 static void show_help(const char *s)
@@ -62,8 +64,8 @@ static void show_help(const char *s)
 "                             Default: processed.hkl).\n"
 "  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
 "\n"
-"      --start-after=<n>     Skip n patterns at the start of the stream.\n"
-"      --stop-after=<n>      Stop after processing n patterns.\n"
+"      --start-after=<n>     Skip <n> crystals at the start of the stream.\n"
+"      --stop-after=<n>      Stop after merging <n> crystals.\n"
 "  -g, --histogram=<h,k,l>   Calculate the histogram of measurements for this\n"
 "                             reflection.\n"
 "  -z, --hist-parameters     Set the range for the histogram and the number of\n"
@@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new,
 }
 
 
-static int merge_pattern(RefList *model, struct image *new, 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)
+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)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
 	double csx, csy, csz;
 
 	if ( reference != NULL ) {
-		scale = scale_intensities(reference, new->reflections, sym);
+		scale = scale_intensities(reference,
+		                          crystal_get_reflections(cr), sym);
 	} else {
 		scale = 1.0;
 	}
 	if ( isnan(scale) ) return 1;
 
-	cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz,
-	                                       &bsx, &bsy, &bsz,
-	                                       &csx, &csy, &csz);
+	cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
+	                                          &bsx, &bsy, &bsz,
+	                                          &csx, &csy, &csz);
 
-	for ( refl = first_refl(new->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
 			yl = h*asy + k*bsy + l*csy;
 			zl = h*asz + k*bsz + l*csz;
 
-			ool = 1.0 / new->lambda;
+			ool = 1.0 / image->lambda;
 			tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+ool);
 			phi = atan2(yl, xl);
 			pa = pow(sin(phi)*sin(tt), 2.0);
@@ -260,60 +263,76 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference,
 }
 
 
-static void merge_all(FILE *fh, RefList *model, RefList *reference,
-                      int config_startafter, int config_stopafter,
-                      const SymOpList *sym,
-                      int n_total_patterns,
-                      double *hist_vals, signed int hist_h,
-                      signed int hist_k, signed int hist_l,
-                      int *hist_i, int config_nopolar, int min_measurements)
+static void display_progress(int n_images, int n_crystals, int n_crystals_used)
+{
+	if ( !isatty(STDERR_FILENO) ) return;
+	if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+	pthread_mutex_lock(&stderr_lock);
+	fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used.",
+	        n_images, n_crystals, n_crystals_used);
+	pthread_mutex_unlock(&stderr_lock);
+
+	fflush(stdout);
+}
+
+
+
+static int merge_all(Stream *st, RefList *model, RefList *reference,
+                     const SymOpList *sym,
+                     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)
 {
 	int rval;
-	int n_patterns = 0;
-	int n_used = 0;
+	int n_images = 0;
+	int n_crystals = 0;
+	int n_crystals_used = 0;
 	Reflection *refl;
 	RefListIterator *iter;
-
-	if ( skip_some_files(fh, config_startafter) ) {
-		ERROR("Failed to skip first %i files.\n", config_startafter);
-		return;
-	}
+	int n_crystals_seen = 0;
 
 	do {
 
 		struct image image;
+		int i;
 
 		image.det = NULL;
 
 		/* Get data from next chunk */
-		rval = read_chunk(fh, &image);
+		rval = read_chunk(st, &image);
 		if ( rval ) break;
 
-		n_patterns++;
+		n_images++;
 
-		if ( (image.reflections != NULL) && (image.indexed_cell) ) {
+		for ( i=0; i<image.n_crystals; i++ ) {
 
 			int r;
+			Crystal *cr = image.crystals[i];
 
-			r = merge_pattern(model, &image, reference, sym,
+			n_crystals_seen++;
+			if ( n_crystals_seen <= start_after ) continue;
+
+			n_crystals++;
+			r = merge_crystal(model, &image, cr, reference, sym,
 			                  hist_vals, hist_h, hist_k, hist_l,
 			                  hist_i, config_nopolar);
 
-			if ( r == 0 ) n_used++;
+			if ( r == 0 ) n_crystals_used++;
+
+			crystal_free(cr);
+
+			if ( n_crystals_used == stop_after ) break;
 
 		}
 
 		free(image.filename);
-		reflist_free(image.reflections);
 		image_feature_list_free(image.features);
-		cell_free(image.indexed_cell);
 
-		progress_bar(n_patterns, n_total_patterns-config_startafter,
-		             "Merging");
+		display_progress(n_images, n_crystals_seen, n_crystals_used);
 
-		if ( config_stopafter ) {
-			if ( n_patterns == config_stopafter ) break;
-		}
+		if ( (stop_after>0) && (n_crystals_used == stop_after) ) break;
 
 	} while ( rval == 0 );
 
@@ -339,7 +358,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference,
 
 	}
 
-	STATUS("%i of the patterns could be used.\n", n_used);
+	return 0;
 }
 
 
@@ -348,14 +367,11 @@ int main(int argc, char *argv[])
 	int c;
 	char *filename = NULL;
 	char *output = NULL;
-	FILE *fh;
+	Stream *st;
 	RefList *model;
 	int config_maxonly = 0;
-	int config_startafter = 0;
-	int config_stopafter = 0;
 	int config_sum = 0;
 	int config_scale = 0;
-	unsigned int n_total_patterns;
 	char *sym_str = NULL;
 	SymOpList *sym;
 	char *histo = NULL;
@@ -369,6 +385,9 @@ int main(int argc, char *argv[])
 	int config_nopolar = 0;
 	char *rval;
 	int min_measurements = 2;
+	int r;
+	int start_after = 0;
+	int stop_after = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -377,8 +396,8 @@ int main(int argc, char *argv[])
 		{"output",             1, NULL,               'o'},
 		{"max-only",           0, &config_maxonly,     1},
 		{"output-every",       1, NULL,               'e'},
-		{"stop-after",         1, NULL,               's'},
-		{"start-after",        1, NULL,               'f'},
+		{"start-after",         1, NULL,              's'},
+		{"stop-after",        1, NULL,                'f'},
 		{"sum",                0, &config_sum,         1},
 		{"scale",              0, &config_scale,       1},
 		{"no-polarisation",    0, &config_nopolar,     1},
@@ -391,7 +410,7 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:e:o:y:g:f:b:z:",
+	while ((c = getopt_long(argc, argv, "hi:e:o:y:g:s:f:z:",
 	                        longopts, NULL)) != -1) {
 
 		switch (c) {
@@ -409,11 +428,21 @@ int main(int argc, char *argv[])
 			break;
 
 			case 's' :
-			config_stopafter = atoi(optarg);
+			errno = 0;
+			start_after = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --start-after.\n");
+				return 1;
+			}
 			break;
 
 			case 'f' :
-			config_startafter = atoi(optarg);
+			errno = 0;
+			stop_after = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --stop-after.\n");
+				return 1;
+			}
 			break;
 
 			case 'y' :
@@ -465,25 +494,11 @@ int main(int argc, char *argv[])
 	free(sym_str);
 
 	/* Open the data stream */
-	if ( strcmp(filename, "-") == 0 ) {
-		fh = stdin;
-	} else {
-		fh = fopen(filename, "r");
-	}
-	free(filename);
-	if ( fh == NULL ) {
-		ERROR("Failed to open input file\n");
-		return 1;
-	}
-
-	/* Count the number of patterns in the file */
-	n_total_patterns = count_patterns(fh);
-	if ( n_total_patterns == 0 ) {
-		ERROR("No patterns to process.\n");
+	st = open_stream_for_read(filename);
+	if ( st == NULL ) {
+		ERROR("Failed to open stream.\n");
 		return 1;
 	}
-	STATUS("There are %i patterns to process\n", n_total_patterns);
-	rewind(fh);
 
 	model = reflist_new();
 
@@ -497,7 +512,8 @@ int main(int argc, char *argv[])
 			return 1;
 		}
 
-		space_for_hist = n_total_patterns * num_equivs(sym, NULL);
+		/* FIXME: This array must grow as necessary. */
+		space_for_hist = 0 * num_equivs(sym, NULL);
 		hist_vals = malloc(space_for_hist * sizeof(double));
 		free(histo);
 		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
@@ -529,40 +545,48 @@ int main(int argc, char *argv[])
 	}
 
 	hist_i = 0;
-	merge_all(fh, model, NULL, config_startafter, config_stopafter,
-	          sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l,
-	          &hist_i, config_nopolar, min_measurements);
-	if ( ferror(fh) ) {
-		ERROR("Stream read error.\n");
+	r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l,
+	              &hist_i, config_nopolar, min_measurements,
+	              start_after, stop_after);
+	fprintf(stderr, "\n");
+	if ( r ) {
+		ERROR("Error while reading stream.\n");
 		return 1;
 	}
-	rewind(fh);
 
 	if ( config_scale ) {
 
 		RefList *reference;
 
-		STATUS("Extra pass for scaling...\n");
+		if ( rewind_stream(st) ) {
 
-		reference = copy_reflist(model);
+			ERROR("Couldn't rewind stream - scaling cannot be "
+			      "performed.\n");
 
-		reflist_free(model);
-		model = reflist_new();
+		} else {
 
-		rewind(fh);
+			int r;
 
-		merge_all(fh, model, reference,
-			  config_startafter, config_stopafter, sym,
-			  n_total_patterns,
-			  hist_vals, hist_h, hist_k, hist_l, &hist_i,
-			  config_nopolar, min_measurements);
+			STATUS("Extra pass for scaling...\n");
 
-		if ( ferror(fh) ) {
-			ERROR("Stream read error.\n");
-			return 1;
-		}
+			reference = copy_reflist(model);
+
+			reflist_free(model);
+			model = reflist_new();
+
+			r = merge_all(st, model, reference, sym,
+				     hist_vals, hist_h, hist_k, hist_l, &hist_i,
+				     config_nopolar, min_measurements,
+				     start_after, stop_after);
+			fprintf(stderr, "\n");
+			if ( r ) {
+				ERROR("Error while reading stream.\n");
+				return 1;
+			}
 
-		reflist_free(reference);
+			reflist_free(reference);
+
+		}
 
 	}
 
@@ -579,7 +603,7 @@ int main(int argc, char *argv[])
 
 	write_reflist(output, model);
 
-	fclose(fh);
+	close_stream(st);
 
 	free(sym);
 	reflist_free(model);