Rename calibrate_detector to sum_stack, add new calibrate_detector

5 files changed, 498 insertions, 351 deletions

diff --git a/.gitignore b/.gitignore
index e2762227..6e1ca3f1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,5 +20,6 @@ src/cubeit
 src/reintegrate
 src/estimate_background
 src/check_hkl
+src/sum_stack
 src/.dirstamp
 *~
diff --git a/Makefile.am b/Makefile.am
index b1ab6bf0..8ecce78b 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -5,7 +5,7 @@ ACLOCAL_AMFLAGS = -I m4
 bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
                src/compare_hkl src/powder_plot src/render_hkl \
                src/calibrate_detector src/partialator src/reintegrate \
-               src/estimate_background src/check_hkl
+               src/estimate_background src/check_hkl src/sum_stack
 
 noinst_PROGRAMS = tests/list_check
 
@@ -75,10 +75,13 @@ src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \
                          src/hdf5-file.c src/image.c src/filters.c \
                          src/thread-pool.c
 
+src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \
+                        src/image.c src/filters.c src/peaks.c src/detector.c \
+                        src/cell.c src/thread-pool.c src/reflist.c
+
 src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \
-                                 src/hdf5-file.c src/image.c src/filters.c \
-                                 src/peaks.c src/detector.c src/cell.c \
-                                 src/thread-pool.c src/reflist.c
+                                 src/hdf5-file.c src/image.c src/thread-pool.c \
+                                 src/detector.c
 
 src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
                           src/utils.c src/detector.c src/peaks.c src/image.c \
diff --git a/Makefile.in b/Makefile.in
index 89fe9a35..e1f68a5c 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -41,7 +41,7 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
 	src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \
 	src/partialator$(EXEEXT) src/reintegrate$(EXEEXT) \
 	src/estimate_background$(EXEEXT) src/check_hkl$(EXEEXT) \
-	$(am__EXEEXT_1) $(am__EXEEXT_2)
+	src/sum_stack$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
 noinst_PROGRAMS = tests/list_check$(EXEEXT)
 TESTS = tests/list_check$(EXEEXT)
 @HAVE_GTK_TRUE@am__append_1 = src/hdfsee
@@ -83,9 +83,8 @@ PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS)
 am__dirstamp = $(am__leading_dot)dirstamp
 am_src_calibrate_detector_OBJECTS = src/calibrate_detector.$(OBJEXT) \
 	src/utils.$(OBJEXT) src/hdf5-file.$(OBJEXT) \
-	src/image.$(OBJEXT) src/filters.$(OBJEXT) src/peaks.$(OBJEXT) \
-	src/detector.$(OBJEXT) src/cell.$(OBJEXT) \
-	src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT)
+	src/image.$(OBJEXT) src/thread-pool.$(OBJEXT) \
+	src/detector.$(OBJEXT)
 src_calibrate_detector_OBJECTS = $(am_src_calibrate_detector_OBJECTS)
 src_calibrate_detector_LDADD = $(LDADD)
 src_calibrate_detector_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
@@ -224,6 +223,14 @@ am_src_render_hkl_OBJECTS = src/render_hkl.$(OBJEXT) \
 src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS)
 src_render_hkl_LDADD = $(LDADD)
 src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
+am_src_sum_stack_OBJECTS = src/sum_stack.$(OBJEXT) src/utils.$(OBJEXT) \
+	src/hdf5-file.$(OBJEXT) src/image.$(OBJEXT) \
+	src/filters.$(OBJEXT) src/peaks.$(OBJEXT) \
+	src/detector.$(OBJEXT) src/cell.$(OBJEXT) \
+	src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT)
+src_sum_stack_OBJECTS = $(am_src_sum_stack_OBJECTS)
+src_sum_stack_LDADD = $(LDADD)
+src_sum_stack_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 am_tests_list_check_OBJECTS = tests/list_check.$(OBJEXT) \
 	src/reflist.$(OBJEXT) src/thread-pool.$(OBJEXT) \
 	src/utils.$(OBJEXT)
@@ -257,7 +264,7 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
 	$(src_partialator_SOURCES) $(src_pattern_sim_SOURCES) \
 	$(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \
 	$(src_reintegrate_SOURCES) $(src_render_hkl_SOURCES) \
-	$(tests_list_check_SOURCES)
+	$(src_sum_stack_SOURCES) $(tests_list_check_SOURCES)
 DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
 	$(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \
 	$(am__src_cubeit_SOURCES_DIST) \
@@ -266,7 +273,8 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
 	$(am__src_indexamajig_SOURCES_DIST) $(src_partialator_SOURCES) \
 	$(am__src_pattern_sim_SOURCES_DIST) $(src_powder_plot_SOURCES) \
 	$(src_process_hkl_SOURCES) $(src_reintegrate_SOURCES) \
-	$(src_render_hkl_SOURCES) $(tests_list_check_SOURCES)
+	$(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \
+	$(tests_list_check_SOURCES)
 RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
 	html-recursive info-recursive install-data-recursive \
 	install-dvi-recursive install-exec-recursive \
@@ -614,10 +622,13 @@ src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \
                          src/hdf5-file.c src/image.c src/filters.c \
                          src/thread-pool.c
 
+src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \
+                        src/image.c src/filters.c src/peaks.c src/detector.c \
+                        src/cell.c src/thread-pool.c src/reflist.c
+
 src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \
-                                 src/hdf5-file.c src/image.c src/filters.c \
-                                 src/peaks.c src/detector.c src/cell.c \
-                                 src/thread-pool.c src/reflist.c
+                                 src/hdf5-file.c src/image.c src/thread-pool.c \
+                                 src/detector.c
 
 src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \
                           src/utils.c src/detector.c src/peaks.c src/image.c \
@@ -782,15 +793,9 @@ src/utils.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
 src/hdf5-file.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/image.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
-src/filters.$(OBJEXT): src/$(am__dirstamp) \
-	src/$(DEPDIR)/$(am__dirstamp)
-src/peaks.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
-src/detector.$(OBJEXT): src/$(am__dirstamp) \
-	src/$(DEPDIR)/$(am__dirstamp)
-src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
 src/thread-pool.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
-src/reflist.$(OBJEXT): src/$(am__dirstamp) \
+src/detector.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibrate_detector_DEPENDENCIES) src/$(am__dirstamp)
 	@rm -f src/calibrate_detector$(EXEEXT)
@@ -798,6 +803,7 @@ src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibra
 src/check_hkl.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/sfac.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
+src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
 src/reflections.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/statistics.$(OBJEXT): src/$(am__dirstamp) \
@@ -816,6 +822,8 @@ src/cubeit.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/render.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
+src/filters.$(OBJEXT): src/$(am__dirstamp) \
+	src/$(DEPDIR)/$(am__dirstamp)
 src/stream.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp)
@@ -842,6 +850,7 @@ src/hdfsee$(EXEEXT): $(src_hdfsee_OBJECTS) $(src_hdfsee_DEPENDENCIES) src/$(am__
 	$(AM_V_CCLD)$(LINK) $(src_hdfsee_OBJECTS) $(src_hdfsee_LDADD) $(LIBS)
 src/indexamajig.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
+src/peaks.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
 src/index.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp)
 src/diffraction.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
@@ -852,6 +861,8 @@ src/templates.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/geometry.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
+src/reflist.$(OBJEXT): src/$(am__dirstamp) \
+	src/$(DEPDIR)/$(am__dirstamp)
 src/diffraction-gpu.$(OBJEXT): src/$(am__dirstamp) \
 	src/$(DEPDIR)/$(am__dirstamp)
 src/cl-utils.$(OBJEXT): src/$(am__dirstamp) \
@@ -895,6 +906,11 @@ src/povray.$(OBJEXT): src/$(am__dirstamp) \
 src/render_hkl$(EXEEXT): $(src_render_hkl_OBJECTS) $(src_render_hkl_DEPENDENCIES) src/$(am__dirstamp)
 	@rm -f src/render_hkl$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(src_render_hkl_OBJECTS) $(src_render_hkl_LDADD) $(LIBS)
+src/sum_stack.$(OBJEXT): src/$(am__dirstamp) \
+	src/$(DEPDIR)/$(am__dirstamp)
+src/sum_stack$(EXEEXT): $(src_sum_stack_OBJECTS) $(src_sum_stack_DEPENDENCIES) src/$(am__dirstamp)
+	@rm -f src/sum_stack$(EXEEXT)
+	$(AM_V_CCLD)$(LINK) $(src_sum_stack_OBJECTS) $(src_sum_stack_LDADD) $(LIBS)
 tests/$(am__dirstamp):
 	@$(MKDIR_P) tests
 	@: > tests/$(am__dirstamp)
@@ -947,6 +963,7 @@ mostlyclean-compile:
 	-rm -f src/sfac.$(OBJEXT)
 	-rm -f src/statistics.$(OBJEXT)
 	-rm -f src/stream.$(OBJEXT)
+	-rm -f src/sum_stack.$(OBJEXT)
 	-rm -f src/symmetry.$(OBJEXT)
 	-rm -f src/templates.$(OBJEXT)
 	-rm -f src/thread-pool.$(OBJEXT)
@@ -994,6 +1011,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sfac.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/statistics.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/stream.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/symmetry.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/templates.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Po@am__quote@
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 417ee36e..25e9850c 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -1,9 +1,9 @@
 /*
- * calibrate-detector.c
+ * calibrate_detector.c
  *
- * Detector calibration
+ * Attempt to refine detector geometry
  *
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -20,306 +20,61 @@
 #include <string.h>
 #include <unistd.h>
 #include <getopt.h>
-#include <pthread.h>
 
 #include "utils.h"
+#include "image.h"
+#include "detector.h"
+#include "index.h"
 #include "hdf5-file.h"
-#include "filters.h"
-#include "peaks.h"
-#include "thread-pool.h"
-
-
-#define INTEGRATION_RADIUS (10)
-
-
-typedef enum
-{
-	SUM_THRESHOLD,
-	SUM_PEAKS
-} SumMethod;
-
-struct sum_args
-{
-	char *filename;
-	int config_cmfilter;
-	int config_noisefilter;
-	double *sum;
-	int w;
-	int h;
-	SumMethod sum_method;
-	double threshold;
-};
-
-
-struct queue_args
-{
-	FILE *fh;
-	char *prefix;
-	int config_cmfilter;
-	int config_noisefilter;
-	double *sum;
-	int w;
-	int h;
-	SumMethod sum_method;
-	double threshold;
-};
 
 
 static void show_help(const char *s)
 {
-	printf("Syntax: %s [options]\n\n", s);
+	printf("Syntax: %s [options] <file.h5>\n\n", s);
 	printf(
-"Calibrate detector geometry from FEL diffraction images.\n"
-"\n"
-"  -h, --help              Display this help message.\n"
-"\n"
-"  -i, --input=<filename>  Specify file containing list of images to process.\n"
-"                           '-' means stdin, which is the default.\n"
-"  -o, --output=<filename> Output filename for summed image in HDF5 format.\n"
-"                           Default: summed.h5.\n"
-"\n"
-"  -p, --intermediate=<P>  Stem of filename for intermediate images.\n"
-"                            The filename stem <p> will be postfixed with a\n"
-"                            hyphen, the current number of patterns processed\n"
-"                            and '.h5'.  Such a pattern will be saved after\n"
-"                            every 1000 input patterns.\n"
-"                            If this option is not specified, no intermediate\n"
-"                            patterns will be saved.\n"
-"\n"
-"  -s, --sum=<method>      Use this method for summation.  Choose from:\n"
-"                           peaks : sum 10px radius circles around peaks.\n"
-"                           threshold : sum thresholded images.\n"
-"  -t, --threshold=<n>     Set the threshold if summing using the 'threshold'\n"
-"                           method.\n"
+"Plot a powder pattern as a 1D graph using the detector geometry.\n"
 "\n"
-"      --filter-cm         Perform common-mode noise subtraction on images\n"
-"                           before proceeding.\n"
-"      --filter-noise      Apply an aggressive noise filter which sets all\n"
-"                           pixels in each 3x3 region to zero if any of them\n"
-"                           have negative values.\n"
-"\n"
-" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
-" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n");
-}
-
-
-static void sum_peaks(struct image *image, double *sum)
-{
-	int x, y, i;
-	int w = image->width;
-	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
-
-	/* FIXME: Get threshold value from command line */
-	search_peaks(image, 800.0, 100000.0);
-
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f = image_get_feature(image->features, i);
-		int xp, yp;
-
-		/* This is not an error. */
-		if ( f == NULL ) continue;
-
-		xp = f->x;
-		yp = f->y;
-
-		for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) {
-		for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) {
-
-			/* Circular mask */
-			if ( x*x + y*y > lim ) continue;
-
-			if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue;
-			if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue;
-
-			float val = image->data[(x+xp)+w*(y+yp)];
-			sum[(x+xp)+w*(y+yp)] += val;
-
-		}
-		}
-
-	}
-}
-
-
-static void sum_threshold(struct image *image, double *sum, double threshold)
-{
-	int x, y;
-
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
-		float val = image->data[x+image->width*y];
-		if ( val > threshold ) {
-			sum[x+image->width*y] += val;
-		}
-	}
-	}
-}
-
-
-static void add_image(void *args, int cookie)
-{
-	struct sum_args *pargs = args;
-	struct hdfile *hdfile;
-	struct image image;
-
-	image.features = NULL;
-	image.data = NULL;
-	image.flags = NULL;
-	image.indexed_cell = NULL;
-	image.filename = pargs->filename;
-	image.det = NULL;
-
-	STATUS("%3i: Processing '%s'\n", cookie, pargs->filename);
-
-	hdfile = hdfile_open(pargs->filename);
-	if ( hdfile == NULL ) {
-		return;
-	} else if ( hdfile_set_first_image(hdfile, "/") ) {
-		ERROR("Couldn't select path\n");
-		hdfile_close(hdfile);
-		return;
-	}
-
-	/* FIXME: Nominal photon energy */
-	hdf5_read(hdfile, &image, 1, 2000.0);
-
-	if ( pargs->config_cmfilter ) {
-		filter_cm(&image);
-	}
-
-	if ( pargs->config_noisefilter ) {
-		filter_noise(&image, NULL);
-	}
-
-	if ( (pargs->w != image.width) || (pargs->h != image.height) ) {
-		ERROR("Wrong image size.\n");
-		goto out;
-	}
-
-	switch ( pargs->sum_method ) {
-
-	case SUM_THRESHOLD :
-		sum_threshold(&image, pargs->sum, pargs->threshold);
-		break;
-
-	case SUM_PEAKS :
-		sum_peaks(&image, pargs->sum);
-		break;
-
-	}
-
-out:
-	free(image.data);
-	image_feature_list_free(image.features);
-	if ( image.flags != NULL ) free(image.flags);
-	hdfile_close(hdfile);
-
-	free(pargs->filename);
-	free(pargs);
-}
-
-
-static void *get_image(void *qp)
-{
-	char line[1024];
-	struct sum_args *pargs;
-	char *rval;
-	struct queue_args *qargs = qp;
-
-	/* Get the next filename */
-	rval = fgets(line, 1023, qargs->fh);
-	if ( rval == NULL ) return NULL;
-
-	pargs = malloc(sizeof(struct sum_args));
-
-	pargs->w = qargs->w;
-	pargs->h = qargs->h;
-	pargs->sum_method = qargs->sum_method;
-	pargs->threshold = qargs->threshold;
-	pargs->config_cmfilter = qargs->config_cmfilter;
-	pargs->config_noisefilter = qargs->config_noisefilter;
-	pargs->sum = qargs->sum;
-
-	chomp(line);
-	pargs->filename = malloc(strlen(qargs->prefix) + strlen(line) + 1);
-	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);
-
-	return pargs;
+"  -h, --help                 Display this help message.\n"
+"  -g. --geometry=<file>      Get detector geometry from file.\n"
+"  -i, --input=<file>         Input filename.\n"
+"\n");
 }
 
 
 int main(int argc, char *argv[])
 {
 	int c;
+	struct image image;
+	int x, y;
+	struct hdfile *hdfile;
 	char *filename = NULL;
-	char *outfile = NULL;
-	FILE *fh;
-	int n_images = 0;
-	int config_cmfilter = 0;
-	int config_noisefilter = 0;
-	char *prefix = NULL;
-	char *sum_str = NULL;
-	char *intermediate = NULL;
-	double threshold = 400.0;
-	SumMethod sum;
-	int nthreads = 1;
-	struct queue_args qargs;
-	int n_done;
-	const int chunk_size = 1000;
+	char *geometry = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
 		{"input",              1, NULL,               'i'},
-		{"output",             1, NULL,               'o'},
-		{"filter-cm",          0, &config_cmfilter,    1},
-		{"filter-noise",       0, &config_noisefilter, 1},
-		{"prefix",             1, NULL,               'x'},
-		{"sum",                1, NULL,               's'},
-		{"intermediate",       1, NULL,               'p'},
-		{"threshold",          1, NULL,               't'},
+		{"geometry",           1, NULL,               'g'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:x:j:o:s:p:t:",
-	                        longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hi:g:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' :
 			show_help(argv[0]);
 			return 0;
 
-		case 'i' :
-			filename = strdup(optarg);
-			break;
-
-		case 'o' :
-			outfile = strdup(optarg);
-			break;
-
-		case 'x' :
-			prefix = strdup(optarg);
-			break;
-
-		case 'j' :
-			nthreads = atoi(optarg);
-			break;
-
-		case 's' :
-			sum_str = strdup(optarg);
-			break;
-
-		case 'p' :
-			intermediate = strdup(optarg);
+		case 0 :
 			break;
 
-		case 't' :
-			threshold = atof(optarg);
+		case 'i' :
+			filename = strdup(optarg);
 			break;
 
-		case 0 :
+		case 'g' :
+			geometry = strdup(optarg);
 			break;
 
 		default :
@@ -329,84 +84,42 @@ int main(int argc, char *argv[])
 	}
 
 	if ( filename == NULL ) {
-		filename = strdup("-");
-	}
-	if ( strcmp(filename, "-") == 0 ) {
-		fh = stdin;
-	} else {
-		fh = fopen(filename, "r");
-	}
-	if ( fh == NULL ) {
-		ERROR("Failed to open input file '%s'\n", filename);
+		ERROR("You must specify the input filename with -i\n");
 		return 1;
 	}
-	free(filename);
 
-	if ( sum_str == NULL ) {
-		STATUS("You didn't specify a summation method, so I'm using"
-		       " the 'peaks' method, which gives the best results.\n");
-		sum = SUM_PEAKS;
-	} else if ( strcmp(sum_str, "peaks") == 0 ) {
-		sum = SUM_PEAKS;
-	} else if ( strcmp(sum_str, "threshold") == 0) {
-		sum = SUM_THRESHOLD;
-	} else {
-		ERROR("Unrecognised summation method '%s'\n", sum_str);
+	if ( geometry == NULL ) {
+		ERROR("You need to specify a geometry file with --geometry\n");
 		return 1;
 	}
-	free(sum_str);
-
-	if ( prefix == NULL ) {
-		prefix = strdup("");
-	}
-
-	if ( outfile == NULL ) {
-		outfile = strdup("summed.h5");
-	}
 
-	if ( nthreads == 0  ) {
-		ERROR("Invalid number of threads.\n");
+	image.det = get_detector_geometry(geometry);
+	if ( image.det == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n", geometry);
 		return 1;
 	}
+	free(geometry);
 
-	qargs.w = 1024;  /* FIXME! */
-	qargs.h = 1024;  /* FIXME! */
-	qargs.sum_method = sum;
-	qargs.threshold = threshold;
-	qargs.config_cmfilter = config_cmfilter;
-	qargs.config_noisefilter = config_noisefilter;
-	qargs.sum = calloc(qargs.w*qargs.h, sizeof(double));
-	qargs.prefix = prefix;
-	qargs.fh = fh;
-
-	do {
-
-		n_done = run_threads(nthreads, add_image, get_image,
-		                     (void *)&qargs, NULL, chunk_size);
+	hdfile = hdfile_open(filename);
+	hdfile_set_image(hdfile, "/data/data");
+	hdf5_read(hdfile, &image, 1, 2000.0);
 
-		n_images += n_done;
+	for ( x=0; x<image.width; x++ ) {
+	for ( y=0; y<image.height; y++ ) {
 
-		/* Write intermediate sum if requested */
-		if ( (intermediate != NULL) && (n_done == chunk_size) ) {
-			char outfile[1024];
-			snprintf(outfile, 1023, "%s-%i.h5",
-			         intermediate, n_images);
-			hdf5_write(outfile, qargs.sum, qargs.w, qargs.h,
-			           H5T_NATIVE_DOUBLE);
-		}
+		double q;
+		int intensity;
+		struct rvec r;
 
-	} while ( n_done == chunk_size );
+		r = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+		q = modulus(r.u, r.v, r.w);
 
-	/* Write the final output */
-	hdf5_write(outfile, qargs.sum, qargs.w, qargs.h, H5T_NATIVE_DOUBLE);
+		intensity = image.data[x + image.width*y];
 
-	free(qargs.sum);
-	free(prefix);
-	free(outfile);
-	if ( intermediate != NULL ) free(intermediate);
-	fclose(fh);
+		printf("%5i\t%5i\t%7.3f\t%7i\n", x, y, q/1.0e9, intensity);
 
-	STATUS("There were %i images.\n", n_images);
+	}
+	}
 
 	return 0;
 }
diff --git a/src/sum_stack.c b/src/sum_stack.c
new file mode 100644
index 00000000..8bd5197e
--- /dev/null
+++ b/src/sum_stack.c
@@ -0,0 +1,412 @@
+/*
+ * sum_stack.c
+ *
+ * Sum a stack of images (e.g. for detector calibration)
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <pthread.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "filters.h"
+#include "peaks.h"
+#include "thread-pool.h"
+
+
+#define INTEGRATION_RADIUS (10)
+
+
+typedef enum
+{
+	SUM_THRESHOLD,
+	SUM_PEAKS
+} SumMethod;
+
+struct sum_args
+{
+	char *filename;
+	int config_cmfilter;
+	int config_noisefilter;
+	double *sum;
+	int w;
+	int h;
+	SumMethod sum_method;
+	double threshold;
+};
+
+
+struct queue_args
+{
+	FILE *fh;
+	char *prefix;
+	int config_cmfilter;
+	int config_noisefilter;
+	double *sum;
+	int w;
+	int h;
+	SumMethod sum_method;
+	double threshold;
+};
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options]\n\n", s);
+	printf(
+"Sum FEL diffraction images.\n"
+"\n"
+"  -h, --help              Display this help message.\n"
+"\n"
+"  -i, --input=<filename>  Specify file containing list of images to process.\n"
+"                           '-' means stdin, which is the default.\n"
+"  -o, --output=<filename> Output filename for summed image in HDF5 format.\n"
+"                           Default: summed.h5.\n"
+"\n"
+"  -p, --intermediate=<P>  Stem of filename for intermediate images.\n"
+"                            The filename stem <p> will be postfixed with a\n"
+"                            hyphen, the current number of patterns processed\n"
+"                            and '.h5'.  Such a pattern will be saved after\n"
+"                            every 1000 input patterns.\n"
+"                            If this option is not specified, no intermediate\n"
+"                            patterns will be saved.\n"
+"\n"
+"  -s, --sum=<method>      Use this method for summation.  Choose from:\n"
+"                           peaks : sum 10px radius circles around peaks.\n"
+"                           threshold : sum thresholded images.\n"
+"  -t, --threshold=<n>     Set the threshold if summing using the 'threshold'\n"
+"                           method.\n"
+"\n"
+"      --filter-cm         Perform common-mode noise subtraction on images\n"
+"                           before proceeding.\n"
+"      --filter-noise      Apply an aggressive noise filter which sets all\n"
+"                           pixels in each 3x3 region to zero if any of them\n"
+"                           have negative values.\n"
+"\n"
+" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
+" -x, --prefix=<p>         Prefix filenames from input file with 'p'.\n");
+}
+
+
+static void sum_peaks(struct image *image, double *sum)
+{
+	int x, y, i;
+	int w = image->width;
+	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
+
+	/* FIXME: Get threshold value from command line */
+	search_peaks(image, 800.0, 100000.0);
+
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+		struct imagefeature *f = image_get_feature(image->features, i);
+		int xp, yp;
+
+		/* This is not an error. */
+		if ( f == NULL ) continue;
+
+		xp = f->x;
+		yp = f->y;
+
+		for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) {
+		for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) {
+
+			/* Circular mask */
+			if ( x*x + y*y > lim ) continue;
+
+			if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue;
+			if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue;
+
+			float val = image->data[(x+xp)+w*(y+yp)];
+			sum[(x+xp)+w*(y+yp)] += val;
+
+		}
+		}
+
+	}
+}
+
+
+static void sum_threshold(struct image *image, double *sum, double threshold)
+{
+	int x, y;
+
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+		float val = image->data[x+image->width*y];
+		if ( val > threshold ) {
+			sum[x+image->width*y] += val;
+		}
+	}
+	}
+}
+
+
+static void add_image(void *args, int cookie)
+{
+	struct sum_args *pargs = args;
+	struct hdfile *hdfile;
+	struct image image;
+
+	image.features = NULL;
+	image.data = NULL;
+	image.flags = NULL;
+	image.indexed_cell = NULL;
+	image.filename = pargs->filename;
+	image.det = NULL;
+
+	STATUS("%3i: Processing '%s'\n", cookie, pargs->filename);
+
+	hdfile = hdfile_open(pargs->filename);
+	if ( hdfile == NULL ) {
+		return;
+	} else if ( hdfile_set_first_image(hdfile, "/") ) {
+		ERROR("Couldn't select path\n");
+		hdfile_close(hdfile);
+		return;
+	}
+
+	/* FIXME: Nominal photon energy */
+	hdf5_read(hdfile, &image, 1, 2000.0);
+
+	if ( pargs->config_cmfilter ) {
+		filter_cm(&image);
+	}
+
+	if ( pargs->config_noisefilter ) {
+		filter_noise(&image, NULL);
+	}
+
+	if ( (pargs->w != image.width) || (pargs->h != image.height) ) {
+		ERROR("Wrong image size.\n");
+		goto out;
+	}
+
+	switch ( pargs->sum_method ) {
+
+	case SUM_THRESHOLD :
+		sum_threshold(&image, pargs->sum, pargs->threshold);
+		break;
+
+	case SUM_PEAKS :
+		sum_peaks(&image, pargs->sum);
+		break;
+
+	}
+
+out:
+	free(image.data);
+	image_feature_list_free(image.features);
+	if ( image.flags != NULL ) free(image.flags);
+	hdfile_close(hdfile);
+
+	free(pargs->filename);
+	free(pargs);
+}
+
+
+static void *get_image(void *qp)
+{
+	char line[1024];
+	struct sum_args *pargs;
+	char *rval;
+	struct queue_args *qargs = qp;
+
+	/* Get the next filename */
+	rval = fgets(line, 1023, qargs->fh);
+	if ( rval == NULL ) return NULL;
+
+	pargs = malloc(sizeof(struct sum_args));
+
+	pargs->w = qargs->w;
+	pargs->h = qargs->h;
+	pargs->sum_method = qargs->sum_method;
+	pargs->threshold = qargs->threshold;
+	pargs->config_cmfilter = qargs->config_cmfilter;
+	pargs->config_noisefilter = qargs->config_noisefilter;
+	pargs->sum = qargs->sum;
+
+	chomp(line);
+	pargs->filename = malloc(strlen(qargs->prefix) + strlen(line) + 1);
+	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);
+
+	return pargs;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *filename = NULL;
+	char *outfile = NULL;
+	FILE *fh;
+	int n_images = 0;
+	int config_cmfilter = 0;
+	int config_noisefilter = 0;
+	char *prefix = NULL;
+	char *sum_str = NULL;
+	char *intermediate = NULL;
+	double threshold = 400.0;
+	SumMethod sum;
+	int nthreads = 1;
+	struct queue_args qargs;
+	int n_done;
+	const int chunk_size = 1000;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"input",              1, NULL,               'i'},
+		{"output",             1, NULL,               'o'},
+		{"filter-cm",          0, &config_cmfilter,    1},
+		{"filter-noise",       0, &config_noisefilter, 1},
+		{"prefix",             1, NULL,               'x'},
+		{"sum",                1, NULL,               's'},
+		{"intermediate",       1, NULL,               'p'},
+		{"threshold",          1, NULL,               't'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:x:j:o:s:p:t:",
+	                        longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' :
+			show_help(argv[0]);
+			return 0;
+
+		case 'i' :
+			filename = strdup(optarg);
+			break;
+
+		case 'o' :
+			outfile = strdup(optarg);
+			break;
+
+		case 'x' :
+			prefix = strdup(optarg);
+			break;
+
+		case 'j' :
+			nthreads = atoi(optarg);
+			break;
+
+		case 's' :
+			sum_str = strdup(optarg);
+			break;
+
+		case 'p' :
+			intermediate = strdup(optarg);
+			break;
+
+		case 't' :
+			threshold = atof(optarg);
+			break;
+
+		case 0 :
+			break;
+
+		default :
+			return 1;
+		}
+
+	}
+
+	if ( filename == NULL ) {
+		filename = strdup("-");
+	}
+	if ( strcmp(filename, "-") == 0 ) {
+		fh = stdin;
+	} else {
+		fh = fopen(filename, "r");
+	}
+	if ( fh == NULL ) {
+		ERROR("Failed to open input file '%s'\n", filename);
+		return 1;
+	}
+	free(filename);
+
+	if ( sum_str == NULL ) {
+		STATUS("You didn't specify a summation method, so I'm using"
+		       " the 'peaks' method, which gives the best results.\n");
+		sum = SUM_PEAKS;
+	} else if ( strcmp(sum_str, "peaks") == 0 ) {
+		sum = SUM_PEAKS;
+	} else if ( strcmp(sum_str, "threshold") == 0) {
+		sum = SUM_THRESHOLD;
+	} else {
+		ERROR("Unrecognised summation method '%s'\n", sum_str);
+		return 1;
+	}
+	free(sum_str);
+
+	if ( prefix == NULL ) {
+		prefix = strdup("");
+	}
+
+	if ( outfile == NULL ) {
+		outfile = strdup("summed.h5");
+	}
+
+	if ( nthreads == 0  ) {
+		ERROR("Invalid number of threads.\n");
+		return 1;
+	}
+
+	qargs.w = 1024;  /* FIXME! */
+	qargs.h = 1024;  /* FIXME! */
+	qargs.sum_method = sum;
+	qargs.threshold = threshold;
+	qargs.config_cmfilter = config_cmfilter;
+	qargs.config_noisefilter = config_noisefilter;
+	qargs.sum = calloc(qargs.w*qargs.h, sizeof(double));
+	qargs.prefix = prefix;
+	qargs.fh = fh;
+
+	do {
+
+		n_done = run_threads(nthreads, add_image, get_image,
+		                     (void *)&qargs, NULL, chunk_size);
+
+		n_images += n_done;
+
+		/* Write intermediate sum if requested */
+		if ( (intermediate != NULL) && (n_done == chunk_size) ) {
+			char outfile[1024];
+			snprintf(outfile, 1023, "%s-%i.h5",
+			         intermediate, n_images);
+			hdf5_write(outfile, qargs.sum, qargs.w, qargs.h,
+			           H5T_NATIVE_DOUBLE);
+		}
+
+	} while ( n_done == chunk_size );
+
+	/* Write the final output */
+	hdf5_write(outfile, qargs.sum, qargs.w, qargs.h, H5T_NATIVE_DOUBLE);
+
+	free(qargs.sum);
+	free(prefix);
+	free(outfile);
+	if ( intermediate != NULL ) free(intermediate);
+	fclose(fh);
+
+	STATUS("There were %i images.\n", n_images);
+
+	return 0;
+}