Merge branch 'tom/pinkindexer'

9 files changed, 698 insertions, 7 deletions

diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index e037abf6..10e0584d 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -12,6 +12,7 @@ pkg_search_module(FFTW fftw3)
 set(HAVE_CURSES ${CURSES_FOUND})
 set(HAVE_FFTW ${FFTW_FOUND})
 set(HAVE_XGANDALF ${XGANDALF_FOUND})
+set(HAVE_PINKINDEXER ${PINKINDEXER_FOUND})
 set(HAVE_FDIP ${FDIP_FOUND})
 
 # Find out where forkpty() is declared
@@ -65,6 +66,7 @@ set(LIBCRYSTFEL_SOURCES
     src/peakfinder8.c
     src/taketwo.c
     src/xgandalf.c
+    src/pinkindexer.c
     src/rational.c
     src/spectrum.c
     ${BISON_symopp_OUTPUTS}
@@ -106,6 +108,7 @@ set(LIBCRYSTFEL_HEADERS
     src/peakfinder8.h
     src/taketwo.h
     src/xgandalf.h
+    src/pinkindexer.h
     src/rational.h
     src/spectrum.h
 )
diff --git a/libcrystfel/config.h.cmake.in b/libcrystfel/config.h.cmake.in
index ea26c1a2..b7ad4c7a 100644
--- a/libcrystfel/config.h.cmake.in
+++ b/libcrystfel/config.h.cmake.in
@@ -4,6 +4,8 @@
 #cmakedefine HAVE_CPU_AFFINITY
 #cmakedefine HAVE_FFTW
 #cmakedefine HAVE_XGANDALF
+#cmakedefine HAVE_NBP
+#cmakedefine HAVE_PINKINDEXER
 #cmakedefine HAVE_FDIP
 #cmakedefine HAVE_CURSES
 
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 201b9041..4d008e65 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -1070,6 +1070,19 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam,
 			}
 		}
 
+	} else if ( strcmp(key, "photon_energy_bandwidth") == 0 ) {
+		if ( beam != NULL ) {
+			double v;
+			char *end;
+			v = strtod(val, &end);
+			if ( (val[0] != '\0') && (end[0] == '\0') ) {
+				beam->bandwidth = v;
+			} else {
+				ERROR("Invalid value for "
+				      "photon_energy_bandwidth\n");
+			}
+		}
+
 	} else if ( strcmp(key, "photon_energy_scale") == 0 ) {
 		if ( beam != NULL ) {
 			beam->photon_energy_scale = atof(val);
@@ -1202,6 +1215,7 @@ struct detector *get_detector_geometry_from_string(const char *string,
 		beam->photon_energy = 0.0;
 		beam->photon_energy_from = NULL;
 		beam->photon_energy_scale = 1.0;
+		beam->bandwidth = 0.00000001;
 	}
 
 	det->n_panels = 0;
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 5b2107d1..fc42a8de 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -332,7 +332,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 			/* Closest point on Ewald sphere.
 			 * Project zl to 0, bit of a hack... */
-			const double zlp0 = zl>0?zl:0;
+			const double zlp0 = zl<0?zl:0;
 			exerr2 = (x-xl)*(x-xl) + (y-yl)*(y-yl) + (z-zl)*(z-zl);
 
 			/* Weighted average between projected lattice point
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 706b06b1..04fb18a0 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -110,6 +110,8 @@ struct beam_params
 	char  *photon_energy_from;  /**< HDF5 dataset name */
 	double photon_energy_scale; /**< Scale factor for photon energy, if it
 	                             *   comes from an image header */
+	double bandwidth;           /**< FWHM bandwidth as a fraction of
+	                             *   wavelength */
 };
 
 
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index f6e81109..049c3ba0 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -58,6 +58,7 @@
 #include "predict-refine.h"
 #include "taketwo.h"
 #include "xgandalf.h"
+#include "pinkindexer.h"
 
 
 /** \file index.h */
@@ -71,6 +72,7 @@ struct _indexingprivate
 
 	struct taketwo_options *ttopts;
 	struct xgandalf_options *xgandalf_opts;
+	struct pinkIndexer_options *pinkIndexer_opts;
 
 	int n_methods;
 	IndexingMethod *methods;
@@ -195,6 +197,10 @@ static char *base_indexer_str(IndexingMethod indm)
 		strcpy(str, "xgandalf");
 		break;
 
+		case INDEXING_PINKINDEXER:
+		strcpy(str, "pinkIndexer");
+		break;
+
 		case INDEXING_SIMULATION :
 		strcpy(str, "simulation");
 		break;
@@ -232,7 +238,9 @@ static char *friendly_indexer_name(IndexingMethod m)
 
 
 static void *prepare_method(IndexingMethod *m, UnitCell *cell,
+                            struct detector *det, struct beam_params *beam,
                             struct xgandalf_options *xgandalf_opts,
+                            struct pinkIndexer_options* pinkIndexer_opts,
                             struct felix_options *felix_opts)
 {
 	char *str;
@@ -277,6 +285,11 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 		priv = xgandalf_prepare(m, cell, xgandalf_opts);
 		break;
 
+		case INDEXING_PINKINDEXER :
+		priv = pinkIndexer_prepare(m, cell, pinkIndexer_opts,
+		                           det, beam);
+		break;
+
 		default :
 		ERROR("Don't know how to prepare indexing method %i\n", *m);
 		break;
@@ -305,10 +318,11 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 
 
 IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
-                                struct detector *det, float *tols,
-                                IndexingFlags flags,
+                                struct detector *det, struct beam_params *beam,
+                                float *tols, IndexingFlags flags,
                                 struct taketwo_options *ttopts,
                                 struct xgandalf_options *xgandalf_opts,
+                                struct pinkIndexer_options *pinkIndexer_opts,
                                 struct felix_options *felix_opts)
 {
 	int i, n;
@@ -333,7 +347,6 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 			ERROR("The indexing method should contain only the method itself and ");
 			ERROR("prior information modifiers ('cell' or 'latt')\n");
 			ERROR("To disable prediction refinement ('norefine'), use --no-refine.\n");
-			ERROR("To check cell axes only ('axes'), use --no-cell-combinations.\n");
 			ERROR("To disable all unit cell checks ('raw'), use --no-check-cell.\n");
 			ERROR("To disable peak alignment check ('bad'), use --no-check-peaks.\n");
 			ERROR("To disable indexing retry ('noretry'), use --no-retry.\n");
@@ -403,11 +416,36 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 		int j;
 
 		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
+		                                          det, beam,
 		                                          xgandalf_opts,
+		                                          pinkIndexer_opts,
 		                                          felix_opts);
 
 		if ( ipriv->engine_private[i] == NULL ) return NULL;
 
+		if ( (methods[i] & INDEXING_METHOD_MASK) == INDEXING_PINKINDEXER ) {
+			if ( n > 1 ) {
+				ERROR("WARNING: Using PinkIndexer at the same "
+				      "time as other indexers is not "
+				      "recommended.\n");
+			}
+
+			if ( flags & INDEXING_CHECK_PEAKS ) {
+				ERROR("WARNING: Setting --no-check-peaks "
+				      "because PinkIndexer is in use.\n");
+			}
+			flags |= INDEXING_CHECK_PEAKS;
+			flags ^= INDEXING_CHECK_PEAKS;
+
+			if ( flags & INDEXING_REFINE ) {
+				ERROR("WARNING: Setting --no-refine because "
+				      "PinkIndexer is in use.\n");
+			}
+			flags |= INDEXING_REFINE;
+			flags ^= INDEXING_REFINE;
+		}
+
+
 		for ( j=0; j<i; j++ ) {
 			if ( methods[i] == methods[j] ) {
 				ERROR("Duplicate indexing method.\n");
@@ -430,6 +468,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 
 	ipriv->ttopts = ttopts;
 	ipriv->xgandalf_opts = xgandalf_opts;
+	ipriv->pinkIndexer_opts = pinkIndexer_opts;
 
 	STATUS("List of indexing methods:\n");
 	for ( i=0; i<n; i++ ) {
@@ -445,6 +484,13 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 }
 
 
+const IndexingMethod *indexing_methods(IndexingPrivate *p, int *n)
+{
+	*n = p->n_methods;
+	return p->methods;
+}
+
+
 void cleanup_indexing(IndexingPrivate *ipriv)
 {
 	int n;
@@ -490,6 +536,10 @@ void cleanup_indexing(IndexingPrivate *ipriv)
 			xgandalf_cleanup(ipriv->engine_private[n]);
 			break;
 
+			case INDEXING_PINKINDEXER :
+			pinkIndexer_cleanup(ipriv->engine_private[n]);
+			break;
+
 			default :
 			ERROR("Don't know how to clean up indexing method %i\n",
 			      ipriv->methods[n]);
@@ -602,6 +652,11 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		r = taketwo_index(image, ipriv->ttopts, mpriv);
 		break;
 
+		case INDEXING_PINKINDEXER :
+		set_last_task(last_task, "indexing:pinkindexer");
+		r = run_pinkIndexer(image, mpriv);
+		break;
+
 		case INDEXING_XGANDALF :
 		set_last_task(last_task, "indexing:xgandalf");
 		r = run_xgandalf(image, mpriv);
@@ -1000,6 +1055,13 @@ IndexingMethod get_indm_from_string_2(const char *str, int *err)
 			method = INDEXING_DEFAULTS_XGANDALF;
 			have_method = 1;
 
+		} else if ( (strcmp(bits[i], "pinkIndexer") == 0)
+		         || (strcmp(bits[i], "pinkindexer") == 0) )
+		{
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_PINKINDEXER;
+			have_method = 1;
+
 		} else if ( strcmp(bits[i], "none") == 0) {
 			if ( have_method ) return warn_method(str);
 			method = INDEXING_NONE;
@@ -1098,9 +1160,10 @@ char *detect_indexing_methods(UnitCell *cell)
 	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) */
+	/* Don't automatically use TakeTwo, Felix or PinkIndexer (yet) */
 	//do_probe(taketwo_probe, cell, methods);
 	//do_probe(felix_probe, cell, methods);
+	//do_probe(pinkIndexer_probe, cell, methods);
 
 	if ( strlen(methods) == 0 ) {
 		free(methods);
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index eb5cf7e4..5e9e40ad 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -63,6 +63,8 @@
 
 #define INDEXING_DEFAULTS_XGANDALF (INDEXING_XGANDALF | INDEXING_USE_CELL_PARAMETERS)
 
+#define INDEXING_DEFAULTS_PINKINDEXER (INDEXING_PINKINDEXER | INDEXING_USE_CELL_PARAMETERS)
+
 /**
  * An enumeration of all the available indexing methods.
  **/
@@ -79,6 +81,7 @@ typedef enum {
 	INDEXING_ASDF = 8,        /**< Use built-in ASDF algorithm */
 	INDEXING_TAKETWO = 9,     /**< Use built-in TakeTwo algorithm */
 	INDEXING_XGANDALF = 10,   /**< Use XGANDALF (via optional library) */
+	INDEXING_PINKINDEXER = 11,/**< Use PinkIndexer (via optional library) */
 
 	INDEXING_ERROR = 255,     /**< Special value for unrecognised indexing
 	                           *   engine */
@@ -146,16 +149,20 @@ extern IndexingMethod get_indm_from_string_2(const char *method, int *err);
 #include "image.h"
 #include "taketwo.h"
 #include "xgandalf.h"
+#include "pinkindexer.h"
 #include "felix.h"
 
-
 extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell,
-                                       struct detector *det, float *ltl,
+                                       struct detector *det,
+                                       struct beam_params *beam, float *ltl,
                                        IndexingFlags flags,
                                        struct taketwo_options *ttopts,
                                        struct xgandalf_options *xgandalf_opts,
+                                       struct pinkIndexer_options *pinkIndexer_opts,
                                        struct felix_options *felix_opts);
 
+extern const IndexingMethod *indexing_methods(IndexingPrivate *p, int *n);
+
 extern char *detect_indexing_methods(UnitCell *cell);
 
 extern void index_pattern(struct image *image, IndexingPrivate *ipriv);
diff --git a/libcrystfel/src/pinkindexer.c b/libcrystfel/src/pinkindexer.c
new file mode 100644
index 00000000..811aa3c8
--- /dev/null
+++ b/libcrystfel/src/pinkindexer.c
@@ -0,0 +1,535 @@
+/*
+ * pinkindexer.c
+ *
+ * Interface to PinkIndexer
+ *
+ * Copyright © 2017-2019 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017-2019 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 "pinkindexer.h"
+
+#ifdef HAVE_PINKINDEXER
+
+#include <stdlib.h>
+
+#include "utils.h"
+#include "cell-utils.h"
+#include "peaks.h"
+
+#include "pinkIndexer/adaptions/crystfel/Lattice.h"
+#include "pinkIndexer/adaptions/crystfel/ExperimentSettings.h"
+#include "pinkIndexer/adaptions/crystfel/PinkIndexer.h"
+
+#define MAX_MULTI_LATTICE_COUNT 8
+
+struct pinkIndexer_private_data {
+	PinkIndexer *pinkIndexer;
+	reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A;
+	float *intensities;
+
+	IndexingMethod indm;
+	UnitCell *cellTemplate;
+	int threadCount;
+	int multi;
+	int min_peaks;
+
+	int no_check_indexed;
+
+	IntegerMatrix *centeringTransformation;
+	LatticeTransform_t latticeReductionTransform;
+};
+
+//static void reduceCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform);
+//static void restoreCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform);
+static void reduceReciprocalCell(UnitCell* cell, LatticeTransform_t* appliedReductionTransform);
+static void restoreReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform);
+static void makeRightHanded(UnitCell* cell);
+static void update_detector(struct detector *det, double xoffs, double yoffs);
+
+int run_pinkIndexer(struct image *image, void *ipriv)
+{
+	struct pinkIndexer_private_data* pinkIndexer_private_data = (struct pinkIndexer_private_data*) ipriv;
+	reciprocalPeaks_1_per_A_t* reciprocalPeaks_1_per_A = &(pinkIndexer_private_data->reciprocalPeaks_1_per_A);
+	float *intensities = pinkIndexer_private_data->intensities;
+
+	int peakCountMax = image_feature_count(image->features);
+	if (peakCountMax < 5) {
+		int goodLatticesCount = 0;
+		return goodLatticesCount;
+	}
+	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->rz * 1e-10;
+		reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount] = f->rx * 1e-10;
+		reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount] = f->ry * 1e-10;
+		intensities[reciprocalPeaks_1_per_A->peakCount] = (float) (f->intensity);
+		reciprocalPeaks_1_per_A->peakCount++;
+	}
+	int indexed = 0;
+	Lattice_t indexedLattice[MAX_MULTI_LATTICE_COUNT];
+	float center_shift[MAX_MULTI_LATTICE_COUNT][2];
+
+	float maxRefinementDisbalance = 0.4;
+
+	do {
+		int peakCount = reciprocalPeaks_1_per_A->peakCount;
+		int matchedPeaksCount = PinkIndexer_indexPattern(pinkIndexer_private_data->pinkIndexer,
+		        &(indexedLattice[indexed]), center_shift[indexed], reciprocalPeaks_1_per_A, intensities,
+		        maxRefinementDisbalance,
+		        pinkIndexer_private_data->threadCount);
+
+		if ((matchedPeaksCount >= 25 && matchedPeaksCount >= peakCount * 0.30)
+		        || matchedPeaksCount >= peakCount * 0.4
+		        || matchedPeaksCount >= 70
+		        || pinkIndexer_private_data->no_check_indexed == 1)
+		                {
+			UnitCell *uc;
+			uc = cell_new();
+
+			Lattice_t *l = &(indexedLattice[indexed]);
+
+			cell_set_reciprocal(uc, l->ay * 1e10, l->az * 1e10, l->ax * 1e10,
+			        l->by * 1e10, l->bz * 1e10, l->bx * 1e10,
+			        l->cy * 1e10, l->cz * 1e10, l->cx * 1e10);
+
+			restoreReciprocalCell(uc, &pinkIndexer_private_data->latticeReductionTransform);
+
+			UnitCell *new_cell_trans = cell_transform_intmat(uc, pinkIndexer_private_data->centeringTransformation);
+			cell_free(uc);
+			uc = new_cell_trans;
+
+			cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(pinkIndexer_private_data->cellTemplate));
+			cell_set_centering(new_cell_trans, cell_get_centering(pinkIndexer_private_data->cellTemplate));
+			cell_set_unique_axis(new_cell_trans, cell_get_unique_axis(pinkIndexer_private_data->cellTemplate));
+
+			if (validate_cell(uc)) {
+				ERROR("pinkIndexer: 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);
+			crystal_set_det_shift(cr, center_shift[indexed][0], center_shift[indexed][1]);
+			update_detector(image->det, center_shift[indexed][0], center_shift[indexed][1]);
+			image_add_crystal(image, cr);
+			indexed++;
+
+		} else {
+			break;
+		}
+	} while (pinkIndexer_private_data->multi
+	        && indexed <= MAX_MULTI_LATTICE_COUNT
+	        && reciprocalPeaks_1_per_A->peakCount >= pinkIndexer_private_data->min_peaks);
+
+	return indexed;
+}
+
+void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+                          struct pinkIndexer_options *pinkIndexer_opts,
+                          struct detector *det, struct beam_params *beam)
+{
+	if ( beam->photon_energy_from != NULL ) {
+		ERROR("For pinkIndexer, the photon_energy must be defined as a "
+		      "constant in the geometry file\n");
+		return NULL;
+	}
+	if ( (det->panels[0].clen_from != NULL) && pinkIndexer_opts->refinement_type ==
+	        REFINEMENT_TYPE_firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed) {
+		ERROR("Using center refinement makes it necessary to have the detector distance fixed in the geometry file!");
+	}
+
+	struct pinkIndexer_private_data* pinkIndexer_private_data = malloc(sizeof(struct pinkIndexer_private_data));
+	allocReciprocalPeaks(&(pinkIndexer_private_data->reciprocalPeaks_1_per_A));
+	pinkIndexer_private_data->intensities = malloc(MAX_PEAK_COUNT_FOR_INDEXER * sizeof(float));
+	pinkIndexer_private_data->indm = *indm;
+	pinkIndexer_private_data->cellTemplate = cell;
+	pinkIndexer_private_data->threadCount = pinkIndexer_opts->thread_count;
+	pinkIndexer_private_data->multi = pinkIndexer_opts->multi;
+	pinkIndexer_private_data->min_peaks = pinkIndexer_opts->min_peaks;
+	pinkIndexer_private_data->no_check_indexed = pinkIndexer_opts->no_check_indexed;
+
+	UnitCell* primitiveCell = uncenter_cell(cell, &pinkIndexer_private_data->centeringTransformation, NULL);
+
+	//reduceCell(primitiveCell, &pinkIndexer_private_data->latticeReductionTransform);
+	reduceReciprocalCell(primitiveCell, &pinkIndexer_private_data->latticeReductionTransform);
+
+	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 lattice = { .ax = asz * 1e-10, .ay = asx * 1e-10, .az = asy * 1e-10,
+	        .bx = bsz * 1e-10, .by = bsx * 1e-10, .bz = bsy * 1e-10,
+	        .cx = csz * 1e-10, .cy = csx * 1e-10, .cz = csy * 1e-10 };
+
+	float detectorDistance_m;
+	if ( det->panels[0].clen_from != NULL ) {
+		detectorDistance_m =  0.25;  /* fake value */
+	} else {
+		detectorDistance_m = det->panels[0].clen + det->panels[0].coffset;
+	}
+
+	float beamEenergy_eV = beam->photon_energy;
+	float nonMonochromaticity = beam->bandwidth;
+	float reflectionRadius_1_per_A;
+	if (pinkIndexer_opts->reflectionRadius < 0) {
+		reflectionRadius_1_per_A = 0.02
+		        * sqrt(lattice.ax * lattice.ax + lattice.ay * lattice.ay + lattice.az * lattice.az);
+	}
+	else {
+		reflectionRadius_1_per_A = pinkIndexer_opts->reflectionRadius;
+	}
+
+	float divergenceAngle_deg = 0.01; //fake
+
+	float tolerance = pinkIndexer_opts->tolerance;
+	Lattice_t sampleReciprocalLattice_1_per_A = lattice;
+	float detectorRadius_m = 0.03; //fake, only for prediction
+	ExperimentSettings* experimentSettings = ExperimentSettings_new(beamEenergy_eV, detectorDistance_m,
+	        detectorRadius_m, divergenceAngle_deg, nonMonochromaticity, sampleReciprocalLattice_1_per_A, tolerance,
+	        reflectionRadius_1_per_A);
+
+	consideredPeaksCount_t consideredPeaksCount = pinkIndexer_opts->considered_peaks_count;
+	angleResolution_t angleResolution = pinkIndexer_opts->angle_resolution;
+	refinementType_t refinementType = pinkIndexer_opts->refinement_type;
+	float maxResolutionForIndexing_1_per_A = pinkIndexer_opts->maxResolutionForIndexing_1_per_A;
+	pinkIndexer_private_data->pinkIndexer = PinkIndexer_new(experimentSettings, consideredPeaksCount, angleResolution,
+	        refinementType,
+	        maxResolutionForIndexing_1_per_A);
+
+	ExperimentSettings_delete(experimentSettings);
+	cell_free(primitiveCell);
+
+	/* Flags that pinkIndexer knows about */
+	*indm &= INDEXING_METHOD_MASK
+	        | INDEXING_USE_CELL_PARAMETERS;
+
+	return pinkIndexer_private_data;
+}
+
+//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 reduceReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	cell_get_reciprocal(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_reciprocal(cell, l.ax, l.ay, l.az,
+	        l.bx, l.by, l.bz,
+	        l.cx, l.cy, l.cz);
+
+	makeRightHanded(cell);
+}
+
+static void restoreReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedReductionTransform)
+{
+
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	cell_get_reciprocal(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_reciprocal(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);
+	}
+}
+
+//hack for electron crystallography while crystal_set_det_shift is not working approprietly
+static void update_detector(struct detector *det, double xoffs, double yoffs)
+{
+	int i;
+
+	for (i = 0; i < det->n_panels; i++) {
+		struct panel *p = &det->panels[i];
+		p->cnx += xoffs * p->res;
+		p->cny += yoffs * p->res;
+	}
+}
+
+void pinkIndexer_cleanup(void *pp)
+{
+	struct pinkIndexer_private_data* pinkIndexer_private_data = (struct pinkIndexer_private_data*) pp;
+
+	freeReciprocalPeaks(pinkIndexer_private_data->reciprocalPeaks_1_per_A);
+	free(pinkIndexer_private_data->intensities);
+	intmat_free(pinkIndexer_private_data->centeringTransformation);
+	PinkIndexer_delete(pinkIndexer_private_data->pinkIndexer);
+}
+
+const char *pinkIndexer_probe(UnitCell *cell)
+{
+	return "pinkIndexer";
+}
+
+#else /* HAVE_PINKINDEXER */
+
+int run_pinkIndexer(struct image *image, void *ipriv)
+{
+	ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n");
+	return 0;
+}
+
+extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+		struct pinkIndexer_options *pinkIndexer_opts,
+		struct detector *det, struct beam_params *beam)
+{
+	ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n");
+	ERROR("To use PINKINDEXER indexing, recompile with PINKINDEXER.\n");
+	return NULL;
+}
+
+void pinkIndexer_cleanup(void *pp)
+{
+}
+
+const char *pinkIndexer_probe(UnitCell *cell)
+{
+	return NULL;
+}
+
+#endif /* HAVE_PINKINDEXER */
+
+static void show_help()
+{
+	printf("Parameters for the PinkIndexer indexing algorithm:\n"
+"     --pinkIndexer-considered-peaks-count\n"
+"                           Considered peaks count, 0 (fewest) to 4 (most)\n"
+"                            Default: 4\n"
+"     --pinkIndexer-angle-resolution\n"
+"                           Angle resolution, 0 (loosest) to 4 (most dense)\n"
+"                            Default: 2\n"
+"     --pinkIndexer-refinement-type\n"
+"                           Refinement type, 0 (none) to 5 (most accurate)\n"
+"                            Default: 1\n"
+"     --pinkIndexer-tolerance\n"
+"                           Relative tolerance of the lattice vectors.\n"
+"                            Default 0.06\n"
+"     --pinkIndexer-reflection-radius\n"
+"                           Radius of the reflections in reciprocal space.\n"
+"                            Specified in 1/A.  Default is 2%% of a*.\n"
+"     --pinkIndexer-max-resolution-for-indexing\n"
+"                           Measured in 1/A\n"
+"     --pinkIndexer-multi   Use pinkIndexers own multi indexing.\n"
+"     --pinkIndexer-thread-count\n"
+"                           Thread count for internal parallelization \n"
+"                            Default: 1\n"
+"     --pinkIndexer-no-check-indexed\n"
+"                           Disable internal check for correct indexing\n"
+"                            solutions\n"
+);
+}
+
+
+static error_t parse_arg(int key, char *arg, struct argp_state *state)
+{
+	float tmp;
+	struct pinkIndexer_options **opts_ptr = state->input;
+
+	switch ( key ) {
+
+		case ARGP_KEY_INIT :
+		*opts_ptr = malloc(sizeof(struct pinkIndexer_options));
+		if ( *opts_ptr == NULL ) return ENOMEM;
+		(*opts_ptr)->considered_peaks_count = 4;
+		(*opts_ptr)->angle_resolution = 2;
+		(*opts_ptr)->refinement_type = 1;
+		(*opts_ptr)->tolerance = 0.06;
+		(*opts_ptr)->maxResolutionForIndexing_1_per_A = +INFINITY;
+		(*opts_ptr)->thread_count = 1;
+		(*opts_ptr)->multi = 0;
+		(*opts_ptr)->no_check_indexed = 0;
+		(*opts_ptr)->min_peaks = 2;
+		(*opts_ptr)->reflectionRadius = -1;
+		break;
+
+		case 1 :
+		show_help();
+		return EINVAL;
+
+		case 2 :
+		if (sscanf(arg, "%u", &(*opts_ptr)->considered_peaks_count) != 1)
+		{
+			ERROR("Invalid value for "
+			      "--pinkIndexer-considered-peaks-count\n");
+			return EINVAL;
+		}
+		break;
+
+		case 3 :
+		if (sscanf(arg, "%u", &(*opts_ptr)->angle_resolution) != 1)
+		{
+			ERROR("Invalid value for "
+			      "--pinkIndexer-angle_resolution\n");
+			return EINVAL;
+		}
+		break;
+
+		case 4 :
+		if (sscanf(arg, "%u", &(*opts_ptr)->refinement_type) != 1)
+		{
+			ERROR("Invalid value for "
+			      "--pinkIndexer-refinement-type\n");
+			return EINVAL;
+		}
+		break;
+
+		case 5 :
+		if (sscanf(arg, "%d", &(*opts_ptr)->thread_count) != 1)
+		{
+			ERROR("Invalid value for --pinkIndexer-thread-count\n");
+			return EINVAL;
+		}
+		break;
+
+		case 6 :
+		if (sscanf(arg, "%f", &(*opts_ptr)->maxResolutionForIndexing_1_per_A) != 1)
+		{
+			ERROR("Invalid value for "
+			      "--pinkIndexer-max-resolution-for-indexing\n");
+			return EINVAL;
+		}
+		break;
+
+		case 7 :
+		if (sscanf(arg, "%f", &(*opts_ptr)->tolerance) != 1)
+		{
+			ERROR("Invalid value for --pinkIndexer-tolerance\n");
+			return EINVAL;
+		}
+		break;
+
+		case 8 :
+		(*opts_ptr)->multi = 1;
+		break;
+
+		case 9 :
+		(*opts_ptr)->no_check_indexed = 1;
+		break;
+
+		case 10 :
+		if (sscanf(arg, "%f", &tmp) != 1) {
+			ERROR("Invalid value for --pinkIndexer-reflection-radius\n");
+			return EINVAL;
+		}
+		(*opts_ptr)->reflectionRadius = tmp / 1e10; /* A^-1 to m^-1 */
+		break;
+
+	}
+
+	return 0;
+}
+
+
+static struct argp_option options[] = {
+
+	{"help-pinkindexer", 1, NULL, OPTION_NO_USAGE,
+	 "Show options for PinkIndexer indexing algorithm", 99},
+
+	{"pinkIndexer-considered-peaks-count", 2, "n", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-cpc", 2, "n", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-angle-resolution", 3, "ang", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-ar", 3, "ang", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-refinement-type", 4, "t", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-rt", 4, "t", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-thread-count", 5, "n", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-tc", 5, "n", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-max-resolution-for-indexing", 6, "res", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-mrfi", 6, "res", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-tolerance", 7, "tol", OPTION_HIDDEN, NULL},
+	{"pinkIndexer-tol", 7, "tol", OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-multi", 8, NULL, OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-no-check-indexed", 9, NULL, OPTION_HIDDEN, NULL},
+
+	{"pinkIndexer-reflection-radius", 10, "r", OPTION_HIDDEN, NULL},
+
+	{0}
+};
+
+
+struct argp pinkIndexer_argp = { options, parse_arg, NULL, NULL, NULL, NULL, NULL };
diff --git a/libcrystfel/src/pinkindexer.h b/libcrystfel/src/pinkindexer.h
new file mode 100644
index 00000000..0169d028
--- /dev/null
+++ b/libcrystfel/src/pinkindexer.h
@@ -0,0 +1,65 @@
+/*
+ * pinkindexer.h
+ *
+ * Interface to PinkIndexer
+ *
+ * Copyright © 2017-2019 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017-2019 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_PINKINDEXER_H_
+#define LIBCRYSTFEL_SRC_PINKINDEXER_H_
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+typedef struct pinkIndexer_options PinkIndexerOptions;
+extern struct argp pinkIndexer_argp;
+
+struct pinkIndexer_options {
+	unsigned int considered_peaks_count;
+	unsigned int angle_resolution;
+	unsigned int refinement_type;
+	float maxResolutionForIndexing_1_per_A;
+	float tolerance;
+	int multi;
+	int thread_count;
+	int min_peaks;
+	int no_check_indexed;
+	float reflectionRadius; /* In m^-1 */
+};
+
+#include <stddef.h>
+#include "index.h"
+
+extern int run_pinkIndexer(struct image *image, void *ipriv);
+
+extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+                                 struct pinkIndexer_options *pinkIndexer_opts,
+                                 struct detector *det, struct beam_params *beam);
+
+extern void pinkIndexer_cleanup(void *pp);
+
+extern const char *pinkIndexer_probe(UnitCell *cell);
+
+#endif /* LIBCRYSTFEL_SRC_PINKINDEXER_H_ */