From 3605e92c4461d4742df22e98e9d1a4bc8e7cbe33 Mon Sep 17 00:00:00 2001 From: Alexandra Tolstikova Date: Wed, 11 Feb 2015 14:00:58 +0100 Subject: Use one fftw_plan for all fft calculations --- libcrystfel/src/asdf.c | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) (limited to 'libcrystfel') 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; -- cgit v1.2.3