diff options
Diffstat (limited to 'libcrystfel/src')
31 files changed, 2485 insertions, 752 deletions
diff --git a/libcrystfel/src/beam-parameters.c b/libcrystfel/src/beam-parameters.c index 4af25261..617f9061 100644 --- a/libcrystfel/src/beam-parameters.c +++ b/libcrystfel/src/beam-parameters.c @@ -8,6 +8,7 @@ * * Authors: * 2010,2012 Thomas White <taw@physics.org> + * 2012 Chunhong Yoon * * This file is part of CrystFEL. * @@ -35,6 +36,7 @@ #include "beam-parameters.h" #include "utils.h" +#include "hdf5-file.h" struct beam_params *get_beam_parameters(const char *filename) @@ -86,7 +88,12 @@ struct beam_params *get_beam_parameters(const char *filename) } else if ( strcmp(bits[0], "beam/radius") == 0 ) { b->beam_radius = atof(bits[2]); } else if ( strcmp(bits[0], "beam/photon_energy") == 0 ) { - b->photon_energy = atof(bits[2]); + if ( strncmp(bits[2], "/", 1) == 0 ) { + b->photon_energy = 0; // 0 means special case + b->photon_energy_from = strdup(bits[2]); + } else { + b->photon_energy = atof(bits[2]); + } } else if ( strcmp(bits[0], "beam/bandwidth") == 0 ) { b->bandwidth = atof(bits[2]); } else if ( strcmp(bits[0], "beam/divergence") == 0 ) { @@ -113,7 +120,7 @@ struct beam_params *get_beam_parameters(const char *filename) ERROR("Invalid or unspecified value for 'beam/radius'.\n"); reject = 1; } - if ( b->photon_energy < 0.0 ) { + if ( b->photon_energy < 0.0 ) { // 0 is ok ERROR("Invalid or unspecified value for" " 'beam/photon_energy'.\n"); reject = 1; @@ -141,3 +148,18 @@ struct beam_params *get_beam_parameters(const char *filename) return b; } + + +void free_beam_parameters(struct beam_params *beam) +{ + free(beam->photon_energy_from); + free(beam); +} + + +void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f) +{ + if ( beam->photon_energy_from != NULL ) { + beam->photon_energy = get_value(f, beam->photon_energy_from); + } +} diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h index 133e041b..de777deb 100644 --- a/libcrystfel/src/beam-parameters.h +++ b/libcrystfel/src/beam-parameters.h @@ -3,11 +3,12 @@ * * 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. * @@ -33,12 +34,16 @@ #include <config.h> #endif +struct beam_params; + +#include "hdf5-file.h" struct beam_params { double fluence; /* photons per pulse */ double beam_radius; /* metres */ double photon_energy; /* eV per photon */ + char *photon_energy_from; /* HDF5 dataset name */ double bandwidth; /* FWHM(wavelength) over wavelength. * Note: current simulation code just uses * a rectangular distribution with this as @@ -50,6 +55,8 @@ struct beam_params extern struct beam_params *get_beam_parameters(const char *filename); +extern void free_beam_parameters(struct beam_params *beam); +extern void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f); #endif /* BEAM_PARAMETERS_H */ 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/cell-utils.h b/libcrystfel/src/cell-utils.h index f92ab22d..8eb08ca7 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -55,7 +55,7 @@ extern int cell_is_sensible(UnitCell *cell); extern int validate_cell(UnitCell *cell); -extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr); +extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t); extern int bravais_lattice(UnitCell *cell); diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index cde25b07..57741b41 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -40,8 +40,15 @@ #include "utils.h" #include "integer_matrix.h" -/* A 3D vector in reciprocal space. - * Note: Heavily abused to serve as a real space vector as well */ +/** + * rvec: + * @u: x component (in reciprocal space) + * @v: y component (in reciprocal space) + * @w: z component (in reciprocal space) + * + * Structure representing a 3D vector in reciprocal space. + * Note: Heavily abused to serve as a real space vector as well. + **/ struct rvec { double u; @@ -49,6 +56,20 @@ struct rvec double w; }; + +/** + * LatticeType: + * @L_TRICLINIC: Triclinic lattice + * @L_MONOCLINIC: Monoclinic lattice + * @L_ORTHORHOMBIC: Orthorhombic lattice + * @L_TETRAGONAL: Tetragonal lattice + * @L_RHOMBOHEDRAL: Rhombohedral lattice + * @L_HEXAGONAL: Hexagonal lattice + * @L_CUBIC: Cubic lattice + * + * An enumeration of the possible lattice types: triclinic, monoclinic, + * orthorhombic, tetragonal, rhombohedral, hexagonal and cubic. + **/ typedef enum { L_TRICLINIC, 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/detector.c b/libcrystfel/src/detector.c index 9cea9a6d..41ebcc9a 100644 --- a/libcrystfel/src/detector.c +++ b/libcrystfel/src/detector.c @@ -462,8 +462,6 @@ void fill_in_values(struct detector *det, struct hdfile *f) if ( p->clen_from != NULL ) { p->clen = get_value(f, p->clen_from) * 1.0e-3; - free(p->clen_from); - p->clen_from = NULL; } p->clen += p->coffset; @@ -693,6 +691,64 @@ static void parse_toplevel(struct detector *det, const char *key, } +/* Test if fs,ss in panel "p" is further {out,in} than {*p_max_d,*p_min_d}, and + * if so update det->furthest_{out,in}_{panel,fs,ss}. */ +static void check_point(struct panel *p, double fs, double ss, + double *p_min_d, double *p_max_d, struct detector *det) +{ + double xs, ys, rx, ry, d; + + xs = fs*p->fsx + ss*p->ssx; + ys = fs*p->fsy + ss*p->ssy; + + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + + d = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + if ( d > *p_max_d ) { + + det->furthest_out_panel = p; + det->furthest_out_fs = fs; + det->furthest_out_ss = ss; + *p_max_d = d; + + } else if ( d < *p_min_d ) { + + det->furthest_in_panel = p; + det->furthest_in_fs = fs; + det->furthest_in_ss = ss; + *p_min_d = d; + + } +} + + +static void find_min_max_d(struct detector *det) +{ + double max_d, min_d; + int i; + + min_d = +INFINITY; + max_d = 0.0; + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p; + double w, h; + + p = &det->panels[i]; + w = p->max_fs - p->min_fs + 1; + h = p->max_ss - p->min_ss + 1; + + check_point(p, 0, 0, &min_d, &max_d, det); + check_point(p, w, 0, &min_d, &max_d, det); + check_point(p, 0, h, &min_d, &max_d, det); + check_point(p, w, h, &min_d, &max_d, det); + + } +} + + struct detector *get_detector_geometry(const char *filename) { FILE *fh; @@ -949,6 +1005,8 @@ out: } + find_min_max_d(det); + if ( reject ) return NULL; fclose(fh); @@ -1071,6 +1129,8 @@ struct detector *simple_geometry(const struct image *image) geom->panels[0].yfs = 0; geom->panels[0].yss = 1; + find_min_max_d(geom); + return geom; } @@ -1129,52 +1189,29 @@ static void check_extents(struct panel p, double *min_x, double *min_y, double largest_q(struct image *image) { - int fs, ss; - double ttm = 0.0; - double qmax = 0.0; - - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; ss++ ) { - - struct rvec q; - double tt; - - q = get_q(image, fs, ss, &tt, 1.0/image->lambda); - - if ( tt > ttm ) { - qmax = modulus(q.u, q.v, q.w); - ttm = tt; - } + struct rvec q; + double tt; - } - } + q = get_q_for_panel(image->det->furthest_out_panel, + image->det->furthest_out_fs, + image->det->furthest_out_ss, + &tt, 1.0/image->lambda); - return qmax; + return modulus(q.u, q.v, q.w); } double smallest_q(struct image *image) { - int fs, ss; - double ttm = +INFINITY; - double qmin = +INFINITY; - for ( fs=0; fs<image->width; fs++ ) { - for ( ss=0; ss<image->height; ss++ ) { - - struct rvec q; - double tt; - - q = get_q(image, fs, ss, &tt, 1.0/image->lambda); - - if ( tt < ttm ) { - qmin = modulus(q.u, q.v, q.w); - ttm = tt; - } + struct rvec q; + double tt; - } - } + q = get_q_for_panel(image->det->furthest_in_panel, + image->det->furthest_in_fs, + image->det->furthest_in_ss, + &tt, 1.0/image->lambda); - return qmin; + return modulus(q.u, q.v, q.w); } diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h index 1cd64716..43bdc42d 100644 --- a/libcrystfel/src/detector.h +++ b/libcrystfel/src/detector.h @@ -106,6 +106,18 @@ struct detector char **rigid_groups; int num_rigid_groups; + /* Location of the pixel furthest away from the beam position, which + * will have the largest value of 2theta regardless of camera length + * and wavelength */ + struct panel *furthest_out_panel; + double furthest_out_fs; + double furthest_out_ss; + + /* As above, but for the smallest 2theta */ + struct panel *furthest_in_panel; + double furthest_in_fs; + double furthest_in_ss; + struct panel defaults; }; diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 43de6c6e..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\n"); - return; + ERROR("Failed to fork for DirAx: %s\n", strerror(errno)); + 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 { @@ -521,9 +588,21 @@ void run_dirax(struct image *image) sval = select(dirax->pty+1, &fds, NULL, NULL, &tv); if ( sval == -1 ) { - int err = errno; - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; + + const int err = errno; + + switch ( err ) { + + case EINTR: + STATUS("Restarting select()\n"); + break; + + default: + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + + } + } else if ( sval != 0 ) { rval = dirax_readable(image, dirax); } else { @@ -531,7 +610,7 @@ void run_dirax(struct image *image) rval = 1; } - } while ( !rval ); + } while ( !rval && !dirax->success ); close(dirax->pty); free(dirax->rbuffer); @@ -541,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/grainspotter.c b/libcrystfel/src/grainspotter.c new file mode 100644 index 00000000..8bdb03aa --- /dev/null +++ b/libcrystfel/src/grainspotter.c @@ -0,0 +1,503 @@ +/* + * grainspotter.c + * + * Invoke GrainSpotter for multi-crystal autoindexing + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White <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 <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <string.h> +#include <unistd.h> +#include <sys/wait.h> +#include <fcntl.h> +#include <assert.h> +#include <sys/ioctl.h> +#include <errno.h> + +#if HAVE_FORKPTY_LINUX +#include <pty.h> +#elif HAVE_FORKPTY_BSD +#include <util.h> +#endif + + +#include "image.h" +#include "utils.h" +#include "peaks.h" +#include "cell.h" +#include "cell-utils.h" +#include "grainspotter.h" + + +#define GRAINSPOTTER_VERBOSE 0 + + +/* Global private data, prepared once */ +struct grainspotter_private +{ + IndexingMethod indm; + UnitCell *cell; +}; + + +/* Data needed during the run of Grainspotter */ +struct grainspotter_data { + + struct grainspotter_private *gp; + + /* Low-level stuff */ + int pty; + pid_t pid; + char *rbuffer; + int rbufpos; + int rbuflen; + +}; + + +static int read_matrix(struct grainspotter_private *gp, struct image *image, + char *filename) +{ + FILE *fh; + int d1; + float d2; + float ubi11, ubi12, ubi13; + float ubi21, ubi22, ubi23; + float ubi31, ubi32, ubi33; + char line[1024]; + int r; + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Can't open '%s'\n", filename); + return 1; + } + + /* Read and discard first line */ + if ( fgets(line, 1024, fh) == NULL ) { + ERROR("Failed to read GFF file.\n"); + return 1; + } + + do { + + Crystal *cr; + UnitCell *cell; + + /* One line per grain */ + if ( fgets(line, 1024, fh) == NULL ) { + ERROR("Failed to read GFF file.\n"); + return 1; + } + + STATUS("'%s'\n", line); + + r = sscanf(line, "%i %f %f %f %f %f %f %f %f %f %f %f %f" + "%f %f %f %f %f %f %f %f %f %f %f", + &d1, &d2, &d2, &d2, &d2, &d2, &d2, &d2, &d2, + &d2, &d2, &d2, &d2, &d2, &d2, + &ubi11, &ubi12, &ubi13, + &ubi21, &ubi22, &ubi23, + &ubi31, &ubi32, &ubi33); + + if ( r != 24 ) { + ERROR("Only %i parameters in GFF file\n", r); + return 1; + } + + cell = cell_new(); + + cell_set_cartesian(cell, ubi12/1e10, ubi13/1e10, ubi11/1e10, + ubi22/1e10, ubi23/1e10, ubi21/1e10, + ubi32/1e10, ubi33/1e10, ubi31/1e10); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, cell); + image_add_crystal(image, cr); + + } while ( !feof(fh) ); + + fclose(fh); + + if ( gp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, image->crystals, + image->n_crystals) ) + { + free_all_crystals(image); + return 0; + } + } + + return 0; +} + + +static void gs_parseline(char *line, struct image *image, + struct grainspotter_data *gs) +{ + #if GRAINSPOTTER_VERBOSE + STATUS("%s\n", line); + #endif +} + + +static int grainspotter_readable(struct image *image, + struct grainspotter_data *gs) +{ + int rval; + int no_string = 0; + + rval = read(gs->pty, gs->rbuffer+gs->rbufpos, gs->rbuflen-gs->rbufpos); + + if ( (rval == -1) || (rval == 0) ) return 1; + + gs->rbufpos += rval; + assert(gs->rbufpos <= gs->rbuflen); + + while ( (!no_string) && (gs->rbufpos > 0) ) { + + int i; + int block_ready = 0; + + /* See if there's a full line in the buffer yet */ + for ( i=0; i<gs->rbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + + if ( (gs->rbuffer[i] == '\r') + && (gs->rbuffer[i+1] == '\n') ) { + block_ready = 1; + break; + } + + } + + if ( block_ready ) { + + unsigned int new_rbuflen; + unsigned int endbit_length; + char *block_buffer = NULL; + + block_buffer = malloc(i+1); + memcpy(block_buffer, gs->rbuffer, i); + block_buffer[i] = '\0'; + + if ( block_buffer[0] == '\r' ) { + memmove(block_buffer, block_buffer+1, i); + } + + gs_parseline(block_buffer, image, gs); + free(block_buffer); + endbit_length = i+2; + + /* Now the block's been parsed, it should be + * forgotten about */ + memmove(gs->rbuffer, + gs->rbuffer + endbit_length, + gs->rbuflen - endbit_length); + + /* Subtract the number of bytes removed */ + gs->rbufpos = gs->rbufpos - endbit_length; + new_rbuflen = gs->rbuflen - endbit_length; + if ( new_rbuflen == 0 ) new_rbuflen = 256; + gs->rbuffer = realloc(gs->rbuffer, new_rbuflen); + gs->rbuflen = new_rbuflen; + + } else { + + if ( gs->rbufpos==gs->rbuflen ) { + + /* More buffer space is needed */ + gs->rbuffer = realloc(gs->rbuffer, + gs->rbuflen + 256); + gs->rbuflen = gs->rbuflen + 256; + /* The new space gets used at the next + * read, shortly... */ + + } + no_string = 1; + + } + + } + + return 0; +} + + +static void write_gve(struct image *image, struct grainspotter_private *gp) +{ + FILE *fh; + int i; + char filename[1024]; + double a, b, c, al, be, ga; + + snprintf(filename, 1023, "xfel-%i.gve", image->id); + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + return; + } + + cell_get_parameters(gp->cell, &a, &b, &c, &al, &be, &ga); + fprintf(fh, "%.6f %.6f %.6f %.6f %.6f %.6f P\n", a*1e10, b*1e10, c*1e10, + rad2deg(al), rad2deg(be), rad2deg(ga)); + fprintf(fh, "# wavelength = %.6f\n", image->lambda*1e10); + fprintf(fh, "# wedge = 0.000000\n"); + fprintf(fh, "# ds h k l\n"); + fprintf(fh, "# xr yr zr xc yc ds eta omega\n"); + + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + fprintf(fh, "%.6f %.6f %.6f 0 0 %.6f %.6f %.6f 0\n", + f->rz/1e10, f->rx/1e10, f->ry/1e10, + modulus(f->rx, f->ry, f->rz)/1e10, /* dstar */ + atan2(f->ry, f->rx), 0.0); /* eta, omega */ + + } + fclose(fh); +} + + +static char *write_ini(struct image *image) +{ + FILE *fh; + char *filename; + double tt; + + filename = malloc(1024); + if ( filename == NULL ) return NULL; + + snprintf(filename, 1023, "xfel-%i.ini", image->id); + + fh = fopen(filename, "w"); + if ( !fh ) { + ERROR("Couldn't open temporary file '%s'\n", filename); + free(filename); + return NULL; + } + + get_q_for_panel(image->det->furthest_out_panel, + image->det->furthest_out_fs, + image->det->furthest_out_ss, + &tt, 1.0/image->lambda); + + fprintf(fh, "spacegroup 96\n"); + fprintf(fh, "!dsrange 0 1.3\n"); + fprintf(fh, "tthrange 0 %.2f\n", 20.0); + fprintf(fh, "etarange 0 360\n"); + fprintf(fh, "domega 1\n"); + fprintf(fh, "omegarange -0.5 0.5\n"); + fprintf(fh, "filespecs xfel-%i.gve xfel-%i.log\n", + image->id, image->id); + fprintf(fh, "cuts 3 0.1 0.5\n"); + fprintf(fh, "eulerstep 8\n"); + fprintf(fh, "uncertainties 0.1 0.2 .5\n"); + fprintf(fh, "nsigmas 2\n"); + fprintf(fh, "minfracg 0.95\n"); + fprintf(fh, "Nhkls_in_indexing 500\n"); + fprintf(fh, "random 10000\n"); + fprintf(fh, "!positionfit\n"); + fprintf(fh, "genhkl\n"); + fprintf(fh, "xfelmode\n"); + + fclose(fh); + + return filename; +} + + +int grainspotter_index(struct image *image, IndexingPrivate *ipriv) +{ + unsigned int opts; + int status; + int rval; + struct grainspotter_data *grainspotter; + struct grainspotter_private *gp = (struct grainspotter_private *)ipriv; + char *ini_filename; + char gff_filename[1024]; + + write_gve(image, gp); + ini_filename = write_ini(image); + + if ( ini_filename == NULL ) { + ERROR("Failed to write ini file for GrainSpotter.\n"); + return 0; + } + + grainspotter = malloc(sizeof(struct grainspotter_data)); + if ( grainspotter == NULL ) { + ERROR("Couldn't allocate memory for GrainSpotter data.\n"); + return 0; + } + + grainspotter->gp = gp; + + snprintf(gff_filename, 1023, "xfel-%i.gff", image->id); + remove(gff_filename); + + grainspotter->pid = forkpty(&grainspotter->pty, NULL, NULL, NULL); + if ( grainspotter->pid == -1 ) { + ERROR("Failed to fork for GrainSpotter: %s\n", strerror(errno)); + return 0; + } + if ( grainspotter->pid == 0 ) { + + /* Child process: invoke GrainSpotter */ + struct termios t; + + /* Turn echo off */ + tcgetattr(STDIN_FILENO, &t); + t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); + tcsetattr(STDIN_FILENO, TCSANOW, &t); + + STATUS("Running GrainSpotter.0.93 '%s'\n", ini_filename); + execlp("GrainSpotter.0.93", "", ini_filename, (char *)NULL); + ERROR("Failed to invoke GrainSpotter.\n"); + _exit(0); + + } + free(ini_filename); + + grainspotter->rbuffer = malloc(256); + grainspotter->rbuflen = 256; + grainspotter->rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(grainspotter->pty, F_GETFL); + fcntl(grainspotter->pty, F_SETFL, opts | O_NONBLOCK); + + do { + + fd_set fds; + struct timeval tv; + int sval; + + FD_ZERO(&fds); + FD_SET(grainspotter->pty, &fds); + + tv.tv_sec = 30; + tv.tv_usec = 0; + + sval = select(grainspotter->pty+1, &fds, NULL, NULL, &tv); + + if ( sval == -1 ) { + + const int err = errno; + + switch ( err ) { + + case EINTR: + STATUS("Restarting select()\n"); + break; + + default: + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + + } + + } else if ( sval != 0 ) { + rval = grainspotter_readable(image, grainspotter); + } else { + ERROR("No response from GrainSpotter..\n"); + rval = 1; + } + + } while ( !rval ); + + close(grainspotter->pty); + free(grainspotter->rbuffer); + waitpid(grainspotter->pid, &status, 0); + + if ( status != 0 ) { + ERROR("GrainSpotter doesn't seem to be working properly.\n"); + free(grainspotter); + return 0; + } + + if ( read_matrix(gp, image, gff_filename) != 0 ) { + free(grainspotter); + return 0; + } + + /* Success! */ + free(grainspotter); + return 1; +} + + + +IndexingPrivate *grainspotter_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl) +{ + struct grainspotter_private *gp; + + if ( cell == NULL ) { + ERROR("GrainSpotter needs a unit cell.\n"); + return NULL; + } + + gp = calloc(1, sizeof(*gp)); + if ( gp == NULL ) return NULL; + + /* Flags that GrainSpotter knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS; + + /* Flags that GrainSpotter requires */ + *indm |= INDEXING_USE_LATTICE_TYPE; + + gp->cell = cell; + gp->indm = *indm; + + return (IndexingPrivate *)gp; +} + + +void grainspotter_cleanup(IndexingPrivate *pp) +{ + struct grainspotter_private *p; + + p = (struct grainspotter_private *)pp; + free(p); +} diff --git a/libcrystfel/src/grainspotter.h b/libcrystfel/src/grainspotter.h new file mode 100644 index 00000000..720fc486 --- /dev/null +++ b/libcrystfel/src/grainspotter.h @@ -0,0 +1,50 @@ +/* + * grainspotter.h + * + * Invoke GrainSpotter for multi-crystal autoindexing + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2013 Thomas White <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 GRAINSPOTTER_H +#define GRAINSPOTTER_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "cell.h" + +extern IndexingPrivate *grainspotter_prepare(IndexingMethod *indm, + UnitCell *cell, + const char *filename, + struct detector *det, + struct beam_params *beam, + float *ltl); + +extern void grainspotter_cleanup(IndexingPrivate *pp); + +extern int grainspotter_index(struct image *image, IndexingPrivate *p); + + +#endif /* GRAINSPOTTER_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 093e20f2..8fe3d69c 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -159,3 +159,31 @@ 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; +} + + +void free_all_crystals(struct image *image) +{ + int i; + + for ( i=0; i<image->n_crystals; i++ ) { + crystal_free(image->crystals[i]); + } +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index afa9e4a7..0d131451 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,7 @@ 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); +extern void free_all_crystals(struct image *image); + #endif /* IMAGE_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index c8a6b08f..20e7cf24 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 * 2013 Cornelius Gati <cornelius.gati@cfel.de> @@ -50,53 +50,41 @@ #include "index.h" #include "index-priv.h" #include "reax.h" +#include "grainspotter.h" #include "geometry.h" #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"; -} +#include "grainspotter.h" 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; + char *str; - 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_XDS : @@ -104,11 +92,45 @@ 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; + + case INDEXING_GRAINSPOTTER : + iprivs[n] = grainspotter_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; + + str = indexer_str(indm[n]); + STATUS("Prepared indexing method %i: %s\n", n, str); + free(str); + + 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; @@ -116,26 +138,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_XDS : @@ -143,7 +165,16 @@ void cleanup_indexing(IndexingPrivate **priv) break; case INDEXING_REAX : - reax_cleanup(priv[n]); + reax_cleanup(privs[n]); + break; + + case INDEXING_GRAINSPOTTER : + grainspotter_cleanup(privs[n]); + break; + + default : + ERROR("Don't know how to clean up indexing method %i\n", + indms[n]); break; } @@ -176,146 +207,229 @@ 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_XDS : + run_XDS(image, cell); + break; + + case INDEXING_REAX : + return reax_index(image, ipriv); + break; + + case INDEXING_GRAINSPOTTER : + return grainspotter_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_XDS : - run_XDS(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; + + case INDEXING_GRAINSPOTTER : + strcpy(str, "grainspotter"); + break; + + default : + ERROR("Unrecognised indexing method %i\n", + indm & INDEXING_METHOD_MASK); + 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[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; + } else if ( strcmp(methods[i], "xds") == 0) { list[i] = INDEXING_XDS; + } 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; } @@ -323,7 +437,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 b2d7553d..d550468b 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,48 +38,79 @@ #endif +#include "beam-parameters.h" #include "cell.h" #include "image.h" #include "detector.h" -/* Indexing methods */ +#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) + +#define INDEXING_DEFAULTS_GRAINSPOTTER (INDEXING_GRAINSPOTTER \ + | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_CHECK_PEAKS) + +/** + * IndexingMethod: + * @INDEXING_NONE: No indexing to be performed + * @INDEXING_DIRAX: Invoke DirAx + * @INDEXING_MOSFLM: Invoke MOSFLM + * @INDEXING_REAX: DPS algorithm using known cell parameters + * + * An enumeration of all the available indexing methods. + **/ typedef enum { - INDEXING_NONE, - INDEXING_DIRAX, - INDEXING_MOSFLM, - INDEXING_XDS, - INDEXING_REAX, -} IndexingMethod; + INDEXING_NONE = 0, -/* Cell reduction methods */ -enum { - CELLR_NONE, - CELLR_REDUCE, - CELLR_COMPARE, - CELLR_COMPARE_AB, -}; + /* The core indexing methods themselves */ + INDEXING_DIRAX = 1, + INDEXING_MOSFLM = 2, + INDEXING_REAX = 3, + INDEXING_GRAINSPOTTER = 4, + INDEXING_XDS = 5, + /* Bits at the top of the IndexingMethod are flags which modify the + * behaviour of the indexer. */ + INDEXING_CHECK_CELL_COMBINATIONS = 256, + INDEXING_CHECK_CELL_AXES = 512, + INDEXING_CHECK_PEAKS = 1024, + INDEXING_USE_LATTICE_TYPE = 2048, -typedef struct _indexingprivate IndexingPrivate; +} IndexingMethod; -extern IndexingPrivate *indexing_private(IndexingMethod indm); +/* This defines the bits in "IndexingMethod" which are used to represent the + * core of the indexing method */ +#define INDEXING_METHOD_MASK (0xff) -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); +/** + * IndexingPrivate: + * + * This is an opaque data structure containing information needed by the + * indexing method. + **/ +typedef void *IndexingPrivate; -extern void index_pattern(struct image *image, UnitCell *cell, - IndexingMethod *indm, int cellr, int verbose, - IndexingPrivate **priv, int config_insane, - const float *ltl); +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, + 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/integer_matrix.c b/libcrystfel/src/integer_matrix.c index 32784f20..96159cc6 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -62,6 +62,8 @@ struct _integermatrix /** * intmat_new: + * @rows: Number of rows that the new matrix is to have + * @cols: Number of columns that the new matrix is to have * * Allocates a new %IntegerMatrix with all elements set to zero. * @@ -406,7 +408,7 @@ void intmat_print(const IntegerMatrix *m) * intmat_is_identity * @m: An %IntegerMatrix * - * Returns true if @m is an identity matrix. + * Returns: true if @m is an identity matrix. * */ int intmat_is_identity(const IntegerMatrix *m) @@ -439,7 +441,7 @@ int intmat_is_identity(const IntegerMatrix *m) * intmat_is_inversion * @m: An %IntegerMatrix * - * Returns true if @m = -I, where I is an identity matrix. + * Returns: true if @m = -I, where I is an identity matrix. * */ int intmat_is_inversion(const IntegerMatrix *m) @@ -473,7 +475,7 @@ int intmat_is_inversion(const IntegerMatrix *m) * @a: An %IntegerMatrix * @b: An %IntegerMatrix * - * Returns true if @a = @b. + * Returns: true if @a = @b. * */ int intmat_equals(const IntegerMatrix *a, const IntegerMatrix *b) diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 38f8478b..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 */ @@ -554,9 +714,9 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); if ( mosflm->pid == -1 ) { - ERROR("Failed to fork for MOSFLM\n"); + 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; @@ -599,10 +763,22 @@ void run_mosflm(struct image *image, UnitCell *cell) sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); - if ( sval == -1 ) { - int err = errno; - ERROR("select() failed: %s\n", strerror(err)); - rval = 1; + if ( sval == -1 ) { + + const int err = errno; + + switch ( err ) { + + case EINTR: + STATUS("Restarting select()\n"); + break; + + default: + ERROR("select() failed: %s\n", strerror(err)); + rval = 1; + + } + } else if ( sval != 0 ) { rval = mosflm_readable(image, mosflm); } else { @@ -610,7 +786,6 @@ void run_mosflm(struct image *image, UnitCell *cell) rval = 1; } - } while ( !rval ); close(mosflm->pty); @@ -621,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 bf6a92f0..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; } @@ -236,6 +248,8 @@ static int integrate_peak(struct image *image, int cfs, int css, if ( p == NULL ) return 1; if ( p->no_index ) return 1; + if ( saturated != NULL ) *saturated = 0; + /* Determine regions where there is expected to be a peak */ p_cfs = cfs - p->min_fs; p_css = css - p->min_ss; /* Panel-relative coordinates */ @@ -293,10 +307,9 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx]; - /* Veto peak if it contains saturation in bg region */ + /* Check if peak contains saturation in bg region */ if ( use_max_adu && (val > p->max_adu) ) { if ( saturated != NULL ) *saturated = 1; - return 1; } bg_tot += val; @@ -353,10 +366,9 @@ static int integrate_peak(struct image *image, int cfs, int css, val = image->data[idx] - bg_mean; - /* Veto peak if it contains saturation */ + /* Check if peak contains saturation */ if ( use_max_adu && (val > p->max_adu) ) { if ( saturated != NULL ) *saturated = 1; - return 1; } pk_counts++; @@ -387,7 +399,8 @@ static int integrate_peak(struct image *image, int cfs, int css, static void search_peaks_in_panel(struct image *image, float threshold, float min_gradient, float min_snr, struct panel *p, - double ir_inn, double ir_mid, double ir_out) + double ir_inn, double ir_mid, double ir_out, + int use_saturated) { int fs, ss, stride; float *data; @@ -400,7 +413,7 @@ 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 nacc = 0; int ncull; @@ -419,13 +432,11 @@ static void search_peaks_in_panel(struct image *image, float threshold, double max; unsigned int did_something; int r; + int saturated; /* 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]; @@ -491,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, NULL); + 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 */ @@ -518,6 +528,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, continue; } + if ( r ) { + /* Bad region - don't detect peak */ + nrej_fail++; + continue; + } + /* Add using "better" coordinates */ image_add_feature(image->features, f_fs, f_ss, image, intensity, NULL); @@ -535,10 +551,13 @@ static void search_peaks_in_panel(struct image *image, float threshold, ncull = 0; } -// STATUS("%i accepted, %i box, %i proximity, %i outside panel, " -// "%i in bad regions, %i with SNR < %g, %i badrow culled.\n", -// nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, -// nrej_snr, min_snr, ncull); + image->num_peaks += nacc; + + //STATUS("%i accepted, %i box, %i proximity, %i outside panel, " + // "%i failed integration, %i with SNR < %g, %i badrow culled, " + // "%i saturated.\n", + // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_fail, + // nrej_snr, min_snr, ncull, nrej_sat); if ( ncull != 0 ) { STATUS("WARNING: %i peaks were badrow culled. This feature" @@ -549,7 +568,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, void search_peaks(struct image *image, float threshold, float min_gradient, - float min_snr, double ir_inn, double ir_mid, double ir_out) + float min_snr, double ir_inn, double ir_mid, double ir_out, + int use_saturated) { int i; @@ -557,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++ ) { @@ -564,35 +586,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient, if ( p->no_index ) continue; search_peaks_in_panel(image, threshold, min_gradient, - min_snr, p, ir_inn, ir_mid, ir_out); + min_snr, p, ir_inn, ir_mid, ir_out, + use_saturated); } } -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); @@ -602,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; } @@ -694,40 +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) +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; + + 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; @@ -794,10 +800,19 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, &intensity, &sigma, ir_inn, ir_mid, ir_out, 1, bgMasks[pnum], &saturated); - if ( saturated ) image->n_saturated++; + if ( saturated ) { + n_saturated++; + if ( !integrate_saturated ) r = 1; + } - /* I/sigma(I) cutoff */ - if ( intensity/sigma < min_snr ) r = 1; + /* I/sigma(I) cutoff + * Rejects reflections below --min-integration-snr, or if the + * SNR is clearly silly. Silly indicates that the intensity + * was zero. */ + snr = fabs(intensity)/sigma; + if ( isnan(snr) || (snr < min_snr) ) { + r = 1; + } /* Record intensity and set redundancy to 1 on success */ if ( r == 0 ) { @@ -808,7 +823,6 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, set_redundancy(refl, 0); } - snr = intensity / sigma; if ( snr > 1.0 ) { if ( first ) { av = snr; @@ -818,27 +832,61 @@ 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++ ) { - free(bgMasks[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; } - free(bgMasks); - image->diffracting_resolution = 0.0; + 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); + } - free(il); + for ( i=0; i<image->det->n_panels; i++ ) { + free(bgMasks[i]); + } + free(bgMasks); } void validate_peaks(struct image *image, double min_snr, - int ir_inn, int ir_mid, int ir_out) + int ir_inn, int ir_mid, int ir_out, int use_saturated) { int i, n; ImageFeatureList *flist; - int n_wtf, n_int, n_dft, n_snr, n_prx; + int n_wtf, n_int, n_dft, n_snr, n_prx, n_sat; flist = image_feature_list_new(); if ( flist == NULL ) return; @@ -846,7 +894,7 @@ void validate_peaks(struct image *image, double min_snr, n = image_feature_count(image->features); /* Loop over peaks, putting each one through the integrator */ - n_wtf = 0; n_int = 0; n_dft = 0; n_snr = 0; n_prx = 0; + n_wtf = 0; n_int = 0; n_dft = 0; n_snr = 0; n_prx = 0; n_sat = 0; for ( i=0; i<n; i++ ) { struct imagefeature *f; @@ -856,6 +904,7 @@ void validate_peaks(struct image *image, double min_snr, double f_fs, f_ss; double intensity, sigma; struct panel *p; + int saturated; f = image_get_feature(image->features, i); if ( f == NULL ) { @@ -871,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, NULL); + 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) ) @@ -903,9 +957,11 @@ void validate_peaks(struct image *image, double min_snr, } //STATUS("HDF5: %i peaks, validated: %i. WTF: %i, integration: %i," - // " drifted: %i, SNR: %i, proximity: %i\n", + // " drifted: %i, SNR: %i, proximity: %i, saturated: %i\n", // n, image_feature_count(flist), - // n_wtf, n_int, n_dft, n_snr, n_prx); + // 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 39fdf54c..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. * @@ -39,21 +39,22 @@ extern void search_peaks(struct image *image, float threshold, float min_gradient, float min_snr, - double ir_inn, double ir_mid, double ir_out); + double ir_inn, double ir_mid, double ir_out, + int use_saturated); extern void integrate_reflections(struct image *image, int use_closer, int bgsub, double min_snr, - double ir_inn, double ir_mid, double ir_out); + double ir_inn, double ir_mid, double ir_out, + 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); extern void validate_peaks(struct image *image, double min_snr, - int ir_inn, int ir_mid, int ir_out); + int ir_inn, int ir_mid, int ir_out, + int use_saturated); #endif /* PEAKS_H */ 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/reflist-utils.c b/libcrystfel/src/reflist-utils.c index f2292929..04096ae4 100644 --- a/libcrystfel/src/reflist-utils.c +++ b/libcrystfel/src/reflist-utils.c @@ -437,6 +437,9 @@ double max_intensity(RefList *list) /** * res_cutoff: * @list: A %RefList + * @cell: A %UnitCell with which to calculate 1/d values for @list + * @min: Minimum acceptable value of 1/d + * @max: Maximum acceptable value of 1/d * * Applies a resolution cutoff to @list, returning the new version and freeing * the old version. 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 db0073c9..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. * @@ -1689,4 +1690,3 @@ const char *symmetry_name(const SymOpList *ops) { return ops->name; } - diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 951d9d03..b75693db 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -100,6 +100,11 @@ extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r); extern size_t notrail(char *s); extern void chomp(char *s); +/** + * AssplodeFlag: + * @ASSPLODE_NONE: Nothing + * @ASSPLODE_DUPS: Don't merge deliminators + **/ typedef enum { ASSPLODE_NONE = 0, ASSPLODE_DUPS = 1<<0 @@ -141,7 +146,8 @@ static inline double modulus_squared(double x, double y, double z) { static inline double distance3d(double x1, double y1, double z1, double x2, double y2, double z2) { - return modulus(x1-x2, y1-y2, z1-z2); + double d = modulus(x1-x2, y1-y2, z1-z2); + return d; } /* Answer in radians */ |