From 74dfd122a3179a719e48ef89a39cbef9ce9dada7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 11 Dec 2010 15:13:28 -0800 Subject: "facetron" is now known as "partialator" --- .gitignore | 2 +- Makefile.am | 17 +- Makefile.in | 77 ++++----- src/facetron.c | 485 ----------------------------------------------------- src/partialator.c | 486 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 535 insertions(+), 532 deletions(-) delete mode 100644 src/facetron.c create mode 100644 src/partialator.c diff --git a/.gitignore b/.gitignore index 2abcd80b..e2762227 100644 --- a/.gitignore +++ b/.gitignore @@ -15,7 +15,7 @@ src/compare_hkl src/powder_plot src/render_hkl src/calibrate_detector -src/facetron +src/partialator src/cubeit src/reintegrate src/estimate_background diff --git a/Makefile.am b/Makefile.am index 0139cbbc..0806dd8c 100644 --- a/Makefile.am +++ b/Makefile.am @@ -4,7 +4,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/facetron src/reintegrate \ + src/calibrate_detector src/partialator src/reintegrate \ src/estimate_background src/check_hkl if HAVE_GTK @@ -69,16 +69,17 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ src/peaks.c src/detector.c src/cell.c \ src/thread-pool.c -src_facetron_SOURCES = src/facetron.c src/cell.c src/hdf5-file.c src/utils.c \ - src/detector.c src/peaks.c src/image.c src/geometry.c \ - src/reflections.c src/stream.c src/thread-pool.c \ - src/beam-parameters.c src/symmetry.c \ - src/post-refinement.c src/hrs-scaling.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 \ + src/geometry.c src/reflections.c src/stream.c \ + src/thread-pool.c src/beam-parameters.c \ + src/symmetry.c src/post-refinement.c \ + src/hrs-scaling.c if HAVE_CAIRO src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ - src/detector.c src/render.c src/filters.c src/image.c \ - src/symmetry.c src/stream.c src/thread-pool.c + src/detector.c src/render.c src/filters.c src/image.c \ + src/symmetry.c src/stream.c src/thread-pool.c endif src_reintegrate_SOURCES = src/reintegrate.c src/cell.c src/hdf5-file.c \ diff --git a/Makefile.in b/Makefile.in index 56ab3939..b96acc49 100644 --- a/Makefile.in +++ b/Makefile.in @@ -39,7 +39,7 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \ src/get_hkl$(EXEEXT) src/indexamajig$(EXEEXT) \ src/compare_hkl$(EXEEXT) src/powder_plot$(EXEEXT) \ src/render_hkl$(EXEEXT) src/calibrate_detector$(EXEEXT) \ - src/facetron$(EXEEXT) src/reintegrate$(EXEEXT) \ + src/partialator$(EXEEXT) src/reintegrate$(EXEEXT) \ src/estimate_background$(EXEEXT) src/check_hkl$(EXEEXT) \ $(am__EXEEXT_1) $(am__EXEEXT_2) @HAVE_GTK_TRUE@am__append_1 = src/hdfsee @@ -120,16 +120,6 @@ src_estimate_background_OBJECTS = \ $(am_src_estimate_background_OBJECTS) src_estimate_background_LDADD = $(LDADD) src_estimate_background_DEPENDENCIES = $(top_builddir)/lib/libgnu.a -am_src_facetron_OBJECTS = src/facetron.$(OBJEXT) src/cell.$(OBJEXT) \ - src/hdf5-file.$(OBJEXT) src/utils.$(OBJEXT) \ - src/detector.$(OBJEXT) src/peaks.$(OBJEXT) src/image.$(OBJEXT) \ - src/geometry.$(OBJEXT) src/reflections.$(OBJEXT) \ - src/stream.$(OBJEXT) src/thread-pool.$(OBJEXT) \ - src/beam-parameters.$(OBJEXT) src/symmetry.$(OBJEXT) \ - src/post-refinement.$(OBJEXT) src/hrs-scaling.$(OBJEXT) -src_facetron_OBJECTS = $(am_src_facetron_OBJECTS) -src_facetron_LDADD = $(LDADD) -src_facetron_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_get_hkl_OBJECTS = src/get_hkl.$(OBJEXT) src/sfac.$(OBJEXT) \ src/cell.$(OBJEXT) src/utils.$(OBJEXT) \ src/reflections.$(OBJEXT) src/symmetry.$(OBJEXT) \ @@ -168,6 +158,16 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS) src_indexamajig_LDADD = $(LDADD) src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a +am_src_partialator_OBJECTS = src/partialator.$(OBJEXT) \ + src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) src/utils.$(OBJEXT) \ + src/detector.$(OBJEXT) src/peaks.$(OBJEXT) src/image.$(OBJEXT) \ + src/geometry.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/stream.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/beam-parameters.$(OBJEXT) src/symmetry.$(OBJEXT) \ + src/post-refinement.$(OBJEXT) src/hrs-scaling.$(OBJEXT) +src_partialator_OBJECTS = $(am_src_partialator_OBJECTS) +src_partialator_LDADD = $(LDADD) +src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am__src_pattern_sim_SOURCES_DIST = src/pattern_sim.c src/diffraction.c \ src/utils.c src/image.c src/cell.c src/hdf5-file.c \ src/detector.c src/sfac.c src/peaks.c src/reflections.c \ @@ -237,17 +237,17 @@ am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY)) am__v_GEN_0 = @echo " GEN " $@; SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \ $(src_compare_hkl_SOURCES) $(src_cubeit_SOURCES) \ - $(src_estimate_background_SOURCES) $(src_facetron_SOURCES) \ - $(src_get_hkl_SOURCES) $(src_hdfsee_SOURCES) \ - $(src_indexamajig_SOURCES) $(src_pattern_sim_SOURCES) \ + $(src_estimate_background_SOURCES) $(src_get_hkl_SOURCES) \ + $(src_hdfsee_SOURCES) $(src_indexamajig_SOURCES) \ + $(src_partialator_SOURCES) $(src_pattern_sim_SOURCES) \ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \ $(src_reintegrate_SOURCES) $(src_render_hkl_SOURCES) DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \ $(am__src_cubeit_SOURCES_DIST) \ - $(src_estimate_background_SOURCES) $(src_facetron_SOURCES) \ - $(src_get_hkl_SOURCES) $(am__src_hdfsee_SOURCES_DIST) \ - $(am__src_indexamajig_SOURCES_DIST) \ + $(src_estimate_background_SOURCES) $(src_get_hkl_SOURCES) \ + $(am__src_hdfsee_SOURCES_DIST) \ + $(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) @@ -585,15 +585,16 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ src/peaks.c src/detector.c src/cell.c \ src/thread-pool.c -src_facetron_SOURCES = src/facetron.c src/cell.c src/hdf5-file.c src/utils.c \ - src/detector.c src/peaks.c src/image.c src/geometry.c \ - src/reflections.c src/stream.c src/thread-pool.c \ - src/beam-parameters.c src/symmetry.c \ - src/post-refinement.c src/hrs-scaling.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 \ + src/geometry.c src/reflections.c src/stream.c \ + src/thread-pool.c src/beam-parameters.c \ + src/symmetry.c src/post-refinement.c \ + src/hrs-scaling.c @HAVE_CAIRO_TRUE@src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ -@HAVE_CAIRO_TRUE@ src/detector.c src/render.c src/filters.c src/image.c \ -@HAVE_CAIRO_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c +@HAVE_CAIRO_TRUE@ src/detector.c src/render.c src/filters.c src/image.c \ +@HAVE_CAIRO_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c src_reintegrate_SOURCES = src/reintegrate.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ @@ -783,21 +784,10 @@ src/estimate_background.$(OBJEXT): src/$(am__dirstamp) \ src/estimate_background$(EXEEXT): $(src_estimate_background_OBJECTS) $(src_estimate_background_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/estimate_background$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_estimate_background_OBJECTS) $(src_estimate_background_LDADD) $(LIBS) -src/facetron.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/geometry.$(OBJEXT): src/$(am__dirstamp) \ +src/get_hkl.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/beam-parameters.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) -src/post-refinement.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/hrs-scaling.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) -src/facetron$(EXEEXT): $(src_facetron_OBJECTS) $(src_facetron_DEPENDENCIES) src/$(am__dirstamp) - @rm -f src/facetron$(EXEEXT) - $(AM_V_CCLD)$(LINK) $(src_facetron_OBJECTS) $(src_facetron_LDADD) $(LIBS) -src/get_hkl.$(OBJEXT): src/$(am__dirstamp) \ - src/$(DEPDIR)/$(am__dirstamp) src/get_hkl$(EXEEXT): $(src_get_hkl_OBJECTS) $(src_get_hkl_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/get_hkl$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_get_hkl_OBJECTS) $(src_get_hkl_LDADD) $(LIBS) @@ -818,6 +808,8 @@ src/mosflm.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/templates.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/geometry.$(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) \ @@ -825,6 +817,15 @@ src/cl-utils.$(OBJEXT): src/$(am__dirstamp) \ src/indexamajig$(EXEEXT): $(src_indexamajig_OBJECTS) $(src_indexamajig_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/indexamajig$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_indexamajig_OBJECTS) $(src_indexamajig_LDADD) $(LIBS) +src/partialator.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/post-refinement.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/hrs-scaling.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) +src/partialator$(EXEEXT): $(src_partialator_OBJECTS) $(src_partialator_DEPENDENCIES) src/$(am__dirstamp) + @rm -f src/partialator$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(src_partialator_OBJECTS) $(src_partialator_LDADD) $(LIBS) src/pattern_sim.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/pattern_sim$(EXEEXT): $(src_pattern_sim_OBJECTS) $(src_pattern_sim_DEPENDENCIES) src/$(am__dirstamp) @@ -868,7 +869,6 @@ mostlyclean-compile: -rm -f src/dirax.$(OBJEXT) -rm -f src/displaywindow.$(OBJEXT) -rm -f src/estimate_background.$(OBJEXT) - -rm -f src/facetron.$(OBJEXT) -rm -f src/filters.$(OBJEXT) -rm -f src/geometry.$(OBJEXT) -rm -f src/get_hkl.$(OBJEXT) @@ -879,6 +879,7 @@ mostlyclean-compile: -rm -f src/index.$(OBJEXT) -rm -f src/indexamajig.$(OBJEXT) -rm -f src/mosflm.$(OBJEXT) + -rm -f src/partialator.$(OBJEXT) -rm -f src/pattern_sim.$(OBJEXT) -rm -f src/peaks.$(OBJEXT) -rm -f src/post-refinement.$(OBJEXT) @@ -913,7 +914,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/dirax.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/displaywindow.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/estimate_background.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/facetron.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/filters.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/geometry.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/get_hkl.Po@am__quote@ @@ -924,6 +924,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/index.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/indexamajig.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/mosflm.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/partialator.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/pattern_sim.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/peaks.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/post-refinement.Po@am__quote@ diff --git a/src/facetron.c b/src/facetron.c deleted file mode 100644 index 5223e070..00000000 --- a/src/facetron.c +++ /dev/null @@ -1,485 +0,0 @@ -/* - * facetron.c - * - * Profile fitting for coherent nanocrystallography - * - * (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 -#include - -#include "utils.h" -#include "hdf5-file.h" -#include "symmetry.h" -#include "reflections.h" -#include "stream.h" -#include "geometry.h" -#include "peaks.h" -#include "thread-pool.h" -#include "beam-parameters.h" -#include "post-refinement.h" -#include "hrs-scaling.h" - - -/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (100) - - -static void show_help(const char *s) -{ - printf("Syntax: %s [options]\n\n", s); - printf( -"Post-refinement and profile fitting for coherent nanocrystallography.\n" -"\n" -" -h, --help Display this help message.\n" -"\n" -" -i, --input= Specify the name of the input 'stream'.\n" -" (must be a file, not e.g. stdin)\n" -" -o, --output= Output filename. Default: facetron.hkl.\n" -" -g. --geometry= Get detector geometry from file.\n" -" -b, --beam= Get beam parameters from file (provides initial\n" -" values for parameters, and nominal wavelengths\n" -" if no per-shot value is found in an HDF5 file.\n" -" -x, --prefix=

Prefix filenames from input file with

.\n" -" --basename Remove the directory parts of the filenames.\n" -" --no-check-prefix Don't attempt to correct the --prefix.\n" -" -y, --symmetry= Merge according to symmetry .\n" -" -n, --iterations= Run cycles of post-refinement.\n" -"\n" -" -j Run analyses in parallel.\n"); -} - - -struct refine_args -{ - const char *sym; - ReflItemList *obs; - double *i_full; - struct image *image; - FILE *graph; - FILE *pgraph; -}; - - -static void refine_image(int mytask, void *tasks) -{ - struct refine_args *all_args = tasks; - struct refine_args *pargs = &all_args[mytask]; - struct image *image = pargs->image; - double nominal_photon_energy = pargs->image->beam->photon_energy; - struct hdfile *hdfile; - struct cpeak *spots; - int n, i; - double dev, last_dev; - - hdfile = hdfile_open(image->filename); - if ( hdfile == NULL ) { - ERROR("Couldn't open '%s'\n", image->filename); - return; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - hdfile_close(hdfile); - return; - } - - if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) { - ERROR("Couldn't read '%s'\n", image->filename); - hdfile_close(hdfile); - return; - } - - double a, b, c, al, be, ga; - cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga); - STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n", - a/1.0e-9, b/1.0e-9, c/1.0e-9, - rad2deg(al), rad2deg(be), rad2deg(ga)); - - spots = find_intersections(image, image->indexed_cell, &n, 0); - dev = +INFINITY; - i = 0; - do { - last_dev = dev; - dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n); - STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev); - i++; - } while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) ); - mean_partial_dev(image, spots, n, pargs->sym, - pargs->i_full, pargs->graph); - if ( pargs->pgraph ) { - fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev); - } - - free(image->data); - if ( image->flags != NULL ) free(image->flags); - hdfile_close(hdfile); - free(spots); - - /* Muppet proofing */ - image->data = NULL; - image->flags = NULL; -} - - -static void refine_all(struct image *images, int n_total_patterns, - struct detector *det, const char *sym, - ReflItemList *obs, double *i_full, int nthreads, - FILE *graph, FILE *pgraph) -{ - struct refine_args *tasks; - int i; - - tasks = malloc(n_total_patterns * sizeof(struct refine_args)); - for ( i=0; ibeam->photon_energy; - - hdfile = hdfile_open(image->filename); - if ( hdfile == NULL ) { - ERROR("Couldn't open '%s'\n", image->filename); - return; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - hdfile_close(hdfile); - return; - } - - if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) { - ERROR("Couldn't read '%s'\n", image->filename); - hdfile_close(hdfile); - return; - } - - /* Figure out which spots should appear in this pattern */ - spots = find_intersections(image, image->indexed_cell, &n, 0); - - /* For each reflection, estimate the partiality */ - for ( j=0; jdata); - if ( image->flags != NULL ) free(image->flags); - hdfile_close(hdfile); - image->cpeaks = spots; - - /* Muppet proofing */ - image->data = NULL; - image->flags = NULL; -} - - -int main(int argc, char *argv[]) -{ - int c; - char *infile = NULL; - char *outfile = NULL; - char *geomfile = NULL; - char *prefix = NULL; - char *sym = NULL; - FILE *fh; - int nthreads = 1; - int config_basename = 0; - int config_checkprefix = 1; - struct detector *det; - double *i_full; - unsigned int *cts; - ReflItemList *obs; - int i; - int n_total_patterns; - struct image *images; - int n_iter = 10; - struct beam_params *beam = NULL; - - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"input", 1, NULL, 'i'}, - {"output", 1, NULL, 'o'}, - {"geometry", 1, NULL, 'g'}, - {"beam", 1, NULL, 'b'}, - {"prefix", 1, NULL, 'x'}, - {"basename", 0, &config_basename, 1}, - {"no-check-prefix", 0, &config_checkprefix, 0}, - {"symmetry", 1, NULL, 'y'}, - {"iterations", 1, NULL, 'n'}, - {0, 0, NULL, 0} - }; - - /* Short options */ - while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:", - longopts, NULL)) != -1) - { - - switch (c) { - case 'h' : - show_help(argv[0]); - return 0; - - case 'i' : - infile = strdup(optarg); - break; - - case 'g' : - geomfile = strdup(optarg); - break; - - case 'x' : - prefix = strdup(optarg); - break; - - case 'j' : - nthreads = atoi(optarg); - break; - - case 'y' : - sym = strdup(optarg); - break; - - case 'o' : - outfile = strdup(optarg); - break; - - case 'n' : - n_iter = atoi(optarg); - break; - - case 'b' : - beam = get_beam_parameters(optarg); - if ( beam == NULL ) { - ERROR("Failed to load beam parameters" - " from '%s'\n", optarg); - return 1; - } - break; - - case 0 : - break; - - default : - return 1; - } - - } - - /* Sanitise input filename and open */ - if ( infile == NULL ) { - infile = strdup("-"); - } - if ( strcmp(infile, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(infile, "r"); - } - if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", infile); - return 1; - } - free(infile); - - /* Sanitise output filename */ - if ( outfile == NULL ) { - outfile = strdup("facetron.hkl"); - } - - /* Sanitise prefix */ - if ( prefix == NULL ) { - prefix = strdup(""); - } else { - if ( config_checkprefix ) { - prefix = check_prefix(prefix); - } - } - - if ( sym == NULL ) sym = strdup("1"); - - /* Get detector geometry */ - det = get_detector_geometry(geomfile); - if ( det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geomfile); - return 1; - } - free(geomfile); - - if ( beam == NULL ) { - ERROR("You must provide a beam parameters file.\n"); - return 1; - } - - /* Prepare for iteration */ - i_full = new_list_intensity(); - obs = new_items(); - - n_total_patterns = count_patterns(fh); - STATUS("There are %i patterns to process\n", n_total_patterns); - - images = malloc(n_total_patterns * sizeof(struct image)); - if ( images == NULL ) { - ERROR("Couldn't allocate memory for images.\n"); - return 1; - } - - /* Fill in what we know about the images so far */ - rewind(fh); - for ( i=0; idivergence; - images[i].bw = beam->bandwidth; - images[i].det = det; - images[i].beam = beam; - images[i].osf = 1.0; - images[i].profile_radius = 0.005e9; - - /* Muppet proofing */ - images[i].data = NULL; - images[i].flags = NULL; - - /* Get reflections from this image */ - integrate_image(&images[i]); - - progress_bar(i, n_total_patterns-1, "Loading pattern data"); - - } - fclose(fh); - free(prefix); - - cts = new_list_count(); - - /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); - scale_intensities(images, n_total_patterns, sym); - - /* Iterate */ - for ( i=0; ipanels); - free(det); - free(beam); - free(cts); - for ( i=0; i + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "utils.h" +#include "hdf5-file.h" +#include "symmetry.h" +#include "reflections.h" +#include "stream.h" +#include "geometry.h" +#include "peaks.h" +#include "thread-pool.h" +#include "beam-parameters.h" +#include "post-refinement.h" +#include "hrs-scaling.h" + + +/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ +#define MAX_CYCLES (100) + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"Scaling and post refinement for coherent nanocrystallography.\n" +"\n" +" -h, --help Display this help message.\n" +"\n" +" -i, --input= Specify the name of the input 'stream'.\n" +" (must be a file, not e.g. stdin)\n" +" -o, --output= Output filename. Default: facetron.hkl.\n" +" -g. --geometry= Get detector geometry from file.\n" +" -b, --beam= Get beam parameters from file, which provides\n" +" initial values for parameters, and nominal\n" +" wavelengths if no per-shot value is found in \n" +" an HDF5 file.\n" +" -x, --prefix=

Prefix filenames from input file with

.\n" +" --basename Remove the directory parts of the filenames.\n" +" --no-check-prefix Don't attempt to correct the --prefix.\n" +" -y, --symmetry= Merge according to symmetry .\n" +" -n, --iterations= Run cycles of scaling and post-refinement.\n" +"\n" +" -j Run analyses in parallel.\n"); +} + + +struct refine_args +{ + const char *sym; + ReflItemList *obs; + double *i_full; + struct image *image; + FILE *graph; + FILE *pgraph; +}; + + +static void refine_image(int mytask, void *tasks) +{ + struct refine_args *all_args = tasks; + struct refine_args *pargs = &all_args[mytask]; + struct image *image = pargs->image; + double nominal_photon_energy = pargs->image->beam->photon_energy; + struct hdfile *hdfile; + struct cpeak *spots; + int n, i; + double dev, last_dev; + + hdfile = hdfile_open(image->filename); + if ( hdfile == NULL ) { + ERROR("Couldn't open '%s'\n", image->filename); + return; + } else if ( hdfile_set_image(hdfile, "/data/data0") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return; + } + + if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) { + ERROR("Couldn't read '%s'\n", image->filename); + hdfile_close(hdfile); + return; + } + + double a, b, c, al, be, ga; + cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga); + STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n", + a/1.0e-9, b/1.0e-9, c/1.0e-9, + rad2deg(al), rad2deg(be), rad2deg(ga)); + + spots = find_intersections(image, image->indexed_cell, &n, 0); + dev = +INFINITY; + i = 0; + do { + last_dev = dev; + dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n); + STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev); + i++; + } while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) ); + mean_partial_dev(image, spots, n, pargs->sym, + pargs->i_full, pargs->graph); + if ( pargs->pgraph ) { + fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev); + } + + free(image->data); + if ( image->flags != NULL ) free(image->flags); + hdfile_close(hdfile); + free(spots); + + /* Muppet proofing */ + image->data = NULL; + image->flags = NULL; +} + + +static void refine_all(struct image *images, int n_total_patterns, + struct detector *det, const char *sym, + ReflItemList *obs, double *i_full, int nthreads, + FILE *graph, FILE *pgraph) +{ + struct refine_args *tasks; + int i; + + tasks = malloc(n_total_patterns * sizeof(struct refine_args)); + for ( i=0; ibeam->photon_energy; + + hdfile = hdfile_open(image->filename); + if ( hdfile == NULL ) { + ERROR("Couldn't open '%s'\n", image->filename); + return; + } else if ( hdfile_set_image(hdfile, "/data/data0") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return; + } + + if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) { + ERROR("Couldn't read '%s'\n", image->filename); + hdfile_close(hdfile); + return; + } + + /* Figure out which spots should appear in this pattern */ + spots = find_intersections(image, image->indexed_cell, &n, 0); + + /* For each reflection, estimate the partiality */ + for ( j=0; jdata); + if ( image->flags != NULL ) free(image->flags); + hdfile_close(hdfile); + image->cpeaks = spots; + + /* Muppet proofing */ + image->data = NULL; + image->flags = NULL; +} + + +int main(int argc, char *argv[]) +{ + int c; + char *infile = NULL; + char *outfile = NULL; + char *geomfile = NULL; + char *prefix = NULL; + char *sym = NULL; + FILE *fh; + int nthreads = 1; + int config_basename = 0; + int config_checkprefix = 1; + struct detector *det; + double *i_full; + unsigned int *cts; + ReflItemList *obs; + int i; + int n_total_patterns; + struct image *images; + int n_iter = 10; + struct beam_params *beam = NULL; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, + {"geometry", 1, NULL, 'g'}, + {"beam", 1, NULL, 'b'}, + {"prefix", 1, NULL, 'x'}, + {"basename", 0, &config_basename, 1}, + {"no-check-prefix", 0, &config_checkprefix, 0}, + {"symmetry", 1, NULL, 'y'}, + {"iterations", 1, NULL, 'n'}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:", + longopts, NULL)) != -1) + { + + switch (c) { + case 'h' : + show_help(argv[0]); + return 0; + + case 'i' : + infile = strdup(optarg); + break; + + case 'g' : + geomfile = strdup(optarg); + break; + + case 'x' : + prefix = strdup(optarg); + break; + + case 'j' : + nthreads = atoi(optarg); + break; + + case 'y' : + sym = strdup(optarg); + break; + + case 'o' : + outfile = strdup(optarg); + break; + + case 'n' : + n_iter = atoi(optarg); + break; + + case 'b' : + beam = get_beam_parameters(optarg); + if ( beam == NULL ) { + ERROR("Failed to load beam parameters" + " from '%s'\n", optarg); + return 1; + } + break; + + case 0 : + break; + + default : + return 1; + } + + } + + /* Sanitise input filename and open */ + if ( infile == NULL ) { + infile = strdup("-"); + } + if ( strcmp(infile, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(infile, "r"); + } + if ( fh == NULL ) { + ERROR("Failed to open input file '%s'\n", infile); + return 1; + } + free(infile); + + /* Sanitise output filename */ + if ( outfile == NULL ) { + outfile = strdup("facetron.hkl"); + } + + /* Sanitise prefix */ + if ( prefix == NULL ) { + prefix = strdup(""); + } else { + if ( config_checkprefix ) { + prefix = check_prefix(prefix); + } + } + + if ( sym == NULL ) sym = strdup("1"); + + /* Get detector geometry */ + det = get_detector_geometry(geomfile); + if ( det == NULL ) { + ERROR("Failed to read detector geometry from '%s'\n", geomfile); + return 1; + } + free(geomfile); + + if ( beam == NULL ) { + ERROR("You must provide a beam parameters file.\n"); + return 1; + } + + /* Prepare for iteration */ + i_full = new_list_intensity(); + obs = new_items(); + + n_total_patterns = count_patterns(fh); + STATUS("There are %i patterns to process\n", n_total_patterns); + + images = malloc(n_total_patterns * sizeof(struct image)); + if ( images == NULL ) { + ERROR("Couldn't allocate memory for images.\n"); + return 1; + } + + /* Fill in what we know about the images so far */ + rewind(fh); + for ( i=0; idivergence; + images[i].bw = beam->bandwidth; + images[i].det = det; + images[i].beam = beam; + images[i].osf = 1.0; + images[i].profile_radius = 0.005e9; + + /* Muppet proofing */ + images[i].data = NULL; + images[i].flags = NULL; + + /* Get reflections from this image */ + integrate_image(&images[i]); + + progress_bar(i, n_total_patterns-1, "Loading pattern data"); + + } + fclose(fh); + free(prefix); + + cts = new_list_count(); + + /* Make initial estimates */ + STATUS("Performing initial scaling.\n"); + scale_intensities(images, n_total_patterns, sym); + + /* Iterate */ + for ( i=0; ipanels); + free(det); + free(beam); + free(cts); + for ( i=0; i