aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/asdf.c72
1 files changed, 51 insertions, 21 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index 53214a12..2f3ddd04 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -17,10 +17,19 @@
#include "asdf.h"
+struct fftw_vars {
+ int N;
+ fftw_plan p;
+ double *in;
+ fftw_complex *out;
+};
+
+
struct asdf_private {
IndexingMethod indm;
float *ltl;
UnitCell *template;
+ struct fftw_vars fftw;
};
@@ -32,6 +41,30 @@ struct tvector {
};
+struct fftw_vars fftw_vars_new()
+{
+ int N = 1024;
+
+ struct fftw_vars fftw;
+
+ fftw.N = N;
+ fftw.in = fftw_malloc(N * sizeof(double));
+ fftw.out = fftw_malloc(N * sizeof(fftw_complex));
+ fftw.p = fftw_plan_dft_r2c_1d(N, fftw.in, fftw.out, FFTW_MEASURE);
+
+ return fftw;
+}
+
+
+static int fftw_vars_free(struct fftw_vars fftw)
+{
+ fftw_free(fftw.in);
+ fftw_free(fftw.out);
+ fftw_destroy_plan(fftw.p);
+ fftw_cleanup();
+}
+
+
struct asdf_cell {
gsl_vector *axes[3];
gsl_vector *reciprocal[3];
@@ -311,7 +344,7 @@ static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3,
static float find_ds_fft(double *projections, int N_projections, double d_max,
- fftw_plan p)
+ struct fftw_vars fftw)
{
int n = N_projections;
double projections_sorted[n];
@@ -320,10 +353,10 @@ static float find_ds_fft(double *projections, int N_projections, double d_max,
int i;
- int N = 1024; // number of points in fft calculation
-
- double in[N];
- fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N);
+ int N = fftw.N; // number of points in fft calculation
+ double *in = fftw.in;
+ fftw_complex *out = fftw.out;
+ fftw_plan p = fftw.p;
for ( i = 0; i < N; i++ ) {
in[i] = 0;
@@ -334,10 +367,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max,
(projections_sorted[n - 1] - projections_sorted[0]) * N)] ++;
}
- 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] -
@@ -356,7 +385,6 @@ static float find_ds_fft(double *projections, int N_projections, double d_max,
double ds = (projections_sorted[n - 1] - projections_sorted[0]) / d;
- fftw_free(out);
return ds;
}
@@ -865,7 +893,7 @@ static double angle_between_gsl(gsl_vector *a, gsl_vector *b)
static int index_refls(gsl_vector **reflections, int N_reflections,
double d_max, double volume_min, double volume_max,
double LevelFit, double IndexFit, int i_max,
- struct asdf_cell *c)
+ struct asdf_cell *c, struct fftw_vars fftw)
{
int i, j, k, l, n;
@@ -906,8 +934,6 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
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]],
@@ -922,7 +948,8 @@ 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, p);
+ ds = find_ds_fft(projections, N_reflections, d_max,
+ fftw);
/* Refine ds, write 1 to fits[i] if reflections[i]
* fits ds */
@@ -952,7 +979,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
if (n > n_max) n_max = n;
}
}
-
+
if ( (i != 0 && i % 10000 == 0) || i == i_max - 1 ) {
/* Sort tvectors by length */
@@ -984,9 +1011,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
free(triplets[i]);
}
free(triplets);
-
- fftw_destroy_plan(p);
-
+
if ( c->n ) return 1;
return 0;
@@ -1035,7 +1060,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
if ( N_reflections == 0 ) return 0;
j = index_refls(reflections, N_reflections, d_max, volume_min,
- volume_max, LevelFit, IndexFit, i_max, &c);
+ volume_max, LevelFit, IndexFit, i_max, &c, dp->fftw);
for ( i = 0; i < N_reflections; i++ ) {
gsl_vector_free(reflections[i]);
@@ -1055,7 +1080,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
gsl_vector_get(c.axes[2], 1) * 1e-10,
gsl_vector_get(c.axes[2], 2) * 1e-10);
- if ( check_cell(dp, image, uc) ) {
+ if ( check_cell(dp, image, uc) ) {
+ asdf_cell_free(c);
cell_free(uc);
return 1;
}
@@ -1063,6 +1089,7 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
cell_free(uc);
}
+ asdf_cell_free(c);
return 0;
}
@@ -1093,7 +1120,9 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
dp->ltl = ltl;
dp->template = cell;
dp->indm = *indm;
-
+
+ dp->fftw = fftw_vars_new();
+
return (IndexingPrivate *)dp;
}
@@ -1102,5 +1131,6 @@ void asdf_cleanup(IndexingPrivate *pp)
{
struct asdf_private *p;
p = (struct asdf_private *)pp;
+ fftw_vars_free(p->fftw);
free(p);
}