Complete indexing functionality

6 files changed, 60 insertions, 27 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 978ab602..9010373c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -16,7 +16,8 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
                       reflections.c
 process_hkl_LDADD = @LIBS@
 
-indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c image.c
+indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c image.c \
+                      intensities.c ewald.c
 indexamajig_LDADD = @LIBS@
 
 if HAVE_GTK
diff --git a/src/dirax.c b/src/dirax.c
index 0c509e06..8606ff44 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -33,6 +33,7 @@
 #include "image.h"
 #include "dirax.h"
 #include "utils.h"
+#include "sfac.h"
 
 
 typedef enum {
@@ -67,10 +68,13 @@ static void dirax_parseline(const char *line, struct image *image)
 		if ( line[i] == 'R' ) rf = 1;
 		if ( (line[i] == 'D') && rf ) {
 			image->dirax_read_cell = 1;
-			if ( image->cell ) {
-				free(image->cell);
+			if ( image->molecule == NULL ) {
+				image->molecule = malloc(sizeof(struct molecule));
+			} else if ( image->molecule->cell ) {
+				free(image->molecule->cell);
+				free(image->molecule);
 			}
-			image->cell = cell_new();
+			image->molecule->cell = cell_new();
 			return;
 		}
 		i++;
@@ -79,24 +83,27 @@ static void dirax_parseline(const char *line, struct image *image)
 	/* Parse unit cell vectors as appropriate */
 	if ( image->dirax_read_cell == 1 ) {
 		/* First row of unit cell values */
-		float x1, x2, x3;
-		sscanf(line, "%f %f %f", &x1, &x2, &x3);
-		cell_set_cartesian_x(image->cell, x1*1e10, x2*1e10, x3*1e10);
+		float x1, x2, x3, d;
+		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &x1, &x2, &x3);
+		cell_set_cartesian_x(image->molecule->cell,
+		                     x1*1e-10, x2*1e-10, x3*1e-10);
 		image->dirax_read_cell++;
 		return;
 	} else if ( image->dirax_read_cell == 2 ) {
 		/* First row of unit cell values */
-		float y1, y2, y3;
-		sscanf(line, "%f %f %f", &y1, &y2, &y3);
-		cell_set_cartesian_y(image->cell, y1*1e10, y2*1e10, y3*1e10);
+		float y1, y2, y3, d;
+		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &y1, &y2, &y3);
+		cell_set_cartesian_y(image->molecule->cell,
+		                     y1*1e-10, y2*1e-10, y3*1e-10);
 		image->dirax_read_cell++;
 		return;
 	} else if ( image->dirax_read_cell == 3 ) {
 		/* First row of unit cell values */
-		float z1, z2, z3;
-		sscanf(line, "%f %f %f", &z1, &z2, &z3);
-		cell_set_cartesian_z(image->cell, z1*1e10, z2*1e10, z3*1e10);
-		STATUS("Read a reciprocal unit cell from DirAx\n");
+		float z1, z2, z3, d;
+		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &z1, &z2, &z3);
+		cell_set_cartesian_z(image->molecule->cell,
+		                     z1*1e-10, z2*1e-10, z3*1e-10);
+		STATUS("Read a direct space unit cell from DirAx\n");
 		/* FIXME: Do something */
 		image->dirax_read_cell = 0;
 		return;
@@ -172,6 +179,7 @@ static void dirax_send_next(struct image *image)
 		default: {
 			image->dirax_step = 0;
 			STATUS("DirAx is idle\n");
+			g_main_loop_quit(image->dirax_ml);
 		}
 
 	}
@@ -371,7 +379,7 @@ static int map_position(struct image *image, double x, double y,
 
 #define PEAK_WINDOW_SIZE (10)
 
-static void search_peaks(struct image *image)
+static void search_peaks(struct image *image, int dump_peaks)
 {
 	FILE *fh;
 	int x, y, width, height;
@@ -478,7 +486,9 @@ static void search_peaks(struct image *image)
 				double rz = 0.0;
 
 				/* Map and record reflection */
-				printf("%i %i\n", x, y);
+				if ( dump_peaks ) {
+					printf("%i %i\n", x, y);
+				}
 
 				image_add_feature(image->features,
 				                  mask_x, mask_y, image, 1.0);
@@ -497,14 +507,13 @@ static void search_peaks(struct image *image)
 }
 
 
-void index_pattern(struct image *image, int no_index)
+void index_pattern(struct image *image, int no_index, int dump_peaks)
 {
 	unsigned int opts;
 	int saved_stderr;
-	GMainLoop *ml;
 
 	/* Do peak search and splurge out 'xfel.drx' */
-	search_peaks(image);
+	search_peaks(image, dump_peaks);
 
 	if ( no_index ) return;
 
@@ -548,8 +557,8 @@ void index_pattern(struct image *image, int no_index)
 	g_io_add_watch(image->dirax, G_IO_IN | G_IO_HUP,
 	               (GIOFunc)dirax_readable, image);
 
-	ml = g_main_loop_new(NULL, FALSE);
-	g_main_loop_run(ml);
+	image->dirax_ml = g_main_loop_new(NULL, FALSE);
+	g_main_loop_run(image->dirax_ml);
 
 	return;
 }
diff --git a/src/dirax.h b/src/dirax.h
index f4b82ba5..4ebad677 100644
--- a/src/dirax.h
+++ b/src/dirax.h
@@ -18,7 +18,7 @@
 #endif
 
 
-extern void index_pattern(struct image *image, int no_index);
+extern void index_pattern(struct image *image, int no_index, int dump_peaks);
 
 
 #endif	/* DIRAX_H */
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 9a65fb16..008aa19f 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -415,13 +415,21 @@ static void load_features_from_file(struct image *image, const char *filename)
 
 	do {
 		char line[1024];
-		int x, y;
-		int r;
+		int x, y, r, d;
+		float df;
 
 		rval = fgets(line, 1023, fh);
 		if ( rval == NULL ) continue;
 		chomp(line);
 
+		/* Try long format (output of pattern_sim --near-bragg) */
+		r = sscanf(line, "%i %i %i %f (at %i,%i)",
+		           &d, &d, &d, &df, &x, &y);
+		if ( r == 6 ) {
+			image_add_feature(image->features, x, y, image, 1.0);
+			continue;
+		}
+
 		r = sscanf(line, "%i %i", &x, &y);
 		if ( r != 2 ) continue;
 
diff --git a/src/image.h b/src/image.h
index 1dcfdc3b..9649fe8e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -94,8 +94,6 @@ struct image {
 	ImageFeatureList	*features;	/* "Experimental" features */
 	ImageFeatureList	*rflist;	/* "Predicted" features */
 
-	UnitCell		*cell;
-
 	/* DirAx auto-indexing low-level stuff */
 	GIOChannel		*dirax;
 	int			dirax_pty;
@@ -103,6 +101,7 @@ struct image {
 	char			*dirax_rbuffer;
 	int			dirax_rbufpos;
 	int			dirax_rbuflen;
+	GMainLoop		*dirax_ml;
 
 	/* DirAx auto-indexing high-level stuff */
 	int			dirax_step;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 1ff6d2e6..41cf5b22 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -24,6 +24,8 @@
 #include "utils.h"
 #include "hdf5-file.h"
 #include "dirax.h"
+#include "intensities.h"
+#include "ewald.h"
 
 
 static void show_help(const char *s)
@@ -195,12 +197,14 @@ int main(int argc, char *argv[])
 	int n_images;
 	int n_hits;
 	int config_noindex = 0;
+	int config_dumpfound = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
 		{"input",              1, NULL,               'i'},
 		{"no-index",           0, &config_noindex,     1},
+		{"dump-found-peaks",   0, &config_dumpfound,   1},
 		{0, 0, NULL, 0}
 	};
 
@@ -257,6 +261,7 @@ int main(int argc, char *argv[])
 		chomp(line);
 
 		image.features = NULL;
+		image.molecule = NULL;
 
 		STATUS("Processing '%s'\n", line);
 
@@ -272,7 +277,18 @@ int main(int argc, char *argv[])
 		fom = image_fom(&image);
 		if ( fom > 0 ) {
 
-			index_pattern(&image, config_noindex);
+			/* Calculate orientation matrix (by magic) */
+			index_pattern(&image, config_noindex, config_dumpfound);
+
+			/* View head-on (unit cell is tilted) */
+			image.orientation.x = 0.0;
+			image.orientation.y = 0.0;
+			image.orientation.z = 0.0;
+			image.orientation.w = 1.0;
+			get_ewald(&image);
+
+			/* Read h,k,l,I */
+			output_intensities(&image);
 
 			n_hits++;