Introduce process_images

8 files changed, 705 insertions, 1 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index bb192be1..50b76850 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = pattern_sim process_hkl get_hkl
+bin_PROGRAMS = pattern_sim process_hkl get_hkl process_images
 
 if HAVE_GTK
 bin_PROGRAMS += hdfsee
@@ -16,6 +16,9 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
                       reflections.c
 process_hkl_LDADD = @LIBS@
 
+process_images_SOURCES = process_images.c hdf5-file.c utils.c dirax.c cell.c
+process_images_LDADD = @LIBS@
+
 if HAVE_GTK
 hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c
 hdfsee_LDADD = @LIBS@
diff --git a/src/cell.c b/src/cell.c
index baf110a6..94c013d1 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -146,6 +146,36 @@ void cell_set_cartesian(UnitCell *cell,
 }
 
 
+void cell_set_cartesian_x(UnitCell *cell, double ax, double bx, double cx)
+{
+	if ( !cell ) return;
+
+	cell->ax = ax;  cell->bx = bx;  cell->cx = cx;
+
+	cell_update_crystallographic(cell);
+}
+
+
+void cell_set_cartesian_y(UnitCell *cell, double ay, double by, double cy)
+{
+	if ( !cell ) return;
+
+	cell->ay = ay;  cell->by = by;  cell->cy = cy;
+
+	cell_update_crystallographic(cell);
+}
+
+
+void cell_set_cartesian_z(UnitCell *cell, double az, double bz, double cz)
+{
+	if ( !cell ) return;
+
+	cell->az = az;  cell->bz = bz;  cell->cz = cz;
+
+	cell_update_crystallographic(cell);
+}
+
+
 UnitCell *cell_new_from_parameters(double a, double b, double c,
                                    double alpha, double beta, double gamma)
 {
diff --git a/src/cell.h b/src/cell.h
index 4f96d870..ba10d638 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -60,6 +60,13 @@ extern void cell_get_cartesian(UnitCell *cell,
                                double *bx, double *by, double *bz,
                                double *cx, double *cy, double *cz);
 
+extern void cell_set_cartesian_x(UnitCell *cell,
+                                double ax, double bx, double cx);
+extern void cell_set_cartesian_y(UnitCell *cell,
+                                 double ay, double by, double cy);
+extern void cell_set_cartesian_z(UnitCell *cell,
+                                 double az, double bz, double cz);
+
 extern void cell_get_reciprocal(UnitCell *cell,
                                double *asx, double *asy, double *asz,
                                double *bsx, double *bsy, double *bsz,
diff --git a/src/dirax.c b/src/dirax.c
new file mode 100644
index 00000000..02837c9d
--- /dev/null
+++ b/src/dirax.c
@@ -0,0 +1,455 @@
+/*
+ * dirax.c
+ *
+ * Invoke the DirAx auto-indexing program
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <glib.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pty.h>
+#include <unistd.h>
+#include <sys/wait.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <sys/ioctl.h>
+#include <termio.h>
+#include <sgtty.h>
+
+
+#include "image.h"
+#include "dirax.h"
+
+
+typedef enum {
+	DIRAX_INPUT_NONE,
+	DIRAX_INPUT_LINE,
+	DIRAX_INPUT_PROMPT
+} DirAxInputType;
+
+
+static void dirax_parseline(const char *line, struct image *image)
+{
+	int i, rf;
+	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);
+
+	if ( strstr(line, "reflections from file") ) {
+		ERROR("DirAx can't understand this data.");
+		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;
+			if ( image->cell ) {
+				free(image->cell);
+			}
+			image->cell = cell_new();
+			return;
+		}
+		i++;
+	}
+
+	/* 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);
+		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);
+		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");
+		/* FIXME: Do something */
+		image->dirax_read_cell = 0;
+		return;
+	}
+
+	image->dirax_read_cell = 0;
+}
+
+
+static void dirax_sendline(const char *line, struct image *image)
+{
+	char *copy;
+	int i;
+
+	write(image->dirax_pty, line, strlen(line));
+
+	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);
+}
+
+
+static void dirax_send_next(struct image *image)
+{
+	switch ( image->dirax_step ) {
+
+		case 1 : {
+			dirax_sendline("\\echo off\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 2 : {
+			dirax_sendline("read xfel.drx\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 3 : {
+			dirax_sendline("dmax 10\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 4 : {
+			dirax_sendline("indexfit 2\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 5 : {
+			dirax_sendline("levelfit 200\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 6 : {
+			dirax_sendline("go\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		case 7 : {
+			dirax_sendline("cell\n", image);
+			image->dirax_step++;
+			break;
+		}
+
+		default: {
+			image->dirax_step = 0;
+			STATUS("DirAx is idle\n");
+		}
+
+	}
+}
+
+
+static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
+                               struct image *image)
+{
+	int rval;
+
+	rval = read(image->dirax_pty, image->dirax_rbuffer+image->dirax_rbufpos,
+			image->dirax_rbuflen-image->dirax_rbufpos);
+
+	if ( (rval == -1) || (rval == 0) ) {
+
+		ERROR("Lost connection to DirAx\n");
+		waitpid(image->dirax_pid, NULL, 0);
+		g_io_channel_shutdown(image->dirax, FALSE, NULL);
+		image->dirax = NULL;
+		return FALSE;
+
+	} else {
+
+		int no_string = 0;
+
+		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 ) {
+					if ( (strncmp(image->dirax_rbuffer+i,
+							"Dirax> ", 7) == 0)
+					  || (strncmp(image->dirax_rbuffer+i,
+					  		"PROMPT:", 7) == 0) ) {
+						block_ready = 1;
+						type = DIRAX_INPUT_PROMPT;
+						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;
+
+				switch ( type ) {
+
+					case DIRAX_INPUT_LINE : {
+
+						char *block_buffer = NULL;
+
+						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;
+
+					}
+
+					default : {
+						ERROR(
+			" Unrecognised input mode (this never happens!)\n");
+						abort();
+					}
+
+				}
+
+				/* 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 TRUE;
+}
+
+
+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;
+
+	if ( image->fmode == FORMULATION_CLEN ) {
+
+		/* Convert pixels to metres */
+		x /= image->resolution;
+		y /= image->resolution;
+		x = x * k / image->camera_len;
+		y = y * k / image->camera_len;
+		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 = x * image->pixel_size;
+		y = y * image->pixel_size;
+		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)
+{
+	FILE *fh;
+	unsigned int opts;
+	int saved_stderr;
+	int x, y;
+	GMainLoop *ml;
+
+	fh = fopen("xfel.drx", "w");
+	if ( !fh ) {
+		ERROR("Couldn't open temporary file xfel.drx\n");
+		return;
+	}
+	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  */
+
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+
+		int val;
+
+		val = image->data[y+image->height*(image->width-1-x)];
+
+		if ( val > 1000 ) {
+
+			double rx, ry, rz;
+
+			/* Map and record reflection */
+			map_position(image, x, y, &rx, &ry, &rz);
+			fprintf(fh, "%10f %10f %10f %8f\n",
+			        rx/1e10, ry/1e10, rz/1e10, 1.0);
+		}
+
+	}
+	}
+
+	fclose(fh);
+
+	saved_stderr = dup(STDERR_FILENO);
+	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);
+
+		/* Reconnect stderr */
+		dup2(saved_stderr, STDERR_FILENO);
+
+		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->dirax = g_io_channel_unix_new(image->dirax_pty);
+	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);
+
+	return;
+}
diff --git a/src/dirax.h b/src/dirax.h
new file mode 100644
index 00000000..2facc4f6
--- /dev/null
+++ b/src/dirax.h
@@ -0,0 +1,24 @@
+/*
+ * dirax.h
+ *
+ * Invoke the DirAx auto-indexing program
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef DIRAX_H
+#define DIRAX_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+extern void index_pattern(struct image *image);
+
+
+#endif	/* DIRAX_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index bd30ff7d..3f3747ca 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -246,8 +246,14 @@ int hdf5_read(struct hdfile *f, struct image *image)
 	image->data = buf;
 	image->height = f->nx;
 	image->width = f->ny;
+
+	/* FIXME: The following are basically made up... */
 	image->x_centre = image->width/2;
 	image->y_centre = image->height/2;
+	image->lambda = ph_en_to_lambda(J_to_eV(1793));
+	image->fmode = FORMULATION_CLEN;
+	image->camera_len = 75.0e-3;  /* 75 mm camera length */
+	image->resolution = 13333.3;  /* 75 micron pixel size */
 
 	return 0;
 }
diff --git a/src/image.h b/src/image.h
index b3fc7e67..65b4c31d 100644
--- a/src/image.h
+++ b/src/image.h
@@ -19,8 +19,10 @@
 
 #include <stdint.h>
 #include <complex.h>
+#include <glib.h>
 
 #include "utils.h"
+#include "cell.h"
 
 
 /* How is the scaling of the image described? */
@@ -92,6 +94,19 @@ struct image {
 	ImageFeatureList	*features;	/* "Experimental" features */
 	ImageFeatureList	*rflist;	/* "Predicted" features */
 
+	UnitCell		*cell;
+
+	/* DirAx auto-indexing low-level stuff */
+	GIOChannel		*dirax;
+	int			dirax_pty;
+	pid_t			dirax_pid;
+	char			*dirax_rbuffer;
+	int			dirax_rbufpos;
+	int			dirax_rbuflen;
+
+	/* DirAx auto-indexing high-level stuff */
+	int			dirax_step;
+	int			dirax_read_cell;
 };
 
 /* An opaque type representing a list of images */
diff --git a/src/process_images.c b/src/process_images.c
new file mode 100644
index 00000000..90f51528
--- /dev/null
+++ b/src/process_images.c
@@ -0,0 +1,164 @@
+/*
+ * process_images.c
+ *
+ * Find hits, index patterns, output hkl+intensity etc.
+ *
+ * (c) 2006-2009 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+
+#include "utils.h"
+#include "hdf5-file.h"
+#include "dirax.h"
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options]\n\n", s);
+	printf(
+"Process and index FEL diffraction images.\n"
+"\n"
+"  -h, --help              Display this help message.\n"
+"\n"
+"  -i, --input=<filename>  Specify file containing list of images to process.\n"
+"                           '-' means stdin, which is the default.\n"
+"\n");
+}
+
+
+static int sum_of_peaks(struct image *image)
+{
+	int x, y;
+	int integr = 0;
+
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+
+		int val;
+
+		val = image->data[y+image->height*(image->width-1-x)];
+
+		if ( val > 1000 ) integr+=val;
+
+	}
+	}
+
+	return integr;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	char *filename = NULL;
+	FILE *fh;
+	char *rval;
+	int n_images;
+	int n_hits;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"input",              1, NULL,               'i'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' : {
+			show_help(argv[0]);
+			return 0;
+		}
+
+		case 'i' : {
+			filename = strdup(optarg);
+			break;
+		}
+
+		case 0 : {
+			break;
+		}
+
+		default : {
+			return 1;
+		}
+		}
+
+	}
+
+	if ( filename == NULL ) {
+		filename = strdup("-");
+	}
+	if ( strcmp(filename, "-") == 0 ) {
+		fh = stdin;
+	} else {
+		fh = fopen(filename, "r");
+	}
+	free(filename);
+	if ( fh == NULL ) {
+		ERROR("Failed to open input file\n");
+		return 1;
+	}
+
+	n_images = 0;
+	n_hits = 0;
+	do {
+
+		char line[1024];
+		struct hdfile *hdfile;
+		struct image image;
+		int integr;
+
+		rval = fgets(line, 1023, fh);
+		chomp(line);
+
+		STATUS("Processing '%s'\n", line);
+
+		hdfile = hdfile_open(line);
+		if ( hdfile == NULL ) {
+			ERROR("Couldn't open file '%s'\n", filename);
+		} else if ( hdfile_set_first_image(hdfile, "/") ) {
+			ERROR("Couldn't select path\n");
+		}
+
+		hdf5_read(hdfile, &image);
+
+		integr = sum_of_peaks(&image);
+		printf("%6i %i\n", n_images, integr);
+		if ( integr > 1e7 ) {
+
+			STATUS("Hit: %s\n", line);
+
+			index_pattern(&image);
+
+			n_hits++;
+
+		}
+
+		n_images++;
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	STATUS("There were %i images.\n", n_images);
+	STATUS("%i hits were found.\n", n_hits);
+
+	return 0;
+}