From 4f11a9f8530178cfb510f11c7bc4b1829bbd0be4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2019 11:40:37 +0100 Subject: Keep track of the "un-centering" matrix, as well as the "centering" This makes it easy to reverse the transformation, if required, which it is when comparing centered cells. --- libcrystfel/src/cell-utils.c | 157 ++++++++++++++++++++++++++++++------------- libcrystfel/src/cell-utils.h | 3 +- libcrystfel/src/rational.c | 25 +++++++ libcrystfel/src/rational.h | 2 + libcrystfel/src/xgandalf.c | 2 +- 5 files changed, 142 insertions(+), 47 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 2bfaa3e1..3a5ef782 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -361,18 +361,45 @@ int bravais_lattice(UnitCell *cell) } +static RationalMatrix *create_rtnl_mtx(signed int a1, signed int a2, + signed int b1, signed int b2, + signed int c1, signed int c2, + signed int d1, signed int d2, + signed int e1, signed int e2, + signed int f1, signed int f2, + signed int g1, signed int g2, + signed int h1, signed int h2, + signed int i1, signed int i2) +{ + RationalMatrix *m = rtnl_mtx_new(3, 3); + if ( m == NULL ) return NULL; + rtnl_mtx_set(m, 0, 0, rtnl(a1, a2)); + rtnl_mtx_set(m, 0, 1, rtnl(b1, b2)); + rtnl_mtx_set(m, 0, 2, rtnl(c1, c2)); + rtnl_mtx_set(m, 1, 0, rtnl(d1, d2)); + rtnl_mtx_set(m, 1, 1, rtnl(e1, e2)); + rtnl_mtx_set(m, 1, 2, rtnl(f1, f2)); + rtnl_mtx_set(m, 2, 0, rtnl(g1, g2)); + rtnl_mtx_set(m, 2, 1, rtnl(h1, h2)); + rtnl_mtx_set(m, 2, 2, rtnl(i1, h2)); + return m; +} + + /* 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) */ + * rarely mP). Store the inverse matrix at pCi */ static IntegerMatrix *centering_transformation(UnitCell *in, char *new_centering, LatticeType *new_latt, - char *new_ua) + char *new_ua, + RationalMatrix **pCi) { LatticeType lt; char ua, cen; - IntegerMatrix *t = NULL; + IntegerMatrix *C = NULL; + RationalMatrix *Ci = NULL; lt = cell_get_lattice_type(in); ua = cell_get_unique_axis(in); @@ -382,13 +409,17 @@ static IntegerMatrix *centering_transformation(UnitCell *in, *new_centering = cen; *new_latt = lt; *new_ua = ua; - t = intmat_identity(3); + C = intmat_identity(3); + Ci = rtnl_mtx_identity(3); } if ( cen == 'I' ) { - t = intmat_create_3x3(0, 1, 1, + C = intmat_create_3x3(0, 1, 1, 1, 0, 1, 1, 1, 0); + Ci = create_rtnl_mtx(-1,2, 1,2, 1,2, + 1,2, -1,2, 1,2, + 1,2, 1,2, -1,2); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; @@ -401,9 +432,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'F' ) { - t = intmat_create_3x3(-1, 1, 1, + C = intmat_create_3x3(-1, 1, 1, 1, -1, 1, 1, 1, -1); + Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + 1,2, 0,1, 1,2, + 1,2, 1,2, 0,1); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; @@ -417,9 +451,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) { /* Obverse setting */ - t = intmat_create_3x3(1, -1, 0, + C = intmat_create_3x3(1, -1, 0, 0, 1, -1, 1, 1, 1); + Ci = create_rtnl_mtx( 2,3, 1,3, 1,3, + -1,3, 1,3, 1,3, + -1,3, -2,3, 1,3); assert(lt == L_HEXAGONAL); assert(ua == 'c'); *new_latt = L_RHOMBOHEDRAL; @@ -428,9 +465,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'A' ) { - t = intmat_create_3x3(1, 0, 0, - 0, 1, -1, - 0, 1, 1); + C = intmat_create_3x3(0, 0, -1, + 1, 1, 0, + 1, -1, 0); + Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + 0,1, 1,2, -1,2, + -1,1, 0,1, 0,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -443,9 +483,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'B' ) { - t = intmat_create_3x3(1, 0, -1, - 0, 1, 0, - 1, 0, 1); + C = intmat_create_3x3(1, 1, 0, + 0, 0, 1, + 1, -1, 0); + Ci = create_rtnl_mtx( 1,2, 0,1, 1,2, + 1,2, 0,1, -1,2, + 0,1, 1,1, 0,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -458,9 +501,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'C' ) { - t = intmat_create_3x3(1, -1, 0, + C = intmat_create_3x3(1, -1, 0, 1, 1, 0, 0, 0, 1); + Ci = create_rtnl_mtx( 1,2, 1,2, 0,1, + -1,2, 1,2, 0,1, + 0,1, 0,1, 1,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -472,45 +518,57 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } } - return t; + *pCi = Ci; + return C; } /** * uncenter_cell: * @in: A %UnitCell - * @t: Location at which to store the transformation which was used. + * @C: Location at which to store the centering transformation + * @Ci: Location at which to store the inverse centering transformation + * + * Turns any cell into a primitive one, e.g. for comparison purposes. * - * 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. + * The transformation which was used is stored at @Ci. The centering + * transformation, which is the transformation you should apply if you want to + * get back the original cell, will be stored at @C. Either or both of these + * can be NULL if you don't need that information. * - * Returns: a primitive version of @in in a conventional (unique axis c) - * setting. + * Returns: a primitive version of @in. * */ -UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **t) +UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi) { - IntegerMatrix *tt; + IntegerMatrix *C; + RationalMatrix *Ci; char new_centering; LatticeType new_latt; char new_ua; UnitCell *out; - tt = centering_transformation(in, &new_centering, &new_latt, &new_ua); - if ( tt == NULL ) return NULL; + C = centering_transformation(in, &new_centering, &new_latt, + &new_ua, &Ci); + if ( C == NULL ) return NULL; - out = cell_transform_intmat_inverse(in, tt); + out = cell_transform_rational(in, Ci); 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; + if ( pC != NULL ) { + *pC = C; } else { - intmat_free(tt); + intmat_free(C); + } + + if ( pCi != NULL ) { + *pCi = Ci; + } else { + rtnl_mtx_free(Ci); } return out; @@ -578,12 +636,12 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, ¢ering); + template = uncenter_cell(template_in, ¢ering, NULL); if ( template == NULL ) return NULL; /* The candidate cell is also uncentered, because it might be centered * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL); + cell = uncenter_cell(cell_in, NULL, NULL); if ( cell == NULL ) return NULL; if ( cell_get_reciprocal(template, &asx, &asy, &asz, @@ -820,11 +878,11 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &to_given_cell); + template = uncenter_cell(template_in, &to_given_cell, NULL); /* The candidate cell is also uncentered, because it might be centered * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL); + cell = uncenter_cell(cell_in, NULL, NULL); /* Get the lengths to match */ if ( cell_get_cartesian(template, &ax, &ay, &az, @@ -1904,8 +1962,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, { UnitCell *cell; UnitCell *reference; - IntegerMatrix *centering_reference; - IntegerMatrix *centering_cell; + IntegerMatrix *CBint; + RationalMatrix *CiA; + RationalMatrix *CB; RationalMatrix *m; double a, b, c, al, be, ga; double av[3], bv[3], cv[3]; @@ -1915,13 +1974,18 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, int ncand_a, ncand_b, ncand_c; int i; int ia, ib; + RationalMatrix *MCiA; + RationalMatrix *CBMCiA; + int ret = 0; /* Actually compare against primitive version of reference */ - reference = uncenter_cell(reference_in, ¢ering_reference); + reference = uncenter_cell(reference_in, &CBint, NULL); if ( reference == NULL ) return 0; + CB = rtnl_mtx_from_intmat(CBint); + intmat_free(CBint); /* Actually compare primitive version of cell */ - cell = uncenter_cell(cell_in, ¢ering_cell); + cell = uncenter_cell(cell_in, NULL, &CiA); if ( cell == NULL ) return 0; /* Get target parameters */ @@ -1955,6 +2019,8 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, } m = rtnl_mtx_new(3, 3); + MCiA = rtnl_mtx_new(3, 3); + CBMCiA = rtnl_mtx_new(3, 3); for ( ia=0; iacols == A->cols); + assert(ans->rows == B->rows); + assert(A->cols == B->rows); + + for ( i=0; irows; i++ ) { + for ( j=0; jcols; j++ ) { + int k; + Rational sum = rtnl_zero(); + for ( k=0; krows; k++ ) { + Rational add; + add = rtnl_mul(rtnl_mtx_get(A, i, k), + rtnl_mtx_get(B, k, j)); + sum = rtnl_add(sum, add); + } + rtnl_mtx_set(ans, i, j, sum); + } + } +} + + void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans) { int i, j; diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 8590e920..9ed20bfd 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -89,6 +89,8 @@ extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); extern RationalMatrix *rtnl_mtx_identity(int rows); extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m); extern void rtnl_mtx_free(RationalMatrix *mtx); +extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, + RationalMatrix *ans); extern void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans); extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c index b122cb4c..51819956 100644 --- a/libcrystfel/src/xgandalf.c +++ b/libcrystfel/src/xgandalf.c @@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, xgandalf_private_data->cellTemplate = cell; - UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation); + UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation, NULL); UnitCell *uc = cell_new_from_cell(primitiveCell); reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform); -- cgit v1.2.3