From facc28d900aeb9c084908c3a21094ee792f78a80 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Oct 2010 17:13:47 +0200 Subject: Add estimate_background --- src/Makefile.am | 5 +- src/Makefile.in | 26 ++++-- src/estimate_background.c | 230 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 253 insertions(+), 8 deletions(-) create mode 100644 src/estimate_background.c (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 2b5c11d0..fc9e3b07 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig compare_hkl \ powder_plot render_hkl calibrate_detector facetron cubeit \ - reintegrate + reintegrate estimate_background if HAVE_GTK bin_PROGRAMS += hdfsee @@ -70,4 +70,7 @@ reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \ thread-pool.c reintegrate_LDADD = @LIBS@ +estimate_background_SOURCES = estimate_background.c stream.c utils.c cell.c +estimate_background_LDADD = @LIBS@ + INCLUDES = "-I$(top_srcdir)/data" diff --git a/src/Makefile.in b/src/Makefile.in index 6c47670f..b4b8d63f 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -36,7 +36,8 @@ 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) \ - reintegrate$(EXEEXT) $(am__EXEEXT_1) + reintegrate$(EXEEXT) estimate_background$(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 @@ -70,6 +71,10 @@ am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \ symmetry.$(OBJEXT) stream.$(OBJEXT) thread-pool.$(OBJEXT) cubeit_OBJECTS = $(am_cubeit_OBJECTS) cubeit_DEPENDENCIES = +am_estimate_background_OBJECTS = estimate_background.$(OBJEXT) \ + stream.$(OBJEXT) utils.$(OBJEXT) cell.$(OBJEXT) +estimate_background_OBJECTS = $(am_estimate_background_OBJECTS) +estimate_background_DEPENDENCIES = am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \ peaks.$(OBJEXT) image.$(OBJEXT) geometry.$(OBJEXT) \ @@ -155,13 +160,14 @@ AM_V_GEN = $(am__v_GEN_$(V)) am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY)) am__v_GEN_0 = @echo " GEN " $@; 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) $(reintegrate_SOURCES) \ - $(render_hkl_SOURCES) + $(cubeit_SOURCES) $(estimate_background_SOURCES) \ + $(facetron_SOURCES) $(get_hkl_SOURCES) $(hdfsee_SOURCES) \ + $(indexamajig_SOURCES) $(pattern_sim_SOURCES) \ + $(powder_plot_SOURCES) $(process_hkl_SOURCES) \ + $(reintegrate_SOURCES) $(render_hkl_SOURCES) DIST_SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \ - $(cubeit_SOURCES) $(facetron_SOURCES) $(get_hkl_SOURCES) \ + $(cubeit_SOURCES) $(estimate_background_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) $(reintegrate_SOURCES) \ @@ -311,6 +317,8 @@ reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \ thread-pool.c reintegrate_LDADD = @LIBS@ +estimate_background_SOURCES = estimate_background.c stream.c utils.c cell.c +estimate_background_LDADD = @LIBS@ INCLUDES = "-I$(top_srcdir)/data" all: all-am @@ -392,6 +400,9 @@ compare_hkl$(EXEEXT): $(compare_hkl_OBJECTS) $(compare_hkl_DEPENDENCIES) cubeit$(EXEEXT): $(cubeit_OBJECTS) $(cubeit_DEPENDENCIES) @rm -f cubeit$(EXEEXT) $(AM_V_CCLD)$(LINK) $(cubeit_OBJECTS) $(cubeit_LDADD) $(LIBS) +estimate_background$(EXEEXT): $(estimate_background_OBJECTS) $(estimate_background_DEPENDENCIES) + @rm -f estimate_background$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(estimate_background_OBJECTS) $(estimate_background_LDADD) $(LIBS) facetron$(EXEEXT): $(facetron_OBJECTS) $(facetron_DEPENDENCIES) @rm -f facetron$(EXEEXT) $(AM_V_CCLD)$(LINK) $(facetron_OBJECTS) $(facetron_LDADD) $(LIBS) @@ -436,6 +447,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/diffraction.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/dirax.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/estimate_background.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/facetron.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/filters.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geometry.Po@am__quote@ diff --git a/src/estimate_background.c b/src/estimate_background.c new file mode 100644 index 00000000..d12bc37e --- /dev/null +++ b/src/estimate_background.c @@ -0,0 +1,230 @@ +/* + * estimate-background.c + * + * Like 'process_hkl', but for signal/noise ratios + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include "utils.h" +#include "statistics.h" +#include "sfac.h" +#include "reflections.h" +#include "symmetry.h" +#include "stream.h" + + +/* Number of bins for plot of resolution shells */ +#define NBINS (10) + + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"Estimate peak SNR from FEL Bragg intensities.\n" +"\n" +" -h, --help Display this help message.\n" +" -i, --input= Specify input filename (\"-\" for stdin).\n" +); +} + + + +int main(int argc, char *argv[]) +{ + int c; + char *filename = NULL; + FILE *fh; + int rval; + double total_vol, vol_per_shell; + double rmin, rmax; + double rmins[NBINS]; + double rmaxs[NBINS]; + double snrs[NBINS]; + unsigned int counts[NBINS]; + UnitCell *real_cell; + double overall_snr; + unsigned int overall_counts; + int i; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"input", 1, NULL, 'i'}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:f:a:r:", + longopts, NULL)) != -1) { + + switch (c) { + case 'h' : + show_help(argv[0]); + return 0; + + case 'i' : + filename = strdup(optarg); + break; + + case 0 : + break; + + default : + return 1; + } + + } + + if ( filename == NULL ) { + ERROR("Please specify filename using the -i option\n"); + return 1; + } + + /* Open the data stream */ + if ( strcmp(filename, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(filename, "r"); + } + free(filename); + if ( fh == NULL ) { + ERROR("Failed to open input file\n"); + return 1; + } + + /* FIXME: Fixed resolution shells */ + rmin = 0.120e9; + rmax = 1.172e9; + + for ( i=0; i %i %i %i %f\n", h, k, l, d/1e9); + + bin = -1; + for ( j=0; jrmins[j]) && (d<=rmaxs[j]) ) { + bin = j; + break; + } + } + if ( bin == -1 ) { + ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9); + continue; + } + + snr = max/bg; + snrs[bin] += snr; + counts[bin]++; + + } while ( !done ); + + } while ( !rval ); + + /* Print out results */ + overall_snr = 0.0; + overall_counts = 0; + for ( i=0; i0 ) { + overall_snr += snrs[i]; + overall_counts += counts[i]; + } + + } + + printf("Overall: %f (%f / %i)\n", overall_snr/(double)overall_counts, + overall_snr, overall_counts); + + + fclose(fh); + + return 0; +} -- cgit v1.2.3