aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/asdf.c
diff options
context:
space:
mode:
authorAlexandra Tolstikova <alexandra.tolstikova@desy.de>2015-02-11 14:00:58 +0100
committerAlexandra Tolstikova <alexandra.tolstikova@desy.de>2015-02-11 14:00:58 +0100
commit3605e92c4461d4742df22e98e9d1a4bc8e7cbe33 (patch)
tree083f1211bc240b378c9770cc7d9f91e5e9d2175b /libcrystfel/src/asdf.c
parentf8d230451e7e0f71d5bf90d5f1079b89d19c8272 (diff)
Use one fftw_plan for all fft calculations
Diffstat (limited to 'libcrystfel/src/asdf.c')
-rw-r--r--libcrystfel/src/asdf.c27
1 files changed, 17 insertions, 10 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index 6ad29645..41973c9a 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -161,7 +161,8 @@ static int cross_product(gsl_vector *vec1, gsl_vector *vec2,
}
/* Returns triple product of three gsl_vectors */
-static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) {
+static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3)
+{
double volume;
gsl_vector *cross = gsl_vector_alloc(3);
@@ -289,7 +290,8 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3,
return 1;
}
-static float find_ds_fft(double *projections, int N_projections, double d_max) {
+static float find_ds_fft(double *projections, int N_projections, double d_max,
+ fftw_plan p) {
int n = N_projections;
double projections_sorted[n];
@@ -300,7 +302,7 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) {
int N = 1024; // number of points in fft calculation
- double *in = malloc(sizeof(double) * N);
+ double in[N];
fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N);
for ( i = 0; i < N; i++ ) {
@@ -311,10 +313,12 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) {
in[(int)((projections_sorted[i] - projections_sorted[0]) /
(projections_sorted[n - 1] - projections_sorted[0]) * N)] ++;
}
-
- fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
- fftw_execute(p);
- fftw_destroy_plan(p);
+
+ if ( !p ) {
+ p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
+ }
+
+ fftw_execute_dft_r2c(p, in, out);
int i_max = (int)(d_max * (projections_sorted[n - 1] -
projections_sorted[0]));
@@ -332,7 +336,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) {
double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d;
- free(in);
fftw_free(out);
return ds;
}
@@ -855,7 +858,9 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
int N_tvectors = 0;
int n_max = 0; // maximum number of reflections fitting one of tvectors
-
+
+ fftw_plan p;
+
for ( i = 0; i < i_max; i++ ) {
if ( calc_normal(reflections[triplets[i][0]],
reflections[triplets[i][1]],
@@ -870,7 +875,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
}
/* Find ds - period in 1d lattice of projections */
- ds = find_ds_fft(projections, N_reflections, d_max);
+ ds = find_ds_fft(projections, N_reflections, d_max, p);
/* Refine ds, write 1 to fits[i] if reflections[i]
* fits ds */
@@ -933,6 +938,8 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
}
free(triplets);
+ fftw_destroy_plan(p);
+
if ( c->n ) return 1;
return 0;