diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | configure.ac | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 5 | ||||
| -rw-r--r-- | src/cell.c | 30 | ||||
| -rw-r--r-- | src/cell.h | 7 | ||||
| -rw-r--r-- | src/dirax.c | 455 | ||||
| -rw-r--r-- | src/dirax.h | 24 | ||||
| -rw-r--r-- | src/hdf5-file.c | 6 | ||||
| -rw-r--r-- | src/image.h | 15 | ||||
| -rw-r--r-- | src/process_images.c | 164 |
11 files changed, 708 insertions, 3 deletions
@@ -13,4 +13,5 @@ src/pattern_sim src/process_hkl src/get_hkl src/hdfsee +src/process_images *~ diff --git a/Makefile.am b/Makefile.am index 8a068240..f6fa5efe 100644 --- a/Makefile.am +++ b/Makefile.am @@ -2,5 +2,5 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h src/relrod.h \ src/utils.h src/diffraction.h src/detector.h src/ewald.h \ src/sfac.h src/intensities.h src/reflections.h src/list_tmp.h \ src/statistics.h src/displaywindow.h src/render.h src/hdfsee.h \ - data/displaywindow.ui + data/displaywindow.ui src/dirax.h SUBDIRS = src data diff --git a/configure.ac b/configure.ac index 7cdf50e5..cce0fc00 100644 --- a/configure.ac +++ b/configure.ac @@ -31,6 +31,6 @@ AC_MSG_WARN([GTK not found. hdfsee will not be built.])) AM_CONDITIONAL([HAVE_GTK], test x$havegtk = xtrue) CFLAGS="$CFLAGS $HDF5_CFLAGS $GTK_CFLAGS" -LIBS="$LIBS $HDF5_LIBS -lm -lz -lgsl -lgslcblas $GTK_LIBS -lgthread-2.0" +LIBS="$LIBS $HDF5_LIBS -lm -lz -lgsl -lgslcblas $GTK_LIBS -lgthread-2.0 -lutil" AC_OUTPUT(Makefile src/Makefile data/Makefile) 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@ @@ -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) { @@ -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; +} |