More work on transformations

8 files changed, 186 insertions, 12 deletions

diff --git a/Makefile.am b/Makefile.am
index 953a04df..f0adc640 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -8,12 +8,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
 
 noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/pr_gradient_check tests/symmetry_check \
-                  tests/centering_check
+                  tests/centering_check tests/transformation_check
 
 TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
         tests/third_merge_check tests/fourth_merge_check \
         tests/integration_check tests/pr_gradient_check \
-        tests/symmetry_check tests/centering_check
+        tests/symmetry_check tests/centering_check tests/transformation_check
 
 EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
               tests/third_merge_check tests/fourth_merge_check
@@ -87,6 +87,10 @@ tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \
 tests_centering_check_SOURCES = tests/centering_check.c libcrystfel/src/cell.c \
                                 libcrystfel/src/cell-utils.c
 
+tests_transformation_check_SOURCES = tests/transformation_check.c \
+                                     libcrystfel/src/cell.c \
+                                     libcrystfel/src/cell-utils.c
+
 INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data
 
 EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \
diff --git a/doc/reference/CrystFEL-docs.sgml b/doc/reference/CrystFEL-docs.sgml
index 89b4832e..e1da1784 100644
--- a/doc/reference/CrystFEL-docs.sgml
+++ b/doc/reference/CrystFEL-docs.sgml
@@ -34,6 +34,7 @@
   <chapter>
     <title>Unit cells</title>
         <xi:include href="xml/unitcell.xml"/>
+        <xi:include href="xml/cell-utils.xml"/>
   </chapter>
 
   <chapter>
diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt
index 284740a1..3ccc1511 100644
--- a/doc/reference/CrystFEL-sections.txt
+++ b/doc/reference/CrystFEL-sections.txt
@@ -88,12 +88,28 @@ cell_set_pointgroup
 cell_set_reciprocal
 cell_set_spacegroup
 <SUBSECTION>
+tfn_combine
+tfn_identity
+tfn_inverse
+tfn_print
+tfn_vector
+tfn_free
+</SECTION>
+
+<SECTION>
+<FILE>cell-utils</FILE>
 cell_rotate
 rotate_cell
 cell_print
 resolution
 match_cell
 match_cell_ab
+cell_is_sensible
+validate_cell
+uncenter_cell
+bravais_lattice
+right_handed
+str_lattice
 </SECTION>
 
 <SECTION>
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 6c711642..9b65cae3 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -47,6 +47,19 @@
 #include "image.h"
 
 
+/**
+ * SECTION:cell-utils
+ * @short_description: Unit cell utilities
+ * @title: Unit cell utilities
+ * @section_id:
+ * @see_also:
+ * @include: "cell-utils.h"
+ * @Image:
+ *
+ * There are some utility functions associated with the core %UnitCell.
+ **/
+
+
 /* Weighting factor of lengths relative to angles */
 #define LWEIGHT (10.0e-9)
 
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index beaa341a..7a62f51c 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -607,33 +607,55 @@ struct _unitcelltransformation
 };
 
 
+/**
+ * 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;
+
 	out = tfn_identity();
-	if ( out == NULL ) return NULL;
+	if ( out == NULL ) {
+		gsl_matrix_free(m);
+		return NULL;
+	}
 
-	perm = gsl_permutation_alloc(t->m->size1);
+	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(t->m->size1, t->m->size2);
+	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(t->m, perm, &s) ) {
+	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(t->m, perm, inv)  ) {
+	if ( gsl_linalg_LU_invert(m, perm, inv)  ) {
 		ERROR("Couldn't invert matrix\n");
 		gsl_permutation_free(perm);
 		return NULL;
@@ -641,6 +663,7 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
 	gsl_permutation_free(perm);
 
 	gsl_matrix_free(out->m);
+	gsl_matrix_free(m);
 	out->m = inv;
 	return out;
 }
@@ -651,7 +674,8 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
  * @cell: A %UnitCell.
  * @t: A %UnitCellTransformation.
  *
- * Applies @t to @cell.
+ * Applies @t to @cell.  Note that the lattice type, centering and unique axis
+ * information will not be preserved.
  *
  * Returns: Transformed copy of @cell.
  *
@@ -675,7 +699,7 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
 	                         &cx, &cy, &cz);
 
 	m = gsl_matrix_alloc(3,3);
-	a = gsl_matrix_alloc(3,3);
+	a = gsl_matrix_calloc(3,3);
 	if ( (m == NULL) || (a == NULL) ) {
 		cell_free(out);
 		return NULL;
@@ -722,7 +746,13 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
  */
 UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t)
 {
-	return cell_transform(cell, tfn_inverse(t));
+	UnitCellTransformation *inv;
+	UnitCell *out;
+
+	inv = tfn_inverse(t);
+	out = cell_transform(cell, inv);
+	tfn_free(inv);
+	return out;
 }
 
 
@@ -769,7 +799,7 @@ void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
 	gsl_matrix *n;
 
 	n = gsl_matrix_alloc(3, 3);
-	a = gsl_matrix_alloc(3, 3);
+	a = gsl_matrix_calloc(3, 3);
 	if ( (n == NULL) || (a == NULL) ) {
 		return;
 	}
@@ -795,6 +825,18 @@ void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
 }
 
 
+/**
+ * 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));
@@ -804,6 +846,13 @@ double *tfn_vector(double a, double b, double c)
 }
 
 
+/**
+ * tfn_print:
+ * @t: A %UnitCellTransformation
+ *
+ * Prints information about @t to stderr.
+ *
+ */
 void tfn_print(UnitCellTransformation *t)
 {
 	gsl_permutation *perm;
@@ -842,3 +891,17 @@ void tfn_print(UnitCellTransformation *t)
 
 	STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s));
 }
+
+
+/**
+ * tfn_free:
+ * @t: A %UnitCellTransformation
+ *
+ * Frees all resources associated with @t.
+ *
+ */
+void tfn_free(UnitCellTransformation *t)
+{
+	gsl_matrix_free(t->m);
+	free(t);
+}
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index ac1bd555..c7b8f8d6 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -146,5 +146,6 @@ extern void tfn_combine(UnitCellTransformation *t,
 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);
 
 #endif	/* CELL_H */
diff --git a/tests/centering_check.c b/tests/centering_check.c
index bfee2806..eca67c85 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -70,7 +70,8 @@ static int check_centering(double a, double b, double c,
 	UnitCellTransformation *t;
 	int fail = 0;
 
-	STATUS("Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n",
+	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,
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
new file mode 100644
index 00000000..771e44fe
--- /dev/null
+++ b/tests/transformation_check.c
@@ -0,0 +1,75 @@
+/*
+ * transformation_check.c
+ *
+ * Check that unit cell transformations work
+ *
+ * Copyright © 2012 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 <cell.h>
+#include <cell-utils.h>
+
+
+int main(int argc, char *argv[])
+{
+	int fail = 0;
+	UnitCell *cell, *cnew, *cback;
+	UnitCellTransformation *tfn, *inv;
+
+	cell = cell_new_from_parameters(50e-10, 55e-10, 70e-10,
+	                                deg2rad(67.0),
+	                                deg2rad(70.0),
+	                                deg2rad(77.0));
+	if ( cell == NULL ) return 1;
+
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+
+	tfn_combine(tfn, tfn_vector(0,1,0),
+	                 tfn_vector(1,0,0),
+	                 tfn_vector(0,0,1));
+
+
+	cell_print(cell);
+	tfn_print(tfn);
+
+	cnew = cell_transform(cell, tfn);
+	cell_print(cnew);
+
+	cback = cell_transform_inverse(cnew, tfn);
+	inv = tfn_inverse(tfn);
+	tfn_print(inv);
+	cell_print(cback);
+
+	cell_get_cartesian(cell, &ax1, &ay1, &az1,
+	                         &by1, &by1, &bz1,
+	                         &cx1, &cy1, &cz1);
+	cell_get_cartesian(cback, &ax, &ay, &az, &by, &by, &bz, &cx, &cy, &cz);
+
+	return fail;
+}