aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/cell-utils.c157
-rw-r--r--libcrystfel/src/cell-utils.h3
-rw-r--r--libcrystfel/src/rational.c25
-rw-r--r--libcrystfel/src/rational.h2
-rw-r--r--libcrystfel/src/xgandalf.c2
-rw-r--r--src/cell_tool.c17
-rw-r--r--tests/centering_check.c9
-rw-r--r--tests/transformation_check.c35
8 files changed, 191 insertions, 59 deletions
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, &centering);
+ template = uncenter_cell(template_in, &centering, 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, &centering_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, &centering_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; ia<ncand_a; ia++ ) {
for ( ib=0; ib<ncand_b; ib++ ) {
@@ -1982,8 +2048,6 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
/* Gamma OK, now look for place for c axis */
for ( ic=0; ic<ncand_c; ic++ ) {
- UnitCell *test2;
-
rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]);
rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]);
@@ -2010,16 +2074,19 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
cell_free(test);
continue;
}
- rtnl_mtx_print(m);
- test2 = cell_transform_intmat(test, centering_reference);
- cell_print(test2);
+
+ /* m is a rational matrix which turns Ap into Bp
+ * We need to return CB.M.CiA, which turns A into B */
+ rtnl_mtx_mtxmult(m, CiA, MCiA);
+ rtnl_mtx_mtxmult(CB, MCiA, CBMCiA);
+ *pmb = CBMCiA;
+ ret = 1;
+
cell_free(test);
- cell_free(test2);
}
}
}
- *pmb = m;
- return 1;
+ return ret;
}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 0face43d..a97f2bfb 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -66,7 +66,8 @@ extern int cell_is_sensible(UnitCell *cell);
extern int validate_cell(UnitCell *cell);
-extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **m);
+extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC,
+ RationalMatrix **pCi);
extern int bravais_lattice(UnitCell *cell);
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index 58178559..4f02c8b1 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -546,6 +546,31 @@ void rtnl_mtx_print(const RationalMatrix *m)
}
+void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
+ RationalMatrix *ans)
+{
+ int i, j;
+
+ assert(ans->cols == A->cols);
+ assert(ans->rows == B->rows);
+ assert(A->cols == B->rows);
+
+ for ( i=0; i<ans->rows; i++ ) {
+ for ( j=0; j<ans->cols; j++ ) {
+ int k;
+ Rational sum = rtnl_zero();
+ for ( k=0; k<A->rows; 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);
diff --git a/src/cell_tool.c b/src/cell_tool.c
index d6c5b6cd..6d0c5f1e 100644
--- a/src/cell_tool.c
+++ b/src/cell_tool.c
@@ -92,7 +92,14 @@ static int comparecells(UnitCell *cell, const char *comparecell,
STATUS("No relationship found between lattices.\n");
return 0;
} else {
+ UnitCell *trans;
STATUS("Relationship found:\n");
+ rtnl_mtx_print(m);
+ STATUS("Transformed version of input unit cell:\n");
+ trans = cell_transform_rational(cell, m);
+ cell_print(trans);
+ cell_free(trans);
+
}
rtnl_mtx_free(m);
@@ -294,15 +301,19 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
static int uncenter(UnitCell *cell, const char *out_file)
{
UnitCell *cnew;
- IntegerMatrix *m;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
- cnew = uncenter_cell(cell, &m);
+ cnew = uncenter_cell(cell, &C, &Ci);
STATUS("------------------> The primitive unit cell:\n");
cell_print(cnew);
STATUS("------------------> The centering transformation:\n");
- intmat_print(m);
+ intmat_print(C);
+
+ STATUS("------------------> The un-centering transformation:\n");
+ rtnl_mtx_print(Ci);
if ( out_file != NULL ) {
FILE *fh = fopen(out_file, "w");
diff --git a/tests/centering_check.c b/tests/centering_check.c
index 774af66b..cc553dd0 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -61,7 +61,8 @@ static int check_centering(double a, double b, double c,
{
UnitCell *cell, *cref;
UnitCell *n;
- IntegerMatrix *t;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
int fail = 0;
int i;
double asx, asy, asz;
@@ -86,12 +87,12 @@ static int check_centering(double a, double b, double c,
cell_free(cref);
check_cell(cell, "Input");
- n = uncenter_cell(cell, &t);
+ n = uncenter_cell(cell, &C, &Ci);
if ( n != NULL ) {
STATUS("The back transformation is:\n");
- intmat_print(t);
+ intmat_print(C);
if ( check_cell(n, "Output") ) fail = 1;
- if ( intmat_is_identity(t) && ((cen=='P') || (cen=='R')) ) {
+ if ( intmat_is_identity(C) && ((cen=='P') || (cen=='R')) ) {
STATUS("Identity, as expected.\n");
} else {
cell_print(n);
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index c5dad16c..06521061 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -224,25 +224,50 @@ static int check_transformation(UnitCell *cell, IntegerMatrix *tfn,
static int check_uncentering(UnitCell *cell)
{
UnitCell *ct;
- IntegerMatrix *tr;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
+ UnitCell *cback;
+ double a[9], b[9];
+ int i;
int fail = 0;
STATUS("-----------------------\n");
STATUS("----> Before transformation:\n");
- cell_print(cell);
+ cell_print_full(cell);
- ct = uncenter_cell(cell, &tr);
+ ct = uncenter_cell(cell, &C, &Ci);
if ( ct == NULL ) return 1;
STATUS("----> The primitive unit cell:\n");
cell_print(ct);
STATUS("----> The matrix to put the centering back:\n");
- intmat_print(tr);
+ intmat_print(C);
+
+ STATUS("----> The recovered centered cell:\n");
+ cback = cell_transform_intmat(ct, C);
+ cell_print(cback);
+
+ cell_get_cartesian(cell, &a[0], &a[1], &a[2],
+ &a[3], &a[4], &a[5],
+ &a[6], &a[7], &a[8]);
+ cell_get_cartesian(cback, &b[0], &b[1], &b[2],
+ &b[3], &b[4], &b[5],
+ &b[6], &b[7], &b[8]);
+ for ( i=0; i<9; i++ ) {
+ if ( fabs(a[i] - b[i]) > 1e-12 ) {
+ fail = 1;
+ }
+ }
+
+ if ( fail ) {
+ ERROR("Original cell not recovered after back transformation\n");
+ }
- fail = check_same_reflections(cell, ct);
+ fail += check_same_reflections(cell, ct);
cell_free(ct);
+ cell_free(cback);
return fail;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-20 02:10:25 +0000