diff options
| -rw-r--r-- | libcrystfel/src/reax.c | 159 |
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; + } } |