Remove sum_stack from core

4 files changed, 4 insertions, 548 deletions

diff --git a/Makefile.am b/Makefile.am
index 53321656..598b8ee0 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/check_hkl src/sum_stack src/partial_sim
+               src/check_hkl src/partial_sim
 
 noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/pr_gradient_check tests/symmetry_check
@@ -58,8 +58,6 @@ src_powder_plot_SOURCES = src/powder_plot.c
 
 src_render_hkl_SOURCES = src/render_hkl.c src/povray.c
 
-src_sum_stack_SOURCES = src/sum_stack.c
-
 src_calibrate_detector_SOURCES = src/calibrate_detector.c
 
 src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
diff --git a/Makefile.in b/Makefile.in
index bffe2b0d..77193584 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -40,8 +40,7 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
 	src/compare_hkl$(EXEEXT) src/powder_plot$(EXEEXT) \
 	src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \
 	src/partialator$(EXEEXT) src/check_hkl$(EXEEXT) \
-	src/sum_stack$(EXEEXT) src/partial_sim$(EXEEXT) \
-	$(am__EXEEXT_1)
+	src/partial_sim$(EXEEXT) $(am__EXEEXT_1)
 noinst_PROGRAMS = tests/list_check$(EXEEXT) \
 	tests/integration_check$(EXEEXT) \
 	tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
@@ -162,11 +161,6 @@ src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS)
 src_render_hkl_LDADD = $(LDADD)
 src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \
 	$(top_builddir)/libcrystfel/libcrystfel.la
-am_src_sum_stack_OBJECTS = src/sum_stack.$(OBJEXT)
-src_sum_stack_OBJECTS = $(am_src_sum_stack_OBJECTS)
-src_sum_stack_LDADD = $(LDADD)
-src_sum_stack_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \
-	$(top_builddir)/libcrystfel/libcrystfel.la
 am__tests_gpu_sim_check_SOURCES_DIST = tests/gpu_sim_check.c
 @HAVE_OPENCL_TRUE@am_tests_gpu_sim_check_OBJECTS =  \
 @HAVE_OPENCL_TRUE@	tests/gpu_sim_check.$(OBJEXT)
@@ -231,7 +225,7 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
 	$(src_partial_sim_SOURCES) $(src_partialator_SOURCES) \
 	$(src_pattern_sim_SOURCES) $(src_powder_plot_SOURCES) \
 	$(src_process_hkl_SOURCES) $(src_render_hkl_SOURCES) \
-	$(src_sum_stack_SOURCES) $(tests_gpu_sim_check_SOURCES) \
+	$(tests_gpu_sim_check_SOURCES) \
 	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
 	$(tests_pr_gradient_check_SOURCES) \
 	$(tests_symmetry_check_SOURCES)
@@ -241,7 +235,7 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \
 	$(src_indexamajig_SOURCES) $(src_partial_sim_SOURCES) \
 	$(am__src_partialator_SOURCES_DIST) $(src_pattern_sim_SOURCES) \
 	$(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \
-	$(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \
+	$(src_render_hkl_SOURCES) \
 	$(am__tests_gpu_sim_check_SOURCES_DIST) \
 	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
 	$(tests_pr_gradient_check_SOURCES) \
@@ -581,7 +575,6 @@ src_compare_hkl_SOURCES = src/compare_hkl.c
 src_check_hkl_SOURCES = src/check_hkl.c
 src_powder_plot_SOURCES = src/powder_plot.c
 src_render_hkl_SOURCES = src/render_hkl.c src/povray.c
-src_sum_stack_SOURCES = src/sum_stack.c
 src_calibrate_detector_SOURCES = src/calibrate_detector.c
 src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
 	src/hrs-scaling.c $(am__append_4)
@@ -802,11 +795,6 @@ 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)
@@ -858,7 +846,6 @@ mostlyclean-compile:
 	-rm -f src/process_hkl.$(OBJEXT)
 	-rm -f src/render_hkl.$(OBJEXT)
 	-rm -f src/scaling-report.$(OBJEXT)
-	-rm -f src/sum_stack.$(OBJEXT)
 	-rm -f tests/gpu_sim_check.$(OBJEXT)
 	-rm -f tests/integration_check.$(OBJEXT)
 	-rm -f tests/list_check.$(OBJEXT)
@@ -885,7 +872,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render_hkl.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/scaling-report.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/gpu_sim_check.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@
diff --git a/README b/README
index cfdb5049..87a128b2 100644
--- a/README
+++ b/README
@@ -52,11 +52,6 @@ In addition, there is also:
 
 	- check_hkl, for determining things like completeness.
 
-	- sum_stack, for summing patterns after peak detection to use
-	             for calibrating your detector.
-
-	- cubeit, for summing patterns in 3D to inspect the peak profiles.
-
 	- calibrate_detector, for refining detector geometry.
 
 	- partial_sim, for calculating partial reflection intensities.
diff --git a/src/sum_stack.c b/src/sum_stack.c
deleted file mode 100644
index a6e7599a..00000000
--- a/src/sum_stack.c
+++ /dev/null
@@ -1,523 +0,0 @@
-/*
- * 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;
-	double min_gradient;
-	double min_snr;
-	char *element;
-};
-
-
-struct queue_args
-{
-	FILE *fh;
-	int config_cmfilter;
-	int config_noisefilter;
-	double *sum;
-	int w;
-	int h;
-	SumMethod sum_method;
-	double threshold;
-	double min_gradient;
-	double min_snr;
-	char *element;
-
-	char *use_this_one_instead;
-	char *prefix;
-	int config_basename;
-};
-
-
-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.  Default: 400 adu\n"
-"      --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
-"                            Default: 100,000.\n"
-"      --min-snr=<n>        Minimum signal-to-noise ration for peaks.\n"
-"                            Default: 5.\n"
-" -e, --image=<element>    Use this image from the HDF5 file.\n"
-"                           Example: /data/data0.\n"
-"                           Default: The first one found.\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, double threshold,
-                      double min_gradient, double min_snr)
-{
-	int i;
-	int w = image->width;
-	int h = image->height;
-	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
-
-	search_peaks(image, threshold, min_gradient, min_snr);
-
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f = image_get_feature(image->features, i);
-		int pfs, pss;
-		int fs, ss;
-
-		/* This is not an error. */
-		if ( f == NULL ) continue;
-
-		pfs = f->fs;
-		pss = f->ss;
-
-		for ( fs=-INTEGRATION_RADIUS; fs<+INTEGRATION_RADIUS; fs++ ) {
-		for ( ss=-INTEGRATION_RADIUS; ss<+INTEGRATION_RADIUS; ss++ ) {
-
-			/* Circular mask */
-			if ( fs*fs + ss*ss > lim ) continue;
-
-			if ( ((fs+pfs)>=w) || ((fs+pfs)<0) ) continue;
-			if ( ((ss+pss)>=h) || ((ss+pss)<0) ) continue;
-
-			float val = image->data[(fs+pfs)+w*(ss+pss)];
-			sum[(fs+pfs)+w*(ss+pss)] += 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;
-
-	if ( pargs->element != NULL ) {
-		int r;
-		r = hdfile_set_image(hdfile, pargs->element);
-		if ( r ) {
-			ERROR("Couldn't select path '%s'\n", pargs->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;
-		}
-
-	}
-
-	hdf5_read(hdfile, &image, 1);
-
-	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, pargs->threshold,
-		          pargs->min_gradient, pargs->min_snr);
-		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;
-	struct sum_args *pargs;
-	char *rval;
-	struct queue_args *qargs = qp;
-
-	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->min_gradient = qargs->min_gradient;
-	pargs->min_snr = qargs->min_snr;
-	pargs->config_cmfilter = qargs->config_cmfilter;
-	pargs->config_noisefilter = qargs->config_noisefilter;
-	pargs->sum = qargs->sum;
-
-	/* Get the next filename */
-	if ( qargs->use_this_one_instead != NULL ) {
-
-		line = qargs->use_this_one_instead;
-		qargs->use_this_one_instead = NULL;
-
-	} else {
-
-		line = malloc(1024*sizeof(char));
-		rval = fgets(line, 1023, qargs->fh);
-		if ( rval == NULL ) {
-			free(pargs);
-			return NULL;
-		}
-		chomp(line);
-
-	}
-
-	if ( qargs->config_basename ) {
-		char *tmp;
-		tmp = safe_basename(line);
-		free(line);
-		line = tmp;
-	}
-
-	pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1);
-	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);
-	free(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;
-	float min_gradient = 100000.0;
-	float min_snr = 5;
-	SumMethod sum;
-	int nthreads = 1;
-	struct queue_args qargs;
-	int n_done;
-	const int chunk_size = 1000;
-	struct hdfile *hdfile;
-	struct image image;
-	char *prepare_line;
-	char prepare_filename[1024];
-	char *rval;
-	int config_basename = 0;
-	char *element = 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'},
-		{"min-gradient",       1, NULL,                4},
-		{"min-snr",            1, NULL,               11},
-		{"basename",           0, &config_basename,    1},
-		{"image",              1, NULL,               'e'},
-		{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 'e' :
-			element = strdup(optarg);
-			break;
-
-		case 4 :
-			min_gradient = strtof(optarg, NULL);
-			break;
-
-		case 11 :
-			min_snr = strtof(optarg, NULL);
-			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;
-	}
-
-	/* Get first filename and use it to set up the summed array */
-	prepare_line = malloc(1024*sizeof(char));
-	rval = fgets(prepare_line, 1023, fh);
-	if ( rval == NULL ) {
-		ERROR("Failed to get filename to prepare indexing.\n");
-		return 1;
-	}
-	chomp(prepare_line);
-	if ( config_basename ) {
-		char *tmp;
-		tmp = safe_basename(prepare_line);
-		free(prepare_line);
-		prepare_line = tmp;
-	}
-	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
-	qargs.use_this_one_instead = prepare_line;
-
-	hdfile = hdfile_open(prepare_filename);
-	if ( hdfile == NULL ) {
-		ERROR("Couldn't open '%s'\n", prepare_filename);
-		return 1;
-	}
-
-	if ( hdf5_read(hdfile, &image, 0) ) {
-		ERROR("Couldn't read '%s'\n", prepare_filename);
-		return 1;
-	}
-
-	hdfile_close(hdfile);
-	qargs.w = image.width;
-	qargs.h = image.height;
-
-	qargs.sum_method = sum;
-	qargs.threshold = threshold;
-	qargs.config_basename = config_basename;
-	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;
-	qargs.element = element;
-
-	do {
-
-		n_done = run_threads(nthreads, add_image, get_image,
-		                     (void *)&qargs, NULL, chunk_size, 0, 0, 0);
-
-		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;
-}