diff options
Diffstat (limited to 'libcrystfel')
33 files changed, 2928 insertions, 602 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 7bb08685..fef2e11f 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -11,7 +11,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ src/xds.c src/integration.c src/predict-refine.c \ src/histogram.c src/events.c src/felix.c \ - src/taketwo.c + src/taketwo.c src/peakfinder8.c if HAVE_FFTW libcrystfel_la_SOURCES += src/asdf.c @@ -32,7 +32,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/xds.h src/predict-refine.h \ src/integration.h src/histogram.h \ src/events.h src/asdf.h src/felix.h \ - src/taketwo.h + src/taketwo.h src/peakfinder8.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c index 49917ad9..0c43fe7a 100644 --- a/libcrystfel/src/asdf.c +++ b/libcrystfel/src/asdf.c @@ -4,12 +4,12 @@ * Alexandra's Superior Direction Finder, or * Algorithm Similar to DirAx, FFT-based * - * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de> - * 2015 Thomas White <taw@physics.org> + * 2015,2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -28,8 +28,6 @@ * */ -#define _ISOC99_SOURCE - #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -1103,7 +1101,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, } -int run_asdf(struct image *image, IndexingPrivate *ipriv) +int run_asdf(struct image *image, void *ipriv) { int i, j; @@ -1204,8 +1202,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *asdf_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct asdf_private *dp; int need_cell = 0; @@ -1234,11 +1232,11 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell, dp->fftw = fftw_vars_new(); - return (IndexingPrivate *)dp; + return (void *)dp; } -void asdf_cleanup(IndexingPrivate *pp) +void asdf_cleanup(void *pp) { struct asdf_private *p; p = (struct asdf_private *)pp; diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h index 1e492a6f..f130d63d 100644 --- a/libcrystfel/src/asdf.h +++ b/libcrystfel/src/asdf.h @@ -4,12 +4,12 @@ * Alexandra's Superior Direction Finder, or * Algorithm Similar to DirAx, FFT-based * - * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de> - * 2015 Thomas White <taw@physics.org> + * 2015,2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -43,33 +43,33 @@ extern "C" { #ifdef HAVE_FFTW -extern int run_asdf(struct image *image, IndexingPrivate *ipriv); +extern int run_asdf(struct image *image, void *ipriv); -extern IndexingPrivate *asdf_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl); +extern void *asdf_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, + float *ltl); -extern void asdf_cleanup(IndexingPrivate *pp); +extern void asdf_cleanup(void *pp); #else /* HAVE_FFTW */ -int run_asdf(struct image *image, IndexingPrivate *ipriv) +int run_asdf(struct image *image, void *ipriv) { ERROR("This copy of CrystFEL was compiled without FFTW support.\n"); return 0; } -IndexingPrivate *asdf_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl) +void *asdf_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, + float *ltl) { ERROR("This copy of CrystFEL was compiled without FFTW support.\n"); ERROR("To use asdf indexing, recompile with FFTW.\n"); return NULL; } -void asdf_cleanup(IndexingPrivate *pp) +void asdf_cleanup(void *pp) { } diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 8fa3b894..6610abc3 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -272,6 +272,12 @@ void cell_print(UnitCell *cell) STATUS("c* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n", csx, csy, csz, modulus(csx, csy, csz)); + STATUS("alpha* = %6.2f deg, beta* = %6.2f deg, " + "gamma* = %6.2f deg\n", + rad2deg(angle_between(bsx, bsy, bsz, csx, csy, csz)), + rad2deg(angle_between(asx, asy, asz, csx, csy, csz)), + rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz))); + STATUS("Cell representation is %s.\n", cell_rep(cell)); } else { diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index ec591e24..3de61073 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -460,10 +460,12 @@ int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, case CELL_REP_RECIP: /* Convert reciprocal -> crystallographic. * Start by converting reciprocal -> cartesian */ - cell_invert(cell->axs, cell->ays, cell->azs, - cell->bxs, cell->bys, cell->bzs, - cell->cxs, cell->cys, cell->czs, - &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + if ( cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) return 1; /* Now convert cartesian -> crystallographic */ *a = modulus(ax, ay, az); diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 42ffbe40..e3b351ff 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -31,8 +31,6 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE #include <stdlib.h> #include <math.h> #include <stdio.h> @@ -530,37 +528,15 @@ int panel_number(struct detector *det, struct panel *p) } -void fill_in_values(struct detector *det, struct hdfile *f, struct event* ev) +void adjust_centering_for_rail(struct panel *p) { - int i; - - for ( i=0; i<det->n_panels; i++ ) { - - double offs; - struct panel *p = &det->panels[i]; - - if ( p->clen_from != NULL ) { - - double val; - int r; + double offs; - r = hdfile_get_value(f, p->clen_from, ev, &val, - H5T_NATIVE_DOUBLE); - if ( r ) { - ERROR("Failed to read '%s'\n", p->clen_from); - } else { - p->clen = val * 1.0e-3; - } - - } - - /* Offset in +z direction from calibrated clen to actual */ - offs = p->clen - p->clen_for_centering; - p->cnx += p->rail_x * offs; - p->cny += p->rail_y * offs; - p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs; - - } + /* Offset in +z direction from calibrated clen to actual */ + offs = p->clen - p->clen_for_centering; + p->cnx += p->rail_x * offs; + p->cny += p->rail_y * offs; + p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs; } @@ -911,7 +887,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key, panel->adu_per_eV = atof(val); } else if ( strcmp(key, "adu_per_photon") == 0 ) { panel->adu_per_photon = atof(val); - STATUS("got adu per photon: %s\n", val); } else if ( strcmp(key, "rigid_group") == 0 ) { add_to_rigid_group(find_or_add_rg(det, val), panel); } else if ( strcmp(key, "clen") == 0 ) { @@ -1079,12 +1054,15 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam, } else if ( strcmp(key, "photon_energy") == 0 ) { if ( beam != NULL ) { - if ( strncmp(val, "/", 1) == 0 ) { + double v; + char *end; + v = strtod(val, &end); + if ( (val[0] != '\0') && (end[0] == '\0') ) { + beam->photon_energy = v; + beam->photon_energy_from = NULL; + } else { beam->photon_energy = 0.0; beam->photon_energy_from = strdup(val); - } else { - beam->photon_energy = atof(val); - beam->photon_energy_from = NULL; } } @@ -1204,7 +1182,7 @@ struct detector *get_detector_geometry_2(const char *filename, int i; int rgi, rgci; int reject = 0; - int path_dim; + int path_dim, mask_path_dim; int dim_dim; int dim_reject = 0; int dim_dim_reject = 0; @@ -1389,6 +1367,7 @@ struct detector *get_detector_geometry_2(const char *filename, } + mask_path_dim = -1; for ( i=0; i<det->n_panels; i++ ) { int panel_mask_dim = 0; @@ -1398,8 +1377,7 @@ struct detector *get_detector_geometry_2(const char *filename, next_instance = det->panels[i].mask; - while(next_instance) - { + while ( next_instance ) { next_instance = strstr(next_instance, "%"); if ( next_instance != NULL ) { next_instance += 1*sizeof(char); @@ -1407,18 +1385,29 @@ struct detector *get_detector_geometry_2(const char *filename, } } - if ( panel_mask_dim != path_dim ) { - dim_reject = 1; + if ( mask_path_dim == -1 ) { + mask_path_dim = panel_mask_dim; + } else { + if ( panel_mask_dim != mask_path_dim ) { + dim_reject = 1; + } } + } } - if ( dim_reject == 1) { + if ( dim_reject == 1 ) { ERROR("All panels' data and mask entries must have the same " "number of placeholders\n"); reject = 1; } + if ( mask_path_dim > path_dim ) { + ERROR("Number of placeholders in mask cannot be larger than " + "for data\n"); + reject = 1; + } + det->path_dim = path_dim; dim_dim_reject = 0; @@ -1993,6 +1982,11 @@ double largest_q(struct image *image) struct rvec q; double tt; + if ( image->det == NULL ) { + ERROR("No detector geometry. assuming detector is infinite!\n"); + return INFINITY; + } + q = get_q_for_panel(image->det->furthest_out_panel, image->det->furthest_out_fs, image->det->furthest_out_ss, diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index 04e3c1ec..a2be2b47 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -3,12 +3,12 @@ * * Detector properties * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2009-2016 Thomas White <taw@physics.org> + * 2009-2017 Thomas White <taw@physics.org> * 2011-2012 Richard Kirian <rkirian@asu.edu> * 2014 Valerio Mariani * 2011 Andrew Aquila @@ -42,7 +42,6 @@ struct rg_collection; struct detector; struct panel; struct badregion; -struct detector; struct beam_params; struct hdfile; struct event; @@ -80,6 +79,47 @@ struct rg_collection }; +/** + * panel: + * @name: Text name for the panel (fixed length array) + * @cnx: Location of corner, in pixels, x coordinate + * @cny: Location of corner, in pixels, y coordinate + * @coffset: The offset to be applied from @clen (which may come from elsewhere) + * @clen: The distance from the interaction point to the corner of the first pixel + * @clen_from: Location to get @clen from, e.g. from HDF5 file + * @mask: Location of mask data + * @mask_file: Filename for mask data + * @satmap: Location of per-pixel saturation map + * @satmap_file: Filename for saturation map + * @res: Resolution of panel in pixels per metre + * @badrow: Readout direction (for filtering out clusters of peaks) + * @no_index: Non-zero if panel is entirely "bad" + * @adu_per_photon: Number of detector intensity units per photon + * @adu_per_eV: Number of detector intensity units per eV of photon energy + * @max_adu: Saturation value + * @dim_structure: Dimension structure + * @fsx: Real-space x-direction of data fast-scan direction + * @fsy: Real-space y-direction of data fast-scan direction + * @fsz: Real-space z-direction of data fast-scan direction + * @ssx: Real-space x-direction of data slow-scan direction + * @ssy: Real-space y-direction of data slow-scan direction + * @ssz: Real-space z-direction of data slow-scan direction + * @rail_x: x direction of camera length "rail" + * @rail_y: y direction of camera length "rail" + * @rail_z: z direction of camera length "rail" + * @clen_for_centering: Value of clen (without coffset) at which beam is centered + * @xfs: Data fast-scan direction of real-space x-direction + * @yfs: Data fast-scan direction of real-space y-direction + * @xss: Data slow-scan direction of real-space x-direction + * @yss: Data slow-scan direction of real-space y-direction + * @orig_min_fs: Minimum fs coordinate of data in file + * @orig_max_fs: Maximum fs coordinate of data in file + * @orig_min_ss: Minimum ss coordinate of data in file (inclusive) + * @orig_max_ss: Maximum ss coordinate of data in file (inclusive) + * @data: Location of data in file + * @w: Width of panel + * @h: Height of panel + */ struct panel { char name[1024]; /* Name for this panel */ @@ -224,9 +264,8 @@ extern void get_pixel_extents(struct detector *det, double *min_x, double *min_y, double *max_x, double *max_y); -extern void fill_in_values(struct detector *det, struct hdfile *f, - struct event* ev); extern void fill_in_adu(struct image *image); +extern void adjust_centering_for_rail(struct panel *p); extern int panel_is_in_rigid_group(const struct rigid_group *rg, struct panel *p); diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 19f35696..e9466a24 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2014,2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -532,7 +532,7 @@ static void write_drx(struct image *image) } -int run_dirax(struct image *image, IndexingPrivate *ipriv) +int run_dirax(struct image *image, void *ipriv) { unsigned int opts; int status; @@ -638,8 +638,8 @@ int run_dirax(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *dirax_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct dirax_private *dp; int need_cell = 0; @@ -670,7 +670,7 @@ IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, } -void dirax_cleanup(IndexingPrivate *pp) +void dirax_cleanup(void *pp) { struct dirax_private *p; p = (struct dirax_private *)pp; diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h index 9776f3f0..96ba6dbc 100644 --- a/libcrystfel/src/dirax.h +++ b/libcrystfel/src/dirax.h @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010,2012-2014 Thomas White <taw@physics.org> + * 2010,2012-2014,2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -39,13 +39,12 @@ extern "C" { #endif -extern int run_dirax(struct image *image, IndexingPrivate *ipriv); +extern int run_dirax(struct image *image, void *ipriv); -extern IndexingPrivate *dirax_prepare(IndexingMethod *indm, - UnitCell *cell, struct detector *det, - float *ltl); +extern void *dirax_prepare(IndexingMethod *indm, + UnitCell *cell, struct detector *det, float *ltl); -extern void dirax_cleanup(IndexingPrivate *pp); +extern void dirax_cleanup(void *pp); #ifdef __cplusplus } diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c index 25771a69..8e4eb861 100644 --- a/libcrystfel/src/events.c +++ b/libcrystfel/src/events.c @@ -3,10 +3,11 @@ * * Event properties * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: + * 2017 Thomas White * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -26,10 +27,8 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE - #include "events.h" +#include "utils.h" #include <hdf5.h> #include <string.h> @@ -585,46 +584,38 @@ char *retrieve_full_path(struct event *ev, const char *data) { int ei ; char *return_value; + char *pholder; return_value = strdup(data); + pholder = strstr(return_value,"%"); + ei = 0; - for ( ei=0; ei<ev->path_length; ei++ ) { + while ( pholder != NULL ) { char *tmp; - tmp = event_path_placeholder_subst(ev->path_entries[ei], - return_value); - free(return_value); - return_value = tmp; - - } - - return return_value; - -} - - -char *partial_event_substitution(struct event *ev, const char *data) -{ - int ei ; - char *return_value; - char *pholder; - - return_value = strdup(data); - pholder = strstr(return_value,"%"); - ei = 0; + /* Check we have enough things to put in the placeholders */ + if ( ei >= ev->path_length ) { + ERROR("Too many placeholders ('%%') in location.\n"); + return NULL; + } - while( pholder != NULL) { + /* Substitute one placeholder */ + tmp = event_path_placeholder_subst(ev->path_entries[ei++], + return_value); - char *tmp_subst_data; + if ( tmp == NULL ) { + ERROR("Couldn't substitute placeholder\n"); + return NULL; + } - tmp_subst_data = event_path_placeholder_subst(ev->path_entries[ei], - return_value); + /* Next time, we will substitute the next part of the path into + * the partially substituted string */ free(return_value); - return_value = strdup(tmp_subst_data); - free(tmp_subst_data); + return_value = tmp; + pholder = strstr(return_value, "%"); - ei += 1; + } return return_value; diff --git a/libcrystfel/src/events.h b/libcrystfel/src/events.h index 7f9c6731..4f717209 100644 --- a/libcrystfel/src/events.h +++ b/libcrystfel/src/events.h @@ -78,7 +78,6 @@ extern char *get_event_string(struct event *ev); extern struct event *get_event_from_event_string(const char *ev_string); extern char *event_path_placeholder_subst(const char *ev_name, const char *data); -extern char *partial_event_substitution(struct event *ev, const char *data); extern char *retrieve_full_path(struct event *ev, const char *data); diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c index 02d523c6..b531e3af 100644 --- a/libcrystfel/src/felix.c +++ b/libcrystfel/src/felix.c @@ -3,8 +3,8 @@ * * Invoke Felix for multi-crystal autoindexing. * - * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2015-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: * 2015 Thomas White <taw@physics.org> @@ -658,9 +658,8 @@ static void parse_options(const char *options, struct felix_private *gp) free(option); } -IndexingPrivate *felix_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl, - const char *options) +void *felix_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl, const char *options) { struct felix_private *gp; diff --git a/libcrystfel/src/felix.h b/libcrystfel/src/felix.h index 40771933..c06e963a 100644 --- a/libcrystfel/src/felix.h +++ b/libcrystfel/src/felix.h @@ -3,12 +3,12 @@ * * Invoke Felix for multi-crystal autoindexing * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White <taw@physics.org> - * 2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de> + * 2010-2013,2017 Thomas White <taw@physics.org> + * 2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de> * * This file is part of CrystFEL. * @@ -36,11 +36,9 @@ #include "cell.h" -extern IndexingPrivate *felix_prepare(IndexingMethod *indm, - UnitCell *cell, - struct detector *det, - float *ltl, - const char *options); +extern void *felix_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl, + const char *options); extern void felix_cleanup(IndexingPrivate *pp); diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 6542b360..13339506 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -44,6 +44,19 @@ #include "utils.h" +/** + * SECTION:hdf5-file + * @short_description: HDF5 file handling + * @title: HDF5 file handling + * @section_id: + * @see_also: + * @include: "hdf5-file.h" + * @Image: + * + * Routines for accessing HDF5 files. + **/ + + struct hdf5_write_location { const char *location; @@ -121,8 +134,7 @@ struct hdfile *hdfile_open(const char *filename) } -int hdfile_set_image(struct hdfile *f, const char *path, - struct panel *p) +int hdfile_set_image(struct hdfile *f, const char *path) { f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT); if ( f->dh < 0 ) { @@ -320,9 +332,34 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line) } -/* Get peaks from HDF5, in "CXI format" (as in "CXIDB") */ -int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe) +/** + * get_peaks_cxi_2: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @fpe: A %filename_plus_event structure specifying the event + * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates + * + * Get peaks from HDF5, in "CXI format" (as in "CXIDB"). The data should be in + * a set of arrays under @p. The number of peaks should be in a 1D array at + * @p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at + * @p/peakXPosRaw and @p/peakYPosRaw respectively (sorry about the naming). The + * first dimension of these arrays should be the event number (as given by + * @fpe). The intensities are expected to be at @p/peakTotalIntensity in a + * similar 2D array. + * + * CrystFEL considers all peak locations to be distances from the corner of the + * detector panel, in pixel units, consistent with its description of detector + * geometry (see 'man crystfel_geometry'). The software which generates the + * CXI files, including Cheetah, may instead consider the peak locations to be + * pixel indices in the data array. In this case, the peak coordinates should + * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero. + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe, int half_pixel_shift) { char path_n[1024]; char path_x[1024]; @@ -338,6 +375,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, float *buf_y; float *buf_i; + double peak_offset = half_pixel_shift ? 0.5 : 0.0; + if ( (fpe != NULL) && (fpe->ev != NULL) && (fpe->ev->dim_entries != NULL) ) { @@ -375,8 +414,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, float fs, ss, val; struct panel *p; - fs = buf_x[pk]; - ss = buf_y[pk]; + fs = buf_x[pk] + peak_offset; + ss = buf_y[pk] + peak_offset; val = buf_i[pk]; p = find_orig_panel(image->det, fs, ss); @@ -395,7 +434,53 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, } -int get_peaks(struct image *image, struct hdfile *f, const char *p) +/** + * get_peaks_cxi: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @fpe: A %filename_plus_event structure specifying the event + * + * This is a wrapper function to preserve API compatibility with older CrystFEL + * versions. Use get_peaks_cxi_2() instead. + * + * This function is equivalent to get_peaks_cxi_2(@image, @f, @p, @fpe, 1). + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe) +{ + return get_peaks_cxi_2(image, f, p, fpe, 1); +} + + +/** + * get_peaks_2: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates + * + * Get peaks from HDF5. The peak list should be located at @p in the HDF5 file, + * a 2D array where the first dimension equals the number of peaks and second + * dimension is three. The first two columns contain the fast scan and slow + * scan coordinates, respectively, of the peaks. The third column contains the + * intensities. + * + * CrystFEL considers all peak locations to be distances from the corner of the + * detector panel, in pixel units, consistent with its description of detector + * geometry (see 'man crystfel_geometry'). The software which generates the + * CXI files, including Cheetah, may instead consider the peak locations to be + * pixel indices in the data array. In this case, the peak coordinates should + * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero. + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks_2(struct image *image, struct hdfile *f, const char *p, + int half_pixel_shift) { hid_t dh, sh; hsize_t size[2]; @@ -405,6 +490,7 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) herr_t r; int tw; char *np; + double peak_offset = half_pixel_shift ? 0.5 : 0.0; if ( image->event != NULL ) { np = retrieve_full_path(image->event, p); @@ -473,8 +559,8 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) float fs, ss, val; struct panel *p; - fs = buf[tw*i+0]; - ss = buf[tw*i+1]; + fs = buf[tw*i+0] + peak_offset; + ss = buf[tw*i+1] + peak_offset; val = buf[tw*i+2]; p = find_orig_panel(image->det, fs, ss); @@ -499,6 +585,26 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p) } +/** + * get_peaks: + * @image: An %image structure + * @f: An %hdfile structure + * @p: The HDF5 path to the peak data + * + * This is a wrapper function to preserve API compatibility with older CrystFEL + * versions. Use get_peaks_2() instead. + * + * This function is equivalent to get_peaks_2(@image, @f, @p, 1). + * + * Returns: non-zero on error, zero otherwise. + * + */ +int get_peaks(struct image *image, struct hdfile *f, const char *p) +{ + return get_peaks_2(image, f, p, 1); +} + + static void cleanup(hid_t fh) { int n_ids, i; @@ -1127,6 +1233,9 @@ static int get_scalar_value(struct hdfile *f, const char *name, void *val, return 1; } + H5Tclose(type); + H5Dclose(dh); + return 0; } @@ -1154,7 +1263,7 @@ static int get_ev_based_value(struct hdfile *f, const char *name, char *subst_name = NULL; if ( ev->path_length != 0 ) { - subst_name = partial_event_substitution(ev, name); + subst_name = retrieve_full_path(ev, name); } else { subst_name = strdup(name); } @@ -1284,8 +1393,10 @@ int hdfile_get_value(struct hdfile *f, const char *name, } -void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, - struct event *ev, struct image *image) +static void hdfile_fill_in_beam_parameters(struct beam_params *beam, + struct hdfile *f, + struct event *ev, + struct image *image) { double eV; @@ -1298,8 +1409,8 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, int r; - r = hdfile_get_value(f, beam->photon_energy_from, ev, &eV, - H5T_NATIVE_DOUBLE); + r = hdfile_get_value(f, beam->photon_energy_from, + ev, &eV, H5T_NATIVE_DOUBLE); if ( r ) { ERROR("Failed to read '%s'\n", beam->photon_energy_from); @@ -1311,6 +1422,36 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f, } +static void hdfile_fill_in_clen(struct detector *det, struct hdfile *f, + struct event *ev) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + + double val; + int r; + + r = hdfile_get_value(f, p->clen_from, ev, &val, + H5T_NATIVE_DOUBLE); + if ( r ) { + ERROR("Failed to read '%s'\n", p->clen_from); + } else { + p->clen = val * 1.0e-3; + } + + } + + adjust_centering_for_rail(p); + + } +} + + int hdf5_read(struct hdfile *f, struct image *image, const char *element, int satcorr) { @@ -1326,7 +1467,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, if ( element == NULL ) { fail = hdfile_set_first_image(f, "/"); } else { - fail = hdfile_set_image(f, element, NULL); + fail = hdfile_set_image(f, element); } if ( fail ) { @@ -1375,7 +1516,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element, if ( image->beam != NULL ) { - fill_in_beam_parameters(image->beam, f, NULL, image); + hdfile_fill_in_beam_parameters(image->beam, f, NULL, image); if ( image->lambda > 1000 ) { /* Error message covers a silly value in the beam file @@ -1634,10 +1775,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( !exists ) { ERROR("Cannot find data for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } - fail = hdfile_set_image(f, panel_full_path, p); + fail = hdfile_set_image(f, panel_full_path); free(panel_full_path); @@ -1654,9 +1798,12 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( !exists ) { ERROR("Cannot find data for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } - fail = hdfile_set_image(f, p->data, p); + fail = hdfile_set_image(f, p->data); } @@ -1664,6 +1811,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( fail ) { ERROR("Couldn't select path for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1673,6 +1823,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, f_count = malloc(hsd->num_dims*sizeof(hsize_t)); if ( (f_offset == NULL) || (f_count == NULL ) ) { ERROR("Failed to allocate offset or count.\n"); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } for ( hsi=0; hsi<hsd->num_dims; hsi++ ) { @@ -1700,6 +1853,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( check < 0 ) { ERROR("Error selecting file dataspace for panel %s\n", p->name); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1713,6 +1869,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, ERROR("Failed to allocate panel %s\n", p->name); free(f_offset); free(f_count); + free(image->dp); + free(image->bad); + free(image->sat); return 1; } for ( i=0; i<p->w*p->h; i++ ) image->sat[pi][i] = INFINITY; @@ -1724,6 +1883,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, p->name); free(f_offset); free(f_count); + for ( i=0; i<=pi; i++ ) { + free(image->dp[i]); + free(image->sat[i]); + } + free(image->dp); + free(image->bad); + free(image->sat); return 1; } @@ -1740,7 +1906,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, if ( load_satmap(f, ev, p, f_offset, f_count, hsd, image->sat[pi]) ) { - ERROR("Failed to laod sat map for panel %s\n", + ERROR("Failed to load sat map for panel %s\n", p->name); } } @@ -1753,13 +1919,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, H5Dclose(f->dh); f->data_open = 0; - fill_in_values(image->det, f, ev); + hdfile_fill_in_clen(image->det, f, ev); if ( satcorr ) debodge_saturation(f, image); if ( image->beam != NULL ) { - fill_in_beam_parameters(image->beam, f, ev, image); + hdfile_fill_in_beam_parameters(image->beam, f, ev, image); if ( (image->lambda > 1.0) || (image->lambda < 1e-20) ) { @@ -1966,7 +2132,7 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name, char *tmp = NULL, *subst_name = NULL; if (ev != NULL && ev->path_length != 0 ) { - subst_name = partial_event_substitution(ev, name); + subst_name = retrieve_full_path(ev, name); } else { subst_name = strdup(name); } @@ -2166,7 +2332,7 @@ int hdfile_set_first_image(struct hdfile *f, const char *group) for ( i=0; i<n; i++ ) { if ( is_image[i] ) { - hdfile_set_image(f, names[i], NULL); + hdfile_set_image(f, names[i]); for ( j=0; j<n; j++ ) free(names[j]); free(is_image); free(is_group); @@ -2484,6 +2650,9 @@ static int check_dims(struct hdfile *hdfile, struct panel *p, struct event *ev, return 1; } + H5Sclose(sh); + H5Dclose(dh); + return 0; } diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h index 8c89eb93..ab53dd2e 100644 --- a/libcrystfel/src/hdf5-file.h +++ b/libcrystfel/src/hdf5-file.h @@ -3,11 +3,11 @@ * * Read/write HDF5 data files * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2015 Thomas White <taw@physics.org> + * 2009-2017 Thomas White <taw@physics.org> * 2014 Valerio Mariani * @@ -67,8 +67,8 @@ extern int hdf5_read2(struct hdfile *f, struct image *image, extern int check_path_existence(hid_t fh, const char *path); extern struct hdfile *hdfile_open(const char *filename); -int hdfile_set_image(struct hdfile *f, const char *path, - struct panel *p); +int hdfile_set_image(struct hdfile *f, const char *path); + extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile, int binning, int16_t *maxp); extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent, @@ -78,9 +78,16 @@ extern void hdfile_close(struct hdfile *f); extern int get_peaks(struct image *image, struct hdfile *f, const char *p); +extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p, + int half_pixel_shift); + extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, struct filename_plus_event *fpe); +extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, + struct filename_plus_event *fpe, + int half_pixel_shift); + extern struct copy_hdf5_field *new_copy_hdf5_field_list(void); extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f); diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 7bb4cede..aad5017c 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -3,12 +3,12 @@ * * Handle images and image features * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de> - * 2011-2016 Thomas White <taw@physics.org> + * 2011-2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -27,14 +27,28 @@ * */ -#define _ISOC99_SOURCE +#include <config.h> + #include <stdlib.h> #include <assert.h> #include <math.h> #include <stdio.h> +#include <hdf5.h> + +#ifdef HAVE_CBFLIB +#ifdef HAVE_CBFLIB_CBF_H +#include <cbflib/cbf.h> +#endif +#ifdef HAVE_CBF_CBF_H +#include <cbf/cbf.h> +#endif +#endif #include "image.h" #include "utils.h" +#include "events.h" +#include "hdf5-file.h" +#include "detector.h" /** * SECTION:image @@ -51,6 +65,14 @@ */ +struct imagefile +{ + enum imagefile_type type; + char *filename; + struct hdfile *hdfile; +}; + + struct _imagefeaturelist { struct imagefeature *features; @@ -300,3 +322,673 @@ void free_all_crystals(struct image *image) free(image->crystals); image->n_crystals = 0; } + + +/**************************** Image field lists *******************************/ + +struct imagefile_field_list +{ + char **fields; + int n_fields; + int max_fields; +}; + + +struct imagefile_field_list *new_imagefile_field_list() +{ + struct imagefile_field_list *n; + + n = calloc(1, sizeof(struct imagefile_field_list)); + if ( n == NULL ) return NULL; + + n->max_fields = 32; + n->fields = malloc(n->max_fields*sizeof(char *)); + if ( n->fields == NULL ) { + free(n); + return NULL; + } + + return n; +} + + +void free_imagefile_field_list(struct imagefile_field_list *n) +{ + int i; + for ( i=0; i<n->n_fields; i++ ) { + free(n->fields[i]); + } + free(n->fields); + free(n); +} + + +void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) +{ + int i; + + /* Already on the list? Don't re-add if so. */ + for ( i=0; i<copyme->n_fields; i++ ) { + if ( strcmp(copyme->fields[i], name) == 0 ) return; + } + + /* Need more space? */ + if ( copyme->n_fields == copyme->max_fields ) { + + char **nfields; + int nmax = copyme->max_fields + 32; + + nfields = realloc(copyme->fields, nmax*sizeof(char *)); + if ( nfields == NULL ) { + ERROR("Failed to allocate space for new HDF5 field.\n"); + return; + } + + copyme->max_fields = nmax; + copyme->fields = nfields; + + } + + copyme->fields[copyme->n_fields] = strdup(name); + if ( copyme->fields[copyme->n_fields] == NULL ) { + ERROR("Failed to add field for copying '%s'\n", name); + return; + } + + copyme->n_fields++; +} + + +/******************************* CBF files ************************************/ + +#ifdef HAVE_CBFLIB + +static char *cbf_strerr(int e) +{ + char *err; + + err = malloc(1024); + if ( err == NULL ) return NULL; + + err[0] = '\0'; + + /* NB Sum of lengths of all strings must be less than 1024 */ + if ( e & CBF_FORMAT ) strcat(err, "Invalid format"); + if ( e & CBF_ALLOC ) strcat(err, "Memory allocation failed"); + if ( e & CBF_ARGUMENT ) strcat(err, "Invalid argument"); + if ( e & CBF_ASCII ) strcat(err, "Value is ASCII"); + if ( e & CBF_BINARY ) strcat(err, "Value is binary"); + if ( e & CBF_BITCOUNT ) strcat(err, "Wrong number of bits"); + if ( e & CBF_ENDOFDATA ) strcat(err, "End of data"); + if ( e & CBF_FILECLOSE ) strcat(err, "File close error"); + if ( e & CBF_FILEOPEN ) strcat(err, "File open error"); + if ( e & CBF_FILEREAD ) strcat(err, "File read error"); + if ( e & CBF_FILETELL ) strcat(err, "File tell error"); + if ( e & CBF_FILEWRITE ) strcat(err, "File write error"); + if ( e & CBF_IDENTICAL ) strcat(err, "Name already exists"); + if ( e & CBF_NOTFOUND ) strcat(err, "Not found"); + if ( e & CBF_OVERFLOW ) strcat(err, "Overflow"); + if ( e & CBF_UNDEFINED ) strcat(err, "Number undefined"); + if ( e & CBF_NOTIMPLEMENTED ) strcat(err, "Not implemented"); + + return err; +} + + +static int unpack_panels(struct image *image, signed int *data, int data_width, + int data_height) +{ + int pi; + + /* FIXME: Load these masks from an HDF5 file, if filenames are + * given in the geometry file */ + uint16_t *flags = NULL; + float *sat = NULL; + + image->dp = malloc(image->det->n_panels * sizeof(float *)); + image->bad = malloc(image->det->n_panels * sizeof(int *)); + image->sat = malloc(image->det->n_panels * sizeof(float *)); + if ( (image->dp == NULL) || (image->bad == NULL) + || (image->sat == NULL) ) + { + ERROR("Failed to allocate panels.\n"); + return 1; + } + + for ( pi=0; pi<image->det->n_panels; pi++ ) { + + struct panel *p; + int fs, ss; + + p = &image->det->panels[pi]; + image->dp[pi] = malloc(p->w*p->h*sizeof(float)); + image->bad[pi] = calloc(p->w*p->h, sizeof(int)); + image->sat[pi] = malloc(p->w*p->h*sizeof(float)); + if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) + || (image->sat[pi] == NULL) ) + { + ERROR("Failed to allocate panel\n"); + return 1; + } + + if ( p->mask != NULL ) { + ERROR("WARNING: Bad pixel masks do not currently work " + "with CBF files\n"); + ERROR(" (bad pixel regions specified in the geometry " + "file will be used, however)\n"); + } + + if ( p->satmap != NULL ) { + ERROR("WARNING: Saturation maps do not currently work " + "with CBF files\n"); + } + + if ( (p->orig_min_fs + p->w > data_width) + || (p->orig_min_ss + p->h > data_height) ) + { + ERROR("Panel %s is outside range of data in CBF file\n", + p->name); + return 1; + } + + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { + + int idx; + int cfs, css; + int bad = 0; + + cfs = fs+p->orig_min_fs; + css = ss+p->orig_min_ss; + idx = cfs + css*data_width; + + image->dp[pi][fs+p->w*ss] = data[idx]; + + if ( sat != NULL ) { + image->sat[pi][fs+p->w*ss] = sat[idx]; + } else { + image->sat[pi][fs+p->w*ss] = INFINITY; + } + + if ( p->no_index ) bad = 1; + + if ( in_bad_region(image->det, p, cfs, css) ) { + bad = 1; + } + + if ( flags != NULL ) { + + int f; + + f = flags[idx]; + + /* Bad if it's missing any of the "good" bits */ + if ( (f & image->det->mask_good) + != image->det->mask_good ) bad = 1; + + /* Bad if it has any of the "bad" bits. */ + if ( f & image->det->mask_bad ) bad = 1; + + } + image->bad[pi][fs+p->w*ss] = bad; + + } + } + + } + + return 0; +} + + +static void cbf_fill_in_beam_parameters(struct beam_params *beam, + struct imagefile *f, + struct image *image) +{ + double eV; + + if ( beam->photon_energy_from == NULL ) { + + /* Explicit value given */ + eV = beam->photon_energy; + + } else { + + ERROR("Can't get photon energy from CBF yet.\n"); + eV = 0.0; + + } + + image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale; +} + + +static void cbf_fill_in_clen(struct detector *det, struct imagefile *f) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + + ERROR("Can't get clen from CBF yet.\n"); + + } + + adjust_centering_for_rail(p); + + } +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + cbf_handle cbfh; + FILE *fh; + int r; + unsigned int compression; + int binary_id, minelement, maxelement, elsigned, elunsigned; + size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; + const char *byteorder; + signed int *data; + + if ( image->det == NULL ) { + ERROR("read_cbf() needs a geometry\n"); + return 1; + } + + if ( cbf_make_handle(&cbfh) ) { + ERROR("Failed to allocate CBF handle\n"); + return 1; + } + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file '%s'\n", f->filename); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + unpack_panels(image, data, dimfast, dimmid); + free(data); + + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + + cbf_free_handle(cbfh); + return 0; +} + + +static float *convert_float(signed int *data, int w, int h) +{ + float *df; + long int i; + + df = malloc(sizeof(float)*w*h); + if ( df == NULL ) return NULL; + + for ( i=0; i<w*h; i++ ) { + df[i] = data[i]; + } + + return df; +} + + +static int read_cbf_simple(struct imagefile *f, struct image *image) +{ + cbf_handle cbfh; + FILE *fh; + int r; + unsigned int compression; + int binary_id, minelement, maxelement, elsigned, elunsigned; + size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; + const char *byteorder; + signed int *data; + + if ( cbf_make_handle(&cbfh) ) { + ERROR("Failed to allocate CBF handle\n"); + return 1; + } + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file '%s'\n", f->filename); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + image->det = simple_geometry(image, dimfast, dimmid); + image->dp = malloc(sizeof(float *)); + if ( image->dp == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + image->dp[0] = convert_float(data, dimfast, dimmid); + if ( image->dp[0] == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + + cbf_free_handle(cbfh); + return 0; +} + +#else /* HAVE_CBFLIB */ + +static int read_cbf_simple(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +#endif /* HAVE_CBFLIB */ + + +/****************************** Image files ***********************************/ + + +static signed int is_cbf_file(const char *filename) +{ + FILE *fh; + char line[1024]; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return -1; + + if ( fgets(line, 1024, fh) == NULL ) return -1; + fclose(fh); + + if ( strstr(line, "CBF") == NULL ) { + return 0; + } + + return 1; +} + + +struct imagefile *imagefile_open(const char *filename) +{ + struct imagefile *f; + + f = malloc(sizeof(struct imagefile)); + if ( f == NULL ) return NULL; + + if ( H5Fis_hdf5(filename) > 0 ) { + + /* This is an HDF5, pass through to HDF5 layer */ + f->type = IMAGEFILE_HDF5; + f->hdfile = hdfile_open(filename); + + if ( f->hdfile == NULL ) { + free(f); + return NULL; + } + + } else if ( is_cbf_file(filename) > 0 ) { + + f->type = IMAGEFILE_CBF; + + } + + f->filename = strdup(filename); + return f; +} + + +int imagefile_read(struct imagefile *f, struct image *image, + struct event *event) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read2(f->hdfile, image, event, 0); + } else if ( f->type == IMAGEFILE_CBF ) { + return read_cbf(f, image); + } else { + ERROR("Unknown file type %i\n", f->type); + return 1; + } +} + + +/* Read a simple file, no multi-event, no prior geometry etc, and + * generate a geometry for it */ +int imagefile_read_simple(struct imagefile *f, struct image *image) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read(f->hdfile, image, NULL, 0); + } else { + return read_cbf_simple(f, image); + } +} + + +enum imagefile_type imagefile_get_type(struct imagefile *f) +{ + assert(f != NULL); + return f->type; +} + + +struct hdfile *imagefile_get_hdfile(struct imagefile *f) +{ + if ( f->type != IMAGEFILE_HDF5 ) { + ERROR("Not an HDF5 file!\n"); + return NULL; + } + + return f->hdfile; +} + + +void imagefile_copy_fields(struct imagefile *f, + const struct imagefile_field_list *copyme, + FILE *fh, struct event *ev) +{ + int i; + + if ( copyme == NULL ) return; + + for ( i=0; i<copyme->n_fields; i++ ) { + + char *val; + char *field; + + field = copyme->fields[i]; + + if ( f->type == IMAGEFILE_HDF5 ) { + val = hdfile_get_string_value(f->hdfile, field, ev); + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + free(val); + + } else { + STATUS("Mock CBF variable\n"); + fprintf(fh, "cbf/%s = %s\n", field, "(FIXME)"); + } + + + } +} + + +void imagefile_close(struct imagefile *f) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + hdfile_close(f->hdfile); + } + free(f->filename); + free(f); +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 9fd9b495..9719bb59 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -3,11 +3,11 @@ * * Handle images and image features * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2009-2016 Thomas White <taw@physics.org> + * 2009-2017 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * @@ -44,6 +44,8 @@ struct detector; struct imagefeature; struct sample; struct image; +struct imagefile; +struct imagefile_field_list; #include "utils.h" #include "cell.h" @@ -51,6 +53,7 @@ struct image; #include "reflist.h" #include "crystal.h" #include "index.h" +#include "events.h" /** * SpectrumType: @@ -67,7 +70,26 @@ typedef enum { SPECTRUM_TWOCOLOUR } SpectrumType; -/* Structure describing a feature in an image */ + +/** + * imagefeature: + * @parent: Image this feature belongs to + * @fs: Fast scan coordinate + * @ss: Slow scan coordinate + * @p: Pointer to panel + * @intensity: Intensity of peak + * @rx: Reciprocal x coordinate in m^-1 + * @ry: Reciprocal y coordinate in m^-1 + * @rz: Reciprocal z coordinate in m^-1 + * @name: Text name for feature + * + * Represents a peak in an image. + * + * Note carefully that the @fs and @ss coordinates are the distances, measured + * in pixels, from the corner of the panel. They are NOT pixel indices. + * If the peak is in the middle of the first pixel, its coordinates would be + * 0.5,0.5. + */ struct imagefeature { struct image *parent; @@ -81,12 +103,21 @@ struct imagefeature { double ry; double rz; - /* Internal use only */ + const char *name; + + /*< private >*/ int valid; +}; - const char *name; + +/* An enum representing the image file formats we can handle */ +enum imagefile_type +{ + IMAGEFILE_HDF5, + IMAGEFILE_CBF }; + /* An opaque type representing a list of image features */ typedef struct _imagefeaturelist ImageFeatureList; @@ -98,44 +129,46 @@ struct sample }; +/** + * beam_params: + * @photon_energy: eV per photon + * @photon_energy_from: HDF5 dataset name + * @photon_energy_scale: Scale factor for photon energy, if it comes from HDF5 + */ struct beam_params { - double photon_energy; /* eV per photon */ - char *photon_energy_from; /* HDF5 dataset name */ - double photon_energy_scale; /* Scale factor for photon energy, if the - * energy is to be from the HDF5 file */ + double photon_energy; + char *photon_energy_from; + double photon_energy_scale; }; /** * image: - * - * <programlisting> - * struct image - * { - * Crystal **crystals; - * int n_crystals; - * IndexingMethod indexed_by; - * - * struct detector *det; - * struct beam_params *beam; - * char *filename; - * const struct copy_hdf5_field *copyme; - * - * int id; - * - * double lambda; - * double div; - * double bw; - * - * int width; - * int height; - * - * long long int num_peaks; - * long long int num_saturated_peaks; - * ImageFeatureList *features; - * }; - * </programlisting> + * @crystals: Array of crystals in the image + * @n_crystals: The number of crystals in the image + * @indexed_by: Indexing method which indexed this pattern + * @det: Detector structure + * @beam: Beam parameters structure + * @filename: Filename for the image file + * @copyme: Fields to copy from the image file to the stream + * @id: ID number of the thread handling this image + * @serial: Serial number for this image + * @lambda: Wavelength + * @div: Divergence + * @bw: Bandwidth + * @num_peaks: The number of peaks + * @num_saturated_peaks: The number of saturated peaks + * @features: The peaks found in the image + * @dp: The image data, by panel + * @bad: The bad pixel mask, array by panel + * @sat: The per-pixel saturation mask, array by panel + * @event: Event ID for the image + * @stuff_from_stream: Items read back from the stream + * @avg_clen: Mean of camera length values for all panels + * @spectrum: Spectrum information + * @nsamples: Number of spectrum samples + * @spectrum_size: SIze of spectrum array * * The field <structfield>data</structfield> contains the raw image data, if it * is currently available. The data might be available throughout the @@ -152,8 +185,8 @@ struct beam_params * returned by the low-level indexing system. <structfield>n_crystals</structfield> * is the number of crystals which were found in the image. * - * <structfield>copyme</structfield> represents a list of HDF5 fields to copy - * to the output stream. + * <structfield>copyme</structfield> represents a list of fields in the image + * file (e.g. HDF5 fields or CBF headers) to copy to the output stream. **/ struct image; @@ -171,7 +204,7 @@ struct image { struct beam_params *beam; /* The nominal beam parameters */ char *filename; struct event *event; - const struct copy_hdf5_field *copyme; + const struct imagefile_field_list *copyme; struct stuff_from_stream *stuff_from_stream; double avg_clen; /* Average camera length extracted @@ -192,8 +225,8 @@ struct image { double bw; /* Bandwidth as a fraction */ /* Detected peaks */ - long long int num_peaks; - long long int num_saturated_peaks; + long long num_peaks; + long long num_saturated_peaks; ImageFeatureList *features; }; @@ -233,6 +266,25 @@ extern void image_add_crystal(struct image *image, Crystal *cryst); extern void remove_flagged_crystals(struct image *image); extern void free_all_crystals(struct image *image); +/* Image files */ +extern struct imagefile *imagefile_open(const char *filename); +extern int imagefile_read(struct imagefile *f, struct image *image, + struct event *event); +extern int imagefile_read_simple(struct imagefile *f, struct image *image); +extern struct hdfile *imagefile_get_hdfile(struct imagefile *f); +extern enum imagefile_type imagefile_get_type(struct imagefile *f); +extern void imagefile_copy_fields(struct imagefile *f, + const struct imagefile_field_list *copyme, + FILE *fh, struct event *ev); +extern void imagefile_close(struct imagefile *f); + +/* Field lists */ +extern struct imagefile_field_list *new_imagefile_field_list(void); +extern void free_imagefile_field_list(struct imagefile_field_list *f); + +extern void add_imagefile_field(struct imagefile_field_list *copyme, + const char *name); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index e7fd879d..bdabb062 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -59,6 +59,17 @@ #include "taketwo.h" +struct _indexingprivate +{ + UnitCell *target_cell; + float tolerance[4]; + + int n_methods; + IndexingMethod *methods; + void **engine_private; +}; + + static int debug_index(struct image *image) { Crystal *cr = crystal_new(); @@ -72,78 +83,119 @@ static int debug_index(struct image *image) } -IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl, - const char *options) +static void *prepare_method(IndexingMethod *m, UnitCell *cell, + struct detector *det, float *ltl, + const char *options) { - int n; - int nm = 0; - IndexingPrivate **iprivs; + char *str; + IndexingMethod in = *m; + void *priv = NULL; - while ( indm[nm] != INDEXING_NONE ) nm++; - iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); + switch ( *m & INDEXING_METHOD_MASK ) { - for ( n=0; n<nm; n++ ) { + case INDEXING_NONE : + priv = "none"; + break; - int i; - IndexingMethod in; - char *str; + case INDEXING_DIRAX : + priv = dirax_prepare(m, cell, det, ltl); + break; - in = indm[n]; + case INDEXING_ASDF : + priv = asdf_prepare(m, cell, det, ltl); + break; - switch ( indm[n] & INDEXING_METHOD_MASK ) { + case INDEXING_MOSFLM : + priv = mosflm_prepare(m, cell, det, ltl); + break; - case INDEXING_DIRAX : - iprivs[n] = dirax_prepare(&indm[n], cell, det, ltl); - break; + case INDEXING_XDS : + priv = xds_prepare(m, cell, det, ltl); + break; - case INDEXING_ASDF : - iprivs[n] = asdf_prepare(&indm[n], cell, det, ltl); - break; + case INDEXING_DEBUG : + priv = (IndexingPrivate *)strdup("Hello!"); + break; - case INDEXING_MOSFLM : - iprivs[n] = mosflm_prepare(&indm[n], cell, det, ltl); - break; + case INDEXING_FELIX : + priv = felix_prepare(m, cell, det, ltl, options); + break; - case INDEXING_XDS : - iprivs[n] = xds_prepare(&indm[n], cell, det, ltl); - break; + case INDEXING_TAKETWO : + priv = taketwo_prepare(m, cell, det, ltl); + break; - case INDEXING_DEBUG : - iprivs[n] = (IndexingPrivate *)strdup("Hello!"); - break; + default : + ERROR("Don't know how to prepare indexing method %i\n", *m); + break; - case INDEXING_FELIX : - iprivs[n] = felix_prepare(&indm[n], cell, det, ltl, - options); - break; + } - case INDEXING_TAKETWO : - iprivs[n] = taketwo_prepare(&indm[n], cell, det, ltl); - break; + str = indexer_str(*m); - default : - ERROR("Don't know how to prepare indexing method %i\n", - indm[n]); - break; + if ( priv == NULL ) { + ERROR("Failed to prepare indexing method %s\n", str); + free(str); + return NULL; + } - } + if ( *m != INDEXING_NONE ) { + STATUS("Prepared indexing method %s\n", str); + } + free(str); - if ( iprivs[n] == NULL ) return NULL; + if ( in != *m ) { + ERROR("Note: flags were altered to take into account " + "the information provided and/or the limitations " + "of the indexing method.\nPlease check the " + "methods listed above carefully.\n"); + } - str = indexer_str(indm[n]); - STATUS("Prepared indexing method %i: %s\n", n, str); - free(str); + return priv; +} + + +IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, + struct detector *det, float *ltl, + int no_refine, const char *options) +{ + int i, n; + char **method_strings; + IndexingPrivate *ipriv; + + /* Parse indexing methods */ + n = assplode(method_list, ",", &method_strings, ASSPLODE_NONE); + + IndexingMethod *methods = malloc(n * sizeof(IndexingMethod)); + if ( methods == NULL ) { + ERROR("Failed to allocate indexing method list\n"); + return NULL; + } - if ( in != indm[n] ) { - ERROR("Note: flags were altered to take into account " - "the information provided and/or the limitations " - "of the indexing method.\nPlease check the " - "methods listed above carefully.\n"); + for ( i=0; i<n; i++ ) { + methods[i] = get_indm_from_string(method_strings[i]); + if ( methods[i] == INDEXING_ERROR ) { + free(methods); + return NULL; } + } + + ipriv = malloc(sizeof(struct _indexingprivate)); + if ( ipriv == NULL ) { + ERROR("Failed to allocate indexing data\n"); + return NULL; + } + + ipriv->engine_private = malloc((n+1) * sizeof(void *)); - for ( i=0; i<n; i++ ) { - if ( indm[i] == indm[n] ) { + for ( i=0; i<n; i++ ) { + + int j; + + ipriv->engine_private[i] = prepare_method(&methods[i], cell, + det, ltl, options); + for ( j=0; j<i; j++ ) { + if ( methods[i] == methods[j] ) { ERROR("Duplicate indexing method.\n"); ERROR("Have you specified some flags which " "aren't accepted by one of your " @@ -152,59 +204,66 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, } } + } + + ipriv->methods = methods; + ipriv->n_methods = n; + if ( cell != NULL ) { + ipriv->target_cell = cell_new_from_cell(cell); + } else { + ipriv->target_cell = NULL; } - iprivs[n] = NULL; + for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i]; - return iprivs; + return ipriv; } -void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) +void cleanup_indexing(IndexingPrivate *ipriv) { - int n = 0; + int n; - if ( indms == NULL ) return; /* Nothing to do */ - if ( privs == NULL ) return; /* Nothing to do */ + if ( ipriv == NULL ) return; /* Nothing to do */ - while ( indms[n] != INDEXING_NONE ) { + for ( n=0; n<ipriv->n_methods; n++ ) { - switch ( indms[n] & INDEXING_METHOD_MASK ) { + switch ( ipriv->methods[n] & INDEXING_METHOD_MASK ) { case INDEXING_NONE : break; case INDEXING_DIRAX : - dirax_cleanup(privs[n]); + dirax_cleanup(ipriv->engine_private[n]); break; case INDEXING_ASDF : - asdf_cleanup(privs[n]); + asdf_cleanup(ipriv->engine_private[n]); break; case INDEXING_MOSFLM : - mosflm_cleanup(privs[n]); + mosflm_cleanup(ipriv->engine_private[n]); break; case INDEXING_XDS : - xds_cleanup(privs[n]); + xds_cleanup(ipriv->engine_private[n]); break; case INDEXING_FELIX : - felix_cleanup(privs[n]); + felix_cleanup(ipriv->engine_private[n]); break; case INDEXING_DEBUG : - free(privs[n]); + free(ipriv->engine_private[n]); break; case INDEXING_TAKETWO : - taketwo_cleanup(privs[n]); + taketwo_cleanup(ipriv->engine_private[n]); break; default : ERROR("Don't know how to clean up indexing method %i\n", - indms[n]); + ipriv->methods[n]); break; } @@ -213,8 +272,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) } - free(indms); - free(privs); + free(ipriv->methods); + free(ipriv->engine_private); + cell_free(ipriv->target_cell); + free(ipriv); } @@ -243,7 +304,7 @@ void map_all_peaks(struct image *image) /* Return non-zero for "success" */ static int try_indexer(struct image *image, IndexingMethod indm, - IndexingPrivate *ipriv) + IndexingPrivate *ipriv, void *mpriv) { int i, r; int n_bad = 0; @@ -254,19 +315,19 @@ static int try_indexer(struct image *image, IndexingMethod indm, return 0; case INDEXING_DIRAX : - r = run_dirax(image, ipriv); + r = run_dirax(image, mpriv); break; case INDEXING_ASDF : - r = run_asdf(image, ipriv); + r = run_asdf(image, mpriv); break; case INDEXING_MOSFLM : - r = run_mosflm(image, ipriv); + r = run_mosflm(image, mpriv); break; case INDEXING_XDS : - r = run_xds(image, ipriv); + r = run_xds(image, mpriv); break; case INDEXING_DEBUG : @@ -274,11 +335,11 @@ static int try_indexer(struct image *image, IndexingMethod indm, break; case INDEXING_FELIX : - r = felix_index(image, ipriv); + r = felix_index(image, mpriv); break; case INDEXING_TAKETWO : - r = taketwo_index(image, ipriv); + r = taketwo_index(image, mpriv); break; default : @@ -296,11 +357,35 @@ static int try_indexer(struct image *image, IndexingMethod indm, crystal_set_mosaicity(cr, 0.0); /* Prediction refinement if requested */ - if ( (indm & INDEXING_REFINE) - && (refine_prediction(image, cr) != 0) ) - { - crystal_set_user_flag(cr, 1); - n_bad++; + if ( indm & INDEXING_REFINE ) { + + UnitCell *out; + + if ( refine_prediction(image, cr) ) { + crystal_set_user_flag(cr, 1); + n_bad++; + continue; + } + + if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS) + || (indm & INDEXING_CHECK_CELL_AXES) ) + { + + /* Check that the cell parameters are still + * within the tolerance */ + out = match_cell(crystal_get_cell(cr), + ipriv->target_cell, 0, + ipriv->tolerance, 0); + + if ( out == NULL ) { + crystal_set_user_flag(cr, 1); + n_bad++; + } + + cell_free(out); + + } + } } @@ -434,20 +519,18 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image) } } -void index_pattern(struct image *image, - IndexingMethod *indms, IndexingPrivate **iprivs) +void index_pattern(struct image *image, IndexingPrivate *ipriv) { - index_pattern_2(image, indms, iprivs, NULL); + index_pattern_2(image, ipriv, NULL); } -void index_pattern_2(struct image *image, IndexingMethod *indms, - IndexingPrivate **iprivs, int *ping) +void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) { int n = 0; ImageFeatureList *orig; - if ( indms == NULL ) return; + if ( ipriv == NULL ) return; map_all_peaks(image); image->crystals = NULL; @@ -455,7 +538,7 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, orig = image->features; - while ( indms[n] != INDEXING_NONE ) { + for ( n=0; n<ipriv->n_methods; n++ ) { int done = 0; int r; @@ -466,10 +549,11 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, do { - r = try_indexer(image, indms[n], iprivs[n]); + r = try_indexer(image, ipriv->methods[n], + ipriv, ipriv->engine_private[n]); success += r; ntry++; - done = finished_retry(indms[n], r, image); + done = finished_retry(ipriv->methods[n], r, image); if ( ntry > 5 ) done = 1; if ( ping != NULL ) (*ping)++; @@ -481,11 +565,14 @@ void index_pattern_2(struct image *image, IndexingMethod *indms, * crystals with a different indexing method) */ if ( success ) break; - n++; + } + if ( n < ipriv->n_methods ) { + image->indexed_by = ipriv->methods[n]; + } else { + image->indexed_by = INDEXING_NONE; } - image->indexed_by = indms[n]; image->features = orig; } @@ -656,100 +743,123 @@ char *indexer_str(IndexingMethod indm) } -IndexingMethod *build_indexer_list(const char *str) +static IndexingMethod warn_method(const char *str) +{ + ERROR("Indexing method must contain exactly one engine name: '%s'\n", + str); + return INDEXING_ERROR; +} + + +IndexingMethod get_indm_from_string(const char *str) { int n, i; - char **methods; - IndexingMethod *list; - int nmeth = 0; + char **bits; + IndexingMethod method = INDEXING_NONE; + int have_method = 0; - n = assplode(str, ",-", &methods, ASSPLODE_NONE); - list = malloc((n+1)*sizeof(IndexingMethod)); + n = assplode(str, "-", &bits, ASSPLODE_NONE); - nmeth = -1; /* So that the first method is #0 */ for ( i=0; i<n; i++ ) { - if ( strcmp(methods[i], "dirax") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_DIRAX; + if ( strcmp(bits[i], "dirax") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_DIRAX; + have_method = 1; - } else if ( strcmp(methods[i], "asdf") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_ASDF; + } else if ( strcmp(bits[i], "asdf") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_ASDF; + have_method = 1; - } else if ( strcmp(methods[i], "mosflm") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_MOSFLM; + } else if ( strcmp(bits[i], "mosflm") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_MOSFLM; + have_method = 1; - } else if ( strcmp(methods[i], "xds") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_XDS; + } else if ( strcmp(bits[i], "xds") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_XDS; + have_method = 1; - } else if ( strcmp(methods[i], "felix") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_FELIX; + } else if ( strcmp(bits[i], "felix") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_FELIX; + have_method = 1; - } else if ( strcmp(methods[i], "taketwo") == 0) { - list[++nmeth] = INDEXING_DEFAULTS_TAKETWO; + } else if ( strcmp(bits[i], "taketwo") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_TAKETWO; + have_method = 1; - } else if ( strcmp(methods[i], "none") == 0) { - list[++nmeth] = INDEXING_NONE; + } else if ( strcmp(bits[i], "none") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_NONE; + have_method = 1; - } else if ( strcmp(methods[i], "simulation") == 0) { - list[++nmeth] = INDEXING_SIMULATION; - return list; + } else if ( strcmp(bits[i], "simulation") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_SIMULATION; + return method; - } else if ( strcmp(methods[i], "debug") == 0) { - list[++nmeth] = INDEXING_DEBUG; - return list; + } else if ( strcmp(bits[i], "debug") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEBUG; + return method; - } else if ( strcmp(methods[i], "raw") == 0) { - list[nmeth] = set_raw(list[nmeth]); + } else if ( strcmp(bits[i], "raw") == 0) { + method = set_raw(method); - } else if ( strcmp(methods[i], "bad") == 0) { - list[nmeth] = set_bad(list[nmeth]); + } else if ( strcmp(bits[i], "bad") == 0) { + method = set_bad(method); - } else if ( strcmp(methods[i], "comb") == 0) { - list[nmeth] = set_comb(list[nmeth]); /* Default */ + } else if ( strcmp(bits[i], "comb") == 0) { + method = set_comb(method); /* Default */ - } else if ( strcmp(methods[i], "axes") == 0) { - list[nmeth] = set_axes(list[nmeth]); + } else if ( strcmp(bits[i], "axes") == 0) { + method = set_axes(method); - } else if ( strcmp(methods[i], "latt") == 0) { - list[nmeth] = set_lattice(list[nmeth]); + } else if ( strcmp(bits[i], "latt") == 0) { + method = set_lattice(method); - } else if ( strcmp(methods[i], "nolatt") == 0) { - list[nmeth] = set_nolattice(list[nmeth]); + } else if ( strcmp(bits[i], "nolatt") == 0) { + method = set_nolattice(method); - } else if ( strcmp(methods[i], "cell") == 0) { - list[nmeth] = set_cellparams(list[nmeth]); + } else if ( strcmp(bits[i], "cell") == 0) { + method = set_cellparams(method); - } else if ( strcmp(methods[i], "nocell") == 0) { - list[nmeth] = set_nocellparams(list[nmeth]); + } else if ( strcmp(bits[i], "nocell") == 0) { + method = set_nocellparams(method); - } else if ( strcmp(methods[i], "retry") == 0) { - list[nmeth] |= INDEXING_RETRY; + } else if ( strcmp(bits[i], "retry") == 0) { + method |= INDEXING_RETRY; - } else if ( strcmp(methods[i], "noretry") == 0) { - list[nmeth] &= ~INDEXING_RETRY; + } else if ( strcmp(bits[i], "noretry") == 0) { + method &= ~INDEXING_RETRY; - } else if ( strcmp(methods[i], "multi") == 0) { - list[nmeth] |= INDEXING_MULTI; + } else if ( strcmp(bits[i], "multi") == 0) { + method |= INDEXING_MULTI; - } else if ( strcmp(methods[i], "nomulti") == 0) { - list[nmeth] &= ~INDEXING_MULTI; + } else if ( strcmp(bits[i], "nomulti") == 0) { + method &= ~INDEXING_MULTI; - } else if ( strcmp(methods[i], "refine") == 0) { - list[nmeth] |= INDEXING_REFINE; + } else if ( strcmp(bits[i], "refine") == 0) { + method |= INDEXING_REFINE; - } else if ( strcmp(methods[i], "norefine") == 0) { - list[nmeth] &= ~INDEXING_REFINE; + } else if ( strcmp(bits[i], "norefine") == 0) { + method &= ~INDEXING_REFINE; } else { ERROR("Bad list of indexing methods: '%s'\n", str); - return NULL; + return INDEXING_ERROR; } - free(methods[i]); + free(bits[i]); } - free(methods); - list[++nmeth] = INDEXING_NONE; + free(bits); + + if ( !have_method ) return warn_method(str); - return list; + return method; } diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 33d568f6..2107ed80 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -81,6 +81,7 @@ * @INDEXING_DEBUG: Results injector for debugging * @INDEXING_ASDF: Use in-built "asdf" indexer * @INDEXING_TAKETWO: Use in-built "taketwo" indexer + * @INDEXING_ERROR: Special value for unrecognised indexing engine name * @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell * axes for agreement with given cell. * @INDEXING_CHECK_CELL_AXES: Check unit cell axes for agreement with given @@ -90,7 +91,7 @@ * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to * guide the indexing process. * @INDEXING_USE_CELL_PARAMETERS: Use the unit cell parameters to guide the - * indexingprocess. + * indexing process. * @INDEXING_RETRY: If the indexer doesn't succeed, delete the weakest peaks * and try again. * @INDEXING_MULTI: If the indexer succeeds, delete the peaks explained by the @@ -115,6 +116,8 @@ typedef enum { INDEXING_ASDF = 8, INDEXING_TAKETWO = 9, + INDEXING_ERROR = 255, /* Unrecognised indexing engine */ + /* Bits at the top of the IndexingMethod are flags which modify the * behaviour of the indexer. */ INDEXING_CHECK_CELL_COMBINATIONS = 256, @@ -143,28 +146,28 @@ extern "C" { * IndexingPrivate: * * This is an opaque data structure containing information needed by the - * indexing method. + * indexing system. **/ -typedef void *IndexingPrivate; +typedef struct _indexingprivate IndexingPrivate; -extern IndexingMethod *build_indexer_list(const char *str); +/* Convert indexing methods to/from text */ extern char *indexer_str(IndexingMethod indm); +extern IndexingMethod get_indm_from_string(const char *method); #include "detector.h" #include "cell.h" #include "image.h" -extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl, - const char *options); +extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell, + struct detector *det, float *ltl, + int no_refine, const char *options); -extern void index_pattern(struct image *image, - IndexingMethod *indms, IndexingPrivate **iprivs); +extern void index_pattern(struct image *image, IndexingPrivate *ipriv); -extern void index_pattern_2(struct image *image, IndexingMethod *indms, - IndexingPrivate **iprivs, int *ping); +extern void index_pattern_2(struct image *image, IndexingPrivate *ipriv, + int *ping); -extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs); +extern void cleanup_indexing(IndexingPrivate *ipriv); #ifdef __cplusplus } diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 05f853eb..a8e6e119 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -725,7 +725,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) } -int run_mosflm(struct image *image, IndexingPrivate *ipriv) +int run_mosflm(struct image *image, void *ipriv) { struct mosflm_data *mosflm; unsigned int opts; @@ -843,8 +843,8 @@ int run_mosflm(struct image *image, IndexingPrivate *ipriv) } -IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct mosflm_private *mp; int need_cell = 0; @@ -911,7 +911,7 @@ IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, } -void mosflm_cleanup(IndexingPrivate *pp) +void mosflm_cleanup(void *pp) { struct mosflm_private *p; p = (struct mosflm_private *)pp; diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h index 572741dd..39ec5390 100644 --- a/libcrystfel/src/mosflm.h +++ b/libcrystfel/src/mosflm.h @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian <rkirian@asu.edu> - * 2012-2014 Thomas White <taw@physics.org> + * 2012-2014,2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -41,12 +41,12 @@ extern "C" { #endif -extern int run_mosflm(struct image *image, IndexingPrivate *ipriv); +extern int run_mosflm(struct image *image, void *ipriv); -extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl); +extern void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl); -extern void mosflm_cleanup(IndexingPrivate *pp); +extern void mosflm_cleanup(void *pp); #ifdef __cplusplus } diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c new file mode 100644 index 00000000..417465df --- /dev/null +++ b/libcrystfel/src/peakfinder8.c @@ -0,0 +1,1195 @@ +/* + * peakfinder8.c + * + * The peakfinder8 algorithm + * + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017 Valerio Mariani <valerio.mariani@desy.de> + * 2017 Anton Barty <anton.barty@desy.de> + * 2017 Oleksandr Yefanov <oleksandr.yefanov@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <stdlib.h> + +#include "peakfinder8.h" + + +// CrystFEL-only block 1 +struct radius_maps +{ + float **r_maps; + int n_rmaps; +}; + + +struct peakfinder_mask +{ + char **masks; + int n_masks; +}; + + +struct peakfinder_panel_data +{ + float **panel_data; + int *panel_h; + int *panel_w; + int num_panels; +}; +// End of CrystFEL-only block 1 + + +struct radial_stats +{ + float *roffset; + float *rthreshold; + float *lthreshold; + float *rsigma; + int *rcount; + int n_rad_bins; +}; + + +struct peakfinder_intern_data +{ + char *pix_in_peak_map; + int *infs; + int *inss; + int *peak_pixels; +}; + + +struct peakfinder_peak_data +{ + int num_found_peaks; + int *npix; + float *com_fs; + float *com_ss; + int *com_index; + float *tot_i; + float *max_i; + float *sigma; + float *snr; +}; + + +// CrystFEL-only block 2 +static struct radius_maps *compute_radius_maps(struct detector *det) +{ + int i, u, iss, ifs; + struct panel p; + struct radius_maps *rm = NULL; + + rm = (struct radius_maps *)malloc(sizeof(struct radius_maps)); + if ( rm == NULL ) { + return NULL; + } + + rm->r_maps = (float **)malloc(det->n_panels*sizeof(float*)); + if ( rm->r_maps == NULL ) { + free(rm); + return NULL; + } + + rm->n_rmaps = det->n_panels; + + for( i=0 ; i<det->n_panels ; i++ ) { + + p = det->panels[i]; + rm->r_maps[i] = (float *)malloc(p.h*p.w*sizeof(float)); + + if ( rm->r_maps[i] == NULL ) { + for ( u = 0; u<i; u++ ) { + free(rm->r_maps[u]); + } + free(rm); + return NULL; + } + + for ( iss=0 ; iss<p.h ; iss++ ) { + for ( ifs=0; ifs<p.w; ifs++ ) { + + int rmi; + int x,y; + + rmi = ifs + p.w * iss; + + x = (p.cnx + ifs * p.fsx + iss * p.ssx); + y = (p.cny + ifs * p.fsy + iss * p.ssy); + + rm->r_maps[i][rmi] = sqrt(x * x + y * y); + } + } + } + return rm; +} + + +static void free_radius_maps(struct radius_maps *r_maps) +{ + int i; + + for ( i=0 ; i<r_maps->n_rmaps ; i++ ) { + free(r_maps->r_maps[i]); + } + free(r_maps->r_maps); + free(r_maps); +} + + +static struct peakfinder_mask *create_peakfinder_mask(struct image *img, + struct radius_maps *rmps, + int min_res, + int max_res) +{ + int i; + struct peakfinder_mask *msk; + + msk = (struct peakfinder_mask *)malloc(sizeof(struct peakfinder_mask)); + msk->masks =(char **) malloc(img->det->n_panels*sizeof(char*)); + msk->n_masks = img->det->n_panels; + for ( i=0; i<img->det->n_panels; i++) { + + struct panel p; + int iss, ifs; + + p = img->det->panels[i]; + + msk->masks[i] = (char *)calloc(p.w*p.h,sizeof(char)); + + for ( iss=0 ; iss<p.h ; iss++ ) { + for ( ifs=0 ; ifs<p.w ; ifs++ ) { + + int idx; + + idx = ifs + iss*p.w; + + if ( rmps->r_maps[i][idx] < max_res + && rmps->r_maps[i][idx] > min_res ) { + + if (! ( ( img->bad != NULL ) + && ( img->bad[i] != NULL ) + && ( img->bad[i][idx] != 0 ) ) ) { + msk->masks[i][idx] = 1; + } + + } + } + } + } + return msk; +} + + +static void free_peakfinder_mask(struct peakfinder_mask * pfmask) +{ + int i; + + for ( i=0 ; i<pfmask->n_masks ; i++ ) { + free(pfmask->masks[i]); + } + free(pfmask->masks); + free(pfmask); +} + + +static struct peakfinder_panel_data *allocate_panel_data(int num_panels) +{ + + struct peakfinder_panel_data *pfdata; + + pfdata = (struct peakfinder_panel_data *)malloc(sizeof(struct peakfinder_panel_data)); + if ( pfdata == NULL ) { + return NULL; + } + + pfdata->panel_h = (int *)malloc(num_panels*sizeof(int)); + if ( pfdata->panel_h == NULL ) { + free(pfdata); + return NULL; + } + + pfdata->panel_w = (int *)malloc(num_panels*sizeof(int)); + if ( pfdata->panel_w == NULL ) { + free(pfdata->panel_h); + free(pfdata); + return NULL; + } + + pfdata->panel_data = (float **)malloc(num_panels*sizeof(float*)); + if ( pfdata->panel_data == NULL ) { + free(pfdata->panel_w); + free(pfdata->panel_h); + free(pfdata); + return NULL; + } + + pfdata->num_panels = num_panels; + + return pfdata; +} + + +static void free_panel_data(struct peakfinder_panel_data *pfdata) +{ + free(pfdata->panel_data); + free(pfdata->panel_w); + free(pfdata->panel_h); + free(pfdata); +} + + +static void compute_num_radial_bins(int w, int h, float *r_map, float *max_r) +{ + int ifs, iss; + int pidx; + + for ( iss=0 ; iss<h ; iss++ ) { + for ( ifs=0 ; ifs<w ; ifs++ ) { + pidx = iss * w + ifs; + if ( r_map[pidx] > *max_r ) { + *max_r = r_map[pidx]; + } + } + } +} +// End of CrystFEL-only block 2 + + +static struct radial_stats* allocate_radial_stats(int num_rad_bins) +{ + struct radial_stats* rstats; + + rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats)); + if ( rstats == NULL ) { + return NULL; + } + + rstats->roffset = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->roffset == NULL ) { + free(rstats); + return NULL; + } + + rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->rthreshold == NULL ) { + free(rstats); + free(rstats->roffset); + return NULL; + } + + rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->lthreshold == NULL ) { + free(rstats); + free(rstats->rthreshold); + free(rstats->roffset); + return NULL; + } + + rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->rsigma == NULL ) { + free(rstats); + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats->lthreshold); + return NULL; + } + + rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int)); + if ( rstats->rcount == NULL ) { + free(rstats); + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats->lthreshold); + free(rstats->rsigma); + return NULL; + } + + rstats->n_rad_bins = num_rad_bins; + + return rstats; +} + + +static void free_radial_stats(struct radial_stats *rstats) +{ + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats->lthreshold); + free(rstats->rsigma); + free(rstats->rcount); + free(rstats); +} + + +static void fill_radial_bins(float *data, + int w, + int h, + float *r_map, + char *mask, + float *rthreshold, + float *lthreshold, + float *roffset, + float *rsigma, + int *rcount) +{ + int iss, ifs; + int pidx; + + int curr_r; + float value; + + for ( iss = 0; iss<h ; iss++ ) { + for ( ifs = 0; ifs<w ; ifs++ ) { + pidx = iss * w + ifs; + if ( mask[pidx] != 0 ) { + curr_r = (int)rint(r_map[pidx]); + value = data[pidx]; + if ( value < rthreshold[curr_r ] && value>lthreshold[curr_r]) { + roffset[curr_r] += value; + rsigma[curr_r] += (value * value); + rcount[curr_r] += 1; + } + } + } + } +} + + +static void compute_radial_stats(float *rthreshold, + float *lthreshold, + float *roffset, + float *rsigma, + int *rcount, + int num_rad_bins, + float min_snr, + float acd_threshold) +{ + int ri; + float this_offset, this_sigma; + + for ( ri=0 ; ri<num_rad_bins ; ri++ ) { + + if ( rcount[ri] == 0 ) { + roffset[ri] = 0; + rsigma[ri] = 0; + rthreshold[ri] = 1e9; + lthreshold[ri] = -1e9; + } else { + this_offset = roffset[ri] / rcount[ri]; + this_sigma = rsigma[ri] / rcount[ri] - (this_offset * this_offset); + if ( this_sigma >= 0 ) { + this_sigma = sqrt(this_sigma); + } + + roffset[ri] = this_offset; + rsigma[ri] = this_sigma; + rthreshold[ri] = roffset[ri] + min_snr*rsigma[ri]; + lthreshold[ri] = roffset[ri] - min_snr*rsigma[ri]; + + if ( rthreshold[ri] < acd_threshold ) { + rthreshold[ri] = acd_threshold; + } + } + } + +} + + +struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) +{ + struct peakfinder_peak_data *pkdata; + + pkdata = (struct peakfinder_peak_data*)malloc(sizeof(struct peakfinder_peak_data)); + if ( pkdata == NULL ) { + return NULL; + } + + pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int)); + if ( pkdata->npix == NULL ) { + free(pkdata); + free(pkdata->npix); + return NULL; + } + + pkdata->com_fs = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->com_fs == NULL ) { + free(pkdata->npix); + free(pkdata); + return NULL; + } + + pkdata->com_ss = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->com_ss == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata); + return NULL; + } + + pkdata->com_index = (int *)malloc(max_num_peaks*sizeof(int)); + if ( pkdata->com_ss == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata); + return NULL; + } + + pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->tot_i == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->com_index); + free(pkdata); + return NULL; + } + + pkdata->max_i = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->max_i == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->com_index); + free(pkdata->tot_i); + free(pkdata); + return NULL; + } + + pkdata->sigma = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->sigma == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->com_index); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata); + return NULL; + } + + pkdata->snr = (float *)malloc(max_num_peaks*sizeof(float)); + if ( pkdata->snr == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->com_index); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata->sigma); + free(pkdata); + return NULL; + } + + return pkdata; +} + + +static void free_peak_data(struct peakfinder_peak_data *pkdata) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata->com_index); + free(pkdata->tot_i); + free(pkdata->max_i); + free(pkdata->sigma); + free(pkdata->snr); + free(pkdata); +} + + +static struct peakfinder_intern_data *allocate_peakfinder_intern_data(int data_size, + int max_pix_count) +{ + + struct peakfinder_intern_data *intern_data; + + intern_data = (struct peakfinder_intern_data *)malloc(sizeof(struct peakfinder_intern_data)); + if ( intern_data == NULL ) { + return NULL; + } + + intern_data->pix_in_peak_map =(char *)calloc(data_size, sizeof(char)); + if ( intern_data->pix_in_peak_map == NULL ) { + free(intern_data); + return NULL; + } + + intern_data->infs =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->infs == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data); + return NULL; + } + + intern_data->inss =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->inss == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data); + return NULL; + } + + intern_data->peak_pixels =(int *)calloc(max_pix_count, sizeof(int)); + if ( intern_data->peak_pixels == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data->inss); + free(intern_data); + return NULL; + } + + return intern_data; +} + + +static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid) +{ + free(pfid->peak_pixels); + free(pfid->pix_in_peak_map); + free(pfid->infs); + free(pfid->inss); + free(pfid); +} + + + +static void peak_search(int p, + struct peakfinder_intern_data *pfinter, + float *copy, char *mask, float *r_map, + float *rthreshold, float *roffset, + int *num_pix_in_peak, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs, float *sum_com_fs, + float *sum_com_ss, float *sum_i, int max_pix_count) +{ + + int k, pi; + int curr_radius; + float curr_threshold; + int curr_fs; + int curr_ss; + float curr_i; + + int search_fs[9] = { 0, -1, 0, 1, -1, 1, -1, 0, 1 }; + int search_ss[9] = { 0, -1, -1, -1, 0, 0, 1, 1, 1 }; + int search_n = 9; + + // Loop through search pattern + for ( k=0; k<search_n; k++ ) { + + if ( (pfinter->infs[p] + search_fs[k]) < 0 ) continue; + if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue; + if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue; + if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue; + + // Neighbour point in big array + curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs; + curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; + + curr_radius = (int)rint(r_map[pi]); + curr_threshold = rthreshold[curr_radius]; + + // Above threshold? + if ( copy[pi] > curr_threshold + && pfinter->pix_in_peak_map[pi] == 0 + && mask[pi] != 0 ) { + + curr_i = copy[pi] - roffset[curr_radius]; + *sum_i += curr_i; + *sum_com_fs += curr_i * ((float)curr_fs); // for center of mass x + *sum_com_ss += curr_i * ((float)curr_ss); // for center of mass y + + pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + search_ss[k]; + pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + search_fs[k]; + pfinter->pix_in_peak_map[pi] = 1; + if ( *num_pix_in_peak < max_pix_count ) { + pfinter->peak_pixels[*num_pix_in_peak] = pi; + } + *num_pix_in_peak = *num_pix_in_peak + 1; + } + } +} + + +static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, + float *copy, float *r_map, + float *rthreshold, float *roffset, + char *pix_in_peak_map, char *mask, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs,float *local_sigma, float *local_offset, + float *background_max_i, int com_idx, + int local_bg_radius) +{ + int ssj, fsi; + float pix_radius; + int curr_fs, curr_ss; + int pi; + int curr_radius; + float curr_threshold; + float curr_i; + + int np_sigma; + int np_counted; + int local_radius; + + float sum_i; + float sum_i_squared; + + ring_width = 2 * local_bg_radius; + + sum_i = 0; + sum_i_squared = 0; + np_sigma = 0; + np_counted = 0; + local_radius = 0; + + for ( ssj = -ring_width ; ssj<ring_width ; ssj++ ) { + for ( fsi = -ring_width ; fsi<ring_width ; fsi++ ) { + + // Within-ASIC check + if ( (com_fs_int + fsi) < 0 ) continue; + if ( (com_fs_int + fsi) >= asic_size_fs ) continue; + if ( (com_ss_int + ssj) < 0 ) continue; + if ( (com_ss_int + ssj) >= asic_size_ss ) + continue; + + // Within outer ring check + pix_radius = sqrt(fsi * fsi + ssj * ssj); + if ( pix_radius>ring_width ) continue; + + // Position of this point in data stream + curr_fs = com_fs_int + fsi + aifs * asic_size_fs; + curr_ss = com_ss_int + ssj + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; + + curr_radius = (int)rint(r_map[pi]); + curr_threshold = rthreshold[curr_radius]; + + // Intensity above background ??? just intensity? + curr_i = copy[pi]; + + // Keep track of value and value-squared for offset and sigma calculation + if ( curr_i < curr_threshold && pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) { + + np_sigma++; + sum_i += curr_i; + sum_i_squared += (curr_i * curr_i); + + if ( curr_i > *background_max_i ) { + *background_max_i = curr_i; + } + } + np_counted += 1; + } + } + + // Calculate local background and standard deviation + if ( np_sigma != 0 ) { + *local_offset = sum_i / np_sigma; + *local_sigma = sum_i_squared / np_sigma - (*local_offset * *local_offset); + if (*local_sigma >= 0) { + *local_sigma = sqrt(*local_sigma); + } else { + *local_sigma = 0.01; + } + } else { + local_radius = (int)rint(r_map[(int)rint(com_idx)]); + *local_offset = roffset[local_radius]; + *local_sigma = 0.01; + } +} + + +static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, + int aiss, int aifs, float *rthreshold, + float *roffset, int *peak_count, + float *copy, struct peakfinder_intern_data *pfinter, + float *r_map, char *mask, int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr, int max_n_peaks) +{ + int pxss, pxfs; + int num_pix_in_peak; + + // Loop over pixels within a module + for ( pxss=1 ; pxss<asic_size_ss-1 ; pxss++ ) { + for ( pxfs=1 ; pxfs<asic_size_fs-1 ; pxfs++ ) { + + float curr_thresh; + int pxidx; + int curr_rad; + + pxidx = (pxss + aiss * asic_size_ss) * num_pix_fs + + pxfs + aifs * asic_size_fs; + + curr_rad = (int)rint(r_map[pxidx]); + curr_thresh = rthreshold[curr_rad]; + + if ( copy[pxidx] > curr_thresh + && pfinter->pix_in_peak_map[pxidx] == 0 + && mask[pxidx] != 0 ) { //??? not sure if needed + + // This might be the start of a new peak - start searching + float sum_com_fs, sum_com_ss; + float sum_i; + float peak_com_fs, peak_com_ss; + float peak_com_fs_int, peak_com_ss_int; + float peak_tot_i, pk_tot_i_raw; + float peak_max_i, pk_max_i_raw; + float peak_snr; + float local_sigma, local_offset; + float background_max_i; + int lt_num_pix_in_pk; + int ring_width; + int peak_idx; + int com_idx; + int p; + + pfinter->infs[0] = pxfs; + pfinter->inss[0] = pxss; + pfinter->peak_pixels[0] = pxidx; + num_pix_in_peak = 0; //y 1; + + sum_i = 0; + sum_com_fs = 0; + sum_com_ss = 0; + + // Keep looping until the pixel count within this peak does not change + do { + lt_num_pix_in_pk = num_pix_in_peak; + + // Loop through points known to be within this peak + for ( p=0; p<=num_pix_in_peak; p++ ) { //changed from 1 to 0 by O.Y. + peak_search(p, + pfinter, copy, mask, + r_map, + rthreshold, + roffset, + &num_pix_in_peak, + asic_size_fs, + asic_size_ss, + aifs, aiss, + num_pix_fs, + &sum_com_fs, + &sum_com_ss, + &sum_i, + max_pix_count); + } + + } while ( lt_num_pix_in_pk != num_pix_in_peak ); + + // Too many or too few pixels means ignore this 'peak'; move on now + if ( num_pix_in_peak < min_pix_count || num_pix_in_peak > max_pix_count ) continue; + + // If for some reason sum_i is 0 - it's better to skip + if ( fabs(sum_i) < 1e-10 ) continue; + + // Calculate center of mass for this peak from initial peak search + peak_com_fs = sum_com_fs / fabs(sum_i); + peak_com_ss = sum_com_ss / fabs(sum_i); + + com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * num_pix_fs; + + peak_com_fs_int = (int)rint(peak_com_fs) - aifs * asic_size_fs; + peak_com_ss_int = (int)rint(peak_com_ss) - aiss * asic_size_ss; + + // Calculate the local signal-to-noise ratio and local background in an annulus around + // this peak (excluding pixels which look like they might be part of another peak) + local_sigma = 0.0; + local_offset = 0.0; + background_max_i = 0.0; + + ring_width = 2 * local_bg_radius; + + search_in_ring(ring_width, peak_com_fs_int, + peak_com_ss_int, + copy, r_map, rthreshold, + roffset, + pfinter->pix_in_peak_map, + mask, asic_size_fs, + asic_size_ss, + aifs, aiss, + num_pix_fs, + &local_sigma, + &local_offset, + &background_max_i, + com_idx, local_bg_radius); + + // Re-integrate (and re-centroid) peak using local background estimates + peak_tot_i = 0; + pk_tot_i_raw = 0; + peak_max_i = 0; + pk_max_i_raw = 0; + sum_com_fs = 0; + sum_com_ss = 0; + + for ( peak_idx = 0 ; + peak_idx < num_pix_in_peak && peak_idx < max_pix_count ; + peak_idx++ ) { + + int curr_idx; + float curr_i; + float curr_i_raw; + int curr_fs, curr_ss; + + curr_idx = pfinter->peak_pixels[peak_idx]; + curr_i_raw = copy[curr_idx]; + curr_i = curr_i_raw - local_offset; + peak_tot_i += curr_i; + pk_tot_i_raw += curr_i_raw; + + // Remember that curr_idx = curr_fs + curr_ss*num_pix_fs + curr_fs = curr_idx % num_pix_fs; + curr_ss = curr_idx / num_pix_fs; + sum_com_fs += curr_i * ((float)curr_fs); + sum_com_ss += curr_i * ((float)curr_ss); + + if ( curr_i_raw > pk_max_i_raw ) pk_max_i_raw = curr_i_raw; + if ( curr_i > peak_max_i ) peak_max_i = curr_i; + } + + + // This CAN happen! Better to skip... + if ( fabs(peak_tot_i) < 1e-10 ) continue; + + peak_com_fs = sum_com_fs / fabs(peak_tot_i); + peak_com_ss = sum_com_ss / fabs(peak_tot_i); + + // Calculate signal-to-noise and apply SNR criteria + if ( fabs(local_sigma) > 1e-10 ) { + peak_snr = peak_tot_i / local_sigma; + } else { + peak_snr = 0; + } + + if (peak_snr < min_snr) continue; + + // Is the maximum intensity in the peak enough above intensity in background region to + // be a peak and not noise? The more pixels there are in the peak, the more relaxed we + // are about this criterion + //f_background_thresh = background_max_i - local_offset; //!!! Ofiget'! If I uncomment + // if (peak_max_i < f_background_thresh) { // these lines the result is + // different! + if (peak_max_i < background_max_i - local_offset) continue; + + // This is a peak? If so, add info to peak list + if ( num_pix_in_peak >= min_pix_count + && num_pix_in_peak <= max_pix_count ) { + + // Bragg peaks in the mask + for ( peak_idx = 0 ; + peak_idx < num_pix_in_peak && + peak_idx < max_pix_count ; + peak_idx++ ) { + pfinter->pix_in_peak_map[pfinter->peak_pixels[peak_idx]] = 2; + } + + int peak_com_idx; + peak_com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * + num_pix_fs; + // Remember peak information + if ( *peak_count < max_n_peaks ) { + + int pidx; + pidx = *peak_count; + + npix[pidx] = num_pix_in_peak; + com_fs[pidx] = peak_com_fs; + com_ss[pidx] = peak_com_ss; + com_index[pidx] = peak_com_idx; + tot_i[pidx] = peak_tot_i; + max_i[pidx] = peak_max_i; + sigma[pidx] = local_sigma; + snr[pidx] = peak_snr; + } + *peak_count += 1; + } + } + } + } +} + + +static int peakfinder8_base(float *roffset, float *rthreshold, + float *data, char *mask, float *r_map, + int asic_size_fs, int num_asics_fs, + int asic_size_ss, int num_asics_ss, + int max_n_peaks, int *num_found_peaks, + int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr, + char* outliersMask) +{ + + int num_pix_fs, num_pix_ss, num_pix_tot; + int aifs, aiss; + int peak_count; + struct peakfinder_intern_data *pfinter; + + num_pix_fs = asic_size_fs * num_asics_fs; + num_pix_ss = asic_size_ss * num_asics_ss; + num_pix_tot = num_pix_fs * num_pix_ss; + + pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count); + if ( pfinter == NULL ) { + return 1; + } + + peak_count = 0; + + // Loop over modules (nxn array) + for ( aiss=0 ; aiss<num_asics_ss ; aiss++ ) { + for ( aifs=0 ; aifs<num_asics_fs ; aifs++ ) { // ??? to change to proper panels need + process_panel(asic_size_fs, asic_size_ss, num_pix_fs, // change copy, mask, r_map + aiss, aifs, rthreshold, roffset, + &peak_count, data, pfinter, r_map, mask, + npix, com_fs, com_ss, com_index, tot_i, + max_i, sigma, snr, min_pix_count, + max_pix_count, local_bg_radius, min_snr, + max_n_peaks); + } + } + *num_found_peaks = peak_count; + + if (outliersMask != NULL) { + memcpy(outliersMask, pfinter->pix_in_peak_map, num_pix_tot*sizeof(char)); + } + + free_peakfinder_intern_data(pfinter); + + return 0; +} + + +int peakfinder8(struct image *img, int max_n_peaks, + float threshold, float min_snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, int min_res, + int max_res, int use_saturated) +{ + struct radius_maps *rmaps; + struct peakfinder_mask *pfmask; + struct peakfinder_panel_data *pfdata; + struct radial_stats *rstats; + struct peakfinder_peak_data *pkdata; + int num_rad_bins; + int pi; + int i, it_counter; + int num_found_peaks; + int remaining_max_num_peaks; + int iterations; + float max_r; + + iterations = 5; + + if ( img-> det == NULL) { + return 1; + } + + rmaps = compute_radius_maps(img->det); + if ( rmaps == NULL ) { + return 1; + } + + pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res); + if ( pfmask == NULL ) { + free_radius_maps(rmaps); + return 1; + } + + pfdata = allocate_panel_data(img->det->n_panels); + if ( pfdata == NULL) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + return 1; + } + + for ( pi=0 ; pi<img->det->n_panels ; pi++ ) { + pfdata->panel_h[pi] = img->det->panels[pi].h; + pfdata->panel_w[pi] = img->det->panels[pi].w; + pfdata->panel_data[pi] = img->dp[pi]; + pfdata->num_panels = img->det->n_panels; + } + + max_r = -1e9; + + for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) { + + compute_num_radial_bins(pfdata->panel_w[pi], + pfdata->panel_h[pi], + rmaps->r_maps[pi], + &max_r); + } + + num_rad_bins = (int)ceil(max_r) + 1; + + rstats = allocate_radial_stats(num_rad_bins); + if ( rstats == NULL ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + return 1; + } + + for ( i=0 ; i<rstats->n_rad_bins ; i++) { + rstats->rthreshold[i] = 1e9; + } + + for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) { + + for ( i=0; i<num_rad_bins; i++ ) { + rstats->roffset[i] = 0; + rstats->rsigma[i] = 0; + rstats->rcount[i] = 0; + } + + for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) { + + fill_radial_bins(pfdata->panel_data[pi], + pfdata->panel_w[pi], + pfdata->panel_h[pi], + rmaps->r_maps[pi], + pfmask->masks[pi], + rstats->rthreshold, + rstats->lthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount); + + } + + compute_radial_stats(rstats->rthreshold, + rstats->lthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount, + num_rad_bins, + min_snr, + threshold); + + } + + pkdata = allocate_peak_data(max_n_peaks); + if ( pkdata == NULL ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } + + remaining_max_num_peaks = max_n_peaks; + + for ( pi=0 ; pi<img->det->n_panels ; pi++) { + + int peaks_to_add; + int pki; + int ret; + + num_found_peaks = 0; + + if ( img->det->panels[pi].no_index ) { + continue; + } + + ret = peakfinder8_base(rstats->roffset, + rstats->rthreshold, + pfdata->panel_data[pi], + pfmask->masks[pi], + rmaps->r_maps[pi], + pfdata->panel_w[pi], 1, + pfdata->panel_h[pi], 1, + max_n_peaks, + &num_found_peaks, + pkdata->npix, + pkdata->com_fs, + pkdata->com_ss, + pkdata->com_index, + pkdata->tot_i, + pkdata->max_i, + pkdata->sigma, + pkdata->snr, + min_pix_count, + max_pix_count, + local_bg_radius, + min_snr, + NULL); + + if ( ret != 0 ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } + + peaks_to_add = num_found_peaks; + + if ( num_found_peaks > remaining_max_num_peaks ) { + peaks_to_add = remaining_max_num_peaks; + } + + remaining_max_num_peaks -= peaks_to_add; + + for ( pki=0 ; pki<peaks_to_add ; pki++ ) { + + struct panel *p; + + p = &img->det->panels[pi]; + + img->num_peaks += 1; + if ( pkdata->max_i[pki] > p->max_adu ) { + img->num_saturated_peaks++; + if ( !use_saturated ) { + continue; + } + } + + image_add_feature(img->features, + pkdata->com_fs[pki]+0.5, + pkdata->com_ss[pki]+0.5, + p, + img, + pkdata->tot_i[pki], + NULL); + } + } + + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + free_peak_data(pkdata); + return 0; +} diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h new file mode 100644 index 00000000..483eebdf --- /dev/null +++ b/libcrystfel/src/peakfinder8.h @@ -0,0 +1,50 @@ +/* + * peakfinder8.h + * + * The peakfinder8 algorithm + * + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017 Valerio Mariani <valerio.mariani@desy.de> + * 2017 Anton Barty <anton.barty@desy.de> + * 2017 Oleksandr Yefanov <oleksandr.yefanov@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef PEAKFINDER8_H +#define PEAKFINDER8_H + +#include "image.h" + +#ifdef __cplusplus +extern "C" { +#endif + +extern int peakfinder8(struct image *img, int max_n_peaks, + float threshold, float min_snr, + int mix_pix_count, int max_pix_count, + int local_bg_radius, int min_res, + int max_res, int use_saturated); + +#ifdef __cplusplus +} +#endif + +#endif // PEAKFINDER8_H diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index cb2a9f5a..28c79538 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -3,7 +3,7 @@ * * Peak search and other image analysis * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * @@ -12,6 +12,7 @@ * 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de> * 2011 Andrew Martin <andrew.martin@desy.de> * 2011 Richard Kirian + * 2017 Valerio Mariani <valerio.mariani@desy.de> * * This file is part of CrystFEL. * @@ -52,6 +53,7 @@ #include "reflist-utils.h" #include "cell-utils.h" #include "geometry.h" +#include "peakfinder8.h" static int cull_peaks_in_panel(struct image *image, struct panel *p) @@ -517,8 +519,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, void search_peaks(struct image *image, float threshold, float min_gradient, - float min_snr, double ir_inn, double ir_mid, double ir_out, - int use_saturated) + float min_snr, double ir_inn, double ir_mid, + double ir_out, int use_saturated) { int i; @@ -541,6 +543,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } +int search_peaks_peakfinder8(struct image *image, int max_n_peaks, + float threshold, float min_snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, int min_res, + int max_res, int use_saturated) +{ + if ( image->features != NULL ) { + image_feature_list_free(image->features); + } + image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; + + return peakfinder8(image, max_n_peaks, threshold, min_snr, + min_pix_count, max_pix_count, + local_bg_radius, min_res, + max_res, use_saturated); +} + + int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { int n_feat = 0; diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index ba724419..a5095127 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -3,11 +3,12 @@ * * Peak search and other image analysis * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2010-2015 Thomas White <taw@physics.org> + * 2017 Valerio Mariani <valerio.mariani@desy.de> * * This file is part of CrystFEL. * @@ -45,9 +46,14 @@ extern "C" { extern int *make_BgMask(struct image *image, struct panel *p, double ir_inn); extern void search_peaks(struct image *image, float threshold, - float min_gradient, float min_snr, - double ir_inn, double ir_mid, double ir_out, - int use_saturated); + float min_gradient, float min_snr, double ir_inn, + double ir_mid, double ir_out, int use_saturated); + +extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks, + float threshold, float min_snr, + int mix_pix_count, int max_pix_count, + int local_bg_radius, int min_res, + int max_res, int use_saturated); extern int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst); diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index 4621f4f4..380a1b8a 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -27,8 +27,6 @@ * */ -#define _ISOC99_SOURCE -#define _GNU_SOURCE #include <math.h> #include <stdio.h> #include <assert.h> diff --git a/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c index 56273fdb..ccf35194 100644 --- a/libcrystfel/src/statistics.c +++ b/libcrystfel/src/statistics.c @@ -30,7 +30,6 @@ #include <config.h> #endif -#define _ISOC99_SOURCE #include <math.h> #include <stdlib.h> #include <gsl/gsl_errno.h> diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index e858172e..fb4b0c70 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -806,8 +806,8 @@ static int write_crystal(Stream *st, Crystal *cr, int include_reflections) } -int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, - int include_peaks, int include_reflections, struct event* ev) +int write_chunk(Stream *st, struct image *i, struct imagefile *imfile, + int include_peaks, int include_reflections, struct event *ev) { int j; char *indexer; @@ -832,7 +832,7 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, fprintf(st->fh, "beam_divergence = %.2e rad\n", i->div); fprintf(st->fh, "beam_bandwidth = %.2e (fraction)\n", i->bw); - copy_hdf5_fields(hdfile, i->copyme, st->fh, ev); + imagefile_copy_fields(imfile, i->copyme, st->fh, ev); if ( i->det != NULL ) { @@ -1196,13 +1196,9 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf) } if ( strncmp(line, "indexed_by = ", 13) == 0 ) { - IndexingMethod *list; - list = build_indexer_list(line+13); - if ( list == NULL ) { + image->indexed_by = get_indm_from_string(line+13); + if ( image->indexed_by == INDEXING_ERROR ) { ERROR("Failed to read indexer list\n"); - } else { - image->indexed_by = list[0]; - free(list); } } diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index ff8628c0..764e3e36 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -3,11 +3,11 @@ * * Stream tools * - * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2017 Thomas White <taw@physics.org> * 2014 Valerio Mariani * 2011 Andrew Aquila * @@ -39,6 +39,7 @@ struct image; struct hdfile; struct event; +struct imagefile; #include "cell.h" #define GEOM_START_MARKER "----- Begin geometry file -----" @@ -106,10 +107,15 @@ extern int read_chunk(Stream *st, struct image *image); extern int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf); -extern int write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, +extern int write_chunk(Stream *st, struct image *image, struct imagefile *imfile, int include_peaks, int include_reflections, struct event *ev); +extern int write_chunk_2(Stream *st, struct image *image, + struct imagefile *imfile, + int include_peaks, int include_reflections, + struct event *ev); + extern void write_command(Stream *st, int argc, char *argv[]); extern void write_geometry_file(Stream *st, const char *geom_filename); extern int rewind_stream(Stream *st); diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c index 461da0cd..90d81a71 100644 --- a/libcrystfel/src/utils.c +++ b/libcrystfel/src/utils.c @@ -387,144 +387,6 @@ int poisson_noise(gsl_rng *rng, double expected) } -/** - * SECTION:quaternion - * @short_description: Simple quaternion handling - * @title: Quaternion - * @section_id: - * @see_also: - * @include: "utils.h" - * @Image: - * - * There is a simple quaternion structure in CrystFEL. At the moment, it is - * only used when simulating patterns, as an argument to cell_rotate() to - * orient the unit cell. - */ - -/** - * quaternion_modulus: - * @q: A %quaternion - * - * If a quaternion represents a pure rotation, its modulus should be unity. - * - * Returns: the modulus of the given quaternion. - **/ -double quaternion_modulus(struct quaternion q) -{ - return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); -} - - -/** - * normalise_quaternion: - * @q: A %quaternion - * - * Rescales the quaternion such that its modulus is unity. - * - * Returns: the normalised version of @q - **/ -struct quaternion normalise_quaternion(struct quaternion q) -{ - double mod; - struct quaternion r; - - mod = quaternion_modulus(q); - - r.w = q.w / mod; - r.x = q.x / mod; - r.y = q.y / mod; - r.z = q.z / mod; - - return r; -} - - -/** - * random_quaternion: - * @rng: A GSL random number generator to use - * - * Returns: a randomly generated, normalised, quaternion. - **/ -struct quaternion random_quaternion(gsl_rng *rng) -{ - struct quaternion q; - - q.w = 2.0*gsl_rng_uniform(rng) - 1.0; - q.x = 2.0*gsl_rng_uniform(rng) - 1.0; - q.y = 2.0*gsl_rng_uniform(rng) - 1.0; - q.z = 2.0*gsl_rng_uniform(rng) - 1.0; - q = normalise_quaternion(q); - - return q; -} - - -/** - * quaternion_valid: - * @q: A %quaternion - * - * Checks if the given quaternion is normalised. - * - * This function performs a nasty floating point comparison of the form - * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be - * relied upon to spot anything other than the most obvious input error. - * - * Returns: 1 if the quaternion is normalised, 0 if not. - **/ -int quaternion_valid(struct quaternion q) -{ - double qmod; - - qmod = quaternion_modulus(q); - - /* Modulus = 1 to within some tolerance? - * Nasty allowance for floating-point accuracy follows... */ - if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; - - return 0; -} - - -/** - * quat_rot - * @q: A vector (in the form of a "struct rvec") - * @z: A %quaternion - * - * Rotates a vector according to a quaternion. - * - * Returns: A rotated version of @p. - **/ -struct rvec quat_rot(struct rvec q, struct quaternion z) -{ - struct rvec res; - double t01, t02, t03, t11, t12, t13, t22, t23, t33; - - t01 = z.w*z.x; - t02 = z.w*z.y; - t03 = z.w*z.z; - t11 = z.x*z.x; - t12 = z.x*z.y; - t13 = z.x*z.z; - t22 = z.y*z.y; - t23 = z.y*z.z; - t33 = z.z*z.z; - - res.u = (1.0 - 2.0 * (t22 + t33)) * q.u - + (2.0 * (t12 + t03)) * q.v - + (2.0 * (t13 - t02)) * q.w; - - res.v = (2.0 * (t12 - t03)) * q.u - + (1.0 - 2.0 * (t11 + t33)) * q.v - + (2.0 * (t01 + t23)) * q.w; - - res.w = (2.0 * (t02 + t13)) * q.u - + (2.0 * (t23 - t01)) * q.v - + (1.0 - 2.0 * (t11 + t22)) * q.w; - - return res; -} - - /* Return non-zero if c is in delims */ static int assplode_isdelim(const char c, const char *delims) { @@ -691,3 +553,141 @@ void utils_fudge_gslcblas() { STATUS("%p\n", cblas_sgemm); } + + +/** + * SECTION:quaternion + * @short_description: Simple quaternion handling + * @title: Quaternion + * @section_id: + * @see_also: + * @include: "utils.h" + * @Image: + * + * There is a simple quaternion structure in CrystFEL. At the moment, it is + * only used when simulating patterns, as an argument to cell_rotate() to + * orient the unit cell. + */ + +/** + * quaternion_modulus: + * @q: A %quaternion + * + * If a quaternion represents a pure rotation, its modulus should be unity. + * + * Returns: the modulus of the given quaternion. + **/ +double quaternion_modulus(struct quaternion q) +{ + return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); +} + + +/** + * normalise_quaternion: + * @q: A %quaternion + * + * Rescales the quaternion such that its modulus is unity. + * + * Returns: the normalised version of @q + **/ +struct quaternion normalise_quaternion(struct quaternion q) +{ + double mod; + struct quaternion r; + + mod = quaternion_modulus(q); + + r.w = q.w / mod; + r.x = q.x / mod; + r.y = q.y / mod; + r.z = q.z / mod; + + return r; +} + + +/** + * random_quaternion: + * @rng: A GSL random number generator to use + * + * Returns: a randomly generated, normalised, quaternion. + **/ +struct quaternion random_quaternion(gsl_rng *rng) +{ + struct quaternion q; + + q.w = 2.0*gsl_rng_uniform(rng) - 1.0; + q.x = 2.0*gsl_rng_uniform(rng) - 1.0; + q.y = 2.0*gsl_rng_uniform(rng) - 1.0; + q.z = 2.0*gsl_rng_uniform(rng) - 1.0; + q = normalise_quaternion(q); + + return q; +} + + +/** + * quaternion_valid: + * @q: A %quaternion + * + * Checks if the given quaternion is normalised. + * + * This function performs a nasty floating point comparison of the form + * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be + * relied upon to spot anything other than the most obvious input error. + * + * Returns: 1 if the quaternion is normalised, 0 if not. + **/ +int quaternion_valid(struct quaternion q) +{ + double qmod; + + qmod = quaternion_modulus(q); + + /* Modulus = 1 to within some tolerance? + * Nasty allowance for floating-point accuracy follows... */ + if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; + + return 0; +} + + +/** + * quat_rot + * @q: A vector (in the form of a "struct rvec") + * @z: A %quaternion + * + * Rotates a vector according to a quaternion. + * + * Returns: A rotated version of @p. + **/ +struct rvec quat_rot(struct rvec q, struct quaternion z) +{ + struct rvec res; + double t01, t02, t03, t11, t12, t13, t22, t23, t33; + + t01 = z.w*z.x; + t02 = z.w*z.y; + t03 = z.w*z.z; + t11 = z.x*z.x; + t12 = z.x*z.y; + t13 = z.x*z.z; + t22 = z.y*z.y; + t23 = z.y*z.z; + t33 = z.z*z.z; + + res.u = (1.0 - 2.0 * (t22 + t33)) * q.u + + (2.0 * (t12 + t03)) * q.v + + (2.0 * (t13 - t02)) * q.w; + + res.v = (2.0 * (t12 - t03)) * q.u + + (1.0 - 2.0 * (t11 + t33)) * q.v + + (2.0 * (t01 + t23)) * q.w; + + res.w = (2.0 * (t02 + t13)) * q.u + + (2.0 * (t23 - t01)) * q.v + + (1.0 - 2.0 * (t11 + t22)) * q.w; + + return res; +} diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 4955f875..a759ff15 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -61,44 +61,10 @@ #define THOMSON_LENGTH (2.81794e-15) -/* ------------------------------ Quaternions ------------------------------- */ - -/** - * quaternion: - * - * <programlisting> - * struct quaternion - * { - * double w - * double x - * double y - * double z - * }; - * </programlisting> - * - * A structure representing a quaternion. - * - **/ -struct quaternion; - -struct quaternion { - double w; - double x; - double y; - double z; -}; - #ifdef __cplusplus extern "C" { #endif -extern struct quaternion normalise_quaternion(struct quaternion q); -extern double quaternion_modulus(struct quaternion q); -extern struct quaternion random_quaternion(gsl_rng *rng); -extern int quaternion_valid(struct quaternion q); -extern struct rvec quat_rot(struct rvec q, struct quaternion z); - - /* --------------------------- Useful functions ----------------------------- */ extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v); @@ -126,17 +92,6 @@ extern double flat_noise(gsl_rng *rng, double expected, double width); extern double gaussian_noise(gsl_rng *rng, double expected, double stddev); extern int poisson_noise(gsl_rng *rng, double expected); -/* Keep these ones inline, to avoid function call overhead */ -static inline struct quaternion invalid_quaternion(void) -{ - struct quaternion quat; - quat.w = 0.0; - quat.x = 0.0; - quat.y = 0.0; - quat.z = 0.0; - return quat; -} - static inline double distance(double x1, double y1, double x2, double y2) { return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); @@ -275,6 +230,47 @@ extern char *safe_basename(const char *in); #define unlikely(x) (x) #endif + +/* ------------------------------ Quaternions ------------------------------- */ + +/** + * quaternion: + * @w: component + * @x: component + * @y: component + * @z: component + * + * A structure representing a quaternion. + * + **/ +struct quaternion; + +struct quaternion { + double w; + double x; + double y; + double z; +}; + +extern struct quaternion normalise_quaternion(struct quaternion q); +extern double quaternion_modulus(struct quaternion q); +extern struct quaternion random_quaternion(gsl_rng *rng); +extern int quaternion_valid(struct quaternion q); +extern struct rvec quat_rot(struct rvec q, struct quaternion z); + +/* Keep these ones inline, to avoid function call overhead */ +static inline struct quaternion invalid_quaternion(void) +{ + struct quaternion quat; + quat.w = 0.0; + quat.x = 0.0; + quat.y = 0.0; + quat.z = 0.0; + return quat; +} + + + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 15826760..8a7daab7 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -3,12 +3,12 @@ * * Invoke xds for crystal autoindexing * - * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2013 Cornelius Gati * * Authors: - * 2010-2016 Thomas White <taw@physics.org> + * 2010-2017 Thomas White <taw@physics.org> * 2013 Cornelius Gati <cornelius.gati@cfel.de> * * This file is part of CrystFEL. @@ -498,7 +498,7 @@ static int write_inp(struct image *image, struct xds_private *xp) } -int run_xds(struct image *image, IndexingPrivate *priv) +int run_xds(struct image *image, void *priv) { unsigned int opts; int status; @@ -619,8 +619,8 @@ int run_xds(struct image *image, IndexingPrivate *priv) } -IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl) +void *xds_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl) { struct xds_private *xp; int need_cell = 0; @@ -687,11 +687,11 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, xp->cell = cell; xp->indm = *indm; - return (IndexingPrivate *)xp; + return xp; } -void xds_cleanup(IndexingPrivate *pp) +void xds_cleanup(void *pp) { struct xds_private *xp; diff --git a/libcrystfel/src/xds.h b/libcrystfel/src/xds.h index a0db2054..094d6d2f 100644 --- a/libcrystfel/src/xds.h +++ b/libcrystfel/src/xds.h @@ -3,13 +3,13 @@ * * Invoke xds for crystal autoindexing * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2013 Cornelius Gati * * Authors: - * 2010-2013 Thomas White <taw@physics.org> - * 2013 Cornelius Gati <cornelius.gati@cfel.de> + * 2010-2013,2017 Thomas White <taw@physics.org> + * 2013 Cornelius Gati <cornelius.gati@cfel.de> * * This file is part of CrystFEL. * @@ -42,12 +42,12 @@ extern "C" { #endif -extern int run_xds(struct image *image, IndexingPrivate *ipriv); +extern int run_xds(struct image *image, void *ipriv); -extern IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell, - struct detector *det, float *ltl); +extern void *xds_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl); -extern void xds_cleanup(IndexingPrivate *pp); +extern void xds_cleanup(void *pp); #ifdef __cplusplus } |