diff options
Diffstat (limited to 'libcrystfel/src/image.c')
| -rw-r--r-- | libcrystfel/src/image.c | 698 |
1 files changed, 695 insertions, 3 deletions
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 7bb4cede..aad5017c 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -3,12 +3,12 @@ * * Handle images and image features * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de> - * 2011-2016 Thomas White <taw@physics.org> + * 2011-2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -27,14 +27,28 @@ * */ -#define _ISOC99_SOURCE +#include <config.h> + #include <stdlib.h> #include <assert.h> #include <math.h> #include <stdio.h> +#include <hdf5.h> + +#ifdef HAVE_CBFLIB +#ifdef HAVE_CBFLIB_CBF_H +#include <cbflib/cbf.h> +#endif +#ifdef HAVE_CBF_CBF_H +#include <cbf/cbf.h> +#endif +#endif #include "image.h" #include "utils.h" +#include "events.h" +#include "hdf5-file.h" +#include "detector.h" /** * SECTION:image @@ -51,6 +65,14 @@ */ +struct imagefile +{ + enum imagefile_type type; + char *filename; + struct hdfile *hdfile; +}; + + struct _imagefeaturelist { struct imagefeature *features; @@ -300,3 +322,673 @@ void free_all_crystals(struct image *image) free(image->crystals); image->n_crystals = 0; } + + +/**************************** Image field lists *******************************/ + +struct imagefile_field_list +{ + char **fields; + int n_fields; + int max_fields; +}; + + +struct imagefile_field_list *new_imagefile_field_list() +{ + struct imagefile_field_list *n; + + n = calloc(1, sizeof(struct imagefile_field_list)); + 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_imagefile_field_list(struct imagefile_field_list *n) +{ + int i; + for ( i=0; i<n->n_fields; i++ ) { + free(n->fields[i]); + } + free(n->fields); + free(n); +} + + +void add_imagefile_field(struct imagefile_field_list *copyme, const char *name) +{ + int i; + + /* Already on the list? Don't re-add if so. */ + for ( i=0; i<copyme->n_fields; i++ ) { + if ( strcmp(copyme->fields[i], name) == 0 ) return; + } + + /* Need more space? */ + if ( copyme->n_fields == copyme->max_fields ) { + + char **nfields; + int nmax = copyme->max_fields + 32; + + nfields = realloc(copyme->fields, nmax*sizeof(char *)); + if ( nfields == NULL ) { + ERROR("Failed to allocate space for new HDF5 field.\n"); + return; + } + + copyme->max_fields = nmax; + copyme->fields = nfields; + + } + + 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++; +} + + +/******************************* CBF files ************************************/ + +#ifdef HAVE_CBFLIB + +static char *cbf_strerr(int e) +{ + char *err; + + err = malloc(1024); + if ( err == NULL ) return NULL; + + err[0] = '\0'; + + /* NB Sum of lengths of all strings must be less than 1024 */ + if ( e & CBF_FORMAT ) strcat(err, "Invalid format"); + if ( e & CBF_ALLOC ) strcat(err, "Memory allocation failed"); + if ( e & CBF_ARGUMENT ) strcat(err, "Invalid argument"); + if ( e & CBF_ASCII ) strcat(err, "Value is ASCII"); + if ( e & CBF_BINARY ) strcat(err, "Value is binary"); + if ( e & CBF_BITCOUNT ) strcat(err, "Wrong number of bits"); + if ( e & CBF_ENDOFDATA ) strcat(err, "End of data"); + if ( e & CBF_FILECLOSE ) strcat(err, "File close error"); + if ( e & CBF_FILEOPEN ) strcat(err, "File open error"); + if ( e & CBF_FILEREAD ) strcat(err, "File read error"); + if ( e & CBF_FILETELL ) strcat(err, "File tell error"); + if ( e & CBF_FILEWRITE ) strcat(err, "File write error"); + if ( e & CBF_IDENTICAL ) strcat(err, "Name already exists"); + if ( e & CBF_NOTFOUND ) strcat(err, "Not found"); + if ( e & CBF_OVERFLOW ) strcat(err, "Overflow"); + if ( e & CBF_UNDEFINED ) strcat(err, "Number undefined"); + if ( e & CBF_NOTIMPLEMENTED ) strcat(err, "Not implemented"); + + return err; +} + + +static int unpack_panels(struct image *image, signed int *data, int data_width, + int data_height) +{ + int pi; + + /* FIXME: Load these masks from an HDF5 file, if filenames are + * given in the geometry file */ + uint16_t *flags = NULL; + float *sat = NULL; + + image->dp = malloc(image->det->n_panels * sizeof(float *)); + image->bad = malloc(image->det->n_panels * sizeof(int *)); + image->sat = malloc(image->det->n_panels * sizeof(float *)); + if ( (image->dp == NULL) || (image->bad == NULL) + || (image->sat == NULL) ) + { + ERROR("Failed to allocate panels.\n"); + return 1; + } + + for ( pi=0; pi<image->det->n_panels; pi++ ) { + + struct panel *p; + int fs, ss; + + p = &image->det->panels[pi]; + image->dp[pi] = malloc(p->w*p->h*sizeof(float)); + image->bad[pi] = calloc(p->w*p->h, sizeof(int)); + image->sat[pi] = malloc(p->w*p->h*sizeof(float)); + if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) + || (image->sat[pi] == NULL) ) + { + ERROR("Failed to allocate panel\n"); + return 1; + } + + if ( p->mask != NULL ) { + ERROR("WARNING: Bad pixel masks do not currently work " + "with CBF files\n"); + ERROR(" (bad pixel regions specified in the geometry " + "file will be used, however)\n"); + } + + if ( p->satmap != NULL ) { + ERROR("WARNING: Saturation maps do not currently work " + "with CBF files\n"); + } + + if ( (p->orig_min_fs + p->w > data_width) + || (p->orig_min_ss + p->h > data_height) ) + { + ERROR("Panel %s is outside range of data in CBF file\n", + p->name); + return 1; + } + + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { + + int idx; + int cfs, css; + int bad = 0; + + cfs = fs+p->orig_min_fs; + css = ss+p->orig_min_ss; + idx = cfs + css*data_width; + + image->dp[pi][fs+p->w*ss] = data[idx]; + + if ( sat != NULL ) { + image->sat[pi][fs+p->w*ss] = sat[idx]; + } else { + image->sat[pi][fs+p->w*ss] = INFINITY; + } + + if ( p->no_index ) bad = 1; + + if ( in_bad_region(image->det, p, cfs, css) ) { + bad = 1; + } + + if ( flags != NULL ) { + + int f; + + f = flags[idx]; + + /* Bad if it's missing any of the "good" bits */ + if ( (f & image->det->mask_good) + != image->det->mask_good ) bad = 1; + + /* Bad if it has any of the "bad" bits. */ + if ( f & image->det->mask_bad ) bad = 1; + + } + image->bad[pi][fs+p->w*ss] = bad; + + } + } + + } + + return 0; +} + + +static void cbf_fill_in_beam_parameters(struct beam_params *beam, + struct imagefile *f, + struct image *image) +{ + double eV; + + if ( beam->photon_energy_from == NULL ) { + + /* Explicit value given */ + eV = beam->photon_energy; + + } else { + + ERROR("Can't get photon energy from CBF yet.\n"); + eV = 0.0; + + } + + image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale; +} + + +static void cbf_fill_in_clen(struct detector *det, struct imagefile *f) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + + struct panel *p = &det->panels[i]; + + if ( p->clen_from != NULL ) { + + ERROR("Can't get clen from CBF yet.\n"); + + } + + adjust_centering_for_rail(p); + + } +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + cbf_handle cbfh; + FILE *fh; + int r; + unsigned int compression; + int binary_id, minelement, maxelement, elsigned, elunsigned; + size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; + const char *byteorder; + signed int *data; + + if ( image->det == NULL ) { + ERROR("read_cbf() needs a geometry\n"); + return 1; + } + + if ( cbf_make_handle(&cbfh) ) { + ERROR("Failed to allocate CBF handle\n"); + return 1; + } + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file '%s'\n", f->filename); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + unpack_panels(image, data, dimfast, dimmid); + free(data); + + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + + cbf_free_handle(cbfh); + return 0; +} + + +static float *convert_float(signed int *data, int w, int h) +{ + float *df; + long int i; + + df = malloc(sizeof(float)*w*h); + if ( df == NULL ) return NULL; + + for ( i=0; i<w*h; i++ ) { + df[i] = data[i]; + } + + return df; +} + + +static int read_cbf_simple(struct imagefile *f, struct image *image) +{ + cbf_handle cbfh; + FILE *fh; + int r; + unsigned int compression; + int binary_id, minelement, maxelement, elsigned, elunsigned; + size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding; + const char *byteorder; + signed int *data; + + if ( cbf_make_handle(&cbfh) ) { + ERROR("Failed to allocate CBF handle\n"); + return 1; + } + + fh = fopen(f->filename, "rb"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", f->filename); + return 1; + } + /* CBFlib calls fclose(fh) when it's ready */ + + if ( cbf_read_widefile(cbfh, fh, 0) ) { + ERROR("Failed to read CBF file '%s'\n", f->filename); + cbf_free_handle(cbfh); + return 1; + } + + /* Select row 0 in data column inside array_data */ + cbf_find_category(cbfh, "array_data"); + cbf_find_column(cbfh, "data"); + cbf_select_row(cbfh, 0); + + /* Get parameters for array read */ + r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id, + &elsize, &elsigned, &elunsigned, + &elements, + &minelement, &maxelement, + &byteorder, + &dimfast, &dimmid, &dimslow, + &padding); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array parameters: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimslow != 0 ) { + ERROR("CBF data array is 3D - don't know what to do with it\n"); + cbf_free_handle(cbfh); + return 1; + } + + if ( dimfast*dimmid*elsize > 10e9 ) { + ERROR("CBF data is far too big (%i x %i x %i bytes).\n", + (int)dimfast, (int)dimmid, (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( elsize != 4 ) { + STATUS("Don't know what to do with element size %i\n", + (int)elsize); + cbf_free_handle(cbfh); + return 1; + } + + if ( strcmp(byteorder, "little_endian") != 0 ) { + STATUS("Don't know what to do with non-little-endian datan\n"); + cbf_free_handle(cbfh); + return 1; + } + + data = malloc(elsize*dimfast*dimmid); + if ( data == NULL ) { + ERROR("Failed to allocate memory for CBF data\n"); + cbf_free_handle(cbfh); + return 1; + } + + r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1, + elsize*dimfast*dimmid, &elread); + if ( r ) { + char *err = cbf_strerr(r); + ERROR("Failed to read CBF array: %s\n", err); + free(err); + cbf_free_handle(cbfh); + return 1; + } + + image->det = simple_geometry(image, dimfast, dimmid); + image->dp = malloc(sizeof(float *)); + if ( image->dp == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + image->dp[0] = convert_float(data, dimfast, dimmid); + if ( image->dp[0] == NULL ) { + ERROR("Failed to allocate dp array\n"); + return 1; + } + + if ( image->beam != NULL ) { + cbf_fill_in_beam_parameters(image->beam, f, image); + if ( image->lambda > 1000 ) { + ERROR("WARNING: Missing or nonsensical wavelength " + "(%e m) for %s.\n", + image->lambda, image->filename); + } + } + cbf_fill_in_clen(image->det, f); + fill_in_adu(image); + + cbf_free_handle(cbfh); + return 0; +} + +#else /* HAVE_CBFLIB */ + +static int read_cbf_simple(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +static int read_cbf(struct imagefile *f, struct image *image) +{ + ERROR("This version of CrystFEL was compiled without CBF support.\n"); + return 1; +} + + +#endif /* HAVE_CBFLIB */ + + +/****************************** Image files ***********************************/ + + +static signed int is_cbf_file(const char *filename) +{ + FILE *fh; + char line[1024]; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return -1; + + if ( fgets(line, 1024, fh) == NULL ) return -1; + fclose(fh); + + if ( strstr(line, "CBF") == NULL ) { + return 0; + } + + return 1; +} + + +struct imagefile *imagefile_open(const char *filename) +{ + struct imagefile *f; + + f = malloc(sizeof(struct imagefile)); + if ( f == NULL ) return NULL; + + if ( H5Fis_hdf5(filename) > 0 ) { + + /* This is an HDF5, pass through to HDF5 layer */ + f->type = IMAGEFILE_HDF5; + f->hdfile = hdfile_open(filename); + + if ( f->hdfile == NULL ) { + free(f); + return NULL; + } + + } else if ( is_cbf_file(filename) > 0 ) { + + f->type = IMAGEFILE_CBF; + + } + + f->filename = strdup(filename); + return f; +} + + +int imagefile_read(struct imagefile *f, struct image *image, + struct event *event) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read2(f->hdfile, image, event, 0); + } else if ( f->type == IMAGEFILE_CBF ) { + return read_cbf(f, image); + } else { + ERROR("Unknown file type %i\n", f->type); + return 1; + } +} + + +/* Read a simple file, no multi-event, no prior geometry etc, and + * generate a geometry for it */ +int imagefile_read_simple(struct imagefile *f, struct image *image) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + return hdf5_read(f->hdfile, image, NULL, 0); + } else { + return read_cbf_simple(f, image); + } +} + + +enum imagefile_type imagefile_get_type(struct imagefile *f) +{ + assert(f != NULL); + return f->type; +} + + +struct hdfile *imagefile_get_hdfile(struct imagefile *f) +{ + if ( f->type != IMAGEFILE_HDF5 ) { + ERROR("Not an HDF5 file!\n"); + return NULL; + } + + return f->hdfile; +} + + +void imagefile_copy_fields(struct imagefile *f, + const struct imagefile_field_list *copyme, + FILE *fh, struct event *ev) +{ + int i; + + if ( copyme == NULL ) return; + + for ( i=0; i<copyme->n_fields; i++ ) { + + char *val; + char *field; + + field = copyme->fields[i]; + + if ( f->type == IMAGEFILE_HDF5 ) { + val = hdfile_get_string_value(f->hdfile, field, ev); + if ( field[0] == '/' ) { + fprintf(fh, "hdf5%s = %s\n", field, val); + } else { + fprintf(fh, "hdf5/%s = %s\n", field, val); + } + free(val); + + } else { + STATUS("Mock CBF variable\n"); + fprintf(fh, "cbf/%s = %s\n", field, "(FIXME)"); + } + + + } +} + + +void imagefile_close(struct imagefile *f) +{ + if ( f->type == IMAGEFILE_HDF5 ) { + hdfile_close(f->hdfile); + } + free(f->filename); + free(f); +} |