aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/peaks.c24
-rw-r--r--libcrystfel/src/peaks.h5
-rw-r--r--libcrystfel/src/reax.c165
3 files changed, 162 insertions, 32 deletions
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 0cbce55f..713dd0a6 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -473,9 +473,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
}
-int peak_sanity_check(struct image *image)
+double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst)
{
-
int i;
int n_feat = 0;
int n_sane = 0;
@@ -483,14 +482,13 @@ int peak_sanity_check(struct image *image)
double bx, by, bz;
double cx, cy, cz;
double min_dist = 0.25;
+ double stot = 0.0;
/* Round towards nearest */
fesetround(1);
/* Cell basis vectors for this image */
- cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
/* Loop over peaks, checking proximity to nearest reflection */
for ( i=0; i<image_feature_count(image->features); i++ ) {
@@ -520,17 +518,25 @@ int peak_sanity_check(struct image *image)
if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist)
&& (fabs(l - ld) < min_dist) )
{
+ double sval;
n_sane++;
+ sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
+ stot += 1.0 - sval;
continue;
}
}
- /* return 0 means fail test, return 1 means pass test */
- // printf("%d out of %d peaks are \"sane\"\n",n_sane,n_feat);
- if ( (float)n_sane / (float)n_feat < 0.5 ) return 0;
+ *pst = stot;
+ return (double)n_sane / (float)n_feat;
+}
+
- return 1;
+int peak_sanity_check(struct image *image)
+{
+ double stot;
+ /* 0 means failed test, 1 means passed test */
+ return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5;
}
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index d52c75e4..ffeb896b 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -28,7 +28,10 @@ extern void integrate_reflections(struct image *image,
int polar, int use_closer, int bgsub,
double min_snr);
-extern int peak_sanity_check(struct image * image);
+extern double peak_lattice_agreement(struct image *image, UnitCell *cell,
+ double *pst);
+
+extern int peak_sanity_check(struct image *image);
/* Exported so it can be poked by integration_check */
extern int integrate_peak(struct image *image, int cfs, int css,
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 80353a2c..ea8f143a 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -838,35 +838,97 @@ static int right_handed(struct rvec a, struct rvec b, struct rvec c)
}
-static int check_twinning(UnitCell *c1, UnitCell *c2)
+struct cell_candidate
+{
+ UnitCell *cell;
+ double fom;
+};
+
+
+struct cell_candidate_list
+{
+ struct cell_candidate *cand;
+ int n_cand;
+};
+
+
+static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
{
int i;
+ int n_dup;
+ const int n_trials = 40;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+
+ cell_get_reciprocal(c1, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ cell_get_cartesian(c2, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ n_dup = 0;
+ for ( i=0; i<n_trials; i++ ) {
- for ( i=0; i<100; i++ ) {
signed int h, k, l;
+ double h2, k2, l2;
+ double rx, ry, rz;
+ double dev;
+ signed int h2i, k2i, l2i;
h = flat_noise(0, 10);
k = flat_noise(0, 10);
l = flat_noise(0, 10);
+ /* Position of this (randomly selected)
+ * reciprocal lattice point */
+ rx = h*asx + k*bsx + l*csx;
+ ry = h*asy + k*bsy + l*csy;
+ rz = h*asz + k*bsz + l*csz;
+
+ /* Indices of this point in the basis of the other cell */
+ h2 = rx*ax + ry*ay + rz*az;
+ k2 = rx*bx + ry*by + rz*bz;
+ l2 = rx*cx + ry*cy + rz*cz;
+
+ h2i = lrint(h2);
+ k2i = lrint(k2);
+ l2i = lrint(l2);
+
+ dev = pow(h2i-h2, 2.0) + pow(k2i-k2, 2.0) + pow(l2i-l2, 2.0);
+
+ if ( verbose ) {
+ STATUS("%3i %3i %3i -> %5.2f %5.2f %5.2f -> "
+ "%3i %3i %3i -> %5.2f\n", h, k, l,
+ h2, k2, l2, h2i, k2i, l2i, dev);
+ }
+
+ if ( dev < 0.1 ) {
+ n_dup++;
+ }
+
+ }
+
+ if ( verbose ) {
+ STATUS("%i duplicates.\n", n_dup);
}
+ if ( n_dup > 10 ) return 1;
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)
+static int twinned(UnitCell *cnew, struct cell_candidate_list *cl)
{
int i;
- for ( i=0; i<ncells; i++ ) {
- if ( check_twinning(cnew, cells[i]) ) return 1;
+ for ( i=0; i<cl->n_cand; i++ ) {
+ if ( check_twinning(cnew, cl->cand[i].cell, 0) ) return 1;
}
return 0;
@@ -885,26 +947,64 @@ static int check_vector_combination(struct dvec *vi, struct dvec *vj,
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);
+ 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);
+ 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 add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
+ double fom)
+{
+ struct cell_candidate cshift;
+ int i, cpos;
+
+ cpos = cl->n_cand;
+ for ( i=0; i<cl->n_cand; i++ ) {
+ if ( fom > cl->cand[i].fom ) {
+ cpos = i;
+ break;
+ }
+ }
+
+ cshift.cell = cnew;
+ cshift.fom = fom;
+
+ for ( i=cpos; i<cl->n_cand; i++ ) {
+
+ struct cell_candidate cshift2;
+ cshift2 = cl->cand[i];
+ cl->cand[i] = cshift;
+ cshift = cshift2;
+
+ }
+
+ if ( cl->n_cand >= MAX_CELL_CANDIDATES ) {
+ /* "cshift" just fell off the end of the list */
+ } else {
+ cl->cand[cl->n_cand++] = cshift;
+ }
+}
+
+
static void assemble_cells_from_candidates(struct image *image,
struct reax_search *s,
UnitCell *cell)
{
int i, j, k;
signed int ti, tj, tk;
- UnitCell *cells[MAX_CELL_CANDIDATES];
- int ncells = 0;
+ struct cell_candidate_list cl;
+
+ cl.cand = calloc(MAX_CELL_CANDIDATES, sizeof(struct cell_candidate));
+ if ( cl.cand == NULL ) {
+ ERROR("Failed to allocate cell candidate list.\n");
+ return;
+ }
+ cl.n_cand = 0;
/* Find candidates for axes 0 and 1 which have the right angle */
for ( i=0; i<s->search[0].n_cand; i++ ) {
@@ -917,6 +1017,7 @@ static void assemble_cells_from_candidates(struct image *image,
struct dvec vi, vj, vk;
struct rvec ai, bi, ci;
UnitCell *cnew;
+ double fom;
vi = s->search[0].cand[i].v;
vj = s->search[1].cand[j].v;
@@ -937,17 +1038,15 @@ static void assemble_cells_from_candidates(struct image *image,
/* We have three vectors with the right angles */
cnew = cell_new_from_direct_axes(ai, bi, ci);
- if ( twinned(cnew, cells, ncells) ) {
+ refine_cell(image, cnew, image->features);
+
+ if ( twinned(cnew, &cl) ) {
cell_free(cnew);
continue;
}
- refine_cell(image, cnew, image->features);
-
- /* FIXME: Rank according to quality of prediction */
- if ( ncells < MAX_CELL_CANDIDATES ) {
- cells[ncells++] = cnew;
- }
+ peak_lattice_agreement(image, cnew, &fom);
+ add_cell_candidate(&cl, cnew, fom);
}
}
@@ -956,11 +1055,33 @@ static void assemble_cells_from_candidates(struct image *image,
}
}
- image->ncells = ncells;
- assert(ncells <= MAX_CELL_CANDIDATES);
- for ( i=0; i<ncells; i++ ) {
- image->candidate_cells[i] = cells[i];
+ for ( i=0; i<cl.n_cand; i++ ) {
+ double a, b, c, al, be, ga;
+ double aA, bA, cA, alA, beA, gaA;
+ int w = 0;
+// STATUS("%i: %f\n", i, cl.cand[i].fom);
+ 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 ( 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;
+ }
+
+ free(cl.cand);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 04:25:54 +0000