aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-09-20 13:41:16 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:42 +0100
commita88d3ef8cd02616ecdb218e1d6c7eecac874af98 (patch)
tree85d1fcd01583a526aa2856ff29d7ec0974d1d9d1
parent636b71d41e8b2cd318e6df24159a6fb65301e5df (diff)
More ReAx work
-rw-r--r--libcrystfel/src/reax.c291
1 files changed, 143 insertions, 148 deletions
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 49231438..4387d039 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -35,6 +35,16 @@
#include "index-priv.h"
+/* Minimum number of standard deviations above the mean a peak must be in the
+ * 1D FT to qualify as a candidate vector */
+#define MIN_SIGMAS (7.0)
+
+
+/* Maximum number of times the angular tolerance that vectors are permitted to
+ * be together before they get merged by squash_vectors() */
+#define INC_TOL_MULTIPLIER (3.0)
+
+
struct dvec
{
double x;
@@ -48,7 +58,7 @@ struct dvec
struct reax_candidate
{
struct dvec v; /* This is the vector for the candidate */
- double strength; /* How much intensity fell into the range */
+ double fom;
};
@@ -151,19 +161,34 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
double tot = 0.0;
double peak = 0.0;
double peak_mod = 0.0;
+ double mean;
+ double sd = 0.0;
int j;
+ int n;
- for ( j=s->search[i].smin; j<=s->search[i].smax; j++ ) {
+ for ( j=0; j<nel/2+1; j++ ) {
- double re, im;
+ double re, im, am;
re = fft_out[j][0];
im = fft_out[j][1];
- peak += sqrt(re*re + im*im);
+ am = sqrt(re*re + im*im);
+
+ tot += am;
+ n++;
+
+ if ( ( j >= s->search[i].smin )
+ && ( j <= s->search[i].smax ) ) {
+ if ( am > peak ) {
+ peak = am;
+ peak_mod = (double)j/(2.0*pmax);
+ }
+ }
}
+ mean = tot/(double)n;
- for ( j=0; j<nel; j++ ) {
+ for ( j=0; j<nel/2+1; j++ ) {
double re, im, am;
@@ -171,19 +196,13 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
im = fft_out[j][1];
am = sqrt(re*re + im*im);
- tot += am;
- if ( (j>=s->search[i].smin)
- && (j<=s->search[i].smax) ) {
- if ( am > peak ) {
- peak = am;
- peak_mod = (double)j/(2.0*pmax);
- }
- }
+ sd += pow(am - mean, 2.0);
}
+ sd = sqrt(sd/(double)n);
/* If sufficiently strong, add to list of candidates */
- if ( peak > 2.0*(tot/nel) ) {
+ if ( peak > mean+MIN_SIGMAS*sd ) {
size_t ns;
struct reax_candidate *cnew;
@@ -197,14 +216,22 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
ERROR("Failed to add candidate.\n");
} else {
+ double fom;
+
+ fom = (peak-mean)/sd;
+
s->search[i].cand = cnew;
s->search[i].cand[cpos].v.x = dir->x*peak_mod;
s->search[i].cand[cpos].v.y = dir->y*peak_mod;
s->search[i].cand[cpos].v.z = dir->z*peak_mod;
s->search[i].cand[cpos].v.th = dir->th;
s->search[i].cand[cpos].v.ph = dir->ph;
+ 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);
+
}
}
@@ -266,15 +293,95 @@ static void fine_search(struct reax_private *p, ImageFeatureList *flist,
}
-static double find_candidates(struct reax_private *p,
- ImageFeatureList *flist, double pmax,
- double *fft_in, fftw_complex *fft_out,
- struct reax_search *s,
- const char *rg, struct detector *det)
+static void squash_vectors(struct reax_search *s, double tol)
+{
+ int i;
+
+ for ( i=0; i<s->n_search; i++ ) {
+
+ struct reax_search_v *sv;
+ struct reax_candidate *new;
+ int j, k;
+ int n_invalid = 0;
+ int n_copied;
+
+ sv = &s->search[i];
+
+ for ( j=0; j<sv->n_cand; j++ ) {
+ for ( k=0; k<sv->n_cand; k++ ) {
+
+ struct reax_candidate *v1, *v2;
+
+ if ( j == k ) continue;
+
+ v1 = &sv->cand[j];
+ v2 = &sv->cand[k];
+
+ if ( angle_between(v1->v.x, v1->v.y, v1->v.z,
+ v2->v.x, v2->v.y, v2->v.z) < tol )
+ {
+ if ( !isnan(v1->fom) && !isnan(v2->fom ) ) {
+ if ( v1->fom > v2->fom ) {
+ v2->fom = NAN;
+ } else {
+ v1->fom = NAN;
+ }
+ n_invalid++;
+ }
+ }
+
+ }
+ }
+
+ new = calloc(sv->n_cand-n_invalid,
+ sizeof(struct reax_candidate));
+ if ( new == NULL ) {
+ ERROR("Failed to allocate memory for squashed"
+ " candidate list.\n");
+ return;
+ }
+
+ n_copied = 0;
+ for ( j=0; j<sv->n_cand; j++ ) {
+ if ( !isnan(sv->cand[j].fom) ) {
+
+ new[n_copied] = sv->cand[j];
+ n_copied++;
+
+ }
+ }
+ assert(sv->n_cand - n_invalid == n_copied);
+
+ free(sv->cand);
+ STATUS("Search vector %i: squashed %i candidates down to %i\n",
+ i, sv->n_cand, n_copied);
+ sv->n_cand = n_copied;
+ sv->cand = new;
+
+ for ( j=0; j<sv->n_cand; j++ ) {
+ STATUS("%i: %+6.2f %+6.2f %+6.2f %.2f\n",
+ j, sv->cand[j].v.x, sv->cand[j].v.y,
+ sv->cand[j].v.z, sv->cand[j].fom);
+ }
+
+ }
+}
+
+
+static void find_candidates(struct reax_private *p,
+ ImageFeatureList *flist, double pmax,
+ double *fft_in, fftw_complex *fft_out,
+ struct reax_search *s,
+ const char *rg, struct detector *det)
{
- double fom = 0.0;
int i;
double th, ph;
+ double fom;
+
+ for ( i=0; i<s->n_search; i++ ) {
+ s->search[i].cand = NULL;
+ s->search[i].n_cand = 0;
+ }
fom = 0.0; th = 0.0; ph = 0.0;
for ( i=0; i<p->n_dir; i++ ) {
@@ -292,19 +399,21 @@ static double find_candidates(struct reax_private *p,
}
+ squash_vectors(s, INC_TOL_MULTIPLIER*p->angular_inc);
+
for ( i=0; i<s->n_search; i++ ) {
struct reax_search_v *sv;
int j;
sv = &s->search[i];
+ STATUS("Search %i: doing fine search for %i candidates\n",
+ i, sv->n_cand);
for ( j=0; j<sv->n_cand; j++ ) {
fine_search(p, flist, pmax, fft_in, fft_out, sv,
&sv->cand[j], rg, det);
}
}
-
- return fom;
}
@@ -337,7 +446,8 @@ static struct reax_search *search_all_axes(UnitCell *cell, double pmax)
s = malloc(3*sizeof(*s));
s->pmax = pmax;
- s->n_search = 1;
+ s->n_search = 3;
+ s->search = malloc(3*sizeof(struct reax_search_v));
smin = 2.0*pmax * amin; smax = 2.0*pmax * amax;
s->search[0].smin = smin; s->search[0].smax = smax;
smin = 2.0*pmax * bmin; smax = 2.0*pmax * bmax;
@@ -349,125 +459,6 @@ static struct reax_search *search_all_axes(UnitCell *cell, double pmax)
}
-/* Refine a direct space vector. From Clegg (1984) */
-static double iterate_refine_vector(double *x, double *y, double *z,
- ImageFeatureList *flist)
-{
- int fi, n, err;
- gsl_matrix *C;
- gsl_vector *A;
- gsl_vector *t;
- gsl_matrix *s_vec;
- gsl_vector *s_val;
- double dtmax;
-
- A = gsl_vector_calloc(3);
- C = gsl_matrix_calloc(3, 3);
-
- n = image_feature_count(flist);
- fesetround(1);
- for ( fi=0; fi<n; fi++ ) {
-
- struct imagefeature *f;
- double val;
- double kn, kno;
- double xv[3];
- int i, j;
-
- f = image_get_feature(flist, fi);
- if ( f == NULL ) continue;
-
- kno = f->rx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */
- kn = nearbyint(kno);
- if ( kn - kno > 0.3 ) continue;
-
- xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz;
-
- for ( i=0; i<3; i++ ) {
-
- val = gsl_vector_get(A, i);
- gsl_vector_set(A, i, val+xv[i]*kn);
-
- for ( j=0; j<3; j++ ) {
- val = gsl_matrix_get(C, i, j);
- gsl_matrix_set(C, i, j, val+xv[i]*xv[j]);
- }
-
- }
-
- }
-
- s_val = gsl_vector_calloc(3);
- s_vec = gsl_matrix_calloc(3, 3);
- err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val);
- if ( err ) {
- ERROR("SVD failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_matrix_free(C);
- gsl_vector_free(A);
- return 0.0;
- }
-
- t = gsl_vector_calloc(3);
- err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t);
- if ( err ) {
- ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_matrix_free(C);
- gsl_vector_free(A);
- gsl_vector_free(t);
- return 0.0;
- }
-
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
-
- dtmax = fabs(*x - gsl_vector_get(t, 0));
- dtmax += fabs(*y - gsl_vector_get(t, 1));
- dtmax += fabs(*z - gsl_vector_get(t, 2));
-
- *x = gsl_vector_get(t, 0);
- *y = gsl_vector_get(t, 1);
- *z = gsl_vector_get(t, 2);
-
- gsl_matrix_free(C);
- gsl_vector_free(A);
-
- return dtmax;
-}
-
-
-static void refine_cell(struct image *image, UnitCell *cell,
- ImageFeatureList *flist)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- int i;
- double sm;
-
- cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- i = 0;
- do {
-
- sm = iterate_refine_vector(&ax, &ay, &az, flist);
- sm += iterate_refine_vector(&bx, &by, &bz, flist);
- sm += iterate_refine_vector(&cx, &cy, &cz, flist);
- i++;
-
- } while ( (sm > 0.001e-9) && (i<10) );
-
- cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
-
- if ( i == 10 ) {
- cell_free(image->indexed_cell);
- image->indexed_cell = NULL;
- }
-}
-
-
static double get_model_phase(double x, double y, double z, ImageFeatureList *f,
int nel, double pmax, double *fft_in,
fftw_complex *fft_out, fftw_plan plan,
@@ -673,6 +664,13 @@ static double max_feature_resolution(ImageFeatureList *flist)
}
+static void assemble_cells_from_candidates(struct image *image,
+ struct reax_search *s,
+ UnitCell *cell)
+{
+}
+
+
void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
{
struct reax_private *p;
@@ -688,22 +686,19 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex));
pmax = max_feature_resolution(image->features);
- s = search_all_axes(cell, pmax);
/* Sanity check */
if ( pmax < 1e4 ) return;
+ s = search_all_axes(cell, pmax);
find_candidates(p, image->features, pmax, fft_in, fft_out, s,
NULL, image->det);
- image->candidate_cells[0] = cell_new();
- /* FIXME: Assemble cell from candidates */
-
// refine_all_rigid_groups(image, image->candidate_cells[0], pmax,
// fft_in, fft_out, p->plan, smin, smax,
// image->det, p);
-// map_all_peaks(image);
- refine_cell(image, image->candidate_cells[0], image->features);
+
+ assemble_cells_from_candidates(image, s, cell);
fftw_free(fft_in);
fftw_free(fft_out);