Use IntegerMatrix for all unit cell transformations

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

1 files changed, 140 insertions, 274 deletions

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;
 }