aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2011-09-06 15:39:02 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:42 +0100
commit636b71d41e8b2cd318e6df24159a6fb65301e5df (patch)
tree925e82dc8c898ff2257b49ba8ac35826d22d314e
parent1e62bb9dd2893e61c953528fe77adab0eafb4f84 (diff)
Initial candidate ReAx stuff
-rw-r--r--libcrystfel/src/reax.c478
1 files changed, 278 insertions, 200 deletions
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index ff1ea582..49231438 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -45,6 +45,31 @@ struct dvec
};
+struct reax_candidate
+{
+ struct dvec v; /* This is the vector for the candidate */
+ double strength; /* How much intensity fell into the range */
+};
+
+
+struct reax_search_v
+{
+ unsigned int smin;
+ unsigned int smax; /* Search for vector in this range */
+
+ struct reax_candidate *cand; /* Candidate vectors go here */
+ int n_cand; /* There are this many candidates */
+};
+
+
+struct reax_search
+{
+ struct reax_search_v *search; /* Search for these vectors */
+ int n_search; /* There are this many vectors to find */
+ double pmax; /* The maximum feature resolution */
+};
+
+
struct reax_private
{
IndexingPrivate base;
@@ -65,14 +90,12 @@ struct reax_private
};
-static double check_dir(struct dvec *dir, ImageFeatureList *flist,
+static void fill_and_transform(struct dvec *dir, ImageFeatureList *flist,
int nel, double pmax, double *fft_in,
fftw_complex *fft_out, fftw_plan plan,
- int smin, int smax,
const char *rg, struct detector *det)
{
int n, i;
- double tot;
for ( i=0; i<nel; i++ ) {
fft_in[i] = 0.0;
@@ -107,19 +130,225 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
}
fftw_execute_dft_r2c(plan, fft_in, fft_out);
+}
+
+
+static double check_dir(struct dvec *dir, ImageFeatureList *flist,
+ int nel, double pmax, double *fft_in,
+ fftw_complex *fft_out, fftw_plan plan,
+ struct reax_search *s,
+ const char *rg, struct detector *det)
+{
+ int i;
+ double tot;
+
+ fill_and_transform(dir, flist, nel, pmax, fft_in, fft_out,
+ plan, rg, det);
tot = 0.0;
- for ( i=smin; i<=smax; i++ ) {
- double re, im;
- re = fft_out[i][0];
- im = fft_out[i][1];
- tot += sqrt(re*re + im*im);
+ for ( i=0; i<s->n_search; i++ ) {
+
+ double tot = 0.0;
+ double peak = 0.0;
+ double peak_mod = 0.0;
+ int j;
+
+ for ( j=s->search[i].smin; j<=s->search[i].smax; j++ ) {
+
+ double re, im;
+
+ re = fft_out[j][0];
+ im = fft_out[j][1];
+ peak += sqrt(re*re + im*im);
+
+ }
+
+ for ( j=0; j<nel; j++ ) {
+
+ double re, im, am;
+
+ re = fft_out[j][0];
+ 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);
+ }
+ }
+
+ }
+
+ /* If sufficiently strong, add to list of candidates */
+ if ( peak > 2.0*(tot/nel) ) {
+
+ size_t ns;
+ struct reax_candidate *cnew;
+ int cpos;
+
+ cpos = s->search[i].n_cand;
+
+ ns = (cpos+1) * sizeof(struct reax_candidate);
+ cnew = realloc(s->search[i].cand, ns);
+ if ( cnew == NULL ) {
+ ERROR("Failed to add candidate.\n");
+ } else {
+
+ 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].n_cand++;
+
+ }
+
+ }
+
}
return tot;
}
+static void fine_search(struct reax_private *p, ImageFeatureList *flist,
+ double pmax, double *fft_in, fftw_complex *fft_out,
+ struct reax_search_v *sv, struct reax_candidate *c,
+ const char *rg, struct detector *det)
+{
+ double fom = 0.0;
+ double th, ph;
+ double inc;
+ struct dvec dir;
+ int i, s;
+ double max;
+
+ inc = p->angular_inc / 100.0;
+
+ for ( th=c->v.th-p->angular_inc; th<=c->v.th+p->angular_inc; th+=inc ) {
+ for ( ph=c->v.ph-p->angular_inc; ph<=c->v.ph+p->angular_inc; ph+=inc ) {
+
+ double new_fom;
+
+ dir.x = cos(ph) * sin(th);
+ dir.y = sin(ph) * sin(th);
+ dir.z = cos(th);
+ dir.th = th;
+ dir.ph = ph;
+
+ fill_and_transform(&dir, flist, p->nel, pmax, fft_in, fft_out,
+ p->plan, rg, det);
+
+ for ( i=sv->smin; i<=sv->smax; i++ ) {
+
+ double re, im, m;
+
+ re = fft_out[i][0];
+ im = fft_out[i][1];
+ m = sqrt(re*re + im*im);
+ if ( m > max ) {
+ max = m;
+ s = i;
+ }
+ }
+
+ if ( new_fom > fom ) {
+ fom = new_fom;
+ c->v = dir;
+ }
+
+ }
+ }
+}
+
+
+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)
+{
+ double fom = 0.0;
+ int i;
+ double th, ph;
+
+ fom = 0.0; th = 0.0; ph = 0.0;
+ for ( i=0; i<p->n_dir; i++ ) {
+
+ double new_fom;
+
+ new_fom = check_dir(&p->directions[i], flist,
+ p->nel, pmax, fft_in, fft_out, p->plan,
+ s, NULL, NULL);
+ if ( new_fom > fom ) {
+ fom = new_fom;
+ th = p->directions[i].th;
+ ph = p->directions[i].ph;
+ }
+
+ }
+
+ for ( i=0; i<s->n_search; i++ ) {
+
+ struct reax_search_v *sv;
+ int j;
+
+ sv = &s->search[i];
+ 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;
+}
+
+
+/* Set up search parameters to look for all three cell axes */
+static struct reax_search *search_all_axes(UnitCell *cell, double pmax)
+{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ double mod_a, mod_b, mod_c;
+ double amin, amax;
+ double bmin, bmax;
+ double cmin, cmax;
+ unsigned int smin, smax;
+ const double ltol = 5.0; /* Direct space axis length tolerance in % */
+ struct reax_search *s;
+
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ mod_a = modulus(ax, ay, az);
+ amin = mod_a * (1.0-ltol/100.0);
+ amax = mod_a * (1.0+ltol/100.0);
+
+ mod_b = modulus(bx, by, bz);
+ bmin = mod_b * (1.0-ltol/100.0);
+ bmax = mod_b * (1.0+ltol/100.0);
+
+ mod_c = modulus(cx, cy, cz);
+ cmin = mod_c * (1.0-ltol/100.0);
+ cmax = mod_c * (1.0+ltol/100.0);
+
+ s = malloc(3*sizeof(*s));
+ s->pmax = pmax;
+ s->n_search = 1;
+ 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;
+ s->search[1].smin = smin; s->search[1].smax = smax;
+ smin = 2.0*pmax * cmin; smax = 2.0*pmax * cmax;
+ s->search[2].smin = smin; s->search[2].smax = smax;
+
+ return s;
+}
+
+
/* Refine a direct space vector. From Clegg (1984) */
static double iterate_refine_vector(double *x, double *y, double *z,
ImageFeatureList *flist)
@@ -239,63 +468,6 @@ static void refine_cell(struct image *image, UnitCell *cell,
}
-static void fine_search(struct reax_private *p, ImageFeatureList *flist,
- int nel, double pmax, double *fft_in,
- fftw_complex *fft_out, fftw_plan plan,
- int smin, int smax, double th_cen, double ph_cen,
- double *x, double *y, double *z)
-{
- double fom = 0.0;
- double th, ph;
- double inc;
- struct dvec dir;
- int i, s;
- double max, modv;
-
- inc = p->angular_inc / 100.0;
-
- for ( th=th_cen-p->angular_inc; th<=th_cen+p->angular_inc; th+=inc ) {
- for ( ph=ph_cen-p->angular_inc; ph<=ph_cen+p->angular_inc; ph+=inc ) {
-
- double new_fom;
-
- dir.x = cos(ph) * sin(th);
- dir.y = sin(ph) * sin(th);
- dir.z = cos(th);
-
- new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out,
- plan, smin, smax, NULL, NULL);
-
- if ( new_fom > fom ) {
- fom = new_fom;
- *x = dir.x; *y = dir.y; *z = dir.z;
- }
-
- }
- }
-
- s = -1;
- max = 0.0;
- for ( i=smin; i<=smax; i++ ) {
-
- double re, im, m;
-
- re = fft_out[i][0];
- im = fft_out[i][1];
- m = sqrt(re*re + im*im);
- if ( m > max ) {
- max = m;
- s = i;
- }
-
- }
- assert(s>0);
-
- modv = (double)s / (2.0*pmax);
- *x *= modv; *y *= modv; *z *= modv;
-}
-
-
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,
@@ -309,8 +481,7 @@ static double get_model_phase(double x, double y, double z, ImageFeatureList *f,
dir.x = x; dir.y = y; dir.z = z;
- check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax,
- rg, det);
+ fill_and_transform(&dir, f, nel, pmax, fft_in, fft_out, plan, rg, det);
s = -1;
max = 0.0;
@@ -336,7 +507,7 @@ static double get_model_phase(double x, double y, double z, ImageFeatureList *f,
static void refine_rigid_group(struct image *image, UnitCell *cell,
- const char *rg, int nel, double pmax,
+ const char *rg, double pmax,
double *fft_in, fftw_complex *fft_out,
fftw_plan plan, int smin, int smax,
struct detector *det, struct reax_private *pr)
@@ -353,8 +524,6 @@ static void refine_rigid_group(struct image *image, UnitCell *cell,
signed int aix, aiy;
signed int bix, biy;
signed int cix, ciy;
- double max;
- int max_i, max_j;
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -362,12 +531,15 @@ static void refine_rigid_group(struct image *image, UnitCell *cell,
mb = modulus(bx, by, bz);
mc = modulus(cx, cy, cz);
- pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, nel, pmax,
- fft_in, fft_out, plan, smin, smax, rg, det);
- phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, nel, pmax,
- fft_in, fft_out, plan, smin, smax, rg, det);
- phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, nel, pmax,
- fft_in, fft_out, plan, smin, smax, rg, det);
+ pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features,
+ pr->nel, pmax, fft_in, fft_out, plan,
+ smin, smax, rg, det);
+ phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features,
+ pr->nel, pmax, fft_in, fft_out, plan,
+ smin, smax, rg, det);
+ phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features,
+ pr->nel, pmax, fft_in, fft_out, plan,
+ smin, smax, rg, det);
for ( i=0; i<det->n_panels; i++ ) {
if ( det->panels[i].rigid_group == rg ) {
@@ -461,7 +633,7 @@ static void refine_rigid_group(struct image *image, UnitCell *cell,
static void refine_all_rigid_groups(struct image *image, UnitCell *cell,
- int nel, double pmax,
+ double pmax,
double *fft_in, fftw_complex *fft_out,
fftw_plan plan, int smin, int smax,
struct detector *det,
@@ -471,67 +643,25 @@ static void refine_all_rigid_groups(struct image *image, UnitCell *cell,
for ( i=0; i<image->det->num_rigid_groups; i++ ) {
refine_rigid_group(image, cell, image->det->rigid_groups[i],
- nel, pmax, fft_in, fft_out, plan, smin, smax,
+ pmax, fft_in, fft_out, plan, smin, smax,
det, p);
}
}
-void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
+static double max_feature_resolution(ImageFeatureList *flist)
{
- int i;
- struct reax_private *p;
- double fom;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- double mod_a, mod_b, mod_c;
- double al, be, ga;
- double th, ph;
- double *fft_in;
- fftw_complex *fft_out;
- int smin, smax;
- double amin, amax;
- double bmin, bmax;
- double cmin, cmax;
double pmax;
- int n;
- const double ltol = 5.0; /* Direct space axis length
- * tolerance in percent */
- const double angtol = deg2rad(1.5); /* Direct space angle tolerance
- * in radians */
-
- assert(pp->indm == INDEXING_REAX);
- p = (struct reax_private *)pp;
-
- fft_in = fftw_malloc(p->nel*sizeof(double));
- fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex));
-
- cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- mod_a = modulus(ax, ay, az);
- amin = mod_a * (1.0-ltol/100.0);
- amax = mod_a * (1.0+ltol/100.0);
-
- mod_b = modulus(bx, by, bz);
- bmin = mod_b * (1.0-ltol/100.0);
- bmax = mod_b * (1.0+ltol/100.0);
-
- mod_c = modulus(cx, cy, cz);
- cmin = mod_c * (1.0-ltol/100.0);
- cmax = mod_c * (1.0+ltol/100.0);
-
- al = angle_between(bx, by, bz, cx, cy, cz);
- be = angle_between(ax, ay, az, cx, cy, cz);
- ga = angle_between(ax, ay, az, bx, by, bz);
+ int i, n;
pmax = 0.0;
- n = image_feature_count(image->features);
+ n = image_feature_count(flist);
for ( i=0; i<n; i++ ) {
struct imagefeature *f;
double val;
- f = image_get_feature(image->features, i);
+ f = image_get_feature(flist, i);
if ( f == NULL ) continue;
val = modulus(f->rx, f->ry, f->rz);
@@ -539,92 +669,40 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
- /* Sanity check */
- if ( pmax < 1e4 ) return;
-
- /* Search for a */
- smin = 2.0*pmax * amin;
- smax = 2.0*pmax * amax;
- fom = 0.0; th = 0.0; ph = 0.0;
- for ( i=0; i<p->n_dir; i++ ) {
-
- double new_fom;
-
- new_fom = check_dir(&p->directions[i], image->features,
- p->nel, pmax, fft_in, fft_out, p->plan,
- smin, smax, NULL, NULL);
- if ( new_fom > fom ) {
- fom = new_fom;
- th = p->directions[i].th;
- ph = p->directions[i].ph;
- }
-
- }
- fine_search(p, image->features, p->nel, pmax, fft_in, fft_out,
- p->plan, smin, smax, th, ph, &ax, &ay, &az);
-
- /* Search for b */
- smin = 2.0*pmax * bmin;
- smax = 2.0*pmax * bmax;
- fom = 0.0; th = 0.0; ph = 0.0;
- for ( i=0; i<p->n_dir; i++ ) {
-
- double new_fom, ang;
-
- ang = angle_between(p->directions[i].x, p->directions[i].y,
- p->directions[i].z, ax, ay, az);
- if ( fabs(ang-ga) > angtol ) continue;
-
- new_fom = check_dir(&p->directions[i], image->features,
- p->nel, pmax, fft_in, fft_out, p->plan,
- smin, smax, NULL, NULL);
- if ( new_fom > fom ) {
- fom = new_fom;
- th = p->directions[i].th;
- ph = p->directions[i].ph;
- }
+ return pmax;
+}
- }
- fine_search(p, image->features, p->nel, pmax, fft_in, fft_out,
- p->plan, smin, smax, th, ph, &bx, &by, &bz);
- /* Search for c */
- smin = 2.0*pmax * cmin;
- smax = 2.0*pmax * cmax;
- fom = 0.0; th = 0.0; ph = 0.0;
- for ( i=0; i<p->n_dir; i++ ) {
+void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
+{
+ struct reax_private *p;
+ double *fft_in;
+ fftw_complex *fft_out;
+ double pmax;
+ struct reax_search *s;
- double new_fom, ang;
+ assert(pp->indm == INDEXING_REAX);
+ p = (struct reax_private *)pp;
- ang = angle_between(p->directions[i].x, p->directions[i].y,
- p->directions[i].z, ax, ay, az);
- if ( fabs(ang-be) > angtol ) continue;
+ fft_in = fftw_malloc(p->nel*sizeof(double));
+ fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex));
- ang = angle_between(p->directions[i].x, p->directions[i].y,
- p->directions[i].z, bx, by, bz);
- if ( fabs(ang-al) > angtol ) continue;
+ pmax = max_feature_resolution(image->features);
+ s = search_all_axes(cell, pmax);
- new_fom = check_dir(&p->directions[i], image->features,
- p->nel, pmax, fft_in, fft_out, p->plan,
- smin, smax, NULL, NULL);
- if ( new_fom > fom ) {
- fom = new_fom;
- th = p->directions[i].th;
- ph = p->directions[i].ph;
- }
+ /* Sanity check */
+ if ( pmax < 1e4 ) return;
- }
- fine_search(p, image->features, p->nel, pmax, fft_in, fft_out,
- p->plan, smin, smax, th, ph, &cx, &cy, &cz);
+ find_candidates(p, image->features, pmax, fft_in, fft_out, s,
+ NULL, image->det);
image->candidate_cells[0] = cell_new();
- cell_set_cartesian(image->candidate_cells[0],
- ax, ay, az, bx, by, bz, cx, cy, cz);
+ /* FIXME: Assemble cell from candidates */
- refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax,
- fft_in, fft_out, p->plan, smin, smax,
- image->det, p);
- map_all_peaks(image);
+// 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);
fftw_free(fft_in);