diff options
| -rw-r--r-- | config.h.in | 3 | ||||
| -rwxr-xr-x | configure | 70 | ||||
| -rw-r--r-- | configure.ac | 15 | ||||
| -rw-r--r-- | src/reax.c | 122 |
4 files changed, 200 insertions, 10 deletions
diff --git a/config.h.in b/config.h.in index ea2b3689..b09f25d8 100644 --- a/config.h.in +++ b/config.h.in @@ -38,6 +38,9 @@ /* Define to 1 if you have the <fcntl.h> header file. */ #undef HAVE_FCNTL_H +/* Define to 1 if FFTW is available. */ +#undef HAVE_FFTW + /* Define to 1 if you have the `forkpty' function. */ #undef HAVE_FORKPTY @@ -616,6 +616,8 @@ HTML_DIR GTKDOC_MKPDF GTKDOC_REBASE GTKDOC_CHECK +HAVE_FFTW_FALSE +HAVE_FFTW_TRUE HAVE_CAIRO_FALSE HAVE_CAIRO_TRUE BUILD_CUBEIT_FALSE @@ -7537,6 +7539,60 @@ fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for fftw_plan_dft_r2c_1d in -lfftw3" >&5 +$as_echo_n "checking for fftw_plan_dft_r2c_1d in -lfftw3... " >&6; } +if test "${ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d+set}" = set; then : + $as_echo_n "(cached) " >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lfftw3 $LIBS" +cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + +/* Override any GCC internal prototype to avoid an error. + Use char because int might match the return type of a GCC + builtin and then its argument prototype would still apply. */ +#ifdef __cplusplus +extern "C" +#endif +char fftw_plan_dft_r2c_1d (); +int +main () +{ +return fftw_plan_dft_r2c_1d (); + ; + return 0; +} +_ACEOF +if ac_fn_c_try_link "$LINENO"; then : + ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d=yes +else + ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d=no +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" >&5 +$as_echo "$ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" >&6; } +if test "x$ac_cv_lib_fftw3_fftw_plan_dft_r2c_1d" = x""yes; then : + + +$as_echo "#define HAVE_FFTW 1" >>confdefs.h + + FFTW_LIBS="-lfftw3" + have_fftw=true + +else + + { $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: ReAx indexing will not be available." >&5 +$as_echo "$as_me: WARNING: ReAx indexing will not be available." >&2;} + have_fftw=false + +fi + + + if test x$have_opencl = xtrue; then HAVE_OPENCL_TRUE= HAVE_OPENCL_FALSE='#' @@ -7576,6 +7632,14 @@ else fi + if test x$have_fftw = xtrue; then + HAVE_FFTW_TRUE= + HAVE_FFTW_FALSE='#' +else + HAVE_FFTW_TRUE='#' + HAVE_FFTW_FALSE= +fi + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for C compiler flag to ignore unused libraries" >&5 @@ -7642,7 +7706,7 @@ CFLAGS="$CFLAGS $GDK_pixbuf_CFLAGS $GDK_pixbuf_2_CFLAGS" LIBS="$LIBS $HDF5_LIBS -lm -lz $GSL_LIBS $GTK_LIBS $OPENCL_LIBS -pthread" LIBS="$LIBS $LIBTIFF_LIBS $libPNG_LIBS $Cairo_LIBS $GDK_pixbuf_LIBS" -LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $LDFLAGS" +LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $FFTW_LIBS $LDFLAGS" @@ -8043,6 +8107,10 @@ if test -z "${HAVE_CAIRO_TRUE}" && test -z "${HAVE_CAIRO_FALSE}"; then as_fn_error $? "conditional \"HAVE_CAIRO\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 fi +if test -z "${HAVE_FFTW_TRUE}" && test -z "${HAVE_FFTW_FALSE}"; then + as_fn_error $? "conditional \"HAVE_FFTW\" was never defined. +Usually this means the macro was only invoked conditionally." "$LINENO" 5 +fi if test -z "${ENABLE_GTK_DOC_TRUE}" && test -z "${ENABLE_GTK_DOC_FALSE}"; then as_fn_error $? "conditional \"ENABLE_GTK_DOC\" was never defined. Usually this means the macro was only invoked conditionally." "$LINENO" 5 diff --git a/configure.ac b/configure.ac index 1b7ccb70..01a89553 100644 --- a/configure.ac +++ b/configure.ac @@ -230,6 +230,18 @@ AC_CHECK_LIB([rt], [clock_gettime], ]) +AC_CHECK_LIB([fftw3], [fftw_plan_dft_r2c_1d], +[ + AC_DEFINE([HAVE_FFTW], [1], + [Define to 1 if FFTW is available.]) + FFTW_LIBS="-lfftw3" + have_fftw=true +], [ + AC_MSG_WARN([ReAx indexing will not be available.]) + have_fftw=false +]) + + dnl Conditionals... AM_CONDITIONAL([HAVE_OPENCL], test x$have_opencl = xtrue) @@ -242,6 +254,7 @@ AM_CONDITIONAL([BUILD_CUBEIT], test x$have_cairo = xtrue \ AM_CONDITIONAL([HAVE_CAIRO], test x$have_cairo = xtrue) +AM_CONDITIONAL([HAVE_FFTW], test x$have_fftw = xtrue) gl_IGNORE_UNUSED_LIBRARIES @@ -252,7 +265,7 @@ CFLAGS="$CFLAGS $GDK_pixbuf_CFLAGS $GDK_pixbuf_2_CFLAGS" LIBS="$LIBS $HDF5_LIBS -lm -lz $GSL_LIBS $GTK_LIBS $OPENCL_LIBS -pthread" LIBS="$LIBS $LIBTIFF_LIBS $libPNG_LIBS $Cairo_LIBS $GDK_pixbuf_LIBS" -LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $LDFLAGS" +LIBS="$LIBS $GDK_pixbuf_2_LIBS $TIMER_LIBS $FFTW_LIBS $LDFLAGS" GTK_DOC_CHECK([1.11],[--flavour no-tmpl]) @@ -18,6 +18,8 @@ #include <stdlib.h> #include <stdio.h> #include <math.h> +#include <assert.h> +#include <fftw3.h> #include "image.h" #include "utils.h" @@ -39,9 +41,116 @@ struct reax_private { IndexingPrivate base; struct dvec *directions; + int n_dir; }; +static double check_dir(struct dvec *dir, ImageFeatureList *flist, double modv, + double pmax, double *fft_in, fftw_complex *fft_out, + fftw_plan plan) +{ + int n, i, v; + double *vals; + + n = image_feature_count(flist); + + vals = malloc(n*sizeof(double)); + v = 0; + for ( i=0; i<n; i++ ) { + + struct imagefeature *f; + + f = image_get_feature(flist, i); + if ( f == NULL ) continue; + + vals[v] = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; + + v++; /* Yuk, yuk, yuk */ + + } + + + + + free(vals); + + return 0.0; +} + + +void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +{ + int i; + struct reax_private *p; + double fom; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double mod_as; + double dx, dy, dz; + int nel, n; + double pmax; + double *fft_in; + fftw_complex *fft_out; + fftw_plan plan; + const double cellmax = 50.0e9; /* 50 nm max cell size */ + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + mod_as = modulus(asx, asy, asz); + + n = image_feature_count(image->features); + pmax = 0.0; + for ( i=0; i<n; i++ ) { + + struct imagefeature *f; + double val; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + val = modulus(f->rx, f->ry, f->rz); + + if ( val > pmax ) pmax = val; + + } + nel = 2.0*pmax*5.0*cellmax; + + fft_in = fftw_malloc(nel*sizeof(double)); + fft_out = fftw_malloc((nel/2 + 1)*sizeof(fftw_complex)); + + plan = fftw_plan_dft_r2c_1d(nel, fft_in, fft_out, FFTW_ESTIMATE); + + /* Search for a* */ + fom = 0.0; dx = 0.0; dy = 0.0; dz = 0.0; + for ( i=0; i<p->n_dir; i++ ) { + + double new_fom; + + new_fom = check_dir(&p->directions[i], image->features, mod_as, + pmax, fft_in, fft_out, plan); + if ( new_fom > fom ) { + fom = new_fom; + dx = p->directions[i].x; + dy = p->directions[i].x; + dz = p->directions[i].x; + } + + } + + fftw_destroy_plan(plan); + fftw_free(fft_in); + fftw_free(fft_out); + + /* No improvement from zero? */ + if ( fom == 0.0 ) return; +} + + IndexingPrivate *reax_prepare() { struct reax_private *priv; @@ -58,12 +167,14 @@ IndexingPrivate *reax_prepare() angular_inc = 0.03; /* From Steller (1997) */ samp = (2.0 * M_PI) / angular_inc; - priv->directions = malloc(samp*samp*sizeof(struct dvec)); + priv->n_dir = samp*samp; + priv->directions = malloc(priv->n_dir*sizeof(struct dvec)); if ( priv == NULL) { free(priv); return NULL; } + /* Generate vectors for 1D Fourier transforms */ for ( ui=0; ui<samp; ui++ ) { for ( vi=0; vi<samp; vi++ ) { @@ -74,8 +185,9 @@ IndexingPrivate *reax_prepare() u = (double)ui/samp; v = (double)vi/samp; + /* Uniform sampling of a hemisphere */ th = 2.0 * M_PI * u; - ph = acos(2.0*v - 1.0); + ph = acos(v); dir = &priv->directions[ui + vi*samp]; @@ -88,9 +200,3 @@ IndexingPrivate *reax_prepare() return (IndexingPrivate *)priv; } - - -void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell) -{ - -} |