aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore1
-rw-r--r--Makefile.am11
-rw-r--r--Makefile.in50
-rw-r--r--src/calibrate_detector.c375
-rw-r--r--src/sum_stack.c412
5 files changed, 498 insertions, 351 deletions
diff --git a/.gitignore b/.gitignore
index e2762227..6e1ca3f1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -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;
+}