aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-08-04 15:39:37 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:36 +0100
commite30a51f72740d95e5468b6ac5a4ea64f4b38f1d0 (patch)
tree449fac36356617b035b692994648f303184f0789 /src
parent4d10776ec53555a236c18dd27d660aa5563d8b64 (diff)
"Walk" the FFT graph to get more accurate axis lengths
Diffstat (limited to 'src')
-rw-r--r--src/reax.c78
1 files changed, 71 insertions, 7 deletions
diff --git a/src/reax.c b/src/reax.c
index d5d16a24..3feb564a 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -89,22 +89,82 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
}
+#define idx_to_m(a) ( 1.0/((2.0*pmax/(double)(a))) )
+
+static void walk_graph(double *x, double *y, double *z, int smin, int smax,
+ fftw_complex *fft_out, int nel, double pmax,
+ double modv_exp)
+{
+ int i, s, mult;
+ double max, modv;
+
+ FILE *fh = fopen("fft.dat", "w");
+ for ( i=0; i<nel/2+1; i++ ) {
+ double re, im, m;
+ re = fft_out[i][0];
+ im = fft_out[i][1];
+ m = sqrt(re*re + im*im);
+ fprintf(fh, "%i %f\n", i, m);
+ }
+ fclose(fh);
+
+ //*x *= modv_exp; *y *= modv_exp; *z *= modv_exp;
+ //return;
+
+ mult = 1;
+ do {
+
+ double new_s;
+
+ 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);
+
+ //STATUS("Estimated axis length:%.5f nm\n",
+ // (idx_to_m(s)/mult)*1e9);
+
+ new_s = (double)s*(mult+1)/mult;
+ smin = new_s - 1;
+ smax = new_s + 1;
+ mult++;
+
+ } while ( mult<5 );
+
+ modv = 2.0*pmax / (double)s;
+ modv *= mult-1;
+ *x *= modv; *y *= modv; *z *= modv;
+
+ //exit(1);
+}
+
+
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 modv,
+ int smin, int smax,
double dx, double dy, double dz,
- double *x, double *y, double *z)
+ double *x, double *y, double *z,
+ double modv_exp)
{
const int fine_samp = 100;
double fom = 0.0;
signed int ui, vi;
+ struct dvec dir;
for ( ui=-fine_samp; ui<fine_samp; ui++ ) {
for ( vi=-fine_samp; vi<fine_samp; vi++ ) {
double u, v;
- struct dvec dir;
double new_fom;
double tx, ty, tz;
@@ -126,11 +186,15 @@ static void fine_search(struct reax_private *p, ImageFeatureList *flist,
if ( new_fom > fom ) {
fom = new_fom;
- *x = dir.x*modv; *y = dir.y*modv; *z = dir.z*modv;
+ *x = dir.x; *y = dir.y; *z = dir.z;
}
}
}
+
+ dir.x = *x; dir.y = *y; dir.z = *z;
+ check_dir(&dir, flist, nel, pmax, fft_in, fft_out, plan, smin, smax);
+ walk_graph(x, y, z, smin, smax, fft_out, nel, pmax, modv_exp);
}
@@ -229,7 +293,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
fine_search(p, image->features, nel, pmax, fft_in, fft_out, plan,
- smin, smax, mod_as, dx, dy, dz, &asx, &asy, &asz);
+ smin, smax, dx, dy, dz, &asx, &asy, &asz, mod_as);
/* Search for b* */
smin = 2.0*pmax / bstmax;
@@ -255,7 +319,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
fine_search(p, image->features, nel, pmax, fft_in, fft_out, plan,
- smin, smax, mod_bs, dx, dy, dz, &bsx, &bsy, &bsz);
+ smin, smax, dx, dy, dz, &bsx, &bsy, &bsz, mod_bs);
/* Search for c* */
smin = 2.0*pmax / cstmax;
@@ -285,7 +349,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
}
fine_search(p, image->features, nel, pmax, fft_in, fft_out, plan,
- smin, smax, mod_cs, dx, dy, dz, &csx, &csy, &csz);
+ smin, smax, dx, dy, dz, &csx, &csy, &csz, mod_cs);
fftw_destroy_plan(plan);
fftw_free(fft_in);