aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/reax.c
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src/reax.c')
-rw-r--r--libcrystfel/src/reax.c89
1 files changed, 67 insertions, 22 deletions
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 7529d12f..3460a4f2 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -3,11 +3,11 @@
*
* A new auto-indexer
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2011-2012 Thomas White <taw@physics.org>
+ * 2011-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -65,6 +65,9 @@
#define MAX_CANDIDATES (1024)
+/* Choose the best solution from this many candidate cells */
+#define MAX_REAX_CELL_CANDIDATES (32)
+
struct dvec
{
double x;
@@ -103,7 +106,7 @@ struct reax_search
struct reax_private
{
- IndexingPrivate base;
+ IndexingMethod indm;
struct dvec *directions;
int n_dir;
double angular_inc;
@@ -907,7 +910,7 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
}
- if ( cl->n_cand >= MAX_CELL_CANDIDATES ) {
+ if ( cl->n_cand >= MAX_REAX_CELL_CANDIDATES ) {
/* "cshift" just fell off the end of the list */
} else {
cl->cand[cl->n_cand++] = cshift;
@@ -915,15 +918,46 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
}
+static int check_cell(struct reax_private *dp, struct image *image,
+ UnitCell *cell)
+{
+ UnitCell *out;
+ Crystal *cr;
+
+ out = cell_new_from_cell(cell);
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to allocate crystal.\n");
+ return 0;
+ }
+
+ crystal_set_cell(cr, out);
+
+ if ( dp->indm & INDEXING_CHECK_PEAKS ) {
+ if ( !peak_sanity_check(image, &cr, 1) ) {
+ crystal_free(cr); /* Frees the cell as well */
+ return 0;
+ }
+ }
+
+ image_add_crystal(image, cr);
+
+ return 1;
+}
+
+
static void assemble_cells_from_candidates(struct image *image,
struct reax_search *s,
- UnitCell *cell)
+ UnitCell *cell,
+ struct reax_private *p)
{
int i, j, k;
signed int ti, tj, tk;
struct cell_candidate_list cl;
- cl.cand = calloc(MAX_CELL_CANDIDATES, sizeof(struct cell_candidate));
+ cl.cand = calloc(MAX_REAX_CELL_CANDIDATES,
+ sizeof(struct cell_candidate));
if ( cl.cand == NULL ) {
ERROR("Failed to allocate cell candidate list.\n");
return;
@@ -967,7 +1001,9 @@ static void assemble_cells_from_candidates(struct image *image,
continue;
}
- peak_lattice_agreement(image, cnew, &fom);
+ /* FIXME! */
+ //peak_lattice_agreement(image, cnew, &fom);
+ fom = 1.0;
add_cell_candidate(&cl, cnew, fom);
}
@@ -985,22 +1021,20 @@ static void assemble_cells_from_candidates(struct image *image,
cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga);
cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA,
&alA, &beA, &gaA);
- if ( (a - aA) > aA/10.0 ) w = 1;
- if ( (b - bA) > bA/10.0 ) w = 1;
- if ( (c - cA) > cA/10.0 ) w = 1;
- if ( (al - alA) > deg2rad(5.0) ) w = 1;
- if ( (be - beA) > deg2rad(5.0) ) w = 1;
- if ( (ga - gaA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(a - aA) > aA/10.0 ) w = 1;
+ if ( fabs(b - bA) > bA/10.0 ) w = 1;
+ if ( fabs(c - cA) > cA/10.0 ) w = 1;
+ if ( fabs(al - alA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(be - beA) > deg2rad(5.0) ) w = 1;
+ if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1;
if ( w ) {
STATUS("This cell is a long way from that sought:\n");
cell_print(cl.cand[i].cell);
}
}
- image->ncells = cl.n_cand;
- assert(image->ncells <= MAX_CELL_CANDIDATES);
for ( i=0; i<cl.n_cand; i++ ) {
- image->candidate_cells[i] = cl.cand[i].cell;
+ if ( check_cell(p, image, cl.cand[i].cell) ) break;
}
free(cl.cand);
@@ -1016,7 +1050,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
struct reax_search *s;
int i;
- assert(pp->indm == INDEXING_REAX);
p = (struct reax_private *)pp;
fft_in = fftw_malloc(p->nel*sizeof(double));
@@ -1039,7 +1072,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
// fft_in, fft_out, p->plan, smin, smax,
// image->det, p);
- assemble_cells_from_candidates(image, s, cell);
+ assemble_cells_from_candidates(image, s, cell, p);
for ( i=0; i<s->n_search; i++ ) {
free(s->search[i].cand);
@@ -1051,16 +1084,27 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
-IndexingPrivate *reax_prepare()
+IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell,
+ const char *filename, struct detector *det,
+ struct beam_params *beam, float *ltl)
{
struct reax_private *p;
int samp;
double th;
+ if ( cell == NULL ) {
+ ERROR("ReAx needs a unit cell.\n");
+ return NULL;
+ }
+
p = calloc(1, sizeof(*p));
if ( p == NULL ) return NULL;
- p->base.indm = INDEXING_REAX;
+ /* Flags that ReAx knows about */
+ *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS;
+
+ /* Flags that ReAx requires */
+ *indm |= INDEXING_USE_LATTICE_TYPE;
p->angular_inc = deg2rad(1.0);
@@ -1122,6 +1166,8 @@ IndexingPrivate *reax_prepare()
p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out,
1, FFTW_MEASURE);
+ p->indm = *indm;
+
return (IndexingPrivate *)p;
}
@@ -1130,7 +1176,6 @@ void reax_cleanup(IndexingPrivate *pp)
{
struct reax_private *p;
- assert(pp->indm == INDEXING_REAX);
p = (struct reax_private *)pp;
free(p->directions);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-15 05:09:54 +0000