added fake mosflm support (runs dirax)

6 files changed, 516 insertions, 5 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 9d0100d5..dfd82d1f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -25,8 +25,8 @@ process_hkl_LDADD = @LIBS@
 
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
                       peaks.c index.c filters.c diffraction.c detector.c \
-                      sfac.c dirax.c reflections.c templates.c symmetry.c \
-                      geometry.c thread-pool.c beam-parameters.c
+                      sfac.c dirax.c mosflm.c reflections.c templates.c \
+                      symmetry.c geometry.c thread-pool.c beam-parameters.c
 indexamajig_LDADD = @LIBS@
 if HAVE_OPENCL
 indexamajig_SOURCES += diffraction-gpu.c cl-utils.c
@@ -76,8 +76,8 @@ endif
 
 reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \
                       peaks.c image.c stream.c \
-                      index.c dirax.c templates.c geometry.c symmetry.c \
-                      thread-pool.c
+                      index.c dirax.c mosflm.c templates.c geometry.c \
+                      symmetry.c thread-pool.c
 reintegrate_LDADD = @LIBS@
 
 estimate_background_SOURCES = estimate_background.c stream.c utils.c cell.c
@@ -87,7 +87,7 @@ INCLUDES = "-I$(top_srcdir)/data"
 
 EXTRA_DIST = cell.h hdf5-file.h image.h  utils.h diffraction.h detector.h \
              sfac.h reflections.h list_tmp.h statistics.h displaywindow.h \
-             render.h hdfsee.h dirax.h peaks.h index.h  filters.h \
+             render.h hdfsee.h dirax.h mosflm.h peaks.h index.h  filters.h \
              diffraction-gpu.h cl-utils.h symmetry.h \
              povray.h index-priv.h geometry.h templates.h render_hkl.h \
              stream.h thread-pool.h beam-parameters.h post-refinement.h \
diff --git a/src/index.c b/src/index.c
index bb725fcb..3e0a09b1 100644
--- a/src/index.c
+++ b/src/index.c
@@ -24,6 +24,7 @@
 #include "utils.h"
 #include "peaks.h"
 #include "dirax.h"
+#include "mosflm.h"
 #include "sfac.h"
 #include "detector.h"
 #include "index.h"
@@ -50,6 +51,8 @@ IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
 		return indexing_private(indm);
 	case INDEXING_DIRAX :
 		return indexing_private(indm);
+	case INDEXING_MOSFLM :
+		return indexing_private(indm);
 	case INDEXING_TEMPLATE :
 		return generate_templates(cell, filename, det,
 		                          nominal_photon_energy);
@@ -67,6 +70,9 @@ void cleanup_indexing(IndexingPrivate *priv)
 	case INDEXING_DIRAX :
 		free(priv);
 		break;
+	case INDEXING_MOSFLM :
+		free(priv);
+		break;
 	case INDEXING_TEMPLATE :
 		free_templates(priv);
 	}
@@ -140,6 +146,9 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
 	case INDEXING_DIRAX :
 		run_dirax(image);
 		break;
+	case INDEXING_MOSFLM :
+		run_mosflm(image);
+		break;
 	case INDEXING_TEMPLATE :
 		match_templates(image, ipriv);
 		break;
diff --git a/src/index.h b/src/index.h
index 81549e00..29be9dcc 100644
--- a/src/index.h
+++ b/src/index.h
@@ -27,6 +27,7 @@
 typedef enum {
 	INDEXING_NONE,
 	INDEXING_DIRAX,
+	INDEXING_MOSFLM,
 	INDEXING_TEMPLATE
 } IndexingMethod;
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 5a33e74d..8d762fa4 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -127,6 +127,7 @@ static void show_help(const char *s)
 "     --indexing=<method>  Use 'method' for indexing.  Choose from:\n"
 "                           none     : no indexing (default)\n"
 "                           dirax    : invoke DirAx\n"
+"                           mosflm   : invoke MOSFLM (DPS)\n"
 "                           template : index by template matching\n"
 " -g. --geometry=<file>    Get detector geometry from file.\n"
 " -b, --beam=<file>        Get beam parameters from file (provides nominal\n"
@@ -788,6 +789,8 @@ int main(int argc, char *argv[])
 		indm = INDEXING_NONE;
 	} else if ( strcmp(indm_str, "dirax") == 0) {
 		indm = INDEXING_DIRAX;
+	} else if ( strcmp(indm_str, "mosflm") == 0) {
+		indm = INDEXING_MOSFLM;
 	} else if ( strcmp(indm_str, "template") == 0) {
 		indm = INDEXING_TEMPLATE;
 	} else {
diff --git a/src/mosflm.c b/src/mosflm.c
new file mode 100644
index 00000000..aed22720
--- /dev/null
+++ b/src/mosflm.c
@@ -0,0 +1,470 @@
+/*
+ * mosflm.c
+ *
+ * Invoke the DPS auto-indexing algorithm through MOSFLM
+ *
+ * This will actuaully run DirAx for now... will be fixed soon...
+ *
+ * (c) 2006-2010 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 <unistd.h>
+#include <sys/wait.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <sys/ioctl.h>
+
+#if HAVE_FORKPTY_LINUX
+#include <pty.h>
+#elif HAVE_FORKPTY_BSD
+#include <util.h>
+#endif
+
+
+#include "image.h"
+#include "mosflm.h"
+#include "utils.h"
+#include "sfac.h"
+#include "peaks.h"
+
+
+
+
+
+
+#define DIRAX_VERBOSE 0
+
+#define MAX_DIRAX_CELL_CANDIDATES (5)
+
+
+typedef enum {
+	DIRAX_INPUT_NONE,
+	DIRAX_INPUT_LINE,
+	DIRAX_INPUT_PROMPT,
+	DIRAX_INPUT_ACL
+} DirAxInputType;
+
+
+static void dirax_parseline(const char *line, struct image *image)
+{
+	int rf, i, di, acl, acl_nh;
+	float d;
+
+	#if DIRAX_VERBOSE
+	char *copy;
+
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='r';
+		if ( copy[i] == '\n' ) copy[i]='\0';
+	}
+	STATUS("DirAx: %s\n", copy);
+	free(copy);
+	#endif
+
+	if ( strstr(line, "reflections from file") ) {
+		ERROR("DirAx can't understand this data.\n");
+		return;
+	}
+
+	/* Is this the first line of a unit cell specification? */
+	rf = 0; i = 0;
+	while ( (i<strlen(line)) && ((line[i] == 'R')
+				|| (line[i] == 'D') || (line[i] == ' ')) ) {
+		if ( line[i] == 'R' ) rf = 1;
+		if ( (line[i] == 'D') && rf ) {
+			image->dirax_read_cell = 1;
+			image->candidate_cells[image->ncells] = cell_new();
+			return;
+		}
+		i++;
+	}
+
+	/* Parse unit cell vectors as appropriate */
+	if ( image->dirax_read_cell == 1 ) {
+		/* First row of unit cell values */
+		float ax, ay, az;
+		int r;
+		r = sscanf(line, "%f %f %f %f %f %f",
+		           &d, &d, &d, &ax, &ay, &az);
+		if ( r != 6 ) {
+			ERROR("Couldn't understand cell line\n");
+			image->dirax_read_cell = 0;
+			free(image->candidate_cells[image->ncells]);
+			return;
+		}
+		cell_set_cartesian_a(image->candidate_cells[image->ncells],
+		                     ax*1e-10, ay*1e-10, az*1e-10);
+		image->dirax_read_cell++;
+		return;
+	} else if ( image->dirax_read_cell == 2 ) {
+		/* Second row of unit cell values */
+		float bx, by, bz;
+		int r;
+		r = sscanf(line, "%f %f %f %f %f %f",
+		           &d, &d, &d, &bx, &by, &bz);
+		if ( r != 6 ) {
+			ERROR("Couldn't understand cell line\n");
+			image->dirax_read_cell = 0;
+			free(image->candidate_cells[image->ncells]);
+			return;
+		}
+		cell_set_cartesian_b(image->candidate_cells[image->ncells],
+		                     bx*1e-10, by*1e-10, bz*1e-10);
+		image->dirax_read_cell++;
+		return;
+	} else if ( image->dirax_read_cell == 3 ) {
+		/* Third row of unit cell values */
+		float cx, cy, cz;
+		int r;
+		r = sscanf(line, "%f %f %f %f %f %f",
+		           &d, &d, &d, &cx, &cy, &cz);
+		if ( r != 6 ) {
+			ERROR("Couldn't understand cell line\n");
+			image->dirax_read_cell = 0;
+			free(image->candidate_cells[image->ncells]);
+			return;
+		}
+		cell_set_cartesian_c(image->candidate_cells[image->ncells++],
+		                     cx*1e-10, cy*1e-10, cz*1e-10);
+		image->dirax_read_cell = 0;
+		return;
+	}
+
+	image->dirax_read_cell = 0;
+
+	if ( sscanf(line, "%i %i %f %f %f %f %f %f %i", &acl, &acl_nh,
+	                               &d, &d, &d, &d, &d, &d, &di) == 9 ) {
+		if ( acl_nh > image->best_acl_nh ) {
+
+			int i, found = 0;
+
+			for ( i=0; i<image->n_acls_tried; i++ ) {
+				if ( image->acls_tried[i] == acl ) found = 1;
+			}
+
+			if ( !found ) {
+				image->best_acl_nh = acl_nh;
+				image->best_acl = acl;
+			}
+
+		}
+	}
+}
+
+
+static void dirax_sendline(const char *line, struct image *image)
+{
+	#if DIRAX_VERBOSE
+	char *copy;
+	int i;
+
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='\0';
+		if ( copy[i] == '\n' ) copy[i]='\0';
+	}
+	STATUS("To DirAx: '%s'\n", copy);
+	free(copy);
+	#endif
+
+	write(image->dirax_pty, line, strlen(line));
+}
+
+
+static void dirax_send_next(struct image *image)
+{
+	char tmp[32];
+
+	switch ( image->dirax_step ) {
+
+	case 1 :
+		dirax_sendline("\\echo off\n", image);
+		break;
+
+	case 2 :
+		snprintf(tmp, 31, "read xfel-%i.drx\n", image->id);
+		dirax_sendline(tmp, image);
+		break;
+
+	case 3 :
+		dirax_sendline("dmax 1000\n", image);
+		break;
+
+	case 4 :
+		dirax_sendline("indexfit 2\n", image);
+		break;
+
+	case 5 :
+		dirax_sendline("levelfit 1000\n", image);
+		break;
+
+	case 6 :
+		dirax_sendline("go\n", image);
+		break;
+
+	case 7 :
+		dirax_sendline("acl\n", image);
+		break;
+
+	case 8 :
+		if ( image->best_acl_nh == 0 ) {
+			STATUS("No more cells to try.\n");
+			/* At this point, DirAx is presenting its ACL prompt
+			 * and waiting for a single number.  Use an extra
+			 * newline to choose automatic ACL selection before
+			 * exiting. */
+			dirax_sendline("\nexit\n", image);
+			break;
+		}
+		snprintf(tmp, 31, "%i\n", image->best_acl);
+		image->acls_tried[image->n_acls_tried++] = image->best_acl;
+		dirax_sendline(tmp, image);
+		break;
+
+	case 9 :
+		dirax_sendline("cell\n", image);
+		break;
+
+	case 10 :
+		if ( image->n_acls_tried == MAX_DIRAX_CELL_CANDIDATES ) {
+			dirax_sendline("exit\n", image);
+		} else {
+			/* Go back round for another cell */
+			image->best_acl_nh = 0;
+			image->dirax_step = 7;
+			dirax_sendline("acl\n", image);
+		}
+		break;
+
+	default:
+		dirax_sendline("exit\n", image);
+		return;
+
+	}
+
+	image->dirax_step++;
+}
+
+
+static int dirax_readable(struct image *image)
+{
+	int rval;
+	int no_string = 0;
+
+	rval = read(image->dirax_pty, image->dirax_rbuffer+image->dirax_rbufpos,
+	            image->dirax_rbuflen-image->dirax_rbufpos);
+
+	if ( (rval == -1) || (rval == 0) ) return 1;
+
+	image->dirax_rbufpos += rval;
+	assert(image->dirax_rbufpos <= image->dirax_rbuflen);
+
+	while ( (!no_string) && (image->dirax_rbufpos > 0) ) {
+
+		int i;
+		int block_ready = 0;
+		DirAxInputType type = DIRAX_INPUT_NONE;
+
+		/* See if there's a full line in the buffer yet */
+		for ( i=0; i<image->dirax_rbufpos-1; i++ ) {
+			/* Means the last value looked at is rbufpos-2 */
+
+			/* Is there a prompt in the buffer? */
+			if ( (i+7 <= image->dirax_rbufpos)
+			  && (!strncmp(image->dirax_rbuffer+i, "Dirax> ", 7)) ) {
+				block_ready = 1;
+				type = DIRAX_INPUT_PROMPT;
+				break;
+			}
+
+			/* How about an ACL prompt? */
+			if ( (i+10 <= image->dirax_rbufpos)
+			  && (!strncmp(image->dirax_rbuffer+i, "acl/auto [", 10)) ) {
+				block_ready = 1;
+				type = DIRAX_INPUT_ACL;
+				break;
+			}
+
+			if ( (image->dirax_rbuffer[i] == '\r')
+			  && (image->dirax_rbuffer[i+1] == '\n') ) {
+				block_ready = 1;
+				type = DIRAX_INPUT_LINE;
+				break;
+			}
+
+		}
+
+		if ( block_ready ) {
+
+			unsigned int new_rbuflen;
+			unsigned int endbit_length;
+			char *block_buffer = NULL;
+
+			switch ( type ) {
+
+			case DIRAX_INPUT_LINE :
+
+				block_buffer = malloc(i+1);
+				memcpy(block_buffer, image->dirax_rbuffer, i);
+				block_buffer[i] = '\0';
+
+				if ( block_buffer[0] == '\r' ) {
+					memmove(block_buffer, block_buffer+1, i);
+				}
+
+				dirax_parseline(block_buffer, image);
+				free(block_buffer);
+				endbit_length = i+2;
+				break;
+
+			case DIRAX_INPUT_PROMPT :
+
+				dirax_send_next(image);
+				endbit_length = i+7;
+				break;
+
+			case DIRAX_INPUT_ACL :
+
+				dirax_send_next(image);
+				endbit_length = i+10;
+				break;
+
+			default :
+
+				/* Obviously, this never happens :) */
+				ERROR("Unrecognised DirAx input mode! "
+				      "I don't know how to understand DirAx\n");
+				return 1;
+
+			}
+
+			/* Now the block's been parsed, it should be
+			 * forgotten about */
+			memmove(image->dirax_rbuffer,
+			        image->dirax_rbuffer + endbit_length,
+			        image->dirax_rbuflen - endbit_length);
+
+			/* Subtract the number of bytes removed */
+			image->dirax_rbufpos = image->dirax_rbufpos
+			                       - endbit_length;
+			new_rbuflen = image->dirax_rbuflen - endbit_length;
+			if ( new_rbuflen == 0 ) new_rbuflen = 256;
+			image->dirax_rbuffer = realloc(image->dirax_rbuffer,
+			                               new_rbuflen);
+			image->dirax_rbuflen = new_rbuflen;
+
+		} else {
+
+			if ( image->dirax_rbufpos==image->dirax_rbuflen ) {
+
+				/* More buffer space is needed */
+				image->dirax_rbuffer = realloc(
+				                    image->dirax_rbuffer,
+				                    image->dirax_rbuflen + 256);
+				image->dirax_rbuflen = image->dirax_rbuflen
+				                                          + 256;
+				/* The new space gets used at the next
+				 * read, shortly... */
+
+			}
+			no_string = 1;
+
+		}
+
+	}
+
+	return 0;
+}
+
+
+
+
+
+void run_mosflm(struct image *image)
+{
+	unsigned int opts;
+	int status;
+	int rval;
+
+	printf("Mosflm is not yet implemented.  Using DirAx for now.\n");
+
+	image->dirax_pid = forkpty(&image->dirax_pty, NULL, NULL, NULL);
+	if ( image->dirax_pid == -1 ) {
+		ERROR("Failed to fork for DirAx\n");
+		return;
+	}
+	if ( image->dirax_pid == 0 ) {
+
+		/* Child process: invoke DirAx */
+		struct termios t;
+
+		/* Turn echo off */
+		tcgetattr(STDIN_FILENO, &t);
+		t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+		tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+		execlp("dirax", "", (char *)NULL);
+		ERROR("Failed to invoke DirAx.\n");
+		_exit(0);
+
+	}
+
+	image->dirax_rbuffer = malloc(256);
+	image->dirax_rbuflen = 256;
+	image->dirax_rbufpos = 0;
+
+	/* Set non-blocking */
+	opts = fcntl(image->dirax_pty, F_GETFL);
+	fcntl(image->dirax_pty, F_SETFL, opts | O_NONBLOCK);
+
+	image->dirax_step = 1;	/* This starts the "initialisation" procedure */
+	image->dirax_read_cell = 0;
+	image->n_acls_tried = 0;
+	image->best_acl_nh = 0;
+
+	do {
+
+		fd_set fds;
+		struct timeval tv;
+		int sval;
+
+		FD_ZERO(&fds);
+		FD_SET(image->dirax_pty, &fds);
+
+		tv.tv_sec = 10;
+		tv.tv_usec = 0;
+
+		sval = select(image->dirax_pty+1, &fds, NULL, NULL, &tv);
+
+		if ( sval == -1 ) {
+			ERROR("select() failed.\n");
+			rval = 1;
+		} else if ( sval != 0 ) {
+			rval = dirax_readable(image);
+		} else {
+			ERROR("No response from DirAx..\n");
+			rval = 1;
+		}
+
+	} while ( !rval );
+
+	close(image->dirax_pty);
+	free(image->dirax_rbuffer);
+	wait(&status);
+
+	return;
+}
diff --git a/src/mosflm.h b/src/mosflm.h
new file mode 100644
index 00000000..8d09d8d5
--- /dev/null
+++ b/src/mosflm.h
@@ -0,0 +1,28 @@
+/*
+ * mosflm.h
+ *
+ * Invoke the DPS auto-indexing algorithm through MOSFLM
+ *
+ * This will actuaully run DirAx for now... will be fixed soon...
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef MOSFLM_H
+#define MOSFLM_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "utils.h"
+
+
+extern void run_mosflm(struct image *image);
+
+
+#endif	/* MOSFLM_H */