aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-08-03 15:14:04 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:35 +0100
commit9b20c8133d46567d4ac90c9ecb9ec9a383285658 (patch)
treea02e1a529d1a08fd425dfd0dc6818e3f951a89ec
parent991b76fd55c9953d1e48886b6bb1d1979b1f54f9 (diff)
Do the FFT and check the results (indexing sort-of works!)
-rw-r--r--src/reax.c115
1 files changed, 105 insertions, 10 deletions
diff --git a/src/reax.c b/src/reax.c
index a50be37e..4264a676 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -45,11 +45,13 @@ struct reax_private
};
-static double check_dir(struct dvec *dir, ImageFeatureList *flist, double modv,
+static double check_dir(struct dvec *dir, ImageFeatureList *flist,
int nel, double pmax, double *fft_in,
- fftw_complex *fft_out, fftw_plan plan)
+ fftw_complex *fft_out, fftw_plan plan,
+ int smin, int smax)
{
int n, i;
+ double tot;
for ( i=0; i<nel; i++ ) {
fft_in[i] = 0.0;
@@ -74,7 +76,15 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist, double modv,
fftw_execute(plan);
- return 0.0;
+ 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);
+ }
+
+ return tot;
}
@@ -86,14 +96,21 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
- double mod_as;
+ double mod_as, mod_bs, mod_cs;
+ double als, bes, gas;
double dx, dy, dz;
int nel, n;
double pmax;
double *fft_in;
fftw_complex *fft_out;
fftw_plan plan;
+ int smin, smax;
+ double astmin, astmax;
+ double bstmin, bstmax;
+ double cstmin, cstmax;
const double cellmax = 50.0e-9; /* 50 nm max cell size */
+ const double ltol = 5.0; /* Cell "wiggle" in percent */
+ const double angtol = deg2rad(1.5); /* Angle tolerance in degrees */
assert(pp->indm == INDEXING_REAX);
p = (struct reax_private *)pp;
@@ -102,6 +119,20 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
mod_as = modulus(asx, asy, asz);
+ astmin = mod_as * (1.0-ltol/100.0);
+ astmax = mod_as * (1.0+ltol/100.0);
+
+ mod_bs = modulus(bsx, bsy, bsz);
+ bstmin = mod_bs * (1.0-ltol/100.0);
+ bstmax = mod_bs * (1.0+ltol/100.0);
+
+ mod_cs = modulus(csx, csy, csz);
+ cstmin = mod_cs * (1.0-ltol/100.0);
+ cstmax = mod_cs * (1.0+ltol/100.0);
+
+ als = angle_between(bsx, bsy, bsz, csx, csy, csz);
+ bes = angle_between(asx, asy, asz, csx, csy, csz);
+ gas = angle_between(asx, asy, asz, bsx, bsy, bsz);
n = image_feature_count(image->features);
pmax = 0.0;
@@ -120,6 +151,13 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
nel = 2.0*pmax*5.0*cellmax;
+ smin = 2.0*pmax / astmax;
+ smax = 2.0*pmax / astmin;
+
+ /* Take the smallest power of 2 greater than "nel"
+ * to make the FFT go fast */
+ nel = pow(2.0, floor((log(nel)/log(2.0))) + 1.0);
+
fft_in = fftw_malloc(nel*sizeof(double));
fft_out = fftw_malloc((nel/2 + 1)*sizeof(fftw_complex));
@@ -131,23 +169,80 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
double new_fom;
- new_fom = check_dir(&p->directions[i], image->features, mod_as,
- nel, pmax, fft_in, fft_out, plan);
+ new_fom = check_dir(&p->directions[i], image->features,
+ nel, pmax, fft_in, fft_out, plan,
+ smin, smax);
+ if ( new_fom > fom ) {
+ fom = new_fom;
+ dx = p->directions[i].x;
+ dy = p->directions[i].y;
+ dz = p->directions[i].z;
+ }
+
+ }
+ asx = dx*mod_as; asy = dy*mod_as; asz = dz*mod_as;
+
+ /* Search for b* */
+ smin = 2.0*pmax / bstmax;
+ smax = 2.0*pmax / bstmin;
+ fom = 0.0; dx = 0.0; dy = 0.0; dz = 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, asx, asy, asz);
+ if ( fabs(ang-gas) > angtol ) continue;
+
+ new_fom = check_dir(&p->directions[i], image->features,
+ nel, pmax, fft_in, fft_out, plan,
+ smin, smax);
+ if ( new_fom > fom ) {
+ fom = new_fom;
+ dx = p->directions[i].x;
+ dy = p->directions[i].y;
+ dz = p->directions[i].z;
+ }
+
+ }
+ bsx = dx*mod_bs; bsy = dy*mod_bs; bsz = dz*mod_bs;
+
+ /* Search for c* */
+ smin = 2.0*pmax / cstmax;
+ smax = 2.0*pmax / cstmin;
+ fom = 0.0; dx = 0.0; dy = 0.0; dz = 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, asx, asy, asz);
+ if ( fabs(ang-bes) > angtol ) continue;
+
+ ang = angle_between(p->directions[i].x, p->directions[i].y,
+ p->directions[i].z, bsx, bsy, bsz);
+ if ( fabs(ang-als) > angtol ) continue;
+
+ new_fom = check_dir(&p->directions[i], image->features,
+ nel, pmax, fft_in, fft_out, plan,
+ smin, smax);
if ( new_fom > fom ) {
fom = new_fom;
dx = p->directions[i].x;
- dy = p->directions[i].x;
- dz = p->directions[i].x;
+ dy = p->directions[i].y;
+ dz = p->directions[i].z;
}
}
+ csx = dx*mod_cs; csy = dy*mod_cs; csz = dz*mod_cs;
fftw_destroy_plan(plan);
fftw_free(fft_in);
fftw_free(fft_out);
- /* No improvement from zero? */
- if ( fom == 0.0 ) return;
+ image->indexed_cell = cell_new();
+ cell_set_reciprocal(image->indexed_cell, asx, asy, asz,
+ bsx, bsy, bsz, csx, csy, csz);
}