First part of stream rework

12 files changed, 177 insertions, 199 deletions

diff --git a/Makefile.am b/Makefile.am
index 5bddb3a8..bdd3f2be 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -44,7 +44,8 @@ endif
 src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
                           src/cell.c src/utils.c src/reflections.c \
                           src/symmetry.c src/stream.c src/beam-parameters.c \
-                          src/thread-pool.c
+                          src/thread-pool.c src/image.c src/detector.c \
+                          src/hdf5-file.c
 
 src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
                           src/cell.c src/image.c src/peaks.c src/index.c \
@@ -52,7 +53,7 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \
                           src/sfac.c src/dirax.c src/mosflm.c \
                           src/reflections.c src/symmetry.c \
                           src/geometry.c src/thread-pool.c \
-                          src/beam-parameters.c src/reflist.c
+                          src/beam-parameters.c src/reflist.c src/stream.c
 if HAVE_OPENCL
 src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c
 endif
diff --git a/Makefile.in b/Makefile.in
index 603c48e8..f37e4218 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -150,7 +150,8 @@ am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \
 	src/filters.c src/diffraction.c src/detector.c src/sfac.c \
 	src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
 	src/geometry.c src/thread-pool.c src/beam-parameters.c \
-	src/reflist.c src/diffraction-gpu.c src/cl-utils.c
+	src/reflist.c src/stream.c src/diffraction-gpu.c \
+	src/cl-utils.c
 @HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \
 @HAVE_OPENCL_TRUE@	src/cl-utils.$(OBJEXT)
 am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
@@ -161,7 +162,7 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \
 	src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \
 	src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \
-	src/reflist.$(OBJEXT) $(am__objects_1)
+	src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1)
 src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS)
 src_indexamajig_LDADD = $(LDADD)
 src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -202,7 +203,9 @@ am_src_process_hkl_OBJECTS = src/process_hkl.$(OBJEXT) \
 	src/sfac.$(OBJEXT) src/statistics.$(OBJEXT) src/cell.$(OBJEXT) \
 	src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \
 	src/symmetry.$(OBJEXT) src/stream.$(OBJEXT) \
-	src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT)
+	src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT) \
+	src/image.$(OBJEXT) src/detector.$(OBJEXT) \
+	src/hdf5-file.$(OBJEXT)
 src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS)
 src_process_hkl_LDADD = $(LDADD)
 src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -618,14 +621,15 @@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \
 src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \
                           src/cell.c src/utils.c src/reflections.c \
                           src/symmetry.c src/stream.c src/beam-parameters.c \
-                          src/thread-pool.c
+                          src/thread-pool.c src/image.c src/detector.c \
+                          src/hdf5-file.c
 
 src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \
 	src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \
 	src/filters.c src/diffraction.c src/detector.c src/sfac.c \
 	src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \
 	src/geometry.c src/thread-pool.c src/beam-parameters.c \
-	src/reflist.c $(am__append_6)
+	src/reflist.c src/stream.c $(am__append_6)
 @BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \
 @BUILD_HDFSEE_TRUE@                     src/hdf5-file.c src/utils.c src/image.c src/filters.c \
 @BUILD_HDFSEE_TRUE@                     src/thread-pool.c src/detector.c
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 1b957b4c..20141522 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -3,7 +3,7 @@
  *
  * Read/write HDF5 data files
  *
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
diff --git a/src/image.h b/src/image.h
index e0c190f2..69c36de8 100644
--- a/src/image.h
+++ b/src/image.h
@@ -85,7 +85,7 @@ struct image {
 	int                     width;
 	int                     height;
 
-	/* Reflections (used for scaling ONLY) */
+	/* Integrated (or about-to-be-integrated) reflections */
 	RefList                 *reflections;
 
 	/* Detected peaks */
diff --git a/src/index.c b/src/index.c
index 554c099d..ddb25eb4 100644
--- a/src/index.c
+++ b/src/index.c
@@ -128,6 +128,8 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
 	int i;
 	int n = 0;
 
+	if ( indm == NULL ) return;
+
 	map_all_peaks(image);
 
 	while ( indm[n] != INDEXING_NONE ) {
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 51461c91..72491450 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -36,6 +36,7 @@
 #include "thread-pool.h"
 #include "beam-parameters.h"
 #include "geometry.h"
+#include "stream.h"
 
 
 enum {
@@ -50,9 +51,8 @@ struct static_index_args
 	UnitCell *cell;
 	int config_cmfilter;
 	int config_noisefilter;
-	int config_dumpfound;
 	int config_verbose;
-	int config_nearbragg;
+	int stream_flags;
 	int config_polar;
 	int config_satcorr;
 	int config_closer;
@@ -109,7 +109,8 @@ static void show_help(const char *s)
 "\n"
 " -i, --input=<filename>   Specify file containing list of images to process.\n"
 "                           '-' means stdin, which is the default.\n"
-" -o, --output=<filename>  Write indexed stream to this file. '-' for stdout.\n"
+" -o, --output=<filename>  Write output stream to this file. '-' for stdout.\n"
+"                           Default: indexamajig.stream\n"
 "\n"
 "     --indexing=<methods> Use 'methods' for indexing.  Provide one or more\n"
 "                           methods separated by commas.  Choose from:\n"
@@ -129,25 +130,19 @@ static void show_help(const char *s)
 "                                    This is the default method.\n"
 "                           hdf5  : Get from /processing/hitfinder/peakinfo\n"
 "                                    in the HDF5 file.\n"
-"\n"
-"\nWith just the above options, this program does not do much of practical use."
-"\nYou should also enable some of the following:\n\n"
-"     --near-bragg         Output a list of reflection intensities to stdout.\n"
-"                           When pixels with fractional indices within 0.1 of\n"
-"                           integer values (the Bragg condition) are found,\n"
-"                           the integral of pixels within a ten pixel radius\n"
-"                           of the nearest-to-Bragg pixel will be reported as\n"
-"                           the intensity.  The centroid of the pixels will\n"
-"                           be given as the coordinates, as well as the h,k,l\n"
-"                           (integer) indices of the reflection.  If a peak\n"
-"                           was located by the initial peak search close to\n"
-"                           the \"near Bragg\" location, its coordinates will\n"
-"                           be taken as the centre instead.\n"
-"     --dump-peaks         Write the results of the peak search to stdout.\n"
-"                           The intensities in this list are from the\n"
-"                           centroid/integration procedure.\n"
-"\n"
-"\nFor more control over the process, you might need:\n\n"
+"\n\n"
+"You can control what information is included in the output stream using\n"
+"' --record=<flags>'.  Possible flags are:\n"
+"pixels         Include a list of sums of pixel values within the\n"
+"                integration domain, correcting for individual pixel\n"
+"                solid angles.\n"
+"integrated     Include a list of reflection intensities, produced by\n"
+"                integrating around predicted peak locations.\n"
+"                The flags 'pixels' and 'integrated' are mutually exclusive.\n"
+"peaks          Include peak locations and intensities from the peak search.\n"
+"peaksifindexed Include peak locations only if the pattern could be indexed.\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"
@@ -196,9 +191,7 @@ static void process_image(void *pp, int cookie)
 	UnitCell *cell = pargs->static_args.cell;
 	int config_cmfilter = pargs->static_args.config_cmfilter;
 	int config_noisefilter = pargs->static_args.config_noisefilter;
-	int config_dumpfound = pargs->static_args.config_dumpfound;
 	int config_verbose = pargs->static_args.config_verbose;
-	int config_nearbragg = pargs->static_args.config_nearbragg;
 	int config_polar = pargs->static_args.config_polar;
 	IndexingMethod *indm = pargs->static_args.indm;
 	const struct beam_params *beam = pargs->static_args.beam;
@@ -274,7 +267,6 @@ static void process_image(void *pp, int cookie)
 		/* Get peaks from HDF5 */
 		if ( get_peaks(&image, hdfile) ) {
 			ERROR("Failed to get peaks from HDF5 file.\n");
-			return;
 		}
 		break;
 	case PEAK_ZAEF :
@@ -288,51 +280,44 @@ static void process_image(void *pp, int cookie)
 	free(image.data);
 	image.data = data_for_measurement;
 
-	if ( config_dumpfound ) {
-		dump_peaks(&image, pargs->static_args.ofh,
-		           pargs->static_args.output_mutex);
-	}
-
-	/* Not indexing? Then there's nothing left to do. */
-	if ( indm == NULL ) goto done;
-
 	/* Calculate orientation matrix (by magic) */
 	index_pattern(&image, cell, indm, pargs->static_args.cellr,
 		      config_verbose, pargs->static_args.ipriv,
 		      pargs->static_args.config_insane);
 
 	/* No cell at this point?  Then we're done. */
-	if ( image.indexed_cell == NULL ) goto done;
-	pargs->indexable = 1;
+	if ( image.indexed_cell != NULL ) pargs->indexable = 1;
 
 	/* Measure intensities if requested */
-	if ( config_nearbragg ) {
 
-		RefList *reflections;
+	/* Do EITHER: */
 
-		image.div = beam->divergence;
-		image.bw = beam->bandwidth;
-		image.profile_radius = 0.0001e9;
+	//image.div = beam->divergence;
+	//image.bw = beam->bandwidth;
+	//image.profile_radius = 0.0001e9;
+	//image.reflections = find_intersections(&image,
+	//                                       image.indexed_cell, 0);
 
-		//reflections = find_intersections(&image, image.indexed_cell,
-		//                                 0);
-		reflections = find_projected_peaks(&image, image.indexed_cell,
-		                                   0, 0.1);
+	image.reflections = find_projected_peaks(&image,
+	                                         image.indexed_cell,
+	                                         0, 0.1);
 
-		output_intensities(&image, image.indexed_cell, reflections,
-		                   pargs->static_args.output_mutex,
-		                   config_polar,
-		                   pargs->static_args.config_closer,
-		                   pargs->static_args.ofh);
+	integrate_reflections(&image, config_polar,
+	                      pargs->static_args.config_closer);
 
-		reflist_free(reflections);
+	/* OR */
 
-	}
+	//image.reflections = integrate_pixels(&image, 0, 0.1,
+	//                                     config_polar);
+
+	pthread_mutex_lock(pargs->static_args.output_mutex);
+	write_chunk(pargs->static_args.ofh, &image,
+	            pargs->static_args.stream_flags);
+	pthread_mutex_unlock(pargs->static_args.output_mutex);
 
 	/* Only free cell if found */
 	cell_free(image.indexed_cell);
 
-done:
 	free(image.data);
 	free(image.flags);
 	image_feature_list_free(image.features);
@@ -737,9 +722,7 @@ int main(int argc, char *argv[])
 	qargs.static_args.cell = cell;
 	qargs.static_args.config_cmfilter = config_cmfilter;
 	qargs.static_args.config_noisefilter = config_noisefilter;
-	qargs.static_args.config_dumpfound = config_dumpfound;
 	qargs.static_args.config_verbose = config_verbose;
-	qargs.static_args.config_nearbragg = config_nearbragg;
 	qargs.static_args.config_polar = config_polar;
 	qargs.static_args.config_satcorr = config_satcorr;
 	qargs.static_args.config_closer = config_closer;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index c823c3a3..948458bc 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -551,15 +551,23 @@ int main(int argc, char *argv[])
 
 		if ( config_nearbragg ) {
 
-			RefList *reflections;
+			/* Do EITHER: */
 
-			reflections = find_projected_peaks(&image, cell,
-			                                   0, 0.1);
+			//image.div = beam->divergence;
+			//image.bw = beam->bandwidth;
+			//image.profile_radius = 0.0001e9;
+			//image.reflections = find_intersections(&image,
+			//                               image.indexed_cell, 0);
 
-			output_intensities(&image, cell, reflections, NULL,
-			                   0, 0, stdout);
+			image.reflections = find_projected_peaks(&image,
+			                            image.indexed_cell, 0, 0.1);
 
-			reflist_free(reflections);
+			integrate_reflections(&image, 0, 0);
+
+			/* OR */
+
+			//image.reflections = integrate_pixels(&image, 0, 0.1,
+			//                                     config_polar);
 
 		}
 
diff --git a/src/peaks.c b/src/peaks.c
index 4bf56ce9..69afc20f 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -416,39 +416,6 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
 }
 
 
-void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex)
-{
-	int i;
-
-	/* Get exclusive access to the output stream if necessary */
-	if ( mutex != NULL ) pthread_mutex_lock(mutex);
-
-	fprintf(ofh, "Peaks from peak search in %s\n", image->filename);
-	fprintf(ofh, "  x/px     y/px   (1/d)/nm^-1    Intensity\n");
-
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f;
-		struct rvec r;
-		double q;
-
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-
-		r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda);
-		q = modulus(r.u, r.v, r.w);
-
-		fprintf(ofh, "%8.3f %8.3f %8.3f    %12.3f\n",
-		       f->x, f->y, q/1.0e9, f->intensity);
-
-	}
-
-	fprintf(ofh, "\n");
-
-	if ( mutex != NULL ) pthread_mutex_unlock(mutex);
-}
-
-
 RefList *find_projected_peaks(struct image *image, UnitCell *cell,
                               int circular_domain, double domain_r)
 {
@@ -600,63 +567,13 @@ int peak_sanity_check(struct image *image, UnitCell *cell,
 }
 
 
-static void output_header(FILE *ofh, UnitCell *cell, struct image *image)
-{
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double a, b, c, al, be, ga;
-
-	fprintf(ofh, "Reflections from indexing in %s\n", image->filename);
-
-	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
-	fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm, %7.5f %7.5f %7.5f deg\n",
-	       a*1.0e9, b*1.0e9, c*1.0e9,
-	       rad2deg(al), rad2deg(be), rad2deg(ga));
-
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-	       asx/1e9, asy/1e9, asz/1e9);
-	fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-	       bsx/1e9, bsy/1e9, bsz/1e9);
-	fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-	       csx/1e9, csy/1e9, csz/1e9);
-
-	if ( image->f0_available ) {
-		fprintf(ofh, "f0 = %7.5f (arbitrary gas detector units)\n",
-		       image->f0);
-	} else {
-		fprintf(ofh, "f0 = invalid\n");
-	}
-
-	fprintf(ofh, "photon_energy_eV = %f\n",
-	        J_to_eV(ph_lambda_to_en(image->lambda)));
-
-}
-
-
-void output_intensities(struct image *image, UnitCell *cell,
-                        RefList *reflections, pthread_mutex_t *mutex, int polar,
-                        int use_closer, FILE *ofh)
+/* Integrate the list of predicted reflections in "image" */
+void integrate_reflections(struct image *image, int polar, int use_closer)
 {
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
 	Reflection *refl;
 	RefListIterator *iter;
 
-	/* Get exclusive access to the output stream if necessary */
-	if ( mutex != NULL ) pthread_mutex_lock(mutex);
-
-	output_header(ofh, cell, image);
-
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-
-	for ( refl = first_refl(reflections, &iter);
+	for ( refl = first_refl(image->reflections, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
 
@@ -666,7 +583,6 @@ void output_intensities(struct image *image, UnitCell *cell,
 		double bg, max;
 		struct panel *p;
 		double px, py;
-		signed int h, k, l;
 
 		get_detector_pos(refl, &px, &py);
 		p = find_panel(image->det, px, py);
@@ -734,23 +650,14 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 		}
 
-		/* Write h,k,l, integrated intensity and centroid coordinates */
-		get_indices(refl, &h, &k, &l);
-		fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n",
-		        h, k, l, intensity, x, y, max, bg);
+		set_int(refl, intensity);
 
 	}
-
-	/* Blank line at end */
-	fprintf(ofh, "\n");
-
-	if ( mutex != NULL ) pthread_mutex_unlock(mutex);
 }
 
 
-void output_pixels(struct image *image, UnitCell *cell,
-                   pthread_mutex_t *mutex, int do_polar,
-                   FILE *ofh, int circular_domain, double domain_r)
+RefList *integrate_pixels(struct image *image, int circular_domain,
+                          double domain_r, int do_polar)
 {
 	int i;
 	double ax, ay, az;
@@ -762,26 +669,25 @@ void output_pixels(struct image *image, UnitCell *cell,
 	double *xmom;
 	double *ymom;
 	ReflItemList *obs;
-
-	/* Get exclusive access to the output stream if necessary */
-	if ( mutex != NULL ) pthread_mutex_lock(mutex);
-
-	output_header(ofh, cell, image);
+	RefList *reflections;
 
 	obs = new_items();
 	intensities = new_list_intensity();
 	xmom = new_list_intensity();
 	ymom = new_list_intensity();
+	reflections = reflist_new();
 
 	/* "Borrow" direction values to get reciprocal lengths */
-	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	cell_get_reciprocal(image->indexed_cell, &ax, &ay, &az,
+	                                         &bx, &by, &bz,
+	                                         &cx, &cy, &cz);
 	aslen = modulus(ax, ay, az);
 	bslen = modulus(bx, by, bz);
 	cslen = modulus(cx, cy, cz);
 
-	cell_get_cartesian(cell, &ax, &ay, &az,
-	                         &bx, &by, &bz,
-	                         &cx, &cy, &cz);
+	cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
+	                                        &bx, &by, &bz,
+	                                        &cx, &cy, &cz);
 	/* For each pixel */
 	fesetround(1);  /* Round towards nearest */
 	for ( fs=0; fs<image->width; fs++ ) {
@@ -890,6 +796,7 @@ void output_pixels(struct image *image, UnitCell *cell,
 		struct refl_item *it;
 		double intensity, xmomv, ymomv;
 		double xp, yp;
+		Reflection *refl;
 
 		it = get_item(obs, i);
 		intensity = lookup_intensity(intensities, it->h, it->k, it->l);
@@ -899,19 +806,16 @@ void output_pixels(struct image *image, UnitCell *cell,
 		xp = xmomv / (double)intensity;
 		yp = ymomv / (double)intensity;
 
-		fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f)\n",
-		       it->h, it->k, it->l, intensity, xp, yp);
+		refl = add_refl(reflections, it->h, it->k, it->l);
+		set_int(refl, intensity);
+		set_detector_pos(refl, 0.0, xp, yp);
 
 	}
 
-	fprintf(ofh, "No peak statistics, because output_pixels() was used.\n");
-	/* Blank line at end */
-	fprintf(ofh, "\n");
-
 	free(xmom);
 	free(ymom);
 	free(intensities);
 	delete_items(obs);
 
-	if ( mutex != NULL ) pthread_mutex_unlock(mutex);
+	return reflections;
 }
diff --git a/src/peaks.h b/src/peaks.h
index 2acb249b..7a6557f1 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -23,21 +23,16 @@
 
 extern void search_peaks(struct image *image, float threshold,
                          float min_gradient);
-extern void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex);
 
-extern void output_intensities(struct image *image, UnitCell *cell,
-                               RefList *reflections,
-                               pthread_mutex_t *mutex, int polar,
-                               int use_closer, FILE *ofh);
-
-extern void output_pixels(struct image *image, UnitCell *cell,
-                          pthread_mutex_t *mutex, int do_polar,
-                          FILE *ofh, int circular_domain, double domain_r);
+extern void integrate_reflections(struct image *image,
+                                  int polar, int use_closer);
 
 extern int peak_sanity_check(struct image *image, UnitCell *cell,
                              int circular_domain, double domain_r);
+
 extern RefList *find_projected_peaks(struct image *image, UnitCell *cell,
                                      int circular_domain, double domain_r);
+
 extern int integrate_peak(struct image *image, int xp, int yp,
                           double *xc, double *yc, double *intensity,
                           double *pbg, double *pmax,
diff --git a/src/reintegrate.c b/src/reintegrate.c
index 67e7b347..fcffd499 100644
--- a/src/reintegrate.c
+++ b/src/reintegrate.c
@@ -140,15 +140,14 @@ static void process_image(void *pg, int cookie)
 		reflections = find_projected_peaks(&image, image.indexed_cell,
 		                                   0, 0.1);
 
-		output_intensities(&image, image.indexed_cell, reflections,
-		                   pargs->output_mutex,
-		                   pargs->config_polar,
-		                   pargs->config_closer,
-		                   pargs->ofh);
-
 		reflist_free(reflections);
 	}
 
+	pthread_mutex_lock(pargs->output_mutex);
+	write_chunk(pargs->ofh, &image, pargs->stream_flags);
+	pthread_mutex_unlock(pargs->output_mutex);
+
+
 	free(image.data);
 	if ( image.flags != NULL ) free(image.flags);
 	hdfile_close(hdfile);
diff --git a/src/stream.c b/src/stream.c
index ffee4a6c..fbe0abad 100644
--- a/src/stream.c
+++ b/src/stream.c
@@ -20,6 +20,7 @@
 
 #include "cell.h"
 #include "utils.h"
+#include "image.h"
 
 
 int count_patterns(FILE *fh)
@@ -73,6 +74,84 @@ static UnitCell *read_orientation_matrix(FILE *fh)
 }
 
 
+static void write_reflections(struct image *image, FILE *ofh)
+{
+}
+
+
+static void write_peaks(struct image *image, FILE *ofh)
+{
+	int i;
+
+	fprintf(ofh, "Peaks from peak search in %s\n", image->filename);
+	fprintf(ofh, "  x/px     y/px   (1/d)/nm^-1    Intensity\n");
+
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+		struct imagefeature *f;
+		struct rvec r;
+		double q;
+
+		f = image_get_feature(image->features, i);
+		if ( f == NULL ) continue;
+
+		r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda);
+		q = modulus(r.u, r.v, r.w);
+
+		fprintf(ofh, "%8.3f %8.3f %8.3f    %12.3f\n",
+		       f->x, f->y, q/1.0e9, f->intensity);
+
+	}
+
+	fprintf(ofh, "\n");
+}
+
+
+
+void write_chunk(FILE *ofh, struct image *image, int flags)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double a, b, c, al, be, ga;
+
+	fprintf(ofh, "----- Begin chunk -----\n");
+
+	fprintf(ofh, "Image filename: %s\n", image->filename);
+
+	cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga);
+	fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm,"
+	             " %7.5f %7.5f %7.5f deg\n",
+	             a*1.0e9, b*1.0e9, c*1.0e9,
+	             rad2deg(al), rad2deg(be), rad2deg(ga));
+
+	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+	        asx/1e9, asy/1e9, asz/1e9);
+	fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+	        bsx/1e9, bsy/1e9, bsz/1e9);
+	fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
+	       csx/1e9, csy/1e9, csz/1e9);
+
+	if ( image->f0_available ) {
+		fprintf(ofh, "I0 = %7.5f (arbitrary gas detector units)\n",
+		       image->f0);
+	} else {
+		fprintf(ofh, "I0 = invalid\n");
+	}
+
+	fprintf(ofh, "photon_energy_eV = %f\n",
+	        J_to_eV(ph_lambda_to_en(image->lambda)));
+
+	write_peaks(image, ofh);
+	write_reflections(image, ofh);
+
+	fprintf(ofh, "----- End chunk -----\n\n");
+}
+
+
 int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev)
 {
 	char line[1024];
diff --git a/src/stream.h b/src/stream.h
index 5865aeca..b5f05b19 100644
--- a/src/stream.h
+++ b/src/stream.h
@@ -17,8 +17,11 @@
 #endif
 
 
+struct image;
+
 extern int count_patterns(FILE *fh);
 extern int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev);
+extern void write_chunk(FILE *ofh, struct image *image, int flags);
 
 
 #endif	/* STREAM_H */