Merge branch 'taketwo'

10 files changed, 1838 insertions, 7 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index eaf812d6..c3189f22 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -402,6 +402,14 @@ Skip the prediction refinement step.
 Display timing data for performance monitoring.
 
 .PD 0
+.IP \fB--taketwo-member-threshold=\fIn\fR
+.IP \fB--taketwo-len-tolerance=\fIn\fR
+.IP \fB--taketwo-angle-tolerance=\fIn\fR
+.IP \fB--taketwo-trace-tolerance=\fIn\fR
+.PD
+These set low-level parameters for the TakeTwo indexing algorithm.  Respectively, the minimum number of vectors in the network before the pattern is considered indexed, the length and angle tolerances (in reciprocal Angstroms and degrees, respectively) and the matrix trace tolerance for considering rotation matrices as equal.
+
+.PD 0
 .IP \fB--felix-options\fR
 .PD
 Specify a comma-separated list of keyword arguments to change the default parameters passed to the felix indexer. For asthetics, this list can also be preceeded and followed by quotations - as in: \fB--felix-options="arg1=10,arg2=600"\fR. A list of parameters which can be modified through this option is detailed below.
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 8af82552..fef2e11f 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -11,7 +11,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
                          src/cell-utils.c src/integer_matrix.c src/crystal.c \
                          src/xds.c src/integration.c src/predict-refine.c \
                          src/histogram.c src/events.c src/felix.c \
-                         src/peakfinder8.c
+                         src/taketwo.c src/peakfinder8.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/asdf.c
@@ -32,7 +32,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \
                                  src/xds.h src/predict-refine.h \
                                  src/integration.h src/histogram.h \
                                  src/events.h src/asdf.h src/felix.h \
-                                 src/peakfinder8.h
+                                 src/taketwo.h src/peakfinder8.h
 
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
 AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 84a14b5c..6610abc3 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -38,7 +38,6 @@
 #include <string.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
 #include <assert.h>
 
 #include "cell.h"
@@ -1438,6 +1437,49 @@ void cell_fudge_gslcblas()
 }
 
 
+UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m)
+{
+	gsl_matrix *c;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	gsl_matrix *res;
+	UnitCell *out;
+
+	cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
+	                        &bsz, &csx, &csy, &csz);
+
+	c = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(c, 0, 0, asx);
+	gsl_matrix_set(c, 1, 0, asy);
+	gsl_matrix_set(c, 2, 0, asz);
+	gsl_matrix_set(c, 0, 1, bsx);
+	gsl_matrix_set(c, 1, 1, bsy);
+	gsl_matrix_set(c, 2, 1, bsz);
+	gsl_matrix_set(c, 0, 2, csx);
+	gsl_matrix_set(c, 1, 2, csy);
+	gsl_matrix_set(c, 2, 2, csz);
+
+	res = gsl_matrix_calloc(3, 3);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
+
+	out = cell_new_from_cell(in);
+	cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0),
+	                         gsl_matrix_get(res, 1, 0),
+	                         gsl_matrix_get(res, 2, 0),
+	                         gsl_matrix_get(res, 0, 1),
+	                         gsl_matrix_get(res, 1, 1),
+	                         gsl_matrix_get(res, 2, 1),
+	                         gsl_matrix_get(res, 0, 2),
+	                         gsl_matrix_get(res, 1, 2),
+	                         gsl_matrix_get(res, 2, 2));
+
+	gsl_matrix_free(res);
+	gsl_matrix_free(c);
+	return out;
+}
+
+
 /**
  * rotate_cell:
  * @in: A %UnitCell to rotate
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 000c863a..efb4b25b 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -35,6 +35,8 @@
 #include <config.h>
 #endif
 
+#include <gsl/gsl_matrix.h>
+
 #include "cell.h"
 
 #ifdef __cplusplus
@@ -47,6 +49,7 @@ extern double resolution(UnitCell *cell,
 extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
 extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
                              double rot);
+extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m);
 
 extern void cell_print(UnitCell *cell);
 
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 446f57a3..b00f13aa 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -56,6 +56,7 @@
 #include "cell-utils.h"
 #include "felix.h"
 #include "predict-refine.h"
+#include "taketwo.h"
 
 
 struct _indexingprivate
@@ -63,6 +64,8 @@ struct _indexingprivate
 	UnitCell *target_cell;
 	float tolerance[4];
 
+	struct taketwo_options *ttopts;
+
 	int n_methods;
 	IndexingMethod *methods;
 	void **engine_private;
@@ -120,6 +123,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 		priv = felix_prepare(m, cell, det, ltl, options);
 		break;
 
+		case INDEXING_TAKETWO :
+		priv = taketwo_prepare(m, cell, det, ltl);
+		break;
+
 		default :
 		ERROR("Don't know how to prepare indexing method %i\n", *m);
 		break;
@@ -152,7 +159,8 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 
 IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
                                 struct detector *det, float *ltl,
-                                int no_refine, const char *options)
+                                int no_refine, const char *options,
+                                struct taketwo_options *ttopts)
 {
 	int i, n;
 	char **method_strings;
@@ -211,6 +219,8 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 	}
 	for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i];
 
+	ipriv->ttopts = ttopts;
+
 	return ipriv;
 }
 
@@ -252,6 +262,10 @@ void cleanup_indexing(IndexingPrivate *ipriv)
 			free(ipriv->engine_private[n]);
 			break;
 
+			case INDEXING_TAKETWO :
+			taketwo_cleanup(ipriv->engine_private[n]);
+			break;
+
 			default :
 			ERROR("Don't know how to clean up indexing method %i\n",
 			      ipriv->methods[n]);
@@ -329,6 +343,10 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		r = felix_index(image, mpriv);
 		break;
 
+		case INDEXING_TAKETWO :
+		r = taketwo_index(image, ipriv->ttopts, mpriv);
+		break;
+
 		default :
 		ERROR("Unrecognised indexing method: %i\n", indm);
 		return 0;
@@ -662,6 +680,10 @@ char *indexer_str(IndexingMethod indm)
 		strcpy(str, "xds");
 		break;
 
+		case INDEXING_TAKETWO :
+		strcpy(str, "taketwo");
+		break;
+
 		case INDEXING_SIMULATION :
 		strcpy(str, "simulation");
 		break;
@@ -770,6 +792,11 @@ IndexingMethod get_indm_from_string(const char *str)
 			method = INDEXING_DEFAULTS_FELIX;
 			have_method = 1;
 
+		} else if ( strcmp(bits[i], "taketwo") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_TAKETWO;
+			have_method = 1;
+
 		} else if ( strcmp(bits[i], "none") == 0) {
 			if ( have_method ) return warn_method(str);
 			method = INDEXING_NONE;
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 8dd7de21..8b3aadeb 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -58,6 +58,11 @@
                                      | INDEXING_USE_CELL_PARAMETERS            \
                                      | INDEXING_RETRY | INDEXING_REFINE)
 
+#define INDEXING_DEFAULTS_TAKETWO (INDEXING_TAKETWO | INDEXING_CHECK_PEAKS     \
+                                     | INDEXING_USE_CELL_PARAMETERS            \
+                                     | INDEXING_USE_LATTICE_TYPE               \
+                                     | INDEXING_RETRY | INDEXING_REFINE)
+
 /* Axis check is needed for XDS, because it likes to permute the axes */
 #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE        \
                                      | INDEXING_USE_CELL_PARAMETERS            \
@@ -75,6 +80,7 @@
  * @INDEXING_SIMULATION: Dummy value
  * @INDEXING_DEBUG: Results injector for debugging
  * @INDEXING_ASDF: Use in-built "asdf" indexer
+ * @INDEXING_TAKETWO: Use in-built "taketwo" indexer
  * @INDEXING_ERROR: Special value for unrecognised indexing engine name
  * @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell
  *   axes for agreement with given cell.
@@ -108,6 +114,7 @@ typedef enum {
 	INDEXING_SIMULATION = 6,
 	INDEXING_DEBUG = 7,
 	INDEXING_ASDF = 8,
+	INDEXING_TAKETWO = 9,
 
 	INDEXING_ERROR = 255,  /* Unrecognised indexing engine */
 
@@ -150,10 +157,13 @@ extern IndexingMethod get_indm_from_string(const char *method);
 #include "detector.h"
 #include "cell.h"
 #include "image.h"
+#include "taketwo.h"
+
 
 extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell,
                                        struct detector *det, float *ltl,
-                                       int no_refine, const char *options);
+                                       int no_refine, const char *options,
+                                       struct taketwo_options *ttopts);
 
 extern void index_pattern(struct image *image, IndexingPrivate *ipriv);
 
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
new file mode 100644
index 00000000..ecd2c447
--- /dev/null
+++ b/libcrystfel/src/taketwo.c
@@ -0,0 +1,1638 @@
+/*
+ * taketwo.c
+ *
+ * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL
+ *
+ * Copyright © 2016 Helen Ginn
+ * Copyright © 2016 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2016 Helen Ginn <helen@strubi.ox.ac.uk>
+ *   2016 Thomas White <taw@physics.org>
+ *
+ * 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 <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <float.h>
+#include <math.h>
+#include <assert.h>
+#include <time.h>
+
+#include "cell-utils.h"
+#include "index.h"
+#include "taketwo.h"
+#include "peaks.h"
+#include "symmetry.h"
+
+/**
+ * spotvec
+ * @obsvec: an observed vector between two spots
+ * @matches: array of matching theoretical vectors from unit cell
+ * @match_num: number of matches
+ * @distance: length of obsvec (do I need this?)
+ * @her_rlp: pointer to first rlp position for difference vec
+ * @his_rlp: pointer to second rlp position for difference vec
+ *
+ * Structure representing 3D vector between two potential Bragg peaks
+ * in reciprocal space, and an array of potential matching theoretical
+ * vectors from unit cell/centering considerations.
+ **/
+struct SpotVec
+{
+	struct rvec obsvec;
+	struct rvec *matches;
+	int match_num;
+	struct rvec *asym_matches;
+	int asym_match_num;
+	double distance;
+	struct rvec *her_rlp;
+	struct rvec *his_rlp;
+};
+
+
+struct taketwo_private
+{
+	IndexingMethod indm;
+	float          *ltl;
+	UnitCell       *cell;
+};
+
+// These rotation symmetry operators
+struct TakeTwoCell
+{
+	UnitCell       *cell;
+	gsl_matrix     **rotSymOps;
+	unsigned int   numOps;
+	struct SpotVec **obs_vecs; // Pointer to an array
+	int obs_vec_count;
+	struct taketwo_options *options;
+	gsl_matrix *twiz1Tmp;
+	gsl_matrix *twiz2Tmp;
+	gsl_vector *vec1Tmp;
+	gsl_vector *vec2Tmp;
+};
+
+
+/* Maximum distance between two rlp sizes to consider info for indexing */
+#define MAX_RECIP_DISTANCE (0.15*1e10)
+
+/* Tolerance for two lengths in reciprocal space to be considered the same */
+#define RECIP_TOLERANCE (0.0010*1e10)
+
+/* Threshold for network members to consider a potential solution */
+#define NETWORK_MEMBER_THRESHOLD (20)
+
+/* Maximum network members (obviously a solution so should stop) */
+#define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 3)
+
+/* Maximum dead ends for a single branch extension during indexing */
+#define MAX_DEAD_ENDS (10)
+
+/* Maximum observed vectors before TakeTwo gives up and deals with
+ * what is already there. */
+#define MAX_OBS_VECTORS 100000
+
+/* Tolerance for two angles to be considered the same */
+#define ANGLE_TOLERANCE (deg2rad(0.6))
+
+/* Tolerance for rot_mats_are_similar */
+#define TRACE_TOLERANCE (deg2rad(3.0))
+
+/* TODO: Multiple lattices */
+
+
+/* ------------------------------------------------------------------------
+ * apologetic function
+ * ------------------------------------------------------------------------*/
+
+void apologise()
+{
+	printf("Error - could not allocate memory for indexing.\n");
+}
+
+/* ------------------------------------------------------------------------
+ * functions concerning aspects of rvec which are very likely to be
+ * incorporated somewhere else in CrystFEL and therefore may need to be
+ * deleted and references connected to a pre-existing function. (Lowest level)
+ * ------------------------------------------------------------------------*/
+
+static struct rvec new_rvec(double new_u, double new_v, double new_w)
+{
+	struct rvec new_rvector;
+	new_rvector.u = new_u;
+	new_rvector.v = new_v;
+	new_rvector.w = new_w;
+
+	return new_rvector;
+}
+
+
+static struct rvec diff_vec(struct rvec from, struct rvec to)
+{
+	struct rvec diff = new_rvec(to.u - from.u,
+				    to.v - from.v,
+				    to.w - from.w);
+
+	return diff;
+}
+
+static double sq_length(struct rvec vec)
+{
+	double sqlength = (vec.u * vec.u + vec.v * vec.v + vec.w * vec.w);
+
+	return sqlength;
+}
+
+
+static double rvec_length(struct rvec vec)
+{
+	return sqrt(sq_length(vec));
+}
+
+
+static void normalise_rvec(struct rvec *vec)
+{
+	double length = rvec_length(*vec);
+	vec->u /= length;
+	vec->v /= length;
+	vec->w /= length;
+}
+
+
+static double rvec_cosine(struct rvec v1, struct rvec v2)
+{
+	double dot_prod = v1.u * v2.u + v1.v * v2.v + v1.w * v2.w;
+	double v1_length = rvec_length(v1);
+	double v2_length = rvec_length(v2);
+
+	double cos_theta = dot_prod / (v1_length * v2_length);
+
+	return cos_theta;
+}
+
+
+static double rvec_angle(struct rvec v1, struct rvec v2)
+{
+	double cos_theta = rvec_cosine(v1, v2);
+	double angle = acos(cos_theta);
+
+	return angle;
+}
+
+
+static struct rvec rvec_cross(struct rvec a, struct rvec b)
+{
+	struct rvec c;
+
+	c.u = a.v*b.w - a.w*b.v;
+	c.v = -(a.u*b.w - a.w*b.u);
+	c.w = a.u*b.v - a.v*b.u;
+
+	return c;
+}
+
+/*
+static void show_rvec(struct rvec r2)
+{
+	struct rvec r = r2;
+	normalise_rvec(&r);
+	STATUS("[ %.3f %.3f %.3f ]\n", r.u, r.v, r.w);
+}
+*/
+
+
+/* ------------------------------------------------------------------------
+ * functions called under the core functions, still specialised (Level 3)
+ * ------------------------------------------------------------------------*/
+
+static void rotation_around_axis(struct rvec c, double th,
+				 gsl_matrix *res)
+{
+	double omc = 1.0 - cos(th);
+	double s = sin(th);
+	gsl_matrix_set(res, 0, 0, cos(th) + c.u*c.u*omc);
+	gsl_matrix_set(res, 0, 1, c.u*c.v*omc - c.w*s);
+	gsl_matrix_set(res, 0, 2, c.u*c.w*omc + c.v*s);
+	gsl_matrix_set(res, 1, 0, c.u*c.v*omc + c.w*s);
+	gsl_matrix_set(res, 1, 1, cos(th) + c.v*c.v*omc);
+	gsl_matrix_set(res, 1, 2, c.v*c.w*omc - c.u*s);
+	gsl_matrix_set(res, 2, 0, c.w*c.u*omc - c.v*s);
+	gsl_matrix_set(res, 2, 1, c.w*c.v*omc + c.u*s);
+	gsl_matrix_set(res, 2, 2, cos(th) + c.w*c.w*omc);
+}
+
+
+/* Rotate vector (vec1) around axis (axis) by angle theta. Find value of
+ * theta for which the angle between (vec1) and (vec2) is minimised. */
+static void closest_rot_mat(struct rvec vec1, struct rvec vec2,
+			    struct rvec axis, gsl_matrix *twizzle)
+{
+	/* Let's have unit vectors */
+	normalise_rvec(&vec1);
+	normalise_rvec(&vec2);
+	normalise_rvec(&axis);
+
+	/* Redeclaring these to try and maintain readability and
+	 * check-ability against the maths I wrote down */
+	double a = vec2.u; double b = vec2.v; double c = vec2.w;
+	double p = vec1.u; double q = vec1.v; double r = vec1.w;
+	double x = axis.u; double y = axis.v; double z = axis.w;
+
+	/* Components in handwritten maths online when I upload it */
+	double A = a*(p*x*x - p + x*y*q + x*z*r) +
+	b*(p*x*y + q*y*y - q + r*y*z) +
+	c*(p*x*z + q*y*z + r*z*z - r);
+
+	double B = a*(y*r - z*q) + b*(p*z - r*x) + c*(q*x - p*y);
+
+	double tan_theta = - B / A;
+	double theta = atan(tan_theta);
+
+	/* Now we have two possible solutions, theta or theta+pi
+	 * and we need to work out which one. This could potentially be
+	 * simplified - do we really need so many cos/sins? maybe check
+	 * the 2nd derivative instead? */
+	double cc = cos(theta);
+	double C = 1 - cc;
+	double s = sin(theta);
+	double occ = -cc;
+	double oC = 1 - occ;
+	double os = -s;
+
+	double pPrime = (x*x*C+cc)*p + (x*y*C-z*s)*q + (x*z*C+y*s)*r;
+	double qPrime = (y*x*C+z*s)*p + (y*y*C+cc)*q + (y*z*C-x*s)*r;
+	double rPrime = (z*x*C-y*s)*p + (z*y*C+x*s)*q + (z*z*C+cc)*r;
+
+	double pDbPrime = (x*x*oC+occ)*p + (x*y*oC-z*os)*q + (x*z*oC+y*os)*r;
+	double qDbPrime = (y*x*oC+z*os)*p + (y*y*oC+occ)*q + (y*z*oC-x*os)*r;
+	double rDbPrime = (z*x*oC-y*os)*p + (z*y*oC+x*os)*q + (z*z*oC+occ)*r;
+
+	double cosAlpha = pPrime * a + qPrime * b + rPrime * c;
+	double cosAlphaOther = pDbPrime * a + qDbPrime * b + rDbPrime * c;
+
+	int addPi = (cosAlphaOther > cosAlpha);
+	double bestAngle = theta + addPi * M_PI;
+
+	/* Don't return an identity matrix which has been rotated by
+	 * theta around "axis", but do assign it to twizzle. */
+	rotation_around_axis(axis, bestAngle, twizzle);
+}
+
+
+static double matrix_trace(gsl_matrix *a)
+{
+	int i;
+	double tr = 0.0;
+
+	assert(a->size1 == a->size2);
+	for ( i=0; i<a->size1; i++ ) {
+		tr += gsl_matrix_get(a, i, i);
+	}
+	return tr;
+}
+
+static char *add_ua(const char *inp, char ua)
+{
+	char *pg = malloc(64);
+	if ( pg == NULL ) return NULL;
+	snprintf(pg, 63, "%s_ua%c", inp, ua);
+	return pg;
+}
+
+
+static char *get_chiral_holohedry(UnitCell *cell)
+{
+	LatticeType lattice = cell_get_lattice_type(cell);
+	char *pg = malloc(64);
+	char *pgout = 0;
+
+	if ( pg == NULL ) return NULL;
+
+	switch (lattice)
+	{
+		case L_TRICLINIC:
+		pg = "1";
+		break;
+
+		case L_MONOCLINIC:
+		pg = "2";
+		break;
+
+		case L_ORTHORHOMBIC:
+		pg = "222";
+		break;
+
+		case L_TETRAGONAL:
+		pg = "422";
+		break;
+
+		case L_RHOMBOHEDRAL:
+		pg = "3_R";
+		break;
+
+		case L_HEXAGONAL:
+		if ( cell_get_centering(cell) == 'H' ) {
+			pg = "3_H";
+		} else {
+			pg = "622";
+		}
+		break;
+
+		case L_CUBIC:
+		pg = "432";
+		break;
+
+		default:
+		pg = "error";
+		break;
+	}
+
+	switch (lattice)
+	{
+		case L_TRICLINIC:
+		case L_ORTHORHOMBIC:
+		case L_RHOMBOHEDRAL:
+		case L_CUBIC:
+		pgout = strdup(pg);
+		break;
+
+		case L_MONOCLINIC:
+		case L_TETRAGONAL:
+		case L_HEXAGONAL:
+		pgout = add_ua(pg, cell_get_unique_axis(cell));
+		break;
+
+		default:
+		break;
+	}
+	
+	return pgout;
+}
+
+
+static SymOpList *sym_ops_for_cell(UnitCell *cell)
+{
+	SymOpList *rawList;
+
+	char *pg = get_chiral_holohedry(cell);
+	rawList = get_pointgroup(pg);
+	free(pg);
+	
+	return rawList;
+}
+
+static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
+  				gsl_matrix *sub, gsl_matrix *mul,
+				double *score, struct TakeTwoCell *cell)
+{
+	double tr;
+
+	gsl_matrix_memcpy(sub, rot1);
+	gsl_matrix_sub(sub, rot2);  /* sub = rot1 - rot2 */
+
+	gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, sub, sub, 0.0, mul);
+
+	tr = matrix_trace(mul);
+	if (score != NULL) *score = tr;
+
+	return (tr < cell->options->trace_tol);
+}
+
+static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
+                                     struct TakeTwoCell *cell)
+{
+	int i;
+	
+	gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+	gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+
+	for (i = 0; i < cell->numOps; i++) {
+		gsl_matrix *testRot = gsl_matrix_alloc(3, 3);
+		gsl_matrix *symOp = cell->rotSymOps[i];
+		
+		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, rot1, symOp,
+		0.0, testRot);
+		
+		if (rot_mats_are_similar(testRot, rot2, sub, mul, NULL, cell)) {
+			gsl_matrix_free(testRot);
+			gsl_matrix_free(sub);
+			gsl_matrix_free(mul);
+			return 1;
+		}
+	
+		gsl_matrix_free(testRot);	
+	}
+	
+	gsl_matrix_free(sub);
+	gsl_matrix_free(mul);
+	
+	return 0;
+}
+
+static void rotation_between_vectors(struct rvec a, struct rvec b,
+				     gsl_matrix *twizzle)
+{
+	double th = rvec_angle(a, b);
+	struct rvec c = rvec_cross(a, b);
+	normalise_rvec(&c);
+	rotation_around_axis(c, th, twizzle);
+}
+
+
+static void rvec_to_gsl(gsl_vector *newVec, struct rvec v)
+{
+	gsl_vector_set(newVec, 0, v.u);
+	gsl_vector_set(newVec, 1, v.v);
+	gsl_vector_set(newVec, 2, v.w);
+}
+
+
+struct rvec gsl_to_rvec(gsl_vector *a)
+{
+	struct rvec v;
+	v.u = gsl_vector_get(a, 0);
+	v.v = gsl_vector_get(a, 1);
+	v.w = gsl_vector_get(a, 2);
+	return v;
+}
+
+
+/* Code me in gsl_matrix language to copy the contents of the function
+ * in cppxfel (IndexingSolution::createSolution). This function is quite
+ * intensive on the number crunching side so simple angle checks are used
+ * to 'pre-scan' vectors beforehand. */
+static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
+				    struct rvec cell1, struct rvec cell2,
+				    struct TakeTwoCell *cell)
+{
+	gsl_matrix *fullMat;
+	rvec_to_gsl(cell->vec1Tmp, cell2);
+
+//	gsl_vector *cell2v = rvec_to_gsl(cell2);
+//	gsl_vector *cell2vr = gsl_vector_calloc(3);
+
+	normalise_rvec(&obs1);
+	normalise_rvec(&obs2);
+	normalise_rvec(&cell1);
+	normalise_rvec(&cell2);
+
+	/* Rotate reciprocal space so that the first simulated vector lines up
+	 * with the observed vector. */
+	rotation_between_vectors(cell1, obs1, cell->twiz1Tmp);
+
+	normalise_rvec(&obs1);
+
+	/* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
+	gsl_blas_dgemv(CblasNoTrans, 1.0, cell->twiz1Tmp, cell->vec1Tmp,
+		       0.0, cell->vec2Tmp);
+
+	/* Now we twirl around the firstAxisUnit until the rotated simulated
+	 * vector matches the second observed vector as closely as possible. */
+
+	closest_rot_mat(gsl_to_rvec(cell->vec2Tmp), obs2, obs1, cell->twiz2Tmp);
+
+	/* We want to apply the first matrix and then the second matrix,
+	 * so we multiply these. */
+	fullMat = gsl_matrix_calloc(3, 3);
+	gsl_blas_dgemm(CblasTrans, CblasTrans, 1.0,
+		       cell->twiz1Tmp, cell->twiz2Tmp, 0.0, fullMat);
+	gsl_matrix_transpose(fullMat);
+
+	return fullMat;
+}
+
+
+static int obs_vecs_share_spot(struct SpotVec *her_obs, struct SpotVec *his_obs)
+{
+	if ( (her_obs->her_rlp == his_obs->her_rlp) ||
+	    (her_obs->her_rlp == his_obs->his_rlp) ||
+	    (her_obs->his_rlp == his_obs->her_rlp) ||
+	    (her_obs->his_rlp == his_obs->his_rlp) ) {
+		return 1;
+	}
+
+	return 0;
+}
+
+
+static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx,
+				   int *members, int num)
+{
+	int i;
+	
+	struct SpotVec *her_obs = &obs_vecs[test_idx];
+
+	for ( i=0; i<num; i++ ) {
+		struct SpotVec *his_obs = &obs_vecs[members[i]];
+
+		int shares = obs_vecs_share_spot(her_obs, his_obs);
+
+		if ( shares ) return 1;
+	}
+
+	return 0;
+}
+
+
+static int obs_vecs_match_angles(struct SpotVec *her_obs,
+                                 struct SpotVec *his_obs,
+                                 int **her_match_idxs, int **his_match_idxs,
+                                 int *match_count, struct TakeTwoCell *cell)
+{
+        int i, j;
+        *match_count = 0;
+
+        double min_angle = deg2rad(2.5);
+        double max_angle = deg2rad(187.5);
+
+        /* calculate angle between observed vectors */
+        double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec);
+
+        /* calculate angle between all potential theoretical vectors */
+
+        for ( i=0; i<her_obs->match_num; i++ ) {
+        for ( j=0; j<his_obs->match_num; j++ ) {
+
+                struct rvec *her_match = &her_obs->matches[i];
+                struct rvec *his_match = &his_obs->matches[j];
+
+                double theory_angle = rvec_angle(*her_match,
+                                                 *his_match);
+
+                /* is this angle a match? */
+
+                double angle_diff = fabs(theory_angle - obs_angle);
+
+                if ( angle_diff < cell->options->angle_tol ) {
+                        // in the case of a brief check only            
+                        if (!her_match_idxs || !his_match_idxs) {
+                                return 1;
+                        }
+                
+                        /* If the angles are too close to 0
+                           or 180, one axis ill-determined */
+                        if (theory_angle < min_angle ||
+                            theory_angle > max_angle) {
+                                continue;
+                        }
+
+                        // check the third vector
+                        
+                        struct rvec theory_diff = diff_vec(*his_match, *her_match);
+                        struct rvec obs_diff = diff_vec(his_obs->obsvec,
+                                                  her_obs->obsvec);
+                        
+                        theory_angle = rvec_angle(*her_match,
+                                                 theory_diff);
+                        obs_angle = rvec_angle(her_obs->obsvec, obs_diff);
+
+                        if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) {
+                                continue;
+                        }
+                        
+                        theory_angle = rvec_angle(*his_match,
+                                                 theory_diff);
+                        obs_angle = rvec_angle(his_obs->obsvec, obs_diff);
+                
+                        if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) {
+                                continue;
+                        }
+                        
+                        size_t new_size = (*match_count + 1) *
+                        sizeof(int);
+                        if (her_match_idxs && his_match_idxs)
+                        {
+                                /* Reallocate the array to fit in another match */
+                                int *temp_hers;
+                                int *temp_his;
+                                temp_hers = realloc(*her_match_idxs, new_size);
+                                temp_his = realloc(*his_match_idxs, new_size);
+
+                                if ( temp_hers == NULL || temp_his == NULL ) {
+                                        apologise();
+                                }
+
+                                (*her_match_idxs) = temp_hers;
+                                (*his_match_idxs) = temp_his;
+
+                                (*her_match_idxs)[*match_count] = i;
+                                (*his_match_idxs)[*match_count] = j;
+                        }
+
+                        (*match_count)++;
+                }
+        }
+        }
+
+        return (*match_count > 0);
+}
+
+/* ------------------------------------------------------------------------
+ * core functions regarding the meat of the TakeTwo algorithm (Level 2)
+ * ------------------------------------------------------------------------*/
+
+static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
+                                    int *obs_members, int *match_members,
+                                    int member_num, struct TakeTwoCell *cell)
+{
+	gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+	gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+	
+	gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2)
+	                                - member_num);
+
+	int i, j, count;
+
+	count = 0;
+	for ( i=0; i<1; i++ ) {
+		for ( j=0; j<member_num; j++ ) {
+			if (i == j) continue;
+			struct SpotVec i_vec = obs_vecs[obs_members[i]];
+			struct SpotVec j_vec = obs_vecs[obs_members[j]];
+
+			struct rvec i_obsvec = i_vec.obsvec;
+			struct rvec j_obsvec = j_vec.obsvec;
+			struct rvec i_cellvec = i_vec.matches[match_members[i]];
+			struct rvec j_cellvec = j_vec.matches[match_members[j]];
+
+			rotations[count] = generate_rot_mat(i_obsvec, j_obsvec,
+							    i_cellvec, j_cellvec,
+							    cell);
+
+			count++;
+		}
+	}
+
+	double min_score = FLT_MAX;
+	int min_rot_index = 0;
+
+	for (i=0; i<count; i++) {
+		double current_score = 0;
+		for (j=0; j<count; j++) {
+			double addition;
+			rot_mats_are_similar(rotations[i], rotations[j],
+					     sub, mul,
+					     &addition, cell);
+
+			current_score += addition;
+		}
+
+		if (current_score < min_score) {
+			min_score = current_score;
+			min_rot_index = i;
+		}
+	}
+
+	gsl_matrix_memcpy(rot, rotations[min_rot_index]);
+
+	for (i=0; i<count; i++) {
+		gsl_matrix_free(rotations[i]);
+	}
+
+	free(rotations);
+	gsl_matrix_free(sub);
+	gsl_matrix_free(mul);
+
+	return 1;
+}
+
+static int weed_duplicate_matches(struct SpotVec *her_obs,
+                                   struct SpotVec *his_obs,
+                                   int **her_match_idxs, int **his_match_idxs,
+                                   int *match_count, struct TakeTwoCell *cell)
+{
+	int num_occupied = 0;
+	gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *));
+	
+	if (old_mats == NULL)
+	{
+		apologise();
+		return 0;
+	}
+
+	signed int i, j;
+	int duplicates = 0;
+	
+	for (i = *match_count - 1; i >= 0; i--) {
+		int her_match = (*her_match_idxs)[i];
+		int his_match = (*his_match_idxs)[i];
+	
+		struct rvec i_obsvec = her_obs->obsvec;
+		struct rvec j_obsvec = his_obs->obsvec;
+		struct rvec i_cellvec = her_obs->matches[her_match];
+		struct rvec j_cellvec = his_obs->matches[his_match];
+
+		gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
+		                                   i_cellvec, j_cellvec, cell);
+
+		int found = 0;
+		
+		for (j = 0; j < num_occupied; j++) {
+			if (old_mats[j] && mat &&
+			    symm_rot_mats_are_similar(old_mats[j], mat, cell))
+			{
+				// we have found a duplicate, so flag as bad.
+				(*her_match_idxs)[i] = -1;
+				(*his_match_idxs)[i] = -1;
+				found = 1;
+	
+				duplicates++;
+
+				gsl_matrix_free(mat);
+				break;
+			}
+		}
+		
+		if (!found) {
+			// we have not found a duplicate, add to list.
+			old_mats[num_occupied] = mat;
+			num_occupied++;
+		}
+	}
+	
+	for (i = 0; i < num_occupied; i++) {
+		if (old_mats[i]) {
+			gsl_matrix_free(old_mats[i]);
+		}
+	}
+
+	free(old_mats);
+	
+	return 1;
+}
+
+static signed int find_next_index(gsl_matrix *rot, int *obs_members,
+				  int *match_members, int start, int member_num,
+				  int *match_found, struct TakeTwoCell *cell)
+{
+	struct SpotVec *obs_vecs = *(cell->obs_vecs);
+	int obs_vec_count = cell->obs_vec_count;
+	gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+	gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+	
+	int i, j, k;
+	
+	for ( i=start; i<obs_vec_count; i++ ) {
+
+		/* first we check for a shared spot - harshest condition */
+		int shared = obs_shares_spot_w_array(obs_vecs, i, obs_members,
+						     member_num);
+
+		if ( !shared ) continue;
+		
+		int skip = 0;
+		for ( j=0; j<member_num && skip == 0; j++ ) {
+			if (i == obs_members[j]) {
+				skip = 1;
+			}
+		}
+		
+		if (skip) {
+			continue;
+		}
+
+		int all_ok = 1;
+		int matched = -1;
+	
+		for ( j=0; j<member_num && all_ok; j++ ) {
+		
+			struct SpotVec *me = &obs_vecs[i];
+			struct SpotVec *you = &obs_vecs[obs_members[j]];
+			struct rvec you_cell = you->matches[match_members[j]];
+
+			struct rvec me_obs = me->obsvec;
+			struct rvec you_obs = you->obsvec;
+
+			int one_is_okay = 0;
+
+			for ( k=0; k<me->match_num; k++ ) {
+				
+				gsl_matrix *test_rot;
+		
+				struct rvec me_cell = me->matches[k];								
+
+				test_rot = generate_rot_mat(me_obs,
+					    you_obs, me_cell, you_cell,
+					    cell);			
+
+				double trace = 0;
+				int ok = rot_mats_are_similar(rot, test_rot,
+						       sub, mul, &trace, cell);
+				
+				gsl_matrix_free(test_rot);
+
+				if (ok) {
+					one_is_okay = 1;
+					
+					if (matched >= 0 && k == matched) {
+						*match_found = k;
+					} else if (matched < 0) {
+						matched = k;
+					} else {
+						one_is_okay = 0;
+						break;
+					}					
+				}
+			}
+			
+			if (!one_is_okay) {
+				all_ok = 0;
+				break;
+			}
+		}
+
+
+		if (all_ok) {
+		
+			for ( j=0; j<member_num; j++ ) {
+			//	STATUS("%i ", obs_members[j]);
+			}
+			//STATUS("%i\n", i);
+			
+			return i;
+		}
+	}
+
+	/* give up. */
+
+	return -1;
+}
+
+
+static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2,
+			int match_idx1, int match_idx2, int *max_members,
+			struct TakeTwoCell *cell)
+{
+	struct SpotVec *obs_vecs = *(cell->obs_vecs);
+	int obs_vec_count = cell->obs_vec_count;
+
+	/* indices of members of the self-consistent network of vectors */
+	int obs_members[MAX_NETWORK_MEMBERS];
+	int match_members[MAX_NETWORK_MEMBERS];
+
+	/* initialise the ones we know already */
+	obs_members[0] = obs_idx1;
+	obs_members[1] = obs_idx2;
+	match_members[0] = match_idx1;
+	match_members[1] = match_idx2;
+	int member_num = 2;
+	*max_members = 2;
+
+	/* counter for dead ends which must not exceed MAX_DEAD_ENDS
+	 * before it is reset in an additional branch */
+	int dead_ends = 0;
+
+	/* we start from 0 */
+	int start = 0;
+
+	while ( 1 ) {
+
+		if (start > obs_vec_count) {
+			return 0;
+		}
+
+		int match_found = -1;
+
+		signed int next_index = find_next_index(rot, obs_members,
+							match_members,
+							start, member_num,
+							&match_found, cell);
+
+		if ( member_num < 2 ) {
+			return 0;
+		}
+
+		if ( next_index < 0 ) {
+			/* If there have been too many dead ends, give up
+			 * on indexing altogether.
+			 **/
+			if ( dead_ends > MAX_DEAD_ENDS ) {
+				break;
+			}
+
+			/* We have not had too many dead ends. Try removing
+			 the last member and continue. */
+			start = obs_members[member_num - 1] + 1;
+			member_num--;
+			dead_ends++;
+
+			continue;
+		}
+
+		/* we have elongated membership - so reset dead_ends counter */
+	//	dead_ends = 0;
+
+		obs_members[member_num] = next_index;
+		match_members[member_num] = match_found;
+
+		member_num++;
+		
+		if (member_num > *max_members) {
+			*max_members = member_num;
+		}
+
+		/* If member_num is high enough, we want to return a yes */
+		if ( member_num > cell->options->member_thresh ) break;
+	}
+
+	finish_solution(rot, obs_vecs, obs_members,
+			  match_members, member_num, cell);
+	
+	return ( member_num > cell->options->member_thresh );
+}
+
+
+static int start_seed(int i, int j, int i_match, int j_match,
+		      gsl_matrix **rotation, int *max_members,
+		      struct TakeTwoCell *cell)
+{
+	struct SpotVec *obs_vecs = *(cell->obs_vecs);
+
+	gsl_matrix *rot_mat;
+	
+	rot_mat = generate_rot_mat(obs_vecs[i].obsvec,
+				   obs_vecs[j].obsvec,
+				   obs_vecs[i].matches[i_match],
+				   obs_vecs[j].matches[j_match],
+				   cell);
+
+	/* Try to expand this rotation matrix to a larger network */
+
+	int success = grow_network(rot_mat, i, j, i_match, j_match, max_members,
+				   cell);
+
+	/* return this matrix and if it was immediately successful */
+	*rotation = rot_mat;
+
+	return success;
+}
+
+static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell)
+{
+	struct SpotVec *obs_vecs = *(cell->obs_vecs);
+	int obs_vec_count = cell->obs_vec_count;
+	
+	/* META: Maximum achieved maximum network membership */
+	int max_max_members = 0;
+	gsl_matrix *best_rotation = NULL;
+
+//	unsigned long start_time = time(NULL);
+
+	/* loop round pairs of vectors to try and find a suitable
+	 * seed to start building a self-consistent network of vectors
+	 */
+	int i, j;
+
+	for ( i=0; i<obs_vec_count; i++ ) {
+		for ( j=0; j<i; j++ ) {
+
+			/** Check to see if there is a shared spot - opportunity
+			 * for optimisation by generating a look-up table
+			 * by spot instead of by vector.
+			 */
+			int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]);
+
+			if ( !shared ) continue;
+
+			/* cell vector index matches stored in i, j and total
+			 * number stored in int matches.
+			 */
+			int *i_idx = 0;
+			int *j_idx = 0;
+			int matches;
+				
+			/* Check to see if any angles match from the cell vectors */
+			obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j],
+					      &i_idx, &j_idx, &matches, cell);
+			if ( matches == 0 ) {
+				free(i_idx);
+				free(j_idx);
+				continue;
+			}
+
+			/* Weed out the duplicate seeds (from symmetric
+			 * reflection pairs)
+			 */
+
+			weed_duplicate_matches(&obs_vecs[i], &obs_vecs[j],
+					       &i_idx, &j_idx, &matches, cell);
+
+			/* We have seeds! Pass each of them through the seed-starter  */
+			/* If a seed has the highest achieved membership, make note...*/
+			int k;
+			for ( k=0; k<matches; k++ ) {
+				if (i_idx[k] < 0 || j_idx[k] < 0) {
+					continue;
+				}
+
+				int max_members = 0;
+
+				int success = start_seed(i, j,
+							 i_idx[k], j_idx[k],
+							 rotation, &max_members,
+							 cell);
+
+				if (success) {
+					free(i_idx); free(j_idx);
+					gsl_matrix_free(best_rotation);
+					return success;
+				} else {
+					if (max_members > max_max_members) {
+						max_max_members = max_members;
+						gsl_matrix_free(best_rotation);
+						best_rotation = *rotation;
+						*rotation = NULL;
+					} else {
+						gsl_matrix_free(*rotation);
+						*rotation = NULL;
+					}
+				}
+			}
+
+			free(i_idx);
+			free(j_idx);
+		}
+	} /* yes this } is meant to be here */
+
+	*rotation = best_rotation;
+	return (best_rotation != NULL);
+}
+
+static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz,
+                           double bsx, double bsy, double bsz,
+                           double csx, double csy, double csz)
+{
+	gsl_matrix_set(mat, 0, 0, asx);
+	gsl_matrix_set(mat, 0, 1, asy);
+	gsl_matrix_set(mat, 0, 2, asz);
+	gsl_matrix_set(mat, 1, 0, bsx);
+	gsl_matrix_set(mat, 1, 1, bsy);
+	gsl_matrix_set(mat, 1, 2, bsz);
+	gsl_matrix_set(mat, 2, 0, csx);
+	gsl_matrix_set(mat, 2, 1, csy);
+	gsl_matrix_set(mat, 2, 2, csz);
+}
+
+static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell)
+{
+	SymOpList *rawList = sym_ops_for_cell(ttCell->cell);
+
+	/* Now we must convert these into rotation matrices rather than hkl
+	 * transformations (affects triclinic, monoclinic, rhombohedral and
+	 * hexagonal space groups only) */
+	 
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+
+	gsl_matrix *recip = gsl_matrix_alloc(3, 3);
+	gsl_matrix *cart = gsl_matrix_alloc(3, 3);
+	cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
+						&bsz, &csx, &csy, &csz);
+
+	set_gsl_matrix(recip, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+	cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
+						&bsz, &csx, &csy, &csz);
+	
+	set_gsl_matrix(cart, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+	int i, j, k;
+	int numOps = num_equivs(rawList, NULL);
+	
+	ttCell->rotSymOps = malloc(numOps * sizeof(gsl_matrix *));
+	ttCell->numOps = numOps;
+	
+	if (ttCell->rotSymOps == NULL) {
+		apologise();
+		return 0;
+	}
+	
+	for (i = 0; i < numOps; i++)
+	{
+		gsl_matrix *symOp = gsl_matrix_alloc(3, 3);
+		IntegerMatrix *op = get_symop(rawList, NULL, i);
+		
+		for (j = 0; j < 3; j++) {
+		for (k = 0; k < 3; k++) {
+			gsl_matrix_set(symOp, j, k, intmat_get(op, j, k));
+		}
+		}
+		
+		gsl_matrix *first = gsl_matrix_calloc(3, 3);
+		gsl_matrix *second = gsl_matrix_calloc(3, 3);
+		
+		/* Each equivalence is of the form:
+	 	   cartesian * symOp * reciprocal.
+	 	   First multiplication: symOp * reciprocal */
+	 	   
+		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
+		               1.0, symOp, recip,
+		               0.0, first);
+		
+		/* Second multiplication: cartesian * first */
+		
+		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
+		               1.0, cart, first,
+		               0.0, second);
+		
+		ttCell->rotSymOps[i] = second;
+		
+		gsl_matrix_free(symOp);
+		gsl_matrix_free(first);
+	}
+
+	gsl_matrix_free(cart);
+	gsl_matrix_free(recip);
+	
+	free_symoplist(rawList);
+	
+	return 1;
+}
+
+
+struct sortme
+{
+	struct rvec v;
+	double dist;
+};
+
+static int sort_func(const void *av, const void *bv)
+{
+	struct sortme *a = (struct sortme *)av;
+	struct sortme *b = (struct sortme *)bv;
+	return a->dist > b->dist;
+}
+
+
+static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
+				  struct SpotVec *obs_vecs, int obs_vec_count,
+				  int is_asymmetric, struct TakeTwoCell *cell)
+{
+	int i, j;
+
+	for ( i=0; i<obs_vec_count; i++ ) {
+
+		int count = 0;
+		struct sortme *for_sort = NULL;
+		
+		for ( j=0; j<cell_vec_count; j++ ) {
+			/* get distance for unit cell vector */
+			double cell_length = rvec_length(cell_vecs[j]);
+			double obs_length = obs_vecs[i].distance;
+
+			/* check if this matches the observed length */
+			double dist_diff = fabs(cell_length - obs_length);
+
+			if ( dist_diff > cell->options->len_tol ) continue;
+
+			/* we have a match, add to array! */
+
+			size_t new_size = (count+1)*sizeof(struct sortme);
+			for_sort = realloc(for_sort, new_size);
+
+			if ( for_sort == NULL ) return 0;
+
+			for_sort[count].v = cell_vecs[j];
+			for_sort[count].dist = dist_diff;
+			count++;
+
+		}
+
+		/* Pointers to relevant things */
+		
+		struct rvec **match_array;
+		int *match_count;
+		
+		if (!is_asymmetric) {
+			match_array = &(obs_vecs[i].matches);
+			match_count = &(obs_vecs[i].match_num);
+		} else {
+			match_array = &(obs_vecs[i].asym_matches);
+			match_count = &(obs_vecs[i].asym_match_num);		
+		}
+
+		/* Sort in order to get most agreeable matches first */
+		qsort(for_sort, count, sizeof(struct sortme), sort_func);
+		*match_array = malloc(count*sizeof(struct rvec));
+		*match_count = count;
+		for ( j=0; j<count; j++ ) {
+			(*match_array)[j] = for_sort[j].v;
+			
+		}
+		free(for_sort);
+	}
+
+	return 1;
+}
+
+static int compare_spot_vecs(const void *av, const void *bv)
+{
+	struct SpotVec *a = (struct SpotVec *)av;
+	struct SpotVec *b = (struct SpotVec *)bv;
+	return a->distance > b->distance;
+}
+
+static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
+			     struct SpotVec **obs_vecs, int *obs_vec_count)
+{
+	int i, j;
+	int count = 0;
+
+	/* maximum distance squared for comparisons */
+	double max_sq_length = pow(MAX_RECIP_DISTANCE, 2);
+
+	for ( i=0; i<rlp_count-1 && count < MAX_OBS_VECTORS; i++ ) {
+		for ( j=i+1; j<rlp_count; j++ ) {
+
+
+			/* calculate difference vector between rlps */
+			struct rvec diff = diff_vec(rlps[i], rlps[j]);
+
+			/* are these two far from each other? */
+			double sqlength = sq_length(diff);
+
+			if ( sqlength > max_sq_length ) continue;
+
+			count++;
+
+			struct SpotVec *temp_obs_vecs;
+			temp_obs_vecs = realloc(*obs_vecs,
+						count*sizeof(struct SpotVec));
+
+			if ( temp_obs_vecs == NULL ) {
+				return 0;
+			} else {
+				*obs_vecs = temp_obs_vecs;
+
+				/* initialise all SpotVec struct members */
+
+				struct SpotVec spot_vec;
+				spot_vec.obsvec = diff;
+				spot_vec.distance = sqrt(sqlength);
+				spot_vec.matches = NULL;
+				spot_vec.match_num = 0;
+				spot_vec.her_rlp = &rlps[i];
+				spot_vec.his_rlp = &rlps[j];
+
+				(*obs_vecs)[count - 1] = spot_vec;
+			}
+		}
+	}
+
+	/* Sort such that the shortest and least error-prone distances
+	   are searched first.
+	 */
+	qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs);
+
+	*obs_vec_count = count;
+
+	return 1;
+}
+
+
+static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
+				struct rvec **asym_vecs, int *vec_count,
+				int *asym_vec_count)
+{
+	double a, b, c, alpha, beta, gamma;
+	int h_max, k_max, l_max;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+			    &bsx, &bsy, &bsz,
+			    &csx, &csy, &csz);
+
+	SymOpList *rawList = sym_ops_for_cell(cell);
+
+	cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+	/* find maximum Miller (h, k, l) indices for a given resolution */
+	h_max = MAX_RECIP_DISTANCE * a;
+	k_max = MAX_RECIP_DISTANCE * b;
+	l_max = MAX_RECIP_DISTANCE * c;
+
+	int h, k, l;
+	int _h, _k, _l;
+	int count = 0;
+	int asym_count = 0;
+
+	for ( h=-h_max; h<=+h_max; h++ ) {
+	for ( k=-k_max; k<=+k_max; k++ ) {
+	for ( l=-l_max; l<=+l_max; l++ ) {
+
+		struct rvec cell_vec;
+
+		/* Exclude systematic absences from centering concerns */
+		if ( forbidden_reflection(cell, h, k, l) ) continue;
+
+		int asymmetric = 0;
+		get_asymm(rawList, h, k, l, &_h, &_k, &_l);
+		
+		if (h == _h && k == _k && l == _l) {
+			asymmetric = 1;
+			asym_count++;
+		}
+		
+		cell_vec.u = h*asx + k*bsx + l*csx;
+		cell_vec.v = h*asy + k*bsy + l*csy;
+		cell_vec.w = h*asz + k*bsz + l*csz;
+
+		/* add this to our array - which may require expanding */
+		count++;
+
+		struct rvec *temp_cell_vecs;
+		temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec));
+		struct rvec *temp_asym_vecs = NULL;
+
+		if (asymmetric) {
+			temp_asym_vecs = realloc(*asym_vecs,
+			                         count*sizeof(struct rvec));
+		}
+
+		if ( temp_cell_vecs == NULL ) {
+			return 0;
+		} else if (asymmetric && temp_asym_vecs == NULL) {
+			return 0;
+		} else {
+			*cell_vecs = temp_cell_vecs;
+			(*cell_vecs)[count - 1] = cell_vec;
+			
+			if (!asymmetric) {
+				continue;
+			}
+			
+			*asym_vecs = temp_asym_vecs;
+			(*asym_vecs)[asym_count - 1] = cell_vec;
+		}
+	}
+	}
+	}
+
+	*vec_count = count;
+	*asym_vec_count = asym_count;
+	
+	free_symoplist(rawList);
+
+
+	return 1;
+}
+
+
+/* ------------------------------------------------------------------------
+ * cleanup functions - called from run_taketwo().
+ * ------------------------------------------------------------------------*/
+ 
+static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs,
+				     int obs_vec_count)
+{
+	int i;
+	for ( i=0; i<obs_vec_count; i++ ) {
+		free(obs_vecs[i].matches);
+		free(obs_vecs[i].asym_matches);
+	}
+
+	free(obs_vecs);
+}
+
+static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell)
+{
+	int i;
+	for ( i=0; i<ttCell->numOps; i++ ) {
+		gsl_matrix_free(ttCell->rotSymOps[i]);
+	}
+
+	free(ttCell->vec1Tmp);
+	free(ttCell->vec2Tmp);
+	free(ttCell->twiz1Tmp);
+	free(ttCell->twiz2Tmp);
+	free(ttCell->rotSymOps);
+}
+
+
+/* ------------------------------------------------------------------------
+ * external functions - top level functions (Level 1)
+ * ------------------------------------------------------------------------*/
+
+/**
+ * @cell: target unit cell
+ * @rlps: spot positions on detector back-projected into recripocal
+ * space depending on detector geometry etc.
+ * @rlp_count: number of rlps in rlps array.
+ * @rot: pointer to be given an assignment (hopefully) if indexing is
+ * successful.
+ **/
+static UnitCell *run_taketwo(UnitCell *cell, struct taketwo_options *opts,
+			     struct rvec *rlps, int rlp_count)
+{
+	int cell_vec_count = 0;
+	int asym_vec_count = 0;
+	struct rvec *cell_vecs = NULL;
+	struct rvec *asym_vecs = NULL;
+	UnitCell *result;
+	int obs_vec_count = 0;
+	struct SpotVec *obs_vecs = NULL;
+	int success = 0;
+	gsl_matrix *solution = NULL;
+	
+	struct TakeTwoCell ttCell;
+	ttCell.cell = cell;
+	ttCell.options = opts;
+	ttCell.rotSymOps = NULL;
+	ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3);
+	ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3);
+	ttCell.vec1Tmp = gsl_vector_calloc(3);
+	ttCell.vec2Tmp = gsl_vector_calloc(3);
+	ttCell.numOps = 0;
+
+	success = generate_rotation_sym_ops(&ttCell);
+
+	success = gen_theoretical_vecs(cell, &cell_vecs, &asym_vecs,
+	                               &cell_vec_count, &asym_vec_count);
+	if ( !success ) return NULL;
+
+	success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count);
+	if ( !success ) return NULL;
+
+	ttCell.obs_vecs = &obs_vecs;
+	ttCell.obs_vec_count = obs_vec_count;
+
+	success = match_obs_to_cell_vecs(asym_vecs, asym_vec_count,
+					 obs_vecs, obs_vec_count, 1, &ttCell);
+
+	success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count,
+					 obs_vecs, obs_vec_count, 0, &ttCell);
+
+	free(cell_vecs);
+	free(asym_vecs);
+
+	if ( !success ) return NULL;
+
+	find_seed(&solution, &ttCell);
+
+	if ( solution == NULL ) {
+		cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
+		return NULL;
+	}
+
+	result = transform_cell_gsl(cell, solution);
+	gsl_matrix_free(solution);
+	cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
+	cleanup_taketwo_cell(&ttCell);
+	
+	return result;
+}
+
+
+/* CrystFEL interface hooks */
+
+int taketwo_index(struct image *image, struct taketwo_options *opts, void *priv)
+{
+	Crystal *cr;
+	UnitCell *cell;
+	struct rvec *rlps;
+	int n_rlps = 0;
+	int i;
+	struct taketwo_private *tp = (struct taketwo_private *)priv;
+
+	/* Replace negative values in options with defaults */
+
+	if (opts->member_thresh < 0) {
+		opts->member_thresh = NETWORK_MEMBER_THRESHOLD;
+	}
+
+	if (opts->len_tol < 0) {
+		opts->len_tol = RECIP_TOLERANCE;
+	}
+	
+	if (opts->angle_tol < 0) {
+		opts->angle_tol = ANGLE_TOLERANCE;
+	}
+
+	if (opts->trace_tol < 0) {
+		opts->trace_tol = TRACE_TOLERANCE;
+	}
+	
+	opts->trace_tol = sqrt(4.0*(1.0-cos(opts->trace_tol)));
+
+	rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec));
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+		struct imagefeature *pk = image_get_feature(image->features, i);
+		if ( pk == NULL ) continue;
+		rlps[n_rlps].u = pk->rx;
+		rlps[n_rlps].v = pk->ry;
+		rlps[n_rlps].w = pk->rz;
+		n_rlps++;
+	}
+	rlps[n_rlps].u = 0.0;
+	rlps[n_rlps].v = 0.0;
+	rlps[n_rlps++].w = 0.0;
+
+	cell = run_taketwo(tp->cell, opts, rlps, n_rlps);
+	free(rlps);
+	if ( cell == NULL ) return 0;
+
+	cr = crystal_new();
+	if ( cr == NULL ) {
+		ERROR("Failed to allocate crystal.\n");
+		return 0;
+	}
+
+	crystal_set_cell(cr, cell);
+
+	if ( tp->indm & INDEXING_CHECK_PEAKS ) {
+		if ( !peak_sanity_check(image, &cr, 1) ) {
+			cell_free(cell);
+			crystal_free(cr);
+		//	STATUS("Rubbish!!\n");
+
+			return 0;
+		} else {
+		//	STATUS("That's good!\n");
+		}
+	}
+
+	image_add_crystal(image, cr);
+
+	return 1;
+}
+
+
+void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell,
+                      struct detector *det, float *ltl)
+{
+	struct taketwo_private *tp;
+
+	/* Flags that TakeTwo knows about */
+	*indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS
+	| INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS
+	| INDEXING_CONTROL_FLAGS;
+
+	if ( !( (*indm & INDEXING_USE_LATTICE_TYPE)
+	       && (*indm & INDEXING_USE_CELL_PARAMETERS)) )
+	{
+		ERROR("TakeTwo indexing requires cell and lattice type "
+		      "information.\n");
+		return NULL;
+	}
+
+	if ( cell == NULL ) {
+		ERROR("TakeTwo indexing requires a unit cell.\n");
+		return NULL;
+	}
+
+	STATUS("*******************************************************************\n");
+	STATUS("*****                    Welcome to TakeTwo                   *****\n");
+	STATUS("*******************************************************************\n");
+	STATUS("      If you use these indexing results, please keep a roof\n");
+	STATUS("           over the author's head by citing this paper.\n\n");
+
+	STATUS("o     o     o     o     o     o     o     o     o     o     o     o\n");
+	STATUS("   o     o     o     o     o     o     o     o     o     o     o   \n");
+	STATUS("o                                                                 o\n");
+	STATUS("   o                      The citation is:                     o   \n");
+	STATUS("o           Ginn et al., Acta Cryst. (2016). D72, 956-965         o\n");
+	STATUS("   o                         Thank you!                        o   \n");
+	STATUS("o                                                                 o\n");
+	STATUS("   o     o     o     o     o     o     o     o     o     o     o   \n");
+	STATUS("o     o     o     o     o     o     o     o     o     o     o     o\n");
+
+
+	STATUS("\n");
+
+	tp = malloc(sizeof(struct taketwo_private));
+	if ( tp == NULL ) return NULL;
+
+	tp->ltl = ltl;
+	tp->cell = cell;
+	tp->indm = *indm;
+
+	return tp;
+}
+
+
+void taketwo_cleanup(IndexingPrivate *pp)
+{
+	struct taketwo_private *tp = (struct taketwo_private *)pp;
+	free(tp);
+}
+
diff --git a/libcrystfel/src/taketwo.h b/libcrystfel/src/taketwo.h
new file mode 100644
index 00000000..a03292a8
--- /dev/null
+++ b/libcrystfel/src/taketwo.h
@@ -0,0 +1,51 @@
+/*
+ * taketwo.h
+ *
+ * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL
+ *
+ * Copyright © 2016-2017 Helen Ginn
+ * Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2016      Helen Ginn <helen@strubi.ox.ac.uk>
+ *   2016-2017 Thomas White <taw@physics.org>
+ *
+ * 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 TAKETWO_H
+#define TAKETWO_H
+
+#include "cell.h"
+#include "index.h"
+
+struct taketwo_options
+{
+	int member_thresh;
+	double len_tol;
+	double angle_tol;
+	double trace_tol;
+};
+
+
+extern void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell,
+                             struct detector *det, float *ltl);
+extern int taketwo_index(struct image *image, struct taketwo_options *opts, void *priv);
+extern void taketwo_cleanup(IndexingPrivate *pp);
+
+#endif /* TAKETWO_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index ce87fde6..d232b66f 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -62,7 +62,7 @@
 #include "reflist-utils.h"
 #include "cell-utils.h"
 #include "integration.h"
-
+#include "taketwo.h"
 #include "im-sandbox.h"
 
 
@@ -171,6 +171,11 @@ static void show_help(const char *s)
 "                          Given as a list of comma separated list of \n"
 "                            indexer specific, key word arguments.\n"
 "                          Example: \"arg1=10,arg2=500\" \n"
+"\nLow-level options for the TakeTwo indexer:\n\n"
+"     --taketwo-member-threshold Minimum number of members in network\n"
+"     --taketwo-len-tolerance    Reciprocal space length tolerance (1/A)\n"
+"     --taketwo-angle-tolerance  Reciprocal space angle tolerance (in degrees)\n"
+"     --taketwo-trace-tolerance  Rotation matrix equivalence tolerance (in degrees)\n"
 );
 }
 
@@ -276,6 +281,10 @@ int main(int argc, char *argv[])
 	iargs.fix_divergence = -1.0;
 	iargs.felix_options = NULL;
 	iargs.profile = 0;
+	iargs.taketwo_opts.member_thresh = -1;
+	iargs.taketwo_opts.len_tol = -1.0;
+	iargs.taketwo_opts.angle_tol = -1.0;
+	iargs.taketwo_opts.trace_tol = -1.0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -344,6 +353,10 @@ int main(int argc, char *argv[])
 	        {"min-res",            1, NULL,               29},
 	        {"max-res",            1, NULL,               30},
 	        {"min-peaks",          1, NULL,               31},
+	        {"taketwo-member-thresh", 1, NULL,            32},
+	        {"taketwo-len-tol",       1, NULL,            33},
+	        {"taketwo-angle-tol",     1, NULL,            34},
+	        {"taketwo-trace-tol",     1, NULL,            35},
 
 		{0, 0, NULL, 0}
 	};
@@ -551,6 +564,42 @@ int main(int argc, char *argv[])
 			iargs.min_peaks = atoi(optarg);
 			break;
 
+			case 32:
+			if ( sscanf(optarg, "%i", &iargs.taketwo_opts.member_thresh) != 1 )
+			{
+				ERROR("Invalid value for --taketwo-member-threshold\n");
+				return 1;
+			}
+			break;
+
+			case 33:
+			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.len_tol) != 1 )
+			{
+				ERROR("Invalid value for --taketwo-len-tolerance\n");
+				return 1;
+			}
+			iargs.taketwo_opts.len_tol *= 1e10;
+			break;
+
+			case 34:
+			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.angle_tol) != 1 )
+			{
+				ERROR("Invalid value for --taketwo-angle-tolerance\n");
+				return 1;
+			}
+			iargs.taketwo_opts.angle_tol = deg2rad(iargs.taketwo_opts.angle_tol);
+			break;
+
+
+			case 35:
+			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.trace_tol) != 1 )
+			{
+				ERROR("Invalid value for --taketwo-trace-tolerance\n");
+				return 1;
+			}
+			iargs.taketwo_opts.trace_tol = deg2rad(iargs.taketwo_opts.trace_tol);
+			break;
+
 			case 0 :
 			break;
 
@@ -803,7 +852,8 @@ int main(int argc, char *argv[])
 
 		iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det,
 		                             iargs.tols, no_refine,
-		                             iargs.felix_options);
+		                             iargs.felix_options,
+		                             &iargs.taketwo_opts);
 		if ( iargs.ipriv == NULL ) {
 			ERROR("Failed to set up indexing system\n");
 			return 1;
diff --git a/src/process_image.h b/src/process_image.h
index ec51c188..3e87b922 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -39,6 +39,7 @@ struct index_args;
 #include "integration.h"
 #include "im-sandbox.h"
 #include "time-accounts.h"
+#include "taketwo.h"
 
 
 enum {
@@ -100,6 +101,7 @@ struct index_args
 	float fix_divergence;
 	char *felix_options;
 	int profile;  /* Whether or not to do wall clock profiling */
+	struct taketwo_options taketwo_opts;
 };