From 4337cafe052c4ad238c969dfa4cb7c7ac52f5e07 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 25 Oct 2018 16:48:36 +0200 Subject: Use IntegerMatrix for all unit cell transformations Get rid of UnitCellTransformation, a thin wrapper which didn't do anything. --- libcrystfel/src/cell-utils.c | 260 ++++++++---------------- libcrystfel/src/cell-utils.h | 3 +- libcrystfel/src/cell.c | 414 +++++++++++++-------------------------- libcrystfel/src/cell.h | 24 +-- libcrystfel/src/integer_matrix.c | 39 ++++ libcrystfel/src/integer_matrix.h | 9 + libcrystfel/src/taketwo.c | 3 +- libcrystfel/src/xgandalf.c | 12 +- 8 files changed, 284 insertions(+), 480 deletions(-) (limited to 'libcrystfel') 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, ¢ering); 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 #include +#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); } -- cgit v1.2.3