Commit bca7a90b authored by Thomas White's avatar Thomas White
Browse files

More ReAx stuff

parent d7ee0d1a
......@@ -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
......
......@@ -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)
{
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment