Merge branch 'tom/xgandalf'

9 files changed, 529 insertions, 11 deletions

diff --git a/AUTHORS b/AUTHORS
index 88d5ec83..e27035d5 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -109,3 +109,4 @@ Contributors
 
 * Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
   peakfinder9
+  xgandalf
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 37e80801..ab581120 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -86,6 +86,10 @@ Invoke XDS, and use its REFIDX procedure to attempt to index the pattern.
 .PD
 Use the TakeTwo algorithm.  See Ginn et al., Acta Cryst. (2016). D72, 956-965.
 
+.IP \fBxgandalf\fR
+.PD
+Invoke XGANDALF - eXtended GrAdient Descent Algorithm for Lattice Finding. Xgandalf must be installed in order to be able to use it.
+
 .PP
 You can add one or more of the following to the above indexing methods, to control what information should be provided to them.  Note that indexamajig performs a series of checks on the indexing results, including checking that the result is consistent with the target unit cell parameters.  To get completely "raw" indexing, you need to disable these checks (see below) \fBand\fR not provide prior information.
 
@@ -397,6 +401,26 @@ The defaults are: \fB--taketwo-member-threshold=20\fR, \fB--taketwo-len-tolernac
 .PD 0
 These set low-level parameters for the Felix indexing algorithm.
 
+.PD 0
+.IP \fB--xgandalf-sampling-pitch=\fIn\fR
+.IP \fB--xgandalf-grad-desc-iterations=\fIn\fR
+.IP \fB--xgandalf-tolerance=\fIn\fR
+.IP \fB--xgandalf-no-deviation-from-provided-cell\fR
+.IP \fB--xgandalf-max-lattice-vector-length=\fIn\fR
+.IP \fB--xgandalf-min-lattice-vector-length=\fIn\fR
+.PD
+These set low-level parameters for the XGANDALF indexing algorithm.
+.IP
+\fB--xgandalf-sampling-pitch\fR selects how dense the reciprocal space is sampled. [0-4]: extremelyLoose to extremelyDense. [5-7]: standardWithSeondaryMillerIndices to extremelyDenseWithSeondaryMillerIndices. Default is 6 (denseWithSeondaryMillerIndices).
+.IP
+\fB--xgandalf-grad-desc-iterations\fR  selects how many gradient descent iterations are performed. [0-5]: veryFew to extremelyMany. Default is 4 (manyMany).
+.IP
+\fB--xgandalf-tolerance\fR relative tolerance of the lattice vectors. Default is 0.02.
+.IP
+\fB--xgandalf-no-deviation-from-provided-cell\fR if a prior unit cell was provided, and this flag is set, the found unit cell will have exactly the same size as the provided one.
+.IP
+\fB--xgandalf-min-lattice-vector-length\fR and \fB--xgandalf-min-lattice-vector-length\fR minimum and maximum possible lattice vector lengths (unit is A). Used for fitting without prior lattice as starting point for gradient descent, so the final minimum lattice vector length can be smaller/highier as min/max. Note: This is valid for the uncentered cell, i.e. the P-cell! Default is 30A and 250A respectively.
+
 .SH INTEGRATION OPTIONS
 .PD 0
 .IP \fB--integration=\fR\fImethod\fR
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index f6f25ec9..25cf66bc 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -68,6 +68,7 @@ set(LIBCRYSTFEL_SOURCES
     src/felix.c
     src/peakfinder8.c
     src/taketwo.c
+    src/xgandalf.c
 )
 
 if (HAVE_FFTW)
@@ -106,6 +107,7 @@ set(LIBCRYSTFEL_HEADERS
     src/felix.h
     src/peakfinder8.h
     src/taketwo.h
+    src/xgandalf.h
 )
 
 add_library(${PROJECT_NAME} SHARED
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index ff84c629..4ff4405a 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;
@@ -945,6 +967,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 +1069,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/xgandalf.c b/libcrystfel/src/xgandalf.c
new file mode 100644
index 00000000..9a2e54ec
--- /dev/null
+++ b/libcrystfel/src/xgandalf.c
@@ -0,0 +1,306 @@
+/*
+ * 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 sampleRealLatticeReduced_A;   //same as cellTemplate
+};
+
+#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);
+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->sampleRealLatticeReduced_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 (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;
+
+	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, NULL);
+		reduceCell(primitiveCell);
+
+		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 sampleRealLatticeReduced_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->sampleRealLatticeReduced_A = sampleRealLatticeReduced_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 sampleRealLatticeReduced_A = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+		xgandalf_private_data->sampleRealLatticeReduced_A = sampleRealLatticeReduced_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);
+}
+
+
+static void reduceCell(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);
+
+	Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+	reduceLattice(&l);
+
+	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 */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 4a790a17..3ad2695f 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -110,13 +110,16 @@ static void show_help(const char *s)
 "                            (peakfinder8 only) Default: 0\n"
 "     --max-res=<n>         Maximum resolution for peak search (in pixels)\n"
 "                            (peakfinder8 only) Default: 1200\n"
-"     --min-snr-biggest-pix=<n>    (peakFinder9 only) min snr of the biggest pixel in "
-"                                  the peak\n"
-"     --min-snr-peak-pix=<n>       (peakFinder9 only) min snr of a peak pixel\n"
-"     --min-sig=<n>                (peakFinder9 only) minimum standard deviation of "
-"                                 the background\n"
-"     --min-peak-over-neighbour=<n>    (peakFinder9 only) just for speed. Biggest pixel"
-"                                     in peak must be n higher than this.\n"
+"     --min-snr-biggest-pix=<n>\n"
+"                           Minimum snr of the biggest pixel in the peak\n"
+"                            (peakfinder9 only)\n"
+"     --min-snr-peak-pix=<n>\n"
+"                           Minimum snr of a peak pixel (peakfinder9 only)\n"
+"     --min-sig=<n>         Minimum standard deviation of the background\n"
+"                            (peakfinder9 only)\n"
+"     --min-peak-over-neighbour=<n>\n"
+"                           Just for speed. Biggest pixel in peak must be n\n"
+"                            higher than this (peakfinder9 only).\n"
 "     --no-use-saturated    Reject saturated peaks\n"
 "     --no-revalidate       Don't re-integrate and check HDF5 peaks\n"
 "     --no-half-pixel-shift\n"
@@ -171,6 +174,26 @@ static void show_help(const char *s)
 "                            Default: 30\n"
 "     --felix-tthrange-min  Minimum 2theta to consider for indexing (degrees)\n"
 "                            Default: 0\n"
+"\n"
+"     --xgandalf-sampling-pitch\n"
+"                           Sampling pitch: 0 (loosest) to 4 (most dense)\n"
+"                            or with secondary Miller indices: 5 (loosest) to\n"
+"                            7 (most dense).  Default: 6\n"
+"     --xgandalf-grad-desc-iterations\n"
+"                           Gradient descent iterations: 0 (few) to 5 (many)\n"
+"                            Default: 4\n"
+"     --xgandalf-tolerance  Relative tolerance of the lattice vectors.\n"
+"                            Default is 0.02\n"
+"     --xgandalf-no-deviation-from-provided-cell\n"
+"                           Force the fitted cell to have the same lattice\n"
+"                            parameters as the provided one\n"
+"     --xgandalf-min-lattice-vector-length\n"
+"                           Minimum possible lattice vector length in A.\n"
+"                            Default: 30 A\n"
+"     --xgandalf-max-lattice-vector-length\n"
+"                           Maximum possible lattice vector length in A.\n"
+"                            Default: 250 A\n"
+"\n"
 "\nIntegration options:\n\n"
 "     --integration=<meth>  Integration method (rings,prof2d)-(cen,nocen)\n"
 "                            Default: rings-nocen\n"
@@ -310,6 +333,12 @@ int main(int argc, char *argv[])
 	iargs.taketwo_opts.len_tol = -1.0;
 	iargs.taketwo_opts.angle_tol = -1.0;
 	iargs.taketwo_opts.trace_tol = -1.0;
+	iargs.xgandalf_opts.sampling_pitch = 6;
+	iargs.xgandalf_opts.grad_desc_iterations = 4;
+	iargs.xgandalf_opts.tolerance = 0.02;
+	iargs.xgandalf_opts.no_deviation_from_provided_cell = 0;
+	iargs.xgandalf_opts.minLatticeVectorLength_A = 30;
+	iargs.xgandalf_opts.maxLatticeVectorLength_A = 250;
 	iargs.felix_opts.ttmin = -1.0;
 	iargs.felix_opts.ttmax = -1.0;
 	iargs.felix_opts.min_visits = 0;
@@ -417,10 +446,22 @@ int main(int argc, char *argv[])
 		{"serial-start",             1, NULL,        344},
 		{"felix-domega",             1, NULL,        345},
 		{"felix-max-internal-angle", 1, NULL,        346},
-		{"min-snr-biggest-pix"              ,1, NULL,347},
-		{"min-snr-peak-pix"                 ,1, NULL,348},
-		{"min-sig"                          ,1, NULL,349},
-		{"min-peak-over-neighbour"          ,1, NULL,350},
+		{"min-snr-biggest-pix",      1, NULL,        347},
+		{"min-snr-peak-pix",         1, NULL,        348},
+		{"min-sig",                  1, NULL,        349},
+		{"min-peak-over-neighbour",  1, NULL,        350},
+		{"xgandalf-sampling-pitch",                  1, NULL, 351},
+		{"xgandalf-sps",                             1, NULL, 351},
+		{"xgandalf-grad-desc-iterations",            1, NULL, 352},
+		{"xgandalf-gdis",                            1, NULL, 352},
+		{"xgandalf-tolerance",                       1, NULL, 353},
+		{"xgandalf-tol",                             1, NULL, 353},
+		{"xgandalf-no-deviation-from-provided-cell", 0, NULL, 354},
+		{"xgandalf-ndfpc",                           0, NULL, 354},
+		{"xgandalf-min-lattice-vector-length",       1, NULL, 355},
+		{"xgandalf-min-lvl",                         1, NULL, 355},
+		{"xgandalf-max-lattice-vector-length",       1, NULL, 356},
+		{"xgandalf-max-lvl",                         1, NULL, 356},
 
 		{0, 0, NULL, 0}
 	};
@@ -763,6 +804,55 @@ int main(int argc, char *argv[])
 			iargs.min_peak_over_neighbour = strtof(optarg, NULL);
 			break;
 
+			case 351:
+			if (sscanf(optarg, "%u", &iargs.xgandalf_opts.sampling_pitch) != 1)
+			{
+				ERROR("Invalid value for --xgandalf-sampling-pitch\n");
+				return 1;
+			}
+			break;
+
+			case 352:
+			if (sscanf(optarg, "%u", &iargs.xgandalf_opts.grad_desc_iterations) != 1)
+			{
+				ERROR("Invalid value for --xgandalf-grad-desc-iterations\n");
+				return 1;
+			}
+			break;
+
+			case 353:
+			if (sscanf(optarg, "%f", &iargs.xgandalf_opts.tolerance) != 1)
+			{
+				ERROR("Invalid value for --xgandalf-tolerance\n");
+				return 1;
+			}
+			break;
+
+			case 354:
+				iargs.xgandalf_opts.no_deviation_from_provided_cell = 1;
+			break;
+
+			case 355:
+			if (sscanf(optarg, "%f",
+					&iargs.xgandalf_opts.minLatticeVectorLength_A) != 1)
+			{
+				ERROR("Invalid value for "
+						"--xgandalf-min-lattice-vector-length\n");
+				return 1;
+			}
+			break;
+
+			case 356:
+			if (sscanf(optarg, "%f",
+					&iargs.xgandalf_opts.maxLatticeVectorLength_A) != 1)
+			{
+				ERROR("Invalid value for "
+						"--xgandalf-max-lattice-vector-length\n");
+				return 1;
+			}
+			break;
+
+
 			case 0 :
 			break;
 
@@ -1060,6 +1150,7 @@ int main(int argc, char *argv[])
 		iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det,
 		                             iargs.tols, flags,
 		                             &iargs.taketwo_opts,
+		                             &iargs.xgandalf_opts,
 		                             &iargs.felix_opts);
 		if ( iargs.ipriv == NULL ) {
 			ERROR("Failed to set up indexing system\n");
diff --git a/src/process_image.h b/src/process_image.h
index 52f063fd..2b58d7bc 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -41,6 +41,7 @@ struct index_args;
 #include "im-sandbox.h"
 #include "time-accounts.h"
 #include "taketwo.h"
+#include "xgandalf.h"
 #include "felix.h"
 
 
@@ -109,6 +110,7 @@ struct index_args
 	int overpredict;
 	int profile;  /* Whether or not to do wall clock profiling */
 	struct taketwo_options taketwo_opts;
+    struct xgandalf_options xgandalf_opts;
 	struct felix_options felix_opts;
 };