aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/cell-utils.c230
-rw-r--r--libcrystfel/src/cell-utils.h4
-rw-r--r--libcrystfel/src/rational.c33
-rw-r--r--libcrystfel/src/rational.h4
-rw-r--r--src/cell_tool.c80
5 files changed, 278 insertions, 73 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index a21ae570..2bfaa3e1 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -1793,3 +1793,233 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
intmat_free(m);
return 0;
}
+
+
+struct cand
+{
+ Rational abc[3];
+ double fom;
+};
+
+
+static int cmpcand(const void *av, const void *bv)
+{
+ const struct cand *a = av;
+ const struct cand *b = bv;
+ return a->fom > b->fom;
+}
+
+
+static Rational *find_candidates(double len, double *a, double *b, double *c,
+ double ltl, int *pncand)
+{
+ Rational *r;
+ struct cand *cands;
+ const int max_cand = 1024;
+ int ncand = 0;
+ Rational *rat;
+ int nrat;
+ int nrej = 0;
+ int ia, ib, ic;
+ int i;
+
+ cands = malloc(max_cand * sizeof(struct cand));
+ if ( cands == NULL ) return NULL;
+
+ rat = rtnl_list(0, 4, -4, 4, &nrat);
+ if ( rat == NULL ) return NULL;
+
+ for ( ia=0; ia<nrat; ia++ ) {
+ for ( ib=0; ib<nrat; ib++ ) {
+ for ( ic=0; ic<nrat; ic++ ) {
+ double vec[3];
+ double abc[3];
+ double veclen;
+ abc[0] = rtnl_as_double(rat[ia]);
+ abc[1] = rtnl_as_double(rat[ib]);
+ abc[2] = rtnl_as_double(rat[ic]);
+ vec[0] = a[0]*abc[0] + b[0]*abc[1] + c[0]*abc[2];
+ vec[1] = a[1]*abc[0] + b[1]*abc[1] + c[1]*abc[2];
+ vec[2] = a[2]*abc[0] + b[2]*abc[1] + c[2]*abc[2];
+ veclen = modulus(vec[0], vec[1], vec[2]);
+ if ( within_tolerance(len, veclen, ltl*100.0) ) {
+ if ( ncand == max_cand ) {
+ nrej++;
+ } else {
+ cands[ncand].abc[0] = rat[ia];
+ cands[ncand].abc[1] = rat[ib];
+ cands[ncand].abc[2] = rat[ic];
+ cands[ncand].fom = fabs(veclen - len);
+ ncand++;
+ }
+ }
+ }
+ }
+ }
+
+ if ( nrej ) {
+ ERROR("WARNING: Too many vector candidates (%i rejected)\n", nrej);
+ }
+
+ /* Sort by difference from reference vector length */
+ qsort(cands, ncand, sizeof(struct cand), cmpcand);
+
+ r = malloc(ncand * 3 * sizeof(Rational));
+ if ( r == 0 ) return NULL;
+
+ for ( i=0; i<ncand; i++ ) {
+ r[3*i+0] = cands[i].abc[0];
+ r[3*i+1] = cands[i].abc[1];
+ r[3*i+2] = cands[i].abc[2];
+ }
+ free(cands);
+
+ *pncand = ncand;
+ return r;
+}
+
+
+/**
+ * compare_reindexed_cell_parameters:
+ * @cell_in: A %UnitCell
+ * @reference_in: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in direct-space axis lengths
+ * @atl: Maximum allowable difference in direct-space angles (in radians)
+ * @pmb: Place to store pointer to matrix
+ *
+ * Compare the @cell_in with @reference_in. If @cell is a derivative lattice
+ * of @reference, within fractional axis length difference @ltl and absolute angle
+ * difference @atl (in radians), this function returns non-zero and stores the
+ * transformation which needs to be applied to @cell_in at @pmb.
+ *
+ * Only the cell parameters will be compared. The relative orientations are
+ * irrelevant.
+ *
+ * Returns: non-zero if the cells match, zero for no match or error.
+ *
+ */
+int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
+ double ltl, double atl,
+ RationalMatrix **pmb)
+{
+ UnitCell *cell;
+ UnitCell *reference;
+ IntegerMatrix *centering_reference;
+ IntegerMatrix *centering_cell;
+ RationalMatrix *m;
+ double a, b, c, al, be, ga;
+ double av[3], bv[3], cv[3];
+ Rational *cand_a;
+ Rational *cand_b;
+ Rational *cand_c;
+ int ncand_a, ncand_b, ncand_c;
+ int i;
+ int ia, ib;
+
+ /* Actually compare against primitive version of reference */
+ reference = uncenter_cell(reference_in, &centering_reference);
+ if ( reference == NULL ) return 0;
+
+ /* Actually compare primitive version of cell */
+ cell = uncenter_cell(cell_in, &centering_cell);
+ if ( cell == NULL ) return 0;
+
+ /* Get target parameters */
+ cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+ cell_get_cartesian(cell, &av[0], &av[1], &av[2],
+ &bv[0], &bv[1], &bv[2],
+ &cv[0], &cv[1], &cv[2]);
+
+ /* Find vectors in 'cell' with lengths close to a, b and c */
+ cand_a = find_candidates(a, av, bv, cv, ltl, &ncand_a);
+ cand_b = find_candidates(b, av, bv, cv, ltl, &ncand_b);
+ cand_c = find_candidates(c, av, bv, cv, ltl, &ncand_c);
+
+ STATUS("Candidates for a: %i\n", ncand_a);
+ for ( i=0; i<10; i++ ) {
+ STATUS("%s %s %s\n", rtnl_format(cand_a[3*i+0]),
+ rtnl_format(cand_a[3*i+1]),
+ rtnl_format(cand_a[3*i+2]));
+ }
+ STATUS("Candidates for b: %i\n", ncand_b);
+ for ( i=0; i<10; i++ ) {
+ STATUS("%s %s %s\n", rtnl_format(cand_b[3*i+0]),
+ rtnl_format(cand_b[3*i+1]),
+ rtnl_format(cand_b[3*i+2]));
+ }
+ STATUS("Candidates for c: %i\n", ncand_c);
+ for ( i=0; i<10; i++ ) {
+ STATUS("%s %s %s\n", rtnl_format(cand_c[3*i+0]),
+ rtnl_format(cand_c[3*i+1]),
+ rtnl_format(cand_c[3*i+2]));
+ }
+
+ m = rtnl_mtx_new(3, 3);
+ for ( ia=0; ia<ncand_a; ia++ ) {
+ for ( ib=0; ib<ncand_b; ib++ ) {
+
+ UnitCell *test;
+ double at, bt, ct, alt, bet, gat;
+ int ic = 0;
+
+ /* Form the matrix using the first candidate for c */
+ 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]);
+ rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+ rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+ rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
+ rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
+ rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+ rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+ /* Check angle between a and b */
+ test = cell_transform_rational(cell, m);
+ cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+ cell_free(test);
+ if ( fabs(gat - ga) > atl ) continue;
+
+ /* 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]);
+ rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+ rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+ rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
+ rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
+ rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+ rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+ if ( rtnl_cmp(rtnl_mtx_det(m),rtnl_zero()) == 0 ) continue;
+
+ test = cell_transform_rational(cell, m);
+ cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+ if ( !right_handed(test) ) {
+ cell_free(test);
+ continue;
+ }
+ if ( fabs(alt - al) > atl ) {
+ cell_free(test);
+ continue;
+ }
+ if ( fabs(bet - be) > atl ) {
+ cell_free(test);
+ continue;
+ }
+ rtnl_mtx_print(m);
+ test2 = cell_transform_intmat(test, centering_reference);
+ cell_print(test2);
+ cell_free(test);
+ cell_free(test2);
+
+ }
+ }
+ }
+
+ *pmb = m;
+ return 1;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index a1a6fd46..0face43d 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -95,6 +95,10 @@ extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a,
double atl,
IntegerMatrix **pmb);
+extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference,
+ double ltl, double atl,
+ RationalMatrix **pmb);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index 5db83165..0f2bcaee 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -226,6 +226,39 @@ char *rtnl_format(Rational rt)
}
+Rational *rtnl_list(signed int num_min, signed int num_max,
+ signed int den_min, signed int den_max,
+ int *pn)
+{
+ signed int num, den;
+ Rational *list;
+ int n = 0;
+
+ list = malloc((1+num_max-num_min)*(1+den_max-den_min)*sizeof(Rational));
+ if ( list == NULL ) return NULL;
+
+ for ( num=num_min; num<=num_max; num++ ) {
+ for ( den=den_min; den<=den_max; den++ ) {
+
+ Rational r = rtnl(num, den);
+
+ /* Denominator zero? */
+ if ( den == 0 ) continue;
+
+ /* Same as last entry? */
+ if ( (n>0) && (rtnl_cmp(list[n-1], r)==0) ) continue;
+
+ /* Can be reduced? */
+ if ( gcd(num, den) != 1 ) continue;
+
+ list[n++] = r;
+ }
+ }
+ *pn = n;
+ return list;
+}
+
+
/**
* SECTION:rational_matrix
* @short_description: Rational matrices
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
index 925f9be1..7983f56f 100644
--- a/libcrystfel/src/rational.h
+++ b/libcrystfel/src/rational.h
@@ -77,6 +77,10 @@ extern Rational rtnl_abs(Rational a);
extern char *rtnl_format(Rational rt);
+extern Rational *rtnl_list(signed int num_min, signed int num_max,
+ signed int den_min, signed int den_max,
+ int *pn);
+
extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols);
extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m);
extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j);
diff --git a/src/cell_tool.c b/src/cell_tool.c
index 03c58f5e..d6c5b6cd 100644
--- a/src/cell_tool.c
+++ b/src/cell_tool.c
@@ -72,14 +72,9 @@ static void show_help(const char *s)
static int comparecells(UnitCell *cell, const char *comparecell,
double ltl, double atl)
{
- signed int i[9];
- int b[9];
- const int maxorder = 2;
UnitCell *cell2;
RationalMatrix *m;
- STATUS("Comparing with: %s\n", comparecell);
-
cell2 = load_cell_from_file(comparecell);
if ( cell2 == NULL ) {
ERROR("Failed to load unit cell from '%s'\n", comparecell);
@@ -92,75 +87,14 @@ static int comparecells(UnitCell *cell, const char *comparecell,
STATUS("------------------> The reference unit cell:\n");
cell_print(cell2);
- STATUS("Comparing cells up to %ix each lattice length.\n", maxorder);
- STATUS("Reciprocal axis length tolerance %f %%\n", ltl);
- STATUS("Reciprocal angle tolerance %f degrees\n", rad2deg(atl));
- STATUS("This will take about 30 seconds. Please wait...\n");
-
- m = rtnl_mtx_new(3, 3);
- for ( i[0]=-maxorder; i[0]<=+maxorder; i[0]++ ) {
- for ( i[1]=-maxorder; i[1]<=+maxorder; i[1]++ ) {
- for ( i[2]=-maxorder; i[2]<=+maxorder; i[2]++ ) {
- for ( i[3]=-maxorder; i[3]<=+maxorder; i[3]++ ) {
- for ( i[4]=-maxorder; i[4]<=+maxorder; i[4]++ ) {
- for ( i[5]=-maxorder; i[5]<=+maxorder; i[5]++ ) {
- for ( i[6]=-maxorder; i[6]<=+maxorder; i[6]++ ) {
- for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) {
- for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) {
- for ( b[0]=0; b[0]<=1; b[0]++ ) {
- for ( b[1]=0; b[1]<=1; b[1]++ ) {
- for ( b[2]=0; b[2]<=1; b[2]++ ) {
- for ( b[3]=0; b[3]<=1; b[3]++ ) {
- for ( b[4]=0; b[4]<=1; b[4]++ ) {
- for ( b[5]=0; b[5]<=1; b[5]++ ) {
- for ( b[6]=0; b[6]<=1; b[6]++ ) {
- for ( b[7]=0; b[7]<=1; b[7]++ ) {
- for ( b[8]=0; b[8]<=1; b[8]++ ) {
-
- UnitCell *nc;
- int j, k;
- int l = 0;
-
- for ( j=0; j<3; j++ ) {
- for ( k=0; k<3; k++ ) {
- if ( b[l] || i[l]==1 || i[l]==-1 ) {
- rtnl_mtx_set(m, j, k, rtnl(i[l], 1));
- } else {
- rtnl_mtx_set(m, j, k, rtnl(1, i[l]));
- }
- l++;
- }
- }
-
- nc = cell_transform_rational(cell, m);
-
- if ( compare_cell_parameters(cell2, nc, ltl, atl) ) {
- STATUS("-----------------------------------------------"
- "-------------------------------------------\n");
- cell_print(nc);
- rtnl_mtx_print(m);
- }
-
- cell_free(nc);
-
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
- }
+ STATUS("------------------> The comparison results:\n");
+ if ( !compare_reindexed_cell_parameters(cell, cell2, ltl, atl, &m) ) {
+ STATUS("No relationship found between lattices.\n");
+ return 0;
+ } else {
+ STATUS("Relationship found:\n");
}
+
rtnl_mtx_free(m);
return 0;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 12:28:35 +0000