aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-09-21 15:31:36 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:59 +0100
commit9e5e6edd61073eb6d19e31dadb325aa88e692155 (patch)
treedbf87dc0fed11d1525fd56a5833380c403ffe826 /src
parentdbd7584014bd221bc3b9af6d3432888b6a3aba91 (diff)
Introduce 'reintegrate' for fast reintegration of data
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am7
-rw-r--r--src/Makefile.in23
-rw-r--r--src/reintegrate.c447
3 files changed, 472 insertions, 5 deletions
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;
+}