Introduce 'reintegrate' for fast reintegration of data

4 files changed, 473 insertions, 5 deletions

diff --git a/.gitignore b/.gitignore
index 2b010e63..7bd9572d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -17,4 +17,5 @@ src/render_hkl
 src/calibrate_detector
 src/facetron
 src/cubeit
+src/reintegrate
 *~
diff --git a/src/Makefile.am b/src/Makefile.am
index cbee354b..fa6a6f41 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,6 @@
 bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig compare_hkl \
-               powder_plot render_hkl calibrate_detector facetron cubeit
+               powder_plot render_hkl calibrate_detector facetron cubeit \
+               reintegrate
 
 if HAVE_GTK
 bin_PROGRAMS += hdfsee
@@ -64,4 +65,8 @@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c diffraction.c \
                  sfac.c render.c filters.c image.c symmetry.c
 cubeit_LDADD = @LIBS@
 
+reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \
+                      peaks.c image.c diffraction.c sfac.c geometry.c
+reintegrate_LDADD = @LIBS@
+
 INCLUDES = "-I$(top_srcdir)/data"
diff --git a/src/Makefile.in b/src/Makefile.in
index 1fce7c2d..81a5c34e 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -36,7 +36,7 @@ bin_PROGRAMS = pattern_sim$(EXEEXT) process_hkl$(EXEEXT) \
 	get_hkl$(EXEEXT) indexamajig$(EXEEXT) compare_hkl$(EXEEXT) \
 	powder_plot$(EXEEXT) render_hkl$(EXEEXT) \
 	calibrate_detector$(EXEEXT) facetron$(EXEEXT) cubeit$(EXEEXT) \
-	$(am__EXEEXT_1)
+	reintegrate$(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
@@ -121,6 +121,12 @@ am_process_hkl_OBJECTS = process_hkl.$(OBJEXT) sfac.$(OBJEXT) \
 	reflections.$(OBJEXT) symmetry.$(OBJEXT)
 process_hkl_OBJECTS = $(am_process_hkl_OBJECTS)
 process_hkl_DEPENDENCIES =
+am_reintegrate_OBJECTS = reintegrate.$(OBJEXT) cell.$(OBJEXT) \
+	hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
+	peaks.$(OBJEXT) image.$(OBJEXT) diffraction.$(OBJEXT) \
+	sfac.$(OBJEXT) geometry.$(OBJEXT)
+reintegrate_OBJECTS = $(am_reintegrate_OBJECTS)
+reintegrate_DEPENDENCIES =
 am_render_hkl_OBJECTS = render_hkl.$(OBJEXT) cell.$(OBJEXT) \
 	reflections.$(OBJEXT) utils.$(OBJEXT) povray.$(OBJEXT) \
 	symmetry.$(OBJEXT) render.$(OBJEXT) hdf5-file.$(OBJEXT) \
@@ -151,12 +157,14 @@ SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \
 	$(cubeit_SOURCES) $(facetron_SOURCES) $(get_hkl_SOURCES) \
 	$(hdfsee_SOURCES) $(indexamajig_SOURCES) \
 	$(pattern_sim_SOURCES) $(powder_plot_SOURCES) \
-	$(process_hkl_SOURCES) $(render_hkl_SOURCES)
+	$(process_hkl_SOURCES) $(reintegrate_SOURCES) \
+	$(render_hkl_SOURCES)
 DIST_SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \
 	$(cubeit_SOURCES) $(facetron_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)
+	$(process_hkl_SOURCES) $(reintegrate_SOURCES) \
+	$(render_hkl_SOURCES)
 ETAGS = etags
 CTAGS = ctags
 DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
@@ -201,7 +209,6 @@ PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@
 PACKAGE_NAME = @PACKAGE_NAME@
 PACKAGE_STRING = @PACKAGE_STRING@
 PACKAGE_TARNAME = @PACKAGE_TARNAME@
-PACKAGE_URL = @PACKAGE_URL@
 PACKAGE_VERSION = @PACKAGE_VERSION@
 PATH_SEPARATOR = @PATH_SEPARATOR@
 PKG_CONFIG = @PKG_CONFIG@
@@ -297,6 +304,10 @@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c diffraction.c \
                  sfac.c render.c filters.c image.c symmetry.c
 
 cubeit_LDADD = @LIBS@
+reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \
+                      peaks.c image.c diffraction.c sfac.c geometry.c
+
+reintegrate_LDADD = @LIBS@
 INCLUDES = "-I$(top_srcdir)/data"
 all: all-am
 
@@ -399,6 +410,9 @@ powder_plot$(EXEEXT): $(powder_plot_OBJECTS) $(powder_plot_DEPENDENCIES)
 process_hkl$(EXEEXT): $(process_hkl_OBJECTS) $(process_hkl_DEPENDENCIES) 
 	@rm -f process_hkl$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(process_hkl_OBJECTS) $(process_hkl_LDADD) $(LIBS)
+reintegrate$(EXEEXT): $(reintegrate_OBJECTS) $(reintegrate_DEPENDENCIES) 
+	@rm -f reintegrate$(EXEEXT)
+	$(AM_V_CCLD)$(LINK) $(reintegrate_OBJECTS) $(reintegrate_LDADD) $(LIBS)
 render_hkl$(EXEEXT): $(render_hkl_OBJECTS) $(render_hkl_DEPENDENCIES) 
 	@rm -f render_hkl$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(render_hkl_OBJECTS) $(render_hkl_LDADD) $(LIBS)
@@ -434,6 +448,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/powder_plot.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/process_hkl.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/reflections.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/reintegrate.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/render.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/render_hkl.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sfac.Po@am__quote@
diff --git a/src/reintegrate.c b/src/reintegrate.c
new file mode 100644
index 00000000..3cb94730
--- /dev/null
+++ b/src/reintegrate.c
@@ -0,0 +1,447 @@
+/*
+ * reintegrate.c
+ *
+ * Like "indexamajig", but skip the indexing step
+ *
+ * (c) 2006-2010 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 <sys/time.h>
+#include <assert.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "symmetry.h"
+#include "peaks.h"
+
+
+#define MAX_THREADS (256)
+
+struct process_args
+{
+	char *filename;
+	int id;
+
+	/* Thread control */
+	pthread_mutex_t control_mutex;  /* Protects the scary stuff below */
+	int start;
+	int finish;
+	int done;
+
+	UnitCell *cell;
+	struct detector *det;
+	pthread_mutex_t *output_mutex;  /* Protects 'stdout' */
+};
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options]\n\n", s);
+	printf(
+"Like 'indexamajig', but skip the indexing step.\n"
+"\n"
+"  -h, --help                 Display this help message.\n"
+"\n"
+"  -i, --input=<filename>     Specify the name of the input 'stream'.\n"
+"                              (must be a file, not e.g. stdin)\n"
+"  -g. --geometry=<file>      Get detector geometry from file.\n"
+"  -x, --prefix=<p>           Prefix filenames from input file with <p>.\n"
+"      --basename             Remove the directory parts of the filenames.\n"
+"      --no-check-prefix      Don't attempt to correct the --prefix.\n"
+"  -j <n>                     Run <n> analyses in parallel.\n");
+}
+
+
+static void process_image(struct process_args *pargs)
+{
+	struct hdfile *hdfile;
+	struct image image;
+
+	image.features = NULL;
+	image.data = NULL;
+	image.flags = NULL;
+	image.indexed_cell = NULL;
+	image.id = pargs->id;
+	image.filename = pargs->filename;
+	image.hits = NULL;
+	image.n_hits = 0;
+	image.det = pargs->det;
+
+	/* 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;
+	} else if ( hdfile_set_first_image(hdfile, "/") ) {
+		ERROR("Couldn't select path\n");
+		hdfile_close(hdfile);
+		return;
+	}
+
+	hdf5_read(hdfile, &image, 1);
+
+	output_intensities(&image, pargs->cell,
+		                   pargs->output_mutex, 1,
+		                   1, 0);
+
+	free(image.data);
+	cell_free(pargs->cell);
+	if ( image.flags != NULL ) free(image.flags);
+	hdfile_close(hdfile);
+}
+
+
+static void *worker_thread(void *pargsv)
+{
+	struct process_args *pargs = pargsv;
+	int finish;
+
+	do {
+
+		int wakeup;
+
+		process_image(pargs);
+
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->done = 1;
+		pargs->start = 0;
+		pthread_mutex_unlock(&pargs->control_mutex);
+
+		/* Go to sleep until told to exit or process next image */
+		do {
+
+			pthread_mutex_lock(&pargs->control_mutex);
+			/* Either of these can result in the thread waking up */
+			wakeup = pargs->start || pargs->finish;
+			finish = pargs->finish;
+			pthread_mutex_unlock(&pargs->control_mutex);
+			usleep(20000);
+
+		} while ( !wakeup );
+
+	} while ( !pargs->finish );
+
+	return NULL;
+}
+
+
+static UnitCell *read_orientation_matrix(FILE *fh)
+{
+	float u, v, w;
+	struct rvec as, bs, cs;
+	UnitCell *cell;
+	char line[1024];
+
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read a-star\n");
+		return NULL;
+	}
+	as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read b-star\n");
+		return NULL;
+	}
+	bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
+	if ( fgets(line, 1023, fh) == NULL ) return NULL;
+	if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) {
+		ERROR("Couldn't read c-star\n");
+		return NULL;
+	}
+	cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
+	cell = cell_new_from_axes(as, bs, cs);
+
+	return cell;
+}
+
+
+static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
+{
+	char line[1024];
+	char *rval = NULL;
+
+	do {
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+
+		chomp(line);
+
+		if ( strncmp(line, "Reflections from indexing", 25) != 0 ) {
+			continue;
+		}
+
+		*filename = strdup(line+29);
+
+		/* Skip two lines (while checking for errors) */
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+
+		*cell = read_orientation_matrix(fh);
+		if ( *cell == NULL ) {
+			STATUS("Got filename but no cell for %s\n", *filename);
+			continue;
+		}
+
+		return 0;
+
+	} while ( rval != NULL );
+
+	return 1;
+}
+
+
+static void integrate_all(int nthreads, struct detector *det, FILE *fh,
+                          int config_basename, const char *prefix)
+{
+	pthread_t workers[MAX_THREADS];
+	struct process_args *worker_args[MAX_THREADS];
+	int worker_active[MAX_THREADS];
+	int i;
+	int rval;
+	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+	/* Initialise worker arguments */
+	for ( i=0; i<nthreads; i++ ) {
+
+		worker_args[i] = malloc(sizeof(struct process_args));
+		worker_args[i]->filename = malloc(1024);
+		worker_active[i] = 0;
+		worker_args[i]->det = det;
+		pthread_mutex_init(&worker_args[i]->control_mutex, NULL);
+		worker_args[i]->output_mutex = &output_mutex;
+
+	}
+
+	/* Start threads off */
+	for ( i=0; i<nthreads; i++ ) {
+
+		struct process_args *pargs;
+		int r;
+		int rval;
+		char *filename;
+		UnitCell *cell;
+
+		pargs = worker_args[i];
+
+		/* Get the next filename */
+		rval = find_chunk(fh, &cell, &filename);
+		if ( rval == 1 ) break;
+		if ( config_basename ) {
+			char *tmp;
+			tmp = basename(filename);
+			free(filename);
+			filename = tmp;
+		}
+		snprintf(pargs->filename, 1023, "%s%s",
+		         prefix, filename);
+		pargs->cell = cell;
+		free(filename);
+
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->done = 0;
+		pargs->start = 1;
+		pargs->finish = 0;
+		pthread_mutex_unlock(&pargs->control_mutex);
+
+		worker_active[i] = 1;
+		r = pthread_create(&workers[i], NULL, worker_thread, pargs);
+		if ( r != 0 ) {
+			worker_active[i] = 0;
+			ERROR("Couldn't start thread %i\n", i);
+		}
+
+	}
+
+	/* Keep threads busy until the end of the data */
+	do {
+
+		int i;
+		rval = 0;
+
+		for ( i=0; i<nthreads; i++ ) {
+
+			struct process_args *pargs;
+			int done;
+			char *filename;
+			UnitCell *cell;
+
+			/* Spend time working, not managing threads */
+			usleep(100000);
+
+			/* Are we using this thread record at all? */
+			if ( !worker_active[i] ) continue;
+
+			/* Has the thread finished yet? */
+			pargs = worker_args[i];
+			pthread_mutex_lock(&pargs->control_mutex);
+			done = pargs->done;
+			pthread_mutex_unlock(&pargs->control_mutex);
+			if ( !done ) continue;
+
+			/* Get the next filename */
+			rval = find_chunk(fh, &cell, &filename);
+			if ( rval == 1 ) break;
+			if ( config_basename ) {
+				char *tmp;
+				tmp = basename(filename);
+				free(filename);
+				filename = tmp;
+			}
+			snprintf(pargs->filename, 1023, "%s%s",
+			         prefix, filename);
+			pargs->cell = cell;
+			free(filename);
+
+			/* Wake the thread up ... */
+			pthread_mutex_lock(&pargs->control_mutex);
+			pargs->done = 0;
+			pargs->start = 1;
+			pthread_mutex_unlock(&pargs->control_mutex);
+
+		}
+
+	} while ( rval == 0 );
+
+	/* Join threads */
+	for ( i=0; i<nthreads; i++ ) {
+
+		if ( !worker_active[i] ) goto free;
+
+		/* Tell the thread to exit */
+		struct process_args *pargs = worker_args[i];
+		pthread_mutex_lock(&pargs->control_mutex);
+		pargs->finish = 1;
+		pthread_mutex_unlock(&pargs->control_mutex);
+
+		/* Wait for it to join */
+		pthread_join(workers[i], NULL);
+
+	free:
+		if ( worker_args[i]->filename != NULL ) {
+			free(worker_args[i]->filename);
+		}
+
+	}
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *infile = NULL;
+	char *geomfile = NULL;
+	FILE *fh;
+	char *prefix = NULL;
+	int nthreads = 1;
+	int config_basename = 0;
+	int config_checkprefix = 1;
+	struct detector *det;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"input",              1, NULL,               'i'},
+		{"geometry",           1, NULL,               'g'},
+		{"prefix",             1, NULL,               'x'},
+		{"basename",           0, &config_basename,    1},
+		{"no-check-prefix",    0, &config_checkprefix, 0},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:g:x:j:",
+	                        longopts, NULL)) != -1)
+	{
+
+		switch (c) {
+		case 'h' :
+			show_help(argv[0]);
+			return 0;
+
+		case 'i' :
+			infile = strdup(optarg);
+			break;
+
+		case 'g' :
+			geomfile = strdup(optarg);
+			break;
+
+		case 'x' :
+			prefix = strdup(optarg);
+			break;
+
+		case 'j' :
+			nthreads = atoi(optarg);
+			break;
+
+		case 0 :
+			break;
+
+		default :
+			return 1;
+		}
+
+	}
+
+	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);
+		return 1;
+	}
+	free(infile);
+
+	if ( prefix == NULL ) {
+		prefix = strdup("");
+	} else {
+		if ( config_checkprefix ) {
+			prefix = check_prefix(prefix);
+		}
+	}
+
+	det = get_detector_geometry(geomfile);
+	if ( det == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n", geomfile);
+		return 1;
+	}
+	free(geomfile);
+
+	rewind(fh);
+	integrate_all(nthreads, det, fh, config_basename, prefix);
+
+	fclose(fh);
+	free(prefix);
+
+	return 0;
+}