aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/CMakeLists.txt14
-rw-r--r--libcrystfel/src/cell-utils.c5
-rw-r--r--libcrystfel/src/crystal.h3
-rw-r--r--libcrystfel/src/detector.c13
-rw-r--r--libcrystfel/src/hdf5-file.c4
-rw-r--r--libcrystfel/src/image.c309
-rw-r--r--libcrystfel/src/image.h30
-rw-r--r--libcrystfel/src/index.c41
-rw-r--r--libcrystfel/src/index.h6
-rw-r--r--libcrystfel/src/mosflm.c57
-rw-r--r--libcrystfel/src/peakfinder8.c10
-rw-r--r--libcrystfel/src/predict-refine.c8
-rw-r--r--libcrystfel/src/predict-refine.h2
-rw-r--r--libcrystfel/src/reflist-utils.c28
-rw-r--r--libcrystfel/src/reflist-utils.h6
-rw-r--r--libcrystfel/src/reflist.c36
-rw-r--r--libcrystfel/src/reflist.h19
-rw-r--r--libcrystfel/src/stream.c5
-rw-r--r--libcrystfel/src/symmetry.c49
-rw-r--r--libcrystfel/src/utils.h2
-rw-r--r--libcrystfel/src/xds.c154
-rw-r--r--libcrystfel/src/xgandalf.c335
-rw-r--r--libcrystfel/src/xgandalf.h58
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, &centering, &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 */