aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2019-08-21 16:27:34 +0200
committerThomas White <taw@physics.org>2019-08-22 17:03:28 +0200
commitbeedf1a19a31ecd83f887ceba6e7f08eba6e930b (patch)
tree4b6604fa35592cb8c635bfddeeee827091d9dc1b /libcrystfel/src
parent0b1db945a12c0e68491412baf322b761511eb016 (diff)
Final implementation of Niggli-based cell comparison
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/cell-utils.c309
1 files changed, 238 insertions, 71 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 961a6554..72225a8c 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -2169,7 +2169,7 @@ static int is_burger(struct g6 g, double eps)
}
-static int is_niggli(struct g6 g, double eps)
+static int UNUSED is_niggli(struct g6 g, double eps)
{
if ( !is_burger(g, eps) ) return 0;
@@ -2188,20 +2188,27 @@ static int is_niggli(struct g6 g, double eps)
}
-static void debug_lattice(struct g6 g, double eps)
+static signed int eps_sign(double v, double eps)
+{
+ return GT(v, 0.0) ? +1 : -1;
+}
+
+
+static void debug_lattice(struct g6 g, double eps, int step)
{
- STATUS("%e %e %e %e %e %e",
+#if 0
+ STATUS("After step %i: %e %e %e %e %e %e --", step,
g.A/1e-0, g.B/1e-0, g.C/1e-0,
g.D/1e-0, g.E/1e-0, g.F/1e-0);
if ( is_burger(g, eps) ) STATUS(" B");
if ( is_niggli(g, eps) ) STATUS(" N");
+ if ( eps_sign(g.D, eps)*eps_sign(g.E, eps)*eps_sign(g.F, eps) > 0 ) {
+ STATUS(" I");
+ } else {
+ STATUS(" II");
+ }
STATUS("\n");
-}
-
-
-static signed int eps_sign(double v, double eps)
-{
- return GT(v, 0.0) ? +1 : -1;
+#endif
}
@@ -2218,6 +2225,23 @@ static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M)
}
+/* Replace the elements of T with those of T*M */
+static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M)
+{
+ int i, j;
+ RationalMatrix *tmp;
+
+ tmp = rtnl_mtx_new(3, 3);
+ rtnl_mtx_mtxmult(T, M, tmp);
+ for ( i=0; i<3; i++ ) {
+ for ( j=0; j<3; j++ ) {
+ rtnl_mtx_set(T, i, j, rtnl_mtx_get(tmp, i, j));
+ }
+ }
+ rtnl_mtx_free(tmp);
+}
+
+
/* Cell volume from G6 components */
static double g6_volume(struct g6 g)
{
@@ -2238,11 +2262,12 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
eps = pow(g6_volume(g), 1.0/3.0) * epsrel;
eps = eps*eps;
- STATUS("eps = %e\n", eps);
T = intmat_identity(3);
M = intmat_new(3, 3);
+ debug_lattice(g, eps, 0);
+
do {
int done_s1s2;
@@ -2250,8 +2275,6 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
do {
done_s1s2 = 1;
- debug_lattice(g, eps);
-
if ( GT(g.A, g.B)
|| (EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E))) )
{
@@ -2266,7 +2289,8 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
intmat_set(M, 2, 2, -1);
mult_in_place(T, M);
- STATUS("step 1\n");
+ debug_lattice(g, eps, 1);
+
} else if ( GT(g.B, g.C)
|| (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) )
@@ -2282,7 +2306,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
intmat_set(M, 2, 1, -1);
mult_in_place(T, M);
- STATUS("step 2\n");
+ debug_lattice(g, eps, 2);
done_s1s2 = 0;
}
@@ -2290,72 +2314,62 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
finished = 0;
- debug_lattice(g, eps);
-
/* K-G paper says g3*g4*g3, which I assume is a misprint */
if ( eps_sign(g.D, eps) * eps_sign(g.E, eps) * eps_sign(g.F, eps) > 0 ) {
intmat_zero(M);
- intmat_set(M, 0, 0, GT(g.D, 0.0) ? 1 : -1);
- intmat_set(M, 1, 1, GT(g.E, 0.0) ? 1 : -1);
- intmat_set(M, 2, 2, GT(g.F, 0.0) ? 1 : -1);
+ intmat_set(M, 0, 0, LT(g.D, 0.0) ? -1 : 1);
+ intmat_set(M, 1, 1, LT(g.E, 0.0) ? -1 : 1);
+ intmat_set(M, 2, 2, LT(g.F, 0.0) ? -1 : 1);
mult_in_place(T, M);
g.D = fabs(g.D);
g.E = fabs(g.E);
g.F = fabs(g.F);
- STATUS("step 3\n");
+ debug_lattice(g, eps, 3);
} else {
- int r, c;
- int have_rc = 0;
-
- intmat_zero(M);
- intmat_set(M, 0, 0, 1);
- intmat_set(M, 1, 1, 1);
- intmat_set(M, 2, 2, 1);
+ int z = 4;
+ signed int ijk[3] = { 1, 1, 1 };
if ( GT(g.D, 0.0) ) {
- intmat_set(M, 0, 0, -1);
+ ijk[0] = -1;
} else if ( !LT(g.D, 0.0) ) {
- r = 0; c = 0;
- have_rc = 1;
+ z = 0;
}
if ( GT(g.E, 0.0) ) {
- intmat_set(M, 1, 1, -1);
+ ijk[1] = -1;
} else if ( !LT(g.D, 0.0) ) {
- r = 1; c = 1;
- have_rc = 1;
- }
+ z = 1;
+ }
if ( GT(g.F, 0.0) ) {
- intmat_set(M, 2, 2, -1);
+ ijk[2] = -1;
} else if ( !LT(g.D, 0.0) ) {
- r = 2; c = 2;
- have_rc = 1;
+ z = 2;
}
- if ( LT(g.D*g.E*g.F, 0.0) ) {
- assert(have_rc);
- if ( have_rc ) {
- intmat_set(M, r, c, -1);
- } else {
- ERROR("Don't have r,c!\n");
- //abort();
+ if ( ijk[0]*ijk[1]*ijk[2] < 0 ) {
+ if ( z == 4 ) {
+ ERROR("No element in reduction step 4");
+ abort();
}
+ ijk[z] = -1;
}
+ intmat_zero(M);
+ intmat_set(M, 0, 0, ijk[0]);
+ intmat_set(M, 1, 1, ijk[1]);
+ intmat_set(M, 2, 2, ijk[2]);
mult_in_place(T, M);
g.D = -fabs(g.D);
g.E = -fabs(g.E);
g.F = -fabs(g.F);
- STATUS("step 4\n");
+ debug_lattice(g, eps, 4);
}
- debug_lattice(g, eps);
-
if ( GT(fabs(g.D), g.B)
|| (EQ(g.B, g.D) && LT(2.0*g.E, g.F))
|| (EQ(g.B, -g.D) && LT(g.F, 0.0)) )
@@ -2373,7 +2387,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
g.D = -2*s*g.B + g.D;
g.E = g.E - s*g.F;
- STATUS("step 5\n");
+ debug_lattice(g, eps, 5);
} else if ( GT(fabs(g.E), g.A)
|| (EQ(g.A, g.E) && LT(2.0*g.D, g.F))
@@ -2392,7 +2406,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
g.D = g.D - s*g.F;
g.E = -2*s*g.A + g.E;
- STATUS("step 6\n");
+ debug_lattice(g, eps, 6);
} else if ( GT(fabs(g.F), g.A)
|| (EQ(g.A, g.F) && LT(2.0*g.D, g.E))
@@ -2411,10 +2425,10 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
g.D = g.D - s*g.E;
g.F = -2*s*g.A + g.F;
- STATUS("step 7\n");
+ debug_lattice(g, eps, 7);
- } else if ( LT(g.A+g.B+g.C+g.D+g.E+g.F, 0.0)
- || ( (EQ(g.A+g.B+g.C+g.D+g.E+g.F, 0.0)
+ } else if ( LT(g.A+g.B+g.D+g.E+g.F, 0.0) /* not g.C */
+ || ( (EQ(g.A+g.B+g.D+g.E+g.F, 0.0)
&& GT(2.0*g.A + 2.0*g.E + g.F, 0.0)) ) )
{
intmat_zero(M);
@@ -2428,7 +2442,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
g.D = 2.0*g.B + g.D + g.F;
g.E = 2.0*g.A + g.E + g.F;
- STATUS("step 8\n");
+ debug_lattice(g, eps, 8);
} else {
finished = 1;
@@ -2436,7 +2450,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
} while ( !finished );
- debug_lattice(g, eps);
+ debug_lattice(g, eps, 99);
assert(is_burger(g, eps));
assert(is_niggli(g, eps));
@@ -2446,6 +2460,96 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
}
+static double cell_diff(UnitCell *cell, double a, double b, double c,
+ double al, double be, double ga)
+{
+ double ta, tb, tc, tal, tbe, tga;
+ double diff = 0.0;
+ cell_get_parameters(cell, &ta, &tb, &tc, &tal, &tbe, &tga);
+ diff += fabs(a - ta);
+ diff += fabs(b - tb);
+ diff += fabs(c - tc);
+ return diff;
+}
+
+
+static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference,
+ IntegerMatrix *RB, IntegerMatrix *CB,
+ const double *tols)
+{
+ IntegerMatrix *m;
+ int i[9];
+ double a, b, c, al, be, ga;
+ double min_dist = +INFINITY;
+ IntegerMatrix *best_m = NULL;
+ RationalMatrix *RBCB;
+ RationalMatrix *rRB;
+ RationalMatrix *rCB;
+
+ m = intmat_new(3, 3);
+ cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+
+ RBCB = rtnl_mtx_new(3, 3);
+ rRB = rtnl_mtx_from_intmat(RB);
+ rCB = rtnl_mtx_from_intmat(CB);
+ rtnl_mtx_mtxmult(rRB, rCB, RBCB);
+ rtnl_mtx_free(rRB);
+ rtnl_mtx_free(rCB);
+
+ 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 *nc;
+ int j, k;
+ int l = 0;
+ signed int det;
+ UnitCell *tmp;
+
+ 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 ) continue;
+
+ tmp = cell_transform_intmat(cell, m);
+ nc = cell_transform_rational(tmp, RBCB);
+ cell_free(tmp);
+
+ if ( compare_cell_parameters(nc, reference, tols) ) {
+ double dist = cell_diff(nc, a, b, c, al, be, ga);
+ if ( dist < min_dist ) {
+ intmat_free(best_m);
+ min_dist = dist;
+ best_m = intmat_copy(m);
+ }
+ }
+
+ cell_free(nc);
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ intmat_free(m);
+ rtnl_mtx_free(RBCB);
+
+ return best_m;
+}
+
+
/**
* \param cell_in: A UnitCell
* \param reference_in: Another UnitCell
@@ -2461,7 +2565,9 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel)
*
* Only the cell parameters will be compared. The relative orientations are
* irrelevant. The tolerances will be applied to the transformed copy of
- * \p cell_in.
+ * \p cell_in, i.e. the version of the input cell which looks similar to
+ * \p reference_in. Subject to the tolerances, the cell will be chosen which
+ * has the lowest total absolute error in unit cell axis lengths.
*
* This is the right function to use for deciding if an indexing solution
* matches a reference cell or not.
@@ -2476,37 +2582,98 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in,
{
UnitCell *cell;
UnitCell *reference;
- IntegerMatrix *CBint;
+ IntegerMatrix *CB;
RationalMatrix *CiA;
- RationalMatrix *CB;
- IntegerMatrix *Mcell;
- IntegerMatrix *Mref;
+ IntegerMatrix *RA;
+ IntegerMatrix *RB;
+ IntegerMatrix *P;
+ UnitCell *cell_reduced;
+ UnitCell *reference_reduced;
+ UnitCell *match;
struct g6 g6cell;
struct g6 g6ref;
+ //STATUS("The input cell:\n");
+ //cell_print(cell_in);
+
+ //STATUS("The reference cell:\n");
+ //cell_print(reference_in);
+
/* 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);
+ reference = uncenter_cell(reference_in, &CB, NULL);
+ if ( reference == NULL ) return NULL;
cell = uncenter_cell(cell_in, NULL, &CiA);
- if ( cell == NULL ) return 0;
+ if ( cell == NULL ) return NULL;
/* Calculate G6 components for both */
g6cell = cell_get_G6(cell);
g6ref = cell_get_G6(reference);
/* Convert both to reduced basis (stably) */
- Mcell = reduce_g6(g6cell, 1e-5);
- Mref = reduce_g6(g6ref, 1e-5);
+ RA = reduce_g6(g6cell, 1e-5);
+ //STATUS("------------------------------------------\n");
+ RB = reduce_g6(g6ref, 1e-5);
+
+ cell_reduced = cell_transform_intmat(cell, RA);
- /* Compare cells (including nearby ones, possibly permuting
- * axes if close to equality) */
- STATUS("Reduced cell:\n");
+ //STATUS("Reduced cell:\n");
+ //cell_print(cell_reduced);
+ //STATUS("Reduced reference:\n");
+ reference_reduced = cell_transform_intmat(reference, RB);
+ //cell_print(reference_reduced);
- intmat_free(Mcell);
- intmat_free(Mref);
+ /* The primitive (non-reduced) cells are no longer needed */
cell_free(reference);
cell_free(cell);
+
+ /* Within tolerance? */
+ P = check_permutations(cell_reduced, reference_in, RB, CB, tols);
+ if ( P != NULL ) {
+
+ match = cell_transform_intmat(cell_reduced, P);
+
+ //STATUS("Original cell transformed to look like reference:\n");
+ //cell_print(match);
+
+ /* Match found. Combined matrix requested? */
+ if ( pmb != NULL ) {
+
+ /* Calculate combined matrix: CB.RiB.RA.CiA */
+ RationalMatrix *rRA = rtnl_mtx_from_intmat(RA);
+ RationalMatrix *rP = rtnl_mtx_from_intmat(P);
+ IntegerMatrix *RiB = intmat_inverse(RB);
+ RationalMatrix *rRiB = rtnl_mtx_from_intmat(RiB);
+ RationalMatrix *rCB = rtnl_mtx_from_intmat(CB);
+ RationalMatrix *comb = rtnl_mtx_identity(3);
+
+ rtnl_mult_in_place(comb, CiA);
+ rtnl_mult_in_place(comb, rRA);
+ rtnl_mult_in_place(comb, rP);
+ rtnl_mult_in_place(comb, rRiB);
+ rtnl_mult_in_place(comb, rCB);
+
+ rtnl_mtx_free(rRA);
+ rtnl_mtx_free(rP);
+ intmat_free(RiB);
+ rtnl_mtx_free(rRiB);
+ rtnl_mtx_free(rCB);
+
+ *pmb = comb;
+
+ }
+
+ } else {
+ match = NULL;
+ }
+
+ rtnl_mtx_free(CiA);
+ intmat_free(RA);
+ intmat_free(RB);
+ intmat_free(CB);
+ intmat_free(P);
+ cell_free(reference_reduced);
+ cell_free(cell_reduced);
+
+ return match;
}