diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/CMakeLists.txt | 14 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.c | 5 | ||||
-rw-r--r-- | libcrystfel/src/crystal.h | 3 | ||||
-rw-r--r-- | libcrystfel/src/detector.c | 13 | ||||
-rw-r--r-- | libcrystfel/src/hdf5-file.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 309 | ||||
-rw-r--r-- | libcrystfel/src/image.h | 30 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 41 | ||||
-rw-r--r-- | libcrystfel/src/index.h | 6 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 57 | ||||
-rw-r--r-- | libcrystfel/src/peakfinder8.c | 10 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 8 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/reflist-utils.c | 28 | ||||
-rw-r--r-- | libcrystfel/src/reflist-utils.h | 6 | ||||
-rw-r--r-- | libcrystfel/src/reflist.c | 36 | ||||
-rw-r--r-- | libcrystfel/src/reflist.h | 19 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 5 | ||||
-rw-r--r-- | libcrystfel/src/symmetry.c | 49 | ||||
-rw-r--r-- | libcrystfel/src/utils.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/xds.c | 154 | ||||
-rw-r--r-- | libcrystfel/src/xgandalf.c | 335 | ||||
-rw-r--r-- | libcrystfel/src/xgandalf.h | 58 |
23 files changed, 915 insertions, 279 deletions
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index f6f25ec9..3d7bec38 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -6,6 +6,7 @@ find_package(XGANDALF) find_package(PINKINDEXER) find_package(NBP) find_package(FDIP) +find_package(ZLIB REQUIRED) pkg_search_module(FFTW fftw3) set(HAVE_CURSES ${CURSES_FOUND}) @@ -68,6 +69,7 @@ set(LIBCRYSTFEL_SOURCES src/felix.c src/peakfinder8.c src/taketwo.c + src/xgandalf.c ) if (HAVE_FFTW) @@ -106,6 +108,7 @@ set(LIBCRYSTFEL_HEADERS src/felix.h src/peakfinder8.h src/taketwo.h + src/xgandalf.h ) add_library(${PROJECT_NAME} SHARED @@ -123,8 +126,9 @@ target_include_directories(${PROJECT_NAME} INTERFACE ${PROJECT_SOURCE_DIR}/src) include_directories(${CMAKE_CURRENT_BINARY_DIR}) add_definitions(-DHAVE_CONFIG_H) -target_include_directories(${PROJECT_NAME} PRIVATE ${HDF5_INCLUDE_DIRS}) -target_link_libraries(${PROJECT_NAME} PRIVATE util ${HDF5_C_LIBRARIES} Threads::Threads GSL::gsl m) +target_include_directories(${PROJECT_NAME} PRIVATE ${HDF5_INCLUDE_DIRS} ${ZLIB_INCLUDE_DIRS}) +target_link_libraries(${PROJECT_NAME} PRIVATE util ${HDF5_C_LIBRARIES} ${ZLIB_LIBRARIES} + Threads::Threads GSL::gsl m) if (XGANDALF_FOUND) target_include_directories(${PROJECT_NAME} PRIVATE ${XGANDALF_INCLUDES}) @@ -162,9 +166,11 @@ if (CURSES_FOUND) endif (CURSES_FOUND) target_compile_options(${PROJECT_NAME} PRIVATE -Wall) +set_target_properties(${PROJECT_NAME} PROPERTIES PUBLIC_HEADER "${LIBCRYSTFEL_HEADERS}") -install (TARGETS libcrystfel - LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) +install(TARGETS libcrystfel + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/crystfel + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) # crystfel.pc configure_file(crystfel.pc.in crystfel.pc) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 5465052c..7b1984bb 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -208,7 +208,10 @@ int right_handed(UnitCell *cell) rh_direct = aCb_dot_c > 0.0; - assert(rh_reciprocal == rh_direct); + if ( rh_reciprocal != rh_direct ) { + ERROR("Whoops, reciprocal and real space handedness are " + "not the same!\n"); + } return rh_direct; } diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h index f4c7477b..cc5754ca 100644 --- a/libcrystfel/src/crystal.h +++ b/libcrystfel/src/crystal.h @@ -36,7 +36,6 @@ #include "cell.h" -#include "reflist.h" /** @@ -47,6 +46,8 @@ **/ typedef struct _crystal Crystal; +#include "reflist.h" + #ifdef __cplusplus extern "C" { #endif diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c index 8d33731a..d4d8b768 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -1429,6 +1429,8 @@ struct detector *get_detector_geometry_2(const char *filename, if ( det->panels[i].dim_structure->dims[di] == HYSL_UNDEFINED ) { dim_dim_reject = 1; + ERROR("Dimension %i for panel %s is undefined.\n", + di, det->panels[i].name); } if ( det->panels[i].dim_structure->dims[di] == HYSL_PLACEHOLDER ) { @@ -1446,18 +1448,23 @@ struct detector *get_detector_geometry_2(const char *filename, } if ( found_ss != 1 ) { - ERROR("Only one slow scan dim coordinate is allowed\n"); + ERROR("Exactly one slow scan dim coordinate is needed " + "(found %i for panel %s)\n", found_ss, + det->panels[i].name); dim_dim_reject = 1; } if ( found_fs != 1 ) { - ERROR("Only one fast scan dim coordinate is allowed\n"); + ERROR("Exactly one fast scan dim coordinate is needed " + "(found %i for panel %s)\n", found_fs, + det->panels[i].name); dim_dim_reject = 1; } if ( panel_dim_dim > 1 ) { ERROR("Maximum one placeholder dim coordinate is " - "allowed\n"); + "allowed (found %i for panel %s)\n", + panel_dim_dim, det->panels[i].name); dim_dim_reject = 1; } diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index 7e3d6b11..f9c2aa97 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -1047,9 +1047,9 @@ int hdf5_write_image(const char *filename, const struct image *image, write_photon_energy(fh, ph_lambda_to_eV(image->lambda), ph_en_loc); - if ( image->spectrum != NULL && image->spectrum_size > 0 ) { + if ( image->spectrum0 != NULL && image->spectrum_size > 0 ) { - write_spectrum(fh, image->spectrum, image->spectrum_size, + write_spectrum(fh, image->spectrum0, image->spectrum_size, image->nsamples); } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 0dda66ed..26d48007 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -34,6 +34,7 @@ #include <math.h> #include <stdio.h> #include <hdf5.h> +#include <zlib.h> #ifdef HAVE_CBFLIB #ifdef HAVE_CBFLIB_CBF_H @@ -439,7 +440,7 @@ static char *cbf_strerr(int e) } -static int unpack_panels(struct image *image, signed int *data, int data_width, +static int unpack_panels(struct image *image, float *data, int data_width, int data_height) { int pi; @@ -587,7 +588,39 @@ static void cbf_fill_in_clen(struct detector *det, struct imagefile *f) } -static int read_cbf(struct imagefile *f, struct image *image) +static float *convert_sint32_float(int32_t *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 float *convert_sint16_float(int16_t *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 float *read_cbf_data(struct imagefile *f, int *w, int *h, cbf_handle *pcbfh) { cbf_handle cbfh; FILE *fh; @@ -596,29 +629,70 @@ static int read_cbf(struct imagefile *f, struct image *image) int binary_id, minelement, maxelement, elsigned, elunsigned; size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; const char *byteorder; - signed int *data; + void *data; + float *dataf; + void *buf = NULL; - if ( image->det == NULL ) { - ERROR("read_cbf() needs a geometry\n"); - return 1; + if ( f->type == IMAGEFILE_CBF ) { + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return NULL; + } + + } else if ( f->type == IMAGEFILE_CBFGZ ) { + + gzFile gzfh; + size_t len, len_read; + const size_t bufinc = 8*1024*1024; /* Allocate buffer in 8Mb chunks */ + size_t bufsz = bufinc; + + gzfh = gzopen(f->filename, "rb"); + if ( gzfh == NULL ) return NULL; + + /* Set larger buffer size for hopefully faster uncompression */ + gzbuffer(gzfh, 128*1024); + + buf = malloc(bufsz); + if ( buf == NULL ) return NULL; + + len = 0; + do { + + len_read = gzread(gzfh, buf+len, bufinc); + if ( len_read == -1 ) return NULL; + len += len_read; + + if ( len_read == bufinc ) { + bufsz += bufinc; + buf = realloc(buf, bufsz); + if ( buf == NULL ) return NULL; + } + + } while ( len_read == bufinc ); + + fh = fmemopen(buf, len, "rb"); + if ( fh == NULL ) return NULL; + + gzclose_r(gzfh); + + } else { + /* Don't know how we ended up here */ + return NULL; } if ( cbf_make_handle(&cbfh) ) { ERROR("Failed to allocate CBF handle\n"); - return 1; + return NULL; } - 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 */ + /* CBFlib will call 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; + return NULL; } /* Select row 0 in data column inside array_data */ @@ -639,40 +713,46 @@ static int read_cbf(struct imagefile *f, struct image *image) ERROR("Failed to read CBF array parameters: %s\n", err); free(err); cbf_free_handle(cbfh); - return 1; + return NULL; } if ( dimslow != 0 ) { ERROR("CBF data array is 3D - don't know what to do with it\n"); cbf_free_handle(cbfh); - return 1; + return NULL; } 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; + return NULL; } - if ( elsize != 4 ) { + if ( (elsize != 4) && (elsize != 2) ) { STATUS("Don't know what to do with element size %i\n", (int)elsize); cbf_free_handle(cbfh); - return 1; + return NULL; + } + + if ( !elsigned ) { + STATUS("Don't know what to do with unsigned data (yet)\n"); + cbf_free_handle(cbfh); + return NULL; } 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; + return NULL; } data = malloc(elsize*dimfast*dimmid); if ( data == NULL ) { ERROR("Failed to allocate memory for CBF data\n"); cbf_free_handle(cbfh); - return 1; + return NULL; } r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, @@ -682,10 +762,44 @@ static int read_cbf(struct imagefile *f, struct image *image) ERROR("Failed to read CBF array: %s\n", err); free(err); cbf_free_handle(cbfh); + return NULL; + } + + if ( elsize == 4 ) { + dataf = convert_sint32_float(data, dimfast, dimmid); + } else if ( elsize == 2 ) { + dataf = convert_sint16_float(data, dimfast, dimmid); + } else { + ERROR("Don't know how to convert element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return NULL; + } + + free(data); + + free(buf); /* Might be NULL */ + + *w = dimfast; + *h = dimmid; + *pcbfh = cbfh; + return dataf; +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + float *data; + int w, h; + cbf_handle cbfh; + + data = read_cbf_data(f, &w, &h, &cbfh); + if ( data == NULL ) { + ERROR("Failed to read CBF data\n"); return 1; } - unpack_panels(image, data, dimfast, dimmid); + unpack_panels(image, data, w, h); free(data); if ( image->beam != NULL ) { @@ -704,126 +818,26 @@ static int read_cbf(struct imagefile *f, struct image *image) } -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) { + float *data; + int w, h; 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); + data = read_cbf_data(f, &w, &h, &cbfh); if ( data == NULL ) { - ERROR("Failed to allocate memory for CBF data\n"); - cbf_free_handle(cbfh); + ERROR("Failed to read CBF data\n"); 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->det = simple_geometry(image, w, h); 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; - } + + image->dp[0] = data; if ( image->beam != NULL ) { cbf_fill_in_beam_parameters(image->beam, f, image); @@ -840,6 +854,7 @@ static int read_cbf_simple(struct imagefile *f, struct image *image) return 0; } + #else /* HAVE_CBFLIB */ static int read_cbf_simple(struct imagefile *f, struct image *image) @@ -881,6 +896,24 @@ signed int is_cbf_file(const char *filename) } +signed int is_cbfgz_file(const char *filename) +{ + gzFile gzfh; + char line[1024]; + + gzfh = gzopen(filename, "rb"); + if ( gzfh == NULL ) return -1; + if ( gzgets(gzfh, line, 1024) == NULL ) return -1; + gzclose_r(gzfh); + + if ( strstr(line, "CBF") == NULL ) { + return 0; + } + + return 1; +} + + struct imagefile *imagefile_open(const char *filename) { struct imagefile *f; @@ -903,6 +936,15 @@ struct imagefile *imagefile_open(const char *filename) f->type = IMAGEFILE_CBF; + } else if ( is_cbfgz_file(filename) ) { + + f->type = IMAGEFILE_CBFGZ; + + } else { + + ERROR("Unrecognised file type: %s\n", filename); + return NULL; + } f->filename = strdup(filename); @@ -917,6 +959,8 @@ int imagefile_read(struct imagefile *f, struct image *image, return hdf5_read2(f->hdfile, image, event, 0); } else if ( f->type == IMAGEFILE_CBF ) { return read_cbf(f, image); + } else if ( f->type == IMAGEFILE_CBFGZ ) { + return read_cbf(f, image); } else { ERROR("Unknown file type %i\n", f->type); return 1; @@ -930,8 +974,13 @@ int imagefile_read_simple(struct imagefile *f, struct image *image) { if ( f->type == IMAGEFILE_HDF5 ) { return hdf5_read(f->hdfile, image, NULL, 0); - } else { + } else if ( f->type == IMAGEFILE_CBF ) { + return read_cbf_simple(f, image); + } else if ( f->type == IMAGEFILE_CBFGZ ) { return read_cbf_simple(f, image); + } else { + ERROR("Unknown file type %i\n", f->type); + return 1; } } diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 40717ab3..cc1bcc5d 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -60,6 +60,7 @@ struct imagefile_field_list; * @SPECTRUM_TOPHAT: A top hat distribution of wavelengths * @SPECTRUM_SASE: A simulated SASE spectrum * @SPECTRUM_TWOCOLOUR: A spectrum containing two peaks + * @SPECTRUM_FROMFILE: An arbitrary spectrum read from input file * * A %SpectrumType represents a type of X-ray energy spectrum to use for * generating simulated data. @@ -67,7 +68,8 @@ struct imagefile_field_list; typedef enum { SPECTRUM_TOPHAT, SPECTRUM_SASE, - SPECTRUM_TWOCOLOUR + SPECTRUM_TWOCOLOUR, + SPECTRUM_FROMFILE } SpectrumType; @@ -114,17 +116,27 @@ struct imagefeature { enum imagefile_type { IMAGEFILE_HDF5, - IMAGEFILE_CBF + IMAGEFILE_CBF, + IMAGEFILE_CBFGZ }; /* An opaque type representing a list of image features */ typedef struct _imagefeaturelist ImageFeatureList; + +struct spectrum +{ + int n; + double *ks; /* 1/m */ + double *weights; +}; + + /* Structure describing a wavelength sample from a spectrum */ struct sample { - double k; + double k; /* 1/m */ double weight; }; @@ -145,9 +157,11 @@ struct beam_params /** * image: + * @hit: Non-zero if the frame was determined to be a "hit" * @crystals: Array of crystals in the image * @n_crystals: The number of crystals in the image * @indexed_by: Indexing method which indexed this pattern + * @n_indexing_tries: Number of times the indexer was tried before indexing * @det: Detector structure * @beam: Beam parameters structure * @filename: Filename for the image file @@ -168,7 +182,7 @@ struct beam_params * @avg_clen: Mean of camera length values for all panels * @spectrum: Spectrum information * @nsamples: Number of spectrum samples - * @spectrum_size: SIze of spectrum array + * @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 @@ -190,9 +204,11 @@ struct image { int **bad; /* Bad pixels by panel */ float **sat; /* Per-pixel saturation values */ + int hit; Crystal **crystals; int n_crystals; IndexingMethod indexed_by; + int n_indexing_tries; struct detector *det; struct beam_params *beam; /* The nominal beam parameters */ @@ -209,7 +225,11 @@ struct image { int serial; /* Monotonically ascending serial * number for this image */ - struct sample *spectrum; + struct spectrum *spectrum; /* Beam spectrum for pink beam data */ + + /* FIXME: These are only used in pattern_sim, which should be changed to + * use the "struct spectrum" from above later */ + struct sample *spectrum0; int nsamples; /* Number of wavelengths */ int spectrum_size; /* Size of "spectrum" */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index ff84c629..1476304d 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -57,6 +57,7 @@ #include "felix.h" #include "predict-refine.h" #include "taketwo.h" +#include "xgandalf.h" struct _indexingprivate @@ -66,6 +67,7 @@ struct _indexingprivate float tolerance[4]; struct taketwo_options *ttopts; + struct xgandalf_options *xgandalf_opts; int n_methods; IndexingMethod *methods; @@ -182,6 +184,10 @@ static char *base_indexer_str(IndexingMethod indm) strcpy(str, "taketwo"); break; + case INDEXING_XGANDALF: + strcpy(str, "xgandalf"); + break; + case INDEXING_SIMULATION : strcpy(str, "simulation"); break; @@ -219,6 +225,7 @@ static char *friendly_indexer_name(IndexingMethod m) static void *prepare_method(IndexingMethod *m, UnitCell *cell, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts) { char *str; @@ -259,6 +266,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, priv = taketwo_prepare(m, cell); break; + case INDEXING_XGANDALF : + priv = xgandalf_prepare(m, cell, xgandalf_opts); + break; + default : ERROR("Don't know how to prepare indexing method %i\n", *m); break; @@ -290,6 +301,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, struct detector *det, float *ltl, IndexingFlags flags, struct taketwo_options *ttopts, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts) { int i, n; @@ -386,6 +398,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, int j; ipriv->engine_private[i] = prepare_method(&methods[i], cell, + xgandalf_opts, felix_opts); if ( ipriv->engine_private[i] == NULL ) return NULL; @@ -411,6 +424,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i]; ipriv->ttopts = ttopts; + ipriv->xgandalf_opts = xgandalf_opts; STATUS("List of indexing methods:\n"); for ( i=0; i<n; i++ ) { @@ -467,6 +481,10 @@ void cleanup_indexing(IndexingPrivate *ipriv) taketwo_cleanup(ipriv->engine_private[n]); break; + case INDEXING_XGANDALF : + xgandalf_cleanup(ipriv->engine_private[n]); + break; + default : ERROR("Don't know how to clean up indexing method %i\n", ipriv->methods[n]); @@ -577,6 +595,10 @@ static int try_indexer(struct image *image, IndexingMethod indm, r = taketwo_index(image, ipriv->ttopts, mpriv); break; + case INDEXING_XGANDALF : + r = run_xgandalf(image, mpriv); + break; + default : ERROR("Unrecognised indexing method: %i\n", indm); return 0; @@ -584,7 +606,11 @@ static int try_indexer(struct image *image, IndexingMethod indm, } /* Stop a really difficult to debug situation in its tracks */ - assert(image->n_crystals - n_before == r); + if ( image->n_crystals - n_before != r ) { + ERROR("Whoops, indexer didn't return the right number " + "of crystals!\n"); + exit(1); + } /* For all the crystals found this time ... */ for ( i=0; i<r; i++ ) { @@ -649,7 +675,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, if ( compare_cells(crystal_get_cell(cr), crystal_get_cell(that_cr), - 0.1, deg2rad(5.0), NULL) ) + 0.1, deg2rad(0.5), NULL) ) { crystal_set_user_flag(cr, 1); } @@ -829,7 +855,10 @@ void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping) /* Stop now if the pattern is indexed (don't try again for more * crystals with a different indexing method) */ - if ( success ) break; + if ( success ) { + image->n_indexing_tries = ntry; + break; + } } @@ -945,6 +974,11 @@ IndexingMethod get_indm_from_string_2(const char *str, int *err) method = INDEXING_DEFAULTS_TAKETWO; have_method = 1; + } else if ( strcmp(bits[i], "xgandalf") == 0) { + if ( have_method ) return warn_method(str); + method = INDEXING_DEFAULTS_XGANDALF; + have_method = 1; + } else if ( strcmp(bits[i], "none") == 0) { if ( have_method ) return warn_method(str); method = INDEXING_NONE; @@ -1042,6 +1076,7 @@ char *detect_indexing_methods(UnitCell *cell) do_probe(dirax_probe, cell, methods); do_probe(asdf_probe, cell, methods); do_probe(xds_probe, cell, methods); + do_probe(xgandalf_probe, cell, methods); /* Don't automatically use TakeTwo or Felix (yet) */ //do_probe(taketwo_probe, cell, methods); //do_probe(felix_probe, cell, methods); diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 8ece9227..2099b4d7 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -56,6 +56,8 @@ #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ | INDEXING_USE_CELL_PARAMETERS) +#define INDEXING_DEFAULTS_XGANDALF (INDEXING_XGANDALF | INDEXING_USE_CELL_PARAMETERS) + /** * IndexingMethod: * @INDEXING_NONE: No indexing to be performed @@ -67,6 +69,7 @@ * @INDEXING_DEBUG: Results injector for debugging * @INDEXING_ASDF: Use in-built "asdf" indexer * @INDEXING_TAKETWO: Use in-built "taketwo" indexer + * @INDEXING_XGANDALF: Invoke XGANDALF * @INDEXING_ERROR: Special value for unrecognised indexing engine name * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to * guide the indexing process. @@ -90,6 +93,7 @@ typedef enum { INDEXING_DEBUG = 7, INDEXING_ASDF = 8, INDEXING_TAKETWO = 9, + INDEXING_XGANDALF = 10, INDEXING_ERROR = 255, /* Unrecognised indexing engine */ @@ -136,6 +140,7 @@ extern IndexingMethod get_indm_from_string_2(const char *method, int *err); #include "cell.h" #include "image.h" #include "taketwo.h" +#include "xgandalf.h" #include "felix.h" @@ -143,6 +148,7 @@ extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell, struct detector *det, float *ltl, IndexingFlags flags, struct taketwo_options *ttopts, + struct xgandalf_options *xgandalf_opts, struct felix_options *felix_opts); extern char *detect_indexing_methods(UnitCell *cell); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 9dc99162..d48dcef2 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -570,18 +570,20 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) cell_get_parameters(mosflm->mp->template, &a, &b, &c, &alpha, &beta, &gamma); cen = cell_get_centering(mosflm->mp->template); + /* Old MOSFLM simply ignores CELL and CENTERING subkeywords. + * So this doesn't do any harm. + * TODO: but still better to show WARNING if MOSFLM is old. */ + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " + "MAXCELL 1000 REFINE " + "CELL %.1f %.1f %.1f %.1f %.1f %.1f " + "CENTERING %c\n", + mosflm->sptfile, a*1e10, b*1e10, c*1e10, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma), + cen); } else { - cen = 'P'; - a = 0; /* Disables prior-cell algorithm in MOSFLM */ + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " + "MAXCELL 1000 REFINE\n", mosflm->sptfile); } - /* Old MOSFLM simply ignores CELL and CENTERING subkeywords. - * So this doesn't do any harm. - * TODO: but still better to show WARNING if MOSFLM is old. */ - snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 MAXCELL 1000 " - "REFINE " - "CELL %.1f %.1f %.1f %.1f %.1f %.1f CENTERING %c\n", - mosflm->sptfile, a*1e10, b*1e10, c*1e10, - rad2deg(alpha), rad2deg(beta), rad2deg(gamma), cen); mosflm_sendline(tmp, mosflm); break; @@ -844,6 +846,15 @@ void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell) *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; + if ( (cell != NULL) && (cell_get_centering(cell) == 'I') + && (cell_get_lattice_type(cell) == L_MONOCLINIC) ) + { + ERROR("WARNING: Mosflm always gives the monoclinic C cell in " + "preference to the monoclinic I cell choice.\n"); + ERROR("To get a higher indexing rate, convert your cell to the " + "monoclinic C cell choice.\n"); + } + mp = malloc(sizeof(struct mosflm_private)); if ( mp == NULL ) return NULL; @@ -875,20 +886,6 @@ static void chop_word(char *s) } -static int file_exists(const char *filename) -{ - struct stat s; - - if ( stat(filename, &s) != 0 ) { - if ( errno == ENOENT ) return 0; - ERROR("Failed to check for %s.\n", filename); - exit(1); - } - - return 1; -} - - const char *mosflm_probe(UnitCell *cell) { pid_t pid; @@ -899,15 +896,6 @@ const char *mosflm_probe(UnitCell *cell) int ok = 0; int l; - /* Mosflm will write mosflm.lp and SUMMARY when we test it, which we are - * are going to delete afterwards. Better check they don't exist first, - * in case they were important. */ - if ( file_exists("mosflm.lp") || file_exists("SUMMARY") ) { - ERROR("Please move or delete mosflm.lp and SUMMARY from the " - "working directory first.\n"); - exit(1); - } - pid = forkpty(&pty, NULL, NULL, NULL); if ( pid == -1 ) { return NULL; @@ -947,9 +935,6 @@ const char *mosflm_probe(UnitCell *cell) close(pty); waitpid(pid, &status, 0); - unlink("mosflm.lp"); - unlink("SUMMARY"); - if ( !ok ) return NULL; if ( cell_has_parameters(cell) ) return "mosflm-cell-nolatt,mosflm-latt-nocell"; diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 39fe8a79..9090d48f 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -298,35 +298,35 @@ static struct radial_stats* allocate_radial_stats(int num_rad_bins) rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rthreshold == NULL ) { - free(rstats); free(rstats->roffset); + free(rstats); return NULL; } rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->lthreshold == NULL ) { - free(rstats); free(rstats->rthreshold); free(rstats->roffset); + free(rstats); 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); + free(rstats); 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); + free(rstats); return NULL; } @@ -434,8 +434,8 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int)); if ( pkdata->npix == NULL ) { - free(pkdata); free(pkdata->npix); + free(pkdata); return NULL; } diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 61a94e05..cccd5d0b 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -319,7 +319,7 @@ static int pair_peaks(struct image *image, Crystal *cr, } -void refine_radius(Crystal *cr, struct image *image) +int refine_radius(Crystal *cr, struct image *image) { int n, n_acc; struct reflpeak *rps; @@ -328,13 +328,13 @@ void refine_radius(Crystal *cr, struct image *image) /* Maximum possible size */ rps = malloc(image_feature_count(image->features) * sizeof(struct reflpeak)); - if ( rps == NULL ) return; + if ( rps == NULL ) return 1; reflist = reflist_new(); n_acc = pair_peaks(image, cr, reflist, rps); if ( n_acc < 3 ) { free(rps); - return; + return 1; } crystal_set_reflections(cr, reflist); update_predictions(cr); @@ -347,6 +347,8 @@ void refine_radius(Crystal *cr, struct image *image) reflist_free(reflist); free(rps); + + return 0; } diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index c763d2ca..fe700f47 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -39,7 +39,7 @@ struct image; extern int refine_prediction(struct image *image, Crystal *cr); -extern void refine_radius(Crystal *cr, struct image *image); +extern int refine_radius(Crystal *cr, struct image *image); #endif /* PREDICT_REFINE_H */ diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c index 54c467b3..70548994 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -3,11 +3,11 @@ * * Utilities to complement the core reflist.c * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -648,6 +648,30 @@ RefList *copy_reflist(RefList *list) } +/** + * free_contribs: + * @list: A %RefList + * + * Goes through @list and frees all the reflection contribution structures. + **/ +void free_contribs(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + struct reflection_contributions *c; + c = get_contributions(refl); + free(c->contribs); + free(c->contrib_crystals); + free(c); + } +} + + static char *full_command_line(int argc, char *argv[]) { int i; diff --git a/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h index f64e9f51..c955491a 100644 --- a/libcrystfel/src/reflist-utils.h +++ b/libcrystfel/src/reflist-utils.h @@ -3,11 +3,11 @@ * * Utilities to complement the core reflist.c * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -69,6 +69,8 @@ extern RefList *res_cutoff(RefList *list, UnitCell *cell, extern RefList *copy_reflist(RefList *list); +extern void free_contribs(RefList *list); + extern void reflist_add_command_and_version(RefList *list, int argcv, char *argv[]); diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index a35aa575..f1d759ea 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2016 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -100,6 +100,9 @@ struct _refldata { double peak; double mean_bg; + /* Contributions */ + struct reflection_contributions *contribs; + /* User-specified temporary values */ double temp1; double temp2; @@ -587,6 +590,18 @@ int get_flag(const Reflection *refl) } +/** + * get_contributions: + * @refl: A %Reflection + * + * Returns: the reflection's contribution list + * + **/ +struct reflection_contributions *get_contributions(const Reflection *refl) +{ + return refl->data.contribs; +} + /********************************** Setters ***********************************/ /** @@ -634,6 +649,7 @@ void set_panel(Reflection *refl, struct panel *p) refl->data.panel = p; } + /** * set_khalf: * @refl: A %Reflection @@ -647,7 +663,6 @@ void set_khalf(Reflection *refl, double khalf) } - /** * set_kpred: * @refl: A %Reflection @@ -848,6 +863,21 @@ void set_flag(Reflection *refl, int flag) } +/** + * set_contributions: + * @refl: A %Reflection + * @contribs: Pointer to the contribution list + * + * Note that the pointer will be stored, not the contents of the structure. + * + **/ +void set_contributions(Reflection *refl, + struct reflection_contributions *contribs) +{ + refl->data.contribs = contribs; +} + + /********************************* Insertion **********************************/ static Reflection *rotate_once(Reflection *refl, int dir) diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 9f494474..30bbaa28 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2011-2017 Thomas White <taw@physics.org> + * 2011-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -38,6 +38,7 @@ #define GET_K(serial) ((((serial) & 0x000ffc00)>>10)-512) #define GET_L(serial) (((serial) & 0x000003ff)-512) + /** * RefList: * @@ -69,10 +70,21 @@ typedef struct _reflection Reflection; **/ typedef struct _reflistiterator RefListIterator; +#include "crystal.h" + #ifdef __cplusplus extern "C" { #endif +/* Structure representing the contributions to a merged reflection */ +struct reflection_contributions +{ + int n_contrib; + int max_contrib; + Reflection **contribs; + Crystal **contrib_crystals; +}; + /* Creation/deletion */ extern RefList *reflist_new(void); extern void reflist_free(RefList *list); @@ -105,6 +117,7 @@ extern double get_phase(const Reflection *refl, int *have_phase); extern double get_peak(const Reflection *refl); extern double get_mean_bg(const Reflection *refl); extern int get_flag(const Reflection *refl); +extern struct reflection_contributions *get_contributions(const Reflection *refl); /* Set */ extern void copy_data(Reflection *to, const Reflection *from); @@ -126,6 +139,8 @@ extern void set_mean_bg(Reflection *refl, double mean_bg); extern void set_symmetric_indices(Reflection *refl, signed int hs, signed int ks, signed int ls); extern void set_flag(Reflection *refl, int flag); +extern void set_contributions(Reflection *refl, + struct reflection_contributions *contribs); /* Insertion */ extern Reflection *add_refl(RefList *list, diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index c11e38df..bcd2627e 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -832,9 +832,13 @@ int write_chunk(Stream *st, struct image *i, struct imagefile *imfile, fprintf(st->fh, "Image serial number: %i\n", i->serial); + fprintf(st->fh, "hit = %i\n", i->hit); indexer = indexer_str(i->indexed_by); fprintf(st->fh, "indexed_by = %s\n", indexer); free(indexer); + if ( i->indexed_by != INDEXING_NONE ) { + fprintf(st->fh, "n_indexing_tries = %i\n", i->n_indexing_tries); + } fprintf(st->fh, "photon_energy_eV = %f\n", J_to_eV(ph_lambda_to_en(i->lambda))); @@ -1355,6 +1359,7 @@ static void read_audit_lines(Stream *st) ERROR("Failed to allocate memory for audit information\n"); return; } + st->audit_info[0] = '\0'; /* Read lines from stream until one of them starts with "-----", * then rewind to the start of that line */ diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 5f3b67f2..06cc368c 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -36,6 +36,7 @@ #include <string.h> #include <math.h> #include <assert.h> +#include <ctype.h> #include "symmetry.h" #include "utils.h" @@ -1656,19 +1657,45 @@ static IntegerMatrix *parse_symmetry_operation(const char *s) int c; size_t cl; - int sign = +1; - int nh = 0; - int nk = 0; - int nl = 0; - + signed int nh = 0; + signed int nk = 0; + signed int nl = 0; + signed int mult = 1; + int ndigit = 0; + signed int sign = +1; + + /* We have one expression something like "-2h+k" */ cl = strlen(els[i]); - for ( c=0; c<cl; c++ ) { - if ( els[i][c] == '-' ) sign = -1; - if ( els[i][c] == '+' ) sign = +1; - if ( els[i][c] == 'h' ) nh += sign; - if ( els[i][c] == 'k' ) nk += sign; - if ( els[i][c] == 'l' ) nl += sign; + + if ( els[i][c] == '-' ) sign *= -1; + if ( els[i][c] == 'h' ) { + nh = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( els[i][c] == 'k' ) { + nk = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( els[i][c] == 'l' ) { + nl = mult*sign; + mult = 1; + ndigit = 0; + sign = +1; + } + if ( isdigit(els[i][c]) ) { + if ( ndigit > 0 ) { + mult *= 10; + mult += els[i][c] - '0'; + } else { + mult *= els[i][c] - '0'; + } + ndigit++; + } } intmat_set(m, i, 0, nh); diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index a759ff15..2c756ad7 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -185,6 +185,8 @@ static inline int within_tolerance(double a, double b, double percent) /* Photon energy (eV) to wavelength (m) */ #define ph_eV_to_lambda(a) ph_en_to_lambda(eV_to_J(a)) +/* Photon energy (eV) to k (1/m) */ +#define ph_eV_to_k(a) ((a)*ELECTRON_CHARGE/PLANCK/C_VACUO) /* ------------------------------ Message macros ---------------------------- */ diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c index 01c12dca..cef619c2 100644 --- a/libcrystfel/src/xds.c +++ b/libcrystfel/src/xds.c @@ -73,18 +73,43 @@ struct xds_private }; +/* Essentially the reverse of spacegroup_for_lattice(), below */ +static int convert_spacegroup_number(int spg, LatticeType *lt, char *cen, + char *ua) +{ + switch ( spg ) { + + case 1: *lt = L_TRICLINIC; *cen = 'P'; *ua = '*'; return 0; + case 3: *lt = L_MONOCLINIC; *cen = 'P'; *ua = 'b'; return 0; + case 5: *lt = L_MONOCLINIC; *cen = 'C'; *ua = 'b'; return 0; + case 16: *lt = L_ORTHORHOMBIC; *cen = 'P'; *ua = '*'; return 0; + case 21: *lt = L_ORTHORHOMBIC; *cen = 'C'; *ua = '*'; return 0; + case 22: *lt = L_ORTHORHOMBIC; *cen = 'F'; *ua = '*'; return 0; + case 23: *lt = L_ORTHORHOMBIC; *cen = 'I'; *ua = '*'; return 0; + case 75: *lt = L_TETRAGONAL; *cen = 'P'; *ua = 'c'; return 0; + case 79: *lt = L_TETRAGONAL; *cen = 'I'; *ua = 'c'; return 0; + case 143: *lt = L_HEXAGONAL; *cen = 'P'; *ua = 'c'; return 0; + case 146: *lt = L_HEXAGONAL; *cen = 'H'; *ua = 'c'; return 0; + case 195: *lt = L_CUBIC; *cen = 'P'; *ua = '*'; return 0; + case 196: *lt = L_CUBIC; *cen = 'F'; *ua = '*'; return 0; + case 197: *lt = L_CUBIC; *cen = 'I'; *ua = '*'; return 0; + default: return 1; + + } +} + + static int read_cell(struct image *image) { FILE * fh; - float axstar, aystar, azstar; - float bxstar, bystar, bzstar; - float cxstar, cystar, czstar; - char asx[11], asy[11], asz[11]; - char bsx[11], bsy[11], bsz[11]; - char csx[11], csy[11], csz[11]; + float ax, ay, az; + float bx, by, bz; + float cx, cy, cz; + int spg; char *rval, line[1024]; - int r; UnitCell *cell; + LatticeType latticetype; + char centering, ua; Crystal *cr; fh = fopen("IDXREF.LP", "r"); @@ -97,74 +122,66 @@ static int read_cell(struct image *image) return 0; } - } while ( strcmp(line, " # COORDINATES OF REC. BASIS VECTOR" - " LENGTH 1/LENGTH\n") != 0 ); + } while ( strcmp(line, " ***** DIFFRACTION PARAMETERS USED AT START OF " + "INTEGRATION *****\n") != 0 ); - /* Free line after chunk */ - rval = fgets(line, 1023, fh); - if ( rval == NULL ) { + /* Find and read space group number */ + do { + rval = fgets(line, 1023, fh); + if ( rval == NULL ) { + fclose(fh); + return 0; + } + } while ( strncmp(line, " SPACE GROUP NUMBER ", 20) != 0 ); + sscanf(line+20, "%i\n", &spg); + + /* Find and read a */ + do { + rval = fgets(line, 1023, fh); + if ( rval == NULL ) { + fclose(fh); + return 0; + } + } while ( strncmp(line, " COORDINATES OF UNIT CELL A-AXIS ", 33) != 0 ); + if ( sscanf(line+33, "%f %f %f\n", &ax, &ay, &az) < 3 ) { fclose(fh); return 0; } - /* Get first vector */ + /* Read b */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - if ( line[4] != '1' ) { - ERROR("No first vector from XDS.\n"); + if ( sscanf(line+33, "%f %f %f\n", &bx, &by, &bz) < 3 ) { + fclose(fh); return 0; } - memcpy(asx, line+7, 10); asx[10] = '\0'; - memcpy(asy, line+17, 10); asy[10] = '\0'; - memcpy(asz, line+27, 10); asz[10] = '\0'; - /* Get second vector */ + /* Read c */ rval = fgets(line, 1023, fh); if ( rval == NULL ) { fclose(fh); return 0; } - if ( line[4] != '2' ) { - ERROR("No second vector from XDS.\n"); + if ( sscanf(line+33, "%f %f %f\n", &cx, &cy, &cz) < 3 ) { + fclose(fh); return 0; } - memcpy(bsx, line+7, 10); bsx[10] = '\0'; - memcpy(bsy, line+17, 10); bsy[10] = '\0'; - memcpy(bsz, line+27, 10); bsz[10] = '\0'; - /* Get third vector */ - rval = fgets(line, 1023, fh); - fclose(fh); - if ( rval == NULL ) return 0; - if ( line[4] != '3' ) return 0; /* No error message this time - * - happens a lot */ - memcpy(csx, line+7, 10); csx[10] = '\0'; - memcpy(csy, line+17, 10); csy[10] = '\0'; - memcpy(csz, line+27, 10); csz[10] = '\0'; - - r = sscanf(asx, "%f", &cxstar); - r += sscanf(asy, "%f", &cystar); - r += sscanf(asz, "%f", &czstar); - r += sscanf(bsx, "%f", &bxstar); - r += sscanf(bsy, "%f", &bystar); - r += sscanf(bsz, "%f", &bzstar); - r += sscanf(csx, "%f", &axstar); - r += sscanf(csy, "%f", &aystar); - r += sscanf(csz, "%f", &azstar); - - if ( r != 9 ) { - STATUS("Fewer than 9 parameters found in NEWMAT file.\n"); + cell = cell_new(); + cell_set_cartesian(cell, + ax*1e-10, ay*1e-10, az*1e-10, + bx*1e-10, by*1e-10, bz*1e-10, + -cx*1e-10, -cy*1e-10, -cz*1e-10); + if ( convert_spacegroup_number(spg, &latticetype, ¢ering, &ua) ) { + ERROR("Failed to convert XDS space group number (%i)\n", spg); return 0; } - - cell = cell_new(); - cell_set_reciprocal(cell, - axstar*10e9, aystar*10e9, azstar*10e9, - bxstar*10e9, bystar*10e9, bzstar*10e9, - -cxstar*10e9, -cystar*10e9, -czstar*10e9); + cell_set_lattice_type(cell, latticetype); + cell_set_centering(cell, centering); + cell_set_unique_axis(cell, ua); cr = crystal_new(); if ( cr == NULL ) { @@ -235,13 +252,15 @@ static const char *spacegroup_for_lattice(UnitCell *cell) switch ( latt ) { case L_TRICLINIC : - g = "1"; + if ( centering == 'P' ) { + g = "1"; + } break; case L_MONOCLINIC : if ( centering == 'P' ) { g = "3"; - } else { + } else if ( centering == 'C' ) { g = "5"; } break; @@ -250,10 +269,10 @@ static const char *spacegroup_for_lattice(UnitCell *cell) if ( centering == 'P' ) { g = "16"; } else if ( centering == 'C' ) { - g = "20"; + g = "21"; } else if ( centering == 'F' ) { g = "22"; - } else { + } else if ( centering == 'I' ) { g = "23"; } break; @@ -261,34 +280,34 @@ static const char *spacegroup_for_lattice(UnitCell *cell) case L_TETRAGONAL : if ( centering == 'P' ) { g = "75"; - } else { + } else if ( centering == 'I' ) { g = "79"; } break; + /* Unfortunately, XDS only does "hexagonal H" */ case L_RHOMBOHEDRAL : + return NULL; + + case L_HEXAGONAL : if ( centering == 'P' ) { g = "143"; - } else { + } + if ( centering == 'H' ) { g = "146"; } break; - case L_HEXAGONAL : - g = "168"; - break; - case L_CUBIC : if ( centering == 'P' ) { g = "195"; } else if ( centering == 'F' ) { g = "196"; - } else { + } else if ( centering == 'I' ) { g = "197"; } break; } - assert(g != NULL); return g; } @@ -348,7 +367,7 @@ static int write_inp(struct image *image, struct xds_private *xp) fprintf(fh, "DETECTOR= CSPAD\n"); fprintf(fh, "MINIMUM_VALID_PIXEL_VALUE= 1\n"); fprintf(fh, "OVERLOAD= 200000000\n"); - fprintf(fh, "INDEX_ERROR= 0.05\n"); + fprintf(fh, "INDEX_ERROR= 0.2\n"); //fprintf(fh, "INDEX_QUALITY= 0.5\n"); fprintf(fh, "REFINE(IDXREF)= CELL ORIENTATION\n"); //fprintf(fh, "MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT= 1\n"); @@ -442,6 +461,11 @@ void *xds_prepare(IndexingMethod *indm, UnitCell *cell) return NULL; } + if ( (cell != NULL) && (spacegroup_for_lattice(cell) == NULL) ) { + ERROR("Don't know how to ask XDS for your cell.\n"); + return NULL; + } + xp = calloc(1, sizeof(*xp)); if ( xp == NULL ) return NULL; diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c new file mode 100644 index 00000000..2434ce2e --- /dev/null +++ b/libcrystfel/src/xgandalf.c @@ -0,0 +1,335 @@ +/* + * xgandalf.c + * + * Interface to XGANDALF indexer + * + * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@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/>. + * + */ + +#include "xgandalf.h" + +#ifdef HAVE_XGANDALF +#include <stdlib.h> + +#include "utils.h" +#include "cell-utils.h" +#include "peaks.h" + +#include "xgandalf/adaptions/crystfel/Lattice.h" +#include "xgandalf/adaptions/crystfel/ExperimentSettings.h" +#include "xgandalf/adaptions/crystfel/IndexerPlain.h" + +struct xgandalf_private_data { + IndexerPlain *indexer; + reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A; + + IndexingMethod indm; + UnitCell *cellTemplate; + Lattice_t sampleRealLattice_A; //same as cellTemplate + UnitCellTransformation *uncenteringTransformation; + LatticeTransform_t latticeReductionTransform; +}; + +#define FAKE_DETECTOR_DISTANCE (0.1) +#define FAKE_DETECTOR_RADIUS (0.1) +#define FAKE_BEAM_ENERGY (1) +#define FAKE_DIVERGENCE_ANGLE_DEG (0.05) +#define FAKE_NON_MONOCHROMATICITY (0.005) +#define FAKE_REFLECTION_RADIUS (0.0001) + +#define MAX_ASSEMBLED_LATTICES_COUNT (10) + +static void reduceCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform); +static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform); +static void makeRightHanded(UnitCell* cell); + +int run_xgandalf(struct image *image, void *ipriv) +{ + struct xgandalf_private_data *xgandalf_private_data = (struct xgandalf_private_data*) ipriv; + reciprocalPeaks_1_per_A_t *reciprocalPeaks_1_per_A = &(xgandalf_private_data->reciprocalPeaks_1_per_A); + + int peakCountMax = image_feature_count(image->features); + reciprocalPeaks_1_per_A->peakCount = 0; + for (int i = 0; i < peakCountMax && i < MAX_PEAK_COUNT_FOR_INDEXER; i++) { + struct imagefeature *f; + f = image_get_feature(image->features, i); + if (f == NULL) { + continue; + } + + reciprocalPeaks_1_per_A->coordinates_x[reciprocalPeaks_1_per_A->peakCount] = f->rx * 1e-10; + reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount] = f->ry * 1e-10; + reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount] = f->rz * 1e-10; + reciprocalPeaks_1_per_A->peakCount++; + } + + Lattice_t assembledLattices[MAX_ASSEMBLED_LATTICES_COUNT]; + int assembledLatticesCount; + IndexerPlain_index(xgandalf_private_data->indexer, + assembledLattices, + &assembledLatticesCount, + MAX_ASSEMBLED_LATTICES_COUNT, + *reciprocalPeaks_1_per_A, + NULL); + + if (assembledLatticesCount > 0) { //no multi-lattice at the moment + assembledLatticesCount = 1; + } + + int goodLatticesCount = assembledLatticesCount; + for (int i = 0; i < assembledLatticesCount && i < 1; i++) { + reorderLattice(&(xgandalf_private_data->sampleRealLattice_A), + &assembledLattices[i]); + + UnitCell *uc; + uc = cell_new(); + + Lattice_t *l = &assembledLattices[i]; + + cell_set_cartesian(uc, l->ax * 1e-10, l->ay * 1e-10, l->az * 1e-10, + l->bx * 1e-10, l->by * 1e-10, l->bz * 1e-10, + l->cx * 1e-10, l->cy * 1e-10, l->cz * 1e-10); + makeRightHanded(uc); + + if(xgandalf_private_data->cellTemplate != NULL){ + restoreCell(uc, &xgandalf_private_data->latticeReductionTransform); + + UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation); + cell_free(uc); + uc = new_cell_trans; + + cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(xgandalf_private_data->cellTemplate)); + cell_set_centering(new_cell_trans, cell_get_centering(xgandalf_private_data->cellTemplate)); + cell_set_unique_axis(new_cell_trans, cell_get_unique_axis(xgandalf_private_data->cellTemplate)); + } + + if (validate_cell(uc)) { + STATUS("Problem with returned cell!\n"); + } + + Crystal *cr = crystal_new(); + if (cr == NULL) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + crystal_set_cell(cr, uc); + image_add_crystal(image, cr); + + } + + return goodLatticesCount; +} + +void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts) +{ + struct xgandalf_private_data *xgandalf_private_data = malloc(sizeof(struct xgandalf_private_data)); + allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A)); + xgandalf_private_data->indm = *indm; + xgandalf_private_data->cellTemplate = NULL; + xgandalf_private_data->uncenteringTransformation = NULL; + + float tolerance = xgandalf_opts->tolerance; + samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch; + gradientDescentIterationsCount_t gradientDescentIterationsCount = xgandalf_opts->grad_desc_iterations; + + if (*indm & INDEXING_USE_CELL_PARAMETERS) { + + xgandalf_private_data->cellTemplate = cell; + + UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation); + + UnitCell *uc = cell_new_from_cell(primitiveCell); + reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform); + + cell_free(uc); + + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + int ret = cell_get_reciprocal(primitiveCell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + if (ret != 0) { + ERROR("cell_get_reciprocal did not finish properly!"); + } + + Lattice_t sampleReciprocalLattice_1_per_A = { + .ax = asx * 1e-10, .ay = asy * 1e-10, .az = asz * 1e-10, + .bx = bsx * 1e-10, .by = bsy * 1e-10, .bz = bsz * 1e-10, + .cx = csx * 1e-10, .cy = csy * 1e-10, .cz = csz * 1e-10 }; + + double ax, ay, az, bx, by, bz, cx, cy, cz; + ret = cell_get_cartesian(primitiveCell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + if (ret != 0) { + ERROR("cell_get_cartesian did not finish properly!"); + } + Lattice_t sampleRealLattice_A = { + .ax = ax * 1e10, .ay = ay * 1e10, .az = az * 1e10, + .bx = bx * 1e10, .by = by * 1e10, .bz = bz * 1e10, + .cx = cx * 1e10, .cy = cy * 1e10, .cz = cz * 1e10 }; + xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A; + + ExperimentSettings *experimentSettings = + ExperimentSettings_new(FAKE_BEAM_ENERGY, + FAKE_DETECTOR_DISTANCE, + FAKE_DETECTOR_RADIUS, + FAKE_DIVERGENCE_ANGLE_DEG, + FAKE_NON_MONOCHROMATICITY, + sampleReciprocalLattice_1_per_A, + tolerance, + FAKE_REFLECTION_RADIUS); + + xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings); + IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer, + samplingPitch); + IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer, + gradientDescentIterationsCount); + + if (xgandalf_opts->no_deviation_from_provided_cell) { + IndexerPlain_setRefineWithExactLattice(xgandalf_private_data->indexer, 1); + } + + ExperimentSettings_delete(experimentSettings); + cell_free(primitiveCell); + + } else { + + Lattice_t sampleRealLattice_A = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + xgandalf_private_data->sampleRealLattice_A = sampleRealLattice_A; + + ExperimentSettings *experimentSettings = + ExperimentSettings_new_nolatt(FAKE_BEAM_ENERGY, + FAKE_DETECTOR_DISTANCE, + FAKE_DETECTOR_RADIUS, + FAKE_DIVERGENCE_ANGLE_DEG, + FAKE_NON_MONOCHROMATICITY, + xgandalf_opts->minLatticeVectorLength_A, + xgandalf_opts->maxLatticeVectorLength_A, + FAKE_REFLECTION_RADIUS); + + xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings); + IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer, + samplingPitch); + IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer, + gradientDescentIterationsCount); + + ExperimentSettings_delete(experimentSettings); + } + + /* Flags that XGANDALF knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_USE_CELL_PARAMETERS; + + return xgandalf_private_data; +} + + +void xgandalf_cleanup(void *pp) +{ + struct xgandalf_private_data *xgandalf_private_data = pp; + + freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A); + IndexerPlain_delete(xgandalf_private_data->indexer); + if(xgandalf_private_data->uncenteringTransformation != NULL){ + tfn_free(xgandalf_private_data->uncenteringTransformation); + } +} + +static void reduceCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz }; + + reduceLattice(&l, appliedReductionTransform); + + cell_set_cartesian(cell, l.ax, l.ay, l.az, + l.bx, l.by, l.bz, + l.cx, l.cy, l.cz); + + makeRightHanded(cell); +} + +static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform) +{ + + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz }; + + restoreLattice(&l, appliedReductionTransform); + + cell_set_cartesian(cell, l.ax, l.ay, l.az, + l.bx, l.by, l.bz, + l.cx, l.cy, l.cz); + + makeRightHanded(cell); +} + +static void makeRightHanded(UnitCell *cell) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + if ( !right_handed(cell) ) { + cell_set_cartesian(cell, -ax, -ay, -az, -bx, -by, -bz, -cx, -cy, -cz); + } +} + + +const char *xgandalf_probe(UnitCell *cell) +{ + return "xgandalf"; +} + +#else + +int run_xgandalf(struct image *image, void *ipriv) +{ + ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n"); + return 0; +} + + +void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts) +{ + ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n"); + ERROR("To use XGANDALF indexing, recompile with XGANDALF.\n"); + return NULL; +} + + +void xgandalf_cleanup(void *pp) +{ +} + + +const char *xgandalf_probe(UnitCell *cell) +{ + return NULL; +} + +#endif // HAVE_XGANDALF diff --git a/libcrystfel/src/xgandalf.h b/libcrystfel/src/xgandalf.h new file mode 100644 index 00000000..1aced417 --- /dev/null +++ b/libcrystfel/src/xgandalf.h @@ -0,0 +1,58 @@ +/* + * xgandalf.h + * + * Interface to XGANDALF indexer + * + * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@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 LIBCRYSTFEL_SRC_XGANDALF_H +#define LIBCRYSTFEL_SRC_XGANDALF_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stddef.h> + +struct xgandalf_options { + unsigned int sampling_pitch; + unsigned int grad_desc_iterations; + float tolerance; + unsigned int no_deviation_from_provided_cell; + float minLatticeVectorLength_A; + float maxLatticeVectorLength_A; +}; + +#include "index.h" + +extern int run_xgandalf(struct image *image, void *ipriv); + +extern void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, + struct xgandalf_options *xgandalf_opts); + +extern void xgandalf_cleanup(void *pp); +extern const char *xgandalf_probe(UnitCell *cell); + + +#endif /* LIBCRYSTFEL_SRC_XGANDALF_H */ |