aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/Makefile.am10
-rw-r--r--src/index.c9
-rw-r--r--src/index.h1
-rw-r--r--src/indexamajig.c3
-rw-r--r--src/mosflm.c470
-rw-r--r--src/mosflm.h28
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 */