diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Makefile.am | 11 | ||||
| -rw-r--r-- | Makefile.in | 50 | ||||
| -rw-r--r-- | src/calibrate_detector.c | 375 | ||||
| -rw-r--r-- | src/sum_stack.c | 412 |
5 files changed, 498 insertions, 351 deletions
@@ -20,5 +20,6 @@ src/cubeit src/reintegrate src/estimate_background src/check_hkl +src/sum_stack src/.dirstamp *~ diff --git a/Makefile.am b/Makefile.am index b1ab6bf0..8ecce78b 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/reintegrate \ - src/estimate_background src/check_hkl + src/estimate_background src/check_hkl src/sum_stack noinst_PROGRAMS = tests/list_check @@ -75,10 +75,13 @@ src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \ src/hdf5-file.c src/image.c src/filters.c \ src/thread-pool.c +src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \ + src/image.c src/filters.c src/peaks.c src/detector.c \ + src/cell.c src/thread-pool.c src/reflist.c + src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ - src/hdf5-file.c src/image.c src/filters.c \ - src/peaks.c src/detector.c src/cell.c \ - src/thread-pool.c src/reflist.c + src/hdf5-file.c src/image.c src/thread-pool.c \ + src/detector.c src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ diff --git a/Makefile.in b/Makefile.in index 89fe9a35..e1f68a5c 100644 --- a/Makefile.in +++ b/Makefile.in @@ -41,7 +41,7 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \ src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \ src/partialator$(EXEEXT) src/reintegrate$(EXEEXT) \ src/estimate_background$(EXEEXT) src/check_hkl$(EXEEXT) \ - $(am__EXEEXT_1) $(am__EXEEXT_2) + src/sum_stack$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2) noinst_PROGRAMS = tests/list_check$(EXEEXT) TESTS = tests/list_check$(EXEEXT) @HAVE_GTK_TRUE@am__append_1 = src/hdfsee @@ -83,9 +83,8 @@ PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS) am__dirstamp = $(am__leading_dot)dirstamp am_src_calibrate_detector_OBJECTS = src/calibrate_detector.$(OBJEXT) \ src/utils.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ - src/image.$(OBJEXT) src/filters.$(OBJEXT) src/peaks.$(OBJEXT) \ - src/detector.$(OBJEXT) src/cell.$(OBJEXT) \ - src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT) + src/image.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/detector.$(OBJEXT) src_calibrate_detector_OBJECTS = $(am_src_calibrate_detector_OBJECTS) src_calibrate_detector_LDADD = $(LDADD) src_calibrate_detector_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -224,6 +223,14 @@ am_src_render_hkl_OBJECTS = src/render_hkl.$(OBJEXT) \ src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS) src_render_hkl_LDADD = $(LDADD) src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a +am_src_sum_stack_OBJECTS = src/sum_stack.$(OBJEXT) src/utils.$(OBJEXT) \ + src/hdf5-file.$(OBJEXT) src/image.$(OBJEXT) \ + src/filters.$(OBJEXT) src/peaks.$(OBJEXT) \ + src/detector.$(OBJEXT) src/cell.$(OBJEXT) \ + src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT) +src_sum_stack_OBJECTS = $(am_src_sum_stack_OBJECTS) +src_sum_stack_LDADD = $(LDADD) +src_sum_stack_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_tests_list_check_OBJECTS = tests/list_check.$(OBJEXT) \ src/reflist.$(OBJEXT) src/thread-pool.$(OBJEXT) \ src/utils.$(OBJEXT) @@ -257,7 +264,7 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \ $(src_partialator_SOURCES) $(src_pattern_sim_SOURCES) \ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \ $(src_reintegrate_SOURCES) $(src_render_hkl_SOURCES) \ - $(tests_list_check_SOURCES) + $(src_sum_stack_SOURCES) $(tests_list_check_SOURCES) DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \ $(am__src_cubeit_SOURCES_DIST) \ @@ -266,7 +273,8 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(am__src_indexamajig_SOURCES_DIST) $(src_partialator_SOURCES) \ $(am__src_pattern_sim_SOURCES_DIST) $(src_powder_plot_SOURCES) \ $(src_process_hkl_SOURCES) $(src_reintegrate_SOURCES) \ - $(src_render_hkl_SOURCES) $(tests_list_check_SOURCES) + $(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \ + $(tests_list_check_SOURCES) RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \ html-recursive info-recursive install-data-recursive \ install-dvi-recursive install-exec-recursive \ @@ -614,10 +622,13 @@ src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \ src/hdf5-file.c src/image.c src/filters.c \ src/thread-pool.c +src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \ + src/image.c src/filters.c src/peaks.c src/detector.c \ + src/cell.c src/thread-pool.c src/reflist.c + src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ - src/hdf5-file.c src/image.c src/filters.c \ - src/peaks.c src/detector.c src/cell.c \ - src/thread-pool.c src/reflist.c + src/hdf5-file.c src/image.c src/thread-pool.c \ + src/detector.c src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ @@ -782,15 +793,9 @@ src/utils.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/hdf5-file.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/image.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/filters.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/peaks.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/detector.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/thread-pool.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) -src/reflist.$(OBJEXT): src/$(am__dirstamp) \ +src/detector.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibrate_detector_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/calibrate_detector$(EXEEXT) @@ -798,6 +803,7 @@ src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibra src/check_hkl.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/sfac.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) +src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/reflections.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/statistics.$(OBJEXT): src/$(am__dirstamp) \ @@ -816,6 +822,8 @@ src/cubeit.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/render.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/filters.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) src/stream.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp) @@ -842,6 +850,7 @@ src/hdfsee$(EXEEXT): $(src_hdfsee_OBJECTS) $(src_hdfsee_DEPENDENCIES) src/$(am__ $(AM_V_CCLD)$(LINK) $(src_hdfsee_OBJECTS) $(src_hdfsee_LDADD) $(LIBS) src/indexamajig.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/peaks.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/index.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) src/diffraction.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) @@ -852,6 +861,8 @@ src/templates.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/geometry.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/reflist.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) src/diffraction-gpu.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/cl-utils.$(OBJEXT): src/$(am__dirstamp) \ @@ -895,6 +906,11 @@ 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) @@ -947,6 +963,7 @@ mostlyclean-compile: -rm -f src/sfac.$(OBJEXT) -rm -f src/statistics.$(OBJEXT) -rm -f src/stream.$(OBJEXT) + -rm -f src/sum_stack.$(OBJEXT) -rm -f src/symmetry.$(OBJEXT) -rm -f src/templates.$(OBJEXT) -rm -f src/thread-pool.$(OBJEXT) @@ -994,6 +1011,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sfac.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/statistics.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/stream.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/symmetry.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/templates.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/thread-pool.Po@am__quote@ diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c index 417ee36e..25e9850c 100644 --- a/src/calibrate_detector.c +++ b/src/calibrate_detector.c @@ -1,9 +1,9 @@ /* - * calibrate-detector.c + * calibrate_detector.c * - * Detector calibration + * Attempt to refine detector geometry * - * (c) 2006-2010 Thomas White <taw@physics.org> + * (c) 2006-2011 Thomas White <taw@physics.org> * * Part of CrystFEL - crystallography with a FEL * @@ -20,306 +20,61 @@ #include <string.h> #include <unistd.h> #include <getopt.h> -#include <pthread.h> #include "utils.h" +#include "image.h" +#include "detector.h" +#include "index.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; -}; - - -struct queue_args -{ - FILE *fh; - char *prefix; - int config_cmfilter; - int config_noisefilter; - double *sum; - int w; - int h; - SumMethod sum_method; - double threshold; -}; static void show_help(const char *s) { - printf("Syntax: %s [options]\n\n", s); + printf("Syntax: %s [options] <file.h5>\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" -" -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.\n" +"Plot a powder pattern as a 1D graph using the detector geometry.\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) -{ - int x, y, i; - int w = image->width; - const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS; - - /* FIXME: Get threshold value from command line */ - search_peaks(image, 800.0, 100000.0); - - for ( i=0; i<image_feature_count(image->features); i++ ) { - - struct imagefeature *f = image_get_feature(image->features, i); - int xp, yp; - - /* This is not an error. */ - if ( f == NULL ) continue; - - xp = f->x; - yp = f->y; - - for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) { - for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) { - - /* Circular mask */ - if ( x*x + y*y > lim ) continue; - - if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; - if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; - - float val = image->data[(x+xp)+w*(y+yp)]; - sum[(x+xp)+w*(y+yp)] += 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; - } else if ( hdfile_set_first_image(hdfile, "/") ) { - ERROR("Couldn't select path\n"); - hdfile_close(hdfile); - return; - } - - /* FIXME: Nominal photon energy */ - hdf5_read(hdfile, &image, 1, 2000.0); - - 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); - 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[1024]; - struct sum_args *pargs; - char *rval; - struct queue_args *qargs = qp; - - /* Get the next filename */ - rval = fgets(line, 1023, qargs->fh); - if ( rval == NULL ) return NULL; - - 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->config_cmfilter = qargs->config_cmfilter; - pargs->config_noisefilter = qargs->config_noisefilter; - pargs->sum = qargs->sum; - - chomp(line); - pargs->filename = malloc(strlen(qargs->prefix) + strlen(line) + 1); - snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line); - - return pargs; +" -h, --help Display this help message.\n" +" -g. --geometry=<file> Get detector geometry from file.\n" +" -i, --input=<file> Input filename.\n" +"\n"); } int main(int argc, char *argv[]) { int c; + struct image image; + int x, y; + struct hdfile *hdfile; 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; - SumMethod sum; - int nthreads = 1; - struct queue_args qargs; - int n_done; - const int chunk_size = 1000; + char *geometry = 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'}, + {"geometry", 1, NULL, 'g'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:x:j:o:s:p:t:", - longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:g:", 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); + case 0 : break; - case 't' : - threshold = atof(optarg); + case 'i' : + filename = strdup(optarg); break; - case 0 : + case 'g' : + geometry = strdup(optarg); break; default : @@ -329,84 +84,42 @@ int main(int argc, char *argv[]) } 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); + ERROR("You must specify the input filename with -i\n"); 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); + if ( geometry == NULL ) { + ERROR("You need to specify a geometry file with --geometry\n"); 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"); + image.det = get_detector_geometry(geometry); + if ( image.det == NULL ) { + ERROR("Failed to read detector geometry from '%s'\n", geometry); return 1; } + free(geometry); - qargs.w = 1024; /* FIXME! */ - qargs.h = 1024; /* FIXME! */ - qargs.sum_method = sum; - qargs.threshold = threshold; - 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; - - do { - - n_done = run_threads(nthreads, add_image, get_image, - (void *)&qargs, NULL, chunk_size); + hdfile = hdfile_open(filename); + hdfile_set_image(hdfile, "/data/data"); + hdf5_read(hdfile, &image, 1, 2000.0); - n_images += n_done; + for ( x=0; x<image.width; x++ ) { + for ( y=0; y<image.height; y++ ) { - /* 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); - } + double q; + int intensity; + struct rvec r; - } while ( n_done == chunk_size ); + r = get_q(&image, x, y, 1, NULL, 1.0/image.lambda); + q = modulus(r.u, r.v, r.w); - /* Write the final output */ - hdf5_write(outfile, qargs.sum, qargs.w, qargs.h, H5T_NATIVE_DOUBLE); + intensity = image.data[x + image.width*y]; - free(qargs.sum); - free(prefix); - free(outfile); - if ( intermediate != NULL ) free(intermediate); - fclose(fh); + printf("%5i\t%5i\t%7.3f\t%7i\n", x, y, q/1.0e9, intensity); - STATUS("There were %i images.\n", n_images); + } + } return 0; } diff --git a/src/sum_stack.c b/src/sum_stack.c new file mode 100644 index 00000000..8bd5197e --- /dev/null +++ b/src/sum_stack.c @@ -0,0 +1,412 @@ +/* + * 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; +}; + + +struct queue_args +{ + FILE *fh; + char *prefix; + int config_cmfilter; + int config_noisefilter; + double *sum; + int w; + int h; + SumMethod sum_method; + double threshold; +}; + + +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.\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) +{ + int x, y, i; + int w = image->width; + const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS; + + /* FIXME: Get threshold value from command line */ + search_peaks(image, 800.0, 100000.0); + + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f = image_get_feature(image->features, i); + int xp, yp; + + /* This is not an error. */ + if ( f == NULL ) continue; + + xp = f->x; + yp = f->y; + + for ( x=-INTEGRATION_RADIUS; x<+INTEGRATION_RADIUS; x++ ) { + for ( y=-INTEGRATION_RADIUS; y<+INTEGRATION_RADIUS; y++ ) { + + /* Circular mask */ + if ( x*x + y*y > lim ) continue; + + if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; + if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; + + float val = image->data[(x+xp)+w*(y+yp)]; + sum[(x+xp)+w*(y+yp)] += 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; + } else if ( hdfile_set_first_image(hdfile, "/") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return; + } + + /* FIXME: Nominal photon energy */ + hdf5_read(hdfile, &image, 1, 2000.0); + + 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); + 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[1024]; + struct sum_args *pargs; + char *rval; + struct queue_args *qargs = qp; + + /* Get the next filename */ + rval = fgets(line, 1023, qargs->fh); + if ( rval == NULL ) return NULL; + + 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->config_cmfilter = qargs->config_cmfilter; + pargs->config_noisefilter = qargs->config_noisefilter; + pargs->sum = qargs->sum; + + chomp(line); + pargs->filename = malloc(strlen(qargs->prefix) + strlen(line) + 1); + snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, 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; + SumMethod sum; + int nthreads = 1; + struct queue_args qargs; + int n_done; + const int chunk_size = 1000; + + /* 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'}, + {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 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; + } + + qargs.w = 1024; /* FIXME! */ + qargs.h = 1024; /* FIXME! */ + qargs.sum_method = sum; + qargs.threshold = threshold; + 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; + + do { + + n_done = run_threads(nthreads, add_image, get_image, + (void *)&qargs, NULL, chunk_size); + + 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; +} |