From d0e9fd4cb3407916e6ceb455a20e01c09847316a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 24 Jan 2019 11:58:22 +0100 Subject: Automatic centering determination after transformation --- libcrystfel/src/cell.c | 191 ++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 172 insertions(+), 19 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index a26e26a6..131e1ac8 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -97,6 +97,20 @@ struct _unitcell { char unique_axis; }; +typedef enum { + CMASK_P = 1<<0, + CMASK_A = 1<<1, + CMASK_B = 1<<2, + CMASK_C = 1<<3, + CMASK_I = 1<<4, + CMASK_F = 1<<5, + CMASK_H = 1<<6, + CMASK_R = 1<<7 +} CenteringMask; + +#define CMASK_ALL (CMASK_P | CMASK_A | CMASK_B | CMASK_C | CMASK_I \ + | CMASK_F | CMASK_H | CMASK_R) + /************************** Setters and Constructors **************************/ @@ -688,28 +702,167 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) } -static char determine_centering(IntegerMatrix *m, char cen) +static int centering_has_point(char cen, Rational *p) +{ + /* First, put the point into the range 0..1 */ + while ( rtnl_cmp(p[0], rtnl_zero()) < 0 ) p[0] = rtnl_add(p[0], rtnl(1, 1)); + while ( rtnl_cmp(p[1], rtnl_zero()) < 0 ) p[1] = rtnl_add(p[1], rtnl(1, 1)); + while ( rtnl_cmp(p[2], rtnl_zero()) < 0 ) p[2] = rtnl_add(p[2], rtnl(1, 1)); + while ( rtnl_cmp(p[0], rtnl(1, 1)) >= 0 ) p[0] = rtnl_sub(p[0], rtnl(1, 1)); + while ( rtnl_cmp(p[1], rtnl(1, 1)) >= 0 ) p[1] = rtnl_sub(p[1], rtnl(1, 1)); + while ( rtnl_cmp(p[2], rtnl(1, 1)) >= 0 ) p[2] = rtnl_sub(p[2], rtnl(1, 1)); + + /* 0,0,0 is present in all centerings */ + if ( (rtnl_cmp(p[0], rtnl_zero()) == 0) + && (rtnl_cmp(p[1], rtnl_zero()) == 0) + && (rtnl_cmp(p[2], rtnl_zero()) == 0) ) return 1; + + /* Only I has 1/2 , 1/2, 1/2 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && (cen == 'I') ) return 1; + + /* A or F has 0 , 1/2, 1/2 */ + if ( (rtnl_cmp(p[0], rtnl_zero()) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && ((cen == 'A') || (cen == 'F')) ) return 1; + + /* B or F has 1/2 , 0 , 1/2 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl_zero()) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && ((cen == 'B') || (cen == 'F')) ) return 1; + + /* C or F has 1/2 , 1/2 , 0 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl_zero()) == 0) + && ((cen == 'C') || (cen == 'F')) ) return 1; + + /* H has 2/3 , 1/3 , 1/3 */ + if ( (rtnl_cmp(p[0], rtnl(2,3)) == 0) + && (rtnl_cmp(p[1], rtnl(1,3)) == 0) + && (rtnl_cmp(p[2], rtnl(1,3)) == 0) + && (cen == 'H') ) return 1; + + /* H has 1/3 , 2/3 , 2/3 */ + if ( (rtnl_cmp(p[0], rtnl(1,3)) == 0) + && (rtnl_cmp(p[1], rtnl(2,3)) == 0) + && (rtnl_cmp(p[2], rtnl(2,3)) == 0) + && (cen == 'H') ) return 1; + + return 0; +} + + +static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc, + char cen) +{ + /* Skip test if this centering isn't even a candidate */ + if ( !(*cmask & c) ) return; + + if ( !centering_has_point(cen, nc) ) { + *cmask |= c; + *cmask ^= c; + } +} + + +/* Check if the point x,y,z in the original cell matches any lattice point + * in the transformed cell */ +static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, + Rational x, Rational y, Rational z) +{ + Rational c[3] = {x, y, z}; + Rational nc[3]; + + /* Transform the lattice point */ + rtnl_mtx_solve(m, c, nc); + + /* Eliminate any centerings which don't include the transformed point */ + maybe_eliminate(CMASK_P, cmask, nc, 'P'); + maybe_eliminate(CMASK_R, cmask, nc, 'R'); + maybe_eliminate(CMASK_A, cmask, nc, 'A'); + maybe_eliminate(CMASK_B, cmask, nc, 'B'); + maybe_eliminate(CMASK_C, cmask, nc, 'C'); + maybe_eliminate(CMASK_I, cmask, nc, 'I'); + maybe_eliminate(CMASK_F, cmask, nc, 'F'); + maybe_eliminate(CMASK_H, cmask, nc, 'H'); +} + + +/* Check if the point x,y,z in the transformed cell matches any lattice point + * in the original cell. If not, eliminate "exclude" from "*mask". */ +static void check_point_bwd(RationalMatrix *m, CenteringMask *mask, + char cen, CenteringMask exclude, + Rational x, Rational y, Rational z) { - Rational c[3]; Rational nc[3]; + Rational c[3] = {x, y, z}; + + rtnl_mtx_mult(m, c, nc); + + if ( !centering_has_point(cen, nc) ) { + *mask |= exclude; + *mask ^= exclude; /* Unset bits */ + } +} + + +static char cmask_decode(CenteringMask mask) +{ + char res[32]; + + res[0] = '\0'; + + if ( mask & CMASK_H ) strcat(res, "H"); + if ( mask & CMASK_F ) strcat(res, "F"); + if ( mask & CMASK_I ) strcat(res, "I"); + if ( mask & CMASK_A ) strcat(res, "A"); + if ( mask & CMASK_B ) strcat(res, "B"); + if ( mask & CMASK_C ) strcat(res, "C"); + if ( mask & CMASK_P ) strcat(res, "P"); + if ( mask & CMASK_R ) strcat(res, "R"); + + STATUS("possible centerings: %s\n", res); + return res[0]; +} + + +static char determine_centering(RationalMatrix *m, char cen) +{ + CenteringMask cmask = CMASK_ALL; + + /* Check whether the current centering can provide all the lattice + * points for the transformed cell. Eliminate any centerings for which + * it can't. */ + check_point_bwd(m, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_bwd(m, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + + /* Check whether the current centering's lattice points will all + * coincide with lattice points in the new centering. Eliminate any + * centerings for which they don't (they give "excess lattice points"). */ + switch ( cen ) { + + case 'P' : + case 'R' : + break; + + case 'C' : + check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + break; + + /* FIXME: Everything else */ + + } - c[0] = rtnl(1, 2); - c[1] = rtnl(1, 2); - c[2] = rtnl_zero(); - intmat_rationalvec_mult(m, c, nc); - STATUS("%s,%s,%s -> %s,%s,%s\n", - rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), - rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); - - c[0] = rtnl_zero(); - c[1] = rtnl_zero(); - c[2] = rtnl_zero(); - intmat_solve_rational(m, nc, c); - STATUS("%s,%s,%s <- %s,%s,%s\n", - rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), - rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); - - return 'A'; + return cmask_decode(cmask); } -- cgit v1.2.3