aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/reax.c159
1 files changed, 154 insertions, 5 deletions
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 2217d6d9..523e8c94 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -45,6 +45,10 @@
#define INC_TOL_MULTIPLIER (3.0)
+/* Maximum number of cells that can be found by the reduction */
+#define MAX_CELLS (48)
+
+
struct dvec
{
double x;
@@ -229,9 +233,6 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
s->search[i].cand[cpos].fom = fom;
s->search[i].n_cand++;
- STATUS("Candidate %.2f %.2f %.2f %.2f sigma\n",
- dir->x, dir->y, dir->z, fom);
-
}
}
@@ -412,6 +413,16 @@ static void fine_search(struct reax_private *p, ImageFeatureList *flist,
}
+static int fom_compare(const void *av, const void *bv)
+{
+ const struct reax_candidate *a = av;
+ const struct reax_candidate *b = bv;
+
+ if ( a->fom > b->fom ) return -1;
+ return +1;
+}
+
+
static void squash_vectors(struct reax_search *s, double tol)
{
int i;
@@ -477,6 +488,9 @@ static void squash_vectors(struct reax_search *s, double tol)
sv->n_cand = n_copied;
sv->cand = new;
+ qsort(sv->cand, sv->n_cand, sizeof(sv->cand[0]),
+ fom_compare);
+
for ( j=0; j<sv->n_cand; j++ ) {
STATUS("%i: %+6.2f %+6.2f %+6.2f nm %.2f\n",
j, sv->cand[j].v.x*1e9,
@@ -792,12 +806,147 @@ static double max_feature_resolution(ImageFeatureList *flist)
}
+static int right_handed(struct rvec a, struct rvec b, struct rvec c)
+{
+ struct rvec aCb;
+ double aCb_dot_c;
+
+ /* "a" cross "b" */
+ aCb.u = a.v*b.w - a.w*b.v;
+ aCb.v = - (a.u*b.w - a.w*b.u);
+ aCb.w = a.u*b.v - a.v*b.u;
+
+ /* "a cross b" dot "c" */
+ aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
+
+ if ( aCb_dot_c > 0.0 ) return 1;
+ return 0;
+}
+
+
+static int check_twinning(UnitCell *c1, UnitCell *c2)
+{
+ int i;
+
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ for ( i=0; i<100; i++ ) {
+ signed int h, k, l;
+
+ h = flat_noise(0, 10);
+ k = flat_noise(0, 10);
+ l = flat_noise(0, 10);
+
+ }
+
+ return 0;
+}
+
+
+/* Return true if "cnew" accounts for more than 25% of the peaks predicted by
+ * any of the "ncells" cells in "cells". */
+static int twinned(UnitCell *cnew, UnitCell **cells, int ncells)
+{
+ int i;
+
+ for ( i=0; i<ncells; i++ ) {
+ if ( check_twinning(cnew, cells[i]) ) return 1;
+ }
+
+ return 0;
+}
+
+
+static int check_vector_combination(struct dvec *vi, struct dvec *vj,
+ struct dvec *vk, UnitCell *cell)
+{
+ double ang;
+ double a, b, c, al, be, ga;
+ const double angtol = deg2rad(2.0);
+
+ cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+
+ ang = angle_between(vi->x, vi->y, vi->z, vj->x, vj->y, vj->z);
+ if ( fabs(ang-ga) > angtol ) return 0;
+
+ ang = angle_between(vi->x, vi->y, vi->z,
+ vk->x, vk->y, vk->z);
+ if ( fabs(ang-be) > angtol ) return 0;
+
+ ang = angle_between(vj->x, vj->y, vj->z,
+ vk->x, vk->y, vk->z);
+ if ( fabs(ang-al) > angtol ) return 0;
+
+ return 1;
+}
+
+
static void assemble_cells_from_candidates(struct image *image,
struct reax_search *s,
UnitCell *cell)
{
- /* FIXME: Implement this */
- image->ncells = 0;
+ int i, j, k;
+ signed int ti, tj, tk;
+ UnitCell *cells[MAX_CELLS];
+ int ncells = 0;
+
+ /* Find candidates for axes 0 and 1 which have the right angle */
+ for ( i=0; i<s->search[0].n_cand; i++ ) {
+ for ( j=0; j<s->search[1].n_cand; j++ ) {
+ for ( k=0; k<s->search[2].n_cand; k++ ) {
+ for ( ti=-1; ti<=1; ti+=2 ) {
+ for ( tj=-1; tj<=1; tj+=2 ) {
+ for ( tk=-1; tk<=1; tk+=2 ) {
+
+ struct dvec vi, vj, vk;
+ struct rvec ai, bi, ci;
+ UnitCell *cnew;
+
+ vi = s->search[0].cand[i].v;
+ vj = s->search[1].cand[j].v;
+ vk = s->search[2].cand[k].v;
+
+ vi.x *= ti; vi.y *= ti; vi.z *= ti;
+ vj.x *= tj; vj.y *= tj; vj.z *= tj;
+ vk.x *= tk; vk.y *= tk; vk.z *= tk;
+
+ if ( !check_vector_combination(&vi, &vj, &vk, cell) ) continue;
+
+ ai.u = vi.x; ai.v = vi.y; ai.w = vi.z;
+ bi.u = vj.x; bi.v = vj.y; bi.w = vj.z;
+ ci.u = vk.x; ci.v = vk.y; ci.w = vk.z;
+
+ if ( !right_handed(ai, bi, ci) ) continue;
+
+ /* We have three vectors with the right angles */
+ cnew = cell_new_from_direct_axes(ai, bi, ci);
+
+ if ( twinned(cnew, cells, ncells) ) {
+ cell_free(cnew);
+ continue;
+ }
+
+ cells[ncells++] = cnew;
+
+ }
+ }
+ }
+ }
+ }
+ }
+
+ image->ncells = ncells;
+ for ( i=0; i<ncells; i++ ) {
+ image->candidate_cells[i] = cells[i];
+ }
+
+ if ( ncells > 0 ) {
+ image->indexed_cell = cells[0];
+ } else {
+ image->indexed_cell = NULL;
+ }
}