diff options
Diffstat (limited to 'libcrystfel/src')
-rw-r--r-- | libcrystfel/src/beam-parameters.h | 8 | ||||
-rw-r--r-- | libcrystfel/src/cell-utils.c | 8 | ||||
-rw-r--r-- | libcrystfel/src/crystal.c | 224 | ||||
-rw-r--r-- | libcrystfel/src/crystal.h | 74 | ||||
-rw-r--r-- | libcrystfel/src/dirax.c | 153 | ||||
-rw-r--r-- | libcrystfel/src/dirax.h | 15 | ||||
-rw-r--r-- | libcrystfel/src/geometry.c | 110 | ||||
-rw-r--r-- | libcrystfel/src/geometry.h | 9 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 18 | ||||
-rw-r--r-- | libcrystfel/src/image.h | 62 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 340 | ||||
-rw-r--r-- | libcrystfel/src/index.h | 73 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 260 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.h | 16 | ||||
-rw-r--r-- | libcrystfel/src/peaks.c | 242 | ||||
-rw-r--r-- | libcrystfel/src/peaks.h | 12 | ||||
-rw-r--r-- | libcrystfel/src/reax.c | 89 | ||||
-rw-r--r-- | libcrystfel/src/reax.h | 17 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 540 | ||||
-rw-r--r-- | libcrystfel/src/stream.h | 43 | ||||
-rw-r--r-- | libcrystfel/src/symmetry.c | 17 |
21 files changed, 1645 insertions, 685 deletions
diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h index 8212811b..de777deb 100644 --- a/libcrystfel/src/beam-parameters.h +++ b/libcrystfel/src/beam-parameters.h @@ -3,11 +3,11 @@ * * Beam parameters * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-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> * 2012 Chunhong Yoon * * This file is part of CrystFEL. @@ -34,6 +34,8 @@ #include <config.h> #endif +struct beam_params; + #include "hdf5-file.h" struct beam_params diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 3ec9a9d6..ee51622b 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -424,14 +424,18 @@ UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) tt = uncentering_transformation(in, &new_centering, &new_latt); if ( tt == NULL ) return NULL; - if ( t != NULL ) *t = tt; - out = cell_transform(in, tt); if ( out == NULL ) return NULL; cell_set_lattice_type(out, new_latt); cell_set_centering(out, new_centering); + if ( t != NULL ) { + *t = tt; + } else { + tfn_free(tt); + } + return out; } diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c new file mode 100644 index 00000000..664d933e --- /dev/null +++ b/libcrystfel/src/crystal.c @@ -0,0 +1,224 @@ +/* + * crystal.c + * + * A class representing a single crystal + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2013 Thomas White <taw@physics.org> + * + * 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 <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "crystal.h" +#include "utils.h" + + +/** + * SECTION:crystal + * @short_description: Crystal + * @title: Crystal + * @section_id: + * @see_also: + * @include: "crystal.h" + * @Image: + * + * This structure represents a single crystal. + */ + + +struct _crystal +{ + /* The image containing the crystal */ + struct image *image; + + /* Information about the crystal */ + UnitCell *cell; + double m; /* Mosaicity in radians */ + double osf; + double profile_radius; + int pr_dud; + double resolution_limit; + + /* Integrated (or about-to-be-integrated) reflections */ + RefList *reflections; + long long int n_saturated; /* Number of overloads */ + + /* User flag, e.g. for "this is a bad crystal". */ + int user_flag; +}; + + +/************************** Constructor / Destructor **************************/ + + +/** + * crystal_new: + * + * Create a new %Crystal. + * + * Returns: the new unit cell, or NULL on failure. + * + */ +Crystal *crystal_new() +{ + Crystal *cryst; + + cryst = malloc(sizeof(Crystal)); + if ( cryst == NULL ) return NULL; + + cryst->cell = NULL; + cryst->reflections = NULL; + cryst->resolution_limit = 0.0; + cryst->n_saturated = 0; + + return cryst; +} + + +/** + * crystal_free: + * @cryst: A %Crystal to free. + * + * Frees a %Crystal, and all internal resources concerning that crystal. + * + */ +void crystal_free(Crystal *cryst) +{ + if ( cryst == NULL ) return; + free(cryst); +} + + +/********************************** Getters ***********************************/ + + +UnitCell *crystal_get_cell(Crystal *cryst) +{ + return cryst->cell; +} + + +double crystal_get_profile_radius(Crystal *cryst) +{ + return cryst->profile_radius; +} + + +RefList *crystal_get_reflections(Crystal *cryst) +{ + return cryst->reflections; +} + + +double crystal_get_resolution_limit(Crystal *cryst) +{ + return cryst->resolution_limit; +} + + +long long int crystal_get_num_saturated_reflections(Crystal *cryst) +{ + return cryst->n_saturated; +} + + +struct image *crystal_get_image(Crystal *cryst) +{ + return cryst->image; +} + + +double crystal_get_osf(Crystal *cryst) +{ + return cryst->osf; +} + + +int crystal_get_user_flag(Crystal *cryst) +{ + return cryst->user_flag; +} + + +double crystal_get_mosaicity(Crystal *cryst) +{ + return cryst->m; +} + + +/********************************** Setters ***********************************/ + + +void crystal_set_cell(Crystal *cryst, UnitCell *cell) +{ + cryst->cell = cell; +} + + +void crystal_set_profile_radius(Crystal *cryst, double r) +{ + cryst->profile_radius = r; +} + + +void crystal_set_reflections(Crystal *cryst, RefList *reflist) +{ + cryst->reflections = reflist; +} + + +void crystal_set_resolution_limit(Crystal *cryst, double res) +{ + cryst->resolution_limit = res; +} + + +void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n) +{ + cryst->n_saturated = n; +} + + +void crystal_set_image(Crystal *cryst, struct image *image) +{ + cryst->image = image; +} + + +void crystal_set_osf(Crystal *cryst, double osf) +{ + cryst->osf = osf; +} + + +void crystal_set_user_flag(Crystal *cryst, int user_flag) +{ + cryst->user_flag = user_flag; +} + + +void crystal_set_mosaicity(Crystal *cryst, double m) +{ + cryst->m = m; +} diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h new file mode 100644 index 00000000..3d0ad9d1 --- /dev/null +++ b/libcrystfel/src/crystal.h @@ -0,0 +1,74 @@ +/* + * crystal.h + * + * A class representing a single crystal + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2013 Thomas White <taw@physics.org> + * + * 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 <http://www.gnu.org/licenses/>. + * + */ + +#ifndef CRYSTAL_H +#define CRYSTAL_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include "cell.h" +#include "reflist.h" + + +/** + * Crystal: + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + **/ +typedef struct _crystal Crystal; + + +extern Crystal *crystal_new(void); +extern void crystal_free(Crystal *cryst); + +extern UnitCell *crystal_get_cell(Crystal *cryst); +extern double crystal_get_profile_radius(Crystal *cryst); +extern RefList *crystal_get_reflections(Crystal *cryst); +extern double crystal_get_resolution_limit(Crystal *cryst); +extern long long int crystal_get_num_saturated_reflections(Crystal *cryst); +extern int crystal_get_user_flag(Crystal *cryst); +extern double crystal_get_osf(Crystal *cryst); +extern struct image *crystal_get_image(Crystal *cryst); +extern double crystal_get_mosaicity(Crystal *cryst); + +extern void crystal_set_cell(Crystal *cryst, UnitCell *cell); +extern void crystal_set_profile_radius(Crystal *cryst, double r); +extern void crystal_set_reflections(Crystal *cryst, RefList *reflist); +extern void crystal_set_resolution_limit(Crystal *cryst, double res); +extern void crystal_set_num_saturated_reflections(Crystal *cryst, + long long int n); +extern void crystal_set_user_flag(Crystal *cryst, int flag); +extern void crystal_set_osf(Crystal *cryst, double osf); +extern void crystal_set_image(Crystal *cryst, struct image *image); +extern void crystal_set_mosaicity(Crystal *cryst, double m); + +#endif /* CRYSTAL_H */ diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 160df719..e6ff36b1 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -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-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -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 */ @@ -83,12 +91,59 @@ struct dirax_data { int read_cell; int best_acl; int best_acl_nh; - int acls_tried[MAX_CELL_CANDIDATES]; + int acls_tried[MAX_DIRAX_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, out); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + cell_free(out); + 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 +174,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 +182,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 +192,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 +211,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 +230,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 +419,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 +433,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 +520,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 +532,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 +569,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 +610,7 @@ void run_dirax(struct image *image) rval = 1; } - } while ( !rval ); + } while ( !rval && !dirax->success ); close(dirax->pty); free(dirax->rbuffer); @@ -553,5 +620,45 @@ 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; + int need_cell = 0; + + if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; + if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; + + if ( need_cell && (cell == NULL) ) { + ERROR("DirAx needs a unit cell for this set of flags.\n"); + return NULL; + } + + /* Flags that DirAx knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS + | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS; + + dp = malloc(sizeof(struct dirax_private)); + if ( dp == NULL ) return NULL; + + dp->ltl = ltl; + dp->template = cell; + dp->indm = *indm; + + return (IndexingPrivate *)dp; +} + + +void dirax_cleanup(IndexingPrivate *pp) +{ + struct dirax_private *p; + p = (struct dirax_private *)pp; + free(p); } diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h index ce3cd4b0..6be8451a 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,15 @@ #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); +extern void dirax_cleanup(IndexingPrivate *pp); #endif /* DIRAX_H */ diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 60fab488..704baa51 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -3,11 +3,11 @@ * * Geometry of diffraction * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -33,6 +33,7 @@ #include <stdlib.h> #include <assert.h> +#include <fenv.h> #include "utils.h" #include "cell.h" @@ -43,6 +44,7 @@ #include "reflist.h" #include "reflist-utils.h" #include "symmetry.h" +#include "geometry.h" static signed int locate_peak(double x, double y, double z, double k, @@ -117,7 +119,7 @@ static double partiality(double rlow, double rhigh, double r) } -static Reflection *check_reflection(struct image *image, +static Reflection *check_reflection(struct image *image, Crystal *cryst, signed int h, signed int k, signed int l, double xl, double yl, double zl) { @@ -131,6 +133,9 @@ static Reflection *check_reflection(struct image *image, double klow, khigh; /* Wavenumber */ Reflection *refl; double cet, cez; + double pr; + + pr = crystal_get_profile_radius(cryst); /* "low" gives the largest Ewald sphere (wavelength short => k large) * "high" gives the smallest Ewald sphere (wavelength long => k small) @@ -152,8 +157,8 @@ static Reflection *check_reflection(struct image *image, rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ if ( (signbit(rlow) == signbit(rhigh)) - && (fabs(rlow) > image->profile_radius) - && (fabs(rhigh) > image->profile_radius) ) return NULL; + && (fabs(rlow) > pr) + && (fabs(rhigh) > pr) ) return NULL; /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the @@ -164,26 +169,26 @@ static Reflection *check_reflection(struct image *image, * et al. (1979). */ clamp_low = 0; clamp_high = 0; - if ( rlow < -image->profile_radius ) { - rlow = -image->profile_radius; + if ( rlow < -pr ) { + rlow = -pr; clamp_low = -1; } - if ( rlow > +image->profile_radius ) { - rlow = +image->profile_radius; + if ( rlow > +pr ) { + rlow = +pr; clamp_low = +1; } - if ( rhigh < -image->profile_radius ) { - rhigh = -image->profile_radius; + if ( rhigh < -pr ) { + rhigh = -pr; clamp_high = -1; } - if ( rhigh > +image->profile_radius ) { - rhigh = +image->profile_radius; + if ( rhigh > +pr ) { + rhigh = +pr; clamp_high = +1; } assert(clamp_low >= clamp_high); /* Calculate partiality */ - part = partiality(rlow, rhigh, image->profile_radius); + part = partiality(rlow, rhigh, pr); /* Locate peak on detector. */ p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det, &xda, &yda); @@ -205,7 +210,7 @@ static Reflection *check_reflection(struct image *image, } -RefList *find_intersections(struct image *image, UnitCell *cell) +RefList *find_intersections(struct image *image, Crystal *cryst) { double ax, ay, az; double bx, by, bz; @@ -217,6 +222,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell) int hmax, kmax, lmax; double mres; signed int h, k, l; + UnitCell *cell; + + cell = crystal_get_cell(cryst); + if ( cell == NULL ) return NULL; reflections = reflist_new(); @@ -262,7 +271,7 @@ RefList *find_intersections(struct image *image, UnitCell *cell) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - refl = check_reflection(image, h, k, l, xl, yl, zl); + refl = check_reflection(image, cryst, h, k, l, xl, yl, zl); if ( refl != NULL ) { add_refl_to_list(refl, reflections); @@ -276,8 +285,68 @@ RefList *find_intersections(struct image *image, UnitCell *cell) } +RefList *select_intersections(struct image *image, Crystal *cryst) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + const double min_dist = 0.25; + RefList *list; + int i; + + /* Round towards nearest */ + fesetround(1); + + /* Cell basis vectors for this image */ + cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az, + &bx, &by, &bz, &cx, &cy, &cz); + + list = reflist_new(); + if ( list == NULL ) return NULL; + + /* Loop over peaks, checking proximity to nearest reflection */ + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f; + struct rvec q; + double h, k, l, hd, kd, ld; + double dsq; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + /* Reciprocal space position of found peak */ + q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); + + /* Decimal and fractional Miller indices of nearest + * reciprocal lattice point */ + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + /* Check distance */ + dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); + + if ( sqrt(dsq) < min_dist ) { + + Reflection *refl; + + refl = add_refl(list, h, k, l); + set_detector_pos(refl, sqrt(dsq), f->fs, f->ss); + + } + + } + + return list; +} + + /* Calculate partialities and apply them to the image's reflections */ -void update_partialities(struct image *image) +void update_partialities(Crystal *cryst) { Reflection *refl; RefListIterator *iter; @@ -285,14 +354,15 @@ void update_partialities(struct image *image) double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; + struct image *image = crystal_get_image(cryst); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); /* Scratch list to give check_reflection() something to add to */ predicted = reflist_new(); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -309,7 +379,7 @@ void update_partialities(struct image *image) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - vals = check_reflection(image, h, k, l, xl, yl, zl); + vals = check_reflection(image, cryst, h, k, l, xl, yl, zl); if ( vals == NULL ) { set_redundancy(refl, 0); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 94b496e6..aecdc28a 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -3,12 +3,12 @@ * * Geometry of diffraction * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2012 Richard Kirian * * This file is part of CrystFEL. @@ -43,9 +43,10 @@ extern "C" { #endif -RefList *find_intersections(struct image *image, UnitCell *cell); +extern RefList *find_intersections(struct image *image, Crystal *cryst); +extern RefList *select_intersections(struct image *image, Crystal *cryst); -void update_partialities(struct image *image); +extern void update_partialities(Crystal *cryst); #ifdef __cplusplus } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 093e20f2..6e4fd40a 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -159,3 +159,21 @@ void image_remove_feature(ImageFeatureList *flist, int idx) { flist->features[idx].valid = 0; } + + +void image_add_crystal(struct image *image, Crystal *cryst) +{ + Crystal **crs; + int n; + + n = image->n_crystals; + crs = realloc(image->crystals, (n+1)*sizeof(Crystal *)); + if ( crs == NULL ) { + ERROR("Failed to allocate memory for crystals.\n"); + return; + } + + crs[n] = cryst; + image->crystals = crs; + image->n_crystals = n+1; +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index afa9e4a7..3739c01f 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -3,11 +3,11 @@ * * Handle images and image features * - * 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: - * 2009-2012 Thomas White <taw@physics.org> + * 2009-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -41,9 +41,8 @@ #include "cell.h" #include "detector.h" #include "reflist.h" - - -#define MAX_CELL_CANDIDATES (32) +#include "crystal.h" +#include "index.h" /* Structure describing a feature in an image */ @@ -79,10 +78,10 @@ typedef struct _imagefeaturelist ImageFeatureList; * uint16_t *flags; * double *twotheta; * - * UnitCell *indexed_cell; - * UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; - * int ncells; - + * Crystal **crystals; + * int n_crystals; + * IndexingMethod indexed_by; + * * struct detector *det; * struct beam_params *beam; * char *filename; @@ -90,22 +89,15 @@ typedef struct _imagefeaturelist ImageFeatureList; * * int id; * - * double m; - * * double lambda; * double div; * double bw; - * double osf; - * double profile_radius; - * int pr_dud; - * double diffracting_resolution; * * int width; * int height; * - * RefList *reflections; - * long long unsigned int n_saturated; - * + * long long int num_peaks; + * long long int num_saturated_peaks; * ImageFeatureList *features; * }; * </programlisting> @@ -121,14 +113,9 @@ typedef struct _imagefeaturelist ImageFeatureList; * by-product of the scattering vector calculation and can be used later for * calculating intensities from differential scattering cross sections. * - * <structfield>candidate_cells</structfield> is an array of unit cells directly - * returned by the low-level indexing system. <structfield>ncells</structfield> - * is the number of candidate unit cells which were found. The maximum number - * of cells which may be returned is <function>MAX_CELL_CANDIDATES</function>. - * <structfield>indexed_cell</structfield> contains the "correct" unit cell - * after cell reduction or matching has been performed. The job of the cell - * reduction is to convert the list of candidate cells into a single indexed - * cell, or <function>NULL</function> on failure. + * <structfield>crystals</structfield> is an array of %Crystal directly + * returned by the low-level indexing system. <structfield>n_crystals</structfield> + * is the number of crystals which were found in the image. * * <structfield>copyme</structfield> represents a list of HDF5 fields to copy * to the output stream. @@ -141,9 +128,9 @@ struct image { uint16_t *flags; double *twotheta; - UnitCell *indexed_cell; - UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; - int ncells; + Crystal **crystals; + int n_crystals; + IndexingMethod indexed_by; struct detector *det; struct beam_params *beam; /* The nominal beam parameters */ @@ -153,26 +140,17 @@ struct image { int id; /* ID number of the thread * handling this image */ - /* Information about the crystal */ - double m; /* Mosaicity in radians */ - /* Per-shot radiation values */ double lambda; /* Wavelength in m */ double div; /* Divergence in radians */ double bw; /* Bandwidth as a fraction */ - double osf; /* Overall scaling factor */ - double profile_radius; /* Radius of reflection */ - int pr_dud; /* Post refinement failed */ - double diffracting_resolution; /* Max 1/d in m^-1 */ int width; int height; - /* Integrated (or about-to-be-integrated) reflections */ - RefList *reflections; - long long int n_saturated; /* Number of overloads */ - /* Detected peaks */ + long long int num_peaks; + long long int num_saturated_peaks; ImageFeatureList *features; }; @@ -196,4 +174,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 86ab8c1a..89fdec4f 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -3,12 +3,12 @@ * * Perform indexing (somehow) * - * 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. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2010-2011 Richard Kirian <rkirian@asu.edu> * 2012 Lorenzo Galli * @@ -53,49 +53,34 @@ #include "cell-utils.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) + struct beam_params *beam, float *ltl) { 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; n<nm; n++ ) { - switch ( indm[n] ) { + int i; + IndexingMethod in; - case INDEXING_NONE : - ERROR("Tried to prepare INDEXING_NONE!\n"); - break; + in = indm[n]; + + switch ( indm[n] & INDEXING_METHOD_MASK ) { case INDEXING_DIRAX : - iprivs[n] = indexing_private(indm[n]); + iprivs[n] = dirax_prepare(&indm[n], cell, filename, + det, beam, ltl); break; case INDEXING_MOSFLM : - iprivs[n] = indexing_private(indm[n]); + iprivs[n] = mosflm_prepare(&indm[n], cell, filename, + det, beam, ltl); break; case INDEXING_GRAINSPOTTER : @@ -103,11 +88,38 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, break; case INDEXING_REAX : - iprivs[n] = reax_prepare(); + iprivs[n] = reax_prepare(&indm[n], cell, filename, + det, beam, ltl); + break; + + default : + ERROR("Don't know how to prepare indexing method %i\n", + indm[n]); break; } + if ( iprivs[n] == NULL ) return NULL; + + STATUS("Prepared indexing method %i: %s\n", + n, indexer_str(indm[n])); + + if ( in != indm[n] ) { + ERROR("Note: flags were altered to take into account " + "the limitations of the indexing method.\n"); + } + + for ( i=0; i<n; i++ ) { + if ( indm[i] == indm[n] ) { + ERROR("Duplicate indexing method.\n"); + ERROR("Have you specified some flags which " + "aren't accepted by one of your " + "chosen indexing methods?\n"); + return NULL; + } + } + + } iprivs[n] = NULL; @@ -115,26 +127,26 @@ 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 ) { + switch ( indms[n] & INDEXING_METHOD_MASK ) { case INDEXING_NONE : - free(priv[n]); break; case INDEXING_DIRAX : - free(priv[n]); + dirax_cleanup(privs[n]); break; case INDEXING_MOSFLM : - free(priv[n]); + mosflm_cleanup(privs[n]); break; case INDEXING_GRAINSPOTTER : @@ -142,7 +154,12 @@ void cleanup_indexing(IndexingPrivate **priv) 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", + indms[n]); break; } @@ -175,146 +192,213 @@ void map_all_peaks(struct image *image) } -void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, - int cellr, int verbose, IndexingPrivate **ipriv, - int config_insane, const float *ltl) +/* 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->indexed_cell = NULL; + image->crystals = NULL; + image->n_crystals = 0; - while ( indm[n] != INDEXING_NONE ) { + while ( indms[n] != INDEXING_NONE ) { - image->ncells = 0; + if ( try_indexer(image, indms[n], iprivs[n]) ) break; + n++; - /* Index as appropriate */ - switch ( indm[n] ) { + } - case INDEXING_NONE : - return; + image->indexed_by = indms[n]; +} - case INDEXING_DIRAX : - run_dirax(image); - break; - case INDEXING_MOSFLM : - run_mosflm(image, cell); - break; +/* Set the indexer flags for "raw mode" ("--cell-reduction=none") */ +static IndexingMethod set_raw(IndexingMethod a) +{ + /* Disable all unit cell checks */ + a &= ~(INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_CELL_AXES); + return a; +} - case INDEXING_GRAINSPOTTER : - run_grainspotter(image, cell); - break; - case INDEXING_REAX : - reax_index(ipriv[n], image, cell); - break; +/* Set the indexer flags for "bad mode" ("--insane) */ +static IndexingMethod set_bad(IndexingMethod a) +{ + /* Disable the peak check */ + return a & ~INDEXING_CHECK_PEAKS; +} - } - if ( image->ncells == 0 ) { - n++; - continue; - } - for ( i=0; i<image->ncells; i++ ) { +/* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */ +static IndexingMethod set_axes(IndexingMethod a) +{ + return (a & ~INDEXING_CHECK_CELL_COMBINATIONS) + | INDEXING_CHECK_CELL_AXES; +} - 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"); - } +/* Set the indexer flags for "combination mode" ("--cell-reduction=reduce") */ +static IndexingMethod set_comb(IndexingMethod a) +{ + return (a & ~INDEXING_CHECK_CELL_AXES) + | INDEXING_CHECK_CELL_COMBINATIONS; +} - /* Match or reduce the cell as appropriate */ - switch ( cellr ) { - case CELLR_NONE : - new_cell = cell_new_from_cell(cand); - break; +/* Set the indexer flags for "use no lattice type information" */ +static IndexingMethod set_nolattice(IndexingMethod a) +{ + return a & ~INDEXING_USE_LATTICE_TYPE; +} - case CELLR_REDUCE : - new_cell = match_cell(cand, cell, verbose, - ltl, 1); - break; - case CELLR_COMPARE : - new_cell = match_cell(cand, cell, verbose, - ltl, 0); - break; +/* Set the indexer flags for "use lattice type information" */ +static IndexingMethod set_lattice(IndexingMethod a) +{ + return a | INDEXING_USE_LATTICE_TYPE; +} - case CELLR_COMPARE_AB : - new_cell = match_cell_ab(cand, cell); - break; - } +char *indexer_str(IndexingMethod indm) +{ + char *str; - /* No cell? Move on to the next candidate */ - if ( new_cell == NULL ) continue; + str = malloc(32); + if ( str == NULL ) { + ERROR("Failed to allocate string.\n"); + return NULL; + } + str[0] = '\0'; - /* Sanity check */ - image->indexed_cell = new_cell; - if ( !config_insane && !peak_sanity_check(image) ) { - cell_free(new_cell); - image->indexed_cell = NULL; - continue; - } + switch ( indm & INDEXING_METHOD_MASK ) { - goto done; /* Success */ + case INDEXING_NONE : + strcpy(str, "none"); + return str; - } + case INDEXING_DIRAX : + strcpy(str, "dirax"); + break; - for ( i=0; i<image->ncells; i++ ) { - cell_free(image->candidate_cells[i]); - image->candidate_cells[i] = NULL; - } + case INDEXING_MOSFLM : + strcpy(str, "mosflm"); + break; - /* Move on to the next indexing method */ - n++; + case INDEXING_REAX : + strcpy(str, "reax"); + break; + + default : + ERROR("Unrecognised indexing method %i\n", indm); + strcpy(str, "(unknown)"); + break; } -done: - for ( i=0; i<image->ncells; i++ ) { - /* May free(NULL) if all algorithms were tried and no success */ - cell_free(image->candidate_cells[i]); + if ( indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + strcat(str, "-comb"); + } else if ( indm & INDEXING_CHECK_CELL_AXES ) { + strcat(str, "-axes"); + } else { + strcat(str, "-raw"); } + + if ( !(indm & INDEXING_CHECK_PEAKS) ) { + strcat(str, "-bad"); + } + + if ( indm & INDEXING_USE_LATTICE_TYPE ) { + strcat(str, "-latt"); + } else { + strcat(str, "-nolatt"); + } + + return str; } -IndexingMethod *build_indexer_list(const char *str, int *need_cell) +IndexingMethod *build_indexer_list(const char *str) { int n, i; char **methods; IndexingMethod *list; - int tmp; - - if ( need_cell == NULL ) need_cell = &tmp; - *need_cell = 0; + int nmeth = 0; - n = assplode(str, ",", &methods, ASSPLODE_NONE); + n = assplode(str, ",-", &methods, ASSPLODE_NONE); list = malloc((n+1)*sizeof(IndexingMethod)); + nmeth = -1; /* So that the first method is #0 */ for ( i=0; i<n; i++ ) { if ( strcmp(methods[i], "dirax") == 0) { - list[i] = INDEXING_DIRAX; + list[++nmeth] = INDEXING_DEFAULTS_DIRAX; + } else if ( strcmp(methods[i], "mosflm") == 0) { - list[i] = INDEXING_MOSFLM; + list[++nmeth] = INDEXING_DEFAULTS_MOSFLM; + } else if ( strcmp(methods[i], "grainspotter") == 0) { - list[i] = INDEXING_GRAINSPOTTER; + list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; + } else if ( strcmp(methods[i], "reax") == 0) { - list[i] = INDEXING_REAX; - *need_cell = 1; + list[++nmeth] = INDEXING_DEFAULTS_REAX; + + } else if ( strcmp(methods[i], "none") == 0) { + list[++nmeth] = INDEXING_NONE; + return list; + + } else if ( strcmp(methods[i], "raw") == 0) { + list[nmeth] = set_raw(list[nmeth]); + + } else if ( strcmp(methods[i], "bad") == 0) { + list[nmeth] = set_bad(list[nmeth]); + + } else if ( strcmp(methods[i], "comb") == 0) { + list[nmeth] = set_comb(list[nmeth]); /* Default */ + + } else if ( strcmp(methods[i], "axes") == 0) { + list[nmeth] = set_axes(list[nmeth]); + + } else if ( strcmp(methods[i], "latt") == 0) { + list[nmeth] = set_lattice(list[nmeth]); + + } else if ( strcmp(methods[i], "nolatt") == 0) { + list[nmeth] = set_nolattice(list[nmeth]); + } else { - ERROR("Unrecognised indexing method '%s'\n", - methods[i]); + ERROR("Bad list of indexing methods: '%s'\n", str); return NULL; } @@ -322,7 +406,7 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) } free(methods); - list[i] = INDEXING_NONE; + list[++nmeth] = INDEXING_NONE; return list; } diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 2df7a311..a47ee9a0 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -3,13 +3,13 @@ * * Perform indexing (somehow) * - * 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. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2010 Richard Kirian * 2012 Lorenzo Galli * @@ -38,11 +38,22 @@ #endif +#include "beam-parameters.h" #include "cell.h" #include "image.h" #include "detector.h" +#define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX | INDEXING_CHECK_PEAKS \ + | INDEXING_CHECK_CELL_COMBINATIONS) + +#define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_CHECK_PEAKS \ + | INDEXING_CHECK_CELL_COMBINATIONS \ + | INDEXING_USE_LATTICE_TYPE) + +#define INDEXING_DEFAULTS_REAX (INDEXING_REAX | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_CHECK_PEAKS) + /** * IndexingMethod: * @INDEXING_NONE: No indexing to be performed @@ -53,21 +64,28 @@ * An enumeration of all the available indexing methods. **/ typedef enum { - INDEXING_NONE, - INDEXING_DIRAX, - INDEXING_MOSFLM, - INDEXING_GRAINSPOTTER, - INDEXING_REAX, -} IndexingMethod; -/* Cell reduction methods */ -enum { - CELLR_NONE, - CELLR_REDUCE, - CELLR_COMPARE, - CELLR_COMPARE_AB, -}; + INDEXING_NONE = 0, + + /* The core indexing methods themselves */ + INDEXING_DIRAX = 1, + INDEXING_MOSFLM = 2, + INDEXING_REAX = 3, + INDEXING_GRAINSPOTTER = 4, + + /* Bits at the top of the IndexingMethod are flags which modify the + * behaviour of the indexer, at the moment just by adding checks. */ + INDEXING_CHECK_CELL_COMBINATIONS = 256, + INDEXING_CHECK_CELL_AXES = 512, + INDEXING_CHECK_PEAKS = 1024, + INDEXING_USE_LATTICE_TYPE = 2048, + +} IndexingMethod; + +/* This defines the bits in "IndexingMethod" which are used to represent the + * core of the indexing method */ +#define INDEXING_METHOD_MASK (0xff) /** @@ -76,24 +94,19 @@ enum { * This is an opaque data structure containing information needed by the * indexing method. **/ -typedef struct _indexingprivate IndexingPrivate; +typedef void *IndexingPrivate; -extern IndexingPrivate *indexing_private(IndexingMethod indm); +extern IndexingMethod *build_indexer_list(const char *str); +extern char *indexer_str(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, - const float *ltl); + const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); -extern void cleanup_indexing(IndexingPrivate **priv); +extern void index_pattern(struct image *image, + IndexingMethod *indms, IndexingPrivate **iprivs); -extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); +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..a03e621d 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * 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. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian <rkirian@asu.edu> - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -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,85 @@ 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 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 ( 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 ( 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, out); + + if ( mp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + cell_free(out); + crystal_free(cr); + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + static void mosflm_parseline(const char *line, struct image *image, struct mosflm_data *dirax) @@ -130,14 +214,54 @@ static void mosflm_parseline(const char *line, struct image *image, } -static int read_newmat(const char *filename, struct image *image) +/* This is the opposite of spacegroup_for_lattice() below. */ +static LatticeType spacegroup_to_lattice(const char *sg) +{ + LatticeType latt; + + if ( sg[1] == '1' ) { + latt = L_TRICLINIC; + } else if ( strncmp(sg+1, "23", 2) == 0 ) { + latt = L_CUBIC; + } else if ( strncmp(sg+1, "222", 3) == 0 ) { + latt = L_ORTHORHOMBIC; + } else if ( sg[1] == '2' ) { + latt = L_MONOCLINIC; + } else if ( sg[1] == '4' ) { + latt = L_TETRAGONAL; + } else if ( sg[1] == '6' ) { + latt = L_HEXAGONAL; + } else if ( sg[1] == '3' ) { + if ( sg[0] == 'H' ) { + latt = L_HEXAGONAL; + } else { + latt = L_RHOMBOHEDRAL; + } + } else { + ERROR("Couldn't understand '%s'\n", sg); + latt = L_TRICLINIC; + } + + return latt; +} + + + +static int read_newmat(struct mosflm_data *mosflm, const char *filename, + struct image *image) { - FILE * fh; + FILE *fh; float asx, asy, asz; float bsx, bsy, bsz; float csx, csy, csz; int n; double c; + UnitCell *cell; + char symm[32]; + char *rval; + int i; + char cen; + LatticeType latt; fh = fopen(filename, "r"); if ( fh == NULL ) { @@ -150,23 +274,54 @@ static int read_newmat(const char *filename, struct image *image) STATUS("Fewer than 9 parameters found in NEWMAT file.\n"); return 1; } + + /* Skip the next six lines */ + for ( i=0; i<6; i++ ) { + char tmp[1024]; + rval = fgets(tmp, 1024, fh); + if ( rval == NULL ) { + ERROR("Failed to read newmat file.\n"); + return 1; + } + } + + rval = fgets(symm, 32, fh); + if ( rval == NULL ) { + ERROR("Failed to read newmat file.\n"); + return 1; + } + fclose(fh); + if ( strncmp(symm, "SYMM ", 5) != 0 ) { + ERROR("Bad 'SYMM' line from MOSFLM.\n"); + return 1; + } + cen = symm[5]; + latt = spacegroup_to_lattice(symm+5); + /* 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); + cell_set_centering(cell, cen); + cell_set_lattice_type(cell, latt); + + if ( check_cell(mosflm->mp, image, cell) ) { + mosflm->success = 1; + mosflm->done = 1; + } - image->ncells = 1; + cell_free(cell); return 0; } @@ -319,7 +474,7 @@ static const char *spacegroup_for_lattice(UnitCell *cell) if ( centering != 'H' ) { g = "6"; } else { - g = "32"; + g = "3"; } break; @@ -343,8 +498,8 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) char tmp[256]; double wavelength; - switch ( mosflm->step ) { - + switch ( mosflm->step ) + { case 1 : mosflm_sendline("DETECTOR ROTATION HORIZONTAL" " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" @@ -352,13 +507,22 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) break; case 2 : - if ( mosflm->target_cell != NULL ) { + if ( (mosflm->mp->indm & INDEXING_USE_LATTICE_TYPE) + && (mosflm->mp->template != NULL) ) + { const char *symm; - symm = spacegroup_for_lattice(mosflm->target_cell); + + 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); mosflm_sendline(tmp, mosflm); + } else { - mosflm_sendline("SYMM P1\n", mosflm); + mosflm_sendline("\n", mosflm); } break; @@ -403,10 +567,9 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) mosflm->finished_ok = 1; break; - default: + default: mosflm_sendline("exit\n", mosflm); return; - } mosflm->step++; @@ -460,8 +623,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) switch ( type ) { - case MOSFLM_INPUT_LINE : - + case MOSFLM_INPUT_LINE : block_buffer = malloc(i+1); memcpy(block_buffer, mosflm->rbuffer, i); block_buffer[i] = '\0'; @@ -475,12 +637,12 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) endbit_length = i+2; break; - case MOSFLM_INPUT_PROMPT : + case MOSFLM_INPUT_PROMPT : mosflm_send_next(image, mosflm); endbit_length = i+7; break; - default : + default : /* Obviously, this never happens :) */ ERROR("Unrecognised MOSFLM input mode! " @@ -527,7 +689,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 +699,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 +716,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,7 +744,11 @@ 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; + rval = 0; do { fd_set fds; @@ -632,8 +796,50 @@ 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; + int need_cell = 0; + + if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; + if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; + if ( *indm & INDEXING_USE_LATTICE_TYPE ) need_cell = 1; + + if ( need_cell && (cell == NULL) ) { + ERROR("MOSFLM needs a unit cell for this set of flags.\n"); + return NULL; + } + + /* Flags that MOSFLM knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS + | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS + | INDEXING_USE_LATTICE_TYPE; + + mp = malloc(sizeof(struct mosflm_private)); + if ( mp == NULL ) return NULL; + + mp->ltl = ltl; + mp->template = cell; + mp->indm = *indm; + + return (IndexingPrivate *)mp; +} + + +void mosflm_cleanup(IndexingPrivate *pp) +{ + struct mosflm_private *p; + p = (struct mosflm_private *)pp; + free(p); } diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h index ba1eb73d..a87232b6 100644 --- a/libcrystfel/src/mosflm.h +++ b/libcrystfel/src/mosflm.h @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * 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. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian <rkirian@asu.edu> - * 2012 Thomas White <taw@physics.org> + * 2012-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -35,10 +35,16 @@ #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); +extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); + +extern void mosflm_cleanup(IndexingPrivate *pp); #endif /* MOSFLM_H */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index dabc3b1a..6c09053b 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -3,12 +3,12 @@ * * Peak search and other image analysis * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de> * 2011 Andrew Martin <andrew.martin@desy.de> * 2011 Richard Kirian @@ -147,24 +147,15 @@ static int cull_peaks(struct image *image) } -/* cfs, css relative to panel origin */ -static int *make_BgMask(struct image *image, struct panel *p, - double ir_out, double ir_inn) +static void add_crystal_to_mask(struct image *image, struct panel *p, + double ir_inn, int w, int h, + int *mask, Crystal *cr) { Reflection *refl; RefListIterator *iter; - int *mask; - int w, h; - - w = p->max_fs - p->min_fs + 1; - h = p->max_ss - p->min_ss + 1; - mask = calloc(w*h, sizeof(int)); - if ( mask == NULL ) return NULL; - - if ( image->reflections == NULL ) return mask; /* Loop over all reflections */ - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -205,6 +196,27 @@ static int *make_BgMask(struct image *image, struct panel *p, } } +} + + +/* cfs, css relative to panel origin */ +static int *make_BgMask(struct image *image, struct panel *p, double ir_inn) +{ + int *mask; + int w, h; + int i; + + w = p->max_fs - p->min_fs + 1; + h = p->max_ss - p->min_ss + 1; + mask = calloc(w*h, sizeof(int)); + if ( mask == NULL ) return NULL; + + if ( image->crystals == NULL ) return mask; + + for ( i=0; i<image->n_crystals; i++ ) { + add_crystal_to_mask(image, p, ir_inn, + w, h, mask, image->crystals[i]); + } return mask; } @@ -401,9 +413,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, int nrej_dis = 0; int nrej_pro = 0; int nrej_fra = 0; - int nrej_bad = 0; + int nrej_fail = 0; int nrej_snr = 0; - int nrej_sat = 0; int nacc = 0; int ncull; @@ -426,9 +437,6 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Overall threshold */ if ( data[fs+stride*ss] < threshold ) continue; - /* Immediate rejection of pixels above max_adu */ - if ( data[fs+stride*ss] > p->max_adu ) continue; - /* Get gradients */ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss]; dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss]; @@ -494,12 +502,11 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Centroid peak and get better coordinates. */ r = integrate_peak(image, mask_fs, mask_ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, &saturated); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); - if ( r ) { - /* Bad region - don't detect peak */ - nrej_bad++; - continue; + if ( saturated ) { + image->num_saturated_peaks++; + if ( !use_saturated ) continue; } /* It is possible for the centroid to fall outside the image */ @@ -521,8 +528,9 @@ static void search_peaks_in_panel(struct image *image, float threshold, continue; } - if ( saturated && !use_saturated ) { - nrej_sat++; + if ( r ) { + /* Bad region - don't detect peak */ + nrej_fail++; continue; } @@ -543,10 +551,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, ncull = 0; } + image->num_peaks += nacc; + //STATUS("%i accepted, %i box, %i proximity, %i outside panel, " - // "%i in bad regions, %i with SNR < %g, %i badrow culled, " + // "%i failed integration, %i with SNR < %g, %i badrow culled, " // "%i saturated.\n", - // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, + // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_fail, // nrej_snr, min_snr, ncull, nrej_sat); if ( ncull != 0 ) { @@ -567,6 +577,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient, image_feature_list_free(image->features); } image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; for ( i=0; i<image->det->n_panels; i++ ) { @@ -581,29 +593,19 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } -double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) +int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { - int i; int n_feat = 0; int n_sane = 0; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double min_dist = 0.25; - double stot = 0.0; - - /* Round towards nearest */ - fesetround(1); - - /* Cell basis vectors for this image */ - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + int i; + const double min_dist = 0.25; - /* Loop over peaks, checking proximity to nearest reflection */ for ( i=0; i<image_feature_count(image->features); i++ ) { struct imagefeature *f; struct rvec q; double h,k,l,hd,kd,ld; + int j; /* Assume all image "features" are genuine peaks */ f = image_get_feature(image->features, i); @@ -613,38 +615,42 @@ double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) /* Reciprocal space position of found peak */ q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); - /* Decimal and fractional Miller indices of nearest - * reciprocal lattice point */ - hd = q.u * ax + q.v * ay + q.w * az; - kd = q.u * bx + q.v * by + q.w * bz; - ld = q.u * cx + q.v * cy + q.w * cz; - h = lrint(hd); - k = lrint(kd); - l = lrint(ld); - - /* Check distance */ - if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist) - && (fabs(l - ld) < min_dist) ) - { - double sval; - n_sane++; - sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); - stot += 1.0 - sval; - continue; - } + 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(crystals[j]), + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + /* Decimal and fractional Miller indices of nearest + * reciprocal lattice point */ + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + /* Check distance */ + if ( (fabs(h - hd) < min_dist) + && (fabs(k - kd) < min_dist) + && (fabs(l - ld) < min_dist) ) + { + n_sane++; + continue; + } - } + } - *pst = stot; - return (double)n_sane / (float)n_feat; -} + } -int peak_sanity_check(struct image *image) -{ - double stot; /* 0 means failed test, 1 means passed test */ - return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5; + return ((double)n_sane / n_feat) >= 0.5; } @@ -705,43 +711,29 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell, } -/* Integrate the list of predicted reflections in "image" */ -void integrate_reflections(struct image *image, int use_closer, int bgsub, - double min_snr, - double ir_inn, double ir_mid, double ir_out, - int integrate_saturated) +static void integrate_crystal(Crystal *cr, struct image *image, int use_closer, + int bgsub, double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int **bgMasks, + int res_cutoff) { + RefList *reflections; struct integr_ind *il; int n, i; double av = 0.0; int first = 1; - int **bgMasks; + int n_saturated = 0; - if ( num_reflections(image->reflections) == 0 ) return; + reflections = crystal_get_reflections(cr); - il = sort_reflections(image->reflections, image->indexed_cell, &n); + if ( num_reflections(reflections) == 0 ) return; + + il = sort_reflections(reflections, crystal_get_cell(cr), &n); if ( il == NULL ) { ERROR("Couldn't sort reflections\n"); return; } - /* Make background masks for all panels */ - bgMasks = calloc(image->det->n_panels, sizeof(int *)); - if ( bgMasks == NULL ) { - ERROR("Couldn't create list of background masks.\n"); - return; - } - for ( i=0; i<image->det->n_panels; i++ ) { - int *mask; - mask = make_BgMask(image, &image->det->panels[i], - ir_out, ir_inn); - if ( mask == NULL ) { - ERROR("Couldn't create background mask.\n"); - return; - } - bgMasks[i] = mask; - } - for ( i=0; i<n; i++ ) { double fs, ss, intensity; @@ -809,7 +801,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, 1, bgMasks[pnum], &saturated); if ( saturated ) { - image->n_saturated++; + n_saturated++; if ( !integrate_saturated ) r = 1; } @@ -840,18 +832,52 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } //STATUS("%5.2f A, %5.2f, av %5.2f\n", // 1e10/il[i].res, snr, av); - //if ( av < 1.0 ) break; + if ( res_cutoff && (av < 1.0) ) break; + } + } + + crystal_set_num_saturated_reflections(cr, n_saturated); + crystal_set_resolution_limit(cr, 0.0); + + free(il); +} + + +/* Integrate the list of predicted reflections in "image" */ +void integrate_reflections(struct image *image, int use_closer, int bgsub, + double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int res_cutoff) +{ + int i; + int **bgMasks; + + /* Make background masks for all panels */ + bgMasks = calloc(image->det->n_panels, sizeof(int *)); + if ( bgMasks == NULL ) { + ERROR("Couldn't create list of background masks.\n"); + return; + } + for ( i=0; i<image->det->n_panels; i++ ) { + int *mask; + mask = make_BgMask(image, &image->det->panels[i], ir_inn); + if ( mask == NULL ) { + ERROR("Couldn't create background mask.\n"); + return; } + bgMasks[i] = mask; + } + + for ( i=0; i<image->n_crystals; i++ ) { + integrate_crystal(image->crystals[i], image, use_closer, + bgsub, min_snr, ir_inn, ir_mid, ir_out, + integrate_saturated, bgMasks, res_cutoff); } for ( i=0; i<image->det->n_panels; i++ ) { free(bgMasks[i]); } free(bgMasks); - - image->diffracting_resolution = 0.0; - - free(il); } @@ -894,12 +920,17 @@ void validate_peaks(struct image *image, double min_snr, r = integrate_peak(image, f->fs, f->ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, &saturated); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); if ( r ) { n_int++; continue; } + if ( saturated ) { + n_sat++; + if ( !use_saturated ) continue; + } + /* It is possible for the centroid to fall outside the image */ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs) || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) @@ -920,11 +951,6 @@ void validate_peaks(struct image *image, double min_snr, continue; } - if ( saturated && !use_saturated ) { - n_sat++; - continue; - } - /* Add using "better" coordinates */ image_add_feature(flist, f_fs, f_ss, image, intensity, NULL); @@ -936,4 +962,6 @@ void validate_peaks(struct image *image, double min_snr, // n_wtf, n_int, n_dft, n_snr, n_prx, n_sat); image_feature_list_free(image->features); image->features = flist; + image->num_saturated_peaks = n_sat; + image->num_peaks = image_feature_count(flist); } diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 3fae13fb..6be728fe 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -3,11 +3,11 @@ * * Peak search and other image analysis * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -45,12 +45,10 @@ extern void search_peaks(struct image *image, float threshold, extern void integrate_reflections(struct image *image, int use_closer, int bgsub, double min_snr, double ir_inn, double ir_mid, double ir_out, - int integrate_saturated); + int integrate_saturated, int res_cutoff); -extern double peak_lattice_agreement(struct image *image, UnitCell *cell, - double *pst); - -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.c b/libcrystfel/src/reax.c index 7529d12f..3460a4f2 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -3,11 +3,11 @@ * * A new auto-indexer * - * 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: - * 2011-2012 Thomas White <taw@physics.org> + * 2011-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -65,6 +65,9 @@ #define MAX_CANDIDATES (1024) +/* Choose the best solution from this many candidate cells */ +#define MAX_REAX_CELL_CANDIDATES (32) + struct dvec { double x; @@ -103,7 +106,7 @@ struct reax_search struct reax_private { - IndexingPrivate base; + IndexingMethod indm; struct dvec *directions; int n_dir; double angular_inc; @@ -907,7 +910,7 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, } - if ( cl->n_cand >= MAX_CELL_CANDIDATES ) { + if ( cl->n_cand >= MAX_REAX_CELL_CANDIDATES ) { /* "cshift" just fell off the end of the list */ } else { cl->cand[cl->n_cand++] = cshift; @@ -915,15 +918,46 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, } +static int check_cell(struct reax_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + out = cell_new_from_cell(cell); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + 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 assemble_cells_from_candidates(struct image *image, struct reax_search *s, - UnitCell *cell) + UnitCell *cell, + struct reax_private *p) { int i, j, k; signed int ti, tj, tk; struct cell_candidate_list cl; - cl.cand = calloc(MAX_CELL_CANDIDATES, sizeof(struct cell_candidate)); + cl.cand = calloc(MAX_REAX_CELL_CANDIDATES, + sizeof(struct cell_candidate)); if ( cl.cand == NULL ) { ERROR("Failed to allocate cell candidate list.\n"); return; @@ -967,7 +1001,9 @@ static void assemble_cells_from_candidates(struct image *image, continue; } - peak_lattice_agreement(image, cnew, &fom); + /* FIXME! */ + //peak_lattice_agreement(image, cnew, &fom); + fom = 1.0; add_cell_candidate(&cl, cnew, fom); } @@ -985,22 +1021,20 @@ static void assemble_cells_from_candidates(struct image *image, cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga); cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA, &alA, &beA, &gaA); - if ( (a - aA) > aA/10.0 ) w = 1; - if ( (b - bA) > bA/10.0 ) w = 1; - if ( (c - cA) > cA/10.0 ) w = 1; - if ( (al - alA) > deg2rad(5.0) ) w = 1; - if ( (be - beA) > deg2rad(5.0) ) w = 1; - if ( (ga - gaA) > deg2rad(5.0) ) w = 1; + if ( fabs(a - aA) > aA/10.0 ) w = 1; + if ( fabs(b - bA) > bA/10.0 ) w = 1; + if ( fabs(c - cA) > cA/10.0 ) w = 1; + if ( fabs(al - alA) > deg2rad(5.0) ) w = 1; + if ( fabs(be - beA) > deg2rad(5.0) ) w = 1; + if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1; if ( w ) { STATUS("This cell is a long way from that sought:\n"); cell_print(cl.cand[i].cell); } } - image->ncells = cl.n_cand; - assert(image->ncells <= MAX_CELL_CANDIDATES); for ( i=0; i<cl.n_cand; i++ ) { - image->candidate_cells[i] = cl.cand[i].cell; + if ( check_cell(p, image, cl.cand[i].cell) ) break; } free(cl.cand); @@ -1016,7 +1050,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) struct reax_search *s; int i; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; fft_in = fftw_malloc(p->nel*sizeof(double)); @@ -1039,7 +1072,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) // fft_in, fft_out, p->plan, smin, smax, // image->det, p); - assemble_cells_from_candidates(image, s, cell); + assemble_cells_from_candidates(image, s, cell, p); for ( i=0; i<s->n_search; i++ ) { free(s->search[i].cand); @@ -1051,16 +1084,27 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) } -IndexingPrivate *reax_prepare() +IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) { struct reax_private *p; int samp; double th; + if ( cell == NULL ) { + ERROR("ReAx needs a unit cell.\n"); + return NULL; + } + p = calloc(1, sizeof(*p)); if ( p == NULL ) return NULL; - p->base.indm = INDEXING_REAX; + /* Flags that ReAx knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS; + + /* Flags that ReAx requires */ + *indm |= INDEXING_USE_LATTICE_TYPE; p->angular_inc = deg2rad(1.0); @@ -1122,6 +1166,8 @@ IndexingPrivate *reax_prepare() p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, 1, FFTW_MEASURE); + p->indm = *indm; + return (IndexingPrivate *)p; } @@ -1130,7 +1176,6 @@ void reax_cleanup(IndexingPrivate *pp) { struct reax_private *p; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; free(p->directions); diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h index 01dfe649..657a2cdf 100644 --- a/libcrystfel/src/reax.h +++ b/libcrystfel/src/reax.h @@ -33,17 +33,26 @@ #include <config.h> #endif +#include "index.h" #include "cell.h" +#include "beam-parameters.h" +#include "detector.h" #ifdef HAVE_FFTW -extern IndexingPrivate *reax_prepare(void); +extern IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl); + 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(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) { return NULL; } @@ -52,7 +61,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) { } diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index b8ff9238..692a860a 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -3,12 +3,12 @@ * * Stream tools * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2011 Andrew Aquila * @@ -37,6 +37,10 @@ #include <stdlib.h> #include <stdio.h> #include <string.h> +#include <assert.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <fcntl.h> #include "cell.h" #include "utils.h" @@ -45,107 +49,17 @@ #include "reflist.h" #include "reflist-utils.h" +#define LATEST_MAJOR_VERSION (2) +#define LATEST_MINOR_VERSION (1) -#define CHUNK_START_MARKER "----- Begin chunk -----" -#define CHUNK_END_MARKER "----- End chunk -----" -#define PEAK_LIST_START_MARKER "Peaks from peak search" -#define PEAK_LIST_END_MARKER "End of peak list" -#define REFLECTION_START_MARKER "Reflections measured after indexing" -/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also - * used to terminate a standalone list of reflections */ -static void exclusive(const char *a, const char *b) +struct _stream { - ERROR("The stream options '%s' and '%s' are mutually exclusive.\n", - a, b); -} - - -int parse_stream_flags(const char *a) -{ - int n, i; - int ret = STREAM_NONE; - char **flags; - - n = assplode(a, ",", &flags, ASSPLODE_NONE); - - for ( i=0; i<n; i++ ) { - - if ( strcmp(flags[i], "integrated") == 0) { - - ret |= STREAM_INTEGRATED; - - } else if ( strcmp(flags[i], "peaks") == 0) { - if ( ret & STREAM_PEAKS_IF_INDEXED ) { - exclusive("peaks", "peaksifindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) { - exclusive("peaks", "peaksifnotindexed"); - return -1; - } - ret |= STREAM_PEAKS; - - } else if ( strcmp(flags[i], "peaksifindexed") == 0) { - if ( ret & STREAM_PEAKS ) { - exclusive("peaks", "peaksifindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) { - exclusive("peaksifnotindexed", - "peaksifindexed"); - return -1; - } - ret |= STREAM_PEAKS_IF_INDEXED; - - } else if ( strcmp(flags[i], "peaksifnotindexed") == 0) { - if ( ret & STREAM_PEAKS ) { - exclusive("peaks", "peaksifnotindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_INDEXED ) { - exclusive("peaksifnotindexed", - "peaksifindexed"); - return -1; - } - ret |= STREAM_PEAKS_IF_NOT_INDEXED; - - } else { - ERROR("Unrecognised stream flag '%s'\n", flags[i]); - return -1; - } - - free(flags[i]); - - } - free(flags); - - return ret; -} - - -int count_patterns(FILE *fh) -{ - char *rval; - - int n_total_patterns = 0; - do { - char line[1024]; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++; - - } while ( rval != NULL ); - - if ( ferror(fh) ) { - ERROR("Read error while counting patterns.\n"); - return 0; - } + FILE *fh; - return n_total_patterns; -} + int major_version; + int minor_version; +}; static int read_peaks(FILE *fh, struct image *image) @@ -215,96 +129,107 @@ static void write_peaks(struct image *image, FILE *ofh) } -void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f) +static void write_crystal(Stream *st, Crystal *cr, int include_reflections) { + UnitCell *cell; + RefList *reflist; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double a, b, c, al, be, ga; - fprintf(ofh, CHUNK_START_MARKER"\n"); + fprintf(st->fh, CRYSTAL_START_MARKER"\n"); - fprintf(ofh, "Image filename: %s\n", i->filename); + cell = crystal_get_cell(cr); + assert(cell != NULL); - if ( i->indexed_cell != NULL ) { + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + fprintf(st->fh, "Cell parameters %7.5f %7.5f %7.5f nm," + " %7.5f %7.5f %7.5f deg\n", + a*1.0e9, b*1.0e9, c*1.0e9, + rad2deg(al), rad2deg(be), rad2deg(ga)); - cell_get_parameters(i->indexed_cell, &a, &b, &c, - &al, &be, &ga); - fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm," - " %7.5f %7.5f %7.5f deg\n", - a*1.0e9, b*1.0e9, c*1.0e9, - rad2deg(al), rad2deg(be), rad2deg(ga)); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + fprintf(st->fh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", + asx/1e9, asy/1e9, asz/1e9); + fprintf(st->fh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + bsx/1e9, bsy/1e9, bsz/1e9); + fprintf(st->fh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + csx/1e9, csy/1e9, csz/1e9); - cell_get_reciprocal(i->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", - asx/1e9, asy/1e9, asz/1e9); - fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - bsx/1e9, bsy/1e9, bsz/1e9); - fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - csx/1e9, csy/1e9, csz/1e9); + reflist = crystal_get_reflections(cr); + if ( reflist != NULL ) { - } else { + fprintf(st->fh, "diffraction_resolution_limit" + " = %.2f nm^-1 or %.2f A\n", + crystal_get_resolution_limit(cr)/1e9, + 1e9 / crystal_get_resolution_limit(cr)); - fprintf(ofh, "No unit cell from indexing.\n"); + fprintf(st->fh, "num_saturated_reflections = %lli\n", + crystal_get_num_saturated_reflections(cr)); } - fprintf(ofh, "photon_energy_eV = %f\n", - J_to_eV(ph_lambda_to_en(i->lambda))); + if ( include_reflections ) { + + if ( reflist != NULL ) { - if ( i->reflections != NULL ) { + fprintf(st->fh, REFLECTION_START_MARKER"\n"); + write_reflections_to_file(st->fh, reflist); + fprintf(st->fh, REFLECTION_END_MARKER"\n"); - fprintf(ofh, "diffraction_resolution_limit" - " = %.2f nm^-1 or %.2f A\n", - i->diffracting_resolution/1e9, - 1e9 / i->diffracting_resolution); + } else { - fprintf(ofh, "num_saturated_reflections" - " = %lli\n", i->n_saturated); + fprintf(st->fh, "No integrated reflections.\n"); + } } + fprintf(st->fh, CRYSTAL_END_MARKER"\n"); +} + + +void write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, + int include_peaks, int include_reflections) +{ + int j; + + fprintf(st->fh, CHUNK_START_MARKER"\n"); + + fprintf(st->fh, "Image filename: %s\n", i->filename); + fprintf(st->fh, "indexed_by = %s\n", indexer_str(i->indexed_by)); + if ( i->det != NULL ) { int j; for ( j=0; j<i->det->n_panels; j++ ) { - fprintf(ofh, "camera_length_%s = %f\n", + fprintf(st->fh, "camera_length_%s = %f\n", i->det->panels[j].name, i->det->panels[j].clen); } } - copy_hdf5_fields(hdfile, i->copyme, ofh); + copy_hdf5_fields(hdfile, i->copyme, st->fh); - if ( (f & STREAM_PEAKS) - || ((f & STREAM_PEAKS_IF_INDEXED) && (i->indexed_cell != NULL)) - || ((f & STREAM_PEAKS_IF_NOT_INDEXED) && (i->indexed_cell == NULL)) ) - { - fprintf(ofh, "\n"); - write_peaks(i, ofh); + fprintf(st->fh, "num_peaks = %lli\n", i->num_peaks); + fprintf(st->fh, "num_saturated_peaks = %lli\n", i->num_saturated_peaks); + if ( include_peaks ) { + write_peaks(i, st->fh); } - if ( f & STREAM_INTEGRATED ) { - - fprintf(ofh, "\n"); - - if ( i->reflections != NULL ) { - - fprintf(ofh, REFLECTION_START_MARKER"\n"); - write_reflections_to_file(ofh, i->reflections); - fprintf(ofh, REFLECTION_END_MARKER"\n"); - - } else { - - fprintf(ofh, "No integrated reflections.\n"); + fprintf(st->fh, "photon_energy_eV = %f\n", + J_to_eV(ph_lambda_to_en(i->lambda))); - } + for ( j=0; j<i->n_crystals; j++ ) { + write_crystal(st, i->crystals[j], include_reflections); } - fprintf(ofh, CHUNK_END_MARKER"\n\n"); + fprintf(st->fh, CHUNK_END_MARKER"\n"); + + fflush(st->fh); } @@ -328,8 +253,7 @@ static int find_start_of_chunk(FILE *fh) } -/* Read the next chunk from a stream and fill in 'image' */ -int read_chunk(FILE *fh, struct image *image) +void read_crystal(Stream *st, struct image *image) { char line[1024]; char *rval = NULL; @@ -337,22 +261,116 @@ int read_chunk(FILE *fh, struct image *image) int have_as = 0; int have_bs = 0; int have_cs = 0; + Crystal *cr; + int n; + Crystal **crystals_new; + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal!\n"); + return; + } + + do { + + float u, v, w; + + rval = fgets(line, 1023, st->fh); + + /* Trouble? */ + if ( rval == NULL ) break; + + chomp(line); + if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) { + as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; + have_as = 1; + } + + if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) { + bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; + have_bs = 1; + } + + if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) { + cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; + have_cs = 1; + } + + if ( have_as && have_bs && have_cs ) { + + UnitCell *cell; + + cell = crystal_get_cell(cr); + + if ( cell != NULL ) { + ERROR("Duplicate cell found in stream!\n"); + ERROR("I'll use the most recent one.\n"); + cell_free(cell); + } + + cell = cell_new_from_reciprocal_axes(as, bs, cs); + crystal_set_cell(cr, cell); + + have_as = 0; have_bs = 0; have_cs = 0; + + } + + if ( strncmp(line, "num_saturated_reflections = ", 28) == 0 ) { + int n = atoi(line+28); + crystal_set_num_saturated_reflections(cr, n); + } + + if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { + + RefList *reflections; + + reflections = read_reflections_from_file(st->fh); + + if ( reflections == NULL ) { + ERROR("Failed while reading reflections\n"); + break; + } + + crystal_set_reflections(cr, reflections); + + } + + if ( strcmp(line, CRYSTAL_END_MARKER) == 0 ) break; + + } while ( 1 ); + + /* Add crystal to the list for this image */ + n = image->n_crystals+1; + crystals_new = realloc(image->crystals, n*sizeof(Crystal *)); + + if ( crystals_new == NULL ) { + ERROR("Failed to expand crystal list!\n"); + } else { + image->crystals = crystals_new; + image->crystals[image->n_crystals++] = cr; + } + +} + + +/* Read the next chunk from a stream and fill in 'image' */ +int read_chunk(Stream *st, struct image *image) +{ + char line[1024]; + char *rval = NULL; int have_filename = 0; int have_ev = 0; - if ( find_start_of_chunk(fh) ) return 1; + if ( find_start_of_chunk(st->fh) ) return 1; image->lambda = -1.0; image->features = NULL; - image->reflections = NULL; - image->indexed_cell = NULL; - image->n_saturated = 0; + image->crystals = NULL; + image->n_crystals = 0; do { - float u, v, w; - - rval = fgets(line, 1023, fh); + rval = fgets(line, 1023, st->fh); /* Trouble? */ if ( rval == NULL ) break; @@ -364,6 +382,11 @@ int read_chunk(FILE *fh, struct image *image) have_filename = 1; } + if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) { + image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19))); + have_ev = 1; + } + if ( strncmp(line, "camera_length_", 14) == 0 ) { if ( image->det != NULL ) { @@ -390,56 +413,19 @@ int read_chunk(FILE *fh, struct image *image) } } - if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) { - as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; - have_as = 1; - } - - if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) { - bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; - have_bs = 1; - } - - if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) { - cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; - have_cs = 1; - } - - if ( have_as && have_bs && have_cs ) { - if ( image->indexed_cell != NULL ) { - ERROR("Duplicate cell found in stream!\n"); - cell_free(image->indexed_cell); - } - image->indexed_cell = cell_new_from_reciprocal_axes(as, - bs, - cs); - have_as = 0; have_bs = 0; have_cs = 0; - } - - if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) { - image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19))); - have_ev = 1; - } - - if ( strncmp(line, "num_saturated_reflections = ", 28) == 0 ) { - image->n_saturated = atoi(line+28); - } - if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) { - if ( read_peaks(fh, image) ) { + if ( read_peaks(st->fh, image) ) { ERROR("Failed while reading peaks\n"); return 1; } } - if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { - image->reflections = read_reflections_from_file(fh); - if ( image->reflections == NULL ) { - ERROR("Failed while reading reflections\n"); - return 1; - } + if ( strcmp(line, CRYSTAL_START_MARKER) == 0 ) { + read_crystal(st, image); } + /* A chunk must have at least a filename and a wavelength, + * otherwise it's incomplete */ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) { if ( have_filename && have_ev ) return 0; ERROR("Incomplete chunk found in input file.\n"); @@ -448,6 +434,10 @@ int read_chunk(FILE *fh, struct image *image) } while ( 1 ); + if ( !feof(st->fh) ) { + ERROR("Error reading stream.\n"); + } + return 1; /* Either error or EOF, don't care because we will complain * on the terminal if it was an error. */ } @@ -457,7 +447,6 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) { int i; - fprintf(ofh, "CrystFEL stream format 2.0\n"); fprintf(ofh, "Command line:"); for ( i=0; i<argc; i++ ) { fprintf(ofh, " %s", argv[i]); @@ -467,25 +456,87 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) } -int skip_some_files(FILE *fh, int n) +Stream *open_stream_for_read(const char *filename) { - char *rval = NULL; - int n_patterns = 0; + Stream *st; - do { + st = malloc(sizeof(struct _stream)); + if ( st == NULL ) return NULL; - char line[1024]; + st->fh = fopen(filename, "r"); + if ( st->fh == NULL ) { + free(st); + return NULL; + } - if ( n_patterns == n ) return 0; + char line[1024]; + char *rval; - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++; + rval = fgets(line, 1023, st->fh); + if ( rval == NULL ) { + ERROR("Failed to read stream version.\n"); + close_stream(st); + return NULL; + } - } while ( rval != NULL ); + if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) { + st->major_version = 2; + st->minor_version = 0; + } else if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) { + st->major_version = 2; + st->minor_version = 1; + } else { + ERROR("Invalid stream, or stream format is too new.\n"); + close_stream(st); + return NULL; + } - return 1; + return st; +} + + +Stream *open_stream_fd_for_write(int fd) +{ + Stream *st; + + st = malloc(sizeof(struct _stream)); + if ( st == NULL ) return NULL; + + st->fh = fdopen(fd, "w"); + if ( st->fh == NULL ) { + free(st); + return NULL; + } + + st->major_version = LATEST_MAJOR_VERSION; + st->minor_version = LATEST_MINOR_VERSION; + + fprintf(st->fh, "CrystFEL stream format %i.%i\n", + st->major_version, st->minor_version); + + return st; +} + + + +Stream *open_stream_for_write(const char *filename) +{ + int fd; + + fd = open(filename, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR); + if ( fd == -1 ) { + ERROR("Failed to open stream.\n"); + return NULL; + } + + return open_stream_fd_for_write(fd); +} + + +void close_stream(Stream *st) +{ + fclose(st->fh); + free(st); } @@ -493,20 +544,53 @@ int is_stream(const char *filename) { FILE *fh; char line[1024]; - char *rval = NULL; + char *rval; + fh = fopen(filename, "r"); - if ( fh == NULL ) { - return -1; - } + if ( fh == NULL ) return 0; + rval = fgets(line, 1023, fh); fclose(fh); - if ( rval == NULL ) { - return -1; - } - if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) { - return 1; - } - else { - return 0; + if ( rval == NULL ) return 0; + + if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) return 1; + if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) return 1; + + return 0; +} + + +void write_line(Stream *st, const char *line) +{ + fprintf(st->fh, "%s\n", line); +} + + +void write_command(Stream *st, int argc, char *argv[]) +{ + int i; + + for ( i=0; i<argc; i++ ) { + if ( i > 0 ) fprintf(st->fh, " "); + fprintf(st->fh, "%s", argv[i]); } + fprintf(st->fh, "\n"); +} + + +/** + * rewind_stream: + * @st: A %Stream + * + * Attempts to set the file pointer for @st to the start of the stream, so that + * later calls to read_chunk() will repeat the sequence of chunks from the + * start. + * + * Programs must not assume that this operation always succeeds! + * + * Returns: non-zero if the stream could not be rewound. + */ +int rewind_stream(Stream *st) +{ + return fseek(st->fh, 0, SEEK_SET); } diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index a544798f..08e680bc 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -3,11 +3,11 @@ * * Stream tools * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Andrew Aquila * * This file is part of CrystFEL. @@ -38,30 +38,31 @@ struct image; struct hdfile; -/* Possible options dictating what goes into the output stream */ -enum -{ - STREAM_NONE = 0, - STREAM_INTEGRATED = 1<<0, - STREAM_PEAKS = 1<<2, - STREAM_PEAKS_IF_INDEXED = 1<<3, - STREAM_PEAKS_IF_NOT_INDEXED = 1<<4, -}; +#define CHUNK_START_MARKER "----- Begin chunk -----" +#define CHUNK_END_MARKER "----- End chunk -----" +#define PEAK_LIST_START_MARKER "Peaks from peak search" +#define PEAK_LIST_END_MARKER "End of peak list" +#define CRYSTAL_START_MARKER "--- Begin crystal" +#define CRYSTAL_END_MARKER "--- End crystal" +#define REFLECTION_START_MARKER "Reflections measured after indexing" +/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also + * used to terminate a standalone list of reflections */ +typedef struct _stream Stream; -extern int count_patterns(FILE *fh); +extern Stream *open_stream_for_read(const char *filename); +extern Stream *open_stream_for_write(const char *filename); +extern Stream *open_stream_fd_for_write(int fd); +extern void close_stream(Stream *st); -extern void write_stream_header(FILE *ofh, int argc, char *argv[]); +extern int read_chunk(Stream *st, struct image *image); +extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, + int include_peaks, int include_reflections); -extern void write_chunk(FILE *ofh, struct image *image, struct hdfile *hdfile, - int flags); - -extern int parse_stream_flags(const char *a); - -extern int read_chunk(FILE *fh, struct image *image); - -extern int skip_some_files(FILE *fh, int n); +extern void write_line(Stream *st, const char *line); +extern void write_command(Stream *st, int argc, char *argv[]); +extern int rewind_stream(Stream *st); extern int is_stream(const char *filename); #endif /* STREAM_H */ diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index a652e413..3c2011d1 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -1069,11 +1069,12 @@ static void do_op(const IntegerMatrix *op, * This function applies the @idx-th symmetry operation from @ops to the * reflection @h, @k, @l, and stores the result at @he, @ke and @le. * - * If you don't mind that the same equivalent might appear twice, simply call - * this function the number of times returned by num_ops(), using the actual - * point group. If repeating the same equivalent twice (for example, if the - * given reflection is a special high-symmetry one), call special_position() - * first to get a "specialised" SymOpList and use that instead. + * Call this function multiple times with idx=0 .. num_equivs(ops, m) to get all + * of the equivalent reflections in turn. + * + * If you don't mind that the same equivalent might appear twice, simply let + * @m = NULL. Otherwise, call new_symopmask() and then special_position() to + * set up a %SymOpMask appropriately. **/ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx, signed int h, signed int k, signed int l, @@ -1133,8 +1134,8 @@ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx, * @k: index of a reflection * @l: index of a reflection * - * This function determines which operations in @ops map the reflection @h, @k, - * @l onto itself, and uses @m to flag the operations in @ops which cause this. + * This function sets up @m to contain information about which operations in + * @ops map the reflection @h, @k, @l onto itself. * **/ void special_position(const SymOpList *ops, SymOpMask *m, @@ -1467,7 +1468,7 @@ static SymOpList *flack_reorder(const SymOpList *source) * operators, which are the symmetry operations which can be added to @target * to generate @source. Only rotations are allowable - no mirrors nor * inversions. - * To count the number of possibilities, use num_ops() on the result. + * To count the number of possibilities, use num_equivs() on the result. * * The algorithm used is "Algorithm A" from Flack (1987), Acta Cryst A43 p564. * |