aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am5
-rw-r--r--src/cell.c30
-rw-r--r--src/cell.h7
-rw-r--r--src/dirax.c455
-rw-r--r--src/dirax.h24
-rw-r--r--src/hdf5-file.c6
-rw-r--r--src/image.h15
-rw-r--r--src/process_images.c164
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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-29 13:01:39 +0000