From efb829adf0c494a3ed632cef472ac1e8e5267765 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 8 Jan 2013 17:22:02 +0100 Subject: Skeleton indexing method for GrainSpotter --- libcrystfel/Makefile.am | 5 +- libcrystfel/src/grainspotter.c | 214 +++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/grainspotter.h | 42 ++++++++ libcrystfel/src/index.c | 15 +++ libcrystfel/src/index.h | 1 + 5 files changed, 275 insertions(+), 2 deletions(-) create mode 100644 libcrystfel/src/grainspotter.c create mode 100644 libcrystfel/src/grainspotter.h diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 1c35e03f..22009223 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -8,7 +8,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ - src/cell-utils.c src/integer_matrix.c + src/cell-utils.c src/integer_matrix.c \ + src/grainspotter.c if HAVE_FFTW libcrystfel_la_SOURCES += src/reax.c @@ -24,7 +25,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ src/render.h src/index.h src/image.h \ src/filters.h src/dirax.h src/mosflm.h \ src/index-priv.h src/reax.h src/cell-utils.h \ - src/integer_matrix.h + src/integer_matrix.h src/grainspotter.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/grainspotter.c b/libcrystfel/src/grainspotter.c new file mode 100644 index 00000000..06f65f08 --- /dev/null +++ b/libcrystfel/src/grainspotter.c @@ -0,0 +1,214 @@ +/* + * grainspotter.c + * + * Invoke GrainSpotter for multi-crystal autoindexing + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if HAVE_FORKPTY_LINUX +#include +#elif HAVE_FORKPTY_BSD +#include +#endif + + +#include "image.h" +#include "grainspotter.h" +#include "utils.h" +#include "peaks.h" + + +#define GRAINSPOTTER_VERBOSE 0 + + +struct grainspotter_data { + + /* Low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + +}; + + +static int grainspotter_readable(struct image *image, + struct grainspotter_data *grainspotter) +{ + int rval; + + rval = read(grainspotter->pty, + grainspotter->rbuffer+grainspotter->rbufpos, + grainspotter->rbuflen-grainspotter->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + /* FIXME! (if needed) */ + //grainspotter->rbufpos += rval; + //assert(grainspotter->rbufpos <= grainspotter->rbuflen); + + return 0; +} + + +static void write_gve(struct image *image) +{ + FILE *fh; + int i; + char filename[1024]; + + snprintf(filename, 1023, "xfel-%i.gve", image->id); + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + return; + } + fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */ + + for ( i=0; ifeatures); i++ ) { + + struct imagefeature *f; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + fprintf(fh, "%10f %10f %10f %8f\n", + f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0); + + } + fclose(fh); +} + + +void run_grainspotter(struct image *image, UnitCell *cell) +{ + unsigned int opts; + int status; + int rval; + struct grainspotter_data *grainspotter; + + write_gve(image); + + grainspotter = malloc(sizeof(struct grainspotter_data)); + if ( grainspotter == NULL ) { + ERROR("Couldn't allocate memory for GrainSpotter data.\n"); + return; + } + + grainspotter->pid = forkpty(&grainspotter->pty, NULL, NULL, NULL); + if ( grainspotter->pid == -1 ) { + ERROR("Failed to fork for GrainSpotter: %s\n", strerror(errno)); + return; + } + if ( grainspotter->pid == 0 ) { + + /* Child process: invoke GrainSpotter */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + execlp("GrainSpotter.0.90", "", (char *)NULL); + ERROR("Failed to invoke GrainSpotter.\n"); + _exit(0); + + } + + grainspotter->rbuffer = malloc(256); + grainspotter->rbuflen = 256; + grainspotter->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(grainspotter->pty, F_GETFL); + fcntl(grainspotter->pty, F_SETFL, opts | O_NONBLOCK); + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(grainspotter->pty, &fds); + + tv.tv_sec = 30; + tv.tv_usec = 0; + + sval = select(grainspotter->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + + const int err = errno; + + switch ( err ) { + + case EINTR: + STATUS("Restarting select()\n"); + break; + + default: + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + + } + + } else if ( sval != 0 ) { + rval = grainspotter_readable(image, grainspotter); + } else { + ERROR("No response from GrainSpotter..\n"); + rval = 1; + } + + } while ( !rval ); + + close(grainspotter->pty); + free(grainspotter->rbuffer); + waitpid(grainspotter->pid, &status, 0); + + if ( status != 0 ) { + ERROR("GrainSpotter doesn't seem to be working properly.\n"); + } + + free(grainspotter); +} diff --git a/libcrystfel/src/grainspotter.h b/libcrystfel/src/grainspotter.h new file mode 100644 index 00000000..1aedf57f --- /dev/null +++ b/libcrystfel/src/grainspotter.h @@ -0,0 +1,42 @@ +/* + * grainspotter.h + * + * Invoke GrainSpotter for multi-crystal autoindexing + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifndef GRAINSPOTTER_H +#define GRAINSPOTTER_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "cell.h" + + +extern void run_grainspotter(struct image *image, UnitCell *cell); + + +#endif /* GRAINSPOTTER_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 7d912902..86ab8c1a 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -48,6 +48,7 @@ #include "index.h" #include "index-priv.h" #include "reax.h" +#include "grainspotter.h" #include "geometry.h" #include "cell-utils.h" @@ -97,6 +98,10 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, iprivs[n] = indexing_private(indm[n]); break; + case INDEXING_GRAINSPOTTER : + iprivs[n] = indexing_private(indm[n]); + break; + case INDEXING_REAX : iprivs[n] = reax_prepare(); break; @@ -132,6 +137,10 @@ void cleanup_indexing(IndexingPrivate **priv) free(priv[n]); break; + case INDEXING_GRAINSPOTTER : + free(priv[n]); + break; + case INDEXING_REAX : reax_cleanup(priv[n]); break; @@ -196,6 +205,10 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, run_mosflm(image, cell); break; + case INDEXING_GRAINSPOTTER : + run_grainspotter(image, cell); + break; + case INDEXING_REAX : reax_index(ipriv[n], image, cell); break; @@ -294,6 +307,8 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) list[i] = INDEXING_DIRAX; } else if ( strcmp(methods[i], "mosflm") == 0) { list[i] = INDEXING_MOSFLM; + } else if ( strcmp(methods[i], "grainspotter") == 0) { + list[i] = INDEXING_GRAINSPOTTER; } else if ( strcmp(methods[i], "reax") == 0) { list[i] = INDEXING_REAX; *need_cell = 1; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 4a724645..2df7a311 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -56,6 +56,7 @@ typedef enum { INDEXING_NONE, INDEXING_DIRAX, INDEXING_MOSFLM, + INDEXING_GRAINSPOTTER, INDEXING_REAX, } IndexingMethod; -- cgit v1.2.3 From 6c970262e0d41252b0ada1a80a6c6293ddcdfec9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 9 Jan 2013 16:00:59 +0100 Subject: GrainSpotter interface (works) --- libcrystfel/src/grainspotter.c | 165 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 152 insertions(+), 13 deletions(-) diff --git a/libcrystfel/src/grainspotter.c b/libcrystfel/src/grainspotter.c index 06f65f08..b7e5ce0a 100644 --- a/libcrystfel/src/grainspotter.c +++ b/libcrystfel/src/grainspotter.c @@ -50,9 +50,10 @@ #include "image.h" -#include "grainspotter.h" #include "utils.h" #include "peaks.h" +#include "cell.h" +#include "cell-utils.h" #define GRAINSPOTTER_VERBOSE 0 @@ -70,30 +71,99 @@ struct grainspotter_data { }; +static int read_matrix(struct image *image) +{ + FILE *fh; + int d1; + float d2; + float ubi11, ubi12, ubi13; + float ubi21, ubi22, ubi23; + float ubi31, ubi32, ubi33; + char filename[1024]; + char line[1024]; + int r; + + snprintf(filename, 1023, "xfel-%i.gff", image->id); + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Can't open '%s'\n", filename); + return 1; + } + + /* Read and discard first line */ + if ( fgets(line, 1024, fh) == NULL ) { + ERROR("Failed to read GFF file.\n"); + return 1; + } + + /* One line per grain */ + if ( fgets(line, 1024, fh) == NULL ) { + ERROR("Failed to read GFF file.\n"); + return 1; + } + + STATUS("'%s'\n", line); + + r = sscanf(line, + "%i %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f", + &d1, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, + &ubi11, &ubi12, &ubi13, &ubi21, &ubi22, &ubi23, &ubi31, &ubi32, &ubi33); + + if ( r != 24 ) { + ERROR("Only %i parameters in GFF file\n", r); + return 1; + } + + fclose(fh); + + image->candidate_cells[0] = cell_new(); + +// cell_set_cartesian(image->candidate_cells[0], +// ubi11/1e10, ubi12/1e10, ubi13/1e10, +// ubi21/1e10, ubi22/1e10, ubi23/1e10, +// ubi31/1e10, ubi32/1e10, ubi33/1e10); + + cell_set_cartesian(image->candidate_cells[0], + ubi12/1e10, ubi13/1e10, ubi11/1e10, + ubi22/1e10, ubi23/1e10, ubi21/1e10, + ubi32/1e10, ubi33/1e10, ubi31/1e10); + + + image->ncells = 1; + cell_print(image->candidate_cells[0]); + + return 0; +} + + static int grainspotter_readable(struct image *image, struct grainspotter_data *grainspotter) { int rval; rval = read(grainspotter->pty, - grainspotter->rbuffer+grainspotter->rbufpos, - grainspotter->rbuflen-grainspotter->rbufpos); + grainspotter->rbuffer, grainspotter->rbuflen); - if ( (rval == -1) || (rval == 0) ) return 1; + if ( rval == -1 ) { + ERROR("Read failed: %s\n", strerror(errno)); + return 1; + } + + grainspotter->rbuffer[rval] = '\0'; - /* FIXME! (if needed) */ - //grainspotter->rbufpos += rval; - //assert(grainspotter->rbufpos <= grainspotter->rbuflen); + printf("GrainSpotter: '%s'\n", grainspotter->rbuffer); return 0; } -static void write_gve(struct image *image) +static void write_gve(struct image *image, UnitCell *cell) { FILE *fh; int i; char filename[1024]; + double a, b, c, al, be, ga; snprintf(filename, 1023, "xfel-%i.gve", image->id); @@ -102,7 +172,14 @@ static void write_gve(struct image *image) ERROR("Couldn't open temporary file '%s'\n", filename); return; } - fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */ + + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + fprintf(fh, "%.6f %.6f %.6f %.6f %.6f %.6f P\n", a*1e10, b*1e10, c*1e10, + rad2deg(al), rad2deg(be), rad2deg(ga)); + fprintf(fh, "# wavelength = %.6f\n", image->lambda*1e10); + fprintf(fh, "# wedge = 0.000000\n"); + fprintf(fh, "# ds h k l\n"); + fprintf(fh, "# xr yr zr xc yc ds eta omega\n"); for ( i=0; ifeatures); i++ ) { @@ -111,11 +188,60 @@ static void write_gve(struct image *image) f = image_get_feature(image->features, i); if ( f == NULL ) continue; - fprintf(fh, "%10f %10f %10f %8f\n", - f->rx/1e10, f->ry/1e10, f->rz/1e10, 1.0); + fprintf(fh, "%.6f %.6f %.6f 0 0 %.6f %.6f %.6f 0\n", + f->rz/1e10, f->rx/1e10, f->ry/1e10, + modulus(f->rx, f->ry, f->rz)/1e10, /* dstar */ + atan2(f->ry, f->rx), 0.0); /* eta, omega */ + + } + fclose(fh); +} + + +static char *write_ini(struct image *image) +{ + FILE *fh; + char *filename; + double tt; + + filename = malloc(1024); + if ( filename == NULL ) return NULL; + snprintf(filename, 1023, "xfel-%i.ini", image->id); + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + free(filename); + return NULL; } + + get_q_for_panel(image->det->furthest_out_panel, + image->det->furthest_out_fs, + image->det->furthest_out_ss, + &tt, 1.0/image->lambda); + + fprintf(fh, "spacegroup 96\n"); + fprintf(fh, "!dsrange 0 1.3\n"); + fprintf(fh, "tthrange 0 %.2f\n", 20.0); + fprintf(fh, "etarange 0 360\n"); + fprintf(fh, "domega 1\n"); + fprintf(fh, "omegarange -0.5 0.5\n"); + fprintf(fh, "filespecs xfel-%i.gve xfel-%i.log\n", + image->id, image->id); + fprintf(fh, "cuts 3 0.1 0.5\n"); + fprintf(fh, "eulerstep 8\n"); + fprintf(fh, "uncertainties 0.1 0.2 .5\n"); + fprintf(fh, "nsigmas 2\n"); + fprintf(fh, "minfracg 0.95\n"); + fprintf(fh, "Nhkls_in_indexing 500\n"); + fprintf(fh, "random 10000\n"); + fprintf(fh, "!positionfit\n"); + fprintf(fh, "genhkl\n"); + fclose(fh); + + return filename; } @@ -125,8 +251,15 @@ void run_grainspotter(struct image *image, UnitCell *cell) int status; int rval; struct grainspotter_data *grainspotter; + char *ini_filename; + + write_gve(image, cell); + ini_filename = write_ini(image); - write_gve(image); + if ( ini_filename == NULL ) { + ERROR("Failed to write ini file for GrainSpotter.\n"); + return; + } grainspotter = malloc(sizeof(struct grainspotter_data)); if ( grainspotter == NULL ) { @@ -149,11 +282,13 @@ void run_grainspotter(struct image *image, UnitCell *cell) t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); tcsetattr(STDIN_FILENO, TCSANOW, &t); - execlp("GrainSpotter.0.90", "", (char *)NULL); + STATUS("Running GrainSpotter.0.90 '%s'\n", ini_filename); + execlp("GrainSpotter.0.90", "", ini_filename, (char *)NULL); ERROR("Failed to invoke GrainSpotter.\n"); _exit(0); } + free(ini_filename); grainspotter->rbuffer = malloc(256); grainspotter->rbuflen = 256; @@ -210,5 +345,9 @@ void run_grainspotter(struct image *image, UnitCell *cell) ERROR("GrainSpotter doesn't seem to be working properly.\n"); } + if ( read_matrix(image) != 0 ) { + ERROR("Failed to read matrix\n"); + } + free(grainspotter); } -- cgit v1.2.3 From 60726ebd9b939a17778f785206c9eb00d5e6c800 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 9 Jan 2013 17:06:40 +0100 Subject: Matrix reading stuff --- libcrystfel/src/grainspotter.c | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/libcrystfel/src/grainspotter.c b/libcrystfel/src/grainspotter.c index b7e5ce0a..6217b810 100644 --- a/libcrystfel/src/grainspotter.c +++ b/libcrystfel/src/grainspotter.c @@ -71,7 +71,7 @@ struct grainspotter_data { }; -static int read_matrix(struct image *image) +static int read_matrix(struct image *image, char *filename) { FILE *fh; int d1; @@ -79,12 +79,9 @@ static int read_matrix(struct image *image) float ubi11, ubi12, ubi13; float ubi21, ubi22, ubi23; float ubi31, ubi32, ubi33; - char filename[1024]; char line[1024]; int r; - snprintf(filename, 1023, "xfel-%i.gff", image->id); - fh = fopen(filename, "r"); if ( fh == NULL ) { ERROR("Can't open '%s'\n", filename); @@ -252,6 +249,7 @@ void run_grainspotter(struct image *image, UnitCell *cell) int rval; struct grainspotter_data *grainspotter; char *ini_filename; + char gff_filename[1024]; write_gve(image, cell); ini_filename = write_ini(image); @@ -267,6 +265,9 @@ void run_grainspotter(struct image *image, UnitCell *cell) return; } + snprintf(gff_filename, 1023, "xfel-%i.gff", image->id); + remove(gff_filename); + grainspotter->pid = forkpty(&grainspotter->pty, NULL, NULL, NULL); if ( grainspotter->pid == -1 ) { ERROR("Failed to fork for GrainSpotter: %s\n", strerror(errno)); @@ -345,7 +346,7 @@ void run_grainspotter(struct image *image, UnitCell *cell) ERROR("GrainSpotter doesn't seem to be working properly.\n"); } - if ( read_matrix(image) != 0 ) { + if ( read_matrix(image, gff_filename) != 0 ) { ERROR("Failed to read matrix\n"); } -- cgit v1.2.3 From 76d3f48ddf33cc82dc2aaa52c43463dd8bb9cc70 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 14 Jan 2013 14:22:54 +0100 Subject: Formatting --- libcrystfel/src/grainspotter.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libcrystfel/src/grainspotter.c b/libcrystfel/src/grainspotter.c index 6217b810..7d0f74ce 100644 --- a/libcrystfel/src/grainspotter.c +++ b/libcrystfel/src/grainspotter.c @@ -214,9 +214,9 @@ static char *write_ini(struct image *image) } get_q_for_panel(image->det->furthest_out_panel, - image->det->furthest_out_fs, - image->det->furthest_out_ss, - &tt, 1.0/image->lambda); + image->det->furthest_out_fs, + image->det->furthest_out_ss, + &tt, 1.0/image->lambda); fprintf(fh, "spacegroup 96\n"); fprintf(fh, "!dsrange 0 1.3\n"); -- cgit v1.2.3