Merge branch 'tom/multicrystal'

Conflicts:
	libcrystfel/src/mosflm.c

16 files changed, 1173 insertions, 958 deletions

diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 6ef12e7a..c4eabaf2 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -1819,7 +1819,6 @@ DisplayWindow *displaywindow_open(const char *filename, const char *peaks,
 
 		dw->hdfile = hdfile_open(filename);
 		if ( dw->hdfile == NULL ) {
-			ERROR("Couldn't open file '%s'\n", filename);
 			free(dw);
 			return NULL;
 		} else {
diff --git a/src/hdfsee.c b/src/hdfsee.c
index 620b684b..36472e46 100644
--- a/src/hdfsee.c
+++ b/src/hdfsee.c
@@ -266,11 +266,7 @@ int main(int argc, char *argv[])
 		                                         ring_radii,
 		                                         n_rings,
 		                                         ring_size);
-		if ( main_window_list[i] == NULL ) {
-			ERROR("Couldn't open display window\n");
-		} else {
-			main_n_windows++;
-		}
+		if ( main_window_list[i] != NULL ) main_n_windows++;
 	}
 
 	if ( main_n_windows == 0 ) return 0;
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index f86ae0bb..0498a532 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -57,7 +57,7 @@
 struct scale_queue_args
 {
 	RefList *reference;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 	double max_shift;
 };
@@ -65,7 +65,7 @@ struct scale_queue_args
 
 struct scale_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	double shift;
 	RefList *reference;
 };
@@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs)
 	wargs = malloc(sizeof(struct scale_worker_args));
 	wargs->reference = qargs->reference;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs)
 static void run_scale_job(void *vwargs, int cookie)
 {
 	struct scale_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *reference = wargs->reference;
 	Reflection *refl;
 	RefListIterator *iter;
 	double num = 0.0;
 	double den = 0.0;
 	double corr;
+	const double osf = crystal_get_osf(cr);
 
-	if ( image->pr_dud ) {
+	if ( crystal_get_user_flag(cr) ) {
 		wargs->shift = 0.0;
 		return;
 	}
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie)
 		Ihl = get_intensity(refl) / get_partiality(refl);
 		esd = get_esd_intensity(refl) / get_partiality(refl);
 
-		num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0);
-		den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0);
+		num += Ih * (Ihl/osf) / pow(esd/osf, 2.0);
+		den += pow(Ih, 2.0)/pow(esd/osf, 2.0);
 
 	}
 
@@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie)
 
 	corr = num / den;
 	if ( !isnan(corr) && !isinf(corr) ) {
-		image->osf *= corr;
+		crystal_set_osf(cr, osf*corr);
 	}
 	wargs->shift = fabs(corr-1.0);
 
@@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs)
 }
 
 
-static double iterate_scale(struct image *images, int n, RefList *reference,
+static double iterate_scale(Crystal **crystals, int n, RefList *reference,
                             int n_threads)
 {
 	struct scale_queue_args qargs;
@@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference,
 
 	qargs.reference = reference;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 	qargs.max_shift = 0.0;
 
 	run_threads(n_threads, run_scale_job, create_scale_job,
@@ -178,14 +179,14 @@ struct merge_queue_args
 {
 	RefList *full;
 	pthread_mutex_t full_lock;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 };
 
 
 struct merge_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	RefList *full;
 	pthread_mutex_t *full_lock;
 };
@@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs)
 	wargs->full = qargs->full;
 	wargs->full_lock = &qargs->full_lock;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs)
 static void run_merge_job(void *vwargs, int cookie)
 {
 	struct merge_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *full = wargs->full;
 	Reflection *refl;
 	RefListIterator *iter;
 	double G;
 
-	if ( image->pr_dud ) return;
+	if ( crystal_get_user_flag(cr)) return;
 
-	G = image->osf;
+	G = crystal_get_osf(cr);
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 }
 
 
-static RefList *lsq_intensities(struct image *images, int n, int n_threads)
+static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads)
 {
 	RefList *full;
 	struct merge_queue_args qargs;
@@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
 
 	qargs.full = full;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 	pthread_mutex_init(&qargs.full_lock, NULL);
 
 	run_threads(n_threads, run_merge_job, create_merge_job,
@@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads)
 struct esd_queue_args
 {
 	RefList *full;
-	struct image *images;
+	Crystal **crystals;
 	int n_started;
 };
 
 
 struct esd_worker_args
 {
-	struct image *image;
+	Crystal *crystal;
 	RefList *full;
 };
 
@@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs)
 	wargs = malloc(sizeof(struct esd_worker_args));
 	wargs->full = qargs->full;
 
-	wargs->image = &qargs->images[qargs->n_started++];
+	wargs->crystal = qargs->crystals[qargs->n_started++];
 
 	return wargs;
 }
@@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs)
 static void run_esd_job(void *vwargs, int cookie)
 {
 	struct esd_worker_args *wargs = vwargs;
-	struct image *image = wargs->image;
+	Crystal *cr = wargs->crystal;
 	RefList *full = wargs->full;
 	Reflection *refl;
 	RefListIterator *iter;
 	double G;
 
-	if ( image->pr_dud ) return;
+	if ( crystal_get_user_flag(cr) ) return;
 
-	G = image->osf;
+	G = crystal_get_osf(cr);
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs)
 }
 
 
-static void calculate_esds(struct image *images, int n, RefList *full,
+static void calculate_esds(Crystal **crystals, int n, RefList *full,
                            int n_threads, int min_red)
 {
 	struct esd_queue_args qargs;
@@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
 
 	qargs.full = full;
 	qargs.n_started = 0;
-	qargs.images = images;
+	qargs.crystals = crystals;
 
 	for ( refl = first_refl(full, &iter);
 	      refl != NULL;
@@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full,
 
 
 /* Scale the stack of images */
-RefList *scale_intensities(struct image *images, int n, RefList *gref,
+RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
                            int n_threads, int noscale)
 {
 	int i;
@@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 	RefList *full = NULL;
 	const int min_redundancy = 3;
 
-	for ( i=0; i<n; i++ ) images[i].osf = 1.0;
+	for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0);
 
 	if ( noscale ) {
-		full = lsq_intensities(images, n, n_threads);
-		calculate_esds(images, n, full, n_threads, min_redundancy);
+		full = lsq_intensities(crystals, n, n_threads);
+		calculate_esds(crystals, n, full, n_threads, min_redundancy);
 		return full;
 	}
 
 	/* No reference -> create an initial list to refine against */
 	if ( gref == NULL ) {
-		full = lsq_intensities(images, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads);
 	}
 
 	/* Iterate */
@@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 			reference = full;
 		}
 
-		max_corr = iterate_scale(images, n, reference, n_threads);
+		max_corr = iterate_scale(crystals, n, reference, n_threads);
 		//STATUS("Scaling iteration %2i: max correction = %5.2f\n",
 		//       i+1, max_corr);
 
 		/* No reference -> generate list for next iteration */
 		if ( gref == NULL ) {
 			reflist_free(full);
-			full = lsq_intensities(images, n, n_threads);
+			full = lsq_intensities(crystals, n, n_threads);
 		}
 
 		//show_scale_factors(images, n);
@@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref,
 	} while ( (max_corr > 0.01) && (i < MAX_CYCLES) );
 
 	if ( gref != NULL ) {
-		full = lsq_intensities(images, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads);
 	} /* else we already did it */
 
-	calculate_esds(images, n, full, n_threads, min_redundancy);
+	calculate_esds(crystals, n, full, n_threads, min_redundancy);
 
 	return full;
 }
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 5425b1ff..80940347 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -35,9 +35,10 @@
 #endif
 
 
-#include "image.h"
+#include "crystal.h"
+#include "reflist.h"
 
-extern RefList *scale_intensities(struct image *images, int n,
+extern RefList *scale_intensities(Crystal **crystals, int n,
                                   RefList *reference, int n_threads,
                                   int noscale);
 
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index c8cceb2b..03c553b7 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -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
@@ -41,12 +41,11 @@
 #include <string.h>
 #include <unistd.h>
 #include <getopt.h>
-#include <hdf5.h>
-#include <gsl/gsl_errno.h>
 #include <pthread.h>
 #include <sys/wait.h>
 #include <fcntl.h>
 #include <signal.h>
+#include <sys/stat.h>
 
 #ifdef HAVE_CLOCK_GETTIME
 #include <time.h>
@@ -54,48 +53,56 @@
 #include <sys/time.h>
 #endif
 
-#include "utils.h"
-#include "hdf5-file.h"
-#include "index.h"
-#include "peaks.h"
-#include "detector.h"
-#include "filters.h"
-#include "thread-pool.h"
-#include "beam-parameters.h"
-#include "geometry.h"
-#include "stream.h"
-#include "reflist-utils.h"
-
 #include "im-sandbox.h"
+#include "process_image.h"
 
 
 /* Write statistics at APPROXIMATELY this interval */
 #define STATS_EVERY_N_SECONDS (5)
 
 
+struct sb_reader
+{
+	pthread_mutex_t lock;
+	int done;
+
+	/* If a worker process dies unexpectedly (e.g. if it segfaults), then
+	 * the pipe for its output can still stay open for a little while while
+	 * its buffer empties.  The number of pipes being read from is therefore
+	 * not necessarily the same as the number of worker processes. */
+	int n_read;
+	FILE **fhs;
+	int *fds;
+
+	/* Final output fd */
+	int ofd;
+};
+
+
 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;
 
 	int n_proc;
 	pid_t *pids;
-	FILE *ofh;
-	FILE **fhs;
 
 	int *running;
 	FILE **result_fhs;
 	int *filename_pipes;
-	int *stream_pipe_read;
 	int *stream_pipe_write;
 	char **last_filename;
+
+	struct sb_reader *reader;
 };
 
 
@@ -120,6 +127,7 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
 {
 	char *line;
 	char *filename;
+	size_t len;
 
 	do {
 
@@ -153,7 +161,13 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
 		line = tmp;
 	}
 
-	filename = malloc(strlen(prefix)+strlen(line)+1);
+	len = strlen(prefix)+strlen(line)+1;
+
+	/* Round the length of the buffer, too keep Valgrind quiet when it gets
+	 * given to write() a bit later on */
+	len += 4 - (len % 4);
+
+	filename = malloc(len);
 
 	snprintf(filename, 1023, "%s%s", prefix, line);
 
@@ -163,187 +177,8 @@ 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,
-                          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;
-
-	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;
-
-	hdfile = hdfile_open(image.filename);
-	if ( hdfile == NULL ) return;
-
-	if ( iargs->element != NULL ) {
-
-		int r;
-		r = hdfile_set_image(hdfile, iargs->element);
-		if ( r ) {
-			ERROR("Couldn't select path '%s'\n", iargs->element);
-			hdfile_close(hdfile);
-			return;
-		}
-
-	} else {
-
-		int r;
-		r = hdfile_set_first_image(hdfile, "/");
-		if ( r ) {
-			ERROR("Couldn't select first path\n");
-			hdfile_close(hdfile);
-			return;
-		}
-
-	}
-
-	check = hdf5_read(hdfile, &image, 1);
-	if ( check ) {
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( (image.width != image.det->max_fs + 1 )
-	  || (image.height != image.det->max_ss + 1))
-	{
-		ERROR("Image size doesn't match geometry size"
-			" - rejecting image.\n");
-		ERROR("Image size: %i,%i.  Geometry size: %i,%i\n",
-		      image.width, image.height,
-		      image.det->max_fs + 1, image.det->max_ss + 1);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	fill_in_values(image.det, hdfile);
-	fill_in_beam_parameters(image.beam, hdfile);
-
-	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
-
-	if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
-		/* Error message covers a silly value in the beam file or in
-		 * the HDF5 file. */
-		ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
-		      image.lambda, image.beam->photon_energy, image.filename);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( config_cmfilter ) {
-		filter_cm(&image);
-	}
-
-	/* Take snapshot of image after CM subtraction but before
-	 * the aggressive noise filter. */
-	data_size = image.width * image.height * sizeof(float);
-	data_for_measurement = malloc(data_size);
-
-	if ( config_noisefilter ) {
-		filter_noise(&image, data_for_measurement);
-	} else {
-		memcpy(data_for_measurement, image.data, data_size);
-	}
-
-	switch ( iargs->peaks ) {
-
-		case PEAK_HDF5:
-		// Get peaks from HDF5
-		if (get_peaks(&image, hdfile,
-			iargs->hdf5_peak_path)) {
-			ERROR("Failed to get peaks from HDF5 file.\n");
-		}
-		if ( !iargs->no_revalidate ) {
-			validate_peaks(&image, iargs->min_int_snr,
-				       iargs->ir_inn, iargs->ir_mid,
-				       iargs->ir_out, iargs->use_saturated);
-		}
-		break;
-
-		case PEAK_ZAEF:
-		search_peaks(&image, iargs->threshold,
-		             iargs->min_gradient, iargs->min_snr,
-		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
-		             iargs->use_saturated);
-		break;
-
-	}
-
-	/* Get rid of noise-filtered version at this point
-	 * - it was strictly for the purposes of peak detection. */
-	free(image.data);
-	image.data = data_for_measurement;
-
-	/* Calculate orientation matrix (by magic) */
-	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);
-
-	if ( image.indexed_cell != NULL ) {
-
-		pargs->indexable = 1;
-
-		if ( iargs->integrate_found ) {
-			image.reflections = select_intersections(&image,
-			                                    image.indexed_cell);
-		} else {
-			image.reflections = find_intersections(&image,
-			                                    image.indexed_cell);
-		}
-
-		if (image.reflections != NULL) {
-			integrate_reflections(&image,
-			                      iargs->config_closer,
-			                      iargs->config_bgsub,
-			                      iargs->min_int_snr,
-			                      iargs->ir_inn,
-			                      iargs->ir_mid,
-			                      iargs->ir_out,
-			                      iargs->integrate_saturated);
-		}
-
-	} else {
-		image.reflections = NULL;
-	}
-
-	write_chunk(ofh, &image, hdfile, iargs->stream_flags);
-	fprintf(ofh, "END\n");
-	fflush(ofh);
-
-	/* Only free cell if found */
-	cell_free(image.indexed_cell);
-
-	reflist_free(image.reflections);
-	free(image.data);
-	if ( image.flags != NULL ) free(image.flags);
-	image_feature_list_free(image.features);
-	hdfile_close(hdfile);
-}
-
-
 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 +224,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 +241,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);
@@ -442,11 +277,19 @@ static time_t get_monotonic_seconds()
 
 #endif
 
+size_t vol = 0;
+
+
+static ssize_t lwrite(int fd, const char *a)
+{
+	size_t l = strlen(a);
+	return write(fd, a, l);
+}
+
 
-static int pump_chunk(FILE *fh, FILE *ofh)
+static int pump_chunk(FILE *fh, int ofd)
 {
 	int chunk_started = 0;
-	int chunk_finished = 0;
 
 	do {
 
@@ -457,65 +300,121 @@ static int pump_chunk(FILE *fh, FILE *ofh)
 		if ( rval == NULL ) {
 
 			if ( feof(fh) ) {
-				/* Process died */
+				/* Whoops, connection lost */
 				if ( chunk_started ) {
 					ERROR("EOF during chunk!\n");
-					fprintf(ofh, "Chunk is unfinished!\n");
-				}
+					lwrite(ofd, "Unfinished chunk!\n");
+					lwrite(ofd, CHUNK_END_MARKER"\n");
+				} /* else normal end of output */
 				return 1;
-			} else {
-				ERROR("fgets() failed: %s\n", strerror(errno));
 			}
-			chunk_finished = 1;
-			continue;
 
-		}
+			ERROR("fgets() failed: %s\n", strerror(errno));
+			return 1;
 
-		if ( strcmp(line, "END\n") == 0 ) {
-			chunk_finished = 1;
-		} else {
-			chunk_started = 1;
-			fprintf(ofh, "%s", line);
 		}
 
-	} while ( !chunk_finished );
+		if ( strcmp(line, "FLUSH\n") == 0 ) break;
+		lwrite(ofd, line);
+
+		if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) break;
+		if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) break;
 
+	} while ( 1 );
 	return 0;
 }
 
 
-static void *run_reader(void *sbv)
+/* Add an fd to the list of pipes to be read from */
+static void add_pipe(struct sb_reader *rd, int fd)
 {
-	struct sandbox *sb = sbv;
-	int done = 0;
+	int *fds_new;
+	FILE **fhs_new;
+	int slot;
 
-	while ( !done ) {
+	pthread_mutex_lock(&rd->lock);
+
+	fds_new = realloc(rd->fds, (rd->n_read+1)*sizeof(int));
+	if ( fds_new == NULL ) {
+		ERROR("Failed to allocate memory for new pipe.\n");
+		return;
+	}
+
+	fhs_new = realloc(rd->fhs, (rd->n_read+1)*sizeof(FILE *));
+	if ( fhs_new == NULL ) {
+		ERROR("Failed to allocate memory for new FH.\n");
+		return;
+	}
+
+	rd->fds = fds_new;
+	rd->fhs = fhs_new;
+	slot = rd->n_read;
+
+	rd->fds[slot] = fd;
+
+	rd->fhs[slot] = fdopen(fd, "r");
+	if ( rd->fhs[slot] == NULL ) {
+		ERROR("Couldn't fdopen() stream!\n");
+		return;
+	}
+
+	rd->n_read++;
+
+	pthread_mutex_unlock(&rd->lock);
+}
+
+
+/* Assumes that the caller is already holding rd->lock! */
+static void remove_pipe(struct sb_reader *rd, int d)
+{
+	int i;
+
+	for ( i=d; i<rd->n_read; i++ ) {
+		if ( i < rd->n_read-1 ) {
+			rd->fds[i] = rd->fds[i+1];
+			rd->fhs[i] = rd->fhs[i+1];
+		} /* else don't bother */
+	}
+
+	rd->n_read--;
+
+	/* We don't bother shrinking the arrays */
+}
+
+
+static void *run_reader(void *rdv)
+{
+	struct sb_reader *rd = rdv;
+
+	while ( 1 ) {
 
 		int r, i;
 		struct timeval tv;
 		fd_set fds;
 		int fdmax;
 
-		tv.tv_sec = 5;
+		/* Exit when:
+		 *     - No fhs left open to read from
+		 * AND - Main thread says "done" */
+		if ( (rd->n_read == 0) && rd->done ) break;
+
+		tv.tv_sec = 1;
 		tv.tv_usec = 0;
 
 		FD_ZERO(&fds);
 		fdmax = 0;
-		lock_sandbox(sb);
-		for ( i=0; i<sb->n_proc; i++ ) {
+		pthread_mutex_lock(&rd->lock);
+		for ( i=0; i<rd->n_read; i++ ) {
 
 			int fd;
 
-			if ( !sb->running[i] ) continue;
-
-			fd = sb->stream_pipe_read[i];
+			fd = rd->fds[i];
 
 			FD_SET(fd, &fds);
 			if ( fd > fdmax ) fdmax = fd;
 
 		}
-
-		unlock_sandbox(sb);
+		pthread_mutex_unlock(&rd->lock);
 
 		r = select(fdmax+1, &fds, NULL, NULL, &tv);
 
@@ -526,38 +425,75 @@ static void *run_reader(void *sbv)
 			continue;
 		}
 
-		if ( r == 0 ) continue; /* Nothing this time.  Try again */
+		pthread_mutex_lock(&rd->lock);
+		for ( i=0; i<rd->n_read; i++ ) {
 
-		lock_sandbox(sb);
-		for ( i=0; i<sb->n_proc; i++ ) {
+			if ( !FD_ISSET(rd->fds[i], &fds) ) {
+				continue;
+			}
 
-			if ( !sb->running[i] ) continue;
+			/* If the chunk cannot be read, assume the connection
+			 * is broken and that the process will die soon. */
+			if ( pump_chunk(rd->fhs[i], rd->ofd) ) {
+				/* remove_pipe() assumes that the caller is
+				 * holding rd->lock ! */
+				remove_pipe(rd, i);
+			}
 
-			if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) continue;
+		}
+		pthread_mutex_unlock(&rd->lock);
+
+	}
+
+	return NULL;
+}
 
-			if ( pump_chunk(sb->fhs[i], sb->ofh) ) {
-				sb->running[i] = 0;
-			}
 
+static int create_temporary_folder(signed int n)
+{
+	int r;
+	char tmp[64];
+	struct stat s;
+
+	if ( n < 0 ) {
+		snprintf(tmp, 63, "indexamajig.%i", getpid());
+	} else {
+		snprintf(tmp, 63, "worker.%i", n);
+	}
+
+	if ( stat(tmp, &s) == -1 ) {
+		if ( errno != ENOENT ) {
+			ERROR("Failed to stat temporary folder.\n");
+			return 1;
 		}
 
-		done = 1;
-		for ( i=0; i<sb->n_proc; i++ ) {
-			if ( sb->running[i] ) done = 0;
+		r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+		if ( r ) {
+			ERROR("Failed to create temporary folder: %s\n",
+			      strerror(errno));
+			return 1;
 		}
-		unlock_sandbox(sb);
 
 	}
 
-	return NULL;
+	r = chdir(tmp);
+	if ( r ) {
+		ERROR("Failed to chdir to temporary folder: %s\n",
+		      strerror(errno));
+		return 1;
+	}
+
+	return 0;
 }
 
 
-static void start_worker_process(struct sandbox *sb, int slot)
+static void start_worker_process(struct sandbox *sb, int slot,
+                                 int argc, char *argv[])
 {
 	pid_t p;
 	int filename_pipe[2];
 	int result_pipe[2];
+	int stream_pipe[2];
 
 	if ( pipe(filename_pipe) == - 1 ) {
 		ERROR("pipe() failed!\n");
@@ -569,6 +505,11 @@ static void start_worker_process(struct sandbox *sb, int slot)
 		return;
 	}
 
+	if ( pipe(stream_pipe) == - 1 ) {
+		ERROR("pipe() failed!\n");
+		return;
+	}
+
 	p = fork();
 	if ( p == -1 ) {
 		ERROR("fork() failed!\n");
@@ -577,7 +518,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;
@@ -592,6 +533,8 @@ static void start_worker_process(struct sandbox *sb, int slot)
 			return;
 		}
 
+		create_temporary_folder(slot);
+
 		/* Free resources which will not be needed by worker */
 		for ( j=0; j<sb->n_proc; j++ ) {
 			if ( (j != slot) && (sb->running[j]) ) {
@@ -600,7 +543,9 @@ static void start_worker_process(struct sandbox *sb, int slot)
 		}
 		for ( j=0; j<sb->n_proc; j++ ) {
 			if ( (j != slot) && (sb->running[j]) ) {
-				fclose(sb->result_fhs[j]);
+				if ( sb->result_fhs[j] != NULL ) {
+					fclose(sb->result_fhs[j]);
+				}
 				close(sb->filename_pipes[j]);
 			}
 		}
@@ -614,10 +559,12 @@ 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(stream_pipe[1]);
+		write_command(st, argc, argv);
+		write_line(st, "FLUSH");
 		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]);
@@ -630,14 +577,11 @@ static void start_worker_process(struct sandbox *sb, int slot)
 	 * and the 'read' end of the result pipe. */
 	sb->pids[slot] = p;
 	sb->running[slot] = 1;
+	add_pipe(sb->reader, stream_pipe[0]);
 	close(filename_pipe[0]);
 	close(result_pipe[1]);
+	close(stream_pipe[1]);
 	sb->filename_pipes[slot] = filename_pipe[1];
-	sb->fhs[slot] = fdopen(sb->stream_pipe_read[slot], "r");
-	if ( sb->fhs[slot] == NULL ) {
-		ERROR("Couldn't fdopen() stream!\n");
-		return;
-	}
 
 	sb->result_fhs[slot] = fdopen(result_pipe[0], "r");
 	if ( sb->result_fhs[slot] == NULL ) {
@@ -685,7 +629,7 @@ static void handle_zombie(struct sandbox *sb)
 				STATUS("Last filename was: %s\n",
 				       sb->last_filename[i]);
 				sb->n_processed++;
-				start_worker_process(sb, i);
+				start_worker_process(sb, i, 0, NULL);
 			}
 
 		}
@@ -697,7 +641,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)
+                    int ofd, int argc, char *argv[])
 {
 	int i;
 	int allDone;
@@ -712,48 +656,41 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		return;
 	}
 
-	sb->n_indexable = 0;
+	sb->reader = calloc(1, sizeof(struct sb_reader));
+	if ( sb->reader == NULL ) {
+		ERROR("Couldn't allocate memory for SB reader.\n");
+		free(sb);
+		return;
+	}
+
+	pthread_mutex_init(&sb->lock, NULL);
+	pthread_mutex_init(&sb->reader->lock, NULL);
+
 	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);
+	sb->reader->fds = NULL;
+	sb->reader->fhs = NULL;
+	sb->reader->ofd = ofd;
 
-	sb->stream_pipe_read = calloc(n_proc, sizeof(int));
 	sb->stream_pipe_write = calloc(n_proc, sizeof(int));
-	if ( sb->stream_pipe_read == NULL ) {
-		ERROR("Couldn't allocate memory for pipes.\n");
-		return;
-	}
 	if ( sb->stream_pipe_write == NULL ) {
 		ERROR("Couldn't allocate memory for pipes.\n");
 		return;
 	}
 
-	for ( i=0; i<n_proc; i++ ) {
-
-		int stream_pipe[2];
-
-		if ( pipe(stream_pipe) == - 1 ) {
-			ERROR("pipe() failed!\n");
-			return;
-		}
-
-		sb->stream_pipe_read[i] = stream_pipe[0];
-		sb->stream_pipe_write[i] = stream_pipe[1];
-
-	}
-
 	lock_sandbox(sb);
 	sb->filename_pipes = calloc(n_proc, sizeof(int));
 	sb->result_fhs = calloc(n_proc, sizeof(FILE *));
 	sb->pids = calloc(n_proc, sizeof(pid_t));
 	sb->running = calloc(n_proc, sizeof(int));
-	sb->fhs = calloc(sb->n_proc, sizeof(FILE *));
 	if ( sb->filename_pipes == NULL ) {
 		ERROR("Couldn't allocate memory for pipes.\n");
 		return;
@@ -776,17 +713,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		ERROR("Couldn't allocate memory for last filename list.\n");
 		return;
 	}
-	if ( sb->fhs == NULL ) {
-		ERROR("Couldn't allocate memory for file handles!\n");
-		return;
-	}
 	unlock_sandbox(sb);
 
-	if ( pthread_create(&reader_thread, NULL, run_reader, (void *)sb) ) {
-		ERROR("Failed to create reader thread.\n");
-		return;
-	}
-
 	if ( pipe(signal_pipe) == -1 ) {
 		ERROR("Failed to create signal pipe.\n");
 		return;
@@ -802,13 +730,23 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		return;
 	}
 
+	if ( create_temporary_folder(-1) ) return;
+
 	/* Fork the right number of times */
 	lock_sandbox(sb);
 	for ( i=0; i<n_proc; i++ ) {
-		start_worker_process(sb, i);
+		start_worker_process(sb, i, argc, argv);
 	}
 	unlock_sandbox(sb);
 
+	/* Start reader thread after forking, so that things are definitely
+	 * "running" */
+	if ( pthread_create(&reader_thread, NULL, run_reader,
+	                    (void *)sb->reader) ) {
+		ERROR("Failed to create reader thread.\n");
+		return;
+	}
+
 	allDone = 0;
 	while ( !allDone ) {
 
@@ -828,9 +766,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 
 			int fd;
 
-			if ( !sb->running[i] ) {
-				continue;
-			}
+			if ( sb->result_fhs[i] == NULL) continue;
 
 			fd = fileno(sb->result_fhs[i]);
 			FD_SET(fd, &fds);
@@ -844,14 +780,11 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 
 		r = select(fdmax+1, &fds, NULL, NULL, &tv);
 		if ( r == -1 ) {
-			if ( errno != EINTR ) {
-				ERROR("select() failed: %s\n", strerror(errno));
-			}
-			continue;
+			if ( errno == EINTR ) continue;
+			ERROR("select() failed: %s\n", strerror(errno));
+			break;
 		}
 
-		if ( r == 0 ) continue; /* No progress this time.  Try again */
-
 		if ( FD_ISSET(signal_pipe[0], &fds) ) {
 
 			char d;
@@ -869,14 +802,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 			int fd;
 			char *eptr;
 
-			if ( !sb->running[i] ) {
-				continue;
-			}
+			if ( sb->result_fhs[i] == NULL ) continue;
 
 			fd = fileno(sb->result_fhs[i]);
-			if ( !FD_ISSET(fd, &fds) ) {
-				continue;
-			}
+			if ( !FD_ISSET(fd, &fds) ) continue;
 
 			rval = fgets(results, 1024, sb->result_fhs[i]);
 			if ( rval == NULL ) {
@@ -884,6 +813,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 					ERROR("fgets() failed: %s\n",
 					      strerror(errno));
 				}
+				sb->result_fhs[i] = NULL;
 				continue;
 			}
 
@@ -895,8 +825,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 +865,16 @@ 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,
+			STATUS("%4i indexable out of %4i processed, "
+			       "%4i crystals so far. "
+			       "%4i images processed since the last message.\n",
+			       sb->n_hadcrystals, sb->n_processed,
+			       sb->n_crystals,
 			       sb->n_processed - sb->n_processed_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;
 
 		}
@@ -947,14 +885,18 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		for ( i=0; i<n_proc; i++ ) {
 			if ( sb->running[i] ) allDone = 0;
 		}
-
 		unlock_sandbox(sb);
 
 	}
 
 	fclose(fh);
 
-	pthread_mutex_destroy(&sb->lock);
+	/* Indicate to the reader thread that we are done */
+	pthread_mutex_lock(&sb->reader->lock);
+	sb->reader->done = 1;
+	pthread_mutex_unlock(&sb->reader->lock);
+
+	pthread_join(reader_thread, NULL);
 
 	for ( i=0; i<n_proc; i++ ) {
 		int status;
@@ -963,21 +905,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 
 	for ( i=0; i<n_proc; i++ ) {
 		close(sb->filename_pipes[i]);
-		fclose(sb->result_fhs[i]);
+		if ( sb->result_fhs[i] != NULL ) fclose(sb->result_fhs[i]);
 	}
 
-	for ( i=0; i<sb->n_proc; i++ ) {
-		fclose(sb->fhs[i]);
-	}
-	free(sb->fhs);
+	free(sb->running);
 	free(sb->filename_pipes);
 	free(sb->result_fhs);
 	free(sb->pids);
-	free(sb->running);
 
-	if ( ofh != stdout ) fclose(ofh);
+	pthread_mutex_destroy(&sb->lock);
 
-	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);
 
+	free(sb);
 }
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 872f84ad..9248b226 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,60 +31,12 @@
  *
  */
 
-enum {
-	PEAK_ZAEF,
-	PEAK_HDF5,
-};
-
-
-/* Information about the indexing process which is common to all patterns */
-struct index_args
-{
-	UnitCell *cell;
-	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;
-	float min_snr;
-	double min_int_snr;
-	struct detector *det;
-	IndexingMethod *indm;
-	IndexingPrivate **ipriv;
-	int peaks;                /* Peak detection method */
-	int cellr;
-	float tols[4];
-	struct beam_params *beam;
-	char *element;
-	char *hdf5_peak_path;
-	double ir_inn;
-	double ir_mid;
-	double ir_out;
-	struct copy_hdf5_field *copyme;
-	int integrate_saturated;
-	int use_saturated;
-	int no_revalidate;
-	int integrate_found;
-};
-
-
-
-
-/* Information about the indexing process for one pattern */
-struct pattern_args
-{
-	/* "Input" */
-	char *filename;
-
-	/* "Output" */
-	int indexable;
-};
+#include "index.h"
+#include "stream.h"
+#include "cell.h"
+#include "process_image.h"
 
 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, int streamfd,
+                           int argc, char *argv[]);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index accfa788..9b868b9f 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,13 +3,13 @@
  *
  * Index patterns, output hkl+intensity etc.
  *
- * 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
@@ -44,6 +44,9 @@
 #include <getopt.h>
 #include <hdf5.h>
 #include <gsl/gsl_errno.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
 
 #ifdef HAVE_CLOCK_GETTIME
 #include <time.h>
@@ -81,11 +84,8 @@ static void show_help(const char *s)
 "                           Default: indexamajig.stream\n"
 "\n"
 "     --indexing=<methods> Use 'methods' for indexing.  Provide one or more\n"
-"                           methods separated by commas.  Choose from:\n"
-"                            none     : no indexing (default)\n"
-"                            dirax    : invoke DirAx\n"
-"                            mosflm   : invoke MOSFLM (DPS)\n"
-"                            reax     : DPS algorithm with known unit cell\n"
+"                           methods separated by commas.\n"
+"                           See 'man indexamajig' for details.\n"
 " -g. --geometry=<file>    Get detector geometry from file.\n"
 " -b, --beam=<file>        Get beam parameters from file (provides nominal\n"
 "                           wavelength value if no per-shot value is found in\n"
@@ -101,28 +101,8 @@ static void show_help(const char *s)
 "     --hdf5-peaks=<p>     Find peaks table in HDF5 file here.\n"
 "                           Default: /processing/hitfinder/peakinfo\n"
 "\n\n"
-"You can control what information is included in the output stream using\n"
-"' --record=<flag1>,<flag2>,<flag3>' and so on.  Possible flags are:\n\n"
-" integrated        Include a list of reflection intensities, produced by\n"
-"                    integrating around predicted peak locations.\n"
-"\n"
-" peaks             Include peak locations and intensities from the peak\n"
-"                    search.\n"
-"\n"
-" peaksifindexed    As 'peaks', but only if the pattern could be indexed.\n"
-"\n"
-" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
-"\n\n"
-"The default is '--record=integrated'.\n"
-"\n\n"
 "For more control over the process, you might need:\n\n"
-"  --cell-reduction=<m>  Use <m> as the cell reduction method. Choose from:\n"
-"                         none    : no matching, just use the raw cell.\n"
-"                         reduce  : full cell reduction.\n"
-"                         compare : match by at most changing the order of\n"
-"                                   the indices.\n"
-"                         compare_ab : compare 'a' and 'b' lengths only.\n"
-"    --tolerance=<tol>   Set the tolerances for cell reduction.\n"
+"    --tolerance=<tol>   Set the tolerances for cell comparison.\n"
 "                          Default: 5,5,5,1.5.\n"
 "    --filter-cm         Perform common-mode noise subtraction on images\n"
 "                         before proceeding.  Intensities will be extracted\n"
@@ -134,7 +114,7 @@ static void show_help(const char *s)
 "    --no-sat-corr       Don't correct values of saturated peaks using a\n"
 "                         table included in the HDF5 file.\n"
 "    --threshold=<n>     Only accept peaks above <n> ADU.  Default: 800.\n"
-"    --min-gradient=<n>  Minimum gradient for Zaefferer peak search.\n"
+"    --min-gradient=<n>  Minimum squared gradient for Zaefferer peak search.\n"
 "                         Default: 100,000.\n"
 "    --min-snr=<n>       Minimum signal-to-noise ratio for peaks.\n"
 "                         Default: 5.\n"
@@ -144,21 +124,20 @@ static void show_help(const char *s)
 "-e, --image=<element>   Use this image from the HDF5 file.\n"
 "                          Example: /data/data0.\n"
 "                          Default: The first one found.\n"
+"    --res-cutoff        Estimate the resolution limit of each pattern, and\n"
+"                         don't integrate reflections further out.\n"
 "\n"
 "\nFor time-resolved stuff, you might want to use:\n\n"
 "     --copy-hdf5-field <f>  Copy the value of field <f> into the stream. You\n"
 "                             can use this option as many times as you need.\n"
 "\n"
-"\nOptions for greater performance or verbosity:\n\n"
-"     --verbose            Be verbose about indexing.\n"
+"\nOptions for greater performance:\n\n"
 " -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
 "\n"
 "\nOptions you probably won't need:\n\n"
 "     --no-check-prefix    Don't attempt to correct the --prefix.\n"
 "     --closer-peak        Don't integrate from the location of a nearby peak\n"
 "                           instead of the predicted spot.  Don't use.\n"
-"     --insane             Don't check that the reduced cell accounts for at\n"
-"                           least 10%% of the located peaks.\n"
 "     --no-bg-sub          Don't subtract local background estimates from\n"
 "                           integrated intensities.\n"
 "     --use-saturated      During the initial peak search, don't reject\n"
@@ -170,124 +149,110 @@ static void show_help(const char *s)
 "     --integrate-found    Skip the spot prediction step, and just integrate\n"
 "                           the intensities of the spots found by the initial\n"
 "                           peak search.\n"
+"     --no-peaks-in-stream Do not record peak search results in the stream.\n"
+"     --no-refls-in-stream Do not record integrated reflections in the stream.\n"
 );
 }
 
 
-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;
+	int ofd;
 	char *rval = NULL;
-	int config_noindex = 0;
-	int config_cmfilter = 0;
-	int config_noisefilter = 0;
-	int config_verbose = 0;
-	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;
-	float min_gradient = 100000.0;
-	float min_snr = 5.0;
-	double min_int_snr = -INFINITY;
-	struct detector *det;
-	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;
 	char prepare_filename[1024];
 	char *use_this_one_instead;
 	struct index_args iargs;
-	struct beam_params *beam = NULL;
-	char *element = NULL;
-	double nominal_photon_energy;
-	int stream_flags = STREAM_INTEGRATED;
-	char *hdf5_peak_path = NULL;
-	struct copy_hdf5_field *copyme;
 	char *intrad = NULL;
-	float ir_inn = 4.0;
-	float ir_mid = 5.0;
-	float ir_out = 7.0;
-	int integrate_saturated = 0;
-	int use_saturated = 0;
-	int no_revalidate = 0;
-	int integrate_found = 0;
-
-	copyme = new_copy_hdf5_field_list();
-	if ( copyme == NULL ) {
+
+	/* Defaults */
+	iargs.cell = NULL;
+	iargs.cmfilter = 0;
+	iargs.noisefilter = 0;
+	iargs.satcorr = 1;
+	iargs.closer = 0;
+	iargs.bgsub = 1;
+	iargs.tols[0] = 5.0;
+	iargs.tols[1] = 5.0;
+	iargs.tols[2] = 5.0;
+	iargs.tols[3] = 1.5;
+	iargs.threshold = 800.0;
+	iargs.min_gradient = 100000.0;
+	iargs.min_snr = 5.0;
+	iargs.min_int_snr = -INFINITY;
+	iargs.det = NULL;
+	iargs.peaks = PEAK_ZAEF;
+	iargs.beam = NULL;
+	iargs.element = NULL;
+	iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
+	iargs.copyme = NULL;
+	iargs.ir_inn = 4.0;
+	iargs.ir_mid = 5.0;
+	iargs.ir_out = 7.0;
+	iargs.use_saturated = 0;
+	iargs.integrate_saturated = 0;
+	iargs.no_revalidate = 0;
+	iargs.integrate_found = 0;
+	iargs.stream_peaks = 1;
+	iargs.stream_refls = 1;
+	iargs.res_cutoff = 0;
+	iargs.copyme = new_copy_hdf5_field_list();
+	if ( iargs.copyme == NULL ) {
 		ERROR("Couldn't allocate HDF5 field list.\n");
 		return 1;
 	}
+	iargs.indm = NULL;  /* No default */
+	iargs.ipriv = NULL;  /* No default */
 
 	/* Long options */
 	const struct option longopts[] = {
+
+		/* Options with long and short versions */
 		{"help",               0, NULL,               'h'},
 		{"input",              1, NULL,               'i'},
 		{"output",             1, NULL,               'o'},
-		{"no-index",           0, &config_noindex,     1},
 		{"indexing",           1, NULL,               'z'},
 		{"geometry",           1, NULL,               'g'},
 		{"beam",               1, NULL,               'b'},
-		{"filter-cm",          0, &config_cmfilter,    1},
-		{"filter-noise",       0, &config_noisefilter, 1},
-		{"verbose",            0, &config_verbose,     1},
 		{"pdb",                1, NULL,               'p'},
 		{"prefix",             1, NULL,               'x'},
-		{"no-sat-corr",        0, &config_satcorr,     0},
-		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
 		{"threshold",          1, NULL,               't'},
-		{"no-check-prefix",    0, &config_checkprefix, 0},
-		{"no-closer-peak",     0, &config_closer,      0},
-		{"closer-peak",        0, &config_closer,      1},
-		{"insane",             0, &config_insane,      1},
 		{"image",              1, NULL,               'e'},
-		{"basename",           0, &config_basename,    1},
-		{"bg-sub",             0, &config_bgsub,       1}, /* Compat */
-		{"no-bg-sub",          0, &config_bgsub,       0},
 
+		/* Long-only options with no arguments */
+		{"filter-cm",          0, &iargs.cmfilter,           1},
+		{"filter-noise",       0, &iargs.noisefilter,        1},
+		{"no-sat-corr",        0, &iargs.satcorr,            0},
+		{"sat-corr",           0, &iargs.satcorr,            1},
+		{"no-check-prefix",    0, &config_checkprefix,       0},
+		{"no-closer-peak",     0, &iargs.closer,             0},
+		{"closer-peak",        0, &iargs.closer,             1},
+		{"basename",           0, &config_basename,          1},
+		{"bg-sub",             0, &iargs.bgsub,              1},
+		{"no-bg-sub",          0, &iargs.bgsub,              0},
+		{"no-peaks-in-stream", 0, &iargs.stream_peaks,       0},
+		{"no-refls-in-stream", 0, &iargs.stream_refls,       0},
+		{"res-cutoff",         0, &iargs.res_cutoff,         1},
+		{"integrate-saturated",0, &iargs.integrate_saturated,1},
+		{"use-saturated",      0, &iargs.use_saturated,      1},
+		{"no-revalidate",      0, &iargs.no_revalidate,      1},
+		{"integrate-found",    0, &iargs.integrate_found,    1},
+
+		/* Long-only options with arguments */
 		{"peaks",              1, NULL,                2},
 		{"cell-reduction",     1, NULL,                3},
 		{"min-gradient",       1, NULL,                4},
@@ -302,10 +267,6 @@ int main(int argc, char *argv[])
 		{"tolerance",          1, NULL,               13},
 		{"int-radius",         1, NULL,               14},
 
-		{"integrate-saturated",0, &integrate_saturated,1},
-		{"use-saturated",      0, &use_saturated,      1},
-		{"no-revalidate",      0, &no_revalidate,      1},
-		{"integrate-found",    0, &integrate_found,    1},
 		{0, 0, NULL, 0}
 	};
 
@@ -344,16 +305,21 @@ int main(int argc, char *argv[])
 			break;
 
 			case 'g' :
-			geometry = strdup(optarg);
+			iargs.det = get_detector_geometry(optarg);
+			if ( iargs.det == NULL ) {
+				ERROR("Failed to read detector geometry from "
+				      "'%s'\n", optarg);
+				return 1;
+			}
 			break;
 
 			case 't' :
-			threshold = strtof(optarg, NULL);
+			iargs.threshold = strtof(optarg, NULL);
 			break;
 
 			case 'b' :
-			beam = get_beam_parameters(optarg);
-			if ( beam == NULL ) {
+			iargs.beam = get_beam_parameters(optarg);
+			if ( iargs.beam == NULL ) {
 				ERROR("Failed to load beam parameters"
 				      " from '%s'\n", optarg);
 				return 1;
@@ -361,7 +327,7 @@ int main(int argc, char *argv[])
 			break;
 
 			case 'e' :
-			element = strdup(optarg);
+			iargs.element = strdup(optarg);
 			break;
 
 			case 2 :
@@ -369,17 +335,24 @@ int main(int argc, char *argv[])
 			break;
 
 			case 3 :
-			scellr = strdup(optarg);
-			break;
+			ERROR("The option '--cell-reduction' is no longer "
+			      "used.\n"
+			      "The complete indexing behaviour is now "
+			      "controlled using '--indexing'.\n"
+			      "See 'man indexamajig' for details of the "
+			      "available methods.\n");
+			return 1;
 
 			case 4 :
-			min_gradient = strtof(optarg, NULL);
+			iargs.min_gradient = strtof(optarg, NULL);
 			break;
 
 			case 5 :
-			stream_flags = parse_stream_flags(optarg);
-			if ( stream_flags < 0 ) return 1;
-			break;
+			ERROR("The option '--record' is no longer used.\n"
+			      "Use '--no-peaks-in-stream' and"
+			      "'--no-refls-in-stream' if you need to control"
+			      "the contents of the stream.\n");
+			return 1;
 
 			case 6 :
 			case 7 :
@@ -389,19 +362,20 @@ int main(int argc, char *argv[])
 			break;
 
 			case 9 :
-			hdf5_peak_path = strdup(optarg);
+			free(iargs.hdf5_peak_path);
+			iargs.hdf5_peak_path = strdup(optarg);
 			break;
 
 			case 10 :
-			add_copy_hdf5_field(copyme, optarg);
+			add_copy_hdf5_field(iargs.copyme, optarg);
 			break;
 
 			case 11 :
-			min_snr = strtof(optarg, NULL);
+			iargs.min_snr = strtof(optarg, NULL);
 			break;
 
 			case 12 :
-			min_int_snr = strtof(optarg, NULL);
+			iargs.min_int_snr = strtof(optarg, NULL);
 			break;
 
 			case 13 :
@@ -440,30 +414,15 @@ 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");
-	}
-
 	if ( speaks == NULL ) {
 		speaks = strdup("zaef");
 		STATUS("You didn't specify a peak detection method.\n");
 		STATUS("I'm using 'zaef' for you.\n");
 	}
 	if ( strcmp(speaks, "zaef") == 0 ) {
-		peaks = PEAK_ZAEF;
+		iargs.peaks = PEAK_ZAEF;
 	} else if ( strcmp(speaks, "hdf5") == 0 ) {
-		peaks = PEAK_HDF5;
+		iargs.peaks = PEAK_HDF5;
 	} else {
 		ERROR("Unrecognised peak detection method '%s'\n", speaks);
 		return 1;
@@ -483,57 +442,29 @@ 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",
-		             &tols[0], &tols[1], &tols[2], &tols[3] );
+		             &iargs.tols[0], &iargs.tols[1],
+		             &iargs.tols[2], &iargs.tols[3]);
 		if ( ttt != 4 ) {
 			ERROR("Invalid parameters for '--tolerance'\n");
 			return 1;
@@ -543,7 +474,8 @@ int main(int argc, char *argv[])
 
 	if ( intrad != NULL ) {
 		int r;
-		r = sscanf(intrad, "%f,%f,%f", &ir_inn, &ir_mid, &ir_out);
+		r = sscanf(intrad, "%f,%f,%f",
+		           &iargs.ir_inn, &iargs.ir_mid, &iargs.ir_out);
 		if ( r != 3 ) {
 			ERROR("Invalid parameters for '--int-radius'\n");
 			return 1;
@@ -555,43 +487,38 @@ int main(int argc, char *argv[])
 		       " probably not appropriate for your patterns.\n");
 	}
 
-	if ( geometry == NULL ) {
-		ERROR("You need to specify a geometry file with --geometry\n");
+	if ( iargs.det == NULL ) {
+		ERROR("You need to provide a geometry file (please read the"
+		      " manual for more details).\n");
 		return 1;
 	}
 
-	det = get_detector_geometry(geometry);
-	if ( det == NULL ) {
-		ERROR("Failed to read detector geometry from '%s'\n", geometry);
+	if ( iargs.beam == NULL ) {
+		ERROR("You need to provide a beam parameters file (please read"
+		      " the manual for more details).\n");
 		return 1;
 	}
-	free(geometry);
 
 	if ( pdb != NULL ) {
-		cell = load_cell_from_pdb(pdb);
-		if ( cell == NULL ) {
+		iargs.cell = load_cell_from_pdb(pdb);
+		if ( iargs.cell == NULL ) {
 			ERROR("Couldn't read unit cell (from %s)\n", pdb);
 			return 1;
 		}
 		free(pdb);
-		cell_print(cell);
+		STATUS("This is what I understood your unit cell to be:\n");
+		cell_print(iargs.cell);
 	} else {
 		STATUS("No unit cell given.\n");
-		cell = NULL;
+		iargs.cell = NULL;
 	}
 
-	write_stream_header(ofh, argc, argv);
-
-	if ( beam != NULL ) {
-		nominal_photon_energy = beam->photon_energy;
-	} else {
-		STATUS("No beam parameters file was given, so I'm taking the"
-		       " nominal photon energy to be 2 keV.\n");
-		ERROR("I'm also going to assume 1 ADU per photon, which is");
-		ERROR(" almost certainly wrong.  Peak sigmas will be"
-		      " incorrect.\n");
-		nominal_photon_energy = 2000.0;
+	ofd = open(outfile, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR);
+	if ( ofd == -1 ) {
+		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);
@@ -613,8 +540,8 @@ int main(int argc, char *argv[])
 
 	/* Prepare the indexer */
 	if ( indm != NULL ) {
-		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
-		                         nominal_photon_energy);
+		ipriv = prepare_indexing(indm, iargs.cell, prepare_filename,
+		                         iargs.det, iargs.beam, iargs.tols);
 		if ( ipriv == NULL ) {
 			ERROR("Failed to prepare indexing.\n");
 			return 1;
@@ -625,43 +552,11 @@ int main(int argc, char *argv[])
 
 	gsl_set_error_handler_off();
 
-	/* Static worker args */
-	iargs.cell = cell;
-	iargs.config_cmfilter = config_cmfilter;
-	iargs.config_noisefilter = config_noisefilter;
-	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];
-	iargs.tols[3] = tols[3];
-	iargs.threshold = threshold;
-	iargs.min_gradient = min_gradient;
-	iargs.min_snr = min_snr;
-	iargs.min_int_snr = min_int_snr;
-	iargs.det = det;
 	iargs.indm = indm;
 	iargs.ipriv = ipriv;
-	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;
-	iargs.ir_mid = ir_mid;
-	iargs.ir_out = ir_out;
-	iargs.use_saturated = use_saturated;
-	iargs.integrate_saturated = integrate_saturated;
-	iargs.no_revalidate = no_revalidate;
-	iargs.integrate_found = integrate_found;
 
 	create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
-	                    use_this_one_instead, ofh);
+	               use_this_one_instead, ofd, argc, argv);
 
 	free(prefix);
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index facf14ef..a4a2f875 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -54,11 +54,12 @@
 #define NBINS 50
 
 
-static void mess_up_cell(UnitCell *cell, double cnoise)
+static void mess_up_cell(Crystal *cr, double cnoise)
 {
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
+	UnitCell *cell = crystal_get_cell(cr);
 
 	//STATUS("Real:\n");
 	//cell_print(cell);
@@ -82,17 +83,18 @@ static void mess_up_cell(UnitCell *cell, double cnoise)
 
 /* For each reflection in "partial", fill in what the intensity would be
  * according to "full" */
-static void calculate_partials(RefList *partial, double osf,
+static void calculate_partials(Crystal *cr,
                                RefList *full, const SymOpList *sym,
                                int random_intensities,
                                pthread_mutex_t *full_lock,
                                unsigned long int *n_ref, double *p_hist,
-                               double *p_max, double max_q, UnitCell *cell)
+                               double *p_max, double max_q)
 {
 	Reflection *refl;
 	RefListIterator *iter;
+	double res;
 
-	for ( refl = first_refl(partial, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -137,16 +139,17 @@ static void calculate_partials(RefList *partial, double osf,
 			}
 		}
 
-		Ip = osf * p * If;
+		Ip = crystal_get_osf(cr) * p * If;
 
-		bin = NBINS*2.0*resolution(cell, h, k, l) / max_q;
+		res = resolution(crystal_get_cell(cr), h, k, l);
+		bin = NBINS*2.0*res/max_q;
 		if ( (bin < NBINS) && (bin>=0) ) {
 			p_hist[bin] += p;
 			n_ref[bin]++;
 			if ( p > p_max[bin] ) p_max[bin] = p;
 		} else {
 			STATUS("Reflection out of histogram range: %e %i %f\n",
-			       resolution(cell, h, k, l), bin,  p);
+			       res, bin,  p);
 		}
 
 		Ip = gaussian_noise(Ip, 100.0);
@@ -206,13 +209,14 @@ struct queue_args
 	unsigned long int n_ref[NBINS];
 	double p_max[NBINS];
 
-	FILE *stream;
+	Stream *stream;
 };
 
 
 struct worker_args
 {
 	struct queue_args *qargs;
+	Crystal *crystal;
 	struct image image;
 
 	/* Histogram for this image */
@@ -243,21 +247,31 @@ static void *create_job(void *vqargs)
 
 static void run_job(void *vwargs, int cookie)
 {
-	double osf;
 	struct quaternion orientation;
 	struct worker_args *wargs = vwargs;
 	struct queue_args *qargs = wargs->qargs;
 	int i;
+	Crystal *cr;
+	RefList *reflections;
 
-	osf = gaussian_noise(1.0, 0.3);
+	cr = crystal_new();
+	if ( cr == NULL ) {
+		ERROR("Failed to create crystal.\n");
+		return;
+	}
+	wargs->crystal = cr;
+	crystal_set_image(cr, &wargs->image);
+
+	crystal_set_osf(cr, gaussian_noise(1.0, 0.3));
+	crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
 
 	/* Set up a random orientation */
 	orientation = random_quaternion();
-	wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation);
+	crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
 
 	snprintf(wargs->image.filename, 255, "dummy.h5");
-	wargs->image.reflections = find_intersections(&wargs->image,
-	                                       wargs->image.indexed_cell);
+	reflections = find_intersections(&wargs->image, cr);
+	crystal_set_reflections(cr, reflections);
 
 	for ( i=0; i<NBINS; i++ ) {
 		wargs->n_ref[i] = 0;
@@ -265,14 +279,16 @@ static void run_job(void *vwargs, int cookie)
 		wargs->p_max[i] = 0.0;
 	}
 
-	calculate_partials(wargs->image.reflections, osf, qargs->full,
+	calculate_partials(cr, qargs->full,
 	                   qargs->sym, qargs->random_intensities,
 	                   &qargs->full_lock,
 	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
-	                   qargs->max_q, wargs->image.indexed_cell);
+	                   qargs->max_q);
 
 	/* Give a slightly incorrect cell in the stream */
-	mess_up_cell(wargs->image.indexed_cell, qargs->cnoise);
+	mess_up_cell(cr, qargs->cnoise);
+
+	image_add_crystal(&wargs->image, cr);
 }
 
 
@@ -282,7 +298,7 @@ static void finalise_job(void *vqargs, void *vwargs)
 	struct queue_args *qargs = vqargs;
 	int i;
 
-	write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED);
+	write_chunk(qargs->stream, &wargs->image, NULL, 0, 1);
 
 	for ( i=0; i<NBINS; i++ ) {
 		qargs->n_ref[i] += wargs->n_ref[i];
@@ -295,8 +311,7 @@ static void finalise_job(void *vqargs, void *vwargs)
 	qargs->n_done++;
 	progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
 
-	reflist_free(wargs->image.reflections);
-	cell_free(wargs->image.indexed_cell);
+	crystal_free(wargs->crystal);
 	free(wargs);
 }
 
@@ -315,7 +330,7 @@ int main(int argc, char *argv[])
 	char *sym_str = NULL;
 	SymOpList *sym;
 	UnitCell *cell = NULL;
-	FILE *ofh;
+	Stream *stream;
 	int n = 2;
 	int random_intensities = 0;
 	char *save_file = NULL;
@@ -496,13 +511,12 @@ int main(int argc, char *argv[])
 		ERROR("You must give a filename for the output.\n");
 		return 1;
 	}
-	ofh = fopen(output_file, "w");
-	if ( ofh == NULL ) {
+	stream = open_stream_for_write(output_file);
+	if ( stream == NULL ) {
 		ERROR("Couldn't open output file '%s'\n", output_file);
 		return 1;
 	}
 	free(output_file);
-	write_stream_header(ofh, argc, argv);
 
 	image.det = det;
 	image.width = det->max_fs;
@@ -511,9 +525,12 @@ int main(int argc, char *argv[])
 	image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
 	image.div = beam->divergence;
 	image.bw = beam->bandwidth;
-	image.profile_radius = beam->profile_radius;
+	image.beam = beam;
 	image.filename = malloc(256);
 	image.copyme = NULL;
+	image.crystals = NULL;
+	image.n_crystals = 0;
+	image.indexed_by = INDEXING_NONE;
 
 	if ( random_intensities ) {
 		full = reflist_new();
@@ -528,7 +545,7 @@ int main(int argc, char *argv[])
 	qargs.random_intensities = random_intensities;
 	qargs.cell = cell;
 	qargs.template_image = &image;
-	qargs.stream = ofh;
+	qargs.stream = stream;
 	qargs.cnoise = cnoise;
 	qargs.max_q = largest_q(&image);
 
@@ -574,7 +591,7 @@ int main(int argc, char *argv[])
 
 	}
 
-	fclose(ofh);
+	close_stream(stream);
 	cell_free(cell);
 	free_detector_geometry(det);
 	free(beam);
diff --git a/src/partialator.c b/src/partialator.c
index c2f6c299..a4af3c18 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
  *
  * Scaling and post refinement for coherent nanocrystallography
  *
- * 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.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -87,16 +87,16 @@ static void show_help(const char *s)
 struct refine_args
 {
 	RefList *full;
-	struct image *image;
+	Crystal *crystal;
 };
 
 
 struct queue_args
 {
-	int n;
+	int n_started;
 	int n_done;
-	int n_total_patterns;
-	struct image *images;
+	Crystal **crystals;
+	int n_crystals;
 	struct refine_args task_defaults;
 };
 
@@ -104,10 +104,9 @@ struct queue_args
 static void refine_image(void *task, int id)
 {
 	struct refine_args *pargs = task;
-	struct image *image = pargs->image;
-	image->id = id;
+	Crystal *cr = pargs->crystal;
 
-	pr_refine(image, pargs->full);
+	pr_refine(cr, pargs->full);
 }
 
 
@@ -119,9 +118,9 @@ static void *get_image(void *vqargs)
 	task = malloc(sizeof(struct refine_args));
 	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
 
-	task->image = &qargs->images[qargs->n];
+	task->crystal = qargs->crystals[qargs->n_started];
 
-	qargs->n++;
+	qargs->n_started++;
 
 	return task;
 }
@@ -133,12 +132,12 @@ static void done_image(void *vqargs, void *task)
 
 	qargs->n_done++;
 
-	progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
+	progress_bar(qargs->n_done, qargs->n_crystals, "Refining");
 	free(task);
 }
 
 
-static void refine_all(struct image *images, int n_total_patterns,
+static void refine_all(Crystal **crystals, int n_crystals,
                        struct detector *det,
                        RefList *full, int nthreads)
 {
@@ -146,19 +145,19 @@ static void refine_all(struct image *images, int n_total_patterns,
 	struct queue_args qargs;
 
 	task_defaults.full = full;
-	task_defaults.image = NULL;
+	task_defaults.crystal = NULL;
 
 	qargs.task_defaults = task_defaults;
-	qargs.n = 0;
+	qargs.n_started = 0;
 	qargs.n_done = 0;
-	qargs.n_total_patterns = n_total_patterns;
-	qargs.images = images;
+	qargs.n_crystals = n_crystals;
+	qargs.crystals = crystals;
 
 	/* Don't have threads which are doing nothing */
-	if ( n_total_patterns < nthreads ) nthreads = n_total_patterns;
+	if ( n_crystals < nthreads ) nthreads = n_crystals;
 
 	run_threads(nthreads, refine_image, get_image, done_image,
-	            &qargs, n_total_patterns, 0, 0, 0);
+	            &qargs, n_crystals, 0, 0, 0);
 }
 
 
@@ -201,13 +200,14 @@ static int select_scalable_reflections(RefList *list, RefList *reference)
 }
 
 
-static void select_reflections_for_refinement(struct image *images, int n,
+static void select_reflections_for_refinement(Crystal **crystals, int n,
                                               RefList *full, int have_reference)
 {
 	int i;
 
 	for ( i=0; i<n; i++ ) {
 
+		RefList *reflist;
 		Reflection *refl;
 		RefListIterator *iter;
 		int n_acc = 0;
@@ -216,9 +216,10 @@ static void select_reflections_for_refinement(struct image *images, int n,
 		int n_fewmatch = 0;
 		int n_ref = 0;
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
 
-		for ( refl = first_refl(images[i].reflections, &iter);
+		reflist = crystal_get_reflections(crystals[i]);
+		for ( refl = first_refl(reflist, &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
@@ -272,6 +273,20 @@ static void select_reflections_for_refinement(struct image *images, int n,
 }
 
 
+static void display_progress(int n_images, int n_crystals)
+{
+	if ( !isatty(STDERR_FILENO) ) return;
+	if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+	pthread_mutex_lock(&stderr_lock);
+	fprintf(stderr, "\r%i images loaded, %i crystals.",
+	        n_images, n_crystals);
+	pthread_mutex_unlock(&stderr_lock);
+
+	fflush(stdout);
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -280,16 +295,15 @@ int main(int argc, char *argv[])
 	char *geomfile = NULL;
 	char *sym_str = NULL;
 	SymOpList *sym;
-	FILE *fh;
 	int nthreads = 1;
 	struct detector *det;
 	int i;
-	int n_total_patterns;
 	struct image *images;
 	int n_iter = 10;
 	struct beam_params *beam = NULL;
 	RefList *full;
-	int n_usable_patterns = 0;
+	int n_images = 0;
+	int n_crystals = 0;
 	int nobs;
 	char *reference_file = NULL;
 	RefList *reference = NULL;
@@ -298,6 +312,8 @@ int main(int argc, char *argv[])
 	char cmdline[1024];
 	SRContext *sr;
 	int noscale = 0;
+	Stream *st;
+	Crystal **crystals;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -386,19 +402,15 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	/* Sanitise input filename and open */
 	if ( infile == NULL ) {
 		infile = strdup("-");
 	}
-	if ( strcmp(infile, "-") == 0 ) {
-		fh = stdin;
-	} else {
-		fh = fopen(infile, "r");
-	}
-	if ( fh == NULL ) {
-		ERROR("Failed to open input file '%s'\n", infile);
+	st = open_stream_for_read(infile);
+	if ( st == NULL ) {
+		ERROR("Failed to open input stream '%s'\n", infile);
 		return 1;
 	}
+	/* Don't free "infile", because it's needed for the scaling report */
 
 	/* Sanitise output filename */
 	if ( outfile == NULL ) {
@@ -438,86 +450,116 @@ int main(int argc, char *argv[])
 
 	}
 
-	n_total_patterns = count_patterns(fh);
-	if ( n_total_patterns == 0 ) {
-		ERROR("No patterns to process.\n");
-		return 1;
-	}
-	STATUS("There are %i patterns to process\n", n_total_patterns);
-
 	gsl_set_error_handler_off();
 
-	images = malloc(n_total_patterns * sizeof(struct image));
-	if ( images == NULL ) {
-		ERROR("Couldn't allocate memory for images.\n");
-		return 1;
-	}
-
 	/* Fill in what we know about the images so far */
-	rewind(fh);
-	nobs = 0;
-	for ( i=0; i<n_total_patterns; i++ ) {
+	n_images = 0;
+	n_crystals = 0;
+	images = NULL;
+	crystals = NULL;
+	do {
 
 		RefList *as;
+		int i;
+		struct image *images_new;
 		struct image *cur;
 
-		cur = &images[n_usable_patterns];
+		images_new = realloc(images, (n_images+1)*sizeof(struct image));
+		if ( images_new == NULL ) {
+			ERROR("Failed to allocate memory for image list.\n");
+			return 1;
+		}
+		images = images_new;
+		cur = &images[n_images];
 
 		cur->det = det;
-
-		if ( read_chunk(fh, cur) != 0 ) {
-			/* Should not happen, because we counted the patterns
-			 * earlier. */
-			ERROR("Failed to read chunk from the input stream.\n");
-			return 1;
+		if ( read_chunk(st, cur) != 0 ) {
+			break;
 		}
 
 		/* Won't be needing this, if it exists */
 		image_feature_list_free(cur->features);
 		cur->features = NULL;
-
-		/* "n_usable_patterns" will not be incremented in this case */
-		if ( cur->indexed_cell == NULL ) continue;
-
-		/* Fill in initial estimates of stuff */
 		cur->div = beam->divergence;
 		cur->bw = beam->bandwidth;
 		cur->width = det->max_fs;
 		cur->height = det->max_ss;
-		cur->osf = 1.0;
-		cur->profile_radius = beam->profile_radius;
-		cur->pr_dud = 0;
-
-		/* Muppet proofing */
 		cur->data = NULL;
 		cur->flags = NULL;
 		cur->beam = NULL;
 
-		/* This is the raw list of reflections */
-		as = asymmetric_indices(cur->reflections, sym);
-		reflist_free(cur->reflections);
-		cur->reflections = as;
+		n_images++;
+
+		for ( i=0; i<cur->n_crystals; i++ ) {
+
+			Crystal *cr;
+			Crystal **crystals_new;
+			RefList *cr_refl;
+
+			crystals_new = realloc(crystals,
+			                      (n_crystals+1)*sizeof(Crystal *));
+			if ( crystals_new == NULL ) {
+				ERROR("Failed to allocate memory for crystal "
+				      "list.\n");
+				return 1;
+			}
+			crystals = crystals_new;
+			crystals[n_crystals] = cur->crystals[i];
+			cr = crystals[n_crystals];
+
+			/* Image pointer will change due to later reallocs */
+			crystal_set_image(cr, NULL);
+
+			/* Fill in initial estimates of stuff */
+			crystal_set_osf(cr, 1.0);
+			crystal_set_profile_radius(cr, beam->profile_radius);
+			crystal_set_user_flag(cr, 0);
+
+			/* This is the raw list of reflections */
+			cr_refl = crystal_get_reflections(cr);
+			as = asymmetric_indices(cr_refl, sym);
+			crystal_set_reflections(cr, as);
+			reflist_free(cr_refl);
 
-		update_partialities(cur);
+			n_crystals++;
 
-		nobs += select_scalable_reflections(cur->reflections,
-		                                    reference);
+		}
+
+		display_progress(n_images, n_crystals);
+
+	} while ( 1 );
+
+	close_stream(st);
+
+	/* Fill in image pointers */
+	nobs = 0;
+	for ( i=0; i<n_images; i++ ) {
+		int j;
+		for ( j=0; j<images[i].n_crystals; j++ ) {
+
+			Crystal *cryst;
+			RefList *as;
 
-		progress_bar(i, n_total_patterns-1, "Loading pattern data");
-		n_usable_patterns++;
+			cryst = images[i].crystals[j];
+			crystal_set_image(cryst, &images[i]);
 
+			/* Now it's safe to do the following */
+			update_partialities(cryst);
+			as = crystal_get_reflections(cryst);
+			nobs += select_scalable_reflections(as, reference);
+
+		}
 	}
-	fclose(fh);
 
 	/* Make initial estimates */
-	STATUS("Performing initial scaling.\n");
+	STATUS("\nPerforming initial scaling.\n");
 	if ( noscale ) STATUS("Scale factors fixed at 1.\n");
-	full = scale_intensities(images, n_usable_patterns, reference,
+	full = scale_intensities(crystals, n_crystals, reference,
 	                         nthreads, noscale);
 
-	sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf",
+	sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf",
 	                  infile, cmdline);
-	sr_iteration(sr, 0, images, n_usable_patterns, full);
+	sr_iteration(sr, 0, crystals, n_crystals, full);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -534,54 +576,55 @@ int main(int argc, char *argv[])
 		}
 
 		/* Refine the geometry of all patterns to get the best fit */
-		select_reflections_for_refinement(images, n_usable_patterns,
+		select_reflections_for_refinement(crystals, n_crystals,
 		                                  comp, have_reference);
-		refine_all(images, n_usable_patterns, det, comp, nthreads);
+		refine_all(crystals, n_crystals, det, comp, nthreads);
 
 		nobs = 0;
-		for ( j=0; j<n_usable_patterns; j++ ) {
+		for ( j=0; j<n_crystals; j++ ) {
 
-			struct image *cur = &images[j];
+			Crystal *cr = crystals[j];
+			RefList *rf = crystal_get_reflections(cr);
 
-			nobs += select_scalable_reflections(cur->reflections,
-			                                    reference);
+			nobs += select_scalable_reflections(rf, reference);
 
 		}
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
-		full = scale_intensities(images, n_usable_patterns,
+		full = scale_intensities(crystals, n_crystals,
 		                         reference, nthreads, noscale);
 
-		sr_iteration(sr, i+1, images, n_usable_patterns, full);
+		sr_iteration(sr, i+1, crystals, n_crystals, full);
 
 	}
 
 	sr_finish(sr);
 
 	n_dud = 0;
-	for ( i=0; i<n_usable_patterns; i++ ) {
-		if ( images[i].pr_dud ) n_dud++;
+	for ( i=0; i<n_crystals; i++ ) {
+		if ( crystal_get_user_flag(crystals[i]) ) n_dud++;
 	}
-	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
+	STATUS("%i crystals could not be refined on the last cycle.\n", n_dud);
 
 	/* Output results */
 	write_reflist(outfile, full);
 
 	/* Clean up */
-	for ( i=0; i<n_usable_patterns; i++ ) {
-		reflist_free(images[i].reflections);
+	for ( i=0; i<n_crystals; i++ ) {
+		reflist_free(crystal_get_reflections(crystals[i]));
+		crystal_free(crystals[i]);
 	}
 	reflist_free(full);
 	free(sym);
 	free(outfile);
 	free_detector_geometry(det);
 	free(beam);
+	free(crystals);
 	if ( reference != NULL ) {
 		reflist_free(reference);
 	}
-	for ( i=0; i<n_usable_patterns; i++ ) {
-		cell_free(images[i].indexed_cell);
+	for ( i=0; i<n_images; i++ ) {
 		free(images[i].filename);
 	}
 	free(images);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 7410d931..1439b148 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * 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.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius)
 
 
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, Reflection *refl, double r)
+double gradient(Crystal *cr, int k, Reflection *refl)
 {
 	double ds, azix, aziy;
 	double ttlow, tthigh, tt;
@@ -103,17 +103,19 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
 	int clamp_low, clamp_high;
 	double klow, khigh;
 	double gr;
+	struct image *image = crystal_get_image(cr);
+	double r = crystal_get_profile_radius(cr);
 
 	get_symmetric_indices(refl, &hs, &ks, &ls);
 
-	cell_get_reciprocal(image->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);
 	xl = hs*asx + ks*bsx + ls*csx;
 	yl = hs*asy + ks*bsy + ls*csy;
 	zl = hs*asz + ks*bsz + ls*csz;
 
-	ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
+	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
 	get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
 
 	klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
@@ -237,8 +239,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
 
 
 /* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(struct image *image, int k, double shift)
+static void apply_shift(Crystal *cr, int k, double shift)
 {
+	double t;
+	struct image *image = crystal_get_image(cr);
+
 	switch ( k ) {
 
 		case REF_DIV :
@@ -251,7 +256,9 @@ static void apply_shift(struct image *image, int k, double shift)
 		break;
 
 		case REF_R :
-		image->profile_radius += shift;
+		t = crystal_get_profile_radius(cr);
+		t += shift;
+		crystal_set_profile_radius(cr, t);
 		break;
 
 		case REF_ASX :
@@ -263,7 +270,7 @@ static void apply_shift(struct image *image, int k, double shift)
 		case REF_CSX :
 		case REF_CSY :
 		case REF_CSZ :
-		apply_cell_shift(image->indexed_cell, k, shift);
+		apply_cell_shift(crystal_get_cell(cr), k, shift);
 		break;
 
 		default :
@@ -357,7 +364,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
 
 
 /* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(struct image *image, const RefList *full)
+static double pr_iterate(Crystal *cr, const RefList *full)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -369,7 +376,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 	double max_shift;
 	int nref = 0;
 
-	reflections = image->reflections;
+	reflections = crystal_get_reflections(cr);
 
 	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
 	v = gsl_vector_calloc(NUM_PARAMS);
@@ -387,6 +394,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 		double p;
 		Reflection *match;
 		double gradients[NUM_PARAMS];
+		const double osf = crystal_get_osf(cr);
 
 		if ( !get_refinable(refl) ) continue;
 
@@ -398,7 +406,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 			      " is it marked as refinable?\n", ha, ka, la);
 			continue;
 		}
-		I_full = image->osf * get_intensity(match);
+		I_full = osf * get_intensity(match);
 
 		/* Actual measurement of this reflection from this pattern? */
 		I_partial = get_intensity(refl);
@@ -412,7 +420,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 			double gr;
-			gr = gradient(image, k, refl, image->profile_radius);
+			gr = gradient(cr, k, refl);
 			gradients[k] = gr;
 		}
 
@@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 				if ( g > k ) continue;
 
 				M_c = gradients[g] * gradients[k];
-				M_c *= w * pow(image->osf * I_full, 2.0);
+				M_c *= w * pow(osf * I_full, 2.0);
 
 				M_curr = gsl_matrix_get(M, k, g);
 				gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 
 	//STATUS("%i reflections went into the equations.\n", nref);
 	if ( nref == 0 ) {
-		image->pr_dud = 1;
+		crystal_set_user_flag(cr, 1);
 		gsl_matrix_free(M);
 		gsl_vector_free(v);
 		return 0.0;
@@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 
 		for ( param=0; param<NUM_PARAMS; param++ ) {
 			double shift = gsl_vector_get(shifts, param);
-			apply_shift(image, param, shift);
+			apply_shift(cr, param, shift);
 			//STATUS("Shift %i: %e\n", param, shift);
 			if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
 		}
@@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full)
 }
 
 
-static double guide_dev(struct image *image, const RefList *full)
+static double guide_dev(Crystal *cr, const RefList *full)
 {
 	double dev = 0.0;
 
@@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full)
 	Reflection *refl;
 	RefListIterator *iter;
 
-	for ( refl = first_refl(image->reflections, &iter);
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
 
@@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full)
 		 * scale_intensities() might not yet have been called, so the
 		 * full version may not have been calculated yet. */
 
-		G = image->osf;
+		G = crystal_get_osf(cr);
 		p = get_partiality(refl);
 		I_partial = get_intensity(refl);
 		I_full = get_intensity(full_version);
@@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full)
 }
 
 
-void pr_refine(struct image *image, const RefList *full)
+void pr_refine(Crystal *cr, const RefList *full)
 {
 	double max_shift, dev;
 	int i;
 	const int verbose = 0;
 
 	if ( verbose ) {
-		dev = guide_dev(image, full);
+		dev = guide_dev(cr, full);
 		STATUS("\n");  /* Deal with progress bar */
 		STATUS("Before iteration:                       dev = %10.5e\n",
 		       dev);
 	}
 
 	i = 0;
-	image->pr_dud = 0;
+	crystal_set_user_flag(cr, 0);
 	do {
 
 		double asx, asy, asz;
@@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full)
 		double csx, csy, csz;
 		double dev;
 
-		cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+		cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
-		max_shift = pr_iterate(image, full);
+		max_shift = pr_iterate(cr, full);
 
-		update_partialities(image);
+		update_partialities(cr);
 
 		if ( verbose ) {
-			dev = guide_dev(image, full);
+			dev = guide_dev(cr, full);
 			STATUS("PR Iteration %2i: max shift = %10.2f"
 			       " dev = %10.5e\n", i+1, max_shift, dev);
 		}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 2223dcdf..fe171882 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -39,6 +39,7 @@
 
 #include "image.h"
 #include "utils.h"
+#include "crystal.h"
 
 
 /* Refineable parameters */
@@ -58,10 +59,10 @@ enum {
 };
 
 
-extern void pr_refine(struct image *image, const RefList *full);
+extern void pr_refine(Crystal *cr, const RefList *full);
 
 /* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(struct image *image, int k, Reflection *refl, double r);
+extern double gradient(Crystal *cr, int k, Reflection *refl);
 
 
 #endif	/* POST_REFINEMENT_H */
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);
diff --git a/src/process_image.c b/src/process_image.c
new file mode 100644
index 00000000..010054c8
--- /dev/null
+++ b/src/process_image.c
@@ -0,0 +1,231 @@
+/*
+ * process_image.c
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <hdf5.h>
+#include <gsl/gsl_errno.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "index.h"
+#include "peaks.h"
+#include "detector.h"
+#include "filters.h"
+#include "thread-pool.h"
+#include "beam-parameters.h"
+#include "geometry.h"
+#include "stream.h"
+#include "reflist-utils.h"
+#include "process_image.h"
+
+
+void process_image(const struct index_args *iargs, struct pattern_args *pargs,
+                   Stream *st, int cookie)
+{
+	float *data_for_measurement;
+	size_t data_size;
+	int check;
+	struct hdfile *hdfile;
+	struct image image;
+	int i;
+	char filename[1024];
+
+	/* Prefix to jump out of temporary folder */
+	if ( pargs->filename[0] != '/' ) {
+		snprintf(filename, 1023, "../../%s", pargs->filename);
+	} else {
+		snprintf(filename, 1023, "%s", pargs->filename);
+	}
+
+	image.features = NULL;
+	image.data = NULL;
+	image.flags = NULL;
+	image.copyme = iargs->copyme;
+	image.id = cookie;
+	image.filename = pargs->filename;  /* Relative to top level */
+	image.beam = iargs->beam;
+	image.det = iargs->det;
+	image.crystals = NULL;
+	image.n_crystals = 0;
+
+	hdfile = hdfile_open(filename);  /* Relative to temporary folder */
+	if ( hdfile == NULL ) return;
+
+	if ( iargs->element != NULL ) {
+
+		int r;
+		r = hdfile_set_image(hdfile, iargs->element);
+		if ( r ) {
+			ERROR("Couldn't select path '%s'\n", iargs->element);
+			hdfile_close(hdfile);
+			return;
+		}
+
+	} else {
+
+		int r;
+		r = hdfile_set_first_image(hdfile, "/");
+		if ( r ) {
+			ERROR("Couldn't select first path\n");
+			hdfile_close(hdfile);
+			return;
+		}
+
+	}
+
+	check = hdf5_read(hdfile, &image, 1);
+	if ( check ) {
+		hdfile_close(hdfile);
+		return;
+	}
+
+	if ( (image.width != image.det->max_fs + 1 )
+	  || (image.height != image.det->max_ss + 1))
+	{
+		ERROR("Image size doesn't match geometry size"
+			" - rejecting image.\n");
+		ERROR("Image size: %i,%i.  Geometry size: %i,%i\n",
+		      image.width, image.height,
+		      image.det->max_fs + 1, image.det->max_ss + 1);
+		hdfile_close(hdfile);
+		return;
+	}
+
+	fill_in_values(image.det, hdfile);
+	fill_in_beam_parameters(image.beam, hdfile);
+
+	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+
+	if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
+		/* Error message covers a silly value in the beam file or in
+		 * the HDF5 file. */
+		ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
+		      image.lambda, image.beam->photon_energy, image.filename);
+		hdfile_close(hdfile);
+		return;
+	}
+
+	if ( iargs->cmfilter ) filter_cm(&image);
+
+	/* Take snapshot of image after CM subtraction but before
+	 * the aggressive noise filter. */
+	data_size = image.width * image.height * sizeof(float);
+	data_for_measurement = malloc(data_size);
+
+	if ( iargs->noisefilter ) {
+		filter_noise(&image, data_for_measurement);
+	} else {
+		memcpy(data_for_measurement, image.data, data_size);
+	}
+
+	switch ( iargs->peaks ) {
+
+		case PEAK_HDF5:
+		// Get peaks from HDF5
+		if (get_peaks(&image, hdfile,
+			iargs->hdf5_peak_path)) {
+			ERROR("Failed to get peaks from HDF5 file.\n");
+		}
+		if ( !iargs->no_revalidate ) {
+			validate_peaks(&image, iargs->min_int_snr,
+				       iargs->ir_inn, iargs->ir_mid,
+				       iargs->ir_out, iargs->use_saturated);
+		}
+		break;
+
+		case PEAK_ZAEF:
+		search_peaks(&image, iargs->threshold,
+		             iargs->min_gradient, iargs->min_snr,
+		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+		             iargs->use_saturated);
+		break;
+
+	}
+
+	/* Get rid of noise-filtered version at this point
+	 * - it was strictly for the purposes of peak detection. */
+	free(image.data);
+	image.data = data_for_measurement;
+
+	/* Index the pattern */
+	index_pattern(&image, iargs->indm, iargs->ipriv);
+
+	pargs->n_crystals = image.n_crystals;
+
+	/* Default beam parameters */
+	image.div = image.beam->divergence;
+	image.bw = image.beam->bandwidth;
+
+	/* Integrate each crystal's diffraction spots */
+	for ( i=0; i<image.n_crystals; i++ ) {
+
+		RefList *reflections;
+
+		/* Set default crystal parameter(s) */
+		crystal_set_profile_radius(image.crystals[i],
+		                           image.beam->profile_radius);
+
+		if ( iargs->integrate_found ) {
+			reflections = select_intersections(&image,
+			                                   image.crystals[i]);
+		} else {
+			reflections = find_intersections(&image,
+			                                 image.crystals[i]);
+		}
+
+		crystal_set_reflections(image.crystals[i], reflections);
+
+	}
+
+	/* Integrate all the crystals at once - need all the crystals so that
+	 * overlaps can be detected. */
+	integrate_reflections(&image, iargs->closer,
+	                              iargs->bgsub,
+	                              iargs->min_int_snr,
+	                              iargs->ir_inn,
+	                              iargs->ir_mid,
+	                              iargs->ir_out,
+	                              iargs->integrate_saturated,
+	                              iargs->res_cutoff);
+
+	write_chunk(st, &image, hdfile,
+	            iargs->stream_peaks, iargs->stream_refls);
+
+	for ( i=0; i<image.n_crystals; i++ ) {
+		crystal_free(image.crystals[i]);
+	}
+
+	free(image.data);
+	if ( image.flags != NULL ) free(image.flags);
+	image_feature_list_free(image.features);
+	hdfile_close(hdfile);
+}
diff --git a/src/process_image.h b/src/process_image.h
new file mode 100644
index 00000000..f2f034d9
--- /dev/null
+++ b/src/process_image.h
@@ -0,0 +1,94 @@
+/*
+ * process_image.h
+ *
+ * The processing pipeline for one image
+ *
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2010-2013 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef PROCESS_IMAGE_H
+#define PROCESS_IMAGE_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+enum {
+	PEAK_ZAEF,
+	PEAK_HDF5,
+};
+
+
+/* Information about the indexing process which is common to all patterns */
+struct index_args
+{
+	UnitCell *cell;
+	int cmfilter;
+	int noisefilter;
+	int satcorr;
+	int closer;
+	int bgsub;
+	float threshold;
+	float min_gradient;
+	float min_snr;
+	double min_int_snr;
+	struct detector *det;
+	IndexingMethod *indm;
+	IndexingPrivate **ipriv;
+	int peaks;                /* Peak detection method */
+	float tols[4];
+	struct beam_params *beam;
+	char *element;
+	char *hdf5_peak_path;
+	float ir_inn;
+	float ir_mid;
+	float ir_out;
+	struct copy_hdf5_field *copyme;
+	int integrate_saturated;
+	int use_saturated;
+	int no_revalidate;
+	int integrate_found;
+	int stream_peaks;
+	int stream_refls;
+	int res_cutoff;
+};
+
+
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+	/* "Input" */
+	char *filename;
+
+	/* "Output" */
+	int n_crystals;
+};
+
+
+extern void process_image(const struct index_args *iargs,
+                          struct pattern_args *pargs, Stream *st,
+                          int cookie);
+
+
+#endif	/* PROCESS_IMAGEs_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index ba425b96..a02a1796 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -3,11 +3,11 @@
  *
  * Write a nice PDF of scaling parameters
  *
- * 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.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height,
 }
 
 
-static void partiality_graph(cairo_t *cr, const struct image *images, int n,
+static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
                              RefList *full)
 {
 	const double g_width = 200.0;
@@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 
 	num_nondud = 0;
 	for ( i=0; i<n; i++ ) {
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		num_nondud++;
 	}
 
@@ -229,10 +229,13 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 
 		Reflection *refl;
 		RefListIterator *iter;
+		Crystal *cryst;
 
-		if ( images[i].pr_dud ) continue;
+		cryst = crystals[i];
 
-		for ( refl = first_refl(images[i].reflections, &iter);
+		if ( crystal_get_user_flag(cryst) ) continue;
+
+		for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
@@ -249,7 +252,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 			if ( get_redundancy(f) < 2 ) continue;
 
 			Ipart = get_intensity(refl);
-			Ifull = images[i].osf * get_intensity(f);
+			Ifull = crystal_get_osf(cryst) * get_intensity(f);
 
 			//if ( Ifull < 10 ) continue;  /* FIXME: Ugh */
 
@@ -267,7 +270,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 					double esd_pobs, esd_Ip, esd_If;
 					esd_Ip = get_esd_intensity(refl);
 					esd_If = get_esd_intensity(f);
-					esd_If *= images[i].osf;
+					esd_If *= crystal_get_osf(cryst);
 					esd_pobs  = pow(esd_Ip/Ipart, 2.0);
 					esd_pobs += pow(esd_If/Ifull, 2.0);
 					esd_pobs = sqrt(esd_pobs);
@@ -313,8 +316,8 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n,
 }
 
 
-static void partiality_histogram(cairo_t *cr, const struct image *images,
-                                 int n, RefList *full, int calc, int backwards)
+static void partiality_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                 RefList *full, int calc, int backwards)
 {
 	int f_max;
 	int i, b;
@@ -344,10 +347,11 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 
 		Reflection *refl;
 		RefListIterator *iter;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		for ( refl = first_refl(images[i].reflections, &iter);
+		for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
 		      refl != NULL;
 		      refl = next_refl(refl, iter) )
 		{
@@ -364,7 +368,7 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 			Ipart = get_intensity(refl);
 			Ifull = get_intensity(f);
 
-			pobs = Ipart/(images[i].osf*Ifull);
+			pobs = Ipart/(crystal_get_osf(cryst)*Ifull);
 			pcalc = get_partiality(refl);
 
 			if ( calc ) {
@@ -426,8 +430,8 @@ static void partiality_histogram(cairo_t *cr, const struct image *images,
 }
 
 
-static void scale_factor_histogram(cairo_t *cr, const struct image *images,
-                                   int n, const char *title)
+static void scale_factor_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                   const char *title)
 {
 	int f_max;
 	int i, b;
@@ -451,8 +455,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 
 	osf_max = 0.0;
 	for ( i=0; i<n; i++ ) {
-		double osf = images[i].osf;
-		if ( images[i].pr_dud ) continue;
+		double osf = crystal_get_osf(crystals[i]);
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		if ( osf > osf_max ) osf_max = osf;
 	}
 	osf_max = ceil(osf_max+osf_max/10000.0);
@@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 
 		for ( i=0; i<n; i++ ) {
 
-			double osf = images[i].osf;
+			double osf = crystal_get_osf(crystals[i]);
 
-			if ( images[i].pr_dud ) continue;
+			if ( crystal_get_user_flag(crystals[i]) ) continue;
 
 			for ( b=0; b<nbins; b++ ) {
 				if ( (osf >= osf_low[b])
@@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images,
 }
 
 
-static void intensity_histogram(cairo_t *cr, const struct image *images,
-                                int n, RefList *full,
+static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n,
+                                RefList *full,
                                 signed int h, signed int k, signed int l)
 {
 	int f_max;
@@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
 	for ( i=0; i<n; i++ ) {
 
 		double osf;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		osf = images[i].osf;
+		osf = crystal_get_osf(cryst);
 
-		for ( f = find_refl(images[i].reflections, h, k, l);
+		for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
 		      f != NULL;
 		      f = next_found_refl(f) )
 		{
@@ -616,12 +621,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images,
 	for ( i=0; i<n; i++ ) {
 
 		double osf;
+		Crystal *cryst = crystals[i];
 
-		if ( images[i].pr_dud ) continue;
+		if ( crystal_get_user_flag(cryst) ) continue;
 
-		osf = images[i].osf;
+		osf = crystal_get_osf(cryst);
 
-		for ( f = find_refl(images[i].reflections, h, k, l);
+		for ( f = find_refl(crystal_get_reflections(cryst), h, k, l);
 		      f != NULL;
 		      f = next_found_refl(f) )
 		{
@@ -727,6 +733,13 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h,
 	Reflection *refl;
 	RefListIterator *iter;
 	int *samples;
+	int j;
+
+	for ( j=0; j<n; j++ ) {
+		h[j] = 0;
+		k[j] = 0;
+		l[j] = 0;
+	}
 
 	samples = calloc(n, sizeof(int));
 
@@ -771,7 +784,7 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h,
 
 
 
-SRContext *sr_titlepage(struct image *images, int n,
+SRContext *sr_titlepage(Crystal **crystals, int n,
                         const char *filename, const char *stream_filename,
                         const char *cmdline)
 {
@@ -805,7 +818,7 @@ SRContext *sr_titlepage(struct image *images, int n,
 }
 
 
-void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
+void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n,
                   RefList *full)
 {
 	int i;
@@ -822,7 +835,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 
 	cairo_save(sr->cr);
 	cairo_translate(sr->cr, 480.0, 350.0);
-	scale_factor_histogram(sr->cr, images, n,
+	scale_factor_histogram(sr->cr, crystals, n,
 	                       "Distribution of overall scale factors");
 	cairo_restore(sr->cr);
 
@@ -830,7 +843,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_save(sr->cr);
 
 	cairo_translate(sr->cr, 70.0, 330.0);
-	partiality_graph(sr->cr, images, n, full);
+	partiality_graph(sr->cr, crystals, n, full);
 
 	cairo_save(sr->cr);
 	cairo_move_to(sr->cr, 0.0, 0.0);
@@ -841,7 +854,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_stroke(sr->cr);
 	cairo_set_dash(sr->cr, NULL, 0, 0.0);
 	cairo_translate(sr->cr, 0.0, -150.0);
-	partiality_histogram(sr->cr, images, n, full, 1, 0);
+	partiality_histogram(sr->cr, crystals, n, full, 1, 0);
 	cairo_restore(sr->cr);
 
 	cairo_save(sr->cr);
@@ -854,7 +867,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 	cairo_set_dash(sr->cr, NULL, 0, 0.0);
 	cairo_translate(sr->cr, 230.0, 200.0);
 	cairo_rotate(sr->cr, -M_PI_2);
-	partiality_histogram(sr->cr, images, n, full, 0, 1);
+	partiality_histogram(sr->cr, crystals, n, full, 0, 1);
 	cairo_restore(sr->cr);
 
 	cairo_restore(sr->cr);
@@ -873,7 +886,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
 
 		cairo_save(sr->cr);
 		cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y);
-		intensity_histogram(sr->cr, images, n, full,
+		intensity_histogram(sr->cr, crystals, n, full,
 		                    sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]);
 		cairo_restore(sr->cr);
 
diff --git a/src/scaling-report.h b/src/scaling-report.h
index b9ac3fb7..5b153377 100644
--- a/src/scaling-report.h
+++ b/src/scaling-report.h
@@ -40,25 +40,25 @@ typedef struct _srcontext SRContext;  /* Opaque */
 
 #ifdef HAVE_CAIRO
 
-extern SRContext *sr_titlepage(struct image *images, int n,
+extern SRContext *sr_titlepage(Crystal **crystals, int n,
                                const char *filename,
                                const char *stream_filename,
                                const char *cmdline);
 
-extern void sr_iteration(SRContext *sr, int iteration, struct image *images,
-                         int n, RefList *full);
+extern void sr_iteration(SRContext *sr, int iteration,
+                         Crystal **crystals, int n, RefList *full);
 
 extern void sr_finish(SRContext *sr);
 
 #else
 
-SRContext *sr_titlepage(struct image *images, int n, const char *filename,
+SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename,
                         const char *stream_filename, const char *cmdline)
 {
 	return NULL;
 }
 
-void sr_iteration(SRContext *sr, int iteration, struct image *images, int n,
+void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n,
                   RefList *full)
 {
 }