diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | src/Makefile.am | 5 | ||||
| -rw-r--r-- | src/Makefile.in | 23 | ||||
| -rw-r--r-- | src/cubeit.c | 208 |
4 files changed, 230 insertions, 7 deletions
@@ -16,4 +16,5 @@ src/powder_plot src/render_hkl src/calibrate_detector src/facetron +src/cubeit *~ diff --git a/src/Makefile.am b/src/Makefile.am index 8f88ead7..bbf3bd21 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 calibrate_detector facetron + powder_plot render_hkl calibrate_detector facetron cubeit if HAVE_GTK bin_PROGRAMS += hdfsee @@ -58,4 +58,7 @@ calibrate_detector_LDADD = @LIBS@ facetron_SOURCES = facetron.c cell.c facetron_LDADD = @LIBS@ +cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c +cubeit_LDADD = @LIBS@ + INCLUDES = "-I$(top_srcdir)/data" diff --git a/src/Makefile.in b/src/Makefile.in index 85dc14ee..6a911fa7 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -35,7 +35,8 @@ POST_UNINSTALL = : 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) $(am__EXEEXT_1) + calibrate_detector$(EXEEXT) facetron$(EXEEXT) cubeit$(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 @@ -64,6 +65,10 @@ am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \ compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS) compare_hkl_DEPENDENCIES = am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) +am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \ + hdf5-file.$(OBJEXT) utils.$(OBJEXT) +cubeit_OBJECTS = $(am_cubeit_OBJECTS) +cubeit_DEPENDENCIES = facetron_OBJECTS = $(am_facetron_OBJECTS) facetron_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ @@ -137,12 +142,12 @@ 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) \ - $(facetron_SOURCES) $(get_hkl_SOURCES) $(hdfsee_SOURCES) \ - $(indexamajig_SOURCES) $(pattern_sim_SOURCES) \ - $(powder_plot_SOURCES) $(process_hkl_SOURCES) \ - $(render_hkl_SOURCES) + $(cubeit_SOURCES) $(facetron_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) \ - $(facetron_SOURCES) $(get_hkl_SOURCES) \ + $(cubeit_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) $(render_hkl_SOURCES) @@ -278,6 +283,8 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \ calibrate_detector_LDADD = @LIBS@ facetron_SOURCES = facetron.c cell.c facetron_LDADD = @LIBS@ +cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c +cubeit_LDADD = @LIBS@ INCLUDES = "-I$(top_srcdir)/data" all: all-am @@ -356,6 +363,9 @@ calibrate_detector$(EXEEXT): $(calibrate_detector_OBJECTS) $(calibrate_detector_ 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) +cubeit$(EXEEXT): $(cubeit_OBJECTS) $(cubeit_DEPENDENCIES) + @rm -f cubeit$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(cubeit_OBJECTS) $(cubeit_LDADD) $(LIBS) facetron$(EXEEXT): $(facetron_OBJECTS) $(facetron_DEPENDENCIES) @rm -f facetron$(EXEEXT) $(AM_V_CCLD)$(LINK) $(facetron_OBJECTS) $(facetron_LDADD) $(LIBS) @@ -391,6 +401,7 @@ distclean-compile: @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@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cubeit.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/detector.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/diffraction-gpu.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/diffraction.Po@am__quote@ diff --git a/src/cubeit.c b/src/cubeit.c new file mode 100644 index 00000000..ba7472aa --- /dev/null +++ b/src/cubeit.c @@ -0,0 +1,208 @@ +/* + * cubeit.c + * + * "Full integration" of diffraction data + * + * (c) 2006-2010 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 <errno.h> + +#include "image.h" +#include "cell.h" +#include "hdf5-file.h" + + +#define MAX_HITS (1024) + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"'Full integration' of diffraction data.\n" +"\n" +" -h, --help Display this help message.\n" +"\n" +" -i, --input=<filename> Specify the name of the input stream.\n" +" Can be '-' for stdin.\n" +); +} + + +static UnitCell *read_orientation_matrix(FILE *fh) +{ + float u, v, w; + struct rvec as, bs, cs; + UnitCell *cell; + char line[1024]; + + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) != 3 ) { + ERROR("Couldn't read a-star\n"); + return NULL; + } + as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) != 3 ) { + ERROR("Couldn't read b-star\n"); + return NULL; + } + bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; + if ( fgets(line, 1023, fh) == NULL ) return NULL; + if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) != 3 ) { + ERROR("Couldn't read c-star\n"); + return NULL; + } + cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; + cell = cell_new_from_axes(as, bs, cs); + + return cell; +} + + +static int find_chunk(FILE *fh, UnitCell **cell, char **filename) +{ + char line[1024]; + char *rval = NULL; + + do { + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + + chomp(line); + + if ( strncmp(line, "Reflections from indexing", 25) != 0 ) { + continue; + } + + *filename = strdup(line+29); + + /* Skip two lines (while checking for errors) */ + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + + *cell = read_orientation_matrix(fh); + if ( *cell == NULL ) { + STATUS("Got filename but no cell for %s\n", *filename); + continue; + } + + return 0; + + } while ( rval != NULL ); + + return 1; +} + + +int main(int argc, char *argv[]) +{ + int c; + char *infile = NULL; + FILE *fh; + UnitCell *cell; + char *filename; + unsigned int angles[180]; + 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:", longopts, NULL)) != -1) { + + switch (c) { + case 'h' : + show_help(argv[0]); + return 0; + + case 'i' : + infile = strdup(optarg); + break; + + case 0 : + break; + + default : + return 1; + } + + } + + if ( infile == NULL ) infile = strdup("-"); + if ( strcmp(infile, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(infile, "r"); + } + if ( fh == NULL ) { + ERROR("Couldn't open input stream '%s'\n", infile); + return ENOENT; + } + + /* Initialise histogram */ + for ( i=0; i<180; i++ ) angles[i] = 0; + + /* Loop over all successfully indexed patterns */ + while ( find_chunk(fh, &cell, &filename) == 0 ) { + + struct image image; + struct hdfile *hdfile; + double ang; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + unsigned int bin; + + //STATUS("Processing '%s'\n", filename); + + #if 0 + hdfile = hdfile_open(filename); + if ( hdfile == NULL ) { + return ENOENT; + } else if ( hdfile_set_image(hdfile, "/data/data0") ) { + ERROR("Couldn't select path\n"); + return ENOENT; + } + + hdf5_read(hdfile, &image, 0); + #endif + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + ang = angle_between(csx, csy, csz, 0.0, 0.0, 1.0); /* 0->pi */ + ang = rad2deg(ang); /* 0->180 deg */ + bin = rint(ang); + angles[bin]++; + + } + + for ( i=0; i<180; i++ ) { + STATUS("%i %i\n", i, angles[i]); + } + + return 0; +} |