Indexing pipeline - "done"!

3 files changed, 121 insertions, 163 deletions

diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index c8cceb2b..825f2d6a 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -73,14 +73,27 @@
 #define STATS_EVERY_N_SECONDS (5)
 
 
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+	/* "Input" */
+	char *filename;
+
+	/* "Output" */
+	int n_crystals;
+};
+
+
 struct sandbox
 {
 	pthread_mutex_t lock;
 
-	int n_indexable;
 	int n_processed;
-	int n_indexable_last_stats;
+	int n_hadcrystals;
+	int n_crystals;
 	int n_processed_last_stats;
+	int n_hadcrystals_last_stats;
+	int n_crystals_last_stats;
 	int t_last_stats;
 
 	struct index_args *iargs;
@@ -164,31 +177,29 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
 
 
 static void process_image(const struct index_args *iargs,
-                          struct pattern_args *pargs, FILE *ofh,
+                          struct pattern_args *pargs, Stream *st,
                           int cookie)
 {
 	float *data_for_measurement;
 	size_t data_size;
-	UnitCell *cell = iargs->cell;
 	int config_cmfilter = iargs->config_cmfilter;
 	int config_noisefilter = iargs->config_noisefilter;
-	int config_verbose = iargs->config_verbose;
 	IndexingMethod *indm = iargs->indm;
 	int check;
 	struct hdfile *hdfile;
 	struct image image;
+	int i;
 
 	image.features = NULL;
 	image.data = NULL;
 	image.flags = NULL;
-	image.indexed_cell = NULL;
 	image.copyme = iargs->copyme;
-	image.reflections = NULL;
-	image.n_saturated = 0;
 	image.id = cookie;
 	image.filename = pargs->filename;
 	image.beam = iargs->beam;
 	image.det = iargs->det;
+	image.crystals = NULL;
+	image.n_crystals = 0;
 
 	hdfile = hdfile_open(image.filename);
 	if ( hdfile == NULL ) return;
@@ -291,28 +302,34 @@ static void process_image(const struct index_args *iargs,
 	free(image.data);
 	image.data = data_for_measurement;
 
-	/* Calculate orientation matrix (by magic) */
+	/* Index the pattern */
+	index_pattern(&image, indm, iargs->ipriv);
+
+	/* Default beam parameters */
 	image.div = image.beam->divergence;
 	image.bw = image.beam->bandwidth;
-	image.profile_radius = image.beam->profile_radius;
 
-	index_pattern(&image, cell, indm, iargs->cellr,
-	              config_verbose, iargs->ipriv,
-	              iargs->config_insane, iargs->tols);
+	/* Integrate each crystal's diffraction spots */
+	for ( i=0; i<image.n_crystals; i++ ) {
 
-	if ( image.indexed_cell != NULL ) {
+		RefList *reflections;
 
-		pargs->indexable = 1;
+		/* Set default crystal parameter(s) */
+		crystal_set_profile_radius(image.crystals[i],
+		                           image.beam->profile_radius);
 
 		if ( iargs->integrate_found ) {
-			image.reflections = select_intersections(&image,
-			                                    image.indexed_cell);
+			reflections = select_intersections(&image,
+			                                   image.crystals[i]);
 		} else {
-			image.reflections = find_intersections(&image,
-			                                    image.indexed_cell);
+			reflections = find_intersections(&image,
+			                                 image.crystals[i]);
 		}
 
-		if (image.reflections != NULL) {
+		crystal_set_reflections(image.crystals[i], reflections);
+
+		if ( reflections != NULL ) {
+
 			integrate_reflections(&image,
 			                      iargs->config_closer,
 			                      iargs->config_bgsub,
@@ -323,18 +340,15 @@ static void process_image(const struct index_args *iargs,
 			                      iargs->integrate_saturated);
 		}
 
-	} else {
-		image.reflections = NULL;
 	}
 
-	write_chunk(ofh, &image, hdfile, iargs->stream_flags);
-	fprintf(ofh, "END\n");
-	fflush(ofh);
+	write_chunk(st, &image, hdfile, iargs->include_peaks,
+	            iargs->include_reflections);
 
-	/* Only free cell if found */
-	cell_free(image.indexed_cell);
+	for ( i=0; i<image.n_crystals; i++ ) {
+		crystal_free(image.crystals[i]);
+	}
 
-	reflist_free(image.reflections);
 	free(image.data);
 	if ( image.flags != NULL ) free(image.flags);
 	image_feature_list_free(image.features);
@@ -343,7 +357,7 @@ static void process_image(const struct index_args *iargs,
 
 
 static void run_work(const struct index_args *iargs,
-                     int filename_pipe, int results_pipe, FILE *ofh,
+                     int filename_pipe, int results_pipe, Stream *st,
                      int cookie)
 {
 	int allDone = 0;
@@ -389,12 +403,12 @@ static void run_work(const struct index_args *iargs,
 		} else {
 
 			pargs.filename = line;
-			pargs.indexable = 0;
+			pargs.n_crystals = 0;
 
-			process_image(iargs, &pargs, ofh, cookie);
+			process_image(iargs, &pargs, st, cookie);
 
 			/* Request another image */
-			c = sprintf(buf, "%i\n", pargs.indexable);
+			c = sprintf(buf, "%i\n", pargs.n_crystals);
 			w = write(results_pipe, buf, c);
 			if ( w < 0 ) {
 				ERROR("write P0\n");
@@ -406,7 +420,7 @@ static void run_work(const struct index_args *iargs,
 
 	}
 
-	cleanup_indexing(iargs->ipriv);
+	cleanup_indexing(iargs->indm, iargs->ipriv);
 	free(iargs->indm);
 	free(iargs->ipriv);
 	free_detector_geometry(iargs->det);
@@ -577,7 +591,7 @@ static void start_worker_process(struct sandbox *sb, int slot)
 
 	if ( p == 0 ) {
 
-		FILE *sfh;
+		Stream *st;
 		int j;
 		struct sigaction sa;
 		int r;
@@ -614,10 +628,10 @@ static void start_worker_process(struct sandbox *sb, int slot)
 		close(filename_pipe[1]);
 		close(result_pipe[0]);
 
-		sfh = fdopen(sb->stream_pipe_write[slot], "w");
+		st = open_stream_fd_for_write(sb->stream_pipe_write[slot]);
 		run_work(sb->iargs, filename_pipe[0], result_pipe[1],
-		         sfh, slot);
-		fclose(sfh);
+		         st, slot);
+		close_stream(st);
 
 		//close(filename_pipe[0]);
 		close(result_pipe[1]);
@@ -697,7 +711,7 @@ static void handle_zombie(struct sandbox *sb)
 
 void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
                     int config_basename, FILE *fh, char *use_this_one_instead,
-                    FILE *ofh)
+                    Stream *st)
 {
 	int i;
 	int allDone;
@@ -712,13 +726,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		return;
 	}
 
-	sb->n_indexable = 0;
 	sb->n_processed = 0;
-	sb->n_indexable_last_stats = 0;
+	sb->n_hadcrystals = 0;
+	sb->n_crystals = 0;
 	sb->n_processed_last_stats = 0;
+	sb->n_hadcrystals_last_stats = 0;
+	sb->n_crystals_last_stats = 0;
 	sb->t_last_stats = get_monotonic_seconds();
 	sb->n_proc = n_proc;
-	sb->ofh = ofh;
 	sb->iargs = iargs;
 
 	pthread_mutex_init(&sb->lock, NULL);
@@ -895,8 +910,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 					ERROR("Invalid result '%s'\n", results);
 				}
 			} else {
-				sb->n_indexable += atoi(results);
+
+				int nc = atoi(results);
+				sb->n_crystals += nc;
+				if ( nc > 0 ) {
+					sb->n_hadcrystals++;
+				}
 				sb->n_processed++;
+
 			}
 
 			/* Send next filename */
@@ -929,14 +950,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		tNow = get_monotonic_seconds();
 		if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) {
 
-			STATUS("%i out of %i indexed so far,"
-			       " %i out of %i since the last message.\n",
-			       sb->n_indexable, sb->n_processed,
-			       sb->n_indexable - sb->n_indexable_last_stats,
-			       sb->n_processed - sb->n_processed_last_stats);
+			STATUS("Total so far: %i images processed, "
+			       "%i had crystals, %i crystals overall.  "
+			       "Since the last message: %i images processed,"
+			       "%i had crystals, %i crystals overall.\n",
+			       sb->n_processed, sb->n_hadcrystals,
+			       sb->n_crystals,
+			       sb->n_processed - sb->n_processed_last_stats,
+			       sb->n_hadcrystals - sb->n_hadcrystals_last_stats,
+			       sb->n_crystals - sb->n_crystals_last_stats);
 
-			sb->n_indexable_last_stats = sb->n_indexable;
 			sb->n_processed_last_stats = sb->n_processed;
+			sb->n_hadcrystals_last_stats = sb->n_hadcrystals;
+			sb->n_crystals_last_stats = sb->n_crystals;
 			sb->t_last_stats = tNow;
 
 		}
@@ -975,9 +1001,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 	free(sb->pids);
 	free(sb->running);
 
-	if ( ofh != stdout ) fclose(ofh);
-
-	STATUS("There were %i images, of which %i could be indexed.\n",
-	       sb->n_processed, sb->n_indexable);
+	STATUS("Final:"
+	       " %i images processed, %i had crystals, %i crystals overall.\n",
+	       sb->n_processed, sb->n_hadcrystals, sb->n_crystals);
 
 }
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 872f84ad..176bfde0 100644
--- a/src/im-sandbox.h
+++ b/src/im-sandbox.h
@@ -3,13 +3,13 @@
  *
  * Sandbox for indexing
  *
- * 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 Richard Kirian
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *   2011      Richard Kirian
  *   2012      Lorenzo Galli
  *   2012      Chunhong Yoon
@@ -31,6 +31,8 @@
  *
  */
 
+#include "stream.h"
+
 enum {
 	PEAK_ZAEF,
 	PEAK_HDF5,
@@ -44,10 +46,8 @@ struct index_args
 	int config_cmfilter;
 	int config_noisefilter;
 	int config_verbose;
-	int stream_flags;         /* What goes into the output? */
 	int config_satcorr;
 	int config_closer;
-	int config_insane;
 	int config_bgsub;
 	float threshold;
 	float min_gradient;
@@ -57,7 +57,6 @@ struct index_args
 	IndexingMethod *indm;
 	IndexingPrivate **ipriv;
 	int peaks;                /* Peak detection method */
-	int cellr;
 	float tols[4];
 	struct beam_params *beam;
 	char *element;
@@ -70,21 +69,11 @@ struct index_args
 	int use_saturated;
 	int no_revalidate;
 	int integrate_found;
+	int include_peaks;
+	int include_reflections;
 };
 
 
-
-
-/* Information about the indexing process for one pattern */
-struct pattern_args
-{
-	/* "Input" */
-	char *filename;
-
-	/* "Output" */
-	int indexable;
-};
-
 extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
                            int config_basename, FILE *fh,
-                           char *use_this_one_instead, FILE *ofh);
+                           char *use_this_one_instead, Stream *st);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 01377de3..3076ea86 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -174,37 +174,13 @@ static void show_help(const char *s)
 }
 
 
-static int parse_cell_reduction(const char *scellr, int *err,
-                                int *reduction_needs_cell)
-{
-	*err = 0;
-	if ( strcmp(scellr, "none") == 0 ) {
-		*reduction_needs_cell = 0;
-		return CELLR_NONE;
-	} else if ( strcmp(scellr, "reduce") == 0) {
-		*reduction_needs_cell = 1;
-		return CELLR_REDUCE;
-	} else if ( strcmp(scellr, "compare") == 0) {
-		*reduction_needs_cell = 1;
-		return CELLR_COMPARE;
-	} else if ( strcmp(scellr, "compare_ab") == 0) {
-		*reduction_needs_cell = 1;
-		return CELLR_COMPARE_AB;
-	} else {
-		*err = 1;
-		*reduction_needs_cell = 0;
-		return CELLR_NONE;
-	}
-}
-
-
 int main(int argc, char *argv[])
 {
 	int c;
 	char *filename = NULL;
 	char *outfile = NULL;
 	FILE *fh;
-	FILE *ofh;
+	Stream *st;
 	char *rval = NULL;
 	int config_noindex = 0;
 	int config_cmfilter = 0;
@@ -213,7 +189,6 @@ int main(int argc, char *argv[])
 	int config_satcorr = 1;
 	int config_checkprefix = 1;
 	int config_closer = 0;
-	int config_insane = 0;
 	int config_bgsub = 1;
 	int config_basename = 0;
 	float threshold = 800.0;
@@ -224,17 +199,13 @@ int main(int argc, char *argv[])
 	char *geometry = NULL;
 	IndexingMethod *indm;
 	IndexingPrivate **ipriv;
-	int indexer_needs_cell;
-	int reduction_needs_cell;
 	char *indm_str = NULL;
 	UnitCell *cell;
 	char *pdb = NULL;
 	char *prefix = NULL;
 	char *speaks = NULL;
-	char *scellr = NULL;
 	char *toler = NULL;
 	float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
-	int cellr;
 	int peaks;
 	int n_proc = 1;
 	char *prepare_line;
@@ -243,7 +214,6 @@ int main(int argc, char *argv[])
 	struct index_args iargs;
 	struct beam_params *beam = NULL;
 	char *element = NULL;
-	int stream_flags = STREAM_INTEGRATED;
 	char *hdf5_peak_path = NULL;
 	struct copy_hdf5_field *copyme;
 	char *intrad = NULL;
@@ -255,6 +225,10 @@ int main(int argc, char *argv[])
 	int no_revalidate = 0;
 	int integrate_found = 0;
 
+	/* For ease of upgrading from old method */
+	char *scellr = NULL;
+	int config_insane = 0;
+
 	copyme = new_copy_hdf5_field_list();
 	if ( copyme == NULL ) {
 		ERROR("Couldn't allocate HDF5 field list.\n");
@@ -290,7 +264,6 @@ int main(int argc, char *argv[])
 		{"peaks",              1, NULL,                2},
 		{"cell-reduction",     1, NULL,                3},
 		{"min-gradient",       1, NULL,                4},
-		{"record",             1, NULL,                5},
 		{"cpus",               1, NULL,                6},
 		{"cpugroup",           1, NULL,                7},
 		{"cpuoffset",          1, NULL,                8},
@@ -301,6 +274,8 @@ int main(int argc, char *argv[])
 		{"tolerance",          1, NULL,               13},
 		{"int-radius",         1, NULL,               14},
 
+		/* FIXME: Add '--no-peaks' and '--no-reflections' */
+
 		{"integrate-saturated",0, &integrate_saturated,1},
 		{"use-saturated",      0, &use_saturated,      1},
 		{"no-revalidate",      0, &no_revalidate,      1},
@@ -376,8 +351,7 @@ int main(int argc, char *argv[])
 			break;
 
 			case 5 :
-			stream_flags = parse_stream_flags(optarg);
-			if ( stream_flags < 0 ) return 1;
+			ERROR("The option '--record' is no longer used.\n");
 			break;
 
 			case 6 :
@@ -425,6 +399,12 @@ int main(int argc, char *argv[])
 
 	}
 
+	if ( scellr != NULL ) {
+		ERROR("Old-style indexing method specification.\n");
+		/* FIXME: Translate to new style */
+		return 1;
+	}
+
 	if ( filename == NULL ) {
 		filename = strdup("-");
 	}
@@ -439,17 +419,6 @@ int main(int argc, char *argv[])
 	}
 	free(filename);
 
-	if ( outfile == NULL ) {
-		ofh = stdout;
-	} else {
-		ofh = fopen(outfile, "w");
-		if ( ofh == NULL ) {
-			ERROR("Failed to open output file '%s'\n", outfile);
-			return 1;
-		}
-		free(outfile);
-	}
-
 	if ( hdf5_peak_path == NULL ) {
 		hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
 	}
@@ -482,53 +451,24 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	if ( (indm_str == NULL) ||
-	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
-		STATUS("Not indexing anything.\n");
-		indexer_needs_cell = 0;
-		reduction_needs_cell = 0;
+	if ( indm_str == NULL ) {
+
+		STATUS("You didn't specify an indexing method, so I  won't try "
+		       " to index anything.\n"
+		       "If that isn't what you wanted, re-run with"
+		       " --indexing=<methods>.\n");
 		indm = NULL;
-		cellr = CELLR_NONE;
+
 	} else {
-		if ( indm_str == NULL ) {
-			STATUS("You didn't specify an indexing method, so I "
-			       " won't try to index anything.\n"
-			       "If that isn't what you wanted, re-run with"
-			       " --indexing=<method>.\n");
-			indm = NULL;
-			indexer_needs_cell = 0;
-		} else {
-			indm = build_indexer_list(indm_str,
-			                          &indexer_needs_cell);
-			if ( indm == NULL ) {
-				ERROR("Invalid indexer list '%s'\n", indm_str);
-				return 1;
-			}
-			free(indm_str);
-		}
 
-		reduction_needs_cell = 0;
-		if ( scellr == NULL ) {
-			STATUS("You didn't specify a cell reduction method, so"
-			       " I'm going to use 'reduce'.\n");
-			cellr = CELLR_REDUCE;
-			reduction_needs_cell = 1;
-		} else {
-			int err;
-			cellr = parse_cell_reduction(scellr, &err,
-			                             &reduction_needs_cell);
-			if ( err ) {
-				ERROR("Unrecognised cell reduction '%s'\n",
-			              scellr);
-				return 1;
-			}
-			free(scellr);
+		indm = build_indexer_list(indm_str);
+		if ( indm == NULL ) {
+			ERROR("Invalid indexer list '%s'\n", indm_str);
+			return 1;
 		}
+		free(indm_str);
 	}
 
-	/* No indexing -> no reduction */
-	if ( indm == NULL ) reduction_needs_cell = 0;
-
 	if ( toler != NULL ) {
 		int ttt;
 		ttt = sscanf(toler, "%f,%f,%f,%f",
@@ -586,7 +526,12 @@ int main(int argc, char *argv[])
 		cell = NULL;
 	}
 
-	write_stream_header(ofh, argc, argv);
+	st = open_stream_for_write(outfile);
+	if ( st == NULL ) {
+		ERROR("Failed to open stream '%s'\n", outfile);
+		return 1;
+	}
+	free(outfile);
 
 	/* Get first filename and use it to set up the indexing */
 	prepare_line = malloc(1024);
@@ -609,7 +554,7 @@ int main(int argc, char *argv[])
 	/* Prepare the indexer */
 	if ( indm != NULL ) {
 		ipriv = prepare_indexing(indm, cell, prepare_filename,
-		                         det, beam);
+		                         det, beam, tols);
 		if ( ipriv == NULL ) {
 			ERROR("Failed to prepare indexing.\n");
 			return 1;
@@ -627,9 +572,7 @@ int main(int argc, char *argv[])
 	iargs.config_verbose = config_verbose;
 	iargs.config_satcorr = config_satcorr;
 	iargs.config_closer = config_closer;
-	iargs.config_insane = config_insane;
 	iargs.config_bgsub = config_bgsub;
-	iargs.cellr = cellr;
 	iargs.tols[0] = tols[0];
 	iargs.tols[1] = tols[1];
 	iargs.tols[2] = tols[2];
@@ -644,7 +587,6 @@ int main(int argc, char *argv[])
 	iargs.peaks = peaks;
 	iargs.beam = beam;
 	iargs.element = element;
-	iargs.stream_flags = stream_flags;
 	iargs.hdf5_peak_path = hdf5_peak_path;
 	iargs.copyme = copyme;
 	iargs.ir_inn = ir_inn;
@@ -654,9 +596,11 @@ int main(int argc, char *argv[])
 	iargs.integrate_saturated = integrate_saturated;
 	iargs.no_revalidate = no_revalidate;
 	iargs.integrate_found = integrate_found;
+	iargs.include_peaks = 1;
+	iargs.include_reflections = 1;  /* FIXME! */
 
 	create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
-	                    use_this_one_instead, ofh);
+	                    use_this_one_instead, st);
 
 	free(prefix);