Merge branch 'tom/newcellcheck'

19 files changed, 1354 insertions, 590 deletions

diff --git a/doc/man/cell_tool.1 b/doc/man/cell_tool.1
index 576943d0..86dd9ae0 100644
--- a/doc/man/cell_tool.1
+++ b/doc/man/cell_tool.1
@@ -51,7 +51,7 @@ There are an infinite number of primitive unit cell for any lattice.  This progr
 .PP
 The program will compare the two cells, and report if \fImy_structure.cell\fR can be made to look similar to \fIreference.cell\fR applying any transformation matrix.
 .PP
-The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles.
+The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles.  If either of the unit cells are centered, a primitive version of will be used \fIfor some of the comparison results\fR.  Read the output carefully.
 
 .SH TRANSFORMING A UNIT CELL
 .PP
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index f6e189ad..d92bd6a9 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -574,404 +574,6 @@ UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi)
 }
 
 
-#define MAX_CAND (1024)
-
-static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c)
-{
-	struct rvec aCb;
-	double aCb_dot_c;
-
-	/* "a" cross "b" */
-	aCb.u = a.v*b.w - a.w*b.v;
-	aCb.v = - (a.u*b.w - a.w*b.u);
-	aCb.w = a.u*b.v - a.v*b.u;
-
-	/* "a cross b" dot "c" */
-	aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
-
-	if ( aCb_dot_c > 0.0 ) return 1;
-	return 0;
-}
-
-
-struct cvec {
-	struct rvec vec;
-	float na;
-	float nb;
-	float nc;
-	float fom;
-};
-
-
-static int same_vector(struct cvec a, struct cvec b)
-{
-	if ( a.na != b.na ) return 0;
-	if ( a.nb != b.nb ) return 0;
-	if ( a.nc != b.nc ) return 0;
-	return 1;
-}
-
-
-/* Attempt to make 'cell' fit into 'template' somehow */
-UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
-                     const float *tols, int reduce)
-{
-	signed int n1l, n2l, n3l;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	int i, j;
-	double lengths[3];
-	double angles[3];
-	struct cvec *cand[3];
-	UnitCell *new_cell = NULL;
-	float best_fom = +999999999.9; /* Large number.. */
-	int ncand[3] = {0,0,0};
-	signed int ilow, ihigh;
-	float angtol = deg2rad(tols[3]);
-	UnitCell *cell;
-	UnitCell *template;
-	IntegerMatrix *centering;
-	UnitCell *new_cell_trans;
-
-	/* "Un-center" the template unit cell to make the comparison easier */
-	template = uncenter_cell(template_in, &centering, NULL);
-	if ( template == NULL ) return NULL;
-
-	/* The candidate cell is also uncentered, because it might be centered
-	 * if it came from (e.g.) MOSFLM */
-	cell = uncenter_cell(cell_in, NULL, NULL);
-	if ( cell == NULL ) return NULL;
-
-	if ( cell_get_reciprocal(template, &asx, &asy, &asz,
-	                         &bsx, &bsy, &bsz,
-	                         &csx, &csy, &csz) ) {
-		ERROR("Couldn't get reciprocal cell for template.\n");
-		cell_free(template);
-		cell_free(cell);
-		intmat_free(centering);
-		return NULL;
-	}
-
-	lengths[0] = modulus(asx, asy, asz);
-	lengths[1] = modulus(bsx, bsy, bsz);
-	lengths[2] = modulus(csx, csy, csz);
-
-	angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
-	angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
-	angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
-
-	cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
-	cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
-	cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
-
-	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                         &bsx, &bsy, &bsz,
-	                         &csx, &csy, &csz) ) {
-		ERROR("Couldn't get reciprocal cell.\n");
-		cell_free(template);
-		cell_free(cell);
-		intmat_free(centering);
-		return NULL;
-	}
-
-	if ( reduce ) {
-		ilow = -2;  ihigh = 4;
-	} else {
-		ilow = 0;  ihigh = 1;
-	}
-
-	/* Negative values mean 1/n, positive means n, zero means zero */
-	for ( n1l=ilow; n1l<=ihigh; n1l++ ) {
-	for ( n2l=ilow; n2l<=ihigh; n2l++ ) {
-	for ( n3l=ilow; n3l<=ihigh; n3l++ ) {
-
-		float n1, n2, n3;
-		signed int b1, b2, b3;
-
-		n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
-		n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
-		n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
-
-		if ( !reduce ) {
-			if ( n1l + n2l + n3l > 1 ) continue;
-		}
-
-		/* 'bit' values can be +1 or -1 */
-		for ( b1=-1; b1<=1; b1+=2 ) {
-		for ( b2=-1; b2<=1; b2+=2 ) {
-		for ( b3=-1; b3<=1; b3+=2 ) {
-
-			double tx, ty, tz;
-			double tlen;
-			int i;
-
-			n1 *= b1;  n2 *= b2;  n3 *= b3;
-
-			tx = n1*asx + n2*bsx + n3*csx;
-			ty = n1*asy + n2*bsy + n3*csy;
-			tz = n1*asz + n2*bsz + n3*csz;
-			tlen = modulus(tx, ty, tz);
-
-			/* Test modulus for agreement with moduli of template */
-			for ( i=0; i<3; i++ ) {
-
-				if ( !within_tolerance(lengths[i], tlen,
-				                       tols[i]) )
-				{
-					continue;
-				}
-
-				if ( ncand[i] == MAX_CAND ) {
-					ERROR("Too many cell candidates - ");
-					ERROR("consider tightening the unit ");
-					ERROR("cell tolerances.\n");
-				} else {
-
-					double fom;
-
-					fom = fabs(lengths[i] - tlen);
-
-					cand[i][ncand[i]].vec.u = tx;
-					cand[i][ncand[i]].vec.v = ty;
-					cand[i][ncand[i]].vec.w = tz;
-					cand[i][ncand[i]].na = n1;
-					cand[i][ncand[i]].nb = n2;
-					cand[i][ncand[i]].nc = n3;
-					cand[i][ncand[i]].fom = fom;
-
-					ncand[i]++;
-
-				}
-			}
-
-		}
-		}
-		}
-
-	}
-	}
-	}
-
-	if ( verbose ) {
-		STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
-	}
-
-	for ( i=0; i<ncand[0]; i++ ) {
-	for ( j=0; j<ncand[1]; j++ ) {
-
-		double ang;
-		int k;
-		float fom1;
-
-		if ( same_vector(cand[0][i], cand[1][j]) ) continue;
-
-		/* Measure the angle between the ith candidate for axis 0
-		 * and the jth candidate for axis 1 */
-		ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
-		                    cand[0][i].vec.w, cand[1][j].vec.u,
-		                    cand[1][j].vec.v, cand[1][j].vec.w);
-
-		/* Angle between axes 0 and 1 should be angle 2 */
-		if ( fabs(ang - angles[2]) > angtol ) continue;
-
-		fom1 = fabs(ang - angles[2]);
-
-		for ( k=0; k<ncand[2]; k++ ) {
-
-			float fom2, fom3;
-
-			if ( same_vector(cand[1][j], cand[2][k]) ) continue;
-
-			/* Measure the angle between the current candidate for
-			 * axis 0 and the kth candidate for axis 2 */
-			ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
-			                    cand[0][i].vec.w, cand[2][k].vec.u,
-			                    cand[2][k].vec.v, cand[2][k].vec.w);
-
-			/* ... it should be angle 1 ... */
-			if ( fabs(ang - angles[1]) > angtol ) continue;
-
-			fom2 = fom1 + fabs(ang - angles[1]);
-
-			/* Finally, the angle between the current candidate for
-			 * axis 1 and the kth candidate for axis 2 */
-			ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
-			                    cand[1][j].vec.w, cand[2][k].vec.u,
-			                    cand[2][k].vec.v, cand[2][k].vec.w);
-
-			/* ... it should be angle 0 ... */
-			if ( fabs(ang - angles[0]) > angtol ) continue;
-
-			/* Unit cell must be right-handed */
-			if ( !right_handed_vec(cand[0][i].vec, cand[1][j].vec,
-			                       cand[2][k].vec) ) continue;
-
-			fom3 = fom2 + fabs(ang - angles[0]);
-			fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
-			                   + cand[2][k].fom);
-
-			if ( fom3 < best_fom ) {
-				if ( new_cell != NULL ) free(new_cell);
-				new_cell = cell_new_from_reciprocal_axes(
-				                 cand[0][i].vec, cand[1][j].vec,
-			                         cand[2][k].vec);
-				best_fom = fom3;
-			}
-
-		}
-
-	}
-	}
-
-	free(cand[0]);
-	free(cand[1]);
-	free(cand[2]);
-
-	cell_free(cell);
-
-	/* Reverse the de-centering transformation */
-	if ( new_cell != NULL ) {
-
-		new_cell_trans = cell_transform_intmat(new_cell, centering);
-		cell_free(new_cell);
-		cell_set_lattice_type(new_cell_trans,
-		                      cell_get_lattice_type(template_in));
-		cell_set_centering(new_cell_trans,
-		                   cell_get_centering(template_in));
-		cell_set_unique_axis(new_cell_trans,
-		                     cell_get_unique_axis(template_in));
-
-		cell_free(template);
-		intmat_free(centering);
-
-		return new_cell_trans;
-
-	} else {
-		cell_free(template);
-		intmat_free(centering);
-		return NULL;
-	}
-}
-
-
-UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	int i;
-	double lengths[3];
-	int used[3];
-	struct rvec real_a, real_b, real_c;
-	struct rvec params[3];
-	double alen, blen;
-	float ltl = 5.0;     /* percent */
-	int have_real_a;
-	int have_real_b;
-	int have_real_c;
-	UnitCell *cell;
-	UnitCell *template;
-	IntegerMatrix *to_given_cell;
-	UnitCell *new_cell;
-	UnitCell *new_cell_trans;
-
-	/* "Un-center" the template unit cell to make the comparison easier */
-	template = uncenter_cell(template_in, &to_given_cell, NULL);
-
-	/* The candidate cell is also uncentered, because it might be centered
-	 * if it came from (e.g.) MOSFLM */
-	cell = uncenter_cell(cell_in, NULL, NULL);
-
-	/* Get the lengths to match */
-	if ( cell_get_cartesian(template, &ax, &ay, &az,
-	                                  &bx, &by, &bz,
-	                                  &cx, &cy, &cz) )
-	{
-		ERROR("Couldn't get cell for template.\n");
-		return NULL;
-	}
-	alen = modulus(ax, ay, az);
-	blen = modulus(bx, by, bz);
-
-	/* Get the lengths from the cell and turn them into anonymous vectors */
-	if ( cell_get_cartesian(cell, &ax, &ay, &az,
-	                              &bx, &by, &bz,
-	                              &cx, &cy, &cz) )
-	{
-		ERROR("Couldn't get cell.\n");
-		return NULL;
-	}
-	lengths[0] = modulus(ax, ay, az);
-	lengths[1] = modulus(bx, by, bz);
-	lengths[2] = modulus(cx, cy, cz);
-	used[0] = 0;  used[1] = 0;  used[2] = 0;
-	params[0].u = ax;  params[0].v = ay;  params[0].w = az;
-	params[1].u = bx;  params[1].v = by;  params[1].w = bz;
-	params[2].u = cx;  params[2].v = cy;  params[2].w = cz;
-
-	real_a.u = 0.0;  real_a.v = 0.0;  real_a.w = 0.0;
-	real_b.u = 0.0;  real_b.v = 0.0;  real_b.w = 0.0;
-	real_c.u = 0.0;  real_c.v = 0.0;  real_c.w = 0.0;
-
-	/* Check each vector against a and b */
-	have_real_a = 0;
-	have_real_b = 0;
-	for ( i=0; i<3; i++ ) {
-		if ( within_tolerance(lengths[i], alen, ltl)
-		     && !used[i] && !have_real_a )
-		{
-			used[i] = 1;
-			memcpy(&real_a, &params[i], sizeof(struct rvec));
-			have_real_a = 1;
-		}
-		if ( within_tolerance(lengths[i], blen, ltl)
-		     && !used[i] && !have_real_b )
-		{
-			used[i] = 1;
-			memcpy(&real_b, &params[i], sizeof(struct rvec));
-			have_real_b = 1;
-		}
-	}
-
-	/* Have we matched both a and b? */
-	if ( !(have_real_a && have_real_b) ) return NULL;
-
-	/* "c" is "the other one" */
-	have_real_c = 0;
-	for ( i=0; i<3; i++ ) {
-		if ( !used[i] ) {
-			memcpy(&real_c, &params[i], sizeof(struct rvec));
-			have_real_c = 1;
-		}
-	}
-
-	if ( !have_real_c ) {
-		ERROR("Huh?  Couldn't find the third vector.\n");
-		ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]);
-		return NULL;
-	}
-
-	/* Flip c if not right-handed */
-	if ( !right_handed_vec(real_a, real_b, real_c) ) {
-		real_c.u = -real_c.u;
-		real_c.v = -real_c.v;
-		real_c.w = -real_c.w;
-	}
-
-	new_cell = cell_new_from_direct_axes(real_a, real_b, real_c);
-
-	 /* Reverse the de-centering transformation */
-	new_cell_trans = cell_transform_intmat_inverse(new_cell, to_given_cell);
-	cell_free(new_cell);
-	cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in));
-	cell_set_centering(new_cell, cell_get_centering(template_in));
-	cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in));
-
-	return new_cell_trans;
-}
-
-
 /* Return sin(theta)/lambda = 1/2d.  Multiply by two if you want 1/d */
 double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
 {
@@ -1669,14 +1271,13 @@ double cell_get_volume(UnitCell *cell)
 
 
 /**
- * \param cell1: A UnitCell
- * \param cell2: Another UnitCell
- * \param ltl: Maximum allowable fractional difference in axis lengths
- * \param atl: Maximum allowable difference in reciprocal angles (in radians)
+ * \param cell: A UnitCell
+ * \param reference: Another UnitCell
+ * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians)
  *
  * Compare the two unit cells.  If the real space parameters match to within
- * fractional difference \p ltl, and the inter-axial angles match within \p atl,
- * and the centering matches, this function returns 1.  Otherwise 0.
+ * the specified tolerances, and the centering matches, this function returns 1.
+ * Otherwise 0.
  *
  * This function considers the cell parameters and centering, but ignores the
  * orientation of the cell.  If you want to compare the orientation as well,
@@ -1685,25 +1286,25 @@ double cell_get_volume(UnitCell *cell)
  * \returns non-zero if the cells match.
  *
  */
-int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
-                            float ltl, float atl)
+int compare_cell_parameters(UnitCell *cell, UnitCell *reference,
+                            const double *tols)
 {
 	double a1, b1, c1, al1, be1, ga1;
 	double a2, b2, c2, al2, be2, ga2;
 
 	/* Centering must match: we don't arbitrarte primitive vs centered,
 	 * different cell choices etc */
-	if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
+	if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0;
 
-	cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1);
-	cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2);
+	cell_get_parameters(cell, &a1, &b1, &c1, &al1, &be1, &ga1);
+	cell_get_parameters(reference, &a2, &b2, &c2, &al2, &be2, &ga2);
 
-	if ( !within_tolerance(a1, a2, ltl*100.0) ) return 0;
-	if ( !within_tolerance(b1, b2, ltl*100.0) ) return 0;
-	if ( !within_tolerance(c1, c2, ltl*100.0) ) return 0;
-	if ( fabs(al1-al2) > atl ) return 0;
-	if ( fabs(be1-be2) > atl ) return 0;
-	if ( fabs(ga1-ga2) > atl ) return 0;
+	if ( !within_tolerance(a1, a2, tols[0]*100.0) ) return 0;
+	if ( !within_tolerance(b1, b2, tols[1]*100.0) ) return 0;
+	if ( !within_tolerance(c1, c2, tols[2]*100.0) ) return 0;
+	if ( fabs(al1-al2) > tols[3] ) return 0;
+	if ( fabs(be1-be2) > tols[4] ) return 0;
+	if ( fabs(ga1-ga2) > tols[5] ) return 0;
 
 	return 1;
 }
@@ -1719,83 +1320,94 @@ static double moduli_check(double ax, double ay, double az,
 
 
 /**
- * \param cell1: A UnitCell
- * \param cell2: Another UnitCell
- * \param ltl: Maximum allowable fractional difference in axis lengths
- * \param atl: Maximum allowable difference in angles (in radians)
+ * \param cell: A UnitCell
+ * \param reference: Another UnitCell
+ * \param tols: Pointer to six tolerance values (see below)
  *
  * Compare the two unit cells.  If the axes match in length (to within
- * fractional difference \p ltl) and the axes are aligned to within \p atl radians,
- * this function returns non-zero.
+ * the specified tolerances), this function returns non-zero.
  *
  * This function compares the orientation of the cell as well as the parameters.
  * If you just want to see if the parameters are the same, use
  * compare_cell_parameters() instead.
  *
- * The cells \p a and \p b must have the same centering.  Otherwise, this function
- * always returns zero.
+ * The comparison is done by checking that the lengths of the unit cell axes are
+ * the same between the two cells, and that the axes have similar directions in
+ * 3D space.  The first three tolerance values are the maximum allowable fractional
+ * differences between the a,b,c axis lengths (respectively) of the two cells.
+ * The last three tolerance values are the maximum allowable angles, in radians,
+ * between the directions of the a,b,c axes of the two cells.
+ *
+ * \p cell and \p reference must have the same centering.  Otherwise, this
+ * function always returns zero.
  *
  * \returns non-zero if the cells match.
  *
  */
-int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2,
-                                            const double ltl, const double atl)
+int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference,
+                                            const double *tols)
 {
 	double ax1, ay1, az1, bx1, by1, bz1, cx1, cy1, cz1;
 	double ax2, ay2, az2, bx2, by2, bz2, cx2, cy2, cz2;
 
-	if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
+	if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0;
 
-	cell_get_cartesian(cell1, &ax1, &ay1, &az1,
-	                          &bx1, &by1, &bz1,
-	                          &cx1, &cy1, &cz1);
+	cell_get_cartesian(cell, &ax1, &ay1, &az1,
+	                         &bx1, &by1, &bz1,
+	                         &cx1, &cy1, &cz1);
 
-	cell_get_cartesian(cell2, &ax2, &ay2, &az2,
-	                          &bx2, &by2, &bz2,
-	                          &cx2, &cy2, &cz2);
+	cell_get_cartesian(reference, &ax2, &ay2, &az2,
+	                              &bx2, &by2, &bz2,
+	                              &cx2, &cy2, &cz2);
 
-	if ( angle_between(ax1, ay1, az1, ax2, ay2, az2) > atl ) return 0;
-	if ( angle_between(bx1, by1, bz1, bx2, by2, bz2) > atl ) return 0;
-	if ( angle_between(cx1, cy1, cz1, cx2, cy2, cz2) > atl ) return 0;
+	if ( angle_between(ax1, ay1, az1, ax2, ay2, az2) > tols[3] ) return 0;
+	if ( angle_between(bx1, by1, bz1, bx2, by2, bz2) > tols[4] ) return 0;
+	if ( angle_between(cx1, cy1, cz1, cx2, cy2, cz2) > tols[5] ) return 0;
 
-	if ( moduli_check(ax1, ay1, az1, ax2, ay2, az2) > ltl ) return 0;
-	if ( moduli_check(bx1, by1, bz1, bx2, by2, bz2) > ltl ) return 0;
-	if ( moduli_check(cx1, cy1, cz1, cx2, cy2, cz2) > ltl ) return 0;
+	if ( moduli_check(ax1, ay1, az1, ax2, ay2, az2) > tols[0] ) return 0;
+	if ( moduli_check(bx1, by1, bz1, bx2, by2, bz2) > tols[1] ) return 0;
+	if ( moduli_check(cx1, cy1, cz1, cx2, cy2, cz2) > tols[2] ) return 0;
 
 	return 1;
 }
 
 
 /**
- * \param a: A UnitCell
- * \param b: Another UnitCell
- * \param ltl: Maximum allowable fractional difference in reciprocal axis lengths
- * \param atl: Maximum allowable difference in reciprocal angles (in radians)
+ * \param cell: A UnitCell
+ * \param reference: Another UnitCell
+ * \param tols: Pointer to six tolerance values (see below)
  * \param pmb: Place to store pointer to matrix
  *
  * Compare the two unit cells.  If, using any permutation of the axes, the
- * axes can be made to match in length (to within fractional difference \p ltl)
- * and the axes aligned to within \p atl radians, this function returns non-zero
- * and stores the transformation to map \p b onto \p a.
+ * axes can be made to match in length and the axes aligned in space, this
+ * function returns non-zero and stores the transformation which must be applied
+ * to \p cell at \p pmb.
+ *
+ * A "permutation" means a transformation represented by a matrix with all
+ * elements equal to +1, 0 or -1, having determinant +1 or -1.  That means that
+ * this function will find the relationship between a left-handed and a right-
+ * handed basis.
  *
  * Note that the orientations of the cells must match, not just the parameters.
  * The comparison is done after reindexing using
- * compare_cell_parameters_and_orientation().
+ * compare_cell_parameters_and_orientation().  See that function for more
+ * details.
  *
- * The cells \p a and \p b must have the same centering.  Otherwise, this function
- * always returns zero.
+ * \p cell and \p reference must have the same centering.  Otherwise, this
+ * function always returns zero.
  *
  * \returns non-zero if the cells match.
  *
  */
-int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
-                                                      double ltl, double atl,
-                                                      IntegerMatrix **pmb)
+int compare_permuted_cell_parameters_and_orientation(UnitCell *cell,
+                                                     UnitCell *reference,
+                                                     const double *tols,
+                                                     IntegerMatrix **pmb)
 {
 	IntegerMatrix *m;
 	int i[9];
 
-	if ( cell_get_centering(a) != cell_get_centering(b) ) return 0;
+	if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0;
 
 	m = intmat_new(3, 3);
 
@@ -1812,16 +1424,20 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
 		UnitCell *nc;
 		int j, k;
 		int l = 0;
+		signed int det;
 
 		for ( j=0; j<3; j++ )
 			for ( k=0; k<3; k++ )
 				intmat_set(m, j, k, i[l++]);
 
-		if ( intmat_det(m) != +1 ) continue;
+		det = intmat_det(m);
+		if ( (det != +1) && (det != -1) ) continue;
 
-		nc = cell_transform_intmat(b, m);
+		nc = cell_transform_intmat(cell, m);
 
-		if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) {
+		if ( compare_cell_parameters_and_orientation(nc, reference,
+		                                             tols) )
+		{
 			if ( pmb != NULL ) *pmb = m;
 			cell_free(nc);
 			return 1;
@@ -1928,31 +1544,18 @@ static Rational *find_candidates(double len, double *a, double *b, double *c,
 }
 
 
-static void g6_components(double *g6, double a, double b, double c,
-                                      double al, double be, double ga)
-{
-	g6[0] = a*a;
-	g6[1] = b*b;
-	g6[2] = c*c;
-	g6[3] = 2.0*b*c*cos(al);
-	g6[4] = 2.0*a*c*cos(be);
-	g6[5] = 2.0*a*b*cos(ga);
-}
-
-
-static double g6_distance(double a1, double b1, double c1,
-                          double al1, double be1, double ga1,
-                          double a2, double b2, double c2,
-                          double al2, double be2, double ga2)
+static double g6_distance(UnitCell *cell1, UnitCell *cell2)
 {
-	double g1[6], g2[6];
-	int i;
+	struct g6 g1, g2;
 	double total = 0.0;
-	g6_components(g1, a1, b1, c1, al1, be1, ga1);
-	g6_components(g2, a2, b2, c2, al2, be2, ga2);
-	for ( i=0; i<6; i++ ) {
-		total += (g1[i]-g2[i])*(g1[i]-g2[i]);
-	}
+	g1 = cell_get_G6(cell1);
+	g2 = cell_get_G6(cell2);
+	total += (g1.A-g2.A)*(g1.A-g2.A);
+	total += (g1.B-g2.B)*(g1.B-g2.B);
+	total += (g1.C-g2.C)*(g1.C-g2.C);
+	total += (g1.D-g2.D)*(g1.D-g2.D);
+	total += (g1.E-g2.E)*(g1.E-g2.E);
+	total += (g1.F-g2.F)*(g1.F-g2.F);
 	return sqrt(total);
 }
 
@@ -1960,15 +1563,15 @@ static double g6_distance(double a1, double b1, double c1,
 /**
  * \param cell_in: A UnitCell
  * \param reference_in: Another UnitCell
- * \param ltl: Maximum allowable fractional difference in direct-space axis lengths
- * \param atl: Maximum allowable difference in direct-space angles (in radians)
+ * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians)
  * \param csl: Non-zero to look for coincidence site lattice relationships
  * \param pmb: Place to store pointer to matrix
  *
  * Compare the \p cell_in with \p reference_in.  If \p cell is a derivative lattice
- * of \p reference, within fractional axis length difference \p ltl and absolute angle
- * difference \p atl (in radians), this function returns non-zero and stores the
- * transformation which needs to be applied to \p cell_in at \p pmb.
+ * of \p reference, within fractional axis length differences \p tols[0..2]
+ * and absolute angle difference \p tols[3..5] (in radians), this function returns
+ * non-zero and stores the transformation which needs to be applied to
+ * \p cell_in at \p pmb.
  *
  * Note that the tolerances will be applied to the primitive unit cell.  If
  * the reference cell is centered, a primitive unit cell will first be calculated.
@@ -1985,12 +1588,16 @@ static double g6_distance(double a1, double b1, double c1,
  * meaning that the lattice points of the transformed version of \p cell_in
  * might not coincide with lattice points of \p reference_in.
  *
+ * This function is used by CrystFEL's cell_tool program to find non-obvious
+ * relationships between crystal lattices.  For most routine comparisons, this
+ * function is probably not the one you need!
+ *
  * \returns non-zero if the cells match, zero for no match or error.
  *
  */
-int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
-                                      double ltl, double atl, int csl,
-                                      RationalMatrix **pmb)
+int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
+                                       const double *tols, int csl,
+                                       RationalMatrix **pmb)
 {
 	UnitCell *cell;
 	UnitCell *reference;
@@ -2005,7 +1612,6 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 	Rational *cand_c;
 	int ncand_a, ncand_b, ncand_c;
 	int ia, ib;
-	RationalMatrix *MCB;
 	RationalMatrix *CiAMCB;
 	double min_dist = +INFINITY;
 
@@ -2026,9 +1632,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 	                         &cv[0], &cv[1], &cv[2]);
 
 	/* Find vectors in 'cell' with lengths close to a, b and c */
-	cand_a = find_candidates(a, av, bv, cv, ltl, csl, &ncand_a);
-	cand_b = find_candidates(b, av, bv, cv, ltl, csl, &ncand_b);
-	cand_c = find_candidates(c, av, bv, cv, ltl, csl, &ncand_c);
+	cand_a = find_candidates(a, av, bv, cv, tols[0], csl, &ncand_a);
+	cand_b = find_candidates(b, av, bv, cv, tols[1], csl, &ncand_b);
+	cand_c = find_candidates(c, av, bv, cv, tols[2], csl, &ncand_c);
 
 	if ( (ncand_a==0) || (ncand_b==0) || (ncand_c==0) ) {
 		*pmb = NULL;
@@ -2040,8 +1646,6 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 	}
 
 	M = rtnl_mtx_new(3, 3);
-	MCB = rtnl_mtx_new(3, 3);
-	CiAMCB = rtnl_mtx_new(3, 3);
 	for ( ia=0; ia<ncand_a; ia++ ) {
 		for ( ib=0; ib<ncand_b; ib++ ) {
 
@@ -2049,6 +1653,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 			double at, bt, ct, alt, bet, gat;
 			double dist;
 			int ic = 0;
+			RationalMatrix *MCB;
 
 			/* Form the matrix using the first candidate for c */
 			rtnl_mtx_set(M, 0, 0, cand_a[3*ia+0]);
@@ -2065,7 +1670,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 			test = cell_transform_rational(cell, M);
 			cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
 			cell_free(test);
-			if ( fabs(gat - ga) > atl ) continue;
+			if ( fabs(gat - ga) > tols[5] ) continue;
 
 			/* Gamma OK, now look for place for c axis */
 			for ( ic=0; ic<ncand_c; ic++ ) {
@@ -2091,21 +1696,21 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 					cell_free(test);
 					continue;
 				}
-				if ( fabs(alt - al) > atl ) {
+				if ( fabs(alt - al) > tols[3] ) {
 					cell_free(test);
 					continue;
 				}
-				if ( fabs(bet - be) > atl ) {
+				if ( fabs(bet - be) > tols[4] ) {
 					cell_free(test);
 					continue;
 				}
 
-				dist = g6_distance(at, bt, ct, alt, bet, gat,
-				                   a, b, c, al, be, ga);
+				dist = g6_distance(test, reference);
 				if ( dist < min_dist ) {
 					min_dist = dist;
-					rtnl_mtx_mtxmult(M, CB, MCB);
-					rtnl_mtx_mtxmult(CiA, MCB, CiAMCB);
+					MCB = rtnlmtx_times_rtnlmtx(M, CB);
+					CiAMCB = rtnlmtx_times_rtnlmtx(CiA, MCB);
+					rtnl_mtx_free(MCB);
 				}
 
 				cell_free(test);
@@ -2115,7 +1720,6 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 	}
 
 	rtnl_mtx_free(M);
-	rtnl_mtx_free(MCB);
 	free(cand_a);
 	free(cand_b);
 	free(cand_c);
@@ -2130,3 +1734,588 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 	*pmb = CiAMCB;
 	return 1;
 }
+
+
+/* Criteria from Grosse-Kunstleve et al., Acta Cryst A60 (2004) p1-6 */
+#define GT(x,y) (y < x - eps)
+#define LT(x,y) GT(y,x)
+#define EQ(x,y) !(GT(x,y) || LT(x,y))
+#define LTE(x,y) !(GT(x,y))
+
+
+static int in_standard_presentation(struct g6 g, double eps)
+{
+	if ( GT(g.A, g.B) ) return 0;
+	if ( GT(g.B, g.C) ) return 0;
+
+	if ( EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E)) ) return 0;
+	if ( EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F)) ) return 0;
+
+	if ( ( GT(g.D, 0.0) && GT(g.E, 0.0) && GT(g.F, 0.0) )
+	  || ( !GT(g.D, 0.0) && !GT(g.E, 0.0) && !GT(g.F, 0.0) ) ) return 1;
+
+	return 0;
+}
+
+
+static int is_burger(struct g6 g, double eps)
+{
+	if ( !in_standard_presentation(g, eps) ) return 0;
+	if ( GT(fabs(g.D), g.B) ) return 0;
+	if ( GT(fabs(g.E), g.A) ) return 0;
+	if ( GT(fabs(g.F), g.A) ) return 0;
+	if ( LT(g.D + g.E + g.F + g.A + g.B, 0.0) ) return 0;
+	return 1;
+}
+
+
+static int UNUSED is_niggli(struct g6 g, double eps)
+{
+	if ( !is_burger(g, eps) ) return 0;
+
+	if ( EQ(g.D, g.B) && GT(g.F, 2.0*g.E) ) return 0;
+	if ( EQ(g.E, g.A) && GT(g.F, 2.0*g.D) ) return 0;
+	if ( EQ(g.F, g.A) && GT(g.E, 2.0*g.D) ) return 0;
+
+	if ( EQ(g.D, -g.B) && !EQ(g.F, 0.0) ) return 0;
+	if ( EQ(g.E, -g.A) && !EQ(g.F, 0.0) ) return 0;
+	if ( EQ(g.F, -g.A) && !EQ(g.E, 0.0) ) return 0;
+
+	if ( EQ(g.D + g.E + g.F + g.A + g.B, 0.0)
+	  && GT(2.0*(g.A + g.E)+g.F, 0.0) ) return 0;
+
+	return 1;
+}
+
+
+static signed int eps_sign(double v, double eps)
+{
+	return GT(v, 0.0) ? +1 : -1;
+}
+
+
+static void debug_lattice(struct g6 g, double eps, int step)
+{
+#if 0
+	STATUS("After step %i:  %e %e %e %e %e %e --", step,
+	       g.A/1e-0, g.B/1e-0, g.C/1e-0,
+	       g.D/1e-0, g.E/1e-0, g.F/1e-0);
+	if ( is_burger(g, eps) ) STATUS(" B");
+	if ( is_niggli(g, eps) ) STATUS(" N");
+	if ( eps_sign(g.D, eps)*eps_sign(g.E, eps)*eps_sign(g.F, eps) > 0 ) {
+		STATUS(" I");
+	} else {
+		STATUS(" II");
+	}
+	STATUS("\n");
+#endif
+}
+
+
+static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M)
+{
+	int i, j;
+	IntegerMatrix *tmp = intmat_times_intmat(T, M);
+	for ( i=0; i<3; i++ ) {
+		for ( j=0; j<3; j++ ) {
+			intmat_set(T, i, j, intmat_get(tmp, i, j));
+		}
+	}
+	intmat_free(tmp);
+}
+
+
+/* Cell volume from G6 components */
+static double g6_volume(struct g6 g)
+{
+	return sqrt(g.A*g.B*g.C
+	       - 0.25*(g.A*g.D*g.D + g.B*g.E*g.E + g.C*g.F*g.F - g.D*g.E*g.F));
+}
+
+
+/* NB The G6 components are passed by value, not reference.
+ * It's the caller's reponsibility to apply the matrix to the cell and
+ * re-calculate the G6 vector, if required. */
+IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
+{
+	IntegerMatrix *T;
+	IntegerMatrix *M;
+	int finished;
+	double eps;
+
+	eps = pow(g6_volume(g), 1.0/3.0) * epsrel;
+	eps = eps*eps;
+
+	T = intmat_identity(3);
+	M = intmat_new(3, 3);
+
+	debug_lattice(g, eps, 0);
+
+	do {
+
+		int done_s1s2;
+
+		do {
+			done_s1s2 = 1;
+
+			if ( GT(g.A, g.B)
+			  || (EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E))) )
+			{
+				/* Swap a and b */
+				double temp;
+				temp = g.A; g.A = g.B;  g.B = temp;
+				temp = g.D; g.D = g.E;  g.E = temp;
+
+				intmat_zero(M);
+				intmat_set(M, 1, 0, -1);
+				intmat_set(M, 0, 1, -1);
+				intmat_set(M, 2, 2, -1);
+				mult_in_place(T, M);
+
+				debug_lattice(g, eps, 1);
+
+			}
+
+			if ( GT(g.B, g.C)
+			  || (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) )
+			{
+				/* Swap b and c */
+				double temp;
+				temp = g.B; g.B = g.C;  g.C = temp;
+				temp = g.E; g.E = g.F;  g.F = temp;
+
+				intmat_zero(M);
+				intmat_set(M, 0, 0, -1);
+				intmat_set(M, 1, 2, -1);
+				intmat_set(M, 2, 1, -1);
+				mult_in_place(T, M);
+
+				debug_lattice(g, eps, 2);
+
+				/* ..."and go to 1." */
+				done_s1s2 = 0;
+			}
+
+		} while ( !done_s1s2 );
+
+		finished = 0;
+
+		/* K-G paper says g3*g4*g3, which I assume is a misprint */
+		if ( eps_sign(g.D, eps) * eps_sign(g.E, eps) * eps_sign(g.F, eps) > 0 ) {
+
+			intmat_zero(M);
+			intmat_set(M, 0, 0, LT(g.D, 0.0) ? -1 : 1);
+			intmat_set(M, 1, 1, LT(g.E, 0.0) ? -1 : 1);
+			intmat_set(M, 2, 2, LT(g.F, 0.0) ? -1 : 1);
+			mult_in_place(T, M);
+
+			g.D = fabs(g.D);
+			g.E = fabs(g.E);
+			g.F = fabs(g.F);
+
+			debug_lattice(g, eps, 3);
+
+		} else  {
+
+			int z = 4;
+			signed int ijk[3] = { 1, 1, 1 };
+
+			if ( GT(g.D, 0.0) ) {
+				ijk[0] = -1;
+			} else if ( !LT(g.D, 0.0) ) {
+				z = 0;
+			}
+			if ( GT(g.E, 0.0) ) {
+				ijk[1] = -1;
+			} else if ( !LT(g.D, 0.0) ) {
+				z = 1;
+			}
+			if ( GT(g.F, 0.0) ) {
+				ijk[2] = -1;
+			} else if ( !LT(g.D, 0.0) ) {
+				z = 2;
+			}
+			if ( ijk[0]*ijk[1]*ijk[2] < 0 ) {
+				if ( z == 4 ) {
+					ERROR("No element in reduction step 4");
+					abort();
+				}
+				ijk[z] = -1;
+			}
+			intmat_zero(M);
+			intmat_set(M, 0, 0, ijk[0]);
+			intmat_set(M, 1, 1, ijk[1]);
+			intmat_set(M, 2, 2, ijk[2]);
+			mult_in_place(T, M);
+
+			g.D = -fabs(g.D);
+			g.E = -fabs(g.E);
+			g.F = -fabs(g.F);
+
+			debug_lattice(g, eps, 4);
+
+		}
+
+		if ( GT(fabs(g.D), g.B)
+		  || (EQ(g.B, g.D) && LT(2.0*g.E, g.F))
+		  || (EQ(g.B, -g.D) && LT(g.F, 0.0)) )
+		{
+			signed int s = g.D > 0.0  ? 1 : -1;
+
+			intmat_zero(M);
+			intmat_set(M, 0, 0, 1);
+			intmat_set(M, 1, 1, 1);
+			intmat_set(M, 2, 2, 1);
+			intmat_set(M, 1, 2, -s);
+			mult_in_place(T, M);
+
+			g.C = g.B + g.C -s*g.D;
+			g.D = -2*s*g.B + g.D;
+			g.E = g.E - s*g.F;
+
+			debug_lattice(g, eps, 5);
+
+		} else if ( GT(fabs(g.E), g.A)
+		         || (EQ(g.A, g.E) && LT(2.0*g.D, g.F))
+		         || (EQ(g.A, -g.E) && LT(g.F, 0.0)) )
+		{
+			signed int s = g.E > 0.0  ? 1 : -1;
+
+			intmat_zero(M);
+			intmat_set(M, 0, 0, 1);
+			intmat_set(M, 1, 1, 1);
+			intmat_set(M, 2, 2, 1);
+			intmat_set(M, 0, 2, -s);
+			mult_in_place(T, M);
+
+			g.C = g.A + g.C -s*g.E;
+			g.D = g.D - s*g.F;
+			g.E = -2*s*g.A + g.E;
+
+			debug_lattice(g, eps, 6);
+
+		} else if ( GT(fabs(g.F), g.A)
+		         || (EQ(g.A, g.F) && LT(2.0*g.D, g.E))
+		         || (EQ(g.A, -g.F) && LT(g.E, 0.0)) )
+		{
+			signed int s = g.F > 0.0  ? 1 : -1;
+
+			intmat_zero(M);
+			intmat_set(M, 0, 0, 1);
+			intmat_set(M, 1, 1, 1);
+			intmat_set(M, 2, 2, 1);
+			intmat_set(M, 0, 1, -s);
+			mult_in_place(T, M);
+
+			g.B = g.A + g.B -s*g.F;
+			g.D = g.D - s*g.E;
+			g.F = -2*s*g.A + g.F;
+
+			debug_lattice(g, eps, 7);
+
+		} else if ( LT(g.A+g.B+g.D+g.E+g.F, 0.0)  /* not g.C */
+		         || ( (EQ(g.A+g.B+g.D+g.E+g.F, 0.0)
+		               && GT(2.0*g.A + 2.0*g.E + g.F, 0.0)) ) )
+		{
+			intmat_zero(M);
+			intmat_set(M, 0, 0, 1);
+			intmat_set(M, 1, 1, 1);
+			intmat_set(M, 2, 2, 1);
+			intmat_set(M, 1, 2, 1);
+			mult_in_place(T, M);
+
+			g.C = g.A+g.B+g.C+g.D+g.E+g.F;
+			g.D = 2.0*g.B + g.D + g.F;
+			g.E = 2.0*g.A + g.E + g.F;
+
+			debug_lattice(g, eps, 8);
+
+		} else {
+			finished = 1;
+		}
+
+	} while ( !finished );
+
+	debug_lattice(g, eps, 99);
+
+	assert(is_burger(g, eps));
+	assert(is_niggli(g, eps));
+
+	intmat_free(M);
+	return T;
+}
+
+
+static double cell_diff(UnitCell *cell, double a, double b, double c,
+                        double al, double be, double ga)
+{
+	double ta, tb, tc, tal, tbe, tga;
+	double diff = 0.0;
+	cell_get_parameters(cell, &ta, &tb, &tc, &tal, &tbe, &tga);
+	diff += fabs(a - ta);
+	diff += fabs(b - tb);
+	diff += fabs(c - tc);
+	return diff;
+}
+
+
+static int random_int(int max)
+{
+	int r;
+	int limit = RAND_MAX;
+	while ( limit % max ) limit--;
+	do {
+		r = rand();
+	} while ( r > limit );
+	return rand() % max;
+}
+
+
+static IntegerMatrix *check_permutations(UnitCell *cell_reduced, UnitCell *reference,
+                                         RationalMatrix *CiARA, IntegerMatrix *RiBCB,
+                                         const double *tols)
+{
+	IntegerMatrix *m;
+	int i[9];
+	double a, b, c, al, be, ga;
+	double min_dist = +INFINITY;
+	int s, sel;
+	IntegerMatrix *best_m[5] = {NULL, NULL, NULL, NULL, NULL};
+	int n_best = 0;
+
+	m = intmat_new(3, 3);
+	cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+
+	for ( i[0]=-1; i[0]<=+1; i[0]++ ) {
+	for ( i[1]=-1; i[1]<=+1; i[1]++ ) {
+	for ( i[2]=-1; i[2]<=+1; i[2]++ ) {
+	for ( i[3]=-1; i[3]<=+1; i[3]++ ) {
+	for ( i[4]=-1; i[4]<=+1; i[4]++ ) {
+	for ( i[5]=-1; i[5]<=+1; i[5]++ ) {
+	for ( i[6]=-1; i[6]<=+1; i[6]++ ) {
+	for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
+	for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
+
+		UnitCell *nc;
+		int j, k;
+		int l = 0;
+		signed int det;
+		UnitCell *tmp;
+
+		for ( j=0; j<3; j++ )
+			for ( k=0; k<3; k++ )
+				intmat_set(m, j, k, i[l++]);
+
+		det = intmat_det(m);
+		if ( det != +1 ) continue;
+
+		tmp = cell_transform_intmat(cell_reduced, m);
+		nc = cell_transform_intmat(tmp, RiBCB);
+		cell_free(tmp);
+
+		if ( compare_cell_parameters(nc, reference, tols) ) {
+			double dist = cell_diff(nc, a, b, c, al, be, ga);
+			if ( dist < min_dist ) {
+
+				/* If the new solution is significantly better,
+				 * dump all the previous ones */
+				for ( s=0; s<n_best; s++ ) {
+					intmat_free(best_m[s]);
+				}
+				min_dist = dist;
+				best_m[0] = intmat_copy(m);
+				n_best = 1;
+
+			} else if ( dist == min_dist ) {
+
+				if ( n_best == 5 ) {
+					ERROR("WARNING: Too many equivalent "
+					      "reindexed lattices\n");
+				} else {
+					/* If the new solution is the same as the
+					* previous one, add it to the list */
+					best_m[n_best++] = intmat_copy(m);
+				}
+
+			}
+		}
+
+		cell_free(nc);
+
+	}
+	}
+	}
+	}
+	}
+	}
+	}
+	}
+	}
+
+	if ( n_best == 0 ) return NULL;
+
+	sel = n_best;
+	if ( n_best == 1 ) {
+
+		/* If there's one solution, choose that one, of course */
+		sel = 0;
+
+	} else {
+
+		/* If one of the solutions results in an identity applied to the
+		 * original cell, choose that one */
+
+		for ( s=0; s<n_best; s++ ) {
+			RationalMatrix *tmp;
+			RationalMatrix *comb;
+			tmp = rtnlmtx_times_intmat(CiARA, best_m[s]);
+			comb = rtnlmtx_times_intmat(tmp, RiBCB);
+			if ( rtnl_mtx_is_identity(comb) ) {
+				sel = s;
+			}
+			rtnl_mtx_free(tmp);
+			rtnl_mtx_free(comb);
+		}
+
+	}
+
+	/* Still undecided?  Choose randomly, to avoid weird distributions
+	 * in the cell parameters */
+	if ( sel == n_best ) {
+		sel = random_int(n_best);
+	}
+
+	/* Free all the others */
+	for ( s=0; s<n_best; s++ ) {
+		if ( s != sel ) intmat_free(best_m[s]);
+	}
+
+	return best_m[sel];
+}
+
+
+/**
+ * \param cell_in: A UnitCell
+ * \param reference_in: Another UnitCell
+ * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians)
+ * \param pmb: Place to store pointer to matrix, or NULL if not needed
+ *
+ * Compare the \p cell_in with \p reference_in.  If they represent the same
+ * lattice, this function returns a copy of \p cell_in transformed to look
+ * similar to \p reference_in.  Otherwise, it returns NULL.
+ *
+ * If \pmb is non-NULL, the transformation which needs to be applied to
+ * \p cell_in will be stored there.
+ *
+ * Only the cell parameters will be compared.  The relative orientations are
+ * irrelevant.  The tolerances will be applied to the transformed copy of
+ * \p cell_in, i.e. the version of the input cell which looks similar to
+ * \p reference_in.  Subject to the tolerances, the cell will be chosen which
+ * has the lowest total absolute error in unit cell axis lengths.
+ *
+ * There will usually be several transformation matrices which produce exactly
+ * the same total absolute error.  If one of the matrices is an identity, that
+ * one will be used.  Otherwise, the matrix will be selected at random from the
+ * possibilities.  This avoids skewed distributions of unit cell parameters,
+ * e.g. the angles always being greater than 90 degrees.
+ *
+ * This is the right function to use for deciding if an indexing solution
+ * matches a reference cell or not.
+ *
+ * \returns A newly allocated UnitCell, or NULL.
+ *
+ */
+UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in,
+                                            UnitCell *reference_in,
+                                            const double *tols,
+                                            RationalMatrix **pmb)
+{
+	UnitCell *cell;
+	UnitCell *reference;
+	IntegerMatrix *CB;
+	RationalMatrix *CiA;
+	IntegerMatrix *RA;
+	IntegerMatrix *RB;
+	IntegerMatrix *RiB;
+	IntegerMatrix *P;
+	RationalMatrix *CiARA;
+	IntegerMatrix *RiBCB;
+	UnitCell *cell_reduced;
+	UnitCell *reference_reduced;
+	UnitCell *match;
+	struct g6 g6cell;
+	struct g6 g6ref;
+
+	//STATUS("The input cell:\n");
+	//cell_print(cell_in);
+
+	//STATUS("The reference cell:\n");
+	//cell_print(reference_in);
+
+	/* Un-center both cells */
+	reference = uncenter_cell(reference_in, &CB, NULL);
+	if ( reference == NULL ) return NULL;
+
+	cell = uncenter_cell(cell_in, NULL, &CiA);
+	if ( cell == NULL ) return NULL;
+
+	/* Calculate G6 components for both */
+	g6cell = cell_get_G6(cell);
+	g6ref = cell_get_G6(reference);
+
+	/* Convert both to reduced basis (stably) */
+	RA = reduce_g6(g6cell, 1e-5);
+	//STATUS("------------------------------------------\n");
+	RB = reduce_g6(g6ref, 1e-5);
+
+	cell_reduced = cell_transform_intmat(cell, RA);
+
+	//STATUS("Reduced cell:\n");
+	//cell_print(cell_reduced);
+	//STATUS("Reduced reference:\n");
+	reference_reduced = cell_transform_intmat(reference, RB);
+	//cell_print(reference_reduced);
+
+	/* The primitive (non-reduced) cells are no longer needed */
+	cell_free(reference);
+	cell_free(cell);
+
+	/* Within tolerance? */
+	RiB = intmat_inverse(RB);
+	RiBCB = intmat_times_intmat(RiB, CB);
+	intmat_free(RiB);
+	CiARA = rtnlmtx_times_intmat(CiA, RA);
+	P = check_permutations(cell_reduced, reference_in, CiARA, RiBCB, tols);
+	if ( P != NULL ) {
+
+		RationalMatrix *tmp;
+		RationalMatrix *comb;
+
+		//STATUS("The best permutation matrix:\n");
+		//intmat_print(P);
+
+		/* Calculate combined matrix: CB.RiB.P.RA.CiA */
+		tmp = rtnlmtx_times_intmat(CiARA, P);
+		comb = rtnlmtx_times_intmat(tmp, RiBCB);
+		rtnl_mtx_free(tmp);
+
+		match = cell_transform_rational(cell_in, comb);
+		//STATUS("Original cell transformed to look like reference:\n");
+		//cell_print(match);
+
+		if ( pmb != NULL ) *pmb = comb;
+
+	} else {
+		match = NULL;
+	}
+
+	rtnl_mtx_free(CiA);
+	intmat_free(RA);
+	intmat_free(RB);
+	intmat_free(CB);
+	intmat_free(P);
+	cell_free(reference_reduced);
+	cell_free(cell_reduced);
+
+	return match;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index be878c10..5b705247 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -59,11 +59,6 @@ extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
 extern void cell_print(UnitCell *cell);
 extern void cell_print_full(UnitCell *cell);
 
-extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
-                            const float *ltl, int reduce);
-
-extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell);
-
 extern UnitCell *load_cell_from_pdb(const char *filename);
 extern UnitCell *load_cell_from_file(const char *filename);
 extern void write_cell(UnitCell *cell, FILE *fh);
@@ -87,24 +82,26 @@ extern int forbidden_reflection(UnitCell *cell,
 
 extern double cell_get_volume(UnitCell *cell);
 
-extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
-                                   float ltl, float atl);
+extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference,
+                                   const double *tols);
 
+extern int compare_cell_parameters_and_orientation(UnitCell *cell,
+                                                   UnitCell *reference,
+                                                   const double *tols);
 
-extern int compare_cell_parameters_and_orientation(UnitCell *cell1,
-                                                   UnitCell *cell2,
-                                                   const double ltl,
-                                                   const double atl);
+extern int compare_permuted_cell_parameters_and_orientation(UnitCell *cell,
+                                                            UnitCell *reference,
+                                                            const double *tols,
+                                                            IntegerMatrix **pmb);
 
-extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a,
-                                                             UnitCell *b,
-                                                             double ltl,
-                                                             double atl,
-                                                             IntegerMatrix **pmb);
+extern int compare_derivative_cell_parameters(UnitCell *cell, UnitCell *reference,
+                                              const double *tols, int csl,
+                                              RationalMatrix **pmb);
 
-extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference,
-                                             double ltl, double atl, int csl,
-                                             RationalMatrix **pmb);
+extern UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in,
+                                                   UnitCell *reference_in,
+                                                   const double *tols,
+                                                   RationalMatrix **pmb);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 6b1d58c7..e785c1c8 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -579,6 +579,21 @@ LatticeType cell_get_lattice_type(UnitCell *cell)
 }
 
 
+struct g6 cell_get_G6(UnitCell *cell)
+{
+	double a, b, c, al, be, ga;
+	struct g6 g;
+	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+	g.A = a*a;
+	g.B = b*b;
+	g.C = c*c;
+	g.D = 2.0*b*c*cos(al);
+	g.E = 2.0*a*c*cos(be);
+	g.F = 2.0*a*b*cos(ga);
+	return g;
+}
+
+
 char cell_get_unique_axis(UnitCell *cell)
 {
 	return cell->unique_axis;
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 4bbc1be2..e12fc209 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -132,6 +132,18 @@ extern void cell_set_reciprocal(UnitCell *cell,
 extern LatticeType cell_get_lattice_type(UnitCell *cell);
 extern void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type);
 
+struct g6
+{
+	double A;
+	double B;
+	double C;
+	double D;
+	double E;
+	double F;
+};
+
+extern struct g6 cell_get_G6(UnitCell *cell);
+
 extern char cell_get_centering(UnitCell *cell);
 extern void cell_set_centering(UnitCell *cell, char centering);
 
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 6f0d1b8c..45b24afb 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -67,7 +67,7 @@ struct _indexingprivate
 {
 	IndexingFlags flags;
 	UnitCell *target_cell;
-	float tolerance[4];
+	double tolerance[6];
 
 	struct taketwo_options *ttopts;
 	struct xgandalf_options *xgandalf_opts;
@@ -315,7 +315,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 
 
 IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
-                                struct detector *det, float *ltl,
+                                struct detector *det, float *tols,
                                 IndexingFlags flags,
                                 struct taketwo_options *ttopts,
                                 struct xgandalf_options *xgandalf_opts,
@@ -439,7 +439,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 	} else {
 		ipriv->target_cell = NULL;
 	}
-	for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i];
+	for ( i=0; i<6; i++ ) ipriv->tolerance[i] = tols[i];
 
 	ipriv->ttopts = ttopts;
 	ipriv->xgandalf_opts = xgandalf_opts;
@@ -542,33 +542,28 @@ void map_all_peaks(struct image *image)
 }
 
 
+/* Return 0 for cell OK, 1 for cell incorrect */
 static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target,
-                      float *tolerance)
+                      double *tolerance)
 {
-	if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS)
-	  || (flags & INDEXING_CHECK_CELL_AXES) )
-	{
-		UnitCell *out;
-		int reduce;
+	UnitCell *out;
+	RationalMatrix *rm;
 
-		if ( flags & INDEXING_CHECK_CELL_COMBINATIONS )
-		{
-			reduce = 1;
-		} else {
-			reduce = 0;
-		}
-
-		out = match_cell(crystal_get_cell(cr),
-		                 target, 0, tolerance, reduce);
-
-		if ( out == NULL ) {
-			return 1;
-		}
+	/* Check at all? */
+	if ( ! ((flags & INDEXING_CHECK_CELL_COMBINATIONS)
+	         || (flags & INDEXING_CHECK_CELL_AXES)) ) return 0;
 
+	if ( compare_reindexed_cell_parameters(crystal_get_cell(cr), target,
+	                                       tolerance, &rm) )
+	{
+		out = cell_transform_rational(crystal_get_cell(cr), rm);
 		cell_free(crystal_get_cell(cr));
 		crystal_set_cell(cr, out);
+		rtnl_mtx_free(rm);
+		return 0;
 	}
-	return 0;
+
+	return 1;
 }
 
 
@@ -697,13 +692,17 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		for ( j=0; j<this_crystal; j++ ) {
 
 			Crystal *that_cr = image->crystals[j];
+			const double tols[] = {0.1, 0.1, 0.1,
+			                       deg2rad(5.0),
+			                       deg2rad(5.0),
+			                       deg2rad(5.0)};
 
 			/* Don't do similarity check against bad crystals */
 			if ( crystal_get_user_flag(that_cr) ) continue;
 
-			if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(cr),
-			                                                       crystal_get_cell(that_cr),
-			                                                       0.1, deg2rad(0.5), NULL) )
+			if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(cr),
+			                                                      crystal_get_cell(that_cr),
+			                                                      tols, NULL) )
 			{
 				crystal_set_user_flag(cr, 1);
 			}
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index c633a620..16aff5bc 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -218,8 +218,8 @@ signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl)
  * \returns A newly allocated \ref IntegerMatrix containing the answer, or NULL on
  * error.
  **/
-IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
-                                  const IntegerMatrix *b)
+IntegerMatrix *intmat_times_intmat(const IntegerMatrix *a,
+                                   const IntegerMatrix *b)
 {
 	unsigned int i, j;
 	IntegerMatrix *ans;
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index 5a365f0a..25c0fb99 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -74,8 +74,8 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i
 extern signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl);
 
 /* Matrix-matrix multiplication */
-extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
-                                         const IntegerMatrix *b);
+extern IntegerMatrix *intmat_times_intmat(const IntegerMatrix *a,
+                                          const IntegerMatrix *b);
 
 /* Inverse */
 extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m);
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index ce286d7a..d59da59c 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -528,15 +528,48 @@ void rtnl_mtx_print(const RationalMatrix *m)
 }
 
 
-void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
-                      RationalMatrix *ans)
+RationalMatrix *rtnlmtx_times_intmat(const RationalMatrix *A,
+                                     const IntegerMatrix *B)
 {
 	int i, j;
+	RationalMatrix *ans;
+	unsigned int B_rows, B_cols;
+
+	intmat_size(B, &B_rows, &B_cols);
+	assert(A->cols == B_rows);
+
+	ans = rtnl_mtx_new(A->rows, B_cols);
+	if ( ans == NULL ) return NULL;
+
+	for ( i=0; i<ans->rows; i++ ) {
+		for ( j=0; j<ans->cols; j++ ) {
+			int k;
+			Rational sum = rtnl_zero();
+			for ( k=0; k<A->rows; k++ ) {
+				Rational add;
+				add = rtnl_mul(rtnl_mtx_get(A, i, k),
+				               rtnl(intmat_get(B, k, j), 1));
+				sum = rtnl_add(sum, add);
+			}
+			rtnl_mtx_set(ans, i, j, sum);
+		}
+	}
+
+	return ans;
+}
+
+
+RationalMatrix *rtnlmtx_times_rtnlmtx(const RationalMatrix *A,
+                                      const RationalMatrix *B)
+{
+	int i, j;
+	RationalMatrix *ans;
 
-	assert(ans->cols == B->cols);
-	assert(ans->rows == A->rows);
 	assert(A->cols == B->rows);
 
+	ans = rtnl_mtx_new(A->rows, B->cols);
+	if ( ans == NULL ) return NULL;
+
 	for ( i=0; i<ans->rows; i++ ) {
 		for ( j=0; j<ans->cols; j++ ) {
 			int k;
@@ -550,6 +583,37 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
 			rtnl_mtx_set(ans, i, j, sum);
 		}
 	}
+
+	return ans;
+}
+
+
+RationalMatrix *intmat_times_rtnlmtx(const IntegerMatrix *A, const RationalMatrix *B)
+{
+	int i, j;
+	RationalMatrix *ans;
+	unsigned int A_rows, A_cols;
+
+	intmat_size(A, &A_rows, &A_cols);
+
+	ans = rtnl_mtx_new(A_rows, B->cols);
+	if ( ans == NULL ) return NULL;
+
+	for ( i=0; i<ans->rows; i++ ) {
+		for ( j=0; j<ans->cols; j++ ) {
+			int k;
+			Rational sum = rtnl_zero();
+			for ( k=0; k<A_rows; k++ ) {
+				Rational add;
+				add = rtnl_mul(rtnl(intmat_get(A, i, k), 1),
+				               rtnl_mtx_get(B, k, j));
+				sum = rtnl_add(sum, add);
+			}
+			rtnl_mtx_set(ans, i, j, sum);
+		}
+	}
+
+	return ans;
 }
 
 
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
index 8681bb06..880bdbeb 100644
--- a/libcrystfel/src/rational.h
+++ b/libcrystfel/src/rational.h
@@ -90,8 +90,13 @@ extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m);
 extern RationalMatrix *rtnl_mtx_identity(int rows);
 extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m);
 extern void rtnl_mtx_free(RationalMatrix *mtx);
-extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
-                             RationalMatrix *ans);
+extern RationalMatrix *rtnlmtx_times_rtnlmtx(const RationalMatrix *A,
+                                             const RationalMatrix *B);
+extern RationalMatrix *rtnlmtx_times_intmat(const RationalMatrix *A,
+                                            const IntegerMatrix *B);
+extern RationalMatrix *intmat_times_rtnlmtx(const IntegerMatrix *a,
+                                            const RationalMatrix *b);
+
 extern int transform_fractional_coords_rtnl(const RationalMatrix *P,
                                             const Rational *ivec,
                                             Rational *ans);
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 6739daaf..e834e357 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -326,7 +326,7 @@ static void expand_ops(SymOpList *s)
 				int k, nk;
 				int found;
 
-				m = intmat_intmat_mult(opi, opj);
+				m = intmat_times_intmat(opi, opj);
 				assert(m != NULL);
 
 				nk = num_ops(s);
@@ -380,13 +380,13 @@ static void transform_ops(SymOpList *s, IntegerMatrix *P)
 
 		IntegerMatrix *r, *f;
 
-		r = intmat_intmat_mult(P, s->ops[i]);
+		r = intmat_times_intmat(P, s->ops[i]);
 		if ( r == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
 		}
 
-		f = intmat_intmat_mult(r, Pi);
+		f = intmat_times_intmat(r, Pi);
 		if ( f == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
@@ -1340,7 +1340,7 @@ static int check_mult(const IntegerMatrix *ans,
 	int val;
 	IntegerMatrix *m;
 
-	m = intmat_intmat_mult(a, b);
+	m = intmat_times_intmat(a, b);
 	assert(m != NULL);
 
 	val = intmat_equals(ans, m);
@@ -1403,7 +1403,7 @@ static int order(const IntegerMatrix *m)
 
 		IntegerMatrix *anew;
 
-		anew = intmat_intmat_mult(m, a);
+		anew = intmat_times_intmat(m, a);
 		assert(anew != NULL);
 		intmat_free(a);
 		a = anew;
diff --git a/src/cell_explorer.c b/src/cell_explorer.c
index 446aec5d..812e16d5 100644
--- a/src/cell_explorer.c
+++ b/src/cell_explorer.c
@@ -2040,6 +2040,9 @@ int main(int argc, char *argv[])
 			}
 
 			w.cells[w.n_cells] = crystal_get_cell(cr);
+			if ( !right_handed(w.cells[w.n_cells]) ) {
+				ERROR("WARNING: Left-handed cell encountered\n");
+			}
 			w.indms[w.n_cells] = cur.indexed_by;
 			w.n_cells++;
 
diff --git a/src/cell_tool.c b/src/cell_tool.c
index cd457f25..83056a6f 100644
--- a/src/cell_tool.c
+++ b/src/cell_tool.c
@@ -67,16 +67,16 @@ static void show_help(const char *s)
 "     --tolerance=<tol>      Set the tolerances for cell comparison.\n"
 "                             Default: 5,1.5 (axis percentage, angle deg).\n"
 "     --highres=n            Resolution limit (Angstroms) for --rings\n"
-"     --csl                  Allow --compare to find coincidence site lattice relationships.\n"
 );
 }
 
 
 static int comparecells(UnitCell *cell, const char *comparecell,
-                        double ltl, double atl, int csl)
+                        double ltl, double atl)
 {
 	UnitCell *cell2;
 	RationalMatrix *m;
+	double tolerance[6];
 
 	cell2 = load_cell_from_file(comparecell);
 	if ( cell2 == NULL ) {
@@ -90,10 +90,17 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 	STATUS("------------------> The reference unit cell:\n");
 	cell_print(cell2);
 
-	STATUS("------------------> The comparison results:\n");
-	if ( !compare_reindexed_cell_parameters(cell, cell2, ltl, atl, csl, &m) ) {
+	tolerance[0] = ltl;
+	tolerance[1] = ltl;
+	tolerance[2] = ltl;
+	tolerance[3] = atl;
+	tolerance[4] = atl;
+	tolerance[5] = atl;
+
+	STATUS("------------------> Reindexed (strictly the same lattice):\n");
+	STATUS("Tolerances applied directly to the unit cells\n");
+	if ( !compare_reindexed_cell_parameters(cell, cell2, tolerance, &m) ) {
 		STATUS("No relationship found between lattices.\n");
-		return 0;
 	} else {
 		UnitCell *trans;
 		STATUS("Relationship found.  To become similar to the reference"
@@ -109,7 +116,46 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 
 	}
 
-	rtnl_mtx_free(m);
+	STATUS("------------------> Derivative lattice  "
+	       "(strictly the same lattice):\n");
+	STATUS("Tolerances applied to primitive versions of the unit cells\n");
+	if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 0, &m) ) {
+		STATUS("No relationship found between lattices.\n");
+	} else {
+		UnitCell *trans;
+		STATUS("Relationship found.  To become similar to the reference"
+		       " cell, the input cell should be transformed by:\n");
+		rtnl_mtx_print(m);
+		STATUS("Transformed version of input unit cell:\n");
+		trans = cell_transform_rational(cell, m);
+		cell_print(trans);
+		cell_free(trans);
+		STATUS("NB transformed cell might not really be triclinic, "
+		       "it's just that I don't (yet) know how to work out what "
+		       "it is.\n");
+		rtnl_mtx_free(m);
+	}
+
+	STATUS("------------------> Coincidence site lattice "
+	       "(not strictly the same lattice):\n");
+	STATUS("Tolerances applied to primitive versions of the unit cells\n");
+	if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 1, &m) ) {
+		STATUS("No relationship found between lattices.\n");
+		return 0;
+	} else {
+		UnitCell *trans;
+		STATUS("Relationship found.  To become similar to the reference"
+		       " cell, the input cell should be transformed by:\n");
+		rtnl_mtx_print(m);
+		STATUS("Transformed version of input unit cell:\n");
+		trans = cell_transform_rational(cell, m);
+		cell_print(trans);
+		cell_free(trans);
+		STATUS("NB transformed cell might not really be triclinic, "
+		       "it's just that I don't (yet) know how to work out what "
+		       "it is.\n");
+		rtnl_mtx_free(m);
+	}
 
 	return 0;
 }
@@ -227,6 +273,14 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 	SymOpList *ops;
 	signed int i[9];
 	const int maxorder = 3;
+	double tolerance[6];
+
+	tolerance[0] = ltl;
+	tolerance[1] = ltl;
+	tolerance[2] = ltl;
+	tolerance[3] = atl;
+	tolerance[4] = atl;
+	tolerance[5] = atl;
 
 	ops = get_pointgroup("1");
 	if ( ops == NULL ) return 1;
@@ -265,7 +319,7 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 
 		nc = cell_transform_intmat(cell, m);
 
-		if ( compare_cell_parameters(cell, nc, ltl, atl) ) {
+		if ( compare_cell_parameters(cell, nc, tolerance) ) {
 			if ( !intmat_is_identity(m) ) add_symop(ops, m);
 			STATUS("-----------------------------------------------"
 			       "-------------------------------------------\n");
@@ -426,7 +480,6 @@ int main(int argc, char *argv[])
 	float highres;
 	double rmax = 1/(2.0e-10);
 	char *trans_str = NULL;
-	int csl = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -445,7 +498,6 @@ int main(int argc, char *argv[])
 
 		{"transform",          1, NULL,                4},
 		{"highres",            1, NULL,                5},
-		{"csl",                0, &csl,                1},
 
 		{0, 0, NULL, 0}
 	};
@@ -572,8 +624,7 @@ int main(int argc, char *argv[])
 	if ( mode == CT_FINDAMBI ) return find_ambi(cell, sym, ltl, atl);
 	if ( mode == CT_UNCENTER ) return uncenter(cell, out_file);
 	if ( mode == CT_RINGS ) return all_rings(cell, sym, rmax);
-	if ( mode == CT_COMPARE ) return comparecells(cell, comparecell,
-	                                              ltl, atl, csl);
+	if ( mode == CT_COMPARE ) return comparecells(cell, comparecell, ltl, atl);
 	if ( mode == CT_TRANSFORM ) return transform(cell, trans_str, out_file);
 	if ( mode == CT_CHOICES ) return cell_choices(cell);
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 739b4969..78f26709 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -135,7 +135,7 @@ static void show_help(const char *s)
 " -p, --pdb=<file>          Unit cell file (PDB or CrystFEL unit cell format)\n"
 "                             Default: 'molecule.pdb'\n"
 "     --tolerance=<tol>     Tolerances for cell comparison\n"
-"                              Default: 5,5,5,1.5\n"
+"                              Default: 5,5,5,1.5,1.5,1.5\n"
 "     --no-check-cell       Don't check lattice parameters against input cell\n"
 "     --no-cell-combinations\n"
 "                           Don't use axis combinations when checking cell\n"
@@ -297,10 +297,12 @@ int main(int argc, char *argv[])
 	iargs.noisefilter = 0;
 	iargs.median_filter = 0;
 	iargs.satcorr = 1;
-	iargs.tols[0] = 5.0;
-	iargs.tols[1] = 5.0;
-	iargs.tols[2] = 5.0;
+	iargs.tols[0] = 0.05;
+	iargs.tols[1] = 0.05;
+	iargs.tols[2] = 0.05;
 	iargs.tols[3] = 1.5;
+	iargs.tols[4] = 1.5;
+	iargs.tols[5] = 1.5;
 	iargs.threshold = 800.0;
 	iargs.min_sq_gradient = 100000.0;
 	iargs.min_snr = 5.0;
@@ -1056,14 +1058,29 @@ int main(int argc, char *argv[])
 	/* Parse unit cell tolerance */
 	if ( toler != NULL ) {
 		int ttt;
-		ttt = sscanf(toler, "%f,%f,%f,%f",
-		             &iargs.tols[0], &iargs.tols[1],
-		             &iargs.tols[2], &iargs.tols[3]);
-		if ( ttt != 4 ) {
-			ERROR("Invalid parameters for '--tolerance'\n");
-			return 1;
+		ttt = sscanf(toler, "%f,%f,%f,%f,%f,%f",
+		             &iargs.tols[0], &iargs.tols[1], &iargs.tols[2],
+		             &iargs.tols[3], &iargs.tols[4], &iargs.tols[5]);
+		if ( ttt != 6 ) {
+			ttt = sscanf(toler, "%f,%f,%f,%f",
+			             &iargs.tols[0], &iargs.tols[1],
+			             &iargs.tols[2], &iargs.tols[3]);
+			if ( ttt != 4 ) {
+				ERROR("Invalid parameters for '--tolerance'\n");
+				return 1;
+			}
+			iargs.tols[4] = iargs.tols[3];
+			iargs.tols[5] = iargs.tols[3];
 		}
 		free(toler);
+
+		/* Percent to fraction */
+		iargs.tols[0] /= 100.0;
+		iargs.tols[1] /= 100.0;
+		iargs.tols[2] /= 100.0;
+		iargs.tols[3] = deg2rad(iargs.tols[3]);
+		iargs.tols[4] = deg2rad(iargs.tols[4]);
+		iargs.tols[5] = deg2rad(iargs.tols[5]);
 	}
 
 	/* Parse integration radii */
diff --git a/src/process_image.h b/src/process_image.h
index 8d6682a6..f8b0f4b1 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -74,7 +74,7 @@ struct index_args
 	struct detector *det;
 	IndexingPrivate *ipriv;
 	int peaks;                /* Peak detection method */
-	float tols[4];
+	float tols[6];
 	struct beam_params *beam;
 	char *hdf5_peak_path;
 	int half_pixel_shift;
diff --git a/src/whirligig.c b/src/whirligig.c
index cc2b531e..fd94a0f2 100644
--- a/src/whirligig.c
+++ b/src/whirligig.c
@@ -297,6 +297,8 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2,
 	IntegerMatrix *m;
 	struct image *i1;
 	struct image *i2;
+	const double tols[] = {0.1, 0.1, 0.1,
+	                       deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)};
 
 	assert(n1 >= 0);
 	assert(n2 >= 0);
@@ -309,9 +311,9 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2,
 	for ( i=0; i<i1->n_crystals; i++ ) {
 	for ( j=0; j<i2->n_crystals; j++ ) {
 
-		if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]),
-		                                                       crystal_get_cell(i2->crystals[j]),
-		                                                       0.1, deg2rad(5.0), &m) )
+		if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]),
+		                                                      crystal_get_cell(i2->crystals[j]),
+		                                                      tols, &m) )
 		{
 			if ( !crystal_used(win, n1, i)
 			  && !crystal_used(win, n2, j) )
@@ -365,6 +367,8 @@ static int try_join(struct window *win, int sn)
 	Crystal *cr;
 	UnitCell *ref;
 	const int sp = win->join_ptr - 1;
+	const double tols[] = {0.1, 0.1, 0.1,
+	                       deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)};
 
 	/* Get the appropriately transformed cell from the last crystal in this
 	 * series */
@@ -374,9 +378,9 @@ static int try_join(struct window *win, int sn)
 	for ( j=0; j<win->img[win->join_ptr].n_crystals; j++ ) {
 		Crystal *cr2;
 		cr2 = win->img[win->join_ptr].crystals[j];
-		if ( compare_reindexed_cell_parameters_and_orientation(ref, crystal_get_cell(cr2),
-		                                                       0.1, deg2rad(5.0),
-		                                                       &win->mat[sn][win->join_ptr]) )
+		if ( compare_permuted_cell_parameters_and_orientation(ref, crystal_get_cell(cr2),
+		                                                      tols,
+		                                                      &win->mat[sn][win->join_ptr]) )
 		{
 			win->ser[sn][win->join_ptr] = j;
 			cell_free(ref);
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 192913ce..bd9959ae 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -92,3 +92,8 @@ add_executable(spectrum_check spectrum_check.c)
 target_include_directories(spectrum_check PRIVATE ${COMMON_INCLUDES})
 target_link_libraries(spectrum_check ${COMMON_LIBRARIES})
 add_test(spectrum_check spectrum_check)
+
+add_executable(cellcompare_check cellcompare_check.c)
+target_include_directories(cellcompare_check PRIVATE ${COMMON_INCLUDES})
+target_link_libraries(cellcompare_check ${COMMON_LIBRARIES})
+add_test(cellcompare_check cellcompare_check)
diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c
new file mode 100644
index 00000000..f674a7f2
--- /dev/null
+++ b/tests/cellcompare_check.c
@@ -0,0 +1,403 @@
+/*
+ * cellcompare_check.c
+ *
+ * Check that unit cell comparison works
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 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/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <gsl/gsl_randist.h>
+
+#include <cell.h>
+#include <cell-utils.h>
+#include <utils.h>
+
+
+static void complain(UnitCell *cell, UnitCell *cref, const char *t, const char *c)
+{
+	STATUS("These cells should %sbe the same%s:\n", t, c);
+	STATUS("Transformed: ----------------------------\n");
+	cell_print_full(cell);
+	STATUS("Original: ---------------------------\n");
+	cell_print_full(cref);
+}
+
+
+static RationalMatrix *random_derivative(gsl_rng *rng)
+{
+	int i, j;
+	RationalMatrix *tr;
+
+	tr = rtnl_mtx_new(3, 3);
+	if ( tr == NULL ) return NULL;
+
+	do {
+		for ( i=0; i<3; i++ ) {
+			for ( j=0; j<3; j++ ) {
+				/* 0..6 inclusive -> -3..3 inclusive */
+				signed int a = gsl_rng_uniform_int(rng, 7) - 3;
+				/* 0..2 inclusive */
+				signed int b = gsl_rng_uniform_int(rng, 3);
+				if ( b == 0 ) {
+					a = 0;
+					b = 1;
+				}
+				rtnl_mtx_set(tr, i, j, rtnl(a, b));
+			}
+		}
+
+	} while ( rtnl_cmp(rtnl_mtx_det(tr), rtnl_zero()) == 0 );
+
+	return tr;
+}
+
+
+static IntegerMatrix *random_permutation(gsl_rng *rng)
+{
+	int i;
+	IntegerMatrix *tr;
+	int v[] = {0, 1, 2};
+
+	tr = intmat_new(3, 3);
+	if ( tr == NULL ) return NULL;
+
+	gsl_ran_shuffle(rng, v, 3, sizeof(int));
+	for ( i=0; i<3; i++ ) {
+		intmat_set(tr, i, v[i], 1);
+	}
+
+	return tr;
+}
+
+
+static int check_ccp(UnitCell *cell, UnitCell *cref, double *tols,
+                     int should_match)
+{
+	const char *a;
+	const char *b;
+
+	a = should_match ? "" : "NOT ";
+	b = " with compare_cell_parameters";
+
+	if ( compare_cell_parameters(cell, cref, tols) != should_match )
+	{
+		complain(cell, cref, a, b);
+		return 1;
+	}
+	return 0;
+}
+
+
+static int check_ccpao(UnitCell *cell, UnitCell *cref, double *tols,
+                       int should_match)
+{
+	const char *a;
+	const char *b;
+
+	a = should_match ? "" : "NOT ";
+	b = " with compare_cell_parameters_and_orientation";
+
+	if ( compare_cell_parameters_and_orientation(cell, cref,
+	                                             tols) != should_match )
+	{
+		complain(cell, cref, a, b);
+		return 1;
+	}
+	return 0;
+}
+
+
+static int check_cpcpao(UnitCell *cell, UnitCell *cref, double *tols,
+                        int should_match)
+{
+	IntegerMatrix *m = NULL;
+	const char *a;
+	const char *b;
+
+	a = should_match ? "" : "NOT ";
+	b = " with compare_permuted_cell_parameters_and_orientation";
+
+	if ( compare_permuted_cell_parameters_and_orientation(cell, cref, tols, &m)
+	                        != should_match )
+	{
+		complain(cell, cref, a, b);
+		intmat_free(m);
+		return 1;
+	}
+	intmat_free(m);
+	return 0;
+}
+
+
+static int check_cdcp(UnitCell *cell, UnitCell *cref, double *tols,
+                      RationalMatrix *tr, int should_match)
+{
+	RationalMatrix *m = NULL;
+	const char *a;
+	const char *b;
+
+	a = should_match ? "" : "NOT ";
+	b = " with compare_derivative_cell_parameters";
+
+	if ( compare_derivative_cell_parameters(cref, cell, tols, 1, &m) != should_match )
+	{
+		complain(cell, cref, a, b);
+		STATUS("The transformation matrix to create the cell was:\n");
+		rtnl_mtx_print(tr);
+		STATUS("Cell comparison says do this to the reference to "
+		       "create the cell:\n");
+		rtnl_mtx_print(m);
+		rtnl_mtx_free(m);
+		return 1;
+	}
+	rtnl_mtx_free(m);
+	return 0;
+}
+
+
+static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols,
+                      RationalMatrix *tr, int should_match)
+{
+	RationalMatrix *m = NULL;
+	const char *a;
+	const char *b;
+	UnitCell *match;
+
+	a = should_match ? "" : "NOT ";
+	b = " with compare_reindexed_cell_parameters";
+
+	match = compare_reindexed_cell_parameters(cell, cref, tols, &m);
+
+	if ( (match != NULL) != should_match )
+	{
+		complain(cell, cref, a, b);
+		STATUS("The transformation matrix to create the cell was:\n");
+		rtnl_mtx_print(tr);
+		STATUS("Cell comparison says do this to the reference to "
+		       "create the cell:\n");
+		rtnl_mtx_print(m);
+		rtnl_mtx_free(m);
+		return 1;
+	}
+	rtnl_mtx_free(m);
+	return 0;
+}
+
+
+static void yaro_test()
+{
+	UnitCell *cell;
+	UnitCell *reference;
+	UnitCell *cmatch;
+	//float tols[] = {5, 5, 5, 1.5};
+	double dtols[] = {0.05, 0.05, 0.05, deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)};
+
+	cell = cell_new_from_parameters(63.24e-10, 63.94e-10, 64.61e-10,
+	                                deg2rad(89.98), deg2rad(89.82), deg2rad(119.87));
+	cell_set_unique_axis(cell, 'c');
+	cell_set_lattice_type(cell, L_HEXAGONAL);
+	cell_set_centering(cell, 'P');
+
+	reference = cell_new_from_parameters(64.7e-10, 64.7e-10, 65.2e-10,
+	                                     deg2rad(90.0), deg2rad(90.0), deg2rad(120.0));
+	cell_set_unique_axis(reference, 'c');
+	cell_set_lattice_type(reference, L_HEXAGONAL);
+	cell_set_centering(reference, 'P');
+
+	STATUS("The cell:\n");
+	cell_print(cell);
+	STATUS("The reference:\n");
+	cell_print(reference);
+	//cmatch = match_cell(cell, reference, 0, tols, 1);
+	//STATUS("The match:\n");
+	//cell_print(cmatch);
+	//cell_free(cmatch);
+
+	RationalMatrix *m = NULL;
+	cmatch = compare_reindexed_cell_parameters(cell, reference, dtols, &m);
+	STATUS("The new match:\n")
+	cell_print(cmatch);
+	STATUS("The matrix:\n");
+	rtnl_mtx_print(m);
+	cell_free(cmatch);
+	rtnl_mtx_free(m);
+
+	cell_free(cell);
+	cell_free(reference);
+}
+
+
+extern IntegerMatrix *reduce_g6(struct g6 g, double epsrel);
+
+int main(int argc, char *argv[])
+{
+	UnitCell *cell, *cref;
+	gsl_rng *rng;
+	int i;
+	const int ntrial = 10;
+	double tols[] = { 0.01, 0.01, 0.01,
+	                  deg2rad(1.0), deg2rad(1.0), deg2rad(1.0) };
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
+	yaro_test();
+
+	cref = cell_new_from_parameters(3e-10, 5.196e-10, 2e-10,
+	                                deg2rad(103.9166666),
+	                                deg2rad(109.4666666),
+	                                deg2rad(134.8833333));
+	cell_set_centering(cref, 'P');
+	if ( cref == NULL ) return 1;
+
+	/* Just rotate cell */
+	for ( i=0; i<ntrial; i++ ) {
+
+		cell = cell_rotate(cref, random_quaternion(rng));
+		if ( cell == NULL ) return 1;
+
+		if ( check_ccp(cell, cref, tols, 1) ) return 1;
+		if ( check_ccpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cpcpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1;
+		if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1;
+
+		cell_free(cell);
+		progress_bar(i+1, ntrial, "Plain rotation");
+	}
+
+	/* Permute axes but don't rotate */
+	for ( i=0; i<ntrial; i++ ) {
+
+		IntegerMatrix *tr;
+
+		tr = random_permutation(rng);
+		cell = cell_transform_intmat(cref, tr);
+
+		if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1;
+		if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1;
+		if ( check_cpcpao(cell, cref, tols, 1) ) return 1;
+		if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1;
+		if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1;
+
+		cell_free(cell);
+		intmat_free(tr);
+		progress_bar(i+1, ntrial, "Axis permutation");
+	}
+
+	/* Rotate cell and permute axes */
+	for ( i=0; i<ntrial; i++ ) {
+
+		IntegerMatrix *tr;
+		UnitCell *cell2;
+
+		cell2 = cell_rotate(cref, random_quaternion(rng));
+		if ( cell2 == NULL ) return 1;
+
+		tr = random_permutation(rng);
+		cell = cell_transform_intmat(cell2, tr);
+		cell_free(cell2);
+
+		if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1;
+		if ( check_ccpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cpcpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1;
+		if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1;
+
+		cell_free(cell);
+		intmat_free(tr);
+		progress_bar(i+1, ntrial, "Rotation with axis permutation");
+	}
+
+	/* Derivative lattice */
+	for ( i=0; i<ntrial; i++ ) {
+
+		RationalMatrix *tr;
+
+		cell = NULL;
+		tr = NULL;
+		do {
+			cell_free(cell);
+			rtnl_mtx_free(tr);
+			tr = random_derivative(rng);
+			cell = cell_transform_rational(cref, tr);
+		} while ( (cell_get_centering(cell) == '?')
+		       || (cell_get_centering(cell) == 'H' ) );
+		/* H centering is no good because it needs a unique axis to
+		 * be specified in order for uncentering in c_r_c_p to work.
+		 * cell_transform_rational doesn't set the unique axis (yet?) */
+
+		if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1;
+		if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1;
+		if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1;
+		if ( check_cdcp(cell, cref, tols, tr, 1) ) return 1;
+		/* check_crcp: Sometimes yes, hard to tell */
+
+		cell_free(cell);
+		rtnl_mtx_free(tr);
+		progress_bar(i+1, ntrial, "Derivative lattice");
+	}
+
+	/* Derivative lattice and rotate */
+	for ( i=0; i<ntrial; i++ ) {
+
+		RationalMatrix *tr;
+		UnitCell *cell2;
+
+		cell2 = cell_rotate(cref, random_quaternion(rng));
+		if ( cell2 == NULL ) return 1;
+
+		cell = NULL;
+		tr = NULL;
+		do {
+			cell_free(cell);
+			rtnl_mtx_free(tr);
+			tr = random_derivative(rng);
+			cell = cell_transform_rational(cell2, tr);
+		} while ( (cell_get_centering(cell) == '?')
+		       || (cell_get_centering(cell) == 'H' ) );  /* See above */
+		cell_free(cell2);
+
+		if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1;
+		if ( check_ccpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cpcpao(cell, cref, tols, 0) ) return 1;
+		if ( check_cdcp(cell, cref, tols, tr, 1) ) return 1;
+		/* check_crcp: Sometimes yes, hard to tell */
+
+		cell_free(cell);
+		rtnl_mtx_free(tr);
+		progress_bar(i+1, ntrial, "Derivative lattice with rotation");
+	}
+
+	cell_free(cref);
+	gsl_rng_free(rng);
+
+	return 0;
+}
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 33f44340..8c8f7fa0 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -396,7 +396,7 @@ int main(int argc, char *argv[])
 	intmat_set_all_3x3(part2,  0,0,1,
 	                           0,1,0,
 	                          -1,0,0);
-	tfn = intmat_intmat_mult(part1, part2);
+	tfn = intmat_times_intmat(part1, part2);
 	fail += check_identity(cell, tfn);
 	intmat_free(part1);
 	intmat_free(part2);