Add pinkIndexer interface

9 files changed, 613 insertions, 3 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 5b51566f..3d96cf0b 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -90,6 +90,11 @@ Use the TakeTwo algorithm.  See Ginn et al., Acta Cryst. (2016). D72, 956-965.
 .PD
 Invoke XGANDALF - eXtended GrAdieNt Descent Algorithm for Lattice Finding. Xgandalf must be installed in order to be able to use it.
 
+.IP \fBpinkIndexer\fR
+.PD
+Invoke pinkIndexer. pinkIndexer must be installed in order to be able to use it. The geometry file should contain the value photon_energy_bandwidth, which sets the bandwidth as (delta energy)/(mean energy).
+
+
 .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.
 
@@ -437,6 +442,38 @@ These set low-level parameters for the XGANDALF indexing algorithm.
 .IP
 \fB--xgandalf-fast-execution\fR Shortcut to set --xgandalf-sampling-pitch=2 --xgandalf-grad-desc-iterations=3
 
+.PD 0
+.IP \fB--pinkIndexer-considered-peaks-count=\fIn\fR
+.IP \fB--pinkIndexer-angle-resolution=\fIn\fR
+.IP \fB--pinkIndexer-refinement-type=\fIn\fR
+.IP \fB--pinkIndexer-tolerance=\fIn\fR
+.IP \fB--pinkIndexer-reflection-radius=\fIn\fR
+.IP \fB--pinkIndexer-max-resolution-for-indexing=\fIn\fR
+.IP \fB--pinkIndexer-multi=\fIn\fR
+.IP \fB--pinkIndexer-thread-count=\fIn\fR
+.IP \fB--pinkIndexer-no-check-indexed=\fIn\fR
+
+.PD
+These set low-level parameters for the XGANDALF indexing algorithm.
+.IP
+\fB--pinkIndexer-considered-peaks-count\fR selects how many peaks are considered for indexing. [0-4] (veryFew to manyMany). Default is 4 (manyMany).
+.IP
+\fB--pinkIndexer-angle-resolution\fR selects how dense the orientation angles of the sample lattice are sampled. [0-4] (extremelyLoose to extremelyDense). Default is 2 (normal).
+.IP
+\fB--pinkIndexer-refinement-type\fR selects the refinement type. 0 = none, 1 = fixedLatticeParameters, 2 = variableLatticeParameters, 3 = firstFixedThenVariableLatticeParameters, 4 = firstFixedThenVariableLatticeParametersMultiSeed, 5 = firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed.
+.IP
+\fB--pinkIndexer-tolerance\fR selects the tolerance of the pinkIndexer (relative tolerance of the lattice vectors). Default is 0.06. For bad geometrys or cell parameters use a high tolerance. For a well known geometry and cell use a small tolerance. Only important for refinement and indexed/not indexed identificaton. Too small tolerance will lead to refining to only a fraction of the peaks and possibly discarding of correctly indexed images. Too high tolerance will lead to bad fitting in presence of multiples or noise and can mark wrongly-indexed patterns as indexed.
+.IP
+\fB--pinkIndexer-reflection-radius\fR sets radius of the reflections in reciprocal space in 1/A. Default is 2%% of a* (which works quiet well for X-rays). Should be chosen much bigger for electrons (~0.002).
+.IP
+\fB--pinkIndexer-max-resolution-for-indexing\fR sets the maximum resolition in 1/A used for indexing. Peaks at high resolution don't add much information, but they add a lot of computation time. Default is infinity. Does not influence the refinement.
+.IP
+\fB--pinkIndexer-multi\fR Use pinkIndexers own multi indexing. Should be combined with the --no-multi flag.
+.IP
+\fB--pinkIndexer-thread-count\fR sets the thread count for internal parallelization. Default is 1. Very useful for small datasets (e.g. for screening). Internal parallelization does not significantly increase the amount of RAM needed, whereas CrystFELs parallelization does. For HPCs typically a mixture of both parallelizations leads to best results.
+.IP
+\fB--pinkIndexer-no-check-indexed\fR Leave the check whether a pattern is indexed completely to CrystFEL. Useful for monochromatic (since CrystFELs prediction model is smarter than the one of pinkIndexer) or in combnation with --no-check-peaks for geometry optimization. This flag is meant to eventually disappear, when the full pink pipeline is implemented.
+
 .SH INTEGRATION OPTIONS
 .PD 0
 .IP \fB--integration=\fR\fImethod\fR
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index 67a20ee5..fdba725e 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -14,6 +14,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
@@ -67,6 +68,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}
@@ -108,6 +110,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/index.c b/libcrystfel/src/index.c
index 4c82e868..cc004d9a 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;
@@ -233,6 +239,7 @@ static char *friendly_indexer_name(IndexingMethod m)
 
 static void *prepare_method(IndexingMethod *m, UnitCell *cell,
                             struct xgandalf_options *xgandalf_opts,
+                            struct pinkIndexer_options* pinkIndexer_opts,
                             struct felix_options *felix_opts)
 {
 	char *str;
@@ -277,6 +284,10 @@ 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);
+		break;
+
 		default :
 		ERROR("Don't know how to prepare indexing method %i\n", *m);
 		break;
@@ -309,6 +320,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
                                 IndexingFlags flags,
                                 struct taketwo_options *ttopts,
                                 struct xgandalf_options *xgandalf_opts,
+                                struct pinkIndexer_options *pinkIndexer_opts,
                                 struct felix_options *felix_opts)
 {
 	int i, n;
@@ -403,7 +415,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 		int j;
 
 		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
-		                                          xgandalf_opts,
+		                                          xgandalf_opts, pinkIndexer_opts,
 		                                          felix_opts);
 
 		if ( ipriv->engine_private[i] == NULL ) return NULL;
@@ -430,6 +442,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++ ) {
@@ -490,6 +503,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]);
@@ -601,6 +618,10 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		r = taketwo_index(image, ipriv->ttopts, mpriv);
 		break;
 
+		case INDEXING_PINKINDEXER :
+		r = run_pinkIndexer(image, mpriv);
+		break;
+
 		case INDEXING_XGANDALF :
 		set_last_task(last_task, "indexing:xgandalf");
 		r = run_xgandalf(image, mpriv);
@@ -999,6 +1020,11 @@ 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) {
+			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;
@@ -1097,9 +1123,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 ec7ed61b..ab38016f 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,14 +149,15 @@ 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,
                                        IndexingFlags flags,
                                        struct taketwo_options *ttopts,
                                        struct xgandalf_options *xgandalf_opts,
+                                       struct pinkIndexer_options *pinkIndexer_opts,
                                        struct felix_options *felix_opts);
 
 extern char *detect_indexing_methods(UnitCell *cell);
diff --git a/libcrystfel/src/pinkindexer.c b/libcrystfel/src/pinkindexer.c
new file mode 100644
index 00000000..2002533f
--- /dev/null
+++ b/libcrystfel/src/pinkindexer.c
@@ -0,0 +1,361 @@
+/*
+ * pinkIndexer.c
+ *
+ *  Created on: Nov 27, 2017
+ *      Author: gevorkov
+ */
+
+#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;
+		printf("\ntotal peaks count %i\n", peakCount);
+		int matchedPeaksCount = PinkIndexer_indexPattern(pinkIndexer_private_data->pinkIndexer,
+		        &(indexedLattice[indexed]), center_shift[indexed], reciprocalPeaks_1_per_A, intensities,
+		        maxRefinementDisbalance,
+		        pinkIndexer_private_data->threadCount);
+
+		printf("center shift = %f %f", center_shift[indexed][0], center_shift[indexed][1]);
+		printf("\nmatchedPeaksCount %i from %i\n", matchedPeaksCount, peakCount);
+
+		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);
+			printf("before restoration\n");
+			cell_print(uc);
+
+			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;
+			printf("after restoration\n");
+			cell_print(uc);
+
+			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)) {
+				printf("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++;
+
+			printf("crystal %i in image added\n", indexed);
+		}
+		else {
+			break;
+		}
+	} while (pinkIndexer_private_data->multi
+	        && indexed <= MAX_MULTI_LATTICE_COUNT
+	        && reciprocalPeaks_1_per_A->peakCount >= pinkIndexer_private_data->min_peaks);
+
+	printf("\n\nfound %i crystals in image\n", indexed);
+	return indexed;
+}
+
+void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
+        struct pinkIndexer_options *pinkIndexer_opts)
+{
+	printf("preparing pink indexer\n");
+
+	if (pinkIndexer_opts->beamEnergy == 0.0 || pinkIndexer_opts->detectorDistance <= 0) {
+		ERROR("ERROR!!!!!! photon_energy and photon_energy_bandwidth must be defined as constants "
+				"in the geometry file for the pinkIndexer!!!!!!");
+	}
+	if (pinkIndexer_opts->beamBandwidth == 0.0) {
+		printf("using default bandwidth of 0.01 for pinkIndexer!");
+		pinkIndexer_opts->beamBandwidth = 0.01;
+	}
+	if (pinkIndexer_opts->detectorDistance == 0.0 && 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);
+
+	printf("reduced cell:\n");
+	cell_print(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 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 = pinkIndexer_opts->detectorDistance;
+	float beamEenergy_eV = pinkIndexer_opts->beamEnergy;
+	float nonMonochromaticity = pinkIndexer_opts->beamBandwidth;
+	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)
+{
+	printf("cleaning up pink indexer\n");
+
+	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)
+{
+	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 */
diff --git a/libcrystfel/src/pinkindexer.h b/libcrystfel/src/pinkindexer.h
new file mode 100644
index 00000000..cd2aaba2
--- /dev/null
+++ b/libcrystfel/src/pinkindexer.h
@@ -0,0 +1,46 @@
+/*
+ * pinkIndexer.h
+ *
+ *  Created on: Nov 27, 2017
+ *      Author: gevorkov
+ */
+
+#ifndef LIBCRYSTFEL_SRC_PINKINDEXER_H_
+#define LIBCRYSTFEL_SRC_PINKINDEXER_H_
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+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 beamEnergy; //in eV
+	float beamBandwidth; //(delta lambda)/lambda
+	float detectorDistance; //in m
+
+	float reflectionRadius; //in 1/A
+};
+
+#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);
+
+extern void pinkIndexer_cleanup(void *pp);
+
+extern const char *pinkIndexer_probe(UnitCell *cell);
+
+#endif /* LIBCRYSTFEL_SRC_PINKINDEXER_H_ */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index bd159a9a..5770051a 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -204,6 +204,32 @@ static void show_help(const char *s)
 "                           All peaks are used for refinement.\n"
 "                            Default: 250\n"
 "\n"
+"     --pinkIndexer-considered-peaks-count  Considered peaks count selector \n"
+"                                           [0-4] (veryFew to manyMany)\n"
+"                                           Default is 4 (manyMany)\n"
+"     --pinkIndexer-angle-resolution     Angle resolution selector \n"
+"                                        [0-4] (extremelyLoose to extremelyDense)\n"
+"                                        Default is 2 (normal)\n"
+"     --pinkIndexer-refinement-type    Refinement type \n"
+"                                      Default is 1\n"
+"               0 = none\n"
+"               1 = fixedLatticeParameters\n"
+"               2 = variableLatticeParameters\n"
+"               3 = firstFixedThenVariableLatticeParameters\n"
+"               4 = firstFixedThenVariableLatticeParametersMultiSeed\n"
+"               5 = firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed\n"
+"     --pinkIndexer-tolerance  Relative tolerance of the lattice vectors.\n"
+"                              Default is 0.06\n"
+"     --pinkIndexer-reflection-radius  radius of the reflections in reciprocal space in 1/A.\n"
+"                              Default is 2%% of a*.\n"
+"                              Should be chosen ~0.002 for electrons.\n"
+"     --pinkIndexer-max-resolution-for-indexing   Measured in 1/A\n"
+"     --pinkIndexer-multi    Use pinkIndexers own multi indexing.\n"
+"     --pinkIndexer-thread-count  thread count for internal parallelization \n"
+"                                 Default is 1\n"
+"     --pinkIndexer-no-check-indexed   Leave the check whether a pattern is indexed completely to CrystFEL\n"
+"									   useful for monochromatic and tweaking."
+"\n"
 "\nIntegration options:\n\n"
 "     --integration=<meth>  Integration method (rings,prof2d)-(cen,nocen)\n"
 "                            Default: rings-nocen\n"
@@ -253,6 +279,7 @@ static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copym
 int main(int argc, char *argv[])
 {
 	int c;
+	unsigned int tmp_enum;
 	char *filename = NULL;
 	char *outfile = NULL;
 	FILE *fh;
@@ -360,6 +387,16 @@ int main(int argc, char *argv[])
 	iargs.xgandalf_opts.minLatticeVectorLength_A = 30;
 	iargs.xgandalf_opts.maxLatticeVectorLength_A = 250;
 	iargs.xgandalf_opts.maxPeaksForIndexing = 250;
+	iargs.pinkIndexer_opts.considered_peaks_count = 4;
+	iargs.pinkIndexer_opts.angle_resolution = 2;
+	iargs.pinkIndexer_opts.refinement_type = 1;
+	iargs.pinkIndexer_opts.tolerance = 0.06;
+	iargs.pinkIndexer_opts.maxResolutionForIndexing_1_per_A = FLT_MAX;
+	iargs.pinkIndexer_opts.thread_count = 1;
+	iargs.pinkIndexer_opts.multi = 0;
+	iargs.pinkIndexer_opts.no_check_indexed = 0;
+	iargs.pinkIndexer_opts.min_peaks = 2;
+	iargs.pinkIndexer_opts.reflectionRadius = -1;
 	iargs.felix_opts.ttmin = -1.0;
 	iargs.felix_opts.ttmax = -1.0;
 	iargs.felix_opts.min_visits = 0;
@@ -491,6 +528,21 @@ int main(int argc, char *argv[])
 		{"min-sq-gradient",    1, NULL,              359}, /* compat */
 		{"xgandalf-fast-execution",                  0, NULL, 360},
 		{"xgandalf-max-peaks",                       1, NULL, 361},
+		{"pinkIndexer-considered-peaks-count",       1, NULL, 362},
+		{"pinkIndexer-cpc",                          1, NULL, 362},
+		{"pinkIndexer-angle-resolution",             1, NULL, 363},
+		{"pinkIndexer-ar",                           1, NULL, 363},
+		{"pinkIndexer-refinement-type",              1, NULL, 364},
+		{"pinkIndexer-rt",                           1, NULL, 364},
+		{"pinkIndexer-thread-count",                 1, NULL, 365},
+		{"pinkIndexer-tc",                           1, NULL, 365},
+		{"pinkIndexer-max-resolution-for-indexing",  1, NULL, 366},
+		{"pinkIndexer-mrfi",                         1, NULL, 366},
+		{"pinkIndexer-tolerance",                    1, NULL, 367},
+		{"pinkIndexer-tol",                          1, NULL, 367},
+		{"pinkIndexer-multi",                        0, NULL, 368},
+		{"pinkIndexer-no-check-indexed",             0, NULL, 369},
+		{"pinkIndexer-reflection-radius",            1, NULL, 370},
 
 		{0, 0, NULL, 0}
 	};
@@ -689,6 +741,7 @@ int main(int argc, char *argv[])
 
 			case 331:
 			iargs.min_peaks = atoi(optarg);
+			iargs.pinkIndexer_opts.min_peaks = iargs.min_peaks;
 			break;
 
 			case 332:
@@ -922,6 +975,74 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 362:
+			if (sscanf(optarg, "%u", &tmp_enum) != 1)
+			{
+				ERROR("Invalid value for "
+						"--pinkIndexer-considered-peaks-count\n");
+				return 1;
+			}
+			iargs.pinkIndexer_opts.considered_peaks_count = tmp_enum;
+			break;
+
+			case 363:
+			if (sscanf(optarg, "%u", &tmp_enum) != 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-angle-resolution \n");
+				return 1;
+			}
+			iargs.pinkIndexer_opts.angle_resolution = tmp_enum;
+			break;
+
+			case 364:
+			if (sscanf(optarg, "%u", &tmp_enum) != 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-refinement-type \n");
+				return 1;
+			}
+			iargs.pinkIndexer_opts.refinement_type = tmp_enum;
+			break;
+
+			case 365:
+			if (sscanf(optarg, "%d", &iargs.pinkIndexer_opts.thread_count) != 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-thread-count \n");
+				return 1;
+			}
+			break;
+
+			case 366:
+			if (sscanf(optarg, "%f", &iargs.pinkIndexer_opts.maxResolutionForIndexing_1_per_A) != 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-max-resolution-for-indexing \n");
+				return 1;
+			}
+			break;
+
+			case 367:
+			if (sscanf(optarg, "%f", &iargs.pinkIndexer_opts.tolerance)	!= 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-tolerance \n");
+				return 1;
+			}
+			break;
+
+			case 368:
+				iargs.pinkIndexer_opts.multi = 1;
+			break;
+
+			case 369:
+				iargs.pinkIndexer_opts.no_check_indexed = 1;
+			break;
+
+			case 370:
+			if (sscanf(optarg, "%f", &iargs.pinkIndexer_opts.reflectionRadius)	!= 1)
+			{
+				ERROR("Invalid value for --pinkIndexer-reflection-radius \n");
+				return 1;
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -1021,6 +1142,12 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 	add_geom_beam_stuff_to_field_list(iargs.copyme, iargs.det, iargs.beam);
+	iargs.pinkIndexer_opts.beamEnergy = iargs.beam->photon_energy;
+	iargs.pinkIndexer_opts.beamBandwidth = iargs.beam->photon_energy_bandwidth;
+	iargs.pinkIndexer_opts.detectorDistance = iargs.det->panels[0].clen;
+	if(iargs.det->panels[0].clen_from != NULL){
+		iargs.pinkIndexer_opts.detectorDistance = 0;
+	}
 
 	/* If no peak path from geometry file, use these (but see later) */
 	if ( iargs.hdf5_peak_path == NULL ) {
@@ -1263,6 +1390,7 @@ int main(int argc, char *argv[])
 		                             iargs.tols, flags,
 		                             &iargs.taketwo_opts,
 		                             &iargs.xgandalf_opts,
+		                             &iargs.pinkIndexer_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 9a58a64f..ced3bcf2 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -46,6 +46,7 @@ struct index_args;
 #include "time-accounts.h"
 #include "taketwo.h"
 #include "xgandalf.h"
+#include "../libcrystfel/src/pinkindexer.h"
 #include "felix.h"
 
 
@@ -117,6 +118,7 @@ struct index_args
 	int profile;  /* Whether or not to do wall clock profiling */
 	struct taketwo_options taketwo_opts;
 	struct xgandalf_options xgandalf_opts;
+    struct pinkIndexer_options pinkIndexer_opts;
 	struct felix_options felix_opts;
 	Spectrum *spectrum;
 	signed int wait_for_file; /* -1 means wait forever */