From 5911a4b627e1c6676c5522a8ecd61a2a834f544e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 12 Sep 2012 17:57:38 +0200 Subject: More work on transformations --- Makefile.am | 8 +++- doc/reference/CrystFEL-docs.sgml | 1 + doc/reference/CrystFEL-sections.txt | 16 ++++++++ libcrystfel/src/cell-utils.c | 13 ++++++ libcrystfel/src/cell.c | 81 ++++++++++++++++++++++++++++++++----- libcrystfel/src/cell.h | 1 + tests/centering_check.c | 3 +- tests/transformation_check.c | 75 ++++++++++++++++++++++++++++++++++ 8 files changed, 186 insertions(+), 12 deletions(-) create mode 100644 tests/transformation_check.c 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 @@ Unit cells + 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 +tfn_combine +tfn_identity +tfn_inverse +tfn_print +tfn_vector +tfn_free + + +
+cell-utils 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
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 + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#include +#include + + +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; +} -- cgit v1.2.3