From 38089071300b8e04ed42236dd08d9055094fb3b8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Nov 2011 12:05:55 +0100 Subject: Introduce "libcrystfel" --- libcrystfel/src/beam-parameters.c | 133 ++++ libcrystfel/src/beam-parameters.h | 40 ++ libcrystfel/src/cell.c | 1164 ++++++++++++++++++++++++++++++++++ libcrystfel/src/cell.h | 107 ++++ libcrystfel/src/detector.c | 1233 +++++++++++++++++++++++++++++++++++++ libcrystfel/src/detector.h | 131 ++++ libcrystfel/src/hdf5-file.c | 817 ++++++++++++++++++++++++ libcrystfel/src/hdf5-file.h | 55 ++ libcrystfel/src/image.c | 144 +++++ libcrystfel/src/image.h | 184 ++++++ libcrystfel/src/list_tmp.h | 106 ++++ libcrystfel/src/reflist.c | 1048 +++++++++++++++++++++++++++++++ libcrystfel/src/reflist.h | 112 ++++ libcrystfel/src/thread-pool.c | 283 +++++++++ libcrystfel/src/thread-pool.h | 71 +++ libcrystfel/src/utils.c | 488 +++++++++++++++ libcrystfel/src/utils.h | 236 +++++++ 17 files changed, 6352 insertions(+) create mode 100644 libcrystfel/src/beam-parameters.c create mode 100644 libcrystfel/src/beam-parameters.h create mode 100644 libcrystfel/src/cell.c create mode 100644 libcrystfel/src/cell.h create mode 100644 libcrystfel/src/detector.c create mode 100644 libcrystfel/src/detector.h create mode 100644 libcrystfel/src/hdf5-file.c create mode 100644 libcrystfel/src/hdf5-file.h create mode 100644 libcrystfel/src/image.c create mode 100644 libcrystfel/src/image.h create mode 100644 libcrystfel/src/list_tmp.h create mode 100644 libcrystfel/src/reflist.c create mode 100644 libcrystfel/src/reflist.h create mode 100644 libcrystfel/src/thread-pool.c create mode 100644 libcrystfel/src/thread-pool.h create mode 100644 libcrystfel/src/utils.c create mode 100644 libcrystfel/src/utils.h (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/beam-parameters.c b/libcrystfel/src/beam-parameters.c new file mode 100644 index 00000000..ac30ad04 --- /dev/null +++ b/libcrystfel/src/beam-parameters.c @@ -0,0 +1,133 @@ +/* + * beam-parameters.c + * + * Beam parameters + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "beam-parameters.h" +#include "utils.h" + + +struct beam_params *get_beam_parameters(const char *filename) +{ + FILE *fh; + struct beam_params *b; + char *rval; + int reject; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + b = calloc(1, sizeof(struct beam_params)); + if ( b == NULL ) return NULL; + + b->fluence = -1.0; + b->beam_radius = -1.0; + b->photon_energy = -1.0; + b->bandwidth = -1.0; + b->divergence = -1.0; + b->dqe = -1.0; + b->adu_per_photon = -1.0; + + do { + + int n1; + char line[1024]; + int i; + char **bits; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; ifluence = atof(bits[2]); + } 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]); + } else if ( strcmp(bits[0], "beam/bandwidth") == 0 ) { + b->bandwidth = atof(bits[2]); + } else if ( strcmp(bits[0], "beam/divergence") == 0 ) { + b->divergence = atof(bits[2]); + } else if ( strcmp(bits[0], "detector/dqe") == 0 ) { + b->dqe = atof(bits[2]); + } else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) { + b->adu_per_photon = atof(bits[2]); + } else { + ERROR("Unrecognised field '%s'\n", bits[0]); + } + + for ( i=0; ifluence < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/fluence'.\n"); + reject = 1; + } + if ( b->beam_radius < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/radius'.\n"); + reject = 1; + } + if ( b->photon_energy < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'beam/photon_energy'.\n"); + reject = 1; + } + if ( b->bandwidth < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/bandwidth'.\n"); + reject = 1; + } + if ( b->divergence < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/divergence'.\n"); + reject = 1; + } + if ( b->dqe < 0.0 ) { + ERROR("Invalid or unspecified value for 'detector/dqe'.\n"); + reject = 1; + } + if ( b->adu_per_photon < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'detector/adu_per_photon'.\n"); + reject = 1; + } + + if ( reject ) { + ERROR("Please fix the above problems with the beam" + " parameters file and try again.\n"); + return NULL; + } + + return b; +} diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h new file mode 100644 index 00000000..411ab449 --- /dev/null +++ b/libcrystfel/src/beam-parameters.h @@ -0,0 +1,40 @@ +/* + * beam-parameters.h + * + * Beam parameters + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef BEAM_PARAMETERS_H +#define BEAM_PARAMETERS_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +struct beam_params +{ + double fluence; /* photons per pulse */ + double beam_radius; /* metres */ + double photon_energy; /* eV per photon */ + double bandwidth; /* FWHM(wavelength) over wavelength. + * Note: current simulation code just uses + * a rectangular distribution with this as + * its (full) width. */ + double divergence; /* divergence (radians) */ + + double dqe; /* Detector DQE (fraction) */ + double adu_per_photon; /* Detector "gain" */ +}; + + +extern struct beam_params *get_beam_parameters(const char *filename); + + +#endif /* BEAM_PARAMETERS_H */ diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c new file mode 100644 index 00000000..c31697bc --- /dev/null +++ b/libcrystfel/src/cell.c @@ -0,0 +1,1164 @@ +/* + * cell.c + * + * Unit Cell Calculations + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include "cell.h" +#include "utils.h" +#include "image.h" + + +/** + * SECTION:unitcell + * @short_description: Unit cell + * @title: UnitCell + * @section_id: + * @see_also: + * @include: "cell.h" + * @Image: + * + * This structure represents a unit cell. + */ + + + +/* Weighting factor of lengths relative to angles */ +#define LWEIGHT (10.0e-9) + +typedef enum { + CELL_REP_CRYST, + CELL_REP_CART, + CELL_REP_RECIP +} CellRepresentation; + +struct _unitcell { + + CellRepresentation rep; + + /* Crystallographic representation */ + double a; /* m */ + double b; /* m */ + double c; /* m */ + double alpha; /* Radians */ + double beta; /* Radians */ + double gamma; /* Radians */ + + /* Cartesian representation */ + double ax; double bx; double cx; + double ay; double by; double cy; + double az; double bz; double cz; + + /* Cartesian representation of reciprocal axes */ + double axs; double bxs; double cxs; + double ays; double bys; double cys; + double azs; double bzs; double czs; + + char *pointgroup; + char *spacegroup; +}; + + +/************************** Setters and Constructors **************************/ + + +/** + * cell_new: + * + * Create a new %UnitCell. + * + * Returns: the new unit cell, or NULL on failure. + * + */ +UnitCell *cell_new() +{ + UnitCell *cell; + + cell = malloc(sizeof(UnitCell)); + if ( cell == NULL ) return NULL; + + cell->a = 1.0; + cell->b = 1.0; + cell->c = 1.0; + cell->alpha = M_PI_2; + cell->beta = M_PI_2; + cell->gamma = M_PI_2; + + cell->rep = CELL_REP_CRYST; + + cell->pointgroup = strdup("1"); + cell->spacegroup = strdup("P 1"); + + return cell; +} + + +/** + * cell_free: + * @cell: A %UnitCell to free. + * + * Frees a %UnitCell, and all internal resources concerning that cell. + * + */ +void cell_free(UnitCell *cell) +{ + if ( cell == NULL ) return; + free(cell->pointgroup); + free(cell->spacegroup); + free(cell); +} + + +void cell_set_parameters(UnitCell *cell, double a, double b, double c, + double alpha, double beta, double gamma) +{ + if ( cell == NULL ) return; + + cell->a = a; + cell->b = b; + cell->c = c; + cell->alpha = alpha; + cell->beta = beta; + cell->gamma = gamma; + + cell->rep = CELL_REP_CRYST; +} + + +void cell_set_cartesian(UnitCell *cell, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz) +{ + if ( cell == NULL ) return; + + cell->ax = ax; cell->ay = ay; cell->az = az; + cell->bx = bx; cell->by = by; cell->bz = bz; + cell->cx = cx; cell->cy = cy; cell->cz = cz; + + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az) +{ + if ( cell == NULL ) return; + cell->ax = ax; cell->ay = ay; cell->az = az; + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz) +{ + if ( cell == NULL ) return; + cell->bx = bx; cell->by = by; cell->bz = bz; + cell->rep = CELL_REP_CART; +} + + +void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz) +{ + if ( cell == NULL ) return; + cell->cx = cx; cell->cy = cy; cell->cz = cz; + cell->rep = CELL_REP_CART; +} + + +UnitCell *cell_new_from_parameters(double a, double b, double c, + double alpha, double beta, double gamma) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell_set_parameters(cell, a, b, c, alpha, beta, gamma); + + return cell; +} + + +UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, + struct rvec cs) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell->axs = as.u; cell->ays = as.v; cell->azs = as.w; + cell->bxs = bs.u; cell->bys = bs.v; cell->bzs = bs.w; + cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w; + + cell->rep = CELL_REP_RECIP; + + return cell; +} + + +UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c) +{ + UnitCell *cell; + + cell = cell_new(); + if ( cell == NULL ) return NULL; + + cell->ax = a.u; cell->ay = a.v; cell->az = a.w; + cell->bx = b.u; cell->by = b.v; cell->bz = b.w; + cell->cx = c.u; cell->cy = c.v; cell->cz = c.w; + + cell->rep = CELL_REP_CART; + + return cell; +} + + +UnitCell *cell_new_from_cell(UnitCell *orig) +{ + UnitCell *new; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + new = cell_new(); + + cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz); + cell_set_pointgroup(new, orig->pointgroup); + cell_set_spacegroup(new, orig->spacegroup); + + return new; +} + + +void cell_set_reciprocal(UnitCell *cell, + double asx, double asy, double asz, + double bsx, double bsy, double bsz, + double csx, double csy, double csz) +{ + if ( cell == NULL ) return; + + cell->axs = asx; cell->ays = asy; cell->azs = asz; + cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz; + cell->cxs = csx; cell->cys = csy; cell->czs = csz; + + cell->rep = CELL_REP_RECIP; +} + + +void cell_set_spacegroup(UnitCell *cell, const char *sym) +{ + free(cell->spacegroup); + cell->spacegroup = strdup(sym); +} + + +void cell_set_pointgroup(UnitCell *cell, const char *sym) +{ + free(cell->pointgroup); + cell->pointgroup = strdup(sym); +} + + +/************************* Getter helper functions ****************************/ + +static int cell_crystallographic_to_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz) +{ + double tmp, V, cosalphastar, cstar; + + /* Firstly: Get a in terms of x, y and z + * +a (cryst) is defined to lie along +x (cart) */ + *ax = cell->a; + *ay = 0.0; + *az = 0.0; + + /* b in terms of x, y and z + * b (cryst) is defined to lie in the xy (cart) plane */ + *bx = cell->b*cos(cell->gamma); + *by = cell->b*sin(cell->gamma); + *bz = 0.0; + + tmp = cos(cell->alpha)*cos(cell->alpha) + + cos(cell->beta)*cos(cell->beta) + + cos(cell->gamma)*cos(cell->gamma) + - 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma); + V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp); + + cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha); + cosalphastar /= sin(cell->beta)*sin(cell->gamma); + + cstar = (cell->a * cell->b * sin(cell->gamma))/V; + + /* c in terms of x, y and z */ + *cx = cell->c*cos(cell->beta); + *cy = -cell->c*sin(cell->beta)*cosalphastar; + *cz = 1.0/cstar; + + return 0; +} + + +/* Why yes, I do enjoy long argument lists...! */ +static int cell_invert(double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz) +{ + int s; + gsl_matrix *m; + gsl_matrix *inv; + gsl_permutation *perm; + + m = gsl_matrix_alloc(3, 3); + if ( m == NULL ) { + ERROR("Couldn't allocate memory for matrix\n"); + return 1; + } + gsl_matrix_set(m, 0, 0, ax); + gsl_matrix_set(m, 0, 1, bx); + gsl_matrix_set(m, 0, 2, cx); + gsl_matrix_set(m, 1, 0, ay); + gsl_matrix_set(m, 1, 1, by); + gsl_matrix_set(m, 1, 2, cy); + gsl_matrix_set(m, 2, 0, az); + gsl_matrix_set(m, 2, 1, bz); + gsl_matrix_set(m, 2, 2, cz); + + /* Invert */ + perm = gsl_permutation_alloc(m->size1); + if ( perm == NULL ) { + ERROR("Couldn't allocate permutation\n"); + gsl_matrix_free(m); + return 1; + } + inv = gsl_matrix_alloc(m->size1, m->size2); + if ( inv == NULL ) { + ERROR("Couldn't allocate inverse\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + if ( gsl_linalg_LU_decomp(m, perm, &s) ) { + ERROR("Couldn't decompose matrix\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + if ( gsl_linalg_LU_invert(m, perm, inv) ) { + ERROR("Couldn't invert matrix\n"); + gsl_matrix_free(m); + gsl_permutation_free(perm); + return 1; + } + gsl_permutation_free(perm); + gsl_matrix_free(m); + + /* Transpose */ + gsl_matrix_transpose(inv); + + *asx = gsl_matrix_get(inv, 0, 0); + *bsx = gsl_matrix_get(inv, 0, 1); + *csx = gsl_matrix_get(inv, 0, 2); + *asy = gsl_matrix_get(inv, 1, 0); + *bsy = gsl_matrix_get(inv, 1, 1); + *csy = gsl_matrix_get(inv, 1, 2); + *asz = gsl_matrix_get(inv, 2, 0); + *bsz = gsl_matrix_get(inv, 2, 1); + *csz = gsl_matrix_get(inv, 2, 2); + + gsl_matrix_free(inv); + + return 0; +} + + +/********************************** Getters ***********************************/ + +int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, + double *alpha, double *beta, double *gamma) +{ + double ax, ay, az, bx, by, bz, cx, cy, cz; + + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Direct response */ + *a = cell->a; + *b = cell->b; + *c = cell->c; + *alpha = cell->alpha; + *beta = cell->beta; + *gamma = cell->gamma; + return 0; + + case CELL_REP_CART: + /* Convert cartesian -> crystallographic */ + *a = modulus(cell->ax, cell->ay, cell->az); + *b = modulus(cell->bx, cell->by, cell->bz); + *c = modulus(cell->cx, cell->cy, cell->cz); + + *alpha = angle_between(cell->bx, cell->by, cell->bz, + cell->cx, cell->cy, cell->cz); + *beta = angle_between(cell->ax, cell->ay, cell->az, + cell->cx, cell->cy, cell->cz); + *gamma = angle_between(cell->ax, cell->ay, cell->az, + cell->bx, cell->by, cell->bz); + return 0; + + case CELL_REP_RECIP: + /* Convert reciprocal -> crystallographic. + * Start by converting reciprocal -> cartesian */ + cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + /* Now convert cartesian -> crystallographic */ + *a = modulus(ax, ay, az); + *b = modulus(bx, by, bz); + *c = modulus(cx, cy, cz); + + *alpha = angle_between(bx, by, bz, cx, cy, cz); + *beta = angle_between(ax, ay, az, cx, cy, cz); + *gamma = angle_between(ax, ay, az, bx, by, bz); + return 0; + } + + return 1; +} + + +int cell_get_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz) +{ + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Convert crystallographic -> cartesian. */ + return cell_crystallographic_to_cartesian(cell, + ax, ay, az, + bx, by, bz, + cx, cy, cz); + + case CELL_REP_CART: + /* Direct response */ + *ax = cell->ax; *ay = cell->ay; *az = cell->az; + *bx = cell->bx; *by = cell->by; *bz = cell->bz; + *cx = cell->cx; *cy = cell->cy; *cz = cell->cz; + return 0; + + case CELL_REP_RECIP: + /* Convert reciprocal -> cartesian */ + return cell_invert(cell->axs, cell->ays, cell->azs, + cell->bxs, cell->bys, cell->bzs, + cell->cxs, cell->cys, cell->czs, + ax, ay, az, bx, by, bz, cx, cy, cz); + + } + + return 1; +} + + +int cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz) +{ + int r; + double ax, ay, az, bx, by, bz, cx, cy, cz; + if ( cell == NULL ) return 1; + + switch ( cell->rep ) { + + case CELL_REP_CRYST: + /* Convert crystallographic -> reciprocal */ + r = cell_crystallographic_to_cartesian(cell, + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + if ( r ) return r; + return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz, + asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + case CELL_REP_CART: + /* Convert cartesian -> reciprocal */ + cell_invert(cell->ax, cell->ay, cell->az, + cell->bx, cell->by, cell->bz, + cell->cx, cell->cy, cell->cz, + asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + return 0; + + case CELL_REP_RECIP: + /* Direct response */ + *asx = cell->axs; *asy = cell->ays; *asz = cell->azs; + *bsx = cell->bxs; *bsy = cell->bys; *bsz = cell->bzs; + *csx = cell->cxs; *csy = cell->cys; *csz = cell->czs; + return 0; + + } + + return 1; +} + + +const char *cell_get_pointgroup(UnitCell *cell) +{ + return cell->pointgroup; +} + + +const char *cell_get_spacegroup(UnitCell *cell) +{ + return cell->spacegroup; +} + + + + + +/********************************* Utilities **********************************/ + +static const char *cell_rep(UnitCell *cell) +{ + switch ( cell->rep ) { + case CELL_REP_CRYST: + return "crystallographic, direct space"; + case CELL_REP_CART: + return "cartesian, direct space"; + case CELL_REP_RECIP: + return "cartesian, reciprocal space"; + } + + return "unknown"; +} + + +/** + * cell_rotate: + * @in: A %UnitCell to rotate + * @quat: A %quaternion + * + * Rotate a %UnitCell using a %quaternion. + * + * Returns: a newly allocated rotated copy of @in. + * + */ +UnitCell *cell_rotate(UnitCell *in, struct quaternion quat) +{ + struct rvec a, b, c; + struct rvec an, bn, cn; + UnitCell *out = cell_new_from_cell(in); + + cell_get_cartesian(in, &a.u, &a.v, &a.w, + &b.u, &b.v, &b.w, + &c.u, &c.v, &c.w); + + an = quat_rot(a, quat); + bn = quat_rot(b, quat); + cn = quat_rot(c, quat); + + cell_set_cartesian(out, an.u, an.v, an.w, + bn.u, bn.v, bn.w, + cn.u, cn.v, cn.w); + + return out; +} + + +void cell_print(UnitCell *cell) +{ + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double a, b, c, alpha, beta, gamma; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + STATUS(" a b c alpha beta gamma\n"); + STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n", + a*1e9, b*1e9, c*1e9, + rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az); + STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz); + STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz); + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + asx, asy, asz, modulus(asx, asy, asz)); + STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + bsx, bsy, bsz, modulus(bsx, bsy, bsz)); + STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + csx, csy, csz, modulus(csx, csy, csz)); + + STATUS("Point group: %s\n", cell_get_pointgroup(cell)); + STATUS("Cell representation is %s.\n", cell_rep(cell)); +} + + +#define MAX_CAND (1024) + +static int right_handed(struct rvec a, struct rvec b, struct rvec c) +{ + struct rvec aCb; + double aCb_dot_c; + + /* "a" cross "b" */ + aCb.u = a.v*b.w - a.w*b.v; + aCb.v = - (a.u*b.w - a.w*b.u); + aCb.w = a.u*b.v - a.v*b.u; + + /* "a cross b" dot "c" */ + aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; + + if ( aCb_dot_c > 0.0 ) return 1; + return 0; +} + + +struct cvec { + struct rvec vec; + float na; + float nb; + float nc; + float fom; +}; + + +static int same_vector(struct cvec a, struct cvec b) +{ + if ( a.na != b.na ) return 0; + if ( a.nb != b.nb ) return 0; + if ( a.nc != b.nc ) return 0; + return 1; +} + + +/* Attempt to make 'cell' fit into 'template' somehow */ +UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, + int reduce) +{ + signed int n1l, n2l, n3l; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + int i, j; + double lengths[3]; + double angles[3]; + struct cvec *cand[3]; + UnitCell *new_cell = NULL; + float best_fom = +999999999.9; /* Large number.. */ + int ncand[3] = {0,0,0}; + signed int ilow, ihigh; + float ltl = 5.0; /* percent */ + float angtol = deg2rad(1.5); + + if ( cell_get_reciprocal(template, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell for template.\n"); + return NULL; + } + + lengths[0] = modulus(asx, asy, asz); + lengths[1] = modulus(bsx, bsy, bsz); + lengths[2] = modulus(csx, csy, csz); + + angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); + angles[1] = angle_between(asx, asy, asz, csx, csy, csz); + angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); + + cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); + cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); + + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal cell.\n"); + return NULL; + } + + if ( reduce ) { + ilow = -2; ihigh = 4; + } else { + ilow = 0; ihigh = 1; + } + + /* Negative values mean 1/n, positive means n, zero means zero */ + for ( n1l=ilow; n1l<=ihigh; n1l++ ) { + for ( n2l=ilow; n2l<=ihigh; n2l++ ) { + for ( n3l=ilow; n3l<=ihigh; n3l++ ) { + + float n1, n2, n3; + signed int b1, b2, b3; + + n1 = (n1l>=0) ? (n1l) : (1.0/n1l); + n2 = (n2l>=0) ? (n2l) : (1.0/n2l); + n3 = (n3l>=0) ? (n3l) : (1.0/n3l); + + if ( !reduce ) { + if ( n1l + n2l + n3l > 1 ) continue; + } + + /* 'bit' values can be +1 or -1 */ + for ( b1=-1; b1<=1; b1+=2 ) { + for ( b2=-1; b2<=1; b2+=2 ) { + for ( b3=-1; b3<=1; b3+=2 ) { + + double tx, ty, tz; + double tlen; + int i; + + n1 *= b1; n2 *= b2; n3 *= b3; + + tx = n1*asx + n2*bsx + n3*csx; + ty = n1*asy + n2*bsy + n3*csy; + tz = n1*asz + n2*bsz + n3*csz; + tlen = modulus(tx, ty, tz); + + /* Test modulus for agreement with moduli of template */ + for ( i=0; i<3; i++ ) { + + if ( !within_tolerance(lengths[i], tlen, ltl) ) + continue; + + cand[i][ncand[i]].vec.u = tx; + cand[i][ncand[i]].vec.v = ty; + cand[i][ncand[i]].vec.w = tz; + cand[i][ncand[i]].na = n1; + cand[i][ncand[i]].nb = n2; + cand[i][ncand[i]].nc = n3; + cand[i][ncand[i]].fom = fabs(lengths[i] - tlen); + if ( ncand[i] == MAX_CAND ) { + ERROR("Too many candidates\n"); + } else { + ncand[i]++; + } + } + + } + } + } + + } + } + } + + if ( verbose ) { + STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); + } + + for ( i=0; i angtol ) continue; + + fom1 = fabs(ang - angles[2]); + + for ( k=0; k angtol ) continue; + + fom2 = fom1 + fabs(ang - angles[1]); + + /* Finally, the angle between the current candidate for + * axis 1 and the kth candidate for axis 2 */ + ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, + cand[1][j].vec.w, cand[2][k].vec.u, + cand[2][k].vec.v, cand[2][k].vec.w); + + /* ... it should be angle 0 ... */ + if ( fabs(ang - angles[0]) > angtol ) continue; + + /* Unit cell must be right-handed */ + if ( !right_handed(cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec) ) continue; + + fom3 = fom2 + fabs(ang - angles[0]); + fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom + + cand[2][k].fom); + + if ( fom3 < best_fom ) { + if ( new_cell != NULL ) free(new_cell); + new_cell = cell_new_from_reciprocal_axes( + cand[0][i].vec, cand[1][j].vec, + cand[2][k].vec); + best_fom = fom3; + } + + } + + } + } + + free(cand[0]); + free(cand[1]); + free(cand[2]); + + return new_cell; +} + + + +UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + int i; + double lengths[3]; + int used[3]; + struct rvec real_a, real_b, real_c; + struct rvec params[3]; + double alen, blen, clen; + float ltl = 5.0; /* percent */ + int have_real_a; + int have_real_b; + int have_real_c; + + /* Get the lengths to match */ + if ( cell_get_cartesian(template, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell for template.\n"); + return NULL; + } + alen = modulus(ax, ay, az); + blen = modulus(bx, by, bz); + clen = modulus(cx, cy, cz); + + /* Get the lengths from the cell and turn them into anonymous vectors */ + if ( cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz) ) + { + ERROR("Couldn't get cell.\n"); + return NULL; + } + lengths[0] = modulus(ax, ay, az); + lengths[1] = modulus(bx, by, bz); + lengths[2] = modulus(cx, cy, cz); + used[0] = 0; used[1] = 0; used[2] = 0; + params[0].u = ax; params[0].v = ay; params[0].w = az; + params[1].u = bx; params[1].v = by; params[1].w = bz; + params[2].u = cx; params[2].v = cy; params[2].w = cz; + + real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; + real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; + real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; + + /* Check each vector against a and b */ + have_real_a = 0; + have_real_b = 0; + for ( i=0; i<3; i++ ) { + if ( within_tolerance(lengths[i], alen, ltl) + && !used[i] && !have_real_a ) + { + used[i] = 1; + memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); + have_real_a = 1; + } + if ( within_tolerance(lengths[i], blen, ltl) + && !used[i] && !have_real_b ) + { + used[i] = 1; + memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); + have_real_b = 1; + } + } + + /* Have we matched both a and b? */ + if ( !(have_real_a && have_real_b) ) return NULL; + + /* "c" is "the other one" */ + have_real_c = 0; + for ( i=0; i<3; i++ ) { + if ( !used[i] ) { + memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); + have_real_c = 1; + } + } + + if ( !have_real_c ) { + ERROR("Huh? Couldn't find the third vector.\n"); + ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); + return NULL; + } + + /* Flip c if not right-handed */ + if ( !right_handed(real_a, real_b, real_c) ) { + real_c.u = -real_c.u; + real_c.v = -real_c.v; + real_c.w = -real_c.w; + } + + return cell_new_from_direct_axes(real_a, real_b, real_c); +} + + +/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ +double resolution(UnitCell *cell, signed int h, signed int k, signed int l) +{ + double a, b, c, alpha, beta, gamma; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha) + - cos(beta)*cos(beta) + - cos(gamma)*cos(gamma) + + 2*cos(alpha)*cos(beta)*cos(gamma) ); + + const double S11 = b*b*c*c*sin(alpha)*sin(alpha); + const double S22 = a*a*c*c*sin(beta)*sin(beta); + const double S33 = a*a*b*b*sin(gamma)*sin(gamma); + const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma)); + const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha)); + const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta)); + + const double brackets = S11*h*h + S22*k*k + S33*l*l + + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l; + const double oneoverdsq = brackets / Vsq; + const double oneoverd = sqrt(oneoverdsq); + + return oneoverd / 2; +} + + +static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym) +{ + char *new = NULL; + + if ( strcmp(sym, "P 1") == 0 ) new = "1"; + if ( strcmp(sym, "P 63") == 0 ) new = "6"; + if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222"; + if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222"; + if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422"; + + if ( new != NULL ) { + cell_set_pointgroup(cell, new); + } else { + ERROR("Can't determine point group for '%s'\n", sym); + } +} + + +UnitCell *load_cell_from_pdb(const char *filename) +{ + FILE *fh; + char *rval; + UnitCell *cell = NULL; + + fh = fopen(filename, "r"); + if ( fh == NULL ) { + ERROR("Couldn't open '%s'\n", filename); + return NULL; + } + + do { + + char line[1024]; + + rval = fgets(line, 1023, fh); + + if ( strncmp(line, "CRYST1", 6) == 0 ) { + + float a, b, c, al, be, ga; + char as[10], bs[10], cs[10]; + char als[8], bes[8], gas[8]; + char *sym; + int r; + + memcpy(as, line+6, 9); as[9] = '\0'; + memcpy(bs, line+15, 9); bs[9] = '\0'; + memcpy(cs, line+24, 9); cs[9] = '\0'; + memcpy(als, line+33, 7); als[7] = '\0'; + memcpy(bes, line+40, 7); bes[7] = '\0'; + memcpy(gas, line+47, 7); gas[7] = '\0'; + + r = sscanf(as, "%f", &a); + r += sscanf(bs, "%f", &b); + r += sscanf(cs, "%f", &c); + r += sscanf(als, "%f", &al); + r += sscanf(bes, "%f", &be); + r += sscanf(gas, "%f", &ga); + + if ( r != 6 ) { + STATUS("Couldn't understand CRYST1 line.\n"); + continue; + } + + cell = cell_new_from_parameters(a*1e-10, + b*1e-10, c*1e-10, + deg2rad(al), + deg2rad(be), + deg2rad(ga)); + + if ( strlen(line) > 65 ) { + sym = strndup(line+55, 10); + notrail(sym); + cell_set_pointgroup_from_pdb(cell, sym); + cell_set_spacegroup(cell, sym); + free(sym); + } else { + cell_set_pointgroup_from_pdb(cell, "P 1"); + cell_set_spacegroup(cell, "P 1"); + ERROR("CRYST1 line without space group.\n"); + } + + break; /* Done */ + } + + } while ( rval != NULL ); + + fclose(fh); + + return cell; +} + + +#ifdef GSL_FUDGE +/* Force the linker to bring in CBLAS to make GSL happy */ +void cell_fudge_gslcblas() +{ + STATUS("%p\n", cblas_sgemm); +} +#endif + + +UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot) +{ + UnitCell *out; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xnew, ynew, znew; + + cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */ + xnew = asx*cos(omega) + asy*sin(omega); + ynew = -asx*sin(omega) + asy*cos(omega); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(omega) + bsy*sin(omega); + ynew = -bsx*sin(omega) + bsy*cos(omega); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(omega) + csy*sin(omega); + ynew = -csx*sin(omega) + csy*cos(omega); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "phi" about +x (not parallel to anything specific) */ + xnew = asx; + ynew = asy*cos(phi) + asz*sin(phi); + znew = -asy*sin(phi) + asz*cos(phi); + asx = xnew; asy = ynew; asz = znew; + xnew = bsx; + ynew = bsy*cos(phi) + bsz*sin(phi); + znew = -bsy*sin(phi) + bsz*cos(phi); + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx; + ynew = csy*cos(phi) + csz*sin(phi); + znew = -csy*sin(phi) + csz*cos(phi); + csx = xnew; csy = ynew; csz = znew; + + /* Rotate by "rot" about the new +z (in-plane rotation) */ + xnew = asx*cos(rot) + asy*sin(rot); + ynew = -asx*sin(rot) + asy*cos(rot); + znew = asz; + asx = xnew; asy = ynew; asz = znew; + xnew = bsx*cos(rot) + bsy*sin(rot); + ynew = -bsx*sin(rot) + bsy*cos(rot); + znew = bsz; + bsx = xnew; bsy = ynew; bsz = znew; + xnew = csx*cos(rot) + csy*sin(rot); + ynew = -csx*sin(rot) + csy*cos(rot); + znew = csz; + csx = xnew; csy = ynew; csz = znew; + + out = cell_new_from_cell(in); + cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + return out; +} + + +int cell_is_sensible(UnitCell *cell) +{ + double a, b, c, al, be, ga; + + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + if ( al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be - ga >= 2.0*M_PI ) return 0; + if ( al - be + ga >= 2.0*M_PI ) return 0; + if ( - al + be + ga >= 2.0*M_PI ) return 0; + if ( al + be + ga <= 0.0 ) return 0; + if ( al + be - ga <= 0.0 ) return 0; + if ( al - be + ga <= 0.0 ) return 0; + if ( - al + be + ga <= 0.0 ) return 0; + if ( isnan(al) ) return 0; + if ( isnan(be) ) return 0; + if ( isnan(ga) ) return 0; + return 1; +} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h new file mode 100644 index 00000000..b5d31fc6 --- /dev/null +++ b/libcrystfel/src/cell.h @@ -0,0 +1,107 @@ +/* + * cell.h + * + * Unit Cell Calculations + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef CELL_H +#define CELL_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "utils.h" + +/* A 3D vector in reciprocal space */ +struct rvec +{ + double u; + double v; + double w; +}; + + +/** + * UnitCell: + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + **/ +typedef struct _unitcell UnitCell; + +extern UnitCell *cell_new(void); +extern UnitCell *cell_new_from_cell(UnitCell *orig); +extern void cell_free(UnitCell *cell); + +/* Lengths in m, angles in radians */ +extern UnitCell *cell_new_from_parameters(double a, double b, double c, + double alpha, double beta, double gamma); + +extern UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs, + struct rvec cs); + +extern UnitCell *cell_new_from_direct_axes(struct rvec as, struct rvec bs, + struct rvec cs); + +extern void cell_set_cartesian(UnitCell *cell, + double ax, double ay, double az, + double bx, double by, double bz, + double cx, double cy, double cz); + +extern void cell_set_parameters(UnitCell *cell, double a, double b, double c, + double alpha, double beta, double gamma); + +extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az); +extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz); +extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz); +extern void cell_set_spacegroup(UnitCell *cell, const char *sym); +extern void cell_set_pointgroup(UnitCell *cell, const char *sym); + + +extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c, + double *alpha, double *beta, double *gamma); + +extern int cell_get_cartesian(UnitCell *cell, + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz); + +extern int cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz); + +extern void cell_set_reciprocal(UnitCell *cell, + double asx, double asy, double asz, + double bsx, double bsy, double bsz, + double csx, double csy, double csz); + +extern const char *cell_get_pointgroup(UnitCell *cell); + +extern const char *cell_get_spacegroup(UnitCell *cell); + +extern double resolution(UnitCell *cell, + signed int h, signed int k, signed int l); + +extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat); +extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, + double rot); + +extern void cell_print(UnitCell *cell); + +extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, + int reduce); + +extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template); + +extern UnitCell *load_cell_from_pdb(const char *filename); + +extern int cell_is_sensible(UnitCell *cell); + +#endif /* CELL_H */ diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c new file mode 100644 index 00000000..54bae1d7 --- /dev/null +++ b/libcrystfel/src/detector.c @@ -0,0 +1,1233 @@ +/* + * detector.c + * + * Detector properties + * + * (c) 2006-2011 Thomas White + * (c) 2011 Rick Kirian + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include +#include +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" +#include "detector.h" +#include "beam-parameters.h" +#include "hdf5-file.h" + + +/** + * SECTION:detector + * @short_description: Detector geometry + * @title: Detector + * @section_id: + * @see_also: + * @include: "detector.h" + * @Image: + * + * This structure represents the detector geometry + */ + + +static int atob(const char *a) +{ + if ( strcasecmp(a, "true") == 0 ) return 1; + if ( strcasecmp(a, "false") == 0 ) return 0; + return atoi(a); +} + + +static int assplode_algebraic(const char *a_orig, char ***pbits) +{ + int len, i; + int nexp; + char **bits; + char *a; + int idx, istr; + + len = strlen(a_orig); + + /* Add plus at start if no sign there already */ + if ( (a_orig[0] != '+') && (a_orig[0] != '-') ) { + len += 1; + a = malloc(len+1); + snprintf(a, len+1, "+%s", a_orig); + a[len] = '\0'; + + } else { + a = strdup(a_orig); + } + + /* Count the expressions */ + nexp = 0; + for ( i=0; i= 0 ) bits[idx][istr] = '\0'; + idx++; + bits[idx] = malloc(len+1); + istr = 0; + } + + if ( !isdigit(ch) && (ch != '.') && (ch != 'x') && (ch != 'y') + && (ch != '+') && (ch != '-') ) + { + ERROR("Invalid character '%C' found.\n", ch); + return 0; + } + + assert(idx >= 0); + bits[idx][istr++] = ch; + + } + if ( idx >= 0 ) bits[idx][istr] = '\0'; + + *pbits = bits; + free(a); + + return nexp; +} + + +/* Parses the scan directions (accounting for possible rotation) + * Assumes all white spaces have been already removed */ +static int dir_conv(const char *a, double *sx, double *sy) +{ + int n; + char **bits; + int i; + + *sx = 0.0; *sy = 0.0; + + n = assplode_algebraic(a, &bits); + + if ( n == 0 ) { + ERROR("Invalid direction '%s'\n", a); + return 1; + } + + for ( i=0; i %5.2fx + %5.2fy\n", a, *sx, *sy); + + return 0; +} + + +struct rvec get_q_for_panel(struct panel *p, double fs, double ss, + double *ttp, double k) +{ + struct rvec q; + double twotheta, r, az; + double rx, ry; + double xs, ys; + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + 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; + + /* Calculate q-vector for this sub-pixel */ + r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + twotheta = atan2(r, p->clen); + az = atan2(ry, rx); + if ( ttp != NULL ) *ttp = twotheta; + + q.u = k * sin(twotheta)*cos(az); + q.v = k * sin(twotheta)*sin(az); + q.w = k * (cos(twotheta) - 1.0); + + return q; +} + + +struct rvec get_q(struct image *image, double fs, double ss, + double *ttp, double k) +{ + struct panel *p; + const unsigned int fsi = fs; + const unsigned int ssi = ss; /* Explicit rounding */ + + /* Determine which panel to use */ + p = find_panel(image->det, fsi, ssi); + assert(p != NULL); + + return get_q_for_panel(p, fs-(double)p->min_fs, ss-(double)p->min_ss, + ttp, k); +} + + +int in_bad_region(struct detector *det, double fs, double ss) +{ + double rx, ry; + struct panel *p; + double xs, ys; + int i; + + /* Determine which panel to use */ + p = find_panel(det, fs, ss); + + /* No panel found -> definitely bad! */ + if ( p == NULL ) return 1; + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + xs = (fs-(double)p->min_fs)*p->fsx + (ss-(double)p->min_ss)*p->ssx; + ys = (fs-(double)p->min_fs)*p->fsy + (ss-(double)p->min_ss)*p->ssy; + + rx = xs + p->cnx; + ry = ys + p->cny; + + for ( i=0; in_bad; i++ ) { + struct badregion *b = &det->bad[i]; + if ( rx < b->min_x ) continue; + if ( rx > b->max_x ) continue; + if ( ry < b->min_y ) continue; + if ( ry > b->max_y ) continue; + return 1; + } + + return 0; +} + + +double get_tt(struct image *image, double fs, double ss) +{ + double r, rx, ry; + struct panel *p; + double xs, ys; + + p = find_panel(image->det, fs, ss); + + /* Convert xs and ys, which are in fast scan/slow scan coordinates, + * to x and y */ + xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx; + ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy; + + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + + r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + return atan2(r, p->clen); +} + + +void record_image(struct image *image, int do_poisson) +{ + int x, y; + double total_energy, energy_density; + double ph_per_e; + double area; + double max_tt = 0.0; + int n_inf1 = 0; + int n_neg1 = 0; + int n_nan1 = 0; + int n_inf2 = 0; + int n_neg2 = 0; + int n_nan2 = 0; + + /* How many photons are scattered per electron? */ + area = M_PI*pow(image->beam->beam_radius, 2.0); + total_energy = image->beam->fluence * ph_lambda_to_en(image->lambda); + energy_density = total_energy / area; + ph_per_e = (image->beam->fluence /area) * pow(THOMSON_LENGTH, 2.0); + STATUS("Fluence = %8.2e photons, " + "Energy density = %5.3f kJ/cm^2, " + "Total energy = %5.3f microJ\n", + image->beam->fluence, energy_density/1e7, total_energy*1e6); + + for ( x=0; xwidth; x++ ) { + for ( y=0; yheight; y++ ) { + + double counts; + double cf; + double intensity, sa; + double pix_area, Lsq; + double xs, ys, rx, ry; + double dsq, proj_area; + struct panel *p; + + intensity = (double)image->data[x + image->width*y]; + if ( isinf(intensity) ) n_inf1++; + if ( intensity < 0.0 ) n_neg1++; + if ( isnan(intensity) ) n_nan1++; + + p = find_panel(image->det, x, y); + + /* Area of one pixel */ + pix_area = pow(1.0/p->res, 2.0); + Lsq = pow(p->clen, 2.0); + + /* Area of pixel as seen from crystal (approximate) */ + proj_area = pix_area * cos(image->twotheta[x + image->width*y]); + + /* Calculate distance from crystal to pixel */ + xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx; + ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy; + rx = (xs + p->cnx) / p->res; + ry = (ys + p->cny) / p->res; + dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + + /* Projected area of pixel divided by distance squared */ + sa = proj_area / (dsq + Lsq); + + if ( do_poisson ) { + counts = poisson_noise(intensity * ph_per_e + * sa * image->beam->dqe ); + } else { + cf = intensity * ph_per_e * sa * image->beam->dqe; + counts = cf; + } + + image->data[x + image->width*y] = counts + * image->beam->adu_per_photon; + + /* Sanity checks */ + if ( isinf(image->data[x+image->width*y]) ) n_inf2++; + if ( isnan(image->data[x+image->width*y]) ) n_nan2++; + if ( image->data[x+image->width*y] < 0.0 ) n_neg2++; + + if ( image->twotheta[x + image->width*y] > max_tt ) { + max_tt = image->twotheta[x + image->width*y]; + } + + } + progress_bar(x, image->width-1, "Post-processing"); + } + + STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", + rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); + + double tt_side = image->twotheta[(image->width/2)+image->width*0]; + STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n", + rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); + + tt_side = image->twotheta[0+image->width*(image->height/2)]; + STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n", + rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9); + + STATUS("Halve the d values to get the voxel size for a synthesis.\n"); + + if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) { + ERROR("WARNING: The raw calculation produced %i negative" + " values, %i infinities and %i NaNs.\n", + n_neg1, n_inf1, n_nan1); + ERROR("WARNING: After processing, there were %i negative" + " values, %i infinities and %i NaNs.\n", + n_neg2, n_inf2, n_nan2); + } +} + + +struct panel *find_panel(struct detector *det, double fs, double ss) +{ + int p; + + for ( p=0; pn_panels; p++ ) { + if ( (fs >= det->panels[p].min_fs) + && (fs <= det->panels[p].max_fs) + && (ss >= det->panels[p].min_ss) + && (ss <= det->panels[p].max_ss) ) { + return &det->panels[p]; + } + } + + return NULL; +} + + +int find_panel_number(struct detector *det, int fs, int ss) +{ + int p; + + for ( p=0; pn_panels; p++ ) { + if ( (fs >= det->panels[p].min_fs) + && (fs <= det->panels[p].max_fs) + && (ss >= det->panels[p].min_ss) + && (ss <= det->panels[p].max_ss) ) return p; + } + + return -1; +} + + +void fill_in_values(struct detector *det, struct hdfile *f) +{ + int i; + + for ( i=0; in_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + 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; + + } +} + + +static struct panel *new_panel(struct detector *det, const char *name) +{ + struct panel *new; + + det->n_panels++; + det->panels = realloc(det->panels, det->n_panels*sizeof(struct panel)); + + new = &det->panels[det->n_panels-1]; + memcpy(new, &det->defaults, sizeof(struct panel)); + + /* Create a new copy of the camera length location if needed */ + if ( new->clen_from != NULL ) { + new->clen_from = strdup(new->clen_from); + } + strcpy(new->name, name); + + return new; +} + + +static struct badregion *new_bad_region(struct detector *det, const char *name) +{ + struct badregion *new; + + det->n_bad++; + det->bad = realloc(det->bad, det->n_bad*sizeof(struct badregion)); + + new = &det->bad[det->n_bad-1]; + new->min_x = NAN; + new->max_x = NAN; + new->min_y = NAN; + new->max_y = NAN; + strcpy(new->name, name); + + return new; +} + + +struct panel *find_panel_by_name(struct detector *det, const char *name) +{ + int i; + + for ( i=0; in_panels; i++ ) { + if ( strcmp(det->panels[i].name, name) == 0 ) { + return &det->panels[i]; + } + } + + return NULL; +} + + +static struct badregion *find_bad_region_by_name(struct detector *det, + const char *name) +{ + int i; + + for ( i=0; in_bad; i++ ) { + if ( strcmp(det->bad[i].name, name) == 0 ) { + return &det->bad[i]; + } + } + + return NULL; +} + + +static char *find_or_add_rg(struct detector *det, const char *name) +{ + int i; + char **new; + char *tmp; + + for ( i=0; inum_rigid_groups; i++ ) { + + if ( strcmp(det->rigid_groups[i], name) == 0 ) { + return det->rigid_groups[i]; + } + + } + + new = realloc(det->rigid_groups, + (1+det->num_rigid_groups)*sizeof(char *)); + if ( new == NULL ) return NULL; + + det->rigid_groups = new; + + tmp = strdup(name); + det->rigid_groups[det->num_rigid_groups] = tmp; + + det->num_rigid_groups++; + + return tmp; +} + + +static int parse_field_for_panel(struct panel *panel, const char *key, + const char *val, struct detector *det) +{ + int reject = 0; + + if ( strcmp(key, "min_fs") == 0 ) { + panel->min_fs = atof(val); + } else if ( strcmp(key, "max_fs") == 0 ) { + panel->max_fs = atof(val); + } else if ( strcmp(key, "min_ss") == 0 ) { + panel->min_ss = atof(val); + } else if ( strcmp(key, "max_ss") == 0 ) { + panel->max_ss = atof(val); + } else if ( strcmp(key, "corner_x") == 0 ) { + panel->cnx = atof(val); + } else if ( strcmp(key, "corner_y") == 0 ) { + panel->cny = atof(val); + } else if ( strcmp(key, "rigid_group") == 0 ) { + panel->rigid_group = find_or_add_rg(det, val); + } else if ( strcmp(key, "clen") == 0 ) { + + char *end; + double v = strtod(val, &end); + if ( end == val ) { + /* This means "fill in later" */ + panel->clen = -1.0; + panel->clen_from = strdup(val); + } else { + panel->clen = v; + panel->clen_from = NULL; + } + + } else if ( strcmp(key, "coffset") == 0) { + panel->coffset = atof(val); + } else if ( strcmp(key, "res") == 0 ) { + panel->res = atof(val); + } else if ( strcmp(key, "peak_sep") == 0 ) { + panel->peak_sep = atof(val); + } else if ( strcmp(key, "integr_radius") == 0 ) { + panel->integr_radius = atof(val); + } else if ( strcmp(key, "badrow_direction") == 0 ) { + panel->badrow = val[0]; /* First character only */ + if ( (panel->badrow != 'x') && (panel->badrow != 'y') + && (panel->badrow != 'f') && (panel->badrow != 's') + && (panel->badrow != '-') ) { + ERROR("badrow_direction must be x, y, f, s or '-'\n"); + ERROR("Assuming '-'\n."); + panel->badrow = '-'; + } + if ( panel->badrow == 'x' ) panel->badrow = 'f'; + if ( panel->badrow == 'y' ) panel->badrow = 's'; + } else if ( strcmp(key, "no_index") == 0 ) { + panel->no_index = atob(val); + } else if ( strcmp(key, "fs") == 0 ) { + if ( dir_conv(val, &panel->fsx, &panel->fsy) != 0 ) { + ERROR("Invalid fast scan direction '%s'\n", val); + reject = 1; + } + } else if ( strcmp(key, "ss") == 0 ) { + if ( dir_conv(val, &panel->ssx, &panel->ssy) != 0 ) { + ERROR("Invalid slow scan direction '%s'\n", val); + reject = 1; + } + } else { + ERROR("Unrecognised field '%s'\n", key); + } + + return reject; +} + + +static int parse_field_bad(struct badregion *panel, const char *key, + const char *val) +{ + int reject = 0; + + if ( strcmp(key, "min_x") == 0 ) { + panel->min_x = atof(val); + } else if ( strcmp(key, "max_x") == 0 ) { + panel->max_x = atof(val); + } else if ( strcmp(key, "min_y") == 0 ) { + panel->min_y = atof(val); + } else if ( strcmp(key, "max_y") == 0 ) { + panel->max_y = atof(val); + } else { + ERROR("Unrecognised field '%s'\n", key); + } + + return reject; +} + + +static void parse_toplevel(struct detector *det, const char *key, + const char *val) +{ + if ( strcmp(key, "mask") == 0 ) { + + det->mask = strdup(val); + + } else if ( strcmp(key, "mask_bad") == 0 ) { + + char *end; + double v = strtod(val, &end); + + if ( end != val ) { + det->mask_bad = v; + } + + } else if ( strcmp(key, "mask_good") == 0 ) { + + char *end; + double v = strtod(val, &end); + + if ( end != val ) { + det->mask_good = v; + } + + } else if ( strcmp(key, "peak_sep") == 0 ) { + det->defaults.peak_sep = atof(val); + } else if ( strcmp(key, "integr_radius") == 0 ) { + det->defaults.integr_radius = atof(val); + } else if ( parse_field_for_panel(&det->defaults, key, val, det) ) { + ERROR("Unrecognised top level field '%s'\n", key); + } +} + + +struct detector *get_detector_geometry(const char *filename) +{ + FILE *fh; + struct detector *det; + char *rval; + char **bits; + int i; + int reject = 0; + int x, y, max_fs, max_ss; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + det = calloc(1, sizeof(struct detector)); + if ( det == NULL ) { + fclose(fh); + return NULL; + } + + det->n_panels = 0; + det->panels = NULL; + det->n_bad = 0; + det->bad = NULL; + det->mask_good = 0; + det->mask_bad = 0; + det->mask = NULL; + det->num_rigid_groups = 0; + det->rigid_groups = NULL; + + /* The default defaults... */ + det->defaults.min_fs = -1; + det->defaults.min_ss = -1; + det->defaults.max_fs = -1; + det->defaults.max_ss = -1; + det->defaults.cnx = NAN; + det->defaults.cny = NAN; + det->defaults.clen = -1.0; + det->defaults.coffset = 0.0; + det->defaults.res = -1.0; + det->defaults.badrow = '-'; + det->defaults.no_index = 0; + det->defaults.peak_sep = 50.0; + det->defaults.integr_radius = 3.0; + det->defaults.fsx = 1.0; + det->defaults.fsy = 0.0; + det->defaults.ssx = 0.0; + det->defaults.ssy = 1.0; + det->defaults.rigid_group = NULL; + strncpy(det->defaults.name, "", 1023); + + do { + + int n1, n2; + char **path; + char line[1024]; + struct badregion *badregion = NULL; + struct panel *panel = NULL; + char *key; + char wholeval[1024]; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + + if ( line[0] == ';' ) continue; + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; in_panels == -1 ) { + ERROR("No panel descriptions in geometry file.\n"); + fclose(fh); + free(det); + return NULL; + } + + max_fs = 0; + max_ss = 0; + for ( i=0; in_panels; i++ ) { + + if ( det->panels[i].min_fs < 0 ) { + ERROR("Please specify the minimum FS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].max_fs < 0 ) { + ERROR("Please specify the maximum FS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].min_ss < 0 ) { + ERROR("Please specify the minimum SS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].max_ss < 0 ) { + ERROR("Please specify the maximum SS coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( isnan(det->panels[i].cnx) ) { + ERROR("Please specify the corner X coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( isnan(det->panels[i].cny) ) { + ERROR("Please specify the corner Y coordinate for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( (det->panels[i].clen < 0.0) + && (det->panels[i].clen_from == NULL) ) { + ERROR("Please specify the camera length for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + if ( det->panels[i].res < 0 ) { + ERROR("Please specify the resolution for" + " panel %s\n", det->panels[i].name); + reject = 1; + } + /* It's OK if the badrow direction is '0' */ + /* It's not a problem if "no_index" is still zero */ + /* The default peak_sep is OK (maybe) */ + /* The default transformation matrix is at least valid */ + + if ( det->panels[i].max_fs > max_fs ) { + max_fs = det->panels[i].max_fs; + } + if ( det->panels[i].max_ss > max_ss ) { + max_ss = det->panels[i].max_ss; + } + + } + + for ( i=0; in_bad; i++ ) { + + if ( isnan(det->bad[i].min_x) ) { + ERROR("Please specify the minimum x coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].min_y) ) { + ERROR("Please specify the minimum y coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].max_x) ) { + ERROR("Please specify the maximum x coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + if ( isnan(det->bad[i].max_y) ) { + ERROR("Please specify the maximum y coordinate for" + " bad region %s\n", det->bad[i].name); + reject = 1; + } + } + + for ( x=0; x<=max_fs; x++ ) { + for ( y=0; y<=max_ss; y++ ) { + if ( find_panel(det, x, y) == NULL ) { + ERROR("Detector geometry invalid: contains gaps.\n"); + reject = 1; + goto out; + } + } + } + +out: + det->max_fs = max_fs; + det->max_ss = max_ss; + + /* Calculate matrix inverse */ + for ( i=0; in_panels; i++ ) { + + struct panel *p; + double d; + + p = &det->panels[i]; + + if ( p->fsx*p->ssy == p->ssx*p->fsy ) { + ERROR("Panel %i transformation singular.\n", i); + reject = 1; + } + + d = (double)p->fsx*p->ssy - p->ssx*p->fsy; + p->xfs = p->ssy / d; + p->yfs = -p->ssx / d; + p->xss = -p->fsy / d; + p->yss = p->fsx / d; + + } + + if ( reject ) return NULL; + + fclose(fh); + + return det; +} + + +void free_detector_geometry(struct detector *det) +{ + int i; + + for ( i=0; inum_rigid_groups; i++ ) { + free(det->rigid_groups[i]); + } + free(det->rigid_groups); + + free(det->panels); + free(det->bad); + free(det->mask); + free(det); +} + + +struct detector *copy_geom(const struct detector *in) +{ + struct detector *out; + int i; + + out = malloc(sizeof(struct detector)); + memcpy(out, in, sizeof(struct detector)); + + if ( in->mask != NULL ) { + out->mask = strdup(in->mask); + } else { + out->mask = NULL; /* = in->mask */ + } + + out->panels = malloc(out->n_panels * sizeof(struct panel)); + memcpy(out->panels, in->panels, out->n_panels * sizeof(struct panel)); + + out->bad = malloc(out->n_bad * sizeof(struct badregion)); + memcpy(out->bad, in->bad, out->n_bad * sizeof(struct badregion)); + + if ( in->rigid_groups != NULL ) { + + out->rigid_groups = malloc(out->num_rigid_groups*sizeof(char *)); + memcpy(out->rigid_groups, in->rigid_groups, + out->num_rigid_groups*sizeof(char *)); + + for ( i=0; inum_rigid_groups; i++ ) { + out->rigid_groups[i] = strdup(in->rigid_groups[i]); + } + + } + + for ( i=0; in_panels; i++ ) { + + struct panel *p; + + p = &out->panels[i]; + + if ( p->clen_from != NULL ) { + /* Make a copy of the clen_from fields unique to this + * copy of the structure. */ + p->clen_from = strdup(p->clen_from); + } + + } + + for ( i=0; inum_rigid_groups; i++ ) { + + int j; + char *rg = in->rigid_groups[i]; + char *rgn = out->rigid_groups[i]; + + for ( j=0; jn_panels; j++ ) { + + if ( in->panels[j].rigid_group == rg ) { + out->panels[j].rigid_group = rgn; + } + + } + + } + + return out; +} + + +struct detector *simple_geometry(const struct image *image) +{ + struct detector *geom; + + geom = calloc(1, sizeof(struct detector)); + + geom->n_panels = 1; + geom->panels = calloc(1, sizeof(struct panel)); + + geom->panels[0].min_fs = 0; + geom->panels[0].max_fs = image->width-1; + geom->panels[0].min_ss = 0; + geom->panels[0].max_ss = image->height-1; + geom->panels[0].cnx = -image->width / 2.0; + geom->panels[0].cny = -image->height / 2.0; + geom->panels[0].rigid_group = NULL; + + geom->panels[0].fsx = 1; + geom->panels[0].fsy = 0; + geom->panels[0].ssx = 0; + geom->panels[0].ssy = 1; + + geom->panels[0].xfs = 1; + geom->panels[0].xss = 0; + geom->panels[0].yfs = 0; + geom->panels[0].yss = 1; + + return geom; +} + + +int reverse_2d_mapping(double x, double y, double *pfs, double *pss, + struct detector *det) +{ + int i; + + for ( i=0; in_panels; i++ ) { + + struct panel *p = &det->panels[i]; + double cx, cy, fs, ss; + + /* Get position relative to corner */ + cx = x - p->cnx; + cy = y - p->cny; + + /* Reverse the transformation matrix */ + fs = cx*p->xfs + cy*p->yfs; + ss = cx*p->xss + cy*p->yss; + + /* In range? */ + if ( fs < 0 ) continue; + if ( ss < 0 ) continue; + if ( fs > (p->max_fs-p->min_fs+1) ) continue; + if ( ss > (p->max_ss-p->min_ss+1) ) continue; + + *pfs = fs + p->min_fs; + *pss = ss + p->min_ss; + return 0; + + } + + return 1; +} + + +static void check_extents(struct panel p, double *min_x, double *min_y, + double *max_x, double *max_y, double fs, double ss) +{ + double xs, ys, rx, ry; + + xs = fs*p.fsx + ss*p.ssx; + ys = fs*p.fsy + ss*p.ssy; + + rx = xs + p.cnx; + ry = ys + p.cny; + + if ( rx > *max_x ) *max_x = rx; + if ( ry > *max_y ) *max_y = ry; + if ( rx < *min_x ) *min_x = rx; + if ( ry < *min_y ) *min_y = ry; +} + + +double largest_q(struct image *image) +{ + int fs, ss; + double ttm = 0.0; + double qmax = 0.0; + + for ( fs=0; fswidth; fs++ ) { + for ( ss=0; ssheight; 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; + } + + } + } + + return qmax; +} + + +double smallest_q(struct image *image) +{ + int fs, ss; + double ttm = +INFINITY; + double qmin = +INFINITY; + for ( fs=0; fswidth; fs++ ) { + for ( ss=0; ssheight; 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; + } + + } + } + + return qmin; +} + + +void get_pixel_extents(struct detector *det, + double *min_x, double *min_y, + double *max_x, double *max_y) +{ + int i; + + *min_x = 0.0; + *max_x = 0.0; + *min_y = 0.0; + *max_y = 0.0; + + /* To determine the maximum extents of the detector, put all four + * corners of each panel through the transformations and watch for the + * biggest */ + + for ( i=0; in_panels; i++ ) { + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + 0.0, + 0.0); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + 0.0, + det->panels[i].max_ss-det->panels[i].min_ss+1); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + det->panels[i].max_fs-det->panels[i].min_fs+1, + 0.0); + + check_extents(det->panels[i], min_x, min_y, max_x, max_y, + det->panels[i].max_fs-det->panels[i].min_fs+1, + det->panels[i].max_ss-det->panels[i].min_ss+1); + + + } +} + + +int write_detector_geometry(const char *filename, struct detector *det) +{ + struct panel *p; + int pi; + FILE *fh; + + if ( filename == NULL ) return 2; + if ( det->n_panels < 1 ) return 3; + + fh = fopen(filename, "w"); + if ( fh == NULL ) return 1; + + for ( pi=0; pin_panels; pi++) { + + p = &(det->panels[pi]); + + if ( p == NULL ) return 4; + + if ( pi > 0 ) fprintf(fh, "\n"); + + fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs); + fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss); + fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs); + fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss); + fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow); + fprintf(fh, "%s/res = %g\n", p->name, p->res); + fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep); + fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from); + fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy); + fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy); + fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx); + fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny); + + if ( p->no_index ) { + fprintf(fh, "%s/no_index = 1\n", p->name); + } /* else don't clutter up the file */ + + if ( p->rigid_group != NULL ) { + fprintf(fh, "%s/rigid_group = %s\n", + p->name, p->rigid_group); + } + + } + fclose(fh); + + return 0; +} diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h new file mode 100644 index 00000000..dd5dccec --- /dev/null +++ b/libcrystfel/src/detector.h @@ -0,0 +1,131 @@ +/* + * detector.h + * + * Detector properties + * + * (c) 2006-2011 Thomas White + * (c) 2011 Rick Kirian + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef DETECTOR_H +#define DETECTOR_H + +struct image; +struct hdfile; + +#include "hdf5-file.h" +#include "image.h" + + +struct panel +{ + char name[1024]; /* Name for this panel */ + + int min_fs; /* Smallest FS value considered to be in the panel */ + int max_fs; /* Largest FS value considered to be in this panel */ + int min_ss; /* ... and so on */ + int max_ss; + double cnx; /* Location of corner (min_fs,min_ss) in pixels */ + double cny; + double coffset; + double clen; /* Camera length in metres */ + char *clen_from; + double res; /* Resolution in pixels per metre */ + char badrow; /* 'x' or 'y' */ + int no_index; /* Don't index peaks in this panel if non-zero */ + double peak_sep; /* Characteristic peak separation */ + double integr_radius; /* Peak integration radius */ + char *rigid_group; /* Rigid group, or -1 for none */ + + double fsx; + double fsy; + double ssx; + double ssy; + + double xfs; + double yfs; + double xss; + double yss; +}; + + +struct badregion +{ + char name[1024]; + double min_x; + double max_x; + double min_y; + double max_y; +}; + + +struct detector +{ + struct panel *panels; + int n_panels; + + int max_fs; + int max_ss; /* Size of overall array needed, minus 1 */ + + struct badregion *bad; + int n_bad; + + char *mask; + unsigned int mask_bad; + unsigned int mask_good; + + char **rigid_groups; + int num_rigid_groups; + + struct panel defaults; +}; + + +extern struct rvec get_q(struct image *image, double fs, double ss, + double *ttp, double k); + +extern struct rvec get_q_for_panel(struct panel *p, double fs, double ss, + double *ttp, double k); + +extern double get_tt(struct image *image, double xs, double ys); + +extern int in_bad_region(struct detector *det, double fs, double ss); + +extern void record_image(struct image *image, int do_poisson); + +extern struct panel *find_panel(struct detector *det, double fs, double ss); + +extern int find_panel_number(struct detector *det, int fs, int ss); + +extern struct detector *get_detector_geometry(const char *filename); + +extern void free_detector_geometry(struct detector *det); + +extern struct detector *simple_geometry(const struct image *image); + +extern void get_pixel_extents(struct detector *det, + double *min_x, double *min_y, + double *max_x, double *max_y); + +extern void fill_in_values(struct detector *det, struct hdfile *f); + +extern struct detector *copy_geom(const struct detector *in); + +extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss, + struct detector *det); + +extern double largest_q(struct image *image); + +extern double smallest_q(struct image *image); + +extern struct panel *find_panel_by_name(struct detector *det, const char *name); + + +#endif /* DETECTOR_H */ diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c new file mode 100644 index 00000000..355e97f1 --- /dev/null +++ b/libcrystfel/src/hdf5-file.c @@ -0,0 +1,817 @@ +/* + * hdf5-file.c + * + * Read/write HDF5 data files + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "image.h" +#include "hdf5-file.h" +#include "utils.h" + + +struct hdfile { + + const char *path; /* Current data path */ + + size_t nx; /* Image width */ + size_t ny; /* Image height */ + + hid_t fh; /* HDF file handle */ + hid_t dh; /* Dataset handle */ + + int data_open; /* True if dh is initialised */ +}; + + +struct hdfile *hdfile_open(const char *filename) +{ + struct hdfile *f; + + f = malloc(sizeof(struct hdfile)); + if ( f == NULL ) return NULL; + + /* Please stop spamming my terminal */ + H5Eset_auto2(H5E_DEFAULT, NULL, NULL); + + f->fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + if ( f->fh < 0 ) { + ERROR("Couldn't open file: %s\n", filename); + free(f); + return NULL; + } + + f->data_open = 0; + + return f; +} + + +int hdfile_set_image(struct hdfile *f, const char *path) +{ + hsize_t size[2]; + hsize_t max_size[2]; + hid_t sh; + + f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT); + if ( f->dh < 0 ) { + ERROR("Couldn't open dataset\n"); + return -1; + } + f->data_open = 1; + + sh = H5Dget_space(f->dh); + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + ERROR("Dataset is not two-dimensional\n"); + return -1; + } + H5Sget_simple_extent_dims(sh, size, max_size); + H5Sclose(sh); + + f->nx = size[0]; + f->ny = size[1]; + + return 0; +} + + +int get_peaks(struct image *image, struct hdfile *f, const char *p) +{ + hid_t dh, sh; + hsize_t size[2]; + hsize_t max_size[2]; + int i; + float *buf; + herr_t r; + int tw; + + dh = H5Dopen2(f->fh, p, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Peak list (%s) not found.\n", p); + return 1; + } + + sh = H5Dget_space(dh); + if ( sh < 0 ) { + H5Dclose(dh); + ERROR("Couldn't get dataspace for peak list.\n"); + return 1; + } + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + ERROR("Peak list has the wrong dimensionality (%i).\n", + H5Sget_simple_extent_ndims(sh)); + H5Sclose(sh); + H5Dclose(dh); + return 1; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + tw = size[1]; + if ( (tw != 3) && (tw != 4) ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Peak list has the wrong dimensions.\n"); + return 1; + } + + buf = malloc(sizeof(float)*size[0]*size[1]); + if ( buf == NULL ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Couldn't reserve memory for the peak list.\n"); + return 1; + } + r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read peak list.\n"); + free(buf); + return 1; + } + + if ( image->features != NULL ) { + image_feature_list_free(image->features); + } + image->features = image_feature_list_new(); + + for ( i=0; idet, fs, ss); + if ( p == NULL ) continue; + if ( p->no_index ) continue; + + image_add_feature(image->features, fs, ss, image, val, NULL); + + } + + free(buf); + H5Sclose(sh); + H5Dclose(dh); + + return 0; +} + + +static void cleanup(hid_t fh) +{ + int n_ids, i; + hid_t ids[256]; + + n_ids = H5Fget_obj_ids(fh, H5F_OBJ_ALL, 256, ids); + for ( i=0; idata_open ) { + H5Dclose(f->dh); + } + cleanup(f->fh); + + H5Fclose(f->fh); + free(f); +} + + +int hdf5_write(const char *filename, const void *data, + int width, int height, int type) +{ + hid_t fh, gh, sh, dh; /* File, group, dataspace and data handles */ + hid_t ph; /* Property list */ + herr_t r; + hsize_t size[2]; + hsize_t max_size[2]; + + fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if ( fh < 0 ) { + ERROR("Couldn't create file: %s\n", filename); + return 1; + } + + gh = H5Gcreate2(fh, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + if ( gh < 0 ) { + ERROR("Couldn't create group\n"); + H5Fclose(fh); + return 1; + } + + /* Note the "swap" here, according to section 3.2.5, + * "C versus Fortran Dataspaces", of the HDF5 user's guide. */ + size[0] = height; + size[1] = width; + max_size[0] = height; + max_size[1] = width; + sh = H5Screate_simple(2, size, max_size); + + /* Set compression */ + ph = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_chunk(ph, 2, size); + H5Pset_deflate(ph, 3); + + dh = H5Dcreate2(gh, "data", type, sh, + H5P_DEFAULT, ph, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't create dataset\n"); + H5Fclose(fh); + return 1; + } + + /* Muppet check */ + H5Sget_simple_extent_dims(sh, size, max_size); + + r = H5Dwrite(dh, type, H5S_ALL, + H5S_ALL, H5P_DEFAULT, data); + if ( r < 0 ) { + ERROR("Couldn't write data\n"); + H5Dclose(dh); + H5Fclose(fh); + return 1; + } + + H5Pclose(ph); + H5Gclose(gh); + H5Dclose(dh); + H5Fclose(fh); + + return 0; +} + + +static double get_wavelength(struct hdfile *f) +{ + herr_t r; + hid_t dh; + double lambda; + int nm = 1; + + dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_nm", H5P_DEFAULT); + if ( dh < 0 ) { + dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_A", H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't get photon wavelength from HDF5 file.\n"); + return -1.0; + } + nm = 0; + } + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &lambda); + H5Dclose(dh); + + if ( r < 0 ) return -1.0; + if ( isnan(lambda) ) return -1.0; + + /* Convert nm -> m */ + if ( nm ) return lambda / 1.0e9; + return lambda / 1.0e10; +} + + +static double get_i0(struct hdfile *f) +{ + herr_t r; + hid_t dh; + double i0; + + dh = H5Dopen2(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT); + if ( dh < 0 ) return -1.0; + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &i0); + H5Dclose(dh); + if ( r < 0 ) return -1.0; + + return i0; +} + + +static void debodge_saturation(struct hdfile *f, struct image *image) +{ + hid_t dh, sh; + hsize_t size[2]; + hsize_t max_size[2]; + int i; + float *buf; + herr_t r; + + dh = H5Dopen2(f->fh, "/processing/hitfinder/peakinfo_saturated", + H5P_DEFAULT); + + if ( dh < 0 ) { + /* This isn't an error */ + return; + } + + sh = H5Dget_space(dh); + if ( sh < 0 ) { + H5Dclose(dh); + ERROR("Couldn't get dataspace for saturation table.\n"); + return; + } + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + H5Sclose(sh); + H5Dclose(dh); + return; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + if ( size[1] != 3 ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Saturation table has the wrong dimensions.\n"); + return; + } + + buf = malloc(sizeof(float)*size[0]*size[1]); + if ( buf == NULL ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Couldn't reserve memory for saturation table.\n"); + return; + } + r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read saturation table.\n"); + free(buf); + return; + } + + for ( i=0; idata[x+image->width*y] = val / 5.0; + image->data[x+1+image->width*y] = val / 5.0; + image->data[x-1+image->width*y] = val / 5.0; + image->data[x+image->width*(y+1)] = val / 5.0; + image->data[x+image->width*(y-1)] = val / 5.0; + + } + + free(buf); + H5Sclose(sh); + H5Dclose(dh); +} + + +int hdf5_read(struct hdfile *f, struct image *image, int satcorr) +{ + herr_t r; + float *buf; + uint16_t *flags; + hid_t mask_dh; + + /* Note the "swap" here, according to section 3.2.5, + * "C versus Fortran Dataspaces", of the HDF5 user's guide. */ + image->width = f->ny; + image->height = f->nx; + + buf = malloc(sizeof(float)*f->nx*f->ny); + + r = H5Dread(f->dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read data\n"); + free(buf); + return 1; + } + image->data = buf; + + if ( (image->det != NULL) && (image->det->mask != NULL) ) { + + mask_dh = H5Dopen2(f->fh, image->det->mask, H5P_DEFAULT); + if ( mask_dh <= 0 ) { + ERROR("Couldn't open flags\n"); + image->flags = NULL; + } else { + flags = malloc(sizeof(uint16_t)*f->nx*f->ny); + r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL, + H5P_DEFAULT, flags); + if ( r < 0 ) { + ERROR("Couldn't read flags\n"); + free(flags); + image->flags = NULL; + } else { + image->flags = flags; + } + H5Dclose(mask_dh); + } + + } + + /* Read wavelength from file */ + image->lambda = get_wavelength(f); + + image->i0 = get_i0(f); + if ( image->i0 < 0.0 ) { + ERROR("Couldn't read incident intensity - using 1.0.\n"); + image->i0 = 1.0; + image->i0_available = 0; + } else { + image->i0_available = 1; + } + + if ( satcorr ) debodge_saturation(f, image); + + return 0; +} + + +static int looks_like_image(hid_t h) +{ + hid_t sh; + hsize_t size[2]; + hsize_t max_size[2]; + + sh = H5Dget_space(h); + if ( sh < 0 ) return 0; + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + return 0; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + if ( ( size[0] > 64 ) && ( size[1] > 64 ) ) return 1; + + return 0; +} + + +double get_value(struct hdfile *f, const char *name) +{ + hid_t dh; + hid_t sh; + hsize_t size; + hsize_t max_size; + hid_t type; + hid_t class; + herr_t r; + double buf; + + dh = H5Dopen2(f->fh, name, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't open data\n"); + return 0.0; + } + + type = H5Dget_type(dh); + class = H5Tget_class(type); + + if ( class != H5T_FLOAT ) { + ERROR("Not a floating point value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + sh = H5Dget_space(dh); + if ( H5Sget_simple_extent_ndims(sh) != 1 ) { + ERROR("Not a scalar value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + H5Sget_simple_extent_dims(sh, &size, &max_size); + if ( size != 1 ) { + ERROR("Not a scalar value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) { + ERROR("Couldn't read value.\n"); + H5Tclose(type); + H5Dclose(dh); + return 0.0; + } + + return buf; +} + + +struct copy_hdf5_field +{ + char **fields; + int n_fields; + int max_fields; +}; + + +struct copy_hdf5_field *new_copy_hdf5_field_list() +{ + struct copy_hdf5_field *n; + + n = calloc(1, sizeof(struct copy_hdf5_field)); + if ( n == NULL ) return NULL; + + n->max_fields = 32; + n->fields = malloc(n->max_fields*sizeof(char *)); + if ( n->fields == NULL ) { + free(n); + return NULL; + } + + return n; +} + + +void free_copy_hdf5_field_list(struct copy_hdf5_field *n) +{ + int i; + for ( i=0; in_fields; i++ ) { + free(n->fields[i]); + } + free(n->fields); + free(n); +} + + +void add_copy_hdf5_field(struct copy_hdf5_field *copyme, + const char *name) +{ + assert(copyme->n_fields < copyme->max_fields); /* FIXME */ + + copyme->fields[copyme->n_fields] = strdup(name); + if ( copyme->fields[copyme->n_fields] == NULL ) { + ERROR("Failed to add field for copying '%s'\n", name); + return; + } + + copyme->n_fields++; +} + + +void copy_hdf5_fields(struct hdfile *f, const struct copy_hdf5_field *copyme, + FILE *fh) +{ + int i; + + if ( copyme == NULL ) return; + + for ( i=0; in_fields; i++ ) { + + char *val; + char *field; + + field = copyme->fields[i]; + val = hdfile_get_string_value(f, field); + + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + + free(val); + + } +} + + +char *hdfile_get_string_value(struct hdfile *f, const char *name) +{ + hid_t dh; + hid_t sh; + hsize_t size; + hsize_t max_size; + hid_t type; + hid_t class; + + dh = H5Dopen2(f->fh, name, H5P_DEFAULT); + if ( dh < 0 ) return NULL; + + type = H5Dget_type(dh); + class = H5Tget_class(type); + + if ( class == H5T_STRING ) { + + herr_t r; + char *tmp; + hid_t sh; + + size = H5Tget_size(type); + tmp = malloc(size+1); + + sh = H5Screate(H5S_SCALAR); + + r = H5Dread(dh, type, sh, sh, H5P_DEFAULT, tmp); + if ( r < 0 ) goto fail; + + /* Two possibilities: + * String is already zero-terminated + * String is not terminated. + * Make sure things are done properly... */ + tmp[size] = '\0'; + chomp(tmp); + + return tmp; + + } + + sh = H5Dget_space(dh); + if ( H5Sget_simple_extent_ndims(sh) != 1 ) goto fail; + + H5Sget_simple_extent_dims(sh, &size, &max_size); + if ( size != 1 ) { + H5Dclose(dh); + goto fail; + } + + switch ( class ) { + case H5T_FLOAT : { + + herr_t r; + double buf; + char *tmp; + + r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) goto fail; + + tmp = malloc(256); + snprintf(tmp, 255, "%f", buf); + + return tmp; + + } + case H5T_INTEGER : { + + herr_t r; + int buf; + char *tmp; + + r = H5Dread(dh, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &buf); + if ( r < 0 ) goto fail; + + tmp = malloc(256); + snprintf(tmp, 255, "%d", buf); + + return tmp; + } + default : { + goto fail; + } + } + +fail: + H5Tclose(type); + H5Dclose(dh); + return NULL; +} + + +char **hdfile_read_group(struct hdfile *f, int *n, const char *parent, + int **p_is_group, int **p_is_image) +{ + hid_t gh; + hsize_t num; + char **res; + int i; + int *is_group; + int *is_image; + H5G_info_t ginfo; + + gh = H5Gopen2(f->fh, parent, H5P_DEFAULT); + if ( gh < 0 ) { + *n = 0; + return NULL; + } + + if ( H5Gget_info(gh, &ginfo) < 0 ) { + /* Whoopsie */ + *n = 0; + return NULL; + } + num = ginfo.nlinks; + *n = num; + if ( num == 0 ) return NULL; + + res = malloc(num*sizeof(char *)); + is_image = malloc(num*sizeof(int)); + is_group = malloc(num*sizeof(int)); + *p_is_image = is_image; + *p_is_group = is_group; + + for ( i=0; i 1 ) { + snprintf(res[i], 255, "%s/%s", parent, buf); + } else { + snprintf(res[i], 255, "%s%s", parent, buf); + } /* ick */ + + is_image[i] = 0; + is_group[i] = 0; + dh = H5Oopen(gh, buf, H5P_DEFAULT); + if ( dh < 0 ) continue; + type = H5Iget_type(dh); + + if ( type == H5I_GROUP ) { + is_group[i] = 1; + } else if ( type == H5I_DATASET ) { + is_image[i] = looks_like_image(dh); + } + H5Oclose(dh); + + } + + return res; +} + + +int hdfile_set_first_image(struct hdfile *f, const char *group) +{ + char **names; + int *is_group; + int *is_image; + int n, i, j; + + names = hdfile_read_group(f, &n, group, &is_group, &is_image); + if ( n == 0 ) return 1; + + for ( i=0; i + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef HDF5_H +#define HDF5_H + +#include +#include + +#include "image.h" + +struct hdfile; + +struct copy_hdf5_field; + + +extern int hdf5_write(const char *filename, const void *data, + int width, int height, int type); + +extern int hdf5_read(struct hdfile *f, struct image *image, int satcorr); + +extern struct hdfile *hdfile_open(const char *filename); +extern int hdfile_set_image(struct hdfile *f, const char *path); +extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile, + int binning, int16_t *maxp); +extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent, + int **p_is_group, int **p_is_image); +extern int hdfile_set_first_image(struct hdfile *f, const char *group); +extern void hdfile_close(struct hdfile *f); + +extern char *hdfile_get_string_value(struct hdfile *f, const char *name); +extern int get_peaks(struct image *image, struct hdfile *f, const char *p); +extern double get_value(struct hdfile *f, const char *name); + +extern struct copy_hdf5_field *new_copy_hdf5_field_list(void); +extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f); +extern void copy_hdf5_fields(struct hdfile *f, + const struct copy_hdf5_field *copyme, FILE *fh); +extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme, + const char *name); + + +#endif /* HDF5_H */ diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c new file mode 100644 index 00000000..75a8af0a --- /dev/null +++ b/libcrystfel/src/image.c @@ -0,0 +1,144 @@ +/* + * image.c + * + * Handle images and image features + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include +#include +#include +#include + +#include "image.h" +#include "utils.h" + +/** + * SECTION:image + * @short_description: Data structure representing an image + * @title: Image + * @section_id: + * @see_also: + * @include: "image.h" + * @Image: + * + * The image structure represents an image, usually one + * frame from a large series of diffraction patterns, which might be from the + * same or different crystals. + */ + + +struct _imagefeaturelist +{ + struct imagefeature *features; + int n_features; +}; + + +void image_add_feature(ImageFeatureList *flist, double fs, double ss, + struct image *parent, double intensity, const char *name) +{ + if ( flist->features ) { + flist->features = realloc(flist->features, + (flist->n_features+1) + *sizeof(struct imagefeature)); + } else { + assert(flist->n_features == 0); + flist->features = malloc(sizeof(struct imagefeature)); + } + + flist->features[flist->n_features].fs = fs; + flist->features[flist->n_features].ss = ss; + flist->features[flist->n_features].intensity = intensity; + flist->features[flist->n_features].parent = parent; + flist->features[flist->n_features].name = name; + flist->features[flist->n_features].valid = 1; + + flist->n_features++; + +} + + +ImageFeatureList *image_feature_list_new() +{ + ImageFeatureList *flist; + + flist = malloc(sizeof(ImageFeatureList)); + + flist->n_features = 0; + flist->features = NULL; + + return flist; +} + + +void image_feature_list_free(ImageFeatureList *flist) +{ + if ( !flist ) return; + + if ( flist->features ) free(flist->features); + free(flist); +} + + +struct imagefeature *image_feature_closest(ImageFeatureList *flist, + double fs, double ss, + double *d, int *idx) +{ + int i; + double dmin = +HUGE_VAL; + int closest = 0; + + for ( i=0; in_features; i++ ) { + + double ds; + + ds = distance(flist->features[i].fs, flist->features[i].ss, + fs, ss); + + if ( ds < dmin ) { + dmin = ds; + closest = i; + } + + } + + if ( dmin < +HUGE_VAL ) { + *d = dmin; + *idx = closest; + return &flist->features[closest]; + } + + *d = +INFINITY; + return NULL; +} + + +int image_feature_count(ImageFeatureList *flist) +{ + if ( flist == NULL ) return 0; + return flist->n_features; +} + + +struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx) +{ + /* Sanity check */ + if ( flist == NULL ) return NULL; + if ( idx > flist->n_features ) return NULL; + + if ( flist->features[idx].valid == 0 ) return NULL; + + return &flist->features[idx]; +} + + +void image_remove_feature(ImageFeatureList *flist, int idx) +{ + flist->features[idx].valid = 0; +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h new file mode 100644 index 00000000..4d4ac2d7 --- /dev/null +++ b/libcrystfel/src/image.h @@ -0,0 +1,184 @@ +/* + * image.h + * + * Handle images and image features + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef IMAGE_H +#define IMAGE_H + +#include +#include +#include + +#include "utils.h" +#include "cell.h" +#include "detector.h" +#include "reflist.h" + + +#define MAX_CELL_CANDIDATES (32) + + +/* Structure describing a feature in an image */ +struct imagefeature { + + struct image *parent; + double fs; + double ss; + double intensity; + + /* Reciprocal space coordinates (m^-1 of course) of this feature */ + double rx; + double ry; + double rz; + + /* Internal use only */ + int valid; + + const char *name; +}; + +/* An opaque type representing a list of image features */ +typedef struct _imagefeaturelist ImageFeatureList; + + +/** + * image: + * + * + * struct image + * { + * float *data; + * uint16_t *flags; + * double *twotheta; + * + * UnitCell *indexed_cell; + * UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; + * int ncells; + + * struct detector *det; + * struct beam_params *beam; + * char *filename; + * const struct copy_hdf5_field *copyme; + * + * int id; + * + * double m; + * + * double lambda; + * double div; + * double bw; + * double i0; + * int i0_available; + * double osf; + * double profile_radius; + * int pr_dud; + * + * int width; + * int height; + * + * RefList *reflections; + * + * ImageFeatureList *features; + * }; + * + * + * The field data contains the raw image data, if it + * is currently available. The data might be available throughout the + * processing of an experimental pattern, but it might not be available when + * simulating, scaling or merging patterns. Similarly, + * flags contains an array of the same dimensions + * as data to contain the bad pixel flags. + * twotheta likewise contains an array of 2*theta + * (scattering angle) values in radians, since these values are generated as a + * by-product of the scattering vector calculation and can be used later for + * calculating intensities from differential scattering cross sections. + * + * candidate_cells is an array of unit cells directly + * returned by the low-level indexing system. ncells + * is the number of candidate unit cells which were found. The maximum number + * of cells which may be returned is MAX_CELL_CANDIDATES. + * indexed_cell 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 NULL on failure. + * + * copyme represents a list of HDF5 fields to copy + * to the output stream. + **/ +struct image; + +struct image { + + float *data; + uint16_t *flags; + double *twotheta; + + UnitCell *indexed_cell; + UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; + int ncells; + + struct detector *det; + struct beam_params *beam; /* The nominal beam parameters */ + char *filename; + const struct copy_hdf5_field *copyme; + + 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 i0; /* Incident intensity */ + int i0_available; /* 0 if f0 wasn't available + * from the input. */ + double osf; /* Overall scaling factor */ + double profile_radius; /* Radius of reflection */ + int pr_dud; /* Post refinement failed */ + + int width; + int height; + + /* Integrated (or about-to-be-integrated) reflections */ + RefList *reflections; + + /* Detected peaks */ + ImageFeatureList *features; + +}; + + +/* Feature lists */ +extern ImageFeatureList *image_feature_list_new(void); + +extern void image_feature_list_free(ImageFeatureList *flist); + +extern void image_add_feature(ImageFeatureList *flist, double x, double y, + struct image *parent, double intensity, + const char *name); + +extern void image_remove_feature(ImageFeatureList *flist, int idx); + +extern struct imagefeature *image_feature_closest(ImageFeatureList *flist, + double fs, double ss, + double *d, int *idx); + +extern int image_feature_count(ImageFeatureList *flist); +extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); + +#endif /* IMAGE_H */ diff --git a/libcrystfel/src/list_tmp.h b/libcrystfel/src/list_tmp.h new file mode 100644 index 00000000..a524b2f9 --- /dev/null +++ b/libcrystfel/src/list_tmp.h @@ -0,0 +1,106 @@ +/* + * Template for creating indexed 3D lists of a given type, usually indexed + * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX. + * + * These are used, for example, for: + * - a list of 'double complex' values for storing structure factors, + * - a list of 'double' values for storing reflection intensities, + * - a list of 'unsigned int' values for counts of some sort. + * + * When LABEL and TYPE are #defined appropriately, including this header + * defines functions such as: + * - new_list_