WIP on cell transformations

5 files changed, 239 insertions, 42 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 508b70e5..776526aa 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -83,7 +83,7 @@ UnitCell *cell_rotate(UnitCell *in, struct quaternion quat)
 }
 
 
-static const char *str_lattice(LatticeType l)
+const char *str_lattice(LatticeType l)
 {
 	switch ( l )
 	{
@@ -253,17 +253,19 @@ int bravais_lattice(UnitCell *cell)
 		if ( lattice == L_HEXAGONAL ) return 1;
 		return 0;
 
+		case 'R' :
+		if ( lattice == L_RHOMBOHEDRAL ) return 1;
+		return 0;
+
 		default :
 		return 0;
 	}
 }
 
 
-/* Turn any cell into a primitive one, e.g. for comparison purposes */
-UnitCell *uncenter_cell(UnitCell *in)
+static UnitCellTransformation *uncentering_transformation(UnitCell *in)
 {
-	UnitCell *cell;
-
+	UnitCellTransformation *t;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
@@ -276,34 +278,29 @@ UnitCell *uncenter_cell(UnitCell *in)
 	LatticeType lt;
 	char ua, cen;
 
-	if ( !bravais_lattice(in) ) {
-		ERROR("Cannot uncenter: not a Bravais lattice.\n");
-		return NULL;
-	}
-
 	cell_get_cartesian(in, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 	lt = cell_get_lattice_type(in);
 	ua = cell_get_unique_axis(in);
 	cen = cell_get_centering(in);
 
-	if ( ua == 'a' ) {
+	if ( (ua == 'a') || (cen == 'A') ) {
 		double tx, ty, tz;
 		tx = ax;  ty = ay;  tz = az;
 		ax = cx;  ay = cy;  az = cz;
-		cx = tx;  cy = ty;  cz = tz;
+		cx = -tx;  cy = -ty;  cz = -tz;
 		if ( cen == 'A' ) cen = 'C';
+		ua = 'c';
 	}
 
-	if ( ua == 'b' ) {
+	if ( (ua == 'b') || (cen == 'B') ) {
 		double tx, ty, tz;
 		tx = bx;  ty = by;  tz = bz;
-		ax = cx;  ay = cy;  az = cz;
-		cx = tx;  cy = ty;  cz = tz;
+		bx = cx;  by = cy;  bz = cz;
+		cx = -tx;  cy = -ty;  cz = -tz;
 		if ( cen == 'B' ) cen = 'C';
+		ua = 'c';
 	}
 
-	cell = cell_new();
-
 	switch ( cen ) {
 
 		case 'P' :
@@ -392,9 +389,39 @@ UnitCell *uncenter_cell(UnitCell *in)
 
 	}
 
-	cell_set_cartesian(cell, nax, nay, naz, nbx, nby, nbz, ncx, ncy, ncz);
 
-	return cell;
+}
+
+
+/**
+ * uncenter_cell:
+ * @in: A %UnitCell
+ * @t: Location at which to store the transformation which was used.
+ *
+ * Turns any cell into a primitive one, e.g. for comparison purposes.  The
+ * transformation which was used is stored at @t, which can be NULL if the
+ * transformation is not required.
+ *
+ * Returns: a primitive version of @in in a conventional (unique axis c)
+ * setting.
+ *
+ */
+UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
+{
+	UnitCell *cell;
+	UnitCellTransformation *tt;
+
+	if ( !bravais_lattice(in) ) {
+		ERROR("Cannot uncenter: not a Bravais lattice.\n");
+		return NULL;
+	}
+
+	tt = uncentering_transformation(in);
+	if ( tt == NULL ) return NULL;
+
+	if ( t != NULL ) *t = tt;
+
+	return cell_transform(in, tt);
 }
 
 
@@ -455,9 +482,15 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	float angtol = deg2rad(tols[3]);
 	UnitCell *cell;
 	UnitCell *template;
+	UnitCellTransformation *to_given_cell;
+	UnitCell *new_cell_trans;
 
-	cell = uncenter_cell(cell_in);
-	template = uncenter_cell(template_in);
+	/* "Un-center" the template unit cell to make the comparison easier */
+	template = uncenter_cell(template_in, &to_given_cell);
+
+	/* The candidate cell is also uncentered, because it might be centered
+	 * if it came from (e.g.) MOSFLM */
+	cell = uncenter_cell(cell_in, NULL);
 
 	if ( cell_get_reciprocal(template, &asx, &asy, &asz,
 	                         &bsx, &bsy, &bsz,
@@ -638,7 +671,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	free(cand[1]);
 	free(cand[2]);
 
-	return new_cell;
+	/* Reverse the de-centering transformation */
+	new_cell_trans = cell_transform_inverse(new_cell, to_given_cell);
+	cell_free(new_cell);
+
+	return new_cell_trans;
 }
 
 
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index b2cb7b67..8acf2f85 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -55,10 +55,12 @@ extern int cell_is_sensible(UnitCell *cell);
 
 extern void validate_cell(UnitCell *cell);
 
-extern UnitCell *uncenter_cell(UnitCell *in);
+extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr);
 
 extern int bravais_lattice(UnitCell *cell);
 
 extern int right_handed(UnitCell *cell);
 
+extern const char *str_lattice(LatticeType l);
+
 #endif	/* CELL_UTILS_H */
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 37e201e3..e2617dfc 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -599,3 +599,57 @@ const char *cell_rep(UnitCell *cell)
 
 	return "unknown";
 }
+
+
+static UnitCellTransformation *inverse_transformation(UnitCellTransformation *t)
+{
+	/* FIXME: Implementation */
+	return NULL;
+}
+
+
+/**
+ * cell_transform:
+ * @cell: A %UnitCell.
+ * @t: A %UnitCellTransformation.
+ *
+ * Applies @t to @cell.
+ *
+ * Returns: Transformed copy of @cell.
+ *
+ */
+UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
+{
+	UnitCell *out;
+
+	if ( t == NULL ) return NULL;
+
+	out = cell_new();
+	if ( out == NULL ) return NULL;
+
+	/* FIXME: Implementation */
+
+	return out;
+}
+
+
+/**
+ * cell_transform:
+ * @cell: A %UnitCell.
+ * @t: A %UnitCellTransformation.
+ *
+ * Applies the inverse of @t to @cell.
+ *
+ * Returns: Transformed copy of @cell.
+ *
+ */
+UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t)
+{
+	return cell_transform(cell, inverse_transformation(t));
+}
+
+
+void cell_transformation_print(UnitCellTransformation *t)
+{
+
+}
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 4e90a8ad..625884e1 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -68,6 +68,15 @@ 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;
+
 extern UnitCell *cell_new(void);
 extern UnitCell *cell_new_from_cell(UnitCell *orig);
 extern void cell_free(UnitCell *cell);
@@ -127,4 +136,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 void cell_transformation_print(UnitCellTransformation *t);
+
 #endif	/* CELL_H */
diff --git a/tests/centering_check.c b/tests/centering_check.c
index 0dc79a5f..96326d21 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -35,46 +35,63 @@
 #include <cell-utils.h>
 
 
-static int check_cell(UnitCell *cell)
+static int check_cell(UnitCell *cell, const char *text)
 {
 	int err = 0;
 
 	if ( !cell_is_sensible(cell) ) {
-		ERROR("Warning: Unit cell parameters are not sensible.\n");
+		ERROR(" %s unit cell parameters are not sensible.\n", text);
 		err = 1;
 	}
 
 	if ( !bravais_lattice(cell) ) {
-		ERROR("Warning: Unit cell is not a conventional Bravais"
-		      " lattice.\n");
+		ERROR(" %s unit cell is not a conventional Bravais"
+		      " lattice.\n", text);
 		err = 1;
 	}
 
 	if ( !right_handed(cell) ) {
-		ERROR("Warning: Unit cell is not right handed.\n");
+		ERROR(" %s unit cell is not right handed.\n", text);
 		err = 1;
 	}
 
+	if ( err ) cell_print(cell);
+
 	return err;
 }
 
 
-static int check_centering()
+static int check_centering(double a, double b, double c,
+                           double al, double be, double ga,
+                           LatticeType latt, char cen, char ua)
 {
 	UnitCell *cell;
 	UnitCell *n;
+	UnitCellTransformation *t;
 	int fail = 0;
 
-	cell = cell_new();
+	STATUS("Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n",
+	       str_lattice(latt), cen, ua, a, b, c, al, be, ga);
+
+	cell = cell_new_from_parameters(a, b, c,
+	                                deg2rad(al), deg2rad(be), deg2rad(ga));
+	cell_set_lattice_type(cell, latt);
+	cell_set_centering(cell, cen);
+	cell_set_unique_axis(cell, ua);
+
+	if ( check_cell(cell, "Input") ) fail = 1;
+	//cell_print(cell);
+	n = uncenter_cell(cell, &t);
+	if ( n != NULL ) {
+		if ( check_cell(n, "Output") ) fail = 1;
+	} else {
+		fail = 1;
+	}
+
+	STATUS("Transformation was:\n");
+	cell_transformation_print(t);
 
-	cell_set_parameters(cell, 10e-10, 20e-10, 30e-10,
-	                    deg2rad(90.0), deg2rad(90.0), deg2rad(100.0));
-	cell_set_centering(cell, 'C');
-	cell_set_lattice_type(cell, L_MONOCLINIC);
-	cell_set_unique_axis(cell, 'c');
-	if ( check_cell(cell) ) fail = 1;
-	n = uncenter_cell(cell);
-	if ( check_cell(n) ) fail = 1;
+	if ( fail ) ERROR("\n");
 
 	return fail;
 }
@@ -86,28 +103,100 @@ int main(int argc, char *argv[])
 	int fail = 0;
 
 	/* Triclinic P */
-	check_centering(10e-10, 20e-10, 30e-10, 100.0, 120.0, 140.0,
-	                L_TRICLINIC, 'P', '*',
-	                );
+	fail += check_centering(50e-10, 55e-10, 70e-10, 67.0, 70.0, 77.0,
+	                        L_TRICLINIC, 'P', '*');
+
 	/* Monoclinic P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
+	                        L_MONOCLINIC, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
+	                        L_MONOCLINIC, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
+	                        L_MONOCLINIC, 'P', 'c');
+
 	/* Monoclinic A */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
+	                        L_MONOCLINIC, 'A', 'a');
+
 	/* Monoclinic B */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
+	                        L_MONOCLINIC, 'B', 'b');
+
 	/* Monoclinic C */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
+	                        L_MONOCLINIC, 'C', 'c');
+
 	/* Orthorhombic P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'P', '*');
+
 	/* Orthorhombic A */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'A', '*');
+
 	/* Orthorhombic B */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'B', '*');
+
 	/* Orthorhombic C */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'C', '*');
+
 	/* Orthorhombic I */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'I', '*');
+
 	/* Orthorhombic F */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'F', '*');
+
 	/* Tetragonal P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'c');
+
 	/* Tetragonal I */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'c');
+
 	/* Rhombohedral R */
+	fail += check_centering(10e-10, 10e-10, 10e-10, 60.0, 60.0, 60.0,
+	                        L_RHOMBOHEDRAL, 'R', '*');
+
 	/* Hexagonal P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 120.0, 90.0, 90.0,
+	                        L_HEXAGONAL, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 120.0, 90.0,
+	                        L_HEXAGONAL, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 120.0,
+	                        L_HEXAGONAL, 'P', 'c');
+
 	/* Hexagonal H (PDB-speak for rhombohedral) */
+	fail += check_centering(40e-10, 20e-10, 20e-10, 120.0, 90.0, 90.0,
+	                        L_HEXAGONAL, 'H', 'a');
+	fail += check_centering(20e-10, 40e-10, 20e-10, 90.0, 120.0, 90.0,
+	                        L_HEXAGONAL, 'H', 'b');
+	fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
+	                        L_HEXAGONAL, 'H', 'c');
+
 	/* Cubic P */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'P', '*');
+
 	/* Cubic I */
-	/* Cubic F */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'I', '*');
 
+	/* Cubic F */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'F', '*');
 
 	return fail;
 }