Make DirAx use optional

8 files changed, 145 insertions, 63 deletions

diff --git a/Makefile.am b/Makefile.am
index 0406c2f4..c4954775 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -2,5 +2,5 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h src/relrod.h \
              src/utils.h src/diffraction.h src/detector.h src/ewald.h \
              src/sfac.h src/intensities.h src/reflections.h src/list_tmp.h \
              src/statistics.h src/displaywindow.h src/render.h src/hdfsee.h \
-             data/displaywindow.ui src/dirax.h src/peaks.h
+             data/displaywindow.ui src/dirax.h src/peaks.h src/index.h
 SUBDIRS = src data
diff --git a/src/Makefile.am b/src/Makefile.am
index 948d2fbf..8dd6af85 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -17,7 +17,7 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
 process_hkl_LDADD = @LIBS@
 
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c image.c \
-                      intensities.c ewald.c peaks.c
+                      intensities.c ewald.c peaks.c index.c
 indexamajig_LDADD = @LIBS@
 
 if HAVE_GTK
diff --git a/src/dirax.c b/src/dirax.c
index 2db591f4..991a3813 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -323,63 +323,13 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 }
 
 
-static int map_position(struct image *image, double x, double y,
-                        double *rx, double *ry, double *rz)
-{
-	/* "Input" space */
-	double d;
-
-	/* Angular description of reflection */
-	double theta, psi, k;
-
-	x -= image->x_centre;
-	y -= image->y_centre;
-	k = 1.0 / image->lambda;
-
-	/* FIXME: Don't process lower CCD for now */
-	if ( y < 0 ) return 0;
-
-	if ( image->fmode == FORMULATION_CLEN ) {
-
-		/* Convert pixels to metres */
-		x /= image->resolution;
-		y /= image->resolution;	/* Convert pixels to metres */
-		d = sqrt((x*x) + (y*y));
-		theta = atan2(d, image->camera_len);
-
-	} else if (image->fmode == FORMULATION_PIXELSIZE ) {
-
-		/* Convert pixels to metres^-1 */
-		x *= image->pixel_size;
-		y *= image->pixel_size;	/* Convert pixels to metres^-1 */
-		d = sqrt((x*x) + (y*y));
-		theta = atan2(d, k);
-
-	} else {
-		ERROR("Unrecognised formulation mode in mapping_scale.\n");
-		return -1;
-	}
-
-	psi = atan2(y, x);
-
-	*rx = k*sin(theta)*cos(psi);
-	*ry = k*sin(theta)*sin(psi);
-	*rz = k - k*cos(theta);
-
-	return 0;
-}
-
-
-void index_pattern(struct image *image, int no_index, int dump_peaks)
+void run_dirax(struct image *image, int no_index)
 {
 	unsigned int opts;
 	int saved_stderr;
 	FILE *fh;
 	int i;
 
-	/* Do peak search and splurge out 'xfel.drx' */
-	search_peaks(image, dump_peaks);
-
 	fh = fopen("xfel.drx", "w");
 	if ( !fh ) {
 		ERROR("Couldn't open temporary file xfel.drx\n");
@@ -390,14 +340,11 @@ void index_pattern(struct image *image, int no_index, int dump_peaks)
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
 		struct imagefeature *f;
-		double rx = 0.0;
-		double ry = 0.0;
-		double rz = 0.0;
 
 		f = image_get_feature(image->features, i);
-		map_position(image, f->x, f->y, &rx, &ry, &rz);
+
 		fprintf(fh, "%10f %10f %10f %8f\n",
-		        rx/1e10, ry/1e10, rz/1e10, 1.0);
+		        f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0);
 
 	}
 	fclose(fh);
diff --git a/src/dirax.h b/src/dirax.h
index 4ebad677..58a0416e 100644
--- a/src/dirax.h
+++ b/src/dirax.h
@@ -18,7 +18,7 @@
 #endif
 
 
-extern void index_pattern(struct image *image, int no_index, int dump_peaks);
+extern void run_dirax(struct image *image, int no_index);
 
 
 #endif	/* DIRAX_H */
diff --git a/src/image.h b/src/image.h
index 9649fe8e..7cf47761 100644
--- a/src/image.h
+++ b/src/image.h
@@ -46,6 +46,11 @@ struct imagefeature {
 	/* Distance between this feature and its partner, if any. */
 	double				partner_d;
 
+	/* Reciprocal space coordinates (m^-1 of course) of this feature */
+	double				rx;
+	double				ry;
+	double				rz;
+
 };
 
 /* An opaque type representing a list of image features */
diff --git a/src/index.c b/src/index.c
new file mode 100644
index 00000000..c8d5dbe6
--- /dev/null
+++ b/src/index.c
@@ -0,0 +1,98 @@
+/*
+ * index.c
+ *
+ * Perform indexing (somehow)
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "dirax.h"
+
+
+static int map_position(struct image *image, double x, double y,
+                        double *rx, double *ry, double *rz)
+{
+	/* "Input" space */
+	double d;
+
+	/* Angular description of reflection */
+	double theta, psi, k;
+
+	x -= image->x_centre;
+	y -= image->y_centre;
+	k = 1.0 / image->lambda;
+
+	/* FIXME: Don't process lower CCD for now */
+	if ( y < 0 ) return 0;
+
+	if ( image->fmode == FORMULATION_CLEN ) {
+
+		/* Convert pixels to metres */
+		x /= image->resolution;
+		y /= image->resolution;	/* Convert pixels to metres */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, image->camera_len);
+
+	} else if (image->fmode == FORMULATION_PIXELSIZE ) {
+
+		/* Convert pixels to metres^-1 */
+		x *= image->pixel_size;
+		y *= image->pixel_size;	/* Convert pixels to metres^-1 */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, k);
+
+	} else {
+		ERROR("Unrecognised formulation mode in mapping_scale.\n");
+		return -1;
+	}
+
+	psi = atan2(y, x);
+
+	*rx = k*sin(theta)*cos(psi);
+	*ry = k*sin(theta)*sin(psi);
+	*rz = k - k*cos(theta);
+
+	return 0;
+}
+
+
+void index_pattern(struct image *image, int no_index, int dump_peaks,
+                   int use_dirax)
+{
+	int i;
+
+	/* Perform 'fine' peak search */
+	search_peaks(image, dump_peaks);
+
+	/* Map positions to 3D. FIXME: Handle lower detector */
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+		struct imagefeature *f;
+
+		f = image_get_feature(image->features, i);
+		map_position(image, f->x, f->y, &f->rx, &f->ry, &f->rz);
+
+	}
+
+	if ( use_dirax ) {
+		run_dirax(image, no_index);
+		return;
+	}
+}
diff --git a/src/index.h b/src/index.h
new file mode 100644
index 00000000..ea0b4485
--- /dev/null
+++ b/src/index.h
@@ -0,0 +1,25 @@
+/*
+ * index.h
+ *
+ * Perform indexing (somehow)
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef INDEX_H
+#define INDEX_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+extern void index_pattern(struct image *image, int no_index, int dump_peaks,
+                          int use_dirax);
+
+
+#endif	/* INDEX_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 7c2e1b5d..7f863a65 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -23,7 +23,7 @@
 
 #include "utils.h"
 #include "hdf5-file.h"
-#include "dirax.h"
+#include "index.h"
 #include "intensities.h"
 #include "ewald.h"
 #include "peaks.h"
@@ -39,8 +39,10 @@ static void show_help(const char *s)
 "\n"
 "  -i, --input=<filename>  Specify file containing list of images to process.\n"
 "                           '-' means stdin, which is the default.\n"
-"      --no-index          Do everything else (including fine peak search),\n"
-"                           but don't invoke the indexing program.\n"
+"      --no-index          Do everything else (including fine peak search and\n"
+"                           writing 'xfel.drx' if DirAx is being used), but\n"
+"                           don't actually index.\n"
+"      --dirax             Use DirAx for indexing.\n"
 "\n");
 }
 
@@ -55,6 +57,7 @@ int main(int argc, char *argv[])
 	int n_hits;
 	int config_noindex = 0;
 	int config_dumpfound = 0;
+	int config_dirax = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -62,6 +65,7 @@ int main(int argc, char *argv[])
 		{"input",              1, NULL,               'i'},
 		{"no-index",           0, &config_noindex,     1},
 		{"dump-found-peaks",   0, &config_dumpfound,   1},
+		{"dirax",              0, &config_dirax,       1},
 		{0, 0, NULL, 0}
 	};
 
@@ -135,7 +139,10 @@ int main(int argc, char *argv[])
 		if ( fom > 0 ) {
 
 			/* Calculate orientation matrix (by magic) */
-			index_pattern(&image, config_noindex, config_dumpfound);
+			index_pattern(&image, config_noindex, config_dumpfound,
+			              config_dirax);
+
+			if ( image.molecule == NULL ) continue;
 
 			/* View head-on (unit cell is tilted) */
 			image.orientation.x = 0.0;