aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/cell-utils.c105
-rw-r--r--libcrystfel/src/cell-utils.h6
-rw-r--r--libcrystfel/src/index.c16
-rw-r--r--tests/cellcompare_check.c79
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;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 17:52:57 +0000