diff options
Diffstat (limited to 'libcrystfel/src')
| -rw-r--r-- | libcrystfel/src/dirax.c | 123 | ||||
| -rw-r--r-- | libcrystfel/src/dirax.h | 15 | ||||
| -rw-r--r-- | libcrystfel/src/image.h | 2 | ||||
| -rw-r--r-- | libcrystfel/src/index.c | 65 | ||||
| -rw-r--r-- | libcrystfel/src/index.h | 8 | ||||
| -rw-r--r-- | libcrystfel/src/mosflm.c | 103 | ||||
| -rw-r--r-- | libcrystfel/src/mosflm.h | 4 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.c | 6 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.h | 3 | ||||
| -rw-r--r-- | libcrystfel/src/reax.h | 16 |
10 files changed, 269 insertions, 76 deletions
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 160df719..0b2debda 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -53,6 +53,7 @@ #include "dirax.h" #include "utils.h" #include "peaks.h" +#include "cell-utils.h" #define DIRAX_VERBOSE 0 @@ -68,6 +69,13 @@ typedef enum { } DirAxInputType; +struct dirax_private { + IndexingMethod indm; + float *ltl; + UnitCell *template; +}; + + struct dirax_data { /* DirAx auto-indexing low-level stuff */ @@ -85,10 +93,56 @@ struct dirax_data { int best_acl_nh; int acls_tried[MAX_CELL_CANDIDATES]; int n_acls_tried; + UnitCell *cur_cell; + int done; + int success; + + struct dirax_private *dp; }; +static int check_cell(struct dirax_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 1); + if ( out == NULL ) return 0; + + } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 0); + if ( out == NULL ) return 0; + + } else { + out = cell_new_from_cell(cell); + } + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, cell); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + + static void dirax_parseline(const char *line, struct image *image, struct dirax_data *dirax) { @@ -119,7 +173,7 @@ static void dirax_parseline(const char *line, struct image *image, if ( line[i] == 'R' ) rf = 1; if ( (line[i] == 'D') && rf ) { dirax->read_cell = 1; - image->candidate_cells[image->ncells] = cell_new(); + dirax->cur_cell = cell_new(); return; } i++; @@ -127,6 +181,7 @@ static void dirax_parseline(const char *line, struct image *image, /* Parse unit cell vectors as appropriate */ if ( dirax->read_cell == 1 ) { + /* First row of unit cell values */ float ax, ay, az; int r; @@ -136,14 +191,16 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_a(image->candidate_cells[image->ncells], + cell_set_cartesian_a(dirax->cur_cell, 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; @@ -153,14 +210,16 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_b(image->candidate_cells[image->ncells], + cell_set_cartesian_b(dirax->cur_cell, 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; @@ -170,13 +229,21 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_c(image->candidate_cells[image->ncells++], + cell_set_cartesian_c(dirax->cur_cell, cx*1e-10, cy*1e-10, cz*1e-10); dirax->read_cell = 0; + + /* Finished reading a cell. Time to check it... */ + if ( check_cell(dirax->dp, image, dirax->cur_cell) ) { + dirax->done = 1; + dirax->success = 1; + } + return; + } dirax->read_cell = 0; @@ -351,8 +418,7 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax) switch ( type ) { - case DIRAX_INPUT_LINE : - + case DIRAX_INPUT_LINE : block_buffer = malloc(i+1); memcpy(block_buffer, dirax->rbuffer, i); block_buffer[i] = '\0'; @@ -366,20 +432,17 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax) endbit_length = i+2; break; - case DIRAX_INPUT_PROMPT : - + case DIRAX_INPUT_PROMPT : dirax_send_next(image, dirax); endbit_length = i+7; break; - case DIRAX_INPUT_ACL : - + case DIRAX_INPUT_ACL : dirax_send_next(image,dirax ); endbit_length = i+10; break; - default : - + default : /* Obviously, this never happens :) */ ERROR("Unrecognised DirAx input mode! " "I don't know how to understand DirAx\n"); @@ -456,7 +519,7 @@ static void write_drx(struct image *image) } -void run_dirax(struct image *image) +int run_dirax(struct image *image, IndexingPrivate *ipriv) { unsigned int opts; int status; @@ -468,13 +531,13 @@ void run_dirax(struct image *image) dirax = malloc(sizeof(struct dirax_data)); if ( dirax == NULL ) { ERROR("Couldn't allocate memory for DirAx data.\n"); - return; + return 0; } dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL); if ( dirax->pid == -1 ) { ERROR("Failed to fork for DirAx: %s\n", strerror(errno)); - return; + return 0; } if ( dirax->pid == 0 ) { @@ -505,6 +568,9 @@ void run_dirax(struct image *image) dirax->read_cell = 0; dirax->n_acls_tried = 0; dirax->best_acl_nh = 0; + dirax->done = 0; + dirax->success = 0; + dirax->dp = (struct dirax_private *)ipriv; do { @@ -543,7 +609,7 @@ void run_dirax(struct image *image) rval = 1; } - } while ( !rval ); + } while ( !rval && !dirax->success ); close(dirax->pty); free(dirax->rbuffer); @@ -553,5 +619,24 @@ void run_dirax(struct image *image) ERROR("DirAx doesn't seem to be working properly.\n"); } + rval = dirax->success; free(dirax); + return rval; +} + + +IndexingPrivate *dirax_prepare(IndexingMethod indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) +{ + struct dirax_private *dp; + + dp = malloc(sizeof(struct dirax_private)); + if ( dp == NULL ) return NULL; + + dp->ltl = ltl; + dp->template = cell; + dp->indm = indm; + + return (IndexingPrivate *)dp; } diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h index ce3cd4b0..09a794b0 100644 --- a/libcrystfel/src/dirax.h +++ b/libcrystfel/src/dirax.h @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010,2012 Thomas White <taw@physics.org> + * 2010,2012-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -33,10 +33,13 @@ #include <config.h> #endif -#include "utils.h" +#include "index.h" +extern int run_dirax(struct image *image, IndexingPrivate *ipriv); -extern void run_dirax(struct image *image); - +extern IndexingPrivate *dirax_prepare(IndexingMethod indm, + UnitCell *cell, const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); #endif /* DIRAX_H */ diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index d54cf43f..7cbe387a 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -175,4 +175,6 @@ extern struct imagefeature *image_feature_closest(ImageFeatureList *flist, extern int image_feature_count(ImageFeatureList *flist); extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); +extern void image_add_crystal(struct image *image, Crystal *cryst); + #endif /* IMAGE_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index fd4384d7..2370e26e 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -61,7 +61,7 @@ static const char *maybes(int n) IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, const char *filename, struct detector *det, - double nominal_photon_energy, float *ltl) + struct beam_params *beam, float *ltl) { int n; int nm = 0; @@ -76,7 +76,8 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, switch ( indm[n] & INDEXING_METHOD_MASK ) { case INDEXING_DIRAX : - iprivs[n] = NULL; + iprivs[n] = dirax_prepare(indm[nm], cell, filename, det, + beam, ltl); break; case INDEXING_MOSFLM : @@ -84,7 +85,7 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, break; case INDEXING_REAX : - iprivs[n] = reax_prepare(cell, filename, det); + iprivs[n] = reax_prepare(cell, filename, det, beam); break; default : @@ -101,35 +102,30 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, } -void cleanup_indexing(IndexingPrivate **priv) +void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) { int n = 0; - if ( priv == NULL ) return; /* Nothing to do */ + if ( indms == NULL ) return; /* Nothing to do */ + if ( privs == NULL ) return; /* Nothing to do */ - while ( priv[n] != NULL ) { + while ( indms[n] != INDEXING_NONE ) { - switch ( priv[n]->indm & INDEXING_METHOD_MASK ) { + switch ( indms[n] & INDEXING_METHOD_MASK ) { case INDEXING_NONE : - free(priv[n]); - break; - case INDEXING_DIRAX : - free(priv[n]); - break; - case INDEXING_MOSFLM : - free(priv[n]); + /* No cleanup */ break; case INDEXING_REAX : - reax_cleanup(priv[n]); + reax_cleanup(privs[n]); break; default : ERROR("Don't know how to clean up indexing method %i\n", - priv[n]->indm); + indms[n]); break; } @@ -162,19 +158,43 @@ void map_all_peaks(struct image *image) } -static void try_indexer(struct image *image, IndexingMethod indm, +/* Return non-zero for "success" */ +static int try_indexer(struct image *image, IndexingMethod indm, IndexingPrivate *ipriv) { + switch ( indm & INDEXING_METHOD_MASK ) { + + case INDEXING_NONE : + break; + + case INDEXING_DIRAX : + return run_dirax(image, ipriv); + break; + + case INDEXING_MOSFLM : + return run_mosflm(image, ipriv); + break; + + case INDEXING_REAX : + return reax_index(image, ipriv); + break; + + default : + ERROR("Unrecognised indexing method: %i\n", indm); + break; + + } + + return 0; } void index_pattern(struct image *image, IndexingMethod *indms, IndexingPrivate **iprivs) { - int i; int n = 0; - if ( indm == NULL ) return; + if ( indms == NULL ) return; map_all_peaks(image); image->crystals = NULL; @@ -182,7 +202,7 @@ void index_pattern(struct image *image, while ( indms[n] != INDEXING_NONE ) { - if ( try_indexer(image, indms[n], iprivs[i]) ) break; + if ( try_indexer(image, indms[n], iprivs[n]) ) break; n++; } @@ -216,7 +236,8 @@ static IndexingMethod set_bad(IndexingMethod a) /* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */ static IndexingMethod set_axes(IndexingMethod a) { - return (a & ~INDEXING_CHECK_COMBINATIONS) | INDEXING_CHECK_CELL_AXES; + return (a & ~INDEXING_CHECK_CELL_COMBINATIONS) + | INDEXING_CHECK_CELL_AXES; } @@ -234,7 +255,7 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) n = assplode(str, ",-", &methods, ASSPLODE_NONE); list = malloc((n+1)*sizeof(IndexingMethod)); - *nmeth = -1; /* So that the first method is #0 */ + nmeth = -1; /* So that the first method is #0 */ for ( i=0; i<n; i++ ) { if ( strcmp(methods[i], "dirax") == 0) { diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index d5b01b43..e0dcb8e0 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -87,15 +87,11 @@ extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, const char *filename, struct detector *det, - double nominal_photon_energy, - float *ltl); + struct beam_params *beam, float *ltl); extern void index_pattern(struct image *image, IndexingMethod *indms, IndexingPrivate **iprivs); -extern void cleanup_indexing(IndexingPrivate **priv); - -/* For indexing methods to call */ -extern IndexingPrivate *indexing_private(IndexingMethod indm); +extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs); #endif /* INDEX_H */ diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 548e3834..8b300012 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -80,6 +80,7 @@ #include "mosflm.h" #include "utils.h" #include "peaks.h" +#include "cell-utils.h" #define MOSFLM_VERBOSE 0 @@ -92,6 +93,14 @@ typedef enum { } MOSFLMInputType; + +struct mosflm_private { + IndexingMethod indm; + float *ltl; + UnitCell *template; +}; + + struct mosflm_data { /* MOSFLM auto-indexing low-level stuff */ @@ -107,10 +116,53 @@ struct mosflm_data { char sptfile[128]; int step; int finished_ok; - UnitCell *target_cell; + int done; + int success; + + struct mosflm_private *mp; }; +static int check_cell(struct mosflm_private *mp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + if ( mp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + + out = match_cell(cell, mp->template, 0, mp->ltl, 1); + if ( out == NULL ) return 0; + + } else if ( mp->indm & INDEXING_CHECK_CELL_AXES ) { + + out = match_cell(cell, mp->template, 0, mp->ltl, 0); + if ( out == NULL ) return 0; + + } else { + out = cell_new_from_cell(cell); + } + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, cell); + + if ( mp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + static void mosflm_parseline(const char *line, struct image *image, struct mosflm_data *dirax) @@ -130,7 +182,8 @@ static void mosflm_parseline(const char *line, struct image *image, } -static int read_newmat(const char *filename, struct image *image) +static int read_newmat(struct mosflm_data *mosflm, const char *filename, + struct image *image) { FILE * fh; float asx, asy, asz; @@ -138,6 +191,7 @@ static int read_newmat(const char *filename, struct image *image) float csx, csy, csz; int n; double c; + UnitCell *cell; fh = fopen(filename, "r"); if ( fh == NULL ) { @@ -155,18 +209,21 @@ static int read_newmat(const char *filename, struct image *image) /* MOSFLM "A" matrix is multiplied by lambda, so fix this */ c = 1.0/image->lambda; - image->candidate_cells[0] = cell_new(); + 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(image->candidate_cells[0], + cell_set_reciprocal(cell, -asy*c, -asz*c, asx*c, -bsy*c, -bsz*c, bsx*c, -csy*c, -csz*c, csx*c); - image->ncells = 1; + if ( check_cell(mosflm->mp, image, cell) ) { + mosflm->success = 1; + mosflm->done = 1; + } return 0; } @@ -352,9 +409,9 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) break; case 2 : - if ( mosflm->target_cell != NULL ) { + if ( mosflm->mp->template != NULL ) { const char *symm; - symm = spacegroup_for_lattice(mosflm->target_cell); + symm = spacegroup_for_lattice(mosflm->mp->template); snprintf(tmp, 255, "SYMM %s\n", symm); mosflm_sendline(tmp, mosflm); } else { @@ -527,7 +584,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) } -void run_mosflm(struct image *image, UnitCell *cell) +int run_mosflm(struct image *image, IndexingPrivate *ipriv) { struct mosflm_data *mosflm; unsigned int opts; @@ -537,11 +594,9 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm = malloc(sizeof(struct mosflm_data)); if ( mosflm == NULL ) { ERROR("Couldn't allocate memory for MOSFLM data.\n"); - return; + return 0; } - mosflm->target_cell = cell; - snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id); write_img(image, mosflm->imagefile); /* Dummy image */ @@ -556,7 +611,7 @@ void run_mosflm(struct image *image, UnitCell *cell) if ( mosflm->pid == -1 ) { ERROR("Failed to fork for MOSFLM: %s\n", strerror(errno)); free(mosflm); - return; + return 0; } if ( mosflm->pid == 0 ) { @@ -584,6 +639,9 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm->step = 1; /* This starts the "initialisation" procedure */ mosflm->finished_ok = 0; + mosflm->mp = (struct mosflm_private *)ipriv; + mosflm->done = 0; + mosflm->success = 0; do { @@ -632,8 +690,27 @@ void run_mosflm(struct image *image, UnitCell *cell) 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); + read_newmat(mosflm, mosflm->newmatfile, image); } + rval = mosflm->success; free(mosflm); + return rval; +} + + +IndexingPrivate *mosflm_prepare(IndexingMethod indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) +{ + struct mosflm_private *mp; + + mp = malloc(sizeof(struct mosflm_private)); + if ( mp == NULL ) return NULL; + + mp->ltl = ltl; + mp->template = cell; + mp->indm = indm; + + return (IndexingPrivate *)mp; } diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h index ba1eb73d..1d296cb7 100644 --- a/libcrystfel/src/mosflm.h +++ b/libcrystfel/src/mosflm.h @@ -35,10 +35,10 @@ #include <config.h> #endif -#include "utils.h" +#include "index.h" -extern void run_mosflm(struct image *image, UnitCell *cell); +extern int run_mosflm(struct image *image, IndexingPrivate *ipriv); #endif /* MOSFLM_H */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index c350388b..08595ef5 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -594,7 +594,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } -int peak_sanity_check(struct image *image) +int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { int n_feat = 0; int n_sane = 0; @@ -616,13 +616,13 @@ int peak_sanity_check(struct image *image) /* Reciprocal space position of found peak */ q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); - for ( j=0; j<image->n_crystals; j++ ) { + for ( j=0; j<n_cryst; j++ ) { double ax, ay, az; double bx, by, bz; double cx, cy, cz; - cell_get_cartesian(crystal_get_cell(image->crystals[j]), + cell_get_cartesian(crystal_get_cell(crystals[j]), &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 40136eca..b00edf00 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -47,7 +47,8 @@ extern void integrate_reflections(struct image *image, double ir_inn, double ir_mid, double ir_out, int integrate_saturated); -extern int peak_sanity_check(struct image *image); +extern int peak_sanity_check(struct image *image, Crystal **crystals, + int n_cryst); extern void estimate_resolution(RefList *list, UnitCell *cell, double *min, double *max); diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h index 01dfe649..af383e29 100644 --- a/libcrystfel/src/reax.h +++ b/libcrystfel/src/reax.h @@ -34,16 +34,24 @@ #endif #include "cell.h" +#include "beam-parameters.h" +#include "detector.h" #ifdef HAVE_FFTW -extern IndexingPrivate *reax_prepare(void); +extern IndexingPrivate *reax_prepare(UnitCell *cell, const char *filename, + struct detector *det, + struct beam_params *beam); + extern void reax_cleanup(IndexingPrivate *pp); -extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell); + +extern int reax_index(struct image *image, IndexingPrivate *p); #else /* HAVE_FFTW */ -static IndexingPrivate *reax_prepare() +static IndexingPrivate *reax_prepare(UnitCell *cell, const char *filename, + struct detector *det, + struct beam_params *beam) { return NULL; } @@ -52,7 +60,7 @@ static void reax_cleanup(IndexingPrivate *pp) { } -static void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell) +static void reax_index(struct image *image, IndexingPrivate *p) { } |