From 25c3d29ed7701cadbb3813878f25b633a7cd7c2d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 13:59:17 +0100 Subject: Move indexing and rendering to libcrystfel --- libcrystfel/src/dirax.c | 529 +++++++++++++++++++++++++++++++ libcrystfel/src/dirax.h | 26 ++ libcrystfel/src/index-priv.h | 29 ++ libcrystfel/src/index.c | 274 ++++++++++++++++ libcrystfel/src/index.h | 63 ++++ libcrystfel/src/mosflm.c | 609 ++++++++++++++++++++++++++++++++++++ libcrystfel/src/mosflm.h | 27 ++ libcrystfel/src/reax.c | 728 +++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/reax.h | 27 ++ libcrystfel/src/render.c | 442 ++++++++++++++++++++++++++ libcrystfel/src/render.h | 55 ++++ 11 files changed, 2809 insertions(+) create mode 100644 libcrystfel/src/dirax.c create mode 100644 libcrystfel/src/dirax.h create mode 100644 libcrystfel/src/index-priv.h create mode 100644 libcrystfel/src/index.c create mode 100644 libcrystfel/src/index.h create mode 100644 libcrystfel/src/mosflm.c create mode 100644 libcrystfel/src/mosflm.h create mode 100644 libcrystfel/src/reax.c create mode 100644 libcrystfel/src/reax.h create mode 100644 libcrystfel/src/render.c create mode 100644 libcrystfel/src/render.h (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c new file mode 100644 index 00000000..bb1cb808 --- /dev/null +++ b/libcrystfel/src/dirax.c @@ -0,0 +1,529 @@ +/* + * dirax.c + * + * Invoke the DirAx auto-indexing program + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#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 "dirax.h" +#include "utils.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; + + +struct dirax_data { + + /* DirAx auto-indexing low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + + /* DirAx auto-indexing high-level stuff */ + int step; + int finished_ok; + int read_cell; + int best_acl; + int best_acl_nh; + int acls_tried[MAX_CELL_CANDIDATES]; + int n_acls_tried; + +}; + + +static void dirax_parseline(const char *line, struct image *image, + struct dirax_data *dirax) +{ + int rf, i, di, acl, acl_nh; + float d; + + #if DIRAX_VERBOSE + char *copy; + + copy = strdup(line); + for ( i=0; iread_cell = 1; + image->candidate_cells[image->ncells] = cell_new(); + return; + } + i++; + } + + /* Parse unit cell vectors as appropriate */ + if ( 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"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_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); + dirax->read_cell++; + return; + } else if ( 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"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_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); + dirax->read_cell++; + return; + } else if ( 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"); + ERROR("'%s'\n", line); + dirax->read_cell = 0; + cell_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); + dirax->read_cell = 0; + return; + } + + 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 > dirax->best_acl_nh ) { + + int i, found = 0; + + for ( i=0; in_acls_tried; i++ ) { + if ( dirax->acls_tried[i] == acl ) found = 1; + } + + if ( !found ) { + dirax->best_acl_nh = acl_nh; + dirax->best_acl = acl; + } + + } + } +} + + +static void dirax_sendline(const char *line, struct dirax_data *dirax) +{ + #if DIRAX_VERBOSE + char *copy; + int i; + + copy = strdup(line); + for ( i=0; ipty, line, strlen(line)); +} + + +static void dirax_send_next(struct image *image, struct dirax_data *dirax) +{ + char tmp[32]; + + switch ( dirax->step ) { + + case 1 : + dirax_sendline("\\echo off\n", dirax); + break; + + case 2 : + snprintf(tmp, 31, "read xfel-%i.drx\n", image->id); + dirax_sendline(tmp, dirax); + break; + + case 3 : + dirax_sendline("dmax 1000\n", dirax); + break; + + case 4 : + dirax_sendline("indexfit 2\n", dirax); + break; + + case 5 : + dirax_sendline("levelfit 1000\n", dirax); + break; + + case 6 : + dirax_sendline("go\n", dirax); + dirax->finished_ok = 1; + break; + + case 7 : + dirax_sendline("acl\n", dirax); + break; + + case 8 : + if ( dirax->best_acl_nh == 0 ) { + /* 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", dirax); + break; + } + snprintf(tmp, 31, "%i\n", dirax->best_acl); + dirax->acls_tried[dirax->n_acls_tried++] = dirax->best_acl; + dirax_sendline(tmp, dirax); + break; + + case 9 : + dirax_sendline("cell\n", dirax); + break; + + case 10 : + if ( dirax->n_acls_tried == MAX_DIRAX_CELL_CANDIDATES ) { + dirax_sendline("exit\n", dirax); + } else { + /* Go back round for another cell */ + dirax->best_acl_nh = 0; + dirax->step = 7; + dirax_sendline("acl\n", dirax); + } + break; + + default: + dirax_sendline("exit\n", dirax); + return; + + } + + dirax->step++; +} + + +static int dirax_readable(struct image *image, struct dirax_data *dirax) +{ + int rval; + int no_string = 0; + + rval = read(dirax->pty, dirax->rbuffer+dirax->rbufpos, + dirax->rbuflen-dirax->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + dirax->rbufpos += rval; + assert(dirax->rbufpos <= dirax->rbuflen); + + while ( (!no_string) && (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; irbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + /* Is there a prompt in the buffer? */ + if ( (i+7 <= dirax->rbufpos) + && (!strncmp(dirax->rbuffer+i, "Dirax> ", 7)) ) { + block_ready = 1; + type = DIRAX_INPUT_PROMPT; + break; + } + + /* How about an ACL prompt? */ + if ( (i+10 <= dirax->rbufpos) + && (!strncmp(dirax->rbuffer+i, "acl/auto [", 10)) ) { + block_ready = 1; + type = DIRAX_INPUT_ACL; + break; + } + + if ( (dirax->rbuffer[i] == '\r') + && (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, dirax->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + dirax_parseline(block_buffer, image, dirax); + free(block_buffer); + endbit_length = i+2; + break; + + case DIRAX_INPUT_PROMPT : + + dirax_send_next(image, dirax); + endbit_length = i+7; + break; + + case DIRAX_INPUT_ACL : + + dirax_send_next(image,dirax ); + 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(dirax->rbuffer, + dirax->rbuffer + endbit_length, + dirax->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + dirax->rbufpos = dirax->rbufpos + - endbit_length; + new_rbuflen = dirax->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + dirax->rbuffer = realloc(dirax->rbuffer, + new_rbuflen); + dirax->rbuflen = new_rbuflen; + + } else { + + if ( dirax->rbufpos==dirax->rbuflen ) { + + /* More buffer space is needed */ + dirax->rbuffer = realloc( + dirax->rbuffer, + dirax->rbuflen + 256); + dirax->rbuflen = dirax->rbuflen + + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + +static void write_drx(struct image *image) +{ + FILE *fh; + int i; + char filename[1024]; + + snprintf(filename, 1023, "xfel-%i.drx", 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_dirax(struct image *image) +{ + unsigned int opts; + int status; + int rval; + struct dirax_data *dirax; + + write_drx(image); + + dirax = malloc(sizeof(struct dirax_data)); + if ( dirax == NULL ) { + ERROR("Couldn't allocate memory for DirAx data.\n"); + return; + } + + dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL); + if ( dirax->pid == -1 ) { + ERROR("Failed to fork for DirAx\n"); + return; + } + if ( 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); + + } + + dirax->rbuffer = malloc(256); + dirax->rbuflen = 256; + dirax->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(dirax->pty, F_GETFL); + fcntl(dirax->pty, F_SETFL, opts | O_NONBLOCK); + + dirax->step = 1; /* This starts the "initialisation" procedure */ + dirax->finished_ok = 0; + dirax->read_cell = 0; + dirax->n_acls_tried = 0; + dirax->best_acl_nh = 0; + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(dirax->pty, &fds); + + tv.tv_sec = 10; + tv.tv_usec = 0; + + sval = select(dirax->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + int err = errno; + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + } else if ( sval != 0 ) { + rval = dirax_readable(image, dirax); + } else { + ERROR("No response from DirAx..\n"); + rval = 1; + } + + } while ( !rval ); + + close(dirax->pty); + free(dirax->rbuffer); + waitpid(dirax->pid, &status, 0); + + if ( dirax->finished_ok == 0 ) { + ERROR("DirAx doesn't seem to be working properly.\n"); + } + + free(dirax); +} diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h new file mode 100644 index 00000000..8c429710 --- /dev/null +++ b/libcrystfel/src/dirax.h @@ -0,0 +1,26 @@ +/* + * dirax.h + * + * Invoke the DirAx auto-indexing program + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef DIRAX_H +#define DIRAX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "utils.h" + + +extern void run_dirax(struct image *image); + + +#endif /* DIRAX_H */ diff --git a/libcrystfel/src/index-priv.h b/libcrystfel/src/index-priv.h new file mode 100644 index 00000000..3d8c8a22 --- /dev/null +++ b/libcrystfel/src/index-priv.h @@ -0,0 +1,29 @@ +/* + * index-priv.h + * + * Indexing private data + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef INDEXPRIV_H +#define INDEXPRIV_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include "index.h" + +struct _indexingprivate +{ + IndexingMethod indm; +}; + + +#endif /* INDEXPRIV_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c new file mode 100644 index 00000000..d5e76c50 --- /dev/null +++ b/libcrystfel/src/index.c @@ -0,0 +1,274 @@ +/* + * index.c + * + * Perform indexing (somehow) + * + * (c) 2006-2011 Thomas White + * (c) 2010 Richard Kirian + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "peaks.h" +#include "dirax.h" +#include "mosflm.h" +#include "detector.h" +#include "index.h" +#include "index-priv.h" +#include "reax.h" +#include "geometry.h" + + +/* Base class constructor for unspecialised indexing private data */ +IndexingPrivate *indexing_private(IndexingMethod indm) +{ + struct _indexingprivate *priv; + priv = calloc(1, sizeof(struct _indexingprivate)); + priv->indm = indm; + return priv; +} + + +static const char *maybes(int n) +{ + if ( n == 1 ) return ""; + return "s"; +} + + +IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + double nominal_photon_energy) +{ + int n; + int nm = 0; + IndexingPrivate **iprivs; + + while ( indm[nm] != INDEXING_NONE ) nm++; + STATUS("Preparing %i indexing method%s.\n", nm, maybes(nm)); + iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); + + for ( n=0; nindm ) { + case INDEXING_NONE : + free(priv[n]); + break; + case INDEXING_DIRAX : + free(priv[n]); + break; + case INDEXING_MOSFLM : + free(priv[n]); + break; + case INDEXING_REAX : + reax_cleanup(priv[n]); + break; + } + + n++; + + } +} + + +void map_all_peaks(struct image *image) +{ + int i, n; + + n = image_feature_count(image->features); + + /* Map positions to 3D */ + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); + f->rx = r.u; f->ry = r.v; f->rz = r.w; + + } +} + + +void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, + int cellr, int verbose, IndexingPrivate **ipriv, + int config_insane) +{ + int i; + int n = 0; + + if ( indm == NULL ) return; + + map_all_peaks(image); + image->indexed_cell = NULL; + + while ( indm[n] != INDEXING_NONE ) { + + image->ncells = 0; + + /* Index as appropriate */ + switch ( indm[n] ) { + case INDEXING_NONE : + return; + case INDEXING_DIRAX : + run_dirax(image); + break; + case INDEXING_MOSFLM : + run_mosflm(image, cell); + break; + case INDEXING_REAX : + reax_index(ipriv[n], image, cell); + break; + } + if ( image->ncells == 0 ) { + n++; + continue; + } + + for ( i=0; incells; i++ ) { + + UnitCell *new_cell = NULL; + UnitCell *cand = image->candidate_cells[i]; + + if ( verbose ) { + STATUS("--------------------\n"); + STATUS("Candidate cell %i (before matching):\n", + i); + cell_print(image->candidate_cells[i]); + STATUS("--------------------\n"); + } + + /* Match or reduce the cell as appropriate */ + switch ( cellr ) { + case CELLR_NONE : + new_cell = cell_new_from_cell(cand); + break; + case CELLR_REDUCE : + new_cell = match_cell(cand, cell, verbose, 1); + break; + case CELLR_COMPARE : + new_cell = match_cell(cand, cell, verbose, 0); + break; + case CELLR_COMPARE_AB : + new_cell = match_cell_ab(cand, cell); + break; + } + + /* No cell? Move on to the next candidate */ + if ( new_cell == NULL ) continue; + + /* Sanity check */ + image->reflections = find_intersections(image, new_cell); + image->indexed_cell = new_cell; + + if ( !config_insane && !peak_sanity_check(image) ) { + cell_free(new_cell); + image->indexed_cell = NULL; + continue; + } + + goto done; /* Success */ + + } + + for ( i=0; incells; i++ ) { + cell_free(image->candidate_cells[i]); + image->candidate_cells[i] = NULL; + } + + /* Move on to the next indexing method */ + n++; + + } + +done: + for ( i=0; incells; i++ ) { + /* May free(NULL) if all algorithms were tried and no success */ + cell_free(image->candidate_cells[i]); + } +} + + +IndexingMethod *build_indexer_list(const char *str, int *need_cell) +{ + int n, i; + char **methods; + IndexingMethod *list; + int tmp; + + if ( need_cell == NULL ) need_cell = &tmp; + *need_cell = 0; + + n = assplode(str, ",", &methods, ASSPLODE_NONE); + list = malloc((n+1)*sizeof(IndexingMethod)); + + for ( i=0; i + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef INDEX_H +#define INDEX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include "cell.h" +#include "image.h" +#include "detector.h" + + +/* Indexing methods */ +typedef enum { + INDEXING_NONE, + INDEXING_DIRAX, + INDEXING_MOSFLM, + INDEXING_REAX, +} IndexingMethod; + + +/* Cell reduction methods */ +enum { + CELLR_NONE, + CELLR_REDUCE, + CELLR_COMPARE, + CELLR_COMPARE_AB, +}; + + +typedef struct _indexingprivate IndexingPrivate; + +extern IndexingPrivate *indexing_private(IndexingMethod indm); + +extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, + const char *filename, + struct detector *det, + double nominal_photon_energy); + +extern void map_all_peaks(struct image *image); + +extern void index_pattern(struct image *image, UnitCell *cell, + IndexingMethod *indm, int cellr, int verbose, + IndexingPrivate **priv, int config_insane); + +extern void cleanup_indexing(IndexingPrivate **priv); + +extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); + +#endif /* INDEX_H */ diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c new file mode 100644 index 00000000..5df5af21 --- /dev/null +++ b/libcrystfel/src/mosflm.c @@ -0,0 +1,609 @@ +/* + * mosflm.c + * + * Invoke the DPS auto-indexing algorithm through MOSFLM + * + * (c) 2010 Richard Kirian + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +/* TODO: + * + * Properly read the newmat file (don't use fscanf-- spaces between numers + * are not guaranteed) + * + * "Success" is indicated by existence of NEWMAT file written by mosflm. + * Better to interact with mosflm directly in order to somehow verify success. + * + * Investigate how these keywords affect mosflms behavior: + * + * MOSAICITY + * DISPERSION + * DIVERGENCE + * POLARISATION + * POSTREF BEAM + * POSTREF USEBEAM OFF + * PREREFINE ON + * EXTRA ON + * POSTREF ON + * + * These did not seem to affect the results by my (Rick's) experience, probably + * because they are only used conjunction with image intensity data, but it's + * worth another look at the documentation. + */ + +#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 "mosflm.h" +#include "utils.h" +#include "peaks.h" + + +#define MOSFLM_VERBOSE 0 + + +typedef enum { + MOSFLM_INPUT_NONE, + MOSFLM_INPUT_LINE, + MOSFLM_INPUT_PROMPT +} MOSFLMInputType; + + +struct mosflm_data { + + /* MOSFLM auto-indexing low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + + /* MOSFLM high-level stuff */ + char newmatfile[128]; + char imagefile[128]; + char sptfile[128]; + int step; + int finished_ok; + UnitCell *target_cell; + +}; + + +static void mosflm_parseline(const char *line, struct image *image, + struct mosflm_data *dirax) +{ + #if MOSFLM_VERBOSE + char *copy; + int i; + + copy = strdup(line); + for ( i=0; ilambda; + + image->candidate_cells[0] = cell_new(); + + cell_set_reciprocal(image->candidate_cells[0], + asz*c, asy*c, asx*c, + bsz*c, bsy*c, bsx*c, + csz*c, csy*c, csx*c); + + image->ncells = 1; + + return 0; +} + + +/* Need to sort mosflm peaks by intensity... */ +struct sptline { + double x; /* x coordinate of peak */ + double y; /* y coordinate of peak */ + double h; /* height of peak */ + double s; /* sigma of peak */ +}; + + +static int compare_vals(const void *ap, const void *bp) +{ + const struct sptline a = *(struct sptline *)ap; + const struct sptline b = *(struct sptline *)bp; + + if ( a.h < b.h ) return 1; + if ( a.h > b.h ) return -1; + return 0; +} + + +/* Write .spt file for mosflm */ +static void write_spt(struct image *image, const char *filename) +{ + FILE *fh; + int i; + double fclen = 67.8; /* fake camera length in mm */ + double fpix = 0.075; /* fake pixel size in mm */ + double pix; + double height = 100.0; + double sigma = 1.0; + int nPeaks = image_feature_count(image->features); + struct sptline *sptlines; + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + return; + } + + fprintf(fh, "%10d %10d %10.8f %10.6f %10.6f\n", 1, 1, fpix, 1.0, 0.0); + fprintf(fh, "%10d %10d\n", 1, 1); + fprintf(fh, "%10.5f %10.5f\n", 0.0, 0.0); + + sptlines = malloc(sizeof(struct sptline)*nPeaks); + + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + p = find_panel(image->det, f->fs, f->ss); + if ( p == NULL ) continue; + + pix = 1000.0/p->res; /* pixel size in mm */ + height = f->intensity; + + xs = (f->fs-p->min_fs)*p->fsx + (f->ss-p->min_ss)*p->ssx; + ys = (f->ss-p->min_fs)*p->fsy + (f->ss-p->min_ss)*p->ssy; + rx = xs + p->cnx; + ry = ys + p->cny; + + sptlines[i].x = ry*pix*fclen/p->clen/1000.0; + sptlines[i].y = -rx*pix*fclen/p->clen/1000.0; + sptlines[i].h = height; + sptlines[i].s = sigma/1000.0; + + } + + qsort(sptlines, nPeaks, sizeof(struct sptline), compare_vals); + + for ( i=0; ipty, line, strlen(line)); +} + + +static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) +{ + char tmp[256]; + char symm[32]; + const char *sg; + double wavelength; + double a, b, c, alpha, beta, gamma; + int i, j; + + switch ( mosflm->step ) { + + case 1 : + mosflm_sendline("DETECTOR ROTATION HORIZONTAL" + " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" + " RECTANGULAR\n", mosflm); + break; + + case 2 : + if ( mosflm->target_cell != NULL ) { + cell_get_parameters(mosflm->target_cell, &a, &b, &c, + &alpha, &beta, &gamma); + snprintf(tmp, 255, + "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", + a*1e10, b*1e10, c*1e10, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + mosflm_sendline(tmp, mosflm); + } else { + mosflm_sendline("# Do nothing\n", mosflm); + } + break; + + case 3 : + if ( mosflm->target_cell != NULL ) { + sg = cell_get_spacegroup(mosflm->target_cell); + /* Remove white space from space group */ + j = 0; + for ( i=0; ilambda*1e10; + snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength); + mosflm_sendline(tmp, mosflm); + break; + + case 7 : + snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile); + mosflm_sendline(tmp, mosflm); + break; + + case 8 : + snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile); + mosflm_sendline(tmp, mosflm); + break; + + case 9 : + snprintf(tmp, 255, "AUTOINDEX DPS FILE %s" + " IMAGE 1 MAXCELL 1000 REFINE\n", + mosflm->sptfile); + + /* "This option is only available if you e-mail Andrew Leslie + * and ask for it." - MOSFLM + * snprintf(tmp, 255, "AUTOINDEX NODISPLAY IMAGE 1 FILE %s\n", + * mosflm->sptfile); */ + mosflm_sendline(tmp, mosflm); + break; + + case 10 : + mosflm_sendline("GO\n", mosflm); + mosflm->finished_ok = 1; + break; + + default: + mosflm_sendline("exit\n", mosflm); + return; + + } + + mosflm->step++; +} + + +static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) +{ + int rval; + int no_string = 0; + + rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos, + mosflm->rbuflen-mosflm->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + mosflm->rbufpos += rval; + assert(mosflm->rbufpos <= mosflm->rbuflen); + + while ( (!no_string) && (mosflm->rbufpos > 0) ) { + + int i; + int block_ready = 0; + MOSFLMInputType type = MOSFLM_INPUT_NONE; + + /* See if there's a full line in the buffer yet */ + for ( i=0; irbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + /* Is there a prompt in the buffer? */ + if ( (i+10 <= mosflm->rbufpos) + && (!strncmp(mosflm->rbuffer+i, "MOSFLM => ", 10)) ) { + block_ready = 1; + type = MOSFLM_INPUT_PROMPT; + break; + } + + if ( (mosflm->rbuffer[i] == '\r') + && (mosflm->rbuffer[i+1] == '\n') ) { + block_ready = 1; + type = MOSFLM_INPUT_LINE; + break; + } + + } + + if ( block_ready ) { + + unsigned int new_rbuflen; + unsigned int endbit_length; + char *block_buffer = NULL; + + switch ( type ) { + + case MOSFLM_INPUT_LINE : + + block_buffer = malloc(i+1); + memcpy(block_buffer, mosflm->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + mosflm_parseline(block_buffer, image, mosflm); + free(block_buffer); + endbit_length = i+2; + break; + + case MOSFLM_INPUT_PROMPT : + + mosflm_send_next(image, mosflm); + endbit_length = i+7; + break; + + default : + + /* Obviously, this never happens :) */ + ERROR("Unrecognised MOSFLM input mode! " + "I don't know how to understand MOSFLM\n"); + return 1; + + } + + /* Now the block's been parsed, it should be + * forgotten about */ + memmove(mosflm->rbuffer, + mosflm->rbuffer + endbit_length, + mosflm->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + mosflm->rbufpos = mosflm->rbufpos + - endbit_length; + new_rbuflen = mosflm->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + mosflm->rbuffer = realloc(mosflm->rbuffer, + new_rbuflen); + mosflm->rbuflen = new_rbuflen; + + } else { + + if ( mosflm->rbufpos==mosflm->rbuflen ) { + + /* More buffer space is needed */ + mosflm->rbuffer = realloc( + mosflm->rbuffer, + mosflm->rbuflen + 256); + mosflm->rbuflen = mosflm->rbuflen + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + +void run_mosflm(struct image *image, UnitCell *cell) +{ + struct mosflm_data *mosflm; + unsigned int opts; + int status; + int rval; + + mosflm = malloc(sizeof(struct mosflm_data)); + if ( mosflm == NULL ) { + ERROR("Couldn't allocate memory for MOSFLM data.\n"); + return; + } + + mosflm->target_cell = cell; + + snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id); + write_img(image, mosflm->imagefile); /* Dummy image */ + + snprintf(mosflm->sptfile, 127, "xfel-%i_001.spt", image->id); + write_spt(image, mosflm->sptfile); + + snprintf(mosflm->newmatfile, 127, "xfel-%i.newmat", image->id); + remove(mosflm->newmatfile); + + mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); + if ( mosflm->pid == -1 ) { + ERROR("Failed to fork for MOSFLM\n"); + free(mosflm); + return; + } + if ( mosflm->pid == 0 ) { + + /* Child process: invoke MOSFLM */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + execlp("ipmosflm", "", (char *)NULL); + ERROR("Failed to invoke MOSFLM.\n"); + _exit(0); + + } + + mosflm->rbuffer = malloc(256); + mosflm->rbuflen = 256; + mosflm->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(mosflm->pty, F_GETFL); + fcntl(mosflm->pty, F_SETFL, opts | O_NONBLOCK); + + mosflm->step = 1; /* This starts the "initialisation" procedure */ + mosflm->finished_ok = 0; + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(mosflm->pty, &fds); + + tv.tv_sec = 10; + tv.tv_usec = 0; + + sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + int err = errno; + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + } else if ( sval != 0 ) { + rval = mosflm_readable(image, mosflm); + } else { + ERROR("No response from MOSFLM..\n"); + rval = 1; + } + + } while ( !rval ); + + close(mosflm->pty); + free(mosflm->rbuffer); + waitpid(mosflm->pid, &status, 0); + + if ( mosflm->finished_ok == 0 ) { + ERROR("MOSFLM doesn't seem to be working properly.\n"); + } else { + /* Read the mosflm NEWMAT file and get cell if found */ + read_newmat(mosflm->newmatfile, image); + } + + free(mosflm); +} diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h new file mode 100644 index 00000000..82b80a44 --- /dev/null +++ b/libcrystfel/src/mosflm.h @@ -0,0 +1,27 @@ +/* + * mosflm.h + * + * Invoke the DPS auto-indexing algorithm through MOSFLM + * + * (c) 2010 Richard Kirian + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef MOSFLM_H +#define MOSFLM_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "utils.h" + + +extern void run_mosflm(struct image *image, UnitCell *cell); + + +#endif /* MOSFLM_H */ diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c new file mode 100644 index 00000000..ff1ea582 --- /dev/null +++ b/libcrystfel/src/reax.c @@ -0,0 +1,728 @@ +/* + * reax.c + * + * A new auto-indexer + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "peaks.h" +#include "cell.h" +#include "index.h" +#include "index-priv.h" + + +struct dvec +{ + double x; + double y; + double z; + double th; + double ph; +}; + + +struct reax_private +{ + IndexingPrivate base; + struct dvec *directions; + int n_dir; + double angular_inc; + + double *fft_in; + fftw_complex *fft_out; + fftw_plan plan; + int nel; + + fftw_complex *r_fft_in; + fftw_complex *r_fft_out; + fftw_plan r_plan; + int ch; + int cw; +}; + + +static double check_dir(struct dvec *dir, ImageFeatureList *flist, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, + const char *rg, struct detector *det) +{ + int n, i; + double tot; + + for ( i=0; ifs, f->ss); + assert(p != NULL); + + if ( p->rigid_group != rg ) continue; + + } + + val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z; + + idx = nel/2 + nel*val/(2.0*pmax); + fft_in[idx]++; + + } + + fftw_execute_dft_r2c(plan, fft_in, fft_out); + + tot = 0.0; + for ( i=smin; i<=smax; i++ ) { + double re, im; + re = fft_out[i][0]; + im = fft_out[i][1]; + tot += sqrt(re*re + im*im); + } + + return tot; +} + + +/* Refine a direct space vector. From Clegg (1984) */ +static double iterate_refine_vector(double *x, double *y, double *z, + ImageFeatureList *flist) +{ + int fi, n, err; + gsl_matrix *C; + gsl_vector *A; + gsl_vector *t; + gsl_matrix *s_vec; + gsl_vector *s_val; + double dtmax; + + A = gsl_vector_calloc(3); + C = gsl_matrix_calloc(3, 3); + + n = image_feature_count(flist); + fesetround(1); + for ( fi=0; firx*(*x) + f->ry*(*y) + f->rz*(*z); /* Sorry ... */ + kn = nearbyint(kno); + if ( kn - kno > 0.3 ) continue; + + xv[0] = f->rx; xv[1] = f->ry; xv[2] = f->rz; + + for ( i=0; i<3; i++ ) { + + val = gsl_vector_get(A, i); + gsl_vector_set(A, i, val+xv[i]*kn); + + for ( j=0; j<3; j++ ) { + val = gsl_matrix_get(C, i, j); + gsl_matrix_set(C, i, j, val+xv[i]*xv[j]); + } + + } + + } + + s_val = gsl_vector_calloc(3); + s_vec = gsl_matrix_calloc(3, 3); + err = gsl_linalg_SV_decomp_jacobi(C, s_vec, s_val); + if ( err ) { + ERROR("SVD failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + gsl_matrix_free(C); + gsl_vector_free(A); + return 0.0; + } + + t = gsl_vector_calloc(3); + err = gsl_linalg_SV_solve(C, s_vec, s_val, A, t); + if ( err ) { + ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + gsl_matrix_free(C); + gsl_vector_free(A); + gsl_vector_free(t); + return 0.0; + } + + gsl_matrix_free(s_vec); + gsl_vector_free(s_val); + + dtmax = fabs(*x - gsl_vector_get(t, 0)); + dtmax += fabs(*y - gsl_vector_get(t, 1)); + dtmax += fabs(*z - gsl_vector_get(t, 2)); + + *x = gsl_vector_get(t, 0); + *y = gsl_vector_get(t, 1); + *z = gsl_vector_get(t, 2); + + gsl_matrix_free(C); + gsl_vector_free(A); + + return dtmax; +} + + +static void refine_cell(struct image *image, UnitCell *cell, + ImageFeatureList *flist) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double sm; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + i = 0; + do { + + sm = iterate_refine_vector(&ax, &ay, &az, flist); + sm += iterate_refine_vector(&bx, &by, &bz, flist); + sm += iterate_refine_vector(&cx, &cy, &cz, flist); + i++; + + } while ( (sm > 0.001e-9) && (i<10) ); + + cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); + + if ( i == 10 ) { + cell_free(image->indexed_cell); + image->indexed_cell = NULL; + } +} + + +static void fine_search(struct reax_private *p, ImageFeatureList *flist, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, double th_cen, double ph_cen, + double *x, double *y, double *z) +{ + double fom = 0.0; + double th, ph; + double inc; + struct dvec dir; + int i, s; + double max, modv; + + inc = p->angular_inc / 100.0; + + for ( th=th_cen-p->angular_inc; th<=th_cen+p->angular_inc; th+=inc ) { + for ( ph=ph_cen-p->angular_inc; ph<=ph_cen+p->angular_inc; ph+=inc ) { + + double new_fom; + + dir.x = cos(ph) * sin(th); + dir.y = sin(ph) * sin(th); + dir.z = cos(th); + + new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out, + plan, smin, smax, NULL, NULL); + + if ( new_fom > fom ) { + fom = new_fom; + *x = dir.x; *y = dir.y; *z = dir.z; + } + + } + } + + s = -1; + max = 0.0; + for ( i=smin; i<=smax; i++ ) { + + double re, im, m; + + re = fft_out[i][0]; + im = fft_out[i][1]; + m = sqrt(re*re + im*im); + if ( m > max ) { + max = m; + s = i; + } + + } + assert(s>0); + + modv = (double)s / (2.0*pmax); + *x *= modv; *y *= modv; *z *= modv; +} + + +static double get_model_phase(double x, double y, double z, ImageFeatureList *f, + int nel, double pmax, double *fft_in, + fftw_complex *fft_out, fftw_plan plan, + int smin, int smax, const char *rg, + struct detector *det) +{ + struct dvec dir; + int s, i; + double max; + double re, im; + + dir.x = x; dir.y = y; dir.z = z; + + check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax, + rg, det); + + s = -1; + max = 0.0; + for ( i=smin; i<=smax; i++ ) { + + double re, im, m; + + re = fft_out[i][0]; + im = fft_out[i][1]; + m = sqrt(re*re + im*im); + if ( m > max ) { + max = m; + s = i; + } + + } + + re = fft_out[s][0]; + im = fft_out[s][1]; + + return atan2(im, re); +} + + +static void refine_rigid_group(struct image *image, UnitCell *cell, + const char *rg, int nel, double pmax, + double *fft_in, fftw_complex *fft_out, + fftw_plan plan, int smin, int smax, + struct detector *det, struct reax_private *pr) +{ + double ax, ay, az, ma; + double bx, by, bz, mb; + double cx, cy, cz, mc; + double pha, phb, phc; + struct panel *p; + int i, j; + fftw_complex *r_fft_in; + fftw_complex *r_fft_out; + double m2m; + signed int aix, aiy; + signed int bix, biy; + signed int cix, ciy; + double max; + int max_i, max_j; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + ma = modulus(ax, ay, az); + mb = modulus(bx, by, bz); + mc = modulus(cx, cy, cz); + + pha = get_model_phase(ax/ma, ay/ma, az/ma, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + phb = get_model_phase(bx/mb, by/mb, bz/mb, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + phc = get_model_phase(cx/mc, cy/mc, cz/mc, image->features, nel, pmax, + fft_in, fft_out, plan, smin, smax, rg, det); + + for ( i=0; in_panels; i++ ) { + if ( det->panels[i].rigid_group == rg ) { + p = &det->panels[i]; + break; + } + } + + r_fft_in = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); + r_fft_out = fftw_malloc(pr->cw*pr->ch*sizeof(fftw_complex)); + for ( i=0; icw; i++ ) { + for ( j=0; jch; j++ ) { + r_fft_in[i+pr->cw*j][0] = 0.0; + r_fft_in[i+pr->cw*j][1] = 0.0; + } + } + + ma = modulus(ax, ay, 0.0); + mb = modulus(bx, by, 0.0); + mc = modulus(cx, cy, 0.0); + m2m = ma; + if ( mb > m2m ) m2m = mb; + if ( mc > m2m ) m2m = mc; + + aix = (pr->cw/2)*ax/m2m; aiy = (pr->ch/2)*ay/m2m; + bix = (pr->cw/2)*bx/m2m; biy = (pr->ch/2)*by/m2m; + cix = (pr->cw/2)*cx/m2m; ciy = (pr->ch/2)*cy/m2m; + + if ( aix < 0 ) aix += pr->cw/2; + if ( bix < 0 ) bix += pr->cw/2; + if ( cix < 0 ) cix += pr->cw/2; + + if ( aiy < 0 ) aiy += pr->ch/2; + if ( biy < 0 ) biy += pr->ch/2; + if ( ciy < 0 ) ciy += pr->ch/2; + + r_fft_in[aix + pr->cw*aiy][0] = cos(pha); + r_fft_in[aix + pr->cw*aiy][1] = sin(pha); + r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][0] = cos(pha); + r_fft_in[pr->cw-aix + pr->cw*(pr->ch-aiy)][1] = -sin(pha); + + r_fft_in[bix + pr->cw*biy][0] = cos(phb); + r_fft_in[bix + pr->cw*biy][1] = sin(phb); + r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][0] = cos(phb); + r_fft_in[pr->cw-bix + pr->cw*(pr->ch-biy)][1] = -sin(phb); + + r_fft_in[cix + pr->cw*ciy][0] = cos(phc); + r_fft_in[cix + pr->cw*ciy][1] = sin(phc); + r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][0] = cos(phc); + r_fft_in[pr->cw-cix + pr->cw*(pr->ch-ciy)][1] = -sin(phc); + + const int tidx = 1; + r_fft_in[tidx][0] = 1.0; + r_fft_in[tidx][1] = 0.0; + +// STATUS("%i %i\n", aix, aiy); +// STATUS("%i %i\n", bix, biy); +// STATUS("%i %i\n", cix, ciy); + + fftw_execute_dft(pr->r_plan, r_fft_in, r_fft_out); + +// max = 0.0; +// FILE *fh = fopen("centering.dat", "w"); +// for ( i=0; icw; i++ ) { +// for ( j=0; jch; j++ ) { +// +// double re, im, am, ph; +// +// re = r_fft_out[i + pr->cw*j][0]; +// im = r_fft_out[i + pr->cw*j][1]; +// am = sqrt(re*re + im*im); +// ph = atan2(im, re); +// +// if ( am > max ) { +// max = am; +// max_i = i; +// max_j = j; +// } +// +// fprintf(fh, "%f ", am); +// +// } +// fprintf(fh, "\n"); +// } +// STATUS("Max at %i, %i\n", max_i, max_j); +// fclose(fh); +// exit(1); + +// STATUS("Offsets for '%s': %.2f, %.2f pixels\n", rg, dx, dy); +} + + +static void refine_all_rigid_groups(struct image *image, UnitCell *cell, + int nel, double pmax, + double *fft_in, fftw_complex *fft_out, + fftw_plan plan, int smin, int smax, + struct detector *det, + struct reax_private *p) +{ + int i; + + for ( i=0; idet->num_rigid_groups; i++ ) { + refine_rigid_group(image, cell, image->det->rigid_groups[i], + nel, pmax, fft_in, fft_out, plan, smin, smax, + det, p); + } +} + + +void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) +{ + int i; + struct reax_private *p; + double fom; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double mod_a, mod_b, mod_c; + double al, be, ga; + double th, ph; + double *fft_in; + fftw_complex *fft_out; + int smin, smax; + double amin, amax; + double bmin, bmax; + double cmin, cmax; + double pmax; + int n; + const double ltol = 5.0; /* Direct space axis length + * tolerance in percent */ + const double angtol = deg2rad(1.5); /* Direct space angle tolerance + * in radians */ + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + fft_in = fftw_malloc(p->nel*sizeof(double)); + fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + mod_a = modulus(ax, ay, az); + amin = mod_a * (1.0-ltol/100.0); + amax = mod_a * (1.0+ltol/100.0); + + mod_b = modulus(bx, by, bz); + bmin = mod_b * (1.0-ltol/100.0); + bmax = mod_b * (1.0+ltol/100.0); + + mod_c = modulus(cx, cy, cz); + cmin = mod_c * (1.0-ltol/100.0); + cmax = mod_c * (1.0+ltol/100.0); + + al = angle_between(bx, by, bz, cx, cy, cz); + be = angle_between(ax, ay, az, cx, cy, cz); + ga = angle_between(ax, ay, az, bx, by, bz); + + pmax = 0.0; + n = image_feature_count(image->features); + for ( i=0; ifeatures, i); + if ( f == NULL ) continue; + + val = modulus(f->rx, f->ry, f->rz); + if ( val > pmax ) pmax = val; + + } + + /* Sanity check */ + if ( pmax < 1e4 ) return; + + /* Search for a */ + smin = 2.0*pmax * amin; + smax = 2.0*pmax * amax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &ax, &ay, &az); + + /* Search for b */ + smin = 2.0*pmax * bmin; + smax = 2.0*pmax * bmax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom, ang; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, ax, ay, az); + if ( fabs(ang-ga) > angtol ) continue; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &bx, &by, &bz); + + /* Search for c */ + smin = 2.0*pmax * cmin; + smax = 2.0*pmax * cmax; + fom = 0.0; th = 0.0; ph = 0.0; + for ( i=0; in_dir; i++ ) { + + double new_fom, ang; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, ax, ay, az); + if ( fabs(ang-be) > angtol ) continue; + + ang = angle_between(p->directions[i].x, p->directions[i].y, + p->directions[i].z, bx, by, bz); + if ( fabs(ang-al) > angtol ) continue; + + new_fom = check_dir(&p->directions[i], image->features, + p->nel, pmax, fft_in, fft_out, p->plan, + smin, smax, NULL, NULL); + if ( new_fom > fom ) { + fom = new_fom; + th = p->directions[i].th; + ph = p->directions[i].ph; + } + + } + fine_search(p, image->features, p->nel, pmax, fft_in, fft_out, + p->plan, smin, smax, th, ph, &cx, &cy, &cz); + + image->candidate_cells[0] = cell_new(); + cell_set_cartesian(image->candidate_cells[0], + ax, ay, az, bx, by, bz, cx, cy, cz); + + refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax, + fft_in, fft_out, p->plan, smin, smax, + image->det, p); + map_all_peaks(image); + refine_cell(image, image->candidate_cells[0], image->features); + + fftw_free(fft_in); + fftw_free(fft_out); + + image->ncells = 1; +} + + +IndexingPrivate *reax_prepare() +{ + struct reax_private *p; + int samp; + double th; + + p = calloc(1, sizeof(*p)); + if ( p == NULL ) return NULL; + + p->base.indm = INDEXING_REAX; + + p->angular_inc = deg2rad(1.7); /* From Steller (1997) */ + + /* Reserve memory, over-estimating the number of directions */ + samp = 2.0*M_PI / p->angular_inc; + p->directions = malloc(samp*samp*sizeof(struct dvec)); + if ( p == NULL) { + free(p); + return NULL; + } + STATUS("Allocated space for %i directions\n", samp*samp); + + /* Generate vectors for 1D Fourier transforms */ + fesetround(1); /* Round to nearest */ + p->n_dir = 0; + for ( th=0.0; thangular_inc ) { + + double ph, phstep, n_phstep; + + n_phstep = 2.0*M_PI*sin(th)/p->angular_inc; + n_phstep = nearbyint(n_phstep); + phstep = 2.0*M_PI/n_phstep; + + for ( ph=0.0; ph<2.0*M_PI; ph+=phstep ) { + + struct dvec *dir; + + assert(p->n_dirdirections[p->n_dir++]; + + dir->x = cos(ph) * sin(th); + dir->y = sin(ph) * sin(th); + dir->z = cos(th); + dir->th = th; + dir->ph = ph; + + } + + } + STATUS("Generated %i directions (angular increment %.3f deg)\n", + p->n_dir, rad2deg(p->angular_inc)); + + p->nel = 1024; + + /* These arrays are not actually used */ + p->fft_in = fftw_malloc(p->nel*sizeof(double)); + p->fft_out = fftw_malloc((p->nel/2 + 1)*sizeof(fftw_complex)); + + p->plan = fftw_plan_dft_r2c_1d(p->nel, p->fft_in, p->fft_out, + FFTW_MEASURE); + + p->cw = 128; p->ch = 128; + + /* Also not used */ + p->r_fft_in = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); + p->r_fft_out = fftw_malloc(p->cw*p->ch*sizeof(fftw_complex)); + + p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, + 1, FFTW_MEASURE); + + return (IndexingPrivate *)p; +} + + +void reax_cleanup(IndexingPrivate *pp) +{ + struct reax_private *p; + + assert(pp->indm == INDEXING_REAX); + p = (struct reax_private *)pp; + + fftw_destroy_plan(p->plan); + fftw_free(p->fft_in); + fftw_free(p->fft_out); + + fftw_destroy_plan(p->r_plan); + fftw_free(p->r_fft_in); + fftw_free(p->r_fft_out); + + free(p); +} diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h new file mode 100644 index 00000000..543cd0d5 --- /dev/null +++ b/libcrystfel/src/reax.h @@ -0,0 +1,27 @@ +/* + * reax.h + * + * A new auto-indexer + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef REAX_H +#define REAX_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "cell.h" + +extern IndexingPrivate *reax_prepare(void); +extern void reax_cleanup(IndexingPrivate *pp); + +extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell); + +#endif /* REAX_H */ diff --git a/libcrystfel/src/render.c b/libcrystfel/src/render.c new file mode 100644 index 00000000..dbaca2fd --- /dev/null +++ b/libcrystfel/src/render.c @@ -0,0 +1,442 @@ +/* + * render.c + * + * Render a high dynamic range buffer in some sensible way + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#ifdef HAVE_GTK +#include +#endif + +#include +#include +#include + +#ifdef HAVE_TIFF +#include +#endif + + +#include "hdf5-file.h" +#include "render.h" +#include "peaks.h" +#include "filters.h" +#include "utils.h" + + +static void render_rgb(double val, double max, + double *rp, double *gp, double *bp) +{ + int s; + double p; + double r, g, b; + + s = val / (max/6); + p = fmod(val, max/6.0); + p /= (max/6.0); + + r = 0.0; g = 0.0; b = 0.0; + + if ( (val < 0.0) ) { + s = 0; + p = 0; + } + if ( (val > max) ) { + s = 6; + } + switch ( s ) { + case 0 : { /* Black to blue */ + r = 0; g = 0; b = p; + break; + } + case 1 : { /* Blue to pink */ + r = p; g = 0; b = 1.0; + break; + } + case 2 : { /* Pink to red */ + r = 1.0; g = 0; b = (1.0-p)*1.0; + break; + } + case 3 : { /* Red to Orange */ + r = 1.0; g = 0.5*p; b = 0; + break; + } + case 4 : { /* Orange to Yellow */ + r = 1.0; g = 0.5 + 0.5*p; b = 0; + break; + } + case 5 : { /* Yellow to White */ + r = 1.0; g = 1.0; b = 1.0*p; + break; + } + case 6 : { /* Pixel has hit the maximum value */ + r = 1.0; g = 1.0; b = 1.0; + break; + } + } + + *rp = r; + *gp = g; + *bp = b; +} + + +static void render_ratio(double val, double max, + double *rp, double *gp, double *bp) +{ + if ( val <= 1.0 ) { + render_rgb(val, 2.0, rp, gp, bp); + } else { + /* Your homework is to simplify this expression */ + val = ((val-1.0)/(max-1.0)) * (max/2.0) + max/2.0; + render_rgb(val, max, rp, gp, bp); + } +} + + +static void render_mono(double val, double max, + double *rp, double *gp, double *bp) +{ + double p; + p = val / max; + if ( val < 0.0 ) p = 0.0; + if ( val > max ) p = 1.0; + *rp = p; + *gp = p; + *bp = p; +} + + +static void render_invmono(double val, double max, + double *rp, double *gp, double *bp) +{ + double p; + p = val / max; + p = 1.0 - p; + if ( val < 0.0 ) p = 1.0; + if ( val > max ) p = 0.0; + *rp = p; + *gp = p; + *bp = p; +} + + +void render_scale(double val, double max, int scale, + double *rp, double *gp, double *bp) +{ + switch ( scale ) { + case SCALE_COLOUR : + render_rgb(val, max, rp, gp, bp); + break; + case SCALE_MONO : + render_mono(val, max, rp, gp, bp); + break; + case SCALE_INVMONO : + render_invmono(val, max, rp, gp, bp); + break; + case SCALE_RATIO : + render_ratio(val, max, rp, gp, bp); + break; + } +} + + +#ifdef HAVE_GTK + +static float *get_binned_panel(struct image *image, int binning, + int min_fs, int max_fs, int min_ss, int max_ss) +{ + float *data; + int x, y; + int w, h; + int fw; + float *in; + + fw = image->width; + in = image->data; + + /* Some pixels might get discarded */ + w = (max_fs - min_fs + 1) / binning; + h = (max_ss - min_ss + 1) / binning; + + data = malloc(w*h*sizeof(float)); + + for ( x=0; xwidth*image->height; i++ ) { + if ( image->data[i] > max ) max = image->data[i]; + } + hdr = get_binned_panel(image, binning, min_fs, max_fs, min_ss, max_ss); + if ( hdr == NULL ) return NULL; + + /* Rendered (colourful) version */ + data = malloc(3*w*h); + if ( data == NULL ) { + free(hdr); + return NULL; + } + + max /= boost; + if ( max <= 6 ) { max = 10; } + /* These x,y coordinates are measured relative to the bottom-left + * corner */ + for ( y=0; ydet->n_panels; + GdkPixbuf **pixbufs; + + pixbufs = calloc(np, sizeof(GdkPixbuf*)); + if ( pixbufs == NULL ) { + *n_pixbufs = 0; + return NULL; + } + + for ( i=0; idet->panels[i].min_fs, + image->det->panels[i].max_fs, + image->det->panels[i].min_ss, + image->det->panels[i].max_ss); + + } + + *n_pixbufs = np; + + return pixbufs; +} + + +GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale) +{ + guchar *data; + size_t x, y; + int max; + + data = malloc(3*w*h); + if ( data == NULL ) return NULL; + + max = h; + + for ( y=0; ywidth); + TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); + TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); + TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP); + TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 32); + TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); + TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); + TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); + TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, + TIFFDefaultStripSize(th, image->width*4)); + + line = _TIFFmalloc(TIFFScanlineSize(th)); + for ( y=0; yheight; y++ ) { + memcpy(line, &image->data[(image->height-1-y)*image->width], + image->width*4); + TIFFWriteScanline(th, line, y, 0); + } + _TIFFfree(line); + + TIFFClose(th); + +#else + STATUS("No TIFF support.\n"); +#endif + return 0; +} + + +int render_tiff_int16(struct image *image, const char *filename, double boost) +{ +#ifdef HAVE_TIFF + TIFF *th; + int16_t *line; + int x, y; + double max; + + th = TIFFOpen(filename, "w"); + if ( th == NULL ) return 1; + + TIFFSetField(th, TIFFTAG_IMAGEWIDTH, image->width); + TIFFSetField(th, TIFFTAG_IMAGELENGTH, image->height); + TIFFSetField(th, TIFFTAG_SAMPLESPERPIXEL, 1); + TIFFSetField(th, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_INT); /* (signed) */ + TIFFSetField(th, TIFFTAG_BITSPERSAMPLE, 16); + TIFFSetField(th, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); + TIFFSetField(th, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT); + TIFFSetField(th, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); + TIFFSetField(th, TIFFTAG_ROWSPERSTRIP, + TIFFDefaultStripSize(th, image->width*4)); + + line = _TIFFmalloc(TIFFScanlineSize(th)); + max = 0.0; + for ( y=0; yheight; y++ ) { + for ( x=0;xwidth; x++ ) { + double val; + val = image->data[x+image->height*y]; + if ( val > max ) max = val; + } + } + max /= 32767.0; + + for ( y=0; yheight; y++ ) { + for ( x=0;xwidth; x++ ) { + + double val; + + val = image->data[x+(image->height-1-y)*image->width]; + val *= ((double)boost/max); + + /* Clamp to 16-bit range, + * and work round inability of most readers to deal + * with signed integers. */ + val += 1000.0; + if ( val > 32767.0 ) val = 32767.0; + if ( val < 0.0 ) val = 0.0; + + line[x] = val; + } + + TIFFWriteScanline(th, line, y, 0); + } + _TIFFfree(line); + + TIFFClose(th); +#else + STATUS("No TIFF support.\n"); +#endif + return 0; +} + +#endif /* HAVE_GTK */ diff --git a/libcrystfel/src/render.h b/libcrystfel/src/render.h new file mode 100644 index 00000000..fbfbbdb2 --- /dev/null +++ b/libcrystfel/src/render.h @@ -0,0 +1,55 @@ +/* + * render.h + * + * Render a high dynamic range buffer in some sensible way + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef RENDER_H +#define RENDER_H + + +#include + +#include "image.h" + +enum { + SCALE_COLOUR, + SCALE_MONO, + SCALE_INVMONO, + SCALE_RATIO +}; + +/* Colour scale lookup */ +extern void render_scale(double val, double max, int scale, + double *rp, double *gp, double *bp); + + +#ifdef HAVE_GTK + +#include + +extern GdkPixbuf **render_panels(struct image *image, + int binning, int scale, double boost, + int *n_pixbufs); + +extern GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale); + +extern int render_tiff_fp(struct image *image, const char *filename); +extern int render_tiff_int16(struct image *image, const char *filename, + double boost); + +#endif /* HAVE_GTK */ + + + +#endif /* RENDER_H */ -- cgit v1.2.3