Use IntegerMatrix for all unit cell transformations

Get rid of UnitCellTransformation, a thin wrapper which didn't do
anything.

12 files changed, 406 insertions, 574 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 39e893cf..507d4cc2 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -361,157 +361,114 @@ int bravais_lattice(UnitCell *cell)
 }
 
 
-static UnitCellTransformation *uncentering_transformation(UnitCell *in,
-                                                          char *new_centering,
-                                                          LatticeType *new_latt)
+/* Given a centered cell @in, return the integer transformation matrix which
+ * turns a primitive cell into @in. Set new_centering and new_latt to the
+ * centering and lattice type of the primitive cell (usually aP, sometimes rR,
+ * rarely mP) */
+static IntegerMatrix *centering_transformation(UnitCell *in,
+                                               char *new_centering,
+                                               LatticeType *new_latt,
+                                               char *new_ua)
 {
-	UnitCellTransformation *t;
-	const double OT = 1.0/3.0;
-	const double TT = 2.0/3.0;
-	const double H = 0.5;
 	LatticeType lt;
 	char ua, cen;
+	IntegerMatrix *t = NULL;
 
 	lt = cell_get_lattice_type(in);
 	ua = cell_get_unique_axis(in);
 	cen = cell_get_centering(in);
 
-	t = tfn_identity();
-	if ( t == NULL ) return NULL;
-
-	if ( ua == 'a' ) {
-		tfn_combine(t, tfn_vector(0,1,0),
-		               tfn_vector(0,0,1),
-		               tfn_vector(1,0,0));
-		if ( lt == L_MONOCLINIC ) {
-			assert(cen != 'A');
-			switch ( cen ) {
-				case 'B' : cen = 'A'; break;
-				case 'C' : cen = 'B'; break;
-				case 'I' : cen = 'I'; break;
-			}
-		}
-	}
-
-	if ( ua == 'b' ) {
-		tfn_combine(t, tfn_vector(0,0,1),
-		               tfn_vector(1,0,0),
-		               tfn_vector(0,1,0));
-		if ( lt == L_MONOCLINIC ) {
-			assert(cen != 'B');
-			switch ( cen ) {
-				case 'C' : cen = 'A'; break;
-				case 'A' : cen = 'B'; break;
-				case 'I' : cen = 'I'; break;
-			}
-		}
-	}
-
-	switch ( cen ) {
-
-		case 'P' :
+	if ( (cen=='P') || (cen=='R') ) {
+		*new_centering = cen;
 		*new_latt = lt;
-		*new_centering = 'P';
-		break;
-
-		case 'R' :
-		*new_latt = L_RHOMBOHEDRAL;
-		*new_centering = 'R';
-		break;
+		*new_ua = ua;
+		t = intmat_identity(3);
+	}
 
-		case 'I' :
-		tfn_combine(t, tfn_vector(-H,H,H),
-		               tfn_vector(H,-H,H),
-		               tfn_vector(H,H,-H));
+	if ( cen == 'I' ) {
+		t = intmat_create_3x3(0, 1, 1,
+		                      1, 0, 1,
+		                      1, 1, 0);
 		if ( lt == L_CUBIC ) {
 			*new_latt = L_RHOMBOHEDRAL;
 			*new_centering = 'R';
+			*new_ua = '*';
 		} else {
-			/* Tetragonal or orthorhombic */
 			*new_latt = L_TRICLINIC;
 			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		break;
+	}
 
-		case 'F' :
-		tfn_combine(t, tfn_vector(0,H,H),
-		               tfn_vector(H,0,H),
-		               tfn_vector(H,H,0));
+	if ( cen == 'F' ) {
+		t = intmat_create_3x3(-1,  1,  1,
+		                       1, -1,  1,
+		                       1,  1, -1);
 		if ( lt == L_CUBIC ) {
 			*new_latt = L_RHOMBOHEDRAL;
 			*new_centering = 'R';
+			*new_ua = '*';
 		} else {
-			assert(lt == L_ORTHORHOMBIC);
 			*new_latt = L_TRICLINIC;
 			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		break;
+	}
 
-		case 'A' :
-		tfn_combine(t, tfn_vector( 1, 0, 0),
-		               tfn_vector( 0, H, H),
-		               tfn_vector( 0,-H, H));
+	if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) {
+		/* Obverse setting */
+		t = intmat_create_3x3(1, -1,  0,
+		                      0,  1, -1,
+		                      1,  1,  1);
+		assert(lt == L_HEXAGONAL);
+		assert(ua == 'c');
+		*new_latt = L_RHOMBOHEDRAL;
+		*new_centering = 'R';
+		*new_ua = '*';
+	}
+
+	if ( cen == 'A' ) {
+		t = intmat_create_3x3(1,  0,  0,
+		                      0,  1, -1,
+		                      0,  1,  1);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'a';
 		} else {
 			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		*new_centering = 'P';
-		break;
+	}
 
-		case 'B' :
-		tfn_combine(t, tfn_vector( H, 0, H),
-		               tfn_vector( 0, 1, 0),
-		               tfn_vector(-H, 0, H));
+	if ( cen == 'B' ) {
+		t = intmat_create_3x3(1,  0, -1,
+		                      0,  1,  0,
+		                      1,  0,  1);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'b';
 		} else {
 			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		*new_centering = 'P';
-		break;
+	}
 
-		case 'C' :
-		tfn_combine(t, tfn_vector( H, H, 0),
-		               tfn_vector(-H, H, 0),
-		               tfn_vector( 0, 0, 1));
+	if ( cen == 'C' ) {
+		t = intmat_create_3x3(1, -1,  0,
+		                      1,  1,  0,
+		                      0,  0,  1);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'c';
 		} else {
 			*new_latt = L_TRICLINIC;
-		}
-		*new_centering = 'P';
-		break;
-
-		case 'H' :
-		/* Obverse setting */
-		tfn_combine(t, tfn_vector(TT,OT,OT),
-		               tfn_vector(-OT,OT,OT),
-		               tfn_vector(-OT,-TT,OT));
-		assert(lt == L_HEXAGONAL);
-		*new_latt = L_RHOMBOHEDRAL;
-		*new_centering = 'R';
-		break;
-
-		default :
-		ERROR("Invalid centering '%c'\n", cell_get_centering(in));
-		return NULL;
-
-	}
-
-	/* Reverse the axis permutation, but only if this was not an H->R
-	 * transformation */
-	if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) {
-		if ( ua == 'a' ) {
-			tfn_combine(t, tfn_vector(0,0,1),
-				       tfn_vector(1,0,0),
-				       tfn_vector(0,1,0));
-		}
-
-		if ( ua == 'b' ) {
-			tfn_combine(t, tfn_vector(0,1,0),
-				       tfn_vector(0,0,1),
-				       tfn_vector(1,0,0));
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
 	}
 
@@ -532,32 +489,28 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
  * setting.
  *
  */
-UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
+UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **t)
 {
-	UnitCellTransformation *tt;
+	IntegerMatrix *tt;
 	char new_centering;
 	LatticeType new_latt;
+	char new_ua;
 	UnitCell *out;
 
-	if ( !bravais_lattice(in) ) {
-		ERROR("Cannot uncenter: not a Bravais lattice.\n");
-		cell_print(in);
-		return NULL;
-	}
-
-	tt = uncentering_transformation(in, &new_centering, &new_latt);
+	tt = centering_transformation(in, &new_centering, &new_latt, &new_ua);
 	if ( tt == NULL ) return NULL;
 
-	out = cell_transform(in, tt);
+	out = cell_transform_inverse(in, tt);
 	if ( out == NULL ) return NULL;
 
 	cell_set_lattice_type(out, new_latt);
 	cell_set_centering(out, new_centering);
+	cell_set_unique_axis(out, new_ua);
 
 	if ( t != NULL ) {
 		*t = tt;
 	} else {
-		tfn_free(tt);
+		intmat_free(tt);
 	}
 
 	return out;
@@ -621,11 +574,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	float angtol = deg2rad(tols[3]);
 	UnitCell *cell;
 	UnitCell *template;
-	UnitCellTransformation *uncentering;
+	IntegerMatrix *centering;
 	UnitCell *new_cell_trans;
 
 	/* "Un-center" the template unit cell to make the comparison easier */
-	template = uncenter_cell(template_in, &uncentering);
+	template = uncenter_cell(template_in, &centering);
 	if ( template == NULL ) return NULL;
 
 	/* The candidate cell is also uncentered, because it might be centered
@@ -639,7 +592,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		ERROR("Couldn't get reciprocal cell for template.\n");
 		cell_free(template);
 		cell_free(cell);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 
@@ -661,7 +614,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		ERROR("Couldn't get reciprocal cell.\n");
 		cell_free(template);
 		cell_free(cell);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 
@@ -823,7 +776,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	/* Reverse the de-centering transformation */
 	if ( new_cell != NULL ) {
 
-		new_cell_trans = cell_transform_inverse(new_cell, uncentering);
+		new_cell_trans = cell_transform(new_cell, centering);
 		cell_free(new_cell);
 		cell_set_lattice_type(new_cell_trans,
 		                      cell_get_lattice_type(template_in));
@@ -833,13 +786,13 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		                     cell_get_unique_axis(template_in));
 
 		cell_free(template);
-		tfn_free(uncentering);
+		intmat_free(centering);
 
 		return new_cell_trans;
 
 	} else {
 		cell_free(template);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 }
@@ -862,7 +815,7 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
 	int have_real_c;
 	UnitCell *cell;
 	UnitCell *template;
-	UnitCellTransformation *to_given_cell;
+	IntegerMatrix *to_given_cell;
 	UnitCell *new_cell;
 	UnitCell *new_cell_trans;
 
@@ -1455,49 +1408,6 @@ 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
@@ -1658,7 +1568,7 @@ int forbidden_reflection(UnitCell *cell,
 	cen = cell_get_centering(cell);
 
 	/* Reflection conditions here must match the transformation matrices
-	 * in uncentering_transformation().  tests/centering_check verifies
+	 * in centering_transformation().  tests/centering_check verifies
 	 * this (amongst other things). */
 
 	if ( cen == 'P' ) return 0;
@@ -1850,7 +1760,6 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
 	for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
 	for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
 
-		UnitCellTransformation *tfn;
 		UnitCell *nc;
 		int j, k;
 		int l = 0;
@@ -1861,17 +1770,14 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
 
 		if ( intmat_det(m) != +1 ) continue;
 
-		tfn = tfn_from_intmat(m);
-		nc = cell_transform(b, tfn);
+		nc = cell_transform(b, m);
 
 		if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) {
 			if ( pmb != NULL ) *pmb = m;
-			tfn_free(tfn);
 			cell_free(nc);
 			return 1;
 		}
 
-		tfn_free(tfn);
 		cell_free(nc);
 
 	}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 47145481..1a4f3bcc 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -49,7 +49,6 @@ 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);
 extern void cell_print_full(UnitCell *cell);
@@ -67,7 +66,7 @@ extern int cell_is_sensible(UnitCell *cell);
 
 extern int validate_cell(UnitCell *cell);
 
-extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t);
+extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **m);
 
 extern int bravais_lattice(UnitCell *cell);
 
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index cc18b49d..242094f6 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -45,6 +45,7 @@
 #include "cell.h"
 #include "utils.h"
 #include "image.h"
+#include "integer_matrix.h"
 
 
 /**
@@ -600,70 +601,88 @@ const char *cell_rep(UnitCell *cell)
 }
 
 
-struct _unitcelltransformation
+UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
 {
-	gsl_matrix *m;
-};
-
-
-/**
- * tfn_inverse:
- * @t: A %UnitCellTransformation.
- *
- * Calculates the inverse of @t.  That is, if you apply cell_transform() to a
- * %UnitCell using @t, and then apply cell_transform() to the result using
- * tfn_inverse(@t) instead of @t, you will recover the same lattice vectors
- * (but note that the lattice type, centering and unique axis information will
- * be lost).
- *
- * Returns: The inverse of @t.
- *
- */
-UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
-{
-	int s;
-	gsl_matrix *m;
-	gsl_matrix *inv;
-	gsl_permutation *perm;
-	UnitCellTransformation *out;
-
-	m = gsl_matrix_alloc(3, 3);
-	if ( m == NULL ) return NULL;
+	gsl_matrix *c;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	gsl_matrix *res;
+	UnitCell *out;
 
-	out = tfn_identity();
-	if ( out == NULL ) {
-		gsl_matrix_free(m);
-		return NULL;
-	}
+	cell_get_cartesian(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, 0, 1, asy);
+	gsl_matrix_set(c, 0, 2, asz);
+	gsl_matrix_set(c, 1, 0, bsx);
+	gsl_matrix_set(c, 1, 1, bsy);
+	gsl_matrix_set(c, 1, 2, bsz);
+	gsl_matrix_set(c, 2, 0, csx);
+	gsl_matrix_set(c, 2, 1, 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_cartesian(out, gsl_matrix_get(res, 0, 0),
+	                        gsl_matrix_get(res, 0, 1),
+	                        gsl_matrix_get(res, 0, 2),
+	                        gsl_matrix_get(res, 1, 0),
+	                        gsl_matrix_get(res, 1, 1),
+	                        gsl_matrix_get(res, 1, 2),
+	                        gsl_matrix_get(res, 2, 0),
+	                        gsl_matrix_get(res, 2, 1),
+	                        gsl_matrix_get(res, 2, 2));
+
+	gsl_matrix_free(res);
+	gsl_matrix_free(c);
+	return out;
+}
 
-	gsl_matrix_memcpy(m, t->m);
 
-	perm = gsl_permutation_alloc(m->size1);
-	if ( perm == NULL ) {
-		ERROR("Couldn't allocate permutation\n");
-		return NULL;
-	}
-	inv = gsl_matrix_alloc(m->size1, m->size2);
-	if ( inv == NULL ) {
-		ERROR("Couldn't allocate inverse\n");
-		gsl_permutation_free(perm);
-		return NULL;
-	}
-	if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
-		ERROR("Couldn't decompose matrix\n");
-		gsl_permutation_free(perm);
-		return NULL;
-	}
-	if ( gsl_linalg_LU_invert(m, perm, inv)  ) {
-		ERROR("Couldn't invert transformation matrix\n");
-		gsl_permutation_free(perm);
-		return NULL;
-	}
-	gsl_permutation_free(perm);
+UnitCell *cell_transform_gsl_reciprocal(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;
 
-	gsl_matrix_free(out->m);
-	gsl_matrix_free(m);
-	out->m = inv;
+	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, 0, 1, asy);
+	gsl_matrix_set(c, 0, 2, asz);
+	gsl_matrix_set(c, 1, 0, bsx);
+	gsl_matrix_set(c, 1, 1, bsy);
+	gsl_matrix_set(c, 1, 2, bsz);
+	gsl_matrix_set(c, 2, 0, csx);
+	gsl_matrix_set(c, 2, 1, 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, 0, 1),
+	                         gsl_matrix_get(res, 0, 2),
+	                         gsl_matrix_get(res, 1, 0),
+	                         gsl_matrix_get(res, 1, 1),
+	                         gsl_matrix_get(res, 1, 2),
+	                         gsl_matrix_get(res, 2, 0),
+	                         gsl_matrix_get(res, 2, 1),
+	                         gsl_matrix_get(res, 2, 2));
+
+	gsl_matrix_free(res);
+	gsl_matrix_free(c);
 	return out;
 }
 
@@ -671,63 +690,40 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
 /**
  * cell_transform:
  * @cell: A %UnitCell.
- * @t: A %UnitCellTransformation.
+ * @t: An %IntegerMatrix.
  *
- * Applies @t to @cell.  Note that the lattice type, centering and unique axis
- * information will not be preserved.
+ * Applies @t to @cell.
  *
  * Returns: Transformed copy of @cell.
  *
  */
-UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
+UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m)
 {
 	UnitCell *out;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	gsl_matrix *m;
-	gsl_matrix *a;
+	gsl_matrix *tm;
 
-	if ( t == NULL ) return NULL;
-
-	out = cell_new_from_cell(cell);
-	if ( out == NULL ) return NULL;
-
-	if ( cell_get_cartesian(out, &ax, &ay, &az,
-	                             &bx, &by, &bz,
-	                             &cx, &cy, &cz) ) return NULL;
+	if ( m == NULL ) return NULL;
 
-	m = gsl_matrix_alloc(3,3);
-	a = gsl_matrix_calloc(3,3);
-	if ( (m == NULL) || (a == NULL) ) {
-		cell_free(out);
+	tm = gsl_matrix_alloc(3,3);
+	if ( tm == NULL ) {
 		return NULL;
 	}
 
-	gsl_matrix_set(m, 0, 0, ax);
-	gsl_matrix_set(m, 0, 1, ay);
-	gsl_matrix_set(m, 0, 2, az);
-	gsl_matrix_set(m, 1, 0, bx);
-	gsl_matrix_set(m, 1, 1, by);
-	gsl_matrix_set(m, 1, 2, bz);
-	gsl_matrix_set(m, 2, 0, cx);
-	gsl_matrix_set(m, 2, 1, cy);
-	gsl_matrix_set(m, 2, 2, cz);
+	gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0));
+	gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1));
+	gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2));
+	gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0));
+	gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1));
+	gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2));
+	gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0));
+	gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1));
+	gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2));
 
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 0.0, a);
+	out = cell_transform_gsl_direct(cell, tm);
 
-	cell_set_cartesian(out, gsl_matrix_get(a, 0, 0),
-	                        gsl_matrix_get(a, 0, 1),
-	                        gsl_matrix_get(a, 0, 2),
-	                        gsl_matrix_get(a, 1, 0),
-	                        gsl_matrix_get(a, 1, 1),
-	                        gsl_matrix_get(a, 1, 2),
-	                        gsl_matrix_get(a, 2, 0),
-	                        gsl_matrix_get(a, 2, 1),
-	                        gsl_matrix_get(a, 2, 2));
+	/* FIXME: Update centering, unique axis, lattice type */
 
-	gsl_matrix_free(a);
-	gsl_matrix_free(m);
+	gsl_matrix_free(tm);
 
 	return out;
 }
@@ -736,197 +732,67 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
 /**
  * cell_transform_inverse:
  * @cell: A %UnitCell.
- * @t: A %UnitCellTransformation.
+ * @m: An %IntegerMatrix
  *
- * Applies the inverse of @t to @cell.
+ * Applies the inverse of @m to @cell.
  *
  * Returns: Transformed copy of @cell.
  *
  */
-UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t)
+UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m)
 {
-	UnitCellTransformation *inv;
 	UnitCell *out;
+	gsl_matrix *tm;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+	int s;
 
-	inv = tfn_inverse(t);
-	out = cell_transform(cell, inv);
-	tfn_free(inv);
-	return out;
-}
-
-
-/**
- * tfn_identity:
- *
- * Returns: A %UnitCellTransformation corresponding to an identity operation.
- *
- */
-UnitCellTransformation *tfn_identity()
-{
-	UnitCellTransformation *tfn;
-
-	tfn = calloc(1, sizeof(UnitCellTransformation));
-	if ( tfn == NULL ) return NULL;
+	if ( m == NULL ) return NULL;
 
-	tfn->m = gsl_matrix_alloc(3, 3);
-	if ( tfn->m == NULL ) {
-		free(tfn);
+	tm = gsl_matrix_alloc(3,3);
+	if ( tm == NULL ) {
 		return NULL;
 	}
 
-	gsl_matrix_set_identity(tfn->m);
-
-	return tfn;
-}
-
-
-
-/**
- * tfn_from_intmat:
- * @m: An %IntegerMatrix
- *
- * Returns: A %UnitCellTransformation corresponding to @m.
- *
- */
-UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m)
-{
-	UnitCellTransformation *tfn;
-	int i, j;
-
-	tfn = tfn_identity();
-	if ( tfn == NULL ) return NULL;
-
-	for ( i=0; i<3; i++ ) {
-	for ( j=0; j<3; j++ ) {
-		gsl_matrix_set(tfn->m, i, j, intmat_get(m, i, j));
-	}
-	}
-
-	return tfn;
-}
-
-
-/**
- * tfn_combine:
- * @t: A %UnitCellTransformation
- * @na: Pointer to three doubles representing naa, nab, nac
- * @nb: Pointer to three doubles representing nba, nbb, nbc
- * @nc: Pointer to three doubles representing nca, ncb, ncc
- *
- * Updates @t such that it represents its previous transformation followed by
- * a new transformation, corresponding to letting a = naa*a + nab*b + nac*c.
- * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c.
- *
- */
-void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
-{
-	gsl_matrix *a;
-	gsl_matrix *n;
-
-	n = gsl_matrix_alloc(3, 3);
-	a = gsl_matrix_calloc(3, 3);
-	if ( (n == NULL) || (a == NULL) ) {
-		return;
+	gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0));
+	gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1));
+	gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2));
+	gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0));
+	gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1));
+	gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2));
+	gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0));
+	gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1));
+	gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2));
+
+	perm = gsl_permutation_alloc(3);
+	if ( perm == NULL ) {
+		ERROR("Couldn't allocate permutation\n");
+		return NULL;
 	}
-	gsl_matrix_set(n, 0, 0, na[0]);
-	gsl_matrix_set(n, 0, 1, na[1]);
-	gsl_matrix_set(n, 0, 2, na[2]);
-	gsl_matrix_set(n, 1, 0, nb[0]);
-	gsl_matrix_set(n, 1, 1, nb[1]);
-	gsl_matrix_set(n, 1, 2, nb[2]);
-	gsl_matrix_set(n, 2, 0, nc[0]);
-	gsl_matrix_set(n, 2, 1, nc[1]);
-	gsl_matrix_set(n, 2, 2, nc[2]);
-
-	free(na);
-	free(nb);
-	free(nc);
-
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 0.0, a);
-
-	gsl_matrix_free(t->m);
-	t->m = a;
-	gsl_matrix_free(n);
-}
-
-
-/**
- * tfn_vector:
- * @a: Amount of "a" to include in new vector
- * @b: Amount of "b" to include in new vector
- * @c: Amount of "c" to include in new vector
- *
- * This is a convenience function to use when sending vectors to tfn_combine():
- *	tfn_combine(tfn, tfn_vector(1,0,0),
- *	                 tfn_vector(0,2,0),
- *	                 tfn_vector(0,0,1));
- *
- */
-double *tfn_vector(double a, double b, double c)
-{
-	double *vec = malloc(3*sizeof(double));
-	if ( vec == NULL ) return NULL;
-	vec[0] = a;  vec[1] = b;  vec[2] = c;
-	return vec;
-}
-
-
-/**
- * tfn_print:
- * @t: A %UnitCellTransformation
- *
- * Prints information about @t to stderr.
- *
- */
-void tfn_print(UnitCellTransformation *t)
-{
-	gsl_permutation *perm;
-	gsl_matrix *lu;
-	int s;
-
-	STATUS("New a = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 0, 0),
-	                                         gsl_matrix_get(t->m, 0, 1),
-	                                         gsl_matrix_get(t->m, 0, 2));
-	STATUS("New b = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 1, 0),
-	                                         gsl_matrix_get(t->m, 1, 1),
-	                                         gsl_matrix_get(t->m, 1, 2));
-	STATUS("New c = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 2, 0),
-	                                         gsl_matrix_get(t->m, 2, 1),
-	                                         gsl_matrix_get(t->m, 2, 2));
-	lu = gsl_matrix_alloc(3, 3);
-	if ( lu == NULL ) {
-		ERROR("Couldn't allocate LU decomposition.\n");
-		return;
+	inv = gsl_matrix_alloc(3, 3);
+	if ( inv == NULL ) {
+		ERROR("Couldn't allocate inverse\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
-
-	gsl_matrix_memcpy(lu, t->m);
-
-	perm = gsl_permutation_alloc(t->m->size1);
-	if ( perm == NULL ) {
-		ERROR("Couldn't allocate permutation.\n");
-		gsl_matrix_free(lu);
-		return;
+	if ( gsl_linalg_LU_decomp(tm, perm, &s) ) {
+		ERROR("Couldn't decompose matrix\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
-	if ( gsl_linalg_LU_decomp(lu, perm, &s) ) {
-		ERROR("LU decomposition failed.\n");
+	if ( gsl_linalg_LU_invert(tm, perm, inv)  ) {
+		ERROR("Couldn't invert transformation matrix\n");
 		gsl_permutation_free(perm);
-		gsl_matrix_free(lu);
-		return;
+		return NULL;
 	}
+	gsl_permutation_free(perm);
 
-	STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s));
-}
+	out = cell_transform_gsl_direct(cell, inv);
 
+	/* FIXME: Update centering, unique axis, lattice type */
 
-/**
- * tfn_free:
- * @t: A %UnitCellTransformation
- *
- * Frees all resources associated with @t.
- *
- */
-void tfn_free(UnitCellTransformation *t)
-{
-	gsl_matrix_free(t->m);
-	free(t);
+	gsl_matrix_free(tm);
+	gsl_matrix_free(inv);
+
+	return out;
 }
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 39e6a1ed..e6ccbdf7 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -91,14 +91,6 @@ typedef enum
 typedef struct _unitcell UnitCell;
 
 
-/**
- * UnitCellTransformation:
- *
- * This opaque data structure represents a tranformation of a unit cell, such
- * as a rotation or a centering operation.
- **/
-typedef struct _unitcelltransformation UnitCellTransformation;
-
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -156,18 +148,10 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
 
 extern const char *cell_rep(UnitCell *cell);
 
-extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t);
-extern UnitCell *cell_transform_inverse(UnitCell *cell,
-                                        UnitCellTransformation *t);
-
-extern UnitCellTransformation *tfn_identity(void);
-extern UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m);
-extern void tfn_combine(UnitCellTransformation *t,
-                        double *na, double *nb, double *nc);
-extern void tfn_print(UnitCellTransformation *t);
-extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t);
-extern double *tfn_vector(double a, double b, double c);
-extern void tfn_free(UnitCellTransformation *t);
+extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m);
+extern UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m);
+extern UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m);
+extern UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index c31072b4..1e616e43 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -342,6 +342,9 @@ static IntegerMatrix *intmat_cofactors(const IntegerMatrix *m)
  *
  * Calculates the inverse of @m.  Inefficiently.
  *
+ * Works only if the inverse of the matrix is also an integer matrix,
+ * i.e. if the determinant of @m is +/- 1.
+ *
  * Returns: the inverse of @m, or NULL on error.
  **/
 IntegerMatrix *intmat_inverse(const IntegerMatrix *m)
@@ -532,3 +535,39 @@ IntegerMatrix *intmat_identity(int size)
 
 	return m;
 }
+
+
+/**
+ * intmat_set_all_3x3
+ * @m: An %IntegerMatrix
+ *
+ * Returns: an identity %IntegerMatrix with side length @size, or NULL on error.
+ *
+ */
+void intmat_set_all_3x3(IntegerMatrix *m,
+                        signed int m11, signed int m12, signed int m13,
+                        signed int m21, signed int m22, signed int m23,
+                        signed int m31, signed int m32, signed int m33)
+{
+	intmat_set(m, 0, 0, m11);
+	intmat_set(m, 0, 1, m12);
+	intmat_set(m, 0, 2, m13);
+
+	intmat_set(m, 1, 0, m21);
+	intmat_set(m, 1, 1, m22);
+	intmat_set(m, 1, 2, m23);
+
+	intmat_set(m, 2, 0, m31);
+	intmat_set(m, 2, 1, m32);
+	intmat_set(m, 2, 2, m33);
+}
+
+
+IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+                                 signed int m21, signed int m22, signed int m23,
+                                 signed int m31, signed int m32, signed int m33)
+{
+	IntegerMatrix *t = intmat_new(3, 3);
+	intmat_set_all_3x3(t, m11, m12, m13, m21, m22, m23, m31, m32, m33);
+	return t;
+}
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index 6616b5e8..c4983aef 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -56,6 +56,15 @@ extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j,
 extern signed int intmat_get(const IntegerMatrix *m,
                              unsigned int i, unsigned int j);
 
+extern void intmat_set_all_3x3(IntegerMatrix *m,
+                               signed int m11, signed int m12, signed int m13,
+                               signed int m21, signed int m22, signed int m23,
+                               signed int m31, signed int m32, signed int m33);
+
+extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+                                        signed int m21, signed int m22, signed int m23,
+                                        signed int m31, signed int m32, signed int m33);
+
 /* Matrix-(int)vector multiplication */
 extern signed int *intmat_intvec_mult(const IntegerMatrix *m,
                                       const signed int *vec);
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 203e5a05..37a29722 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -98,6 +98,7 @@
 #include <assert.h>
 #include <time.h>
 
+#include "cell.h"
 #include "cell-utils.h"
 #include "index.h"
 #include "taketwo.h"
@@ -2058,7 +2059,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 	tp->numPrevs++;
 
 	/* Prepare the solution for CrystFEL friendliness */
-	result = transform_cell_gsl(cell, solution);
+	result = cell_transform_gsl_reciprocal(cell, solution);
 	cleanup_taketwo_cell(&ttCell);
 
 	return result;
diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c
index 5f31df33..e6618560 100644
--- a/libcrystfel/src/xgandalf.c
+++ b/libcrystfel/src/xgandalf.c
@@ -46,7 +46,7 @@ struct xgandalf_private_data {
 	IndexingMethod indm;
 	UnitCell *cellTemplate;
 	Lattice_t sampleRealLattice_A;   //same as cellTemplate
-	UnitCellTransformation *uncenteringTransformation;
+	IntegerMatrix *centeringTransformation;
 	LatticeTransform_t latticeReductionTransform;
 };
 
@@ -114,7 +114,7 @@ int run_xgandalf(struct image *image, void *ipriv)
 		if(xgandalf_private_data->cellTemplate != NULL){
 			restoreCell(uc, &xgandalf_private_data->latticeReductionTransform);
 
-			UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation);
+			UnitCell *new_cell_trans = cell_transform(uc, xgandalf_private_data->centeringTransformation);
 			cell_free(uc);
 			uc = new_cell_trans;
 
@@ -147,7 +147,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
 	allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A));
 	xgandalf_private_data->indm = *indm;
 	xgandalf_private_data->cellTemplate = NULL;
-	xgandalf_private_data->uncenteringTransformation = NULL;
+	xgandalf_private_data->centeringTransformation = NULL;
 
 	float tolerance = xgandalf_opts->tolerance;
 	samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch;
@@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
 
 		xgandalf_private_data->cellTemplate = cell;
 
-		UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation);
+		UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation);
 
 		UnitCell *uc = cell_new_from_cell(primitiveCell);
 		reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform);
@@ -250,8 +250,8 @@ void xgandalf_cleanup(void *pp)
 
 	freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A);
 	IndexerPlain_delete(xgandalf_private_data->indexer);
-	if(xgandalf_private_data->uncenteringTransformation != NULL){
-		tfn_free(xgandalf_private_data->uncenteringTransformation);
+	if(xgandalf_private_data->centeringTransformation != NULL){
+		intmat_free(xgandalf_private_data->centeringTransformation);
 	}
 	free(xgandalf_private_data);
 }
diff --git a/src/cell_tool.c b/src/cell_tool.c
index e8e60f30..12a00192 100644
--- a/src/cell_tool.c
+++ b/src/cell_tool.c
@@ -105,7 +105,6 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 	for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) {
 	for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) {
 
-		UnitCellTransformation *tfn;
 		UnitCell *nc;
 		IntegerMatrix *m;
 		int j, k;
@@ -120,8 +119,7 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 
 		if ( intmat_det(m) < 1 ) continue;
 
-		tfn = tfn_from_intmat(m);
-		nc = cell_transform(cell, tfn);
+		nc = cell_transform(cell, m);
 
 		if ( compare_cell_parameters(cell2, nc, ltl, atl) ) {
 			STATUS("-----------------------------------------------"
@@ -131,7 +129,6 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 		}
 
 		intmat_free(m);
-		tfn_free(tfn);
 		cell_free(nc);
 
 	}
@@ -283,7 +280,6 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 	for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) {
 	for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) {
 
-		UnitCellTransformation *tfn;
 		UnitCell *nc;
 		IntegerMatrix *m;
 		int j, k;
@@ -298,8 +294,7 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 
 		if ( intmat_det(m) != +1 ) continue;
 
-		tfn = tfn_from_intmat(m);
-		nc = cell_transform(cell, tfn);
+		nc = cell_transform(cell, m);
 
 		if ( compare_cell_parameters(cell, nc, ltl, atl) ) {
 			if ( !intmat_is_identity(m) ) add_symop(ops, m);
@@ -311,7 +306,6 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 			intmat_free(m);
 		}
 
-		tfn_free(tfn);
 		cell_free(nc);
 
 	}
@@ -345,15 +339,15 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
 static int uncenter(UnitCell *cell, const char *out_file)
 {
 	UnitCell *cnew;
-	UnitCellTransformation *trans;
+	IntegerMatrix *m;
 
-	cnew = uncenter_cell(cell, &trans);
+	cnew = uncenter_cell(cell, &m);
 
 	STATUS("------------------> The primitive unit cell:\n");
 	cell_print(cnew);
 
-	STATUS("------------------> The decentering transformation:\n");
-	tfn_print(trans);
+	STATUS("------------------> The centering transformation:\n");
+	intmat_print(m);
 
 	if ( out_file != NULL ) {
 		FILE *fh = fopen(out_file, "w");
diff --git a/src/whirligig.c b/src/whirligig.c
index 54ade40a..7d50141a 100644
--- a/src/whirligig.c
+++ b/src/whirligig.c
@@ -364,15 +364,12 @@ static int try_join(struct window *win, int sn)
 	int j;
 	Crystal *cr;
 	UnitCell *ref;
-	UnitCellTransformation *tfn;
 	const int sp = win->join_ptr - 1;
 
 	/* Get the appropriately transformed cell from the last crystal in this
 	 * series */
-	tfn = tfn_from_intmat(win->mat[sn][sp]);
 	cr = win->img[sp].crystals[win->ser[sn][sp]];
-	ref = cell_transform(crystal_get_cell(cr), tfn);
-	tfn_free(tfn);
+	ref = cell_transform(crystal_get_cell(cr), win->mat[sn][sp]);
 
 	for ( j=0; j<win->img[win->join_ptr].n_crystals; j++ ) {
 		Crystal *cr2;
diff --git a/tests/centering_check.c b/tests/centering_check.c
index 887311dd..774af66b 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -61,7 +61,7 @@ static int check_centering(double a, double b, double c,
 {
 	UnitCell *cell, *cref;
 	UnitCell *n;
-	UnitCellTransformation *t;
+	IntegerMatrix *t;
 	int fail = 0;
 	int i;
 	double asx, asy, asz;
@@ -88,12 +88,17 @@ static int check_centering(double a, double b, double c,
 	check_cell(cell, "Input");
 	n = uncenter_cell(cell, &t);
 	if ( n != NULL ) {
-		STATUS("Transformation was:\n");
-		tfn_print(t);
+		STATUS("The back transformation is:\n");
+		intmat_print(t);
 		if ( check_cell(n, "Output") ) fail = 1;
-		if ( !fail ) cell_print(n);
+		if ( intmat_is_identity(t) && ((cen=='P') || (cen=='R')) ) {
+			STATUS("Identity, as expected.\n");
+		} else {
+			cell_print(n);
+		}
 	} else {
-		fail = 1;
+		STATUS("*************  Could not uncenter.\n");
+		return 1;
 	}
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
@@ -299,10 +304,6 @@ int main(int argc, char *argv[])
 	fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
 	                        L_HEXAGONAL, 'H', 'c', rng);
 
-	/* Cubic P */
-	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_CUBIC, 'P', '*', rng);
-
 	/* Cubic I */
 	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
 	                        L_CUBIC, 'I', '*', rng);
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 2e7e7378..2eb01aa6 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -37,6 +37,7 @@
 
 #include <cell.h>
 #include <cell-utils.h>
+#include <symmetry.h>
 
 
 #define MAX_REFLS (10*1024)
@@ -134,17 +135,44 @@ static int compare_rvecs(struct rvec *a, int na, struct rvec *b, int nb)
 }
 
 
-static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
+static int check_same_reflections(UnitCell *cell, UnitCell *cnew)
+{
+	struct rvec *vecs;
+	struct rvec *tvecs;
+	int na, nb;
+
+	/* Check that the two cells predict the same reflections */
+	vecs = all_refls(cell, 1e9, &na);
+	tvecs = all_refls(cnew, 1e9, &nb);
+	if ( compare_rvecs(vecs, na, tvecs, nb)
+	    || compare_rvecs(tvecs, nb, vecs, na) )
+	{
+		ERROR("******* Transformed cell didn't predict the same reflections\n");
+		//printf("---\n");
+		//for ( i=0; i<na; i++ ) {
+		//	printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w);
+		//}
+		//printf("---\n");
+		//for ( i=0; i<nb; i++ ) {
+		//	printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w);
+		//}
+		return 1;
+	} else {
+		STATUS("The cells predict the same reflections.\n");
+	}
+	free(vecs);
+	free(tvecs);
+	return 0;
+}
+
+
+static int check_transformation(UnitCell *cell, IntegerMatrix *tfn,
                                 int pred_test, UnitCell *ct)
 {
 	UnitCell *cnew, *cback;
-	UnitCellTransformation *inv;
 	double a[9], b[9];
 	int i;
 	int fail = 0;
-	struct rvec *vecs;
-	struct rvec *tvecs;
-	int na, nb;
 
 	STATUS("-----------------------\n");
 	if ( ct == NULL ) {
@@ -154,32 +182,15 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
 	}
 	cback = cell_transform_inverse(cnew, tfn);
 
+	STATUS("----> Before transformation:\n");
 	cell_print(cell);
-	tfn_print(tfn);
+	STATUS("----> The transformation matrix:\n");
+	intmat_print(tfn);
+	STATUS("----> After transformation:\n");
 	cell_print(cnew);
 
 	if ( pred_test ) {
-		/* Check that the two cells predict the same reflections */
-		vecs = all_refls(cell, 1e9, &na);
-		tvecs = all_refls(cnew, 1e9, &nb);
-		if ( compare_rvecs(vecs, na, tvecs, nb)
-		  || compare_rvecs(tvecs, nb, vecs, na) )
-		{
-			ERROR("Transformed cell didn't predict the same reflections\n");
-			//printf("---\n");
-			//for ( i=0; i<na; i++ ) {
-			//	printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w);
-			//}
-			//printf("---\n");
-			//for ( i=0; i<nb; i++ ) {
-			//	printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w);
-			//}
-			return 1;
-		} else {
-			STATUS("The cells predict the same reflections.\n");
-		}
-		free(vecs);
-		free(tvecs);
+		check_same_reflections(cell, cnew);
 	} else {
 		STATUS("Cells not expected to predict the same reflections.\n");
 	}
@@ -193,18 +204,17 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
 	                          &b[6], &b[7], &b[8]);
 	for ( i=0; i<9; i++ ) {
 		if ( !tolerance(a[i], b[i]) ) {
+			//STATUS("%e %e\n", a[i], b[i]);
 			fail = 1;
-			STATUS("%e %e\n", a[i], b[i]);
 		}
 	}
 
 	if ( fail ) {
-		ERROR("Original cell not recovered after transformation:\n");
-		cell_print(cell);
-		tfn_print(tfn);
-		inv = tfn_inverse(tfn);
-		tfn_print(inv);
+		ERROR("******* Original cell not recovered after transformation:\n");
+		STATUS("----> After transformation and transformation back:\n");
 		cell_print(cback);
+	} else {
+		STATUS("The original cell was recovered after inverse transform.\n");
 	}
 
 	return fail;
@@ -214,22 +224,48 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
 static int check_uncentering(UnitCell *cell)
 {
 	UnitCell *ct;
-	UnitCellTransformation *tr;
+	IntegerMatrix *tr;
+	int fail = 0;
+
+	STATUS("-----------------------\n");
+
+	STATUS("----> Before transformation:\n");
+	cell_print(cell);
 
 	ct = uncenter_cell(cell, &tr);
-	return check_transformation(cell, tr, 1, ct);
+	if ( ct == NULL ) return 1;
+
+	STATUS("----> The primitive unit cell:\n");
+	cell_print(ct);
+
+	STATUS("----> The matrix to put the centering back:\n");
+	intmat_print(tr);
+
+	fail = check_same_reflections(cell, ct);
+	cell_free(ct);
+
+	return fail;
 }
 
 
-static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
+static int check_identity(UnitCell *cell, IntegerMatrix *tfn)
 {
 	UnitCell *cnew;
 	double a[9], b[9];
 	int i;
 	int fail = 0;
 
+	STATUS("-----------------------\n");
+
 	cnew = cell_transform(cell, tfn);
 
+	STATUS("----> Before identity transformation:\n");
+	cell_print(cell);
+	STATUS("----> The identity transformation matrix:\n");
+	intmat_print(tfn);
+	STATUS("----> After identity transformation:\n");
+	cell_print(cnew);
+
 	cell_get_cartesian(cell, &a[0], &a[1], &a[2],
 	                         &a[3], &a[4], &a[5],
 	                         &a[6], &a[7], &a[8]);
@@ -239,14 +275,14 @@ static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
 	for ( i=0; i<9; i++ ) {
 		if ( !within_tolerance(a[i], b[i], 0.1) ) {
 			fail = 1;
-			STATUS("%e %e\n", a[i], b[i]);
+			//STATUS("%e %e\n", a[i], b[i]);
 		}
 	}
 
 	if ( fail ) {
-		ERROR("Original cell not recovered after transformation:\n");
+		ERROR("Original cell not recovered after identity transformation:\n");
 		cell_print(cell);
-		tfn_print(tfn);
+		intmat_print(tfn);
 		cell_print(cnew);
 	}
 
@@ -258,7 +294,8 @@ int main(int argc, char *argv[])
 {
 	int fail = 0;
 	UnitCell *cell, *cref;
-	UnitCellTransformation *tfn;
+	IntegerMatrix *tfn;
+	IntegerMatrix *part1, *part2;
 	gsl_rng *rng;
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
@@ -273,53 +310,52 @@ int main(int argc, char *argv[])
 	if ( cell == NULL ) return 1;
 	cell_free(cref);
 
+	tfn = intmat_identity(3);
+
 	/* Permutation of axes */
-	tfn = tfn_identity();
 	if ( tfn == NULL ) return 1;
-	tfn_combine(tfn, tfn_vector(0,1,0),
-	                 tfn_vector(0,0,1),
-	                 tfn_vector(1,0,0));
+	intmat_set_all_3x3(tfn, 0,1,0,
+	                        0,0,1,
+	                        1,0,0);
 	fail += check_transformation(cell, tfn, 1, NULL);
-	tfn_free(tfn);
 
 	/* Doubling of cell in one direction */
-	tfn = tfn_identity();
 	if ( tfn == NULL ) return 1;
-	tfn_combine(tfn, tfn_vector(2,0,0),
-	                 tfn_vector(0,1,0),
-	                 tfn_vector(0,0,1));
+	intmat_set_all_3x3(tfn, 2,0,0,
+	                        0,1,0,
+	                        0,0,1);
 	fail += check_transformation(cell, tfn, 0, NULL);
-	tfn_free(tfn);
 
-	/* Diagonal supercell */
-	tfn = tfn_identity();
+	/* Shearing */
 	if ( tfn == NULL ) return 1;
-	tfn_combine(tfn, tfn_vector(1,1,0),
-	                 tfn_vector(0,1,0),
-	                 tfn_vector(0,0,1));
+	intmat_set_all_3x3(tfn, 1,1,0,
+	                        0,1,0,
+	                        0,0,1);
 	fail += check_transformation(cell, tfn, 1, NULL);
-	tfn_free(tfn);
 
 	/* Crazy */
-	tfn = tfn_identity();
 	if ( tfn == NULL ) return 1;
-	tfn_combine(tfn, tfn_vector(1,1,0),
-	                 tfn_vector(0,1,0),
-	                 tfn_vector(0,1,1));
+	intmat_set_all_3x3(tfn, 1,1,0,
+	                        0,1,0,
+	                        0,1,1);
 	fail += check_transformation(cell, tfn, 0, NULL);
-	tfn_free(tfn);
 
 	/* Identity in two parts */
-	tfn = tfn_identity();
+	part1 = intmat_identity(3);
+	part2 = intmat_identity(3);
 	if ( tfn == NULL ) return 1;
-	tfn_combine(tfn, tfn_vector(0,0,1),
-	                 tfn_vector(0,1,0),
-	                 tfn_vector(-1,0,0));
-	tfn_combine(tfn,  tfn_vector(0,0,-1),
-	                  tfn_vector(0,1,0),
-	                  tfn_vector(1,0,0));
+	intmat_set_all_3x3(part1, 0,0,1,
+	                        0,1,0,
+	                        -1,0,0);
+	intmat_set_all_3x3(part2,  0,0,-1,
+	                         0,1,0,
+	                         1,0,0);
+	tfn = intmat_intmat_mult(part1, part2);
 	fail += check_identity(cell, tfn);
-	tfn_free(tfn);
+	intmat_free(part1);
+	intmat_free(part2);
+
+	intmat_free(tfn);
 
 	/* Check some uncentering transformations */
 	cref = cell_new_from_parameters(50e-10, 50e-10, 50e-10,