diff options
Diffstat (limited to 'libcrystfel/src/image-hdf5.c')
-rw-r--r-- | libcrystfel/src/image-hdf5.c | 334 |
1 files changed, 334 insertions, 0 deletions
diff --git a/libcrystfel/src/image-hdf5.c b/libcrystfel/src/image-hdf5.c index 146918ee..653057eb 100644 --- a/libcrystfel/src/image-hdf5.c +++ b/libcrystfel/src/image-hdf5.c @@ -1670,3 +1670,337 @@ int is_hdf5_file(const char *filename) return ( (strcmp(ext, ".h5") == 0) || (strcmp(ext, ".cxi") == 0) ); } + + +/***************************** Writing *****************************/ + +struct hdf5_write_location { + + const char *location; + int n_panels; + int *panel_idxs; + + int max_ss; + int max_fs; + +}; + + +static void add_panel_to_location(struct hdf5_write_location *loc, + struct panel_template *p, int pi) +{ + int *new_panel_idxs; + + new_panel_idxs = realloc(loc->panel_idxs, + (loc->n_panels+1)*sizeof(int)); + if ( new_panel_idxs == NULL ) { + ERROR("Error while managing write location list.\n"); + return; + } + loc->panel_idxs = new_panel_idxs; + loc->panel_idxs[loc->n_panels] = pi; + loc->n_panels += 1; + if ( p->orig_max_fs > loc->max_fs ) { + loc->max_fs = p->orig_max_fs; + } + if ( p->orig_max_ss > loc->max_ss ) { + loc->max_ss = p->orig_max_ss; + } +} + + +static void add_panel_location(struct panel_template *p, + const char *p_location, int pi, + struct hdf5_write_location **plocations, + int *pnum_locations) +{ + int li; + int num_locations = *pnum_locations; + struct hdf5_write_location *locations = *plocations; + int done = 0; + + /* Does this HDF5 path already exist in the location list? + * If so, add the new panel to it (with a unique index, we hope) */ + for ( li=0; li<num_locations; li++ ) { + if ( strcmp(p_location, locations[li].location) == 0 ) { + add_panel_to_location(&locations[li], p, pi); + done = 1; + } + } + + /* If not, add a new location to ths list */ + if ( !done ) { + + struct hdf5_write_location *new_locations; + size_t nsz; + + nsz = (num_locations+1)*sizeof(struct hdf5_write_location); + new_locations = realloc(locations, nsz); + if ( new_locations == NULL ) { + ERROR("Failed to grow location list.\n"); + return; + } + locations = new_locations; + + locations[num_locations].max_ss = p->orig_max_ss; + locations[num_locations].max_fs = p->orig_max_fs; + locations[num_locations].location = p_location; + locations[num_locations].panel_idxs = malloc(sizeof(int)); + if ( locations[num_locations].panel_idxs == NULL ) { + ERROR("Failed to allocate single idx (!)\n"); + return; + } + locations[num_locations].panel_idxs[0] = pi; + locations[num_locations].n_panels = 1; + + num_locations += 1; + + } + + *plocations = locations; + *pnum_locations = num_locations; +} + + +static struct hdf5_write_location *make_location_list(const DataTemplate *dtempl, + int *pnum_locations) +{ + int pi; + struct hdf5_write_location *locations = NULL; + int num_locations = 0; + + for ( pi=0; pi<dtempl->n_panels; pi++ ) { + + struct panel_template *p; + const char *p_location; + + p = &dtempl->panels[pi]; + + assert(p->data != NULL); + p_location = p->data; + + add_panel_location(p, p_location, pi, + &locations, &num_locations); + + } + + *pnum_locations = num_locations; + return locations; +} + + +static void write_location(hid_t fh, const DataTemplate *dtempl, + float **dp, + struct hdf5_write_location *loc) +{ + hid_t sh, dh, ph; + hid_t dh_dataspace; + hsize_t size[2]; + int pi; + + /* Note the "swap" here, according to section 3.2.5, + * "C versus Fortran Dataspaces", of the HDF5 user's guide. */ + size[0] = loc->max_ss+1; + size[1] = loc->max_fs+1; + sh = H5Screate_simple(2, size, NULL); + + ph = H5Pcreate(H5P_LINK_CREATE); + H5Pset_create_intermediate_group(ph, 1); + + dh = H5Dcreate2(fh, loc->location, H5T_NATIVE_FLOAT, sh, + ph, H5P_DEFAULT, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't create dataset\n"); + H5Fclose(fh); + return; + } + + H5Sget_simple_extent_dims(sh, size, NULL); + + for ( pi=0; pi<loc->n_panels; pi++ ) { + + hsize_t f_offset[2], f_count[2], dims[2]; + hid_t memspace; + struct panel_template *p; + int r; + + p = &dtempl->panels[loc->panel_idxs[pi]]; + + f_offset[0] = p->orig_min_ss; + f_offset[1] = p->orig_min_fs; + f_count[0] = p->orig_max_ss - p->orig_min_ss +1; + f_count[1] = p->orig_max_fs - p->orig_min_fs +1; + + dh_dataspace = H5Dget_space(dh); + r = H5Sselect_hyperslab(dh_dataspace, H5S_SELECT_SET, + f_offset, NULL, f_count, NULL); + if ( r < 0 ) { + ERROR("Error selecting file dataspace " + "for panel %s\n", p->name); + H5Pclose(ph); + H5Dclose(dh); + H5Sclose(dh_dataspace); + H5Sclose(sh); + H5Fclose(fh); + return; + } + + dims[0] = PANEL_HEIGHT(p); + dims[1] = PANEL_WIDTH(p); + memspace = H5Screate_simple(2, dims, NULL); + + r = H5Dwrite(dh, H5T_NATIVE_FLOAT, memspace, dh_dataspace, + H5P_DEFAULT, dp[loc->panel_idxs[pi]]); + if ( r < 0 ) { + ERROR("Couldn't write data\n"); + H5Pclose(ph); + H5Dclose(dh); + H5Sclose(dh_dataspace); + H5Sclose(memspace); + H5Sclose(sh); + H5Fclose(fh); + return; + } + + H5Sclose(dh_dataspace); + H5Sclose(memspace); + } + H5Pclose(ph); + H5Sclose(sh); + H5Dclose(dh); +} + + +static void write_wavelength(hid_t fh, double wl, + const DataTemplate *dtempl) +{ + hid_t ph, sh, dh; + hsize_t size1d[1]; + int r; + + ph = H5Pcreate(H5P_LINK_CREATE); + H5Pset_create_intermediate_group(ph, 1); + + size1d[0] = 1; + sh = H5Screate_simple(1, size1d, NULL); + + dh = H5Dcreate2(fh, "/wavelength", H5T_NATIVE_DOUBLE, sh, + ph, H5S_ALL, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Couldn't create dataset for photon energy.\n"); + return; + } + r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, &wl); + if ( r < 0 ) { + ERROR("Couldn't write photon energy.\n"); + /* carry on */ + } + + H5Pclose(ph); + H5Dclose(dh); +} + + +static void write_spectrum(hid_t fh, Spectrum *s, + const DataTemplate *dtempl) +{ + herr_t r; + double *arr; + int i; + hid_t sh, dh, ph; + double kmin, kmax, step; + const hsize_t n = 1024; + + ph = H5Pcreate(H5P_LINK_CREATE); + H5Pset_create_intermediate_group(ph, 1); + + arr = malloc(n*sizeof(double)); + if ( arr == NULL ) { + ERROR("Failed to allocate memory for spectrum.\n"); + return; + } + + /* Save the wavelength values */ + spectrum_get_range(s, &kmin, &kmax); + step = (kmax-kmin)/n; + for ( i=0; i<n; i++ ) { + arr[i] = 1.0e10/(kmin+i*step); + } + + sh = H5Screate_simple(1, &n, NULL); + + dh = H5Dcreate2(fh, "/spectrum/wavelengths_A", H5T_NATIVE_DOUBLE, + sh, ph, H5S_ALL, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Failed to create dataset for spectrum wavelengths.\n"); + return; + } + r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL, + H5S_ALL, H5P_DEFAULT, arr); + if ( r < 0 ) { + ERROR("Failed to write spectrum wavelengths.\n"); + return; + } + H5Dclose(dh); + + /* Save the probability density values */ + for ( i=0; i<n; i++ ) { + arr[i] = spectrum_get_density_at_k(s, kmin+i*step); + } + + dh = H5Dcreate2(fh, "/spectrum/pdf", H5T_NATIVE_DOUBLE, sh, + H5P_DEFAULT, H5S_ALL, H5P_DEFAULT); + if ( dh < 0 ) { + ERROR("Failed to create dataset for spectrum p.d.f.\n"); + return; + } + r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL, + H5S_ALL, H5P_DEFAULT, arr); + if ( r < 0 ) { + ERROR("Failed to write spectrum p.d.f.\n"); + return; + } + + H5Dclose(dh); + H5Pclose(ph); + free(arr); +} + + +int image_hdf5_write(const struct image *image, + const DataTemplate *dtempl, + const char *filename) +{ + hid_t fh; + int li; + struct hdf5_write_location *locations; + int num_locations; + + if ( dtempl == NULL ) return 1; + + fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if ( fh < 0 ) { + ERROR("Couldn't create file: %s\n", filename); + return 1; + } + + locations = make_location_list(dtempl, &num_locations); + + for ( li=0; li<num_locations; li++ ) { + write_location(fh, dtempl, image->dp, &locations[li]); + } + + write_wavelength(fh, image->lambda, dtempl); + + if ( image->spectrum != NULL ) { + write_spectrum(fh, image->spectrum, dtempl); + } + + H5Fclose(fh); + for ( li=0; li<num_locations; li ++ ) { + free(locations[li].panel_idxs); + } + free(locations); + return 0; +} |