More ReAx stuff

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
 
diff --git a/configure b/configure
index e47de008..569d2a6b 100755
--- a/configure
+++ b/configure
@@ -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])
diff --git a/src/reax.c b/src/reax.c
index 77ca124f..20e8c94c 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -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)
-{
-
-}