From e3f4046056cf92ce5e40e5e715cbcd53df7b0621 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Aug 2019 14:07:31 +0200 Subject: Framework for new unit cell comparison function --- libcrystfel/src/cell-utils.c | 105 +++++++++++++++---------------------------- libcrystfel/src/cell-utils.h | 6 +-- libcrystfel/src/index.c | 16 +------ tests/cellcompare_check.c | 79 +++++++++++++++++++------------- 4 files changed, 87 insertions(+), 119 deletions(-) diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 141eaaab..20ede189 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2341,10 +2341,9 @@ IntegerMatrix *reduce_g6(double *g, double eps) * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * \param pmb: Place to store pointer to matrix * - * Compare \p cell_in with \p reference_in. If they are the same, - * within \p tolerance, taking into account possible permutations of the axes, - * this function returns non-zero and stores the transformation which needs to - * be applied to \p cell_in at \p pmb. + * Compare the \p cell_in with \p reference_in. If they represent the same + * lattice, this function returns non-zero and stores the transformation which + * needs to be applied to \p cell_in at \p pmb. * * Subject to the tolerances, this function will find the transformation which * gives the best match to the reference cell, using the Euclidian norm in @@ -2356,76 +2355,44 @@ IntegerMatrix *reduce_g6(double *g, double eps) * \returns non-zero if the cells match, zero for no match or error. * */ -int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, IntegerMatrix **pmb) +int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, + double *tolerance, RationalMatrix **pmb) { - IntegerMatrix *m; - signed int i[9]; - double a, b, c, al, be, ga; - double min_dist = +INFINITY; - - m = intmat_new(3, 3); - cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga); - *pmb = NULL; - - for ( i[0]=-1; i[0]<=+1; i[0]++ ) { - for ( i[1]=-1; i[1]<=+1; i[1]++ ) { - for ( i[2]=-1; i[2]<=+1; i[2]++ ) { - for ( i[3]=-1; i[3]<=+1; i[3]++ ) { - for ( i[4]=-1; i[4]<=+1; i[4]++ ) { - for ( i[5]=-1; i[5]<=+1; i[5]++ ) { - for ( i[6]=-1; i[6]<=+1; i[6]++ ) { - for ( i[7]=-1; i[7]<=+1; i[7]++ ) { - for ( i[8]=-1; i[8]<=+1; i[8]++ ) { - - UnitCell *test; - int j, k; - int l = 0; - signed int det; - - for ( j=0; j<3; j++ ) - for ( k=0; k<3; k++ ) - intmat_set(m, j, k, i[l++]); - - det = intmat_det(m); - if ( (det != +1) && (det != -1) ) continue; - - test = cell_transform_intmat(cell, m); - - if ( compare_cell_parameters(reference, test, tolerance) ) { - - double at, bt, ct, alt, bet, gat; - double dist; + UnitCell *cell; + UnitCell *reference; + IntegerMatrix *CBint; + RationalMatrix *CiA; + RationalMatrix *CB; + IntegerMatrix *Mcell; + IntegerMatrix *Mref; + double eps; + double G6cell[6]; + double G6ref[6]; - cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat); - dist = g6_distance(at, bt, ct, alt, bet, gat, - a, b, c, al, be, ga); - if ( dist < min_dist ) { - min_dist = dist; - intmat_free(*pmb); - *pmb = intmat_copy(m); - } + /* Un-center both cells */ + reference = uncenter_cell(reference_in, &CBint, NULL); + if ( reference == NULL ) return 0; + CB = rtnl_mtx_from_intmat(CBint); + intmat_free(CBint); - } + cell = uncenter_cell(cell_in, NULL, &CiA); + if ( cell == NULL ) return 0; - cell_free(test); + /* Calculate G6 components for both */ + g6_components(G6cell, cell); + g6_components(G6ref, reference); - } - } - } - } - } - } - } - } - } + /* Convert both to reduced basis (stably) */ + eps = pow(cell_get_volume(cell), 1.0/3.0) * 1e-5; + Mcell = reduce_g6(G6cell, eps); + eps = pow(cell_get_volume(reference), 1.0/3.0) * 1e-5; + Mref = reduce_g6(G6ref, eps); - intmat_free(m); + /* Compare cells (including nearby ones, possibly permuting + * axes if close to equality) */ - if ( isinf(min_dist) ) { - return 0; - } - - /* Solution found */ - return 1; + intmat_free(Mcell); + intmat_free(Mref); + cell_free(reference); + cell_free(cell); } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index fa577269..e46d860b 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -90,9 +90,6 @@ extern double cell_get_volume(UnitCell *cell); extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance); -extern int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, IntegerMatrix **pmb); - extern int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, const double ltl, @@ -108,6 +105,9 @@ extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference double *tolerance, int csl, RationalMatrix **pmb); +extern int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, + double *tolerance, RationalMatrix **pmb); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index b07dd067..aa0cfb6a 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -547,26 +547,13 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, double *tolerance) { UnitCell *out; - IntegerMatrix *im; RationalMatrix *rm; /* Check at all? */ if ( ! ((flags & INDEXING_CHECK_CELL_COMBINATIONS) || (flags & INDEXING_CHECK_CELL_AXES)) ) return 0; - if ( compare_permuted_cell_parameters(crystal_get_cell(cr), target, - tolerance, &im) ) - { - out = cell_transform_intmat(crystal_get_cell(cr), im); - cell_free(crystal_get_cell(cr)); - crystal_set_cell(cr, out); - intmat_free(im); - return 0; - } - - if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS ) - && compare_reindexed_cell_parameters(crystal_get_cell(cr), target, - tolerance, 0, &rm) ) + if ( compare_lattices(crystal_get_cell(cr), target, tolerance, &rm) ) { out = cell_transform_rational(crystal_get_cell(cr), rm); cell_free(crystal_get_cell(cr)); @@ -576,7 +563,6 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, } return 1; - } diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 48b7a9c5..0e98276c 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -116,29 +116,6 @@ static int check_ccp(UnitCell *cell, UnitCell *cref, double *tols, } -static int check_cpcp(UnitCell *cell, UnitCell *cref, double *tols, - int should_match) -{ - IntegerMatrix *m = NULL; - const char *a; - const char *b; - - a = should_match ? "" : "NOT "; - b = " with compare_permuted_cell_parameters"; - - if ( compare_permuted_cell_parameters(cell, cref, tols, &m) != should_match ) - { - complain(cell, cref, a, b); - STATUS("Matrix was:\n"); - intmat_print(m); - intmat_free(m); - return 1; - } - intmat_free(m); - return 0; -} - - static int check_ccpao(UnitCell *cell, UnitCell *cref, double *tols, int should_match) { @@ -206,6 +183,33 @@ static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, } +static void yaro_test() +{ + UnitCell *cell; + UnitCell *reference; + UnitCell *cmatch; + float tols[] = {5, 5, 5, 1.5}; + + cell = cell_new_from_parameters(64.24e-10, 63.94e-10, 64.61e-10, + deg2rad(89.98), deg2rad(89.82), deg2rad(119.87)); + cell_set_unique_axis(cell, 'c'); + cell_set_lattice_type(cell, L_HEXAGONAL); + + reference = cell_new_from_parameters(64.7e-10, 64.7e-10, 65.2e-10, + deg2rad(90.0), deg2rad(90.0), deg2rad(120.0)); + cell_set_unique_axis(reference, 'c'); + cell_set_lattice_type(reference, L_HEXAGONAL); + + cell_print(cell); + cell_print(reference); + cmatch = match_cell(cell, reference, 0, tols, 1); + cell_print(cmatch); +} + + +extern IntegerMatrix *reduce_g6(double *g, double eps); +extern void g6_components(double *g6, UnitCell *cell); + int main(int argc, char *argv[]) { UnitCell *cell, *cref; @@ -216,13 +220,29 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10, - deg2rad(67.0), - deg2rad(70.0), - deg2rad(77.0)); + yaro_test(); + + cref = cell_new_from_parameters(3e-0, 5.196e-0, 2e-0, + deg2rad(103.9166666), + deg2rad(109.4666666), + deg2rad(134.8833333)); cell_set_centering(cref, 'P'); if ( cref == NULL ) return 1; + STATUS("The test cell:\n"); + cell_print(cref); + double g[6]; + g6_components(g, cref); + STATUS("G6: %e %e %e %e %e %e\n", g[0], g[1], g[2], g[3], g[4], g[5]); + double eps = pow(cell_get_volume(cref), 1.0/3.0) * 1e-5; + //eps = 1e-27; + IntegerMatrix *M = reduce_g6(g, eps); + STATUS("The transformation to reduce:\n"); + intmat_print(M); + STATUS("The reduced cell:\n"); + cell = cell_transform_intmat(cref, M); + cell_print(cell); + /* Just rotate cell */ for ( i=0; i<100; i++ ) { @@ -230,7 +250,6 @@ int main(int argc, char *argv[]) if ( cell == NULL ) return 1; if ( check_ccp(cell, cref, tols, 1) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -248,7 +267,6 @@ int main(int argc, char *argv[]) cell = cell_transform_intmat(cref, tr); if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, 1) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -272,7 +290,6 @@ int main(int argc, char *argv[]) cell_free(cell2); if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -301,7 +318,6 @@ int main(int argc, char *argv[]) * cell_transform_rational doesn't set the unique axis (yet?) */ if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; @@ -332,7 +348,6 @@ int main(int argc, char *argv[]) cell_free(cell2); if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; -- cgit v1.2.3