diff options
| -rw-r--r-- | Makefile.am | 4 | ||||
| -rw-r--r-- | Makefile.in | 20 | ||||
| -rw-r--r-- | README | 5 | ||||
| -rw-r--r-- | src/sum_stack.c | 523 |
4 files changed, 4 insertions, 548 deletions
diff --git a/Makefile.am b/Makefile.am index 53321656..598b8ee0 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/check_hkl src/sum_stack src/partial_sim + src/check_hkl src/partial_sim noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_gradient_check tests/symmetry_check @@ -58,8 +58,6 @@ src_powder_plot_SOURCES = src/powder_plot.c src_render_hkl_SOURCES = src/render_hkl.c src/povray.c -src_sum_stack_SOURCES = src/sum_stack.c - src_calibrate_detector_SOURCES = src/calibrate_detector.c src_partialator_SOURCES = src/partialator.c src/post-refinement.c \ diff --git a/Makefile.in b/Makefile.in index bffe2b0d..77193584 100644 --- a/Makefile.in +++ b/Makefile.in @@ -40,8 +40,7 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \ src/compare_hkl$(EXEEXT) src/powder_plot$(EXEEXT) \ src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \ src/partialator$(EXEEXT) src/check_hkl$(EXEEXT) \ - src/sum_stack$(EXEEXT) src/partial_sim$(EXEEXT) \ - $(am__EXEEXT_1) + src/partial_sim$(EXEEXT) $(am__EXEEXT_1) noinst_PROGRAMS = tests/list_check$(EXEEXT) \ tests/integration_check$(EXEEXT) \ tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \ @@ -162,11 +161,6 @@ src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS) src_render_hkl_LDADD = $(LDADD) src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ $(top_builddir)/libcrystfel/libcrystfel.la -am_src_sum_stack_OBJECTS = src/sum_stack.$(OBJEXT) -src_sum_stack_OBJECTS = $(am_src_sum_stack_OBJECTS) -src_sum_stack_LDADD = $(LDADD) -src_sum_stack_DEPENDENCIES = $(top_builddir)/lib/libgnu.a \ - $(top_builddir)/libcrystfel/libcrystfel.la am__tests_gpu_sim_check_SOURCES_DIST = tests/gpu_sim_check.c @HAVE_OPENCL_TRUE@am_tests_gpu_sim_check_OBJECTS = \ @HAVE_OPENCL_TRUE@ tests/gpu_sim_check.$(OBJEXT) @@ -231,7 +225,7 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \ $(src_partial_sim_SOURCES) $(src_partialator_SOURCES) \ $(src_pattern_sim_SOURCES) $(src_powder_plot_SOURCES) \ $(src_process_hkl_SOURCES) $(src_render_hkl_SOURCES) \ - $(src_sum_stack_SOURCES) $(tests_gpu_sim_check_SOURCES) \ + $(tests_gpu_sim_check_SOURCES) \ $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \ $(tests_pr_gradient_check_SOURCES) \ $(tests_symmetry_check_SOURCES) @@ -241,7 +235,7 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(src_indexamajig_SOURCES) $(src_partial_sim_SOURCES) \ $(am__src_partialator_SOURCES_DIST) $(src_pattern_sim_SOURCES) \ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \ - $(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \ + $(src_render_hkl_SOURCES) \ $(am__tests_gpu_sim_check_SOURCES_DIST) \ $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \ $(tests_pr_gradient_check_SOURCES) \ @@ -581,7 +575,6 @@ src_compare_hkl_SOURCES = src/compare_hkl.c src_check_hkl_SOURCES = src/check_hkl.c src_powder_plot_SOURCES = src/powder_plot.c src_render_hkl_SOURCES = src/render_hkl.c src/povray.c -src_sum_stack_SOURCES = src/sum_stack.c src_calibrate_detector_SOURCES = src/calibrate_detector.c src_partialator_SOURCES = src/partialator.c src/post-refinement.c \ src/hrs-scaling.c $(am__append_4) @@ -802,11 +795,6 @@ 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) @@ -858,7 +846,6 @@ mostlyclean-compile: -rm -f src/process_hkl.$(OBJEXT) -rm -f src/render_hkl.$(OBJEXT) -rm -f src/scaling-report.$(OBJEXT) - -rm -f src/sum_stack.$(OBJEXT) -rm -f tests/gpu_sim_check.$(OBJEXT) -rm -f tests/integration_check.$(OBJEXT) -rm -f tests/list_check.$(OBJEXT) @@ -885,7 +872,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/scaling-report.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/sum_stack.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/gpu_sim_check.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@ @@ -52,11 +52,6 @@ In addition, there is also: - check_hkl, for determining things like completeness. - - sum_stack, for summing patterns after peak detection to use - for calibrating your detector. - - - cubeit, for summing patterns in 3D to inspect the peak profiles. - - calibrate_detector, for refining detector geometry. - partial_sim, for calculating partial reflection intensities. diff --git a/src/sum_stack.c b/src/sum_stack.c deleted file mode 100644 index a6e7599a..00000000 --- a/src/sum_stack.c +++ /dev/null @@ -1,523 +0,0 @@ -/* - * 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; - double min_gradient; - double min_snr; - char *element; -}; - - -struct queue_args -{ - FILE *fh; - int config_cmfilter; - int config_noisefilter; - double *sum; - int w; - int h; - SumMethod sum_method; - double threshold; - double min_gradient; - double min_snr; - char *element; - - char *use_this_one_instead; - char *prefix; - int config_basename; -}; - - -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. Default: 400 adu\n" -" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n" -" Default: 100,000.\n" -" --min-snr=<n> Minimum signal-to-noise ration for peaks.\n" -" Default: 5.\n" -" -e, --image=<element> Use this image from the HDF5 file.\n" -" Example: /data/data0.\n" -" Default: The first one found.\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, double threshold, - double min_gradient, double min_snr) -{ - int i; - int w = image->width; - int h = image->height; - const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS; - - search_peaks(image, threshold, min_gradient, min_snr); - - for ( i=0; i<image_feature_count(image->features); i++ ) { - - struct imagefeature *f = image_get_feature(image->features, i); - int pfs, pss; - int fs, ss; - - /* This is not an error. */ - if ( f == NULL ) continue; - - pfs = f->fs; - pss = f->ss; - - for ( fs=-INTEGRATION_RADIUS; fs<+INTEGRATION_RADIUS; fs++ ) { - for ( ss=-INTEGRATION_RADIUS; ss<+INTEGRATION_RADIUS; ss++ ) { - - /* Circular mask */ - if ( fs*fs + ss*ss > lim ) continue; - - if ( ((fs+pfs)>=w) || ((fs+pfs)<0) ) continue; - if ( ((ss+pss)>=h) || ((ss+pss)<0) ) continue; - - float val = image->data[(fs+pfs)+w*(ss+pss)]; - sum[(fs+pfs)+w*(ss+pss)] += 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; - - if ( pargs->element != NULL ) { - int r; - r = hdfile_set_image(hdfile, pargs->element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", pargs->element); - hdfile_close(hdfile); - return; - } - } else { - int r; - r = hdfile_set_first_image(hdfile, "/"); - if ( r ) { - ERROR("Couldn't select first path\n"); - hdfile_close(hdfile); - return; - } - - } - - hdf5_read(hdfile, &image, 1); - - 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, pargs->threshold, - pargs->min_gradient, pargs->min_snr); - 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; - struct sum_args *pargs; - char *rval; - struct queue_args *qargs = qp; - - 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->min_gradient = qargs->min_gradient; - pargs->min_snr = qargs->min_snr; - pargs->config_cmfilter = qargs->config_cmfilter; - pargs->config_noisefilter = qargs->config_noisefilter; - pargs->sum = qargs->sum; - - /* Get the next filename */ - if ( qargs->use_this_one_instead != NULL ) { - - line = qargs->use_this_one_instead; - qargs->use_this_one_instead = NULL; - - } else { - - line = malloc(1024*sizeof(char)); - rval = fgets(line, 1023, qargs->fh); - if ( rval == NULL ) { - free(pargs); - return NULL; - } - chomp(line); - - } - - if ( qargs->config_basename ) { - char *tmp; - tmp = safe_basename(line); - free(line); - line = tmp; - } - - pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1); - snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line); - free(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; - float min_gradient = 100000.0; - float min_snr = 5; - SumMethod sum; - int nthreads = 1; - struct queue_args qargs; - int n_done; - const int chunk_size = 1000; - struct hdfile *hdfile; - struct image image; - char *prepare_line; - char prepare_filename[1024]; - char *rval; - int config_basename = 0; - char *element = 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'}, - {"min-gradient", 1, NULL, 4}, - {"min-snr", 1, NULL, 11}, - {"basename", 0, &config_basename, 1}, - {"image", 1, NULL, 'e'}, - {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 'e' : - element = strdup(optarg); - break; - - case 4 : - min_gradient = strtof(optarg, NULL); - break; - - case 11 : - min_snr = strtof(optarg, NULL); - 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; - } - - /* Get first filename and use it to set up the summed array */ - prepare_line = malloc(1024*sizeof(char)); - rval = fgets(prepare_line, 1023, fh); - if ( rval == NULL ) { - ERROR("Failed to get filename to prepare indexing.\n"); - return 1; - } - chomp(prepare_line); - if ( config_basename ) { - char *tmp; - tmp = safe_basename(prepare_line); - free(prepare_line); - prepare_line = tmp; - } - snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line); - qargs.use_this_one_instead = prepare_line; - - hdfile = hdfile_open(prepare_filename); - if ( hdfile == NULL ) { - ERROR("Couldn't open '%s'\n", prepare_filename); - return 1; - } - - if ( hdf5_read(hdfile, &image, 0) ) { - ERROR("Couldn't read '%s'\n", prepare_filename); - return 1; - } - - hdfile_close(hdfile); - qargs.w = image.width; - qargs.h = image.height; - - qargs.sum_method = sum; - qargs.threshold = threshold; - qargs.config_basename = config_basename; - 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; - qargs.element = element; - - do { - - n_done = run_threads(nthreads, add_image, get_image, - (void *)&qargs, NULL, chunk_size, 0, 0, 0); - - 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; -} |