Add calibrate_detector program for working out detector geometry

3 files changed, 377 insertions, 8 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 9dbf685a..7a440ed8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,5 @@
 bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig compare_hkl \
-               powder_plot render_hkl
+               powder_plot render_hkl calibrate_detector
 
 if HAVE_GTK
 bin_PROGRAMS += hdfsee
@@ -49,4 +49,8 @@ powder_plot_LDADD = @LIBS@
 render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c
 render_hkl_LDADD = @LIBS@
 
+calibrate_detector_SOURCES = calibrate-detector.c utils.c hdf5-file.c image.c \
+                             filters.c
+calibrate_detector_LDADD = @LIBS@
+
 INCLUDES = "-I$(top_srcdir)/data"
diff --git a/src/Makefile.in b/src/Makefile.in
index 28146741..3c5e4179 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -34,7 +34,8 @@ PRE_UNINSTALL = :
 POST_UNINSTALL = :
 bin_PROGRAMS = pattern_sim$(EXEEXT) process_hkl$(EXEEXT) \
 	get_hkl$(EXEEXT) indexamajig$(EXEEXT) compare_hkl$(EXEEXT) \
-	powder_plot$(EXEEXT) render_hkl$(EXEEXT) $(am__EXEEXT_1)
+	powder_plot$(EXEEXT) render_hkl$(EXEEXT) \
+	calibrate_detector$(EXEEXT) $(am__EXEEXT_1)
 @HAVE_GTK_TRUE@am__append_1 = hdfsee
 @HAVE_OPENCL_TRUE@am__append_2 = diffraction-gpu.c cl-utils.c
 @HAVE_OPENCL_TRUE@am__append_3 = diffraction-gpu.c cl-utils.c
@@ -51,6 +52,11 @@ CONFIG_CLEAN_VPATH_FILES =
 @HAVE_GTK_TRUE@am__EXEEXT_1 = hdfsee$(EXEEXT)
 am__installdirs = "$(DESTDIR)$(bindir)"
 PROGRAMS = $(bin_PROGRAMS)
+am_calibrate_detector_OBJECTS = calibrate-detector.$(OBJEXT) \
+	utils.$(OBJEXT) hdf5-file.$(OBJEXT) image.$(OBJEXT) \
+	filters.$(OBJEXT)
+calibrate_detector_OBJECTS = $(am_calibrate_detector_OBJECTS)
+calibrate_detector_DEPENDENCIES =
 am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \
 	cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \
 	statistics.$(OBJEXT)
@@ -124,12 +130,13 @@ am__v_CCLD_0 = @echo "  CCLD  " $@;
 AM_V_GEN = $(am__v_GEN_$(V))
 am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY))
 am__v_GEN_0 = @echo "  GEN   " $@;
-SOURCES = $(compare_hkl_SOURCES) $(get_hkl_SOURCES) $(hdfsee_SOURCES) \
-	$(indexamajig_SOURCES) $(pattern_sim_SOURCES) \
-	$(powder_plot_SOURCES) $(process_hkl_SOURCES) \
-	$(render_hkl_SOURCES)
-DIST_SOURCES = $(compare_hkl_SOURCES) $(get_hkl_SOURCES) \
-	$(am__hdfsee_SOURCES_DIST) $(am__indexamajig_SOURCES_DIST) \
+SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \
+	$(get_hkl_SOURCES) $(hdfsee_SOURCES) $(indexamajig_SOURCES) \
+	$(pattern_sim_SOURCES) $(powder_plot_SOURCES) \
+	$(process_hkl_SOURCES) $(render_hkl_SOURCES)
+DIST_SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \
+	$(get_hkl_SOURCES) $(am__hdfsee_SOURCES_DIST) \
+	$(am__indexamajig_SOURCES_DIST) \
 	$(am__pattern_sim_SOURCES_DIST) $(powder_plot_SOURCES) \
 	$(process_hkl_SOURCES) $(render_hkl_SOURCES)
 ETAGS = etags
@@ -255,6 +262,10 @@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
 powder_plot_LDADD = @LIBS@
 render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c
 render_hkl_LDADD = @LIBS@
+calibrate_detector_SOURCES = calibrate-detector.c utils.c hdf5-file.c image.c \
+                             filters.c
+
+calibrate_detector_LDADD = @LIBS@
 INCLUDES = "-I$(top_srcdir)/data"
 all: all-am
 
@@ -327,6 +338,9 @@ uninstall-binPROGRAMS:
 
 clean-binPROGRAMS:
 	-test -z "$(bin_PROGRAMS)" || rm -f $(bin_PROGRAMS)
+calibrate_detector$(EXEEXT): $(calibrate_detector_OBJECTS) $(calibrate_detector_DEPENDENCIES) 
+	@rm -f calibrate_detector$(EXEEXT)
+	$(AM_V_CCLD)$(LINK) $(calibrate_detector_OBJECTS) $(calibrate_detector_LDADD) $(LIBS)
 compare_hkl$(EXEEXT): $(compare_hkl_OBJECTS) $(compare_hkl_DEPENDENCIES) 
 	@rm -f compare_hkl$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(compare_hkl_OBJECTS) $(compare_hkl_LDADD) $(LIBS)
@@ -358,6 +372,7 @@ mostlyclean-compile:
 distclean-compile:
 	-rm -f *.tab.c
 
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/calibrate-detector.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cell.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cl-utils.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/compare_hkl.Po@am__quote@
diff --git a/src/calibrate-detector.c b/src/calibrate-detector.c
new file mode 100644
index 00000000..15188a61
--- /dev/null
+++ b/src/calibrate-detector.c
@@ -0,0 +1,350 @@
+/*
+ * calibrate-detector.c
+ *
+ * Detector calibration
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#define _GNU_SOURCE 1
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <pthread.h>
+#include <sys/time.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "filters.h"
+
+
+#define MAX_THREADS (96)
+
+struct process_args
+{
+	char *filename;
+	int id;
+	int config_cmfilter;
+	int config_noisefilter;
+	float *sum;
+	int w;
+	int h;
+};
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options]\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"
+"      --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"
+"  -j <n>                  Run <n> analyses in parallel.  Default 1.\n");
+}
+
+
+static void *process_image(void *pargsv)
+{
+	struct process_args *pargs = pargsv;
+	struct hdfile *hdfile;
+	struct image image;
+	int x, y;
+
+	image.features = NULL;
+	image.data = NULL;
+	image.indexed_cell = NULL;
+	image.id = pargs->id;
+	image.filename = pargs->filename;
+	image.hits = NULL;
+	image.n_hits = 0;
+	image.det = NULL;
+
+	/* View head-on (unit cell is tilted) */
+	image.orientation.w = 1.0;
+	image.orientation.x = 0.0;
+	image.orientation.y = 0.0;
+	image.orientation.z = 0.0;
+
+	STATUS("Processing '%s'\n", pargs->filename);
+
+	hdfile = hdfile_open(pargs->filename);
+	if ( hdfile == NULL ) {
+		return NULL;
+	} else if ( hdfile_set_first_image(hdfile, "/") ) {
+		ERROR("Couldn't select path\n");
+		hdfile_close(hdfile);
+		return NULL;
+	}
+
+	hdf5_read(hdfile, &image);
+
+	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;
+	}
+
+	for ( x=0; x<image.width; x++ ) {
+	for ( y=0; y<image.height; y++ ) {
+		pargs->sum[x+pargs->w*y] += image.data[x+image.width*y];
+	}
+	}
+
+out:
+	free(image.data);
+	free(image.flags);
+	hdfile_close(hdfile);
+
+	return NULL;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *filename = NULL;
+	char *outfile = NULL;
+	FILE *fh;
+	char *rval = NULL;
+	int n_images;
+	int config_cmfilter = 0;
+	int config_noisefilter = 0;
+	char *prefix = NULL;
+	int nthreads = 1;
+	pthread_t workers[MAX_THREADS];
+	struct process_args *worker_args[MAX_THREADS];
+	int worker_active[MAX_THREADS];
+	int i;
+	const int w = 1024;  /* FIXME! */
+	const int h = 1024;  /* FIXME! */
+
+	/* 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'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:x:j:o:",
+	                        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 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 ( prefix == NULL ) {
+		prefix = strdup("");
+	}
+
+	if ( outfile == NULL ) {
+		outfile = strdup("summed.h5");
+	}
+
+	if ( (nthreads == 0) || (nthreads > MAX_THREADS) ) {
+		ERROR("Invalid number of threads.\n");
+		return 1;
+	}
+
+	/* Initialise worker arguments */
+	for ( i=0; i<nthreads; i++ ) {
+
+		worker_args[i] = malloc(sizeof(struct process_args));
+		worker_args[i]->filename = malloc(1024);
+		worker_args[i]->sum = calloc(w*h, sizeof(float));
+		worker_active[i] = 0;
+
+		worker_args[i]->w = w;
+		worker_args[i]->h = h;
+
+	}
+
+	n_images = 0;
+
+	/* Initially, fire off the full number of threads */
+	for ( i=0; i<nthreads; i++ ) {
+
+		char line[1024];
+		struct process_args *pargs;
+		int r;
+
+		pargs = worker_args[i];
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		chomp(line);
+		snprintf(pargs->filename, 1023, "%s%s", prefix, line);
+
+		n_images++;
+
+		pargs->config_cmfilter = config_cmfilter;
+		pargs->config_noisefilter = config_noisefilter;
+
+		worker_active[i] = 1;
+		r = pthread_create(&workers[i], NULL, process_image, pargs);
+		if ( r != 0 ) {
+			worker_active[i] = 0;
+			ERROR("Couldn't start thread %i\n", i);
+		}
+
+	}
+
+	/* Start new threads as old ones finish */
+	do {
+
+		int i;
+
+		for ( i=0; i<nthreads; i++ ) {
+
+			char line[1024];
+			int r;
+			struct process_result *result = NULL;
+			struct timespec t;
+			struct timeval tv;
+			struct process_args *pargs;
+
+			if ( !worker_active[i] ) continue;
+
+			pargs = worker_args[i];
+
+			gettimeofday(&tv, NULL);
+			t.tv_sec = tv.tv_sec;
+			t.tv_nsec = tv.tv_usec * 1000 + 20000;
+
+			r = pthread_timedjoin_np(workers[i], (void *)&result,
+			                         &t);
+			if ( r != 0 ) continue; /* Not ready yet */
+
+			worker_active[i] = 0;
+
+			rval = fgets(line, 1023, fh);
+			if ( rval == NULL ) break;
+			chomp(line);
+			snprintf(pargs->filename, 1023, "%s%s", prefix, line);
+
+			worker_active[i] = 1;
+			r = pthread_create(&workers[i], NULL, process_image,
+			                   pargs);
+			if ( r != 0 ) {
+				worker_active[i] = 0;
+				ERROR("Couldn't start thread %i\n", i);
+			}
+
+			n_images++;
+		}
+
+	} while ( rval != NULL );
+
+	/* Catch all remaining threads */
+	for ( i=0; i<nthreads; i++ ) {
+
+		struct process_result *result = NULL;
+
+		if ( !worker_active[i] ) goto free;
+
+		pthread_join(workers[i], (void *)&result);
+
+		worker_active[i] = 0;
+
+	free:
+		if ( worker_args[i]->filename != NULL ) {
+			free(worker_args[i]->filename);
+		}
+		free(worker_args[i]);
+
+	}
+
+	/* Sum the individual sums */
+	for ( i=1; i<nthreads; i++ ) {
+		int x, y;
+		for ( x=0; x<w; x++ ) {
+		for ( y=0; y<h; y++ ) {
+			float val = worker_args[i]->sum[x+w*y];
+			worker_args[0]->sum[x+w*y] += val;
+		}
+		}
+	}
+
+	hdf5_write(outfile, worker_args[0]->sum, w, h, H5T_NATIVE_FLOAT);
+
+	free(prefix);
+	free(outfile);
+	fclose(fh);
+
+	STATUS("There were %i images.\n", n_images);
+
+	return 0;
+}