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/hdf5-file.c | 817 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 817 insertions(+) create mode 100644 libcrystfel/src/hdf5-file.c (limited to 'libcrystfel/src/hdf5-file.c') 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