/* * mosflm.c * * Invoke the DPS auto-indexing algorithm through MOSFLM * * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian * 2010-2018 Thomas White * 2014 Takanori Nakane * * 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 . * */ /* 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 #include #ifdef HAVE_FORKPTY_PTY_H #include #endif #ifdef HAVE_FORKPTY_UTIL_H #include #endif #ifdef HAVE_CLOCK_GETTIME #include #else #include #endif #include "image.h" #include "mosflm.h" #include "utils.h" #include "peaks.h" #include "cell-utils.h" #define MOSFLM_VERBOSE 0 #define FAKE_CLEN (0.1) typedef enum { MOSFLM_INPUT_NONE, MOSFLM_INPUT_LINE, MOSFLM_INPUT_PROMPT } MOSFLMInputType; struct mosflm_private { IndexingMethod indm; UnitCell *template; }; 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; int done; int success; struct mosflm_private *mp; }; static int check_cell(struct mosflm_private *mp, struct image *image, UnitCell *cell) { Crystal *cr; /* If we sent lattice information, make sure that we got back what we * asked for, not (e.g.) some "H" version of a rhombohedral R cell */ if ( mp->indm & INDEXING_USE_LATTICE_TYPE ) { LatticeType latt_m, latt_r; char cen_m, cen_r; /* What we asked for */ latt_r = cell_get_lattice_type(mp->template); cen_r = cell_get_centering(mp->template); /* What we got back */ latt_m = cell_get_lattice_type(cell); cen_m = cell_get_centering(cell); if ( latt_r != latt_m ) { ERROR("Lattice type produced by MOSFLM (%i) does not " "match what was requested (%i). " "Please report this.\n", latt_m, latt_r); return 0; } if ( (latt_m != L_MONOCLINIC) && (cen_r != cen_m) ) { ERROR("Centering produced by MOSFLM (%c) does not " "match what was requested (%c). " "Please report this.\n", cen_m, cen_r); return 0; } /* If it's monoclinic, see the warning in mosflm_prepare() */ } cr = crystal_new(); if ( cr == NULL ) { ERROR("Failed to allocate crystal.\n"); return 0; } crystal_set_cell(cr, cell); image_add_crystal(image, cr); return 1; } static void mosflm_parseline(const char *line, struct image *image, struct mosflm_data *dirax) { if ( MOSFLM_VERBOSE || (strncmp(line, "Invocation:", 11) == 0) ) { char *copy; int i; copy = strdup(line); for ( i=0; ilambda; cell = cell_new(); /* The relationship between the coordinates in the spot file and the * resulting matrix is diabolically complicated. This transformation * seems to work, but is not derived by working through all the * transformations. */ cell_set_reciprocal(cell, -asy*c, -asz*c, asx*c, -bsy*c, -bsz*c, bsx*c, -csy*c, -csz*c, csx*c); cell_set_centering(cell, cen); cell_set_lattice_type(cell, latt); cell_set_unique_axis(cell, ua); //STATUS("My cell:\n"); //cell_print(cell); if ( check_cell(mosflm->mp, image, cell) ) { mosflm->success = 1; mosflm->done = 1; } return 0; } /* Write .spt file for mosflm */ static void write_spt(struct image *image, const char *filename) { FILE *fh; int i; int n; fh = fopen(filename, "w"); if ( !fh ) { ERROR("Couldn't open temporary file '%s'\n", filename); return; } /* Number of pixels in x, number of pixels in y, pixel size (mm), * YSCALE, OMEGA */ fputs("1 1 0.0 1.0 0.0\n", fh); /* INVERTX, ISWUNG */ fputs("0 1\n", fh); /* XBEAM, YBEAM */ fputs("0.0 0.0\n", fh); n = image_feature_count(image->features); for ( i=0; ifeatures, i); if ( f == NULL ) continue; ttx = angle_between_2d(0.0, 1.0, f->rx, 1.0/image->lambda + f->rz); tty = angle_between_2d(0.0, 1.0, f->ry, 1.0/image->lambda + f->rz); if ( f->rx < 0.0 ) ttx *= -1.0; if ( f->ry < 0.0 ) tty *= -1.0; x = -tan(ttx)*FAKE_CLEN; y = tan(tty)*FAKE_CLEN; fprintf(fh, "%10.2f %10.2f 0.0 0.0 1000.0 10.0\n", x*1e3, y*1e3); } fputs("-999.0 -999.0 -999.0 -999.0 -999.0 -999.0\n", fh); fclose(fh); } /* Write a dummy 1x1 pixel image file for MOSFLM. Without post refinement, * MOSFLM will ignore this, but it must be present. */ static void write_img(struct image *image, const char *filename) { FILE *fh; unsigned short int *intimage; intimage = malloc(sizeof(unsigned short int)); intimage[0] = 1; fh = fopen(filename, "w"); if ( !fh ) { ERROR("Couldn't open temporary file '%s'\n", filename); return; } /* Write header */ fprintf(fh, "{\nHEADER_BYTES=512;\n"); fprintf(fh, "BYTE_ORDER=little_endian;\n"); fprintf(fh, "TYPE=unsigned_short;\n"); fprintf(fh, "DIM=2;\n"); fprintf(fh, "SIZE1=1;\n"); fprintf(fh, "SIZE2=1;\n"); fprintf(fh, "}\n"); /* Header padding */ while ( ftell(fh) < 512 ) fprintf(fh," "); fwrite(intimage, sizeof(unsigned short int), 1, fh); free(intimage); fclose(fh); } static void mosflm_sendline(const char *line, struct mosflm_data *mosflm) { #if MOSFLM_VERBOSE char *copy; int i; copy = strdup(line); for ( i=0; ipty, line, strlen(line)) == -1 ) { ERROR("write() to MOSFLM failed: %s\n", strerror(errno)); } } /* Turn what we know about the unit cell into something which we can give to * MOSFLM to make it give us only indexing results compatible with the cell. */ static char *spacegroup_for_lattice(UnitCell *cell) { LatticeType latt; char centering; char *g = NULL; char *result; char ua; latt = cell_get_lattice_type(cell); centering = cell_get_centering(cell); ua = cell_get_unique_axis(cell); switch ( latt ) { case L_TRICLINIC : g = "1"; break; case L_MONOCLINIC : switch ( ua ) { case 'a' : g = "211"; break; case 'b' : g = "121"; break; case 'c' : g = "112"; break; } break; case L_ORTHORHOMBIC : g = "222"; break; case L_TETRAGONAL : g = "4"; break; case L_RHOMBOHEDRAL : g = "3"; break; case L_HEXAGONAL : if ( centering != 'H' ) { g = "6"; } else { g = "3"; } break; case L_CUBIC : g = "23"; break; } assert(g != NULL); result = malloc(32); if ( result == NULL ) return NULL; snprintf(result, 31, "%c%s", centering, g); return result; } static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) { char tmp[256]; char cen; double wavelength; double a = 0, b = 0, c = 0, alpha = 0, beta = 0, gamma = 0; switch ( mosflm->step ) { case 1 : /* Backwards-compatible workaround for different Mosflm behaviour * in version 7.2.2 */ mosflm_sendline("DETECTOR ADSC\n", mosflm); break; case 2 : mosflm_sendline("DETECTOR ROTATION HORIZONTAL" " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" " RECTANGULAR\n", mosflm); break; case 3 : if ( (mosflm->mp->indm & INDEXING_USE_LATTICE_TYPE) && (mosflm->mp->template != NULL) ) { char *symm; if ( cell_get_lattice_type(mosflm->mp->template) == L_RHOMBOHEDRAL ) { mosflm_sendline("CRYSTAL R\n", mosflm); } symm = spacegroup_for_lattice(mosflm->mp->template); snprintf(tmp, 255, "SYMM %s\n", symm); //STATUS("Asking MOSFLM for '%s'\n", symm); free(symm); mosflm_sendline(tmp, mosflm); } else { mosflm_sendline("\n", mosflm); } break; case 4 : snprintf(tmp, 255, "DISTANCE %f\n", FAKE_CLEN*1e3); mosflm_sendline(tmp, mosflm); break; case 5 : mosflm_sendline("BEAM 0.0 0.0\n", mosflm); break; case 6 : wavelength = image->lambda*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 : if ( mosflm->mp->indm & INDEXING_USE_CELL_PARAMETERS ) { cell_get_parameters(mosflm->mp->template, &a, &b, &c, &alpha, &beta, &gamma); cen = cell_get_centering(mosflm->mp->template); /* Old MOSFLM simply ignores CELL and CENTERING subkeywords. * So this doesn't do any harm. * TODO: but still better to show WARNING if MOSFLM is old. */ snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " "MAXCELL 1000 REFINE " "CELL %.1f %.1f %.1f %.1f %.1f %.1f " "CENTERING %c\n", mosflm->sptfile, a*1e10, b*1e10, c*1e10, rad2deg(alpha), rad2deg(beta), rad2deg(gamma), cen); } else { snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1 " "MAXCELL 1000 REFINE\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; } int run_mosflm(struct image *image, void *ipriv) { 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 0; } snprintf(mosflm->imagefile, 127, "xfel_001.img"); write_img(image, mosflm->imagefile); /* Dummy image */ snprintf(mosflm->sptfile, 127, "xfel_001.spt"); write_spt(image, mosflm->sptfile); snprintf(mosflm->newmatfile, 127, "xfel.newmat"); remove(mosflm->newmatfile); mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); if ( mosflm->pid == -1 ) { ERROR("Failed to fork for MOSFLM: %s\n", strerror(errno)); free(mosflm); return 0; } 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("mosflm", "mosflm", (char *)NULL); execlp("ipmosflm", "ipmosflm", (char *)NULL); ERROR("Invocation: Failed to invoke MOSFLM: %s\n", strerror(errno)); _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; mosflm->mp = (struct mosflm_private *)ipriv; mosflm->done = 0; mosflm->success = 0; rval = 0; do { fd_set fds; struct timeval tv; int sval; FD_ZERO(&fds); FD_SET(mosflm->pty, &fds); tv.tv_sec = 30; tv.tv_usec = 0; sval = select(mosflm->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 = 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, mosflm->newmatfile, image); } rval = mosflm->success; free(mosflm); return rval; } void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell) { struct mosflm_private *mp; if ( mosflm_probe(cell) == NULL ) { ERROR("Mosflm does not appear to run properly.\n"); ERROR("Please check your Mosflm installation.\n"); return NULL; } /* Flags that MOSFLM knows about */ *indm &= INDEXING_METHOD_MASK | INDEXING_USE_LATTICE_TYPE | INDEXING_USE_CELL_PARAMETERS; if ( (cell_get_centering(cell) == 'I') && (cell_get_lattice_type(cell) == L_MONOCLINIC) ) { ERROR("WARNING: Mosflm always gives the monoclinic C cell in " "preference to the monoclinic I cell choice.\n"); ERROR("To get a higher indexing rate, convert your cell to the " "monoclinic C cell choice.\n"); } mp = malloc(sizeof(struct mosflm_private)); if ( mp == NULL ) return NULL; mp->template = cell; mp->indm = *indm; return (IndexingPrivate *)mp; } void mosflm_cleanup(void *pp) { struct mosflm_private *p; p = (struct mosflm_private *)pp; free(p); } static void chop_word(char *s) { int i; size_t l = strlen(s); for ( i=0; i