From 652043a185e214754dcc42c864a748bf91255282 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 12 Jun 2010 09:12:34 -0700 Subject: Add calibrate_detector program for working out detector geometry --- src/Makefile.am | 6 +- src/Makefile.in | 29 +++- src/calibrate-detector.c | 350 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 377 insertions(+), 8 deletions(-) create mode 100644 src/calibrate-detector.c (limited to 'src') 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 + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#define _GNU_SOURCE 1 +#include +#include +#include +#include +#include +#include +#include +#include + +#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= Specify file containing list of images to process.\n" +" '-' means stdin, which is the default.\n" +" -o, --output= 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 Run 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; xsum[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; ifilename = 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; ifilename, 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; ifilename, 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; ifilename != NULL ) { + free(worker_args[i]->filename); + } + free(worker_args[i]); + + } + + /* Sum the individual sums */ + for ( i=1; isum[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; +} -- cgit v1.2.3