diff options
author | Thomas White <taw@physics.org> | 2020-05-20 22:42:49 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2020-07-29 18:42:57 +0200 |
commit | 23ea67dc03ac19f7a1457ecfdc8d5ee9cac68632 (patch) | |
tree | c9cf8cbc343b842c3da3803d67343f73071b96dd | |
parent | 302de26924528b31a2320c90fd944224674bd835 (diff) |
Port indexamajig to new API, part II
This also removes a large chunk of legacy code
-rw-r--r-- | libcrystfel/src/hdf5-file.c | 1772 | ||||
-rw-r--r-- | libcrystfel/src/hdf5-file.h | 38 | ||||
-rw-r--r-- | libcrystfel/src/image.c | 439 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/pinkindexer.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/pinkindexer.h | 3 | ||||
-rw-r--r-- | libcrystfel/src/stream.c | 14 | ||||
-rw-r--r-- | src/im-zmq.c | 23 | ||||
-rw-r--r-- | src/im-zmq.h | 2 | ||||
-rw-r--r-- | src/process_image.c | 3 |
10 files changed, 28 insertions, 2276 deletions
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c index ec0d77f2..1b156c54 100644 --- a/libcrystfel/src/hdf5-file.c +++ b/libcrystfel/src/hdf5-file.c @@ -57,34 +57,6 @@ struct hdf5_write_location { }; -int split_group_and_object(const char *path, char **group, char **object) -{ - const char *sep; - const char *store; - - sep = path; - store = sep; - sep = strpbrk(sep + 1, "/"); - if ( sep != NULL ) { - while ( 1 ) { - store = sep; - sep = strpbrk(sep + 1, "/"); - if ( sep == NULL ) { - break; - } - } - } - if ( store == path ) { - *group = NULL; - *object = strdup(path); - } else { - *group = strndup(path, store - path); - *object = strdup(store+1); - } - return 0; -}; - - struct hdfile { const char *path; /* Current data path */ @@ -134,564 +106,6 @@ int hdfile_set_image(struct hdfile *f, const char *path) } -static int read_peak_count(struct hdfile *f, char *path, int line, - int *num_peaks) -{ - - hid_t dh, sh, mh; - hsize_t size[1]; - hsize_t max_size[1]; - hsize_t offset[1], count[1]; - hsize_t m_offset[1], m_count[1], dimmh[1]; - - - int tw, r; - - dh = H5Dopen2(f->fh, path, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Data block %s not found.\n", path); - return 1; - } - - sh = H5Dget_space(dh); - if ( sh < 0 ) { - H5Dclose(dh); - ERROR("Couldn't get dataspace for data.\n"); - return 1; - } - - if ( H5Sget_simple_extent_ndims(sh) != 1 ) { - ERROR("Data block %s has the wrong dimensionality (%i).\n", - path, H5Sget_simple_extent_ndims(sh)); - H5Sclose(sh); - H5Dclose(dh); - return 1; - } - - H5Sget_simple_extent_dims(sh, size, max_size); - - tw = size[0]; - - if ( line > tw-1 ) { - H5Sclose(sh); - H5Dclose(dh); - ERROR("Data block %s does not contain data for required event.\n", - path); - return 1; - } - - offset[0] = line; - count[0] = 1; - - r = H5Sselect_hyperslab(sh, H5S_SELECT_SET, - offset, NULL, count, NULL); - if ( r < 0 ) { - ERROR("Error selecting file dataspace " - "for data block %s\n", path); - H5Dclose(dh); - H5Sclose(sh); - return 1; - } - - m_offset[0] = 0; - m_count[0] = 1; - dimmh[0] = 1; - mh = H5Screate_simple(1, dimmh, NULL); - r = H5Sselect_hyperslab(mh, H5S_SELECT_SET, - m_offset, NULL, m_count, NULL); - if ( r < 0 ) { - ERROR("Error selecting memory dataspace " - "for data block %s\n", path); - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return 1; - } - - r = H5Dread(dh, H5T_NATIVE_INT, mh, - sh, H5P_DEFAULT, num_peaks); - if ( r < 0 ) { - ERROR("Couldn't read data for block %s, line %i\n", path, line); - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return 1; - } - - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return 0; -} - - - -static float *read_hdf5_data(struct hdfile *f, char *path, int line) -{ - - hid_t dh, sh, mh; - hsize_t size[2]; - hsize_t max_size[2]; - hsize_t offset[2], count[2]; - hsize_t m_offset[2], m_count[2], dimmh[2]; - float *buf; - int tw, r; - - dh = H5Dopen2(f->fh, path, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Data block (%s) not found.\n", path); - return NULL; - } - - sh = H5Dget_space(dh); - if ( sh < 0 ) { - H5Dclose(dh); - ERROR("Couldn't get dataspace for data.\n"); - return NULL; - } - - if ( H5Sget_simple_extent_ndims(sh) != 2 ) { - ERROR("Data block %s has the wrong dimensionality (%i).\n", - path, H5Sget_simple_extent_ndims(sh)); - H5Sclose(sh); - H5Dclose(dh); - return NULL; - } - - H5Sget_simple_extent_dims(sh, size, max_size); - - tw = size[0]; - if ( line> tw-1 ) { - H5Sclose(sh); - H5Dclose(dh); - ERROR("Data block %s does not contain data for required event.\n", - path); - return NULL; - } - - offset[0] = line; - offset[1] = 0; - count[0] = 1; - count[1] = size[1]; - - r = H5Sselect_hyperslab(sh, H5S_SELECT_SET, offset, NULL, count, NULL); - if ( r < 0 ) { - ERROR("Error selecting file dataspace " - "for data block %s\n", path); - H5Dclose(dh); - H5Sclose(sh); - return NULL; - } - - m_offset[0] = 0; - m_offset[1] = 0; - m_count[0] = 1; - m_count[1] = size[1]; - dimmh[0] = 1; - dimmh[1] = size[1]; - - mh = H5Screate_simple(2, dimmh, NULL); - r = H5Sselect_hyperslab(mh, H5S_SELECT_SET, - m_offset, NULL, m_count, NULL); - if ( r < 0 ) { - ERROR("Error selecting memory dataspace " - "for data block %s\n", path); - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return NULL; - } - - buf = malloc(size[1]*sizeof(float)); - if ( buf == NULL ) return NULL; - r = H5Dread(dh, H5T_NATIVE_FLOAT, mh, sh, H5P_DEFAULT, buf); - if ( r < 0 ) { - ERROR("Couldn't read data for block %s, line %i\n", path, line); - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return NULL; - } - - H5Dclose(dh); - H5Sclose(sh); - H5Sclose(mh); - return buf; -} - - -/** - * \param image: An \ref image structure - * \param f: An \ref hdfile structure - * \param p: The HDF5 path to the peak data - * \param fpe: A \ref filename_plus_event structure specifying the event - * \param half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates - * - * Get peaks from HDF5, in "CXI format" (as in "CXIDB"). The data should be in - * a set of arrays under \p p. The number of peaks should be in a 1D array at - * \p p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at - * \p p/peakXPosRaw and \p p/peakYPosRaw respectively (sorry about the naming). The - * first dimension of these arrays should be the event number (as given by - * \p fpe). The intensities are expected to be at \p p/peakTotalIntensity in a - * similar 2D array. - * - * CrystFEL considers all peak locations to be distances from the corner of the - * detector panel, in pixel units, consistent with its description of detector - * geometry (see 'man crystfel_geometry'). The software which generates the - * CXI files, including Cheetah, may instead consider the peak locations to be - * pixel indices in the data array. In this case, the peak coordinates should - * have 0.5 added to them. This will be done if \p half_pixel_shift is non-zero. - * - * \returns Non-zero on error, zero otherwise. - * - */ -int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe, int half_pixel_shift) -{ - char path_n[1024]; - char path_x[1024]; - char path_y[1024]; - char path_i[1024]; - int r; - int pk; - - int line = 0; - int num_peaks; - - float *buf_x; - float *buf_y; - float *buf_i; - - double peak_offset = half_pixel_shift ? 0.5 : 0.0; - - if ( (fpe != NULL) && (fpe->ev != NULL) - && (fpe->ev->dim_entries != NULL) ) - { - line = fpe->ev->dim_entries[0]; - } else { - ERROR("CXI format peak list format selected," - "but file has no event structure"); - return 1; - } - - snprintf(path_n, 1024, "%s/nPeaks", p); - snprintf(path_x, 1024, "%s/peakXPosRaw", p); - snprintf(path_y, 1024, "%s/peakYPosRaw", p); - snprintf(path_i, 1024, "%s/peakTotalIntensity", p); - - r = read_peak_count(f, path_n, line, &num_peaks); - if ( r != 0 ) return 1; - - buf_x = read_hdf5_data(f, path_x, line); - if ( r != 0 ) return 1; - - buf_y = read_hdf5_data(f, path_y, line); - if ( r != 0 ) return 1; - - buf_i = read_hdf5_data(f, path_i, line); - if ( r != 0 ) return 1; - - if ( image->features != NULL ) { - image_feature_list_free(image->features); - } - image->features = image_feature_list_new(); - - for ( pk=0; pk<num_peaks; pk++ ) { - - float fs, ss, val; - struct panel *p; - int pn; - - fs = buf_x[pk] + peak_offset; - ss = buf_y[pk] + peak_offset; - val = buf_i[pk]; - - p = find_orig_panel(image->det, fs, ss); - if ( p == NULL ) continue; - if ( p->no_index ) continue; - - /* Convert coordinates to panel-relative */ - fs = fs - p->orig_min_fs; - ss = ss - p->orig_min_ss; - - pn = panel_number(image->det, p); - - image_add_feature(image->features, fs, ss, pn, - image, val, NULL); - - } - - return 0; -} - - -/** - * \param image: An \ref image structure - * \param f: An \ref hdfile structure - * \param p: The HDF5 path to the peak data - * \param fpe: A \ref filename_plus_event structure specifying the event - * - * This is a wrapper function to preserve API compatibility with older CrystFEL - * versions. Use \ref get_peaks_cxi_2 instead. - * - * This function is equivalent to get_peaks_cxi_2(\p image, \p f, \p p, \p fpe, 1). - * - * \returns Non-zero on error, zero otherwise. - * - */ -int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe) -{ - return get_peaks_cxi_2(image, f, p, fpe, 1); -} - - -/** - * \param image: An \ref image structure - * \param f: An \ref hdfile structure - * \param p: The HDF5 path to the peak data - * \param half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates - * - * Get peaks from HDF5. The peak list should be located at \p p in the HDF5 file, - * a 2D array where the first dimension equals the number of peaks and second - * dimension is three. The first two columns contain the fast scan and slow - * scan coordinates, respectively, of the peaks. The third column contains the - * intensities. - * - * CrystFEL considers all peak locations to be distances from the corner of the - * detector panel, in pixel units, consistent with its description of detector - * geometry (see 'man crystfel_geometry'). The software which generates the - * CXI files, including Cheetah, may instead consider the peak locations to be - * pixel indices in the data array. In this case, the peak coordinates should - * have 0.5 added to them. This will be done if \p half_pixel_shift is non-zero. - * - * \returns Non-zero on error, zero otherwise. - * - */ -int get_peaks_2(struct image *image, struct hdfile *f, const char *p, - int half_pixel_shift) -{ - hid_t dh, sh; - hsize_t size[2]; - hsize_t max_size[2]; - int i; - float *buf; - herr_t r; - int tw; - char *np; - double peak_offset = half_pixel_shift ? 0.5 : 0.0; - - if ( image->event != NULL ) { - np = retrieve_full_path(image->event, p); - } else { - np = strdup(p); - } - - dh = H5Dopen2(f->fh, np, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Peak list (%s) not found.\n", np); - return 1; - } - - sh = H5Dget_space(dh); - if ( sh < 0 ) { - H5Dclose(dh); - ERROR("Couldn't get dataspace for peak list.\n"); - free(np); - 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); - free(np); - 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"); - free(np); - 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"); - free(np); - 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); - free(np); - return 1; - } - - if ( image->features != NULL ) { - image_feature_list_free(image->features); - } - image->features = image_feature_list_new(); - - for ( i=0; i<size[0]; i++ ) { - - float fs, ss, val; - struct panel *p; - int pn; - - fs = buf[tw*i+0] + peak_offset; - ss = buf[tw*i+1] + peak_offset; - val = buf[tw*i+2]; - - p = find_orig_panel(image->det, fs, ss); - if ( p == NULL ) continue; - if ( p->no_index ) continue; - - /* Convert coordinates to panel-relative */ - fs = fs - p->orig_min_fs; - ss = ss - p->orig_min_ss; - - pn = panel_number(image->det, p); - - image_add_feature(image->features, fs, ss, pn, - image, val, NULL); - - } - - free(buf); - free(np); - H5Sclose(sh); - H5Dclose(dh); - - return 0; -} - - -/** - * \param image: An \ref image structure - * \param f: An \ref hdfile structure - * \param p: The HDF5 path to the peak data - * - * This is a wrapper function to preserve API compatibility with older CrystFEL - * versions. Use \ref get_peaks_2 instead. - * - * This function is equivalent to \ref get_peaks_2(\p image, \p f, \p p, 1). - * - * \returns Non-zero on error, zero otherwise. - * - */ -int get_peaks(struct image *image, struct hdfile *f, const char *p) -{ - return get_peaks_2(image, f, p, 1); -} - - -static void cleanup(hid_t fh) -{ - int n_ids, i; - hid_t ids[2048]; - - n_ids = H5Fget_obj_ids(fh, H5F_OBJ_ALL, 2048, ids); - - for ( i=0; i<n_ids; i++ ) { - - hid_t id; - H5I_type_t type; - - id = ids[i]; - - type = H5Iget_type(id); - - if ( type == H5I_GROUP ) H5Gclose(id); - if ( type == H5I_DATASET ) H5Dclose(id); - if ( type == H5I_DATATYPE ) H5Tclose(id); - if ( type == H5I_DATASPACE ) H5Sclose(id); - if ( type == H5I_ATTR ) H5Aclose(id); - - } - -} - - -void hdfile_close(struct hdfile *f) -{ - - if ( f->data_open ) { - H5Dclose(f->dh); - } - - cleanup(f->fh); - - H5Fclose(f->fh); - - free(f); -} - - -/* Deprecated */ -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 */ - herr_t r; - hsize_t 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; - sh = H5Screate_simple(2, size, NULL); - - dh = H5Dcreate2(gh, "data", type, sh, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Couldn't create dataset\n"); - H5Fclose(fh); - return 1; - } - - /* Muppet check */ - H5Sget_simple_extent_dims(sh, size, NULL); - - 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; - } - H5Dclose(dh); - H5Gclose(gh); - H5Fclose(fh); - - return 0; -} - - static void add_panel_to_location(struct hdf5_write_location *loc, struct panel *p, int pi) { @@ -1030,149 +444,6 @@ int hdf5_write_image(const char *filename, const struct image *image, } -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; i<size[0]; i++ ) { - - unsigned int fs, ss; - float val; - struct panel *p; - signed int pn; - - fs = buf[3*i+0]; - ss = buf[3*i+1]; - val = buf[3*i+2]; - - /* Turn "original" position into "panel" position */ - pn = find_orig_panel_number(image->det, fs, ss); - if ( pn == -1 ) { - ERROR("Failed to find panel!\n"); - continue; - } - p = &image->det->panels[pn]; - - image->dp[pn][fs+p->w*ss] = val/5.0; - image->dp[pn][fs+1+p->w*ss] = val/5.0; - image->dp[pn][fs-1+p->w*ss] = val/5.0; - image->dp[pn][fs+p->w*(ss-1)] = val/5.0; - image->dp[pn][fs+p->w*(ss+1)] = val/5.0; - - } - - free(buf); - H5Sclose(sh); - H5Dclose(dh); -} - - -static int *make_badmask(int *flags, struct detector *det, float *data, - struct panel *p) -{ - int *badmap; - int fs, ss; - - badmap = malloc(p->w*p->h*sizeof(int)); - if ( badmap == NULL ) { - ERROR("Failed to allocate bad mask for panel %s\n", - p->name); - return NULL; - } - - /* Defaults, then bad pixels arising from bad regions or panels */ - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - int bad = 0; - - if ( p->no_index ) bad = 1; - - if ( in_bad_region(det, p, fs, ss) ) { - bad = 1; - } - - badmap[fs+p->w*ss] = bad; - } - } - - /* Bad pixels from mask */ - if ( flags != NULL ) { - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - int f = flags[fs+p->w*ss]; - int bad = badmap[fs+p->w*ss]; - float val = data[fs+p->w*ss]; - - /* Bad if it's missing any of the "good" bits */ - if ( (f & det->mask_good) != det->mask_good ) bad = 1; - - /* Bad if it has any of the "bad" bits. */ - if ( f & det->mask_bad ) bad = 1; - - /* Bad if pixel value is NaN or inf */ - if ( isnan(val) || isinf(val) ) bad = 1; - - badmap[fs+p->w*ss] = bad; - - } - } - } - - return badmap; -} - - int hdfile_get_value(struct hdfile *f, const char *name, struct event *ev, void *val, hid_t memtype) { @@ -1315,638 +586,6 @@ int hdfile_get_value(struct hdfile *f, const char *name, struct event *ev, } -static void hdfile_fill_in_beam_parameters(struct beam_params *beam, - struct hdfile *f, - struct event *ev, - struct image *image) -{ - double eV; - - if ( beam->photon_energy_from == NULL ) { - - /* Explicit value given */ - eV = beam->photon_energy; - - } else { - - int r; - - r = hdfile_get_value(f, beam->photon_energy_from, - ev, &eV, H5T_NATIVE_DOUBLE); - if ( r ) { - ERROR("Failed to read '%s'\n", - beam->photon_energy_from); - } - - } - - image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale; -} - - -static void hdfile_fill_in_clen(struct detector *det, struct hdfile *f, - struct event *ev) -{ - int i; - - for ( i=0; i<det->n_panels; i++ ) { - - struct panel *p = &det->panels[i]; - - if ( p->clen_from != NULL ) { - - double val; - int r; - - r = hdfile_get_value(f, p->clen_from, ev, &val, - H5T_NATIVE_DOUBLE); - if ( r ) { - ERROR("Failed to read '%s'\n", p->clen_from); - } else { - p->clen = val * 1.0e-3; - } - - } - - adjust_centering_for_rail(p); - - } -} - - -int hdf5_read(struct hdfile *f, struct image *image, const char *element, - int satcorr) -{ - herr_t r; - float *buf; - int fail; - hsize_t *size; - hsize_t *max_size; - hid_t sh; - int sh_dim; - int w, h; - - if ( element == NULL ) { - fail = hdfile_set_first_image(f, "/"); - } else { - fail = hdfile_set_image(f, element); - } - - if ( fail ) { - ERROR("Couldn't select path\n"); - return 1; - } - - sh = H5Dget_space(f->dh); - sh_dim = H5Sget_simple_extent_ndims(sh); - - if ( sh_dim != 2 ) { - ERROR("Dataset is not two-dimensional\n"); - return -1; - } - - size = malloc(sh_dim*sizeof(hsize_t)); - max_size = malloc(sh_dim*sizeof(hsize_t)); - H5Sget_simple_extent_dims(sh, size, max_size); - H5Sclose(sh); - w = size[1]; - h = size[0]; - free(size); - free(max_size); - - buf = malloc(sizeof(float)*w*h); - 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; - } - - if ( image->det != NULL ) { - ERROR("WARNING: hdf5_read() called with geometry structure.\n"); - } - image->det = simple_geometry(image, w, h); - image->dp = malloc(sizeof(double *)); - if ( image->dp == NULL ) { - ERROR("Failed to allocate memory for image data!\n"); - return 1; - } - image->dp[0] = buf; - - if ( satcorr ) debodge_saturation(f, image); - - if ( image->beam != NULL ) { - - hdfile_fill_in_beam_parameters(image->beam, f, NULL, image); - - if ( image->lambda > 1000 ) { - /* Error message covers a silly value in the beam file - * or in the HDF5 file. */ - ERROR("WARNING: Missing or nonsensical wavelength " - "(%e m) for %s.\n", - image->lambda, image->filename); - } - - } - - fill_in_adu(image); - - return 0; -} - - -static hsize_t *first_two_dims(hsize_t *in, struct dim_structure *ds) -{ - int i, j; - hsize_t *out = malloc(2*sizeof(hsize_t)); - - if ( out == NULL ) return NULL; - - j = 0; - for ( i=0; i<ds->num_dims; i++ ) { - if ( (ds->dims[i] == HYSL_FS) || (ds->dims[i] == HYSL_SS) ) { - out[j++] = in[i]; - } - } - return out; -} - - -static int load_satmap(struct hdfile *f, struct event *ev, struct panel *p, - hsize_t *in_f_offset, hsize_t *in_f_count, - struct dim_structure *dim_struct, - float *satmap) -{ - char *loc; /* Sat map location after possible substitution */ - hid_t satmap_dataspace, satmap_dh; - int exists; - int check, r; - hid_t memspace; - hsize_t dimsm[2]; - hid_t fh; - hsize_t *f_offset, *f_count; - - if ( p->satmap_file != NULL ) { - - fh = H5Fopen(p->satmap_file, H5F_ACC_RDONLY, H5P_DEFAULT); - if ( fh < 0 ) { - ERROR("Couldn't open satmap file '%s'\n", p->satmap_file); - return 1; - } - - /* If we have an external map file, we assume it to be a simple - * 2D job */ - f_offset = first_two_dims(in_f_offset, dim_struct); - f_count = first_two_dims(in_f_count, dim_struct); - - } else { - - /* Otherwise, we assume it has the same dimensions as the - * image data itself */ - fh = f->fh; - f_offset = in_f_offset; - f_count = in_f_count; - } - - if ( ev != NULL ) { - loc = retrieve_full_path(ev, p->satmap); - } else { - loc = strdup(p->satmap); - } - - exists = check_path_existence(fh, loc); - if ( !exists ) { - ERROR("Cannot find satmap for panel %s\n", p->name); - goto err; - } - - satmap_dh = H5Dopen2(fh, loc, H5P_DEFAULT); - if ( satmap_dh <= 0 ) { - ERROR("Couldn't open satmap for panel %s\n", p->name); - goto err; - } - - satmap_dataspace = H5Dget_space(satmap_dh); - check = H5Sselect_hyperslab(satmap_dataspace, H5S_SELECT_SET, - f_offset, NULL, f_count, NULL); - if ( check < 0 ) { - ERROR("Error selecting satmap dataspace for panel %s\n", - p->name); - goto err; - } - - dimsm[0] = p->h; - dimsm[1] = p->w; - memspace = H5Screate_simple(2, dimsm, NULL); - if ( check < 0 ) { - ERROR("Error selecting satmap memory dataspace for panel %s\n", - p->name); - goto err; - } - - r = H5Dread(satmap_dh, H5T_NATIVE_FLOAT, memspace, - satmap_dataspace, H5P_DEFAULT, satmap); - if ( r < 0 ) { - ERROR("Couldn't read satmap for panel %s\n", p->name); - goto err; - } - - H5Sclose(satmap_dataspace); - H5Dclose(satmap_dh); - free(loc); - - return 0; - -err: - if ( p->satmap_file != NULL ) H5Fclose(fh); - free(loc); - return 1; -} - - - -static int load_mask(struct hdfile *f, struct event *ev, struct panel *p, - int *flags, - hsize_t *in_f_offset, hsize_t *in_f_count, - struct dim_structure *dim_struct) -{ - char *mask; /* Mask location after possible substitution */ - hid_t mask_dataspace, mask_dh; - int exists; - int check, r; - hid_t memspace; - hsize_t dimsm[2]; - hid_t fh; - hsize_t *f_offset, *f_count; - - if ( p->mask_file != NULL ) { - - fh = H5Fopen(p->mask_file, H5F_ACC_RDONLY, H5P_DEFAULT); - if ( fh < 0 ) { - ERROR("Couldn't open mask file '%s'\n", p->mask_file); - return 1; - } - - /* If we have an external map file, we assume it to be a simple - * 2D job */ - f_offset = first_two_dims(in_f_offset, dim_struct); - f_count = first_two_dims(in_f_count, dim_struct); - - } else { - fh = f->fh; - f_offset = in_f_offset; - f_count = in_f_count; - } - - if ( ev != NULL ) { - mask = retrieve_full_path(ev, p->mask); - } else { - mask = strdup(p->mask); - } - - exists = check_path_existence(fh, mask); - if ( !exists ) { - ERROR("Cannot find flags for panel %s\n", p->name); - goto err; - } - - mask_dh = H5Dopen2(fh, mask, H5P_DEFAULT); - if ( mask_dh <= 0 ) { - ERROR("Couldn't open flags for panel %s\n", p->name); - goto err; - } - - mask_dataspace = H5Dget_space(mask_dh); - check = H5Sselect_hyperslab(mask_dataspace, H5S_SELECT_SET, - f_offset, NULL, f_count, NULL); - if ( check < 0 ) { - ERROR("Error selecting mask dataspace for panel %s\n", p->name); - goto err; - } - - dimsm[0] = p->h; - dimsm[1] = p->w; - memspace = H5Screate_simple(2, dimsm, NULL); - if ( check < 0 ) { - ERROR("Error selecting memory dataspace for panel %s\n", p->name); - goto err; - } - - r = H5Dread(mask_dh, H5T_NATIVE_INT, memspace, - mask_dataspace, H5P_DEFAULT, flags); - if ( r < 0 ) { - ERROR("Couldn't read flags for panel %s\n", p->name); - goto err; - } - - H5Sclose(mask_dataspace); - H5Dclose(mask_dh); - free(mask); - - return 0; - -err: - if ( p->mask_file != NULL ) H5Fclose(fh); - free(mask); - return 1; -} - - -int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev, - int satcorr) -{ - herr_t r; - int pi; - int i; - - if ( image->det == NULL ) { - ERROR("Geometry not available\n"); - return 1; - } - - 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 data arrays.\n"); - return 1; - } - - for ( pi=0; pi<image->det->n_panels; pi++ ) { - - hsize_t *f_offset, *f_count; - int hsi; - struct dim_structure *hsd; - herr_t check; - hid_t dataspace, memspace; - int fail; - struct panel *p; - hsize_t dims[2]; - - p = &image->det->panels[pi]; - - if ( ev != NULL ) { - - int exists; - char *panel_full_path; - - panel_full_path = retrieve_full_path(ev, p->data); - - exists = check_path_existence(f->fh, panel_full_path); - if ( !exists ) { - ERROR("Cannot find data for panel %s\n", - p->name); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - - fail = hdfile_set_image(f, panel_full_path); - - free(panel_full_path); - - } else { - - if ( p->data == NULL ) { - - fail = hdfile_set_first_image(f, "/"); - - } else { - - int exists; - exists = check_path_existence(f->fh, p->data); - if ( !exists ) { - ERROR("Cannot find data for panel %s\n", - p->name); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - fail = hdfile_set_image(f, p->data); - - } - - } - if ( fail ) { - ERROR("Couldn't select path for panel %s\n", - p->name); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - - /* Determine where to read the data from in the file */ - hsd = image->det->panels[pi].dim_structure; - f_offset = malloc(hsd->num_dims*sizeof(hsize_t)); - f_count = malloc(hsd->num_dims*sizeof(hsize_t)); - if ( (f_offset == NULL) || (f_count == NULL ) ) { - ERROR("Failed to allocate offset or count.\n"); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - for ( hsi=0; hsi<hsd->num_dims; hsi++ ) { - - if ( hsd->dims[hsi] == HYSL_FS ) { - f_offset[hsi] = p->orig_min_fs; - f_count[hsi] = p->orig_max_fs - p->orig_min_fs+1; - } else if ( hsd->dims[hsi] == HYSL_SS ) { - f_offset[hsi] = p->orig_min_ss; - f_count[hsi] = p->orig_max_ss - p->orig_min_ss+1; - } else if (hsd->dims[hsi] == HYSL_PLACEHOLDER ) { - f_offset[hsi] = ev->dim_entries[0]; - f_count[hsi] = 1; - } else { - f_offset[hsi] = hsd->dims[hsi]; - f_count[hsi] = 1; - } - - } - - /* Set up dataspace for file */ - dataspace = H5Dget_space(f->dh); - check = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, - f_offset, NULL, f_count, NULL); - if ( check < 0 ) { - ERROR("Error selecting file dataspace for panel %s\n", - p->name); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - - dims[0] = p->h; - dims[1] = p->w; - memspace = H5Screate_simple(2, dims, NULL); - - image->dp[pi] = malloc(p->w*p->h*sizeof(float)); - image->sat[pi] = malloc(p->w*p->h*sizeof(float)); - if ( (image->dp[pi] == NULL) || (image->sat[pi] == NULL) ) { - ERROR("Failed to allocate panel %s\n", p->name); - free(f_offset); - free(f_count); - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - for ( i=0; i<p->w*p->h; i++ ) image->sat[pi][i] = INFINITY; - - r = H5Dread(f->dh, H5T_NATIVE_FLOAT, memspace, dataspace, - H5P_DEFAULT, image->dp[pi]); - if ( r < 0 ) { - ERROR("Couldn't read data for panel %s\n", - p->name); - free(f_offset); - free(f_count); - for ( i=0; i<=pi; i++ ) { - free(image->dp[i]); - free(image->sat[i]); - } - free(image->dp); - free(image->bad); - free(image->sat); - return 1; - } - - if ( p->mask != NULL ) { - int *flags = malloc(p->w*p->h*sizeof(int)); - if ( !load_mask(f, ev, p, flags, f_offset, f_count, hsd) ) { - image->bad[pi] = make_badmask(flags, image->det, - image->dp[pi], p); - } else { - image->bad[pi] = make_badmask(NULL, image->det, - image->dp[pi], p); - } - free(flags); - } else { - image->bad[pi] = make_badmask(NULL, image->det, - image->dp[pi], p); - } - - if ( p->satmap != NULL ) { - if ( load_satmap(f, ev, p, f_offset, f_count, hsd, - image->sat[pi]) ) - { - ERROR("Failed to load sat map for panel %s\n", - p->name); - } - } - - H5Sclose(dataspace); - free(f_offset); - free(f_count); - - } - - H5Dclose(f->dh); - f->data_open = 0; - hdfile_fill_in_clen(image->det, f, ev); - - if ( satcorr ) debodge_saturation(f, image); - - if ( image->beam != NULL ) { - - hdfile_fill_in_beam_parameters(image->beam, f, ev, image); - - if ( (image->lambda > 1.0) || (image->lambda < 1e-20) ) { - - ERROR("WARNING: Nonsensical wavelength (%e m) value " - "for file: %s, event: %s.\n", - image->lambda, image->filename, - get_event_string(image->event)); - } - - } - - fill_in_adu(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 ) { - H5Sclose(sh); - return 0; - } - - H5Sget_simple_extent_dims(sh, size, max_size); - H5Sclose(sh); - - if ( ( size[0] > 64 ) && ( size[1] > 64 ) ) return 1; - - return 0; -} - - -int hdfile_is_scalar(struct hdfile *f, const char *name, int verbose) -{ - hid_t dh; - hid_t sh; - hsize_t size[3]; - hid_t type; - int ndims; - int i; - int check; - - check = check_path_existence(f->fh, name); - if ( check == 0 ) { - ERROR("No such scalar field '%s'\n", name); - return 0; - } - - dh = H5Dopen2(f->fh, name, H5P_DEFAULT); - type = H5Dget_type(dh); - - /* Get the dimensionality. We have to cope with scalars expressed as - * arrays with all dimensions 1, as well as zero-d arrays. */ - sh = H5Dget_space(dh); - ndims = H5Sget_simple_extent_ndims(sh); - if ( ndims > 3 ) { - if ( verbose ) { - ERROR("Too many dimensions (%i).\n", ndims); - } - H5Sclose(sh); - H5Tclose(type); - H5Dclose(dh); - return 0; - } - - /* Check that the size in all dimensions is 1 */ - H5Sget_simple_extent_dims(sh, size, NULL); - H5Sclose(sh); - H5Tclose(type); - H5Dclose(dh); - for ( i=0; i<ndims; i++ ) { - if ( size[i] != 1 ) { - if ( verbose ) { - ERROR("%s not a scalar value (ndims=%i," - "size[%i]=%i)\n", - name, ndims, i, (int)size[i]); - } - return 0; - } - } - - return 1; -} - - struct copy_hdf5_field { char **fields; @@ -2302,124 +941,6 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name, } -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<num; i++ ) { - - char buf[256]; - hid_t dh; - H5I_type_t type; - - H5Lget_name_by_idx(gh, ".", H5_INDEX_NAME, H5_ITER_NATIVE, - i, buf, 255, H5P_DEFAULT); - res[i] = malloc(512); - if ( strlen(parent) > 1 ) { - snprintf(res[i], 511, "%s/%s", parent, buf); - } else { - snprintf(res[i], 511, "%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); - - } - - H5Gclose(gh); - - 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<n; i++ ) { - - if ( is_image[i] ) { - hdfile_set_image(f, names[i]); - for ( j=0; j<n; j++ ) free(names[j]); - free(is_image); - free(is_group); - free(names); - return 0; - } else if ( is_group[i] ) { - if ( !hdfile_set_first_image(f, names[i]) ) { - for ( j=0; j<n; j++ ) free(names[j]); - free(is_image); - free(is_group); - free(names); - return 0; - } - } - - } - - for ( j=0; j<n; j++ ) free(names[j]); - free(is_image); - free(is_group); - free(names); - - return 1; -} - - -struct parse_params { - struct hdfile *hdfile; - int path_dim; - const char *path; - struct event *curr_event; - struct event_list *ev_list; - int top_level; -}; - - int check_path_existence(hid_t fh, const char *path) { @@ -2497,296 +1018,3 @@ int check_path_existence(hid_t fh, const char *path) return 1; } - - -static herr_t parse_file_event_structure(hid_t loc_id, char *name, - const H5L_info_t *info, - struct parse_params *pp) - -{ - char *substituted_path; - char *ph_loc; - char *truncated_path; - htri_t check; - herr_t herrt_iterate, herrt_info; - struct H5O_info_t object_info; - - if ( !pp->top_level ) { - - int fail_push; - - fail_push = push_path_entry_to_event(pp->curr_event, name); - if ( fail_push ) { - return -1; - } - - substituted_path = event_path_placeholder_subst(name, pp->path); - - } else { - substituted_path = strdup(pp->path); - } - - if ( pp->top_level == 1 ) { - pp->top_level = 0; - } - - truncated_path = strdup(substituted_path); - ph_loc = strstr(substituted_path,"%"); - if ( ph_loc != NULL) { - truncated_path[ph_loc-substituted_path] = '\0'; - } - - herrt_iterate = 0; - herrt_info = 0; - - check = check_path_existence(pp->hdfile->fh, truncated_path); - if ( check == 0 ) { - pop_path_entry_from_event(pp->curr_event); - return 0; - } else { - - herrt_info = H5Oget_info_by_name(pp->hdfile->fh, truncated_path, - &object_info, H5P_DEFAULT); - if ( herrt_info < 0 ) { - free(truncated_path); - free(substituted_path); - return -1; - } - - if ( pp->curr_event->path_length == pp->path_dim - && object_info.type == H5O_TYPE_DATASET ) - { - - int fail_append; - - fail_append = append_event_to_event_list(pp->ev_list, - pp->curr_event); - if ( fail_append ) { - free(truncated_path); - free(substituted_path); - return -1; - } - - pop_path_entry_from_event(pp->curr_event); - return 0; - - } else { - - pp->path = substituted_path; - - if ( object_info.type == H5O_TYPE_GROUP ) { - - herrt_iterate = H5Literate_by_name(pp->hdfile->fh, - truncated_path, H5_INDEX_NAME, - H5_ITER_NATIVE, NULL, - (H5L_iterate_t)parse_file_event_structure, - (void *)pp, H5P_DEFAULT); - } - } - } - - pop_path_entry_from_event(pp->curr_event); - - free(truncated_path); - free(substituted_path); - - return herrt_iterate; -} - - -static int fill_paths(struct hdfile *hdfile, struct detector *det, int pi, - struct event_list *master_el) -{ - struct parse_params pparams; - struct event *empty_event; - struct event_list *panel_ev_list; - int ei; - int check; - - empty_event = initialize_event(); - panel_ev_list = initialize_event_list(); - if ( (empty_event == NULL) || (panel_ev_list == NULL) ) - { - ERROR("Failed to allocate memory for event list.\n"); - return 1; - } - - pparams.path = det->panels[pi].data; - pparams.hdfile = hdfile; - pparams.path_dim = det->path_dim; - pparams.curr_event = empty_event; - pparams.top_level = 1; - pparams.ev_list = panel_ev_list; - - check = parse_file_event_structure(hdfile->fh, NULL, NULL, &pparams); - if ( check < 0 ) { - free_event(empty_event); - free_event_list(panel_ev_list); - return 1; - } - - for ( ei=0; ei<panel_ev_list->num_events; ei++ ) { - - int fail_add; - - fail_add = add_non_existing_event_to_event_list(master_el, - panel_ev_list->events[ei]); - if ( fail_add ) { - free_event(empty_event); - free_event_list(panel_ev_list); - return 1; - } - - } - - free_event(empty_event); - free_event_list(panel_ev_list); - - return 0; -} - - -static int check_dims(struct hdfile *hdfile, struct panel *p, struct event *ev, - struct event_list *events, int *global_path_dim) -{ - char *full_panel_path; - hid_t dh; - hid_t sh; - int dims; - hsize_t *size; - hsize_t *max_size; - int hsdi; - int panel_path_dim = 0; - struct dim_structure *panel_dim_structure; - - /* Get the full path for this panel in this event */ - full_panel_path = retrieve_full_path(ev, p->data); - - dh = H5Dopen2(hdfile->fh, full_panel_path, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Error opening '%s'\n", full_panel_path); - ERROR("Failed to enumerate events. " - "Check your geometry file.\n"); - return 1; - } - - sh = H5Dget_space(dh); - dims = H5Sget_simple_extent_ndims(sh); - size = malloc(dims*sizeof(hsize_t)); - max_size = malloc(dims*sizeof(hsize_t)); - if ( (size==NULL) || (max_size==NULL) ) { - ERROR("Failed to allocate memory for dimensions\n"); - return 1; - } - - dims = H5Sget_simple_extent_dims(sh, size, max_size); - - panel_dim_structure = p->dim_structure; - for ( hsdi=0; hsdi<panel_dim_structure->num_dims; hsdi++ ) { - if ( panel_dim_structure->dims[hsdi] == HYSL_PLACEHOLDER ) { - panel_path_dim = size[hsdi]; - break; - } - } - - if ( *global_path_dim == -1 ) { - - *global_path_dim = panel_path_dim; - - } else if ( panel_path_dim != *global_path_dim ) { - - ERROR("All panels must have the same number of frames\n"); - ERROR("Panel %s has %i frames in one dimension, but the first " - "panel has %i.\n", - p->name, panel_path_dim, *global_path_dim); - free(size); - free(max_size); - return 1; - } - - H5Sclose(sh); - H5Dclose(dh); - - return 0; -} - - -struct event_list *fill_event_list(struct hdfile *hdfile, struct detector *det) -{ - struct event_list *master_el; - - master_el = initialize_event_list(); - if ( master_el == NULL ) { - ERROR("Failed to allocate event list.\n"); - return NULL; - } - - /* First expand any placeholders in the HDF5 paths */ - if ( det->path_dim != 0 ) { - int pi; - for ( pi=0; pi<det->n_panels; pi++ ) { - if ( fill_paths(hdfile, det, pi, master_el) ) { - ERROR("Failed to enumerate paths.\n"); - return NULL; - } - } - } - - /* Now enumerate the placeholder dimensions */ - if ( det->dim_dim > 0 ) { - - struct event_list *master_el_with_dims; - int evi; - - /* If there were no HDF5 path placeholders, add a dummy event */ - if ( master_el->num_events == 0 ) { - struct event *empty_ev; - empty_ev = initialize_event(); - append_event_to_event_list(master_el, empty_ev); - free(empty_ev); - } - - master_el_with_dims = initialize_event_list(); - - /* For each event so far, expand the dimensions */ - for ( evi=0; evi<master_el->num_events; evi++ ) { - - int pi; - int global_path_dim = -1; - int mlwd; - - /* Check the dimensionality of each panel */ - for ( pi=0; pi<det->n_panels; pi++ ) { - if ( check_dims(hdfile, &det->panels[pi], - master_el->events[evi], - master_el_with_dims, - &global_path_dim) ) - { - ERROR("Failed to enumerate dims.\n"); - return NULL; - } - } - - /* Add this dimensionality to all events */ - for ( mlwd=0; mlwd<global_path_dim; mlwd++ ) { - - struct event *mlwd_ev; - - mlwd_ev = copy_event(master_el->events[evi]); - push_dim_entry_to_event(mlwd_ev, mlwd); - append_event_to_event_list(master_el_with_dims, - mlwd_ev); - free_event(mlwd_ev); - } - - } - - free_event_list(master_el); - return master_el_with_dims; - - } else { - - return master_el; - - } -} diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h index 99f231e1..eab7b55d 100644 --- a/libcrystfel/src/hdf5-file.h +++ b/libcrystfel/src/hdf5-file.h @@ -42,7 +42,6 @@ struct event_list; #include "image.h" #include "events.h" -struct hdfile; struct copy_hdf5_field; #include "image.h" @@ -54,56 +53,21 @@ extern "C" { /** * \file hdf5-file.h - * HDF5 abstraction layer + * HDF5 utility functions */ -extern int hdf5_write(const char *filename, const void *data, - int width, int height, int type); - extern int hdf5_write_image(const char *filename, const struct image *image, char *element); -extern int hdf5_read(struct hdfile *f, struct image *image, - const char* element, int satcorr); - -extern int hdf5_read2(struct hdfile *f, struct image *image, - struct event *ev, int satcorr); - extern int check_path_existence(hid_t fh, const char *path); -extern struct hdfile *hdfile_open(const char *filename); -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 int get_peaks(struct image *image, struct hdfile *f, const char *p); - -extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p, - int half_pixel_shift); - -extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe); - -extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p, - struct filename_plus_event *fpe, - int half_pixel_shift); - 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, struct event *ev); extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme, const char *name); -extern struct event_list *fill_event_list(struct hdfile* hdfile, - struct detector* det); - extern int hdfile_get_value(struct hdfile *f, const char *name, struct event *ev, void *val, hid_t memtype); extern int hdfile_is_scalar(struct hdfile *f, const char *name, int verbose); diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 3ba71118..c7d46978 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -48,14 +48,6 @@ /** \file image.h */ -struct imagefile -{ - enum imagefile_type type; - char *filename; - struct hdfile *hdfile; -}; - - struct _imagefeaturelist { struct imagefeature *features; @@ -331,233 +323,6 @@ void free_all_crystals(struct image *image) } -/**************************** 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 ************************************/ - -static int unpack_panels(struct image *image, float *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 ( isnan(data[idx]) || isinf(data[idx]) ) 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 void add_out(float val, float *data_out, int nmemb_out, int *outpos, int *nrej) { @@ -1017,75 +782,6 @@ static float *read_cbf_data(const char *filename, int gz, int *w, int *h) } -static int read_cbf(struct imagefile *f, struct image *image) -{ - float *data; - int w, h; - - data = read_cbf_data(f->filename, f->type == IMAGEFILE_CBF, - &w, &h); - if ( data == NULL ) { - ERROR("Failed to read CBF data\n"); - return 1; - } - - unpack_panels(image, data, w, h); - 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); - - return 0; -} - - -static int read_cbf_simple(struct imagefile *f, struct image *image) -{ - float *data; - int w, h; - - data = read_cbf_data(f->filename, f->type == IMAGEFILE_CBF, - &w, &h); - if ( data == NULL ) { - ERROR("Failed to read CBF data\n"); - return 1; - } - - image->det = simple_geometry(image, w, h); - image->dp = malloc(sizeof(float *)); - if ( image->dp == NULL ) { - ERROR("Failed to allocate dp array\n"); - return 1; - } - - image->dp[0] = 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); - - return 0; -} - - -/****************************** Image files ***********************************/ - - signed int is_cbf_file(const char *filename) { FILE *fh; @@ -1123,141 +819,6 @@ signed int is_cbfgz_file(const char *filename) } -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; - - } else if ( is_cbfgz_file(filename) ) { - - f->type = IMAGEFILE_CBFGZ; - - } else { - - ERROR("Unrecognised file type: %s\n", filename); - return NULL; - - } - - 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 if ( f->type == IMAGEFILE_CBFGZ ) { - 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 if ( f->type == IMAGEFILE_CBF ) { - return read_cbf_simple(f, image); - } else if ( f->type == IMAGEFILE_CBFGZ ) { - return read_cbf_simple(f, image); - } else { - ERROR("Unknown file type %i\n", f->type); - return 1; - } -} - - -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 == NULL ) return NULL; - - 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); -} - - /************************** New API (DataTemplate) ****************************/ struct image *image_new() diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index c6a0a561..47d2f414 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -234,7 +234,7 @@ static char *friendly_indexer_name(IndexingMethod m) static void *prepare_method(IndexingMethod *m, UnitCell *cell, - struct detector *det, struct beam_params *beam, + const DataTemplate *dtempl, struct xgandalf_options *xgandalf_opts, struct pinkIndexer_options* pinkIndexer_opts, struct felix_options *felix_opts, @@ -284,7 +284,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, case INDEXING_PINKINDEXER : priv = pinkIndexer_prepare(m, cell, pinkIndexer_opts, - det, beam); + dtempl); break; default : @@ -413,7 +413,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, int j; ipriv->engine_private[i] = prepare_method(&methods[i], cell, - det, beam, + dtempl, xgandalf_opts, pinkIndexer_opts, felix_opts, diff --git a/libcrystfel/src/pinkindexer.c b/libcrystfel/src/pinkindexer.c index 9cf649c1..39913284 100644 --- a/libcrystfel/src/pinkindexer.c +++ b/libcrystfel/src/pinkindexer.c @@ -166,7 +166,7 @@ int run_pinkIndexer(struct image *image, void *ipriv) void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell, struct pinkIndexer_options *pinkIndexer_opts, - struct detector *det, struct beam_params *beam) + const DataTemplate *dtempl) { if ( beam->photon_energy_from != NULL && pinkIndexer_opts->customPhotonEnergy <= 0) { ERROR("For pinkIndexer, the photon_energy must be defined as a " @@ -378,7 +378,7 @@ int run_pinkIndexer(struct image *image, void *ipriv) extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell, struct pinkIndexer_options *pinkIndexer_opts, - struct detector *det, struct beam_params *beam) + const DataTemplate *dtempl) { ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n"); ERROR("To use PINKINDEXER indexing, recompile with PINKINDEXER.\n"); diff --git a/libcrystfel/src/pinkindexer.h b/libcrystfel/src/pinkindexer.h index f79f4331..9d3e2496 100644 --- a/libcrystfel/src/pinkindexer.h +++ b/libcrystfel/src/pinkindexer.h @@ -54,12 +54,13 @@ struct pinkIndexer_options { #include <stddef.h> #include "index.h" +#include "datatemplate.h" extern int run_pinkIndexer(struct image *image, void *ipriv); extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell, struct pinkIndexer_options *pinkIndexer_opts, - struct detector *det, struct beam_params *beam); + const DataTemplate *dtempl); extern void pinkIndexer_cleanup(void *pp); diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index ccf1a632..8255a711 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -857,10 +857,7 @@ int write_chunk(Stream *st, struct image *i, fprintf(st->fh, CHUNK_START_MARKER"\n"); fprintf(st->fh, "Image filename: %s\n", i->filename); - if ( i->event != NULL ) { - fprintf(st->fh, "Event: %s\n", get_event_string(i->event)); - } - + fprintf(st->fh, "Event: %s\n", i->ev); fprintf(st->fh, "Image serial number: %i\n", i->serial); fprintf(st->fh, "hit = %i\n", i->hit); @@ -1087,8 +1084,7 @@ static void read_crystal(Stream *st, struct image *image, StreamReadFlags srf) if ( reflist == NULL ) { ERROR("Failed while reading reflections\n"); ERROR("Filename = %s\n", image->filename); - ERROR("Event = %s\n", - get_event_string(image->event)); + ERROR("Event = %s\n", image->ev); break; } @@ -1200,7 +1196,7 @@ static int read_and_store_field(struct image *image, const char *line) /** * Read the next chunk from a stream and fill in 'image' */ -int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf) +int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf) { char line[1024]; char *rval = NULL; @@ -1213,7 +1209,7 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf) image->features = NULL; image->crystals = NULL; image->n_crystals = 0; - image->event = NULL; + image->ev = NULL; image->stuff_from_stream = NULL; if ( (srf & STREAM_READ_REFLECTIONS) || (srf & STREAM_READ_UNITCELL) ) { @@ -1238,7 +1234,7 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf) } if ( strncmp(line, "Event: ", 7) == 0 ) { - image->event = get_event_from_event_string(line+7); + image->ev = strdup(line+7); } if ( strncmp(line, "indexed_by = ", 13) == 0 ) { diff --git a/src/im-zmq.c b/src/im-zmq.c index f8743dbb..39f1d207 100644 --- a/src/im-zmq.c +++ b/src/im-zmq.c @@ -46,6 +46,7 @@ #include "hdf5-file.h" #include "utils.h" #include "im-zmq.h" +#include "datatemplate_priv.h" struct im_zmq @@ -408,19 +409,19 @@ static double *find_msgpack_data(msgpack_object *obj, int *width, int *height) } -static double *zero_array(DataTemplate *dtempl, int *dw, int *dh) +static double *zero_array(const DataTemplate *dtempl, int *dw, int *dh) { int max_fs = 0; int max_ss = 0; int pi; double *data; - for ( pi=0; pi<det->n_panels; pi++ ) { - if ( det->panels[pi].orig_max_fs > max_fs ) { - max_fs = det->panels[pi].orig_max_fs; + for ( pi=0; pi<dtempl->n_panels; pi++ ) { + if ( dtempl->panels[pi].orig_max_fs > max_fs ) { + max_fs = dtempl->panels[pi].orig_max_fs; } - if ( det->panels[pi].orig_max_ss > max_ss ) { - max_ss = det->panels[pi].orig_max_ss; + if ( dtempl->panels[pi].orig_max_ss > max_ss ) { + max_ss = dtempl->panels[pi].orig_max_ss; } } @@ -456,25 +457,25 @@ struct image *unpack_msgpack_data(msgpack_object *obj, if ( obj == NULL ) { ERROR("No MessagePack object!\n"); - return 1; + return NULL; } if ( !no_image_data ) { data = find_msgpack_data(obj, &data_width, &data_height); if ( data == NULL ) { ERROR("No image data in MessagePack object.\n"); - return 1; + return NULL; } } else { data = zero_array(dtempl, &data_width, &data_height); } image = image_new(); - if ( image == NULL ) return 1; + if ( image == NULL ) return NULL; if ( unpack_slab(image, data, data_width, data_height) ) { ERROR("Failed to unpack data slab.\n"); - return 1; + return NULL; } im_zmq_fill_in_beam_parameters(image->beam, image); @@ -485,5 +486,5 @@ struct image *unpack_msgpack_data(msgpack_object *obj, im_zmq_fill_in_clen(image->det); fill_in_adu(image); - return 0; + return image; } diff --git a/src/im-zmq.h b/src/im-zmq.h index 270a7c3a..ddf3f172 100644 --- a/src/im-zmq.h +++ b/src/im-zmq.h @@ -55,7 +55,7 @@ extern ImageFeatureList *get_peaks_msgpack(msgpack_object *obj, int half_pixel_shift); extern struct image *unpack_msgpack_data(msgpack_object *obj, - DataTemplate *dtempl, + const DataTemplate *dtempl, int no_image_data); #else /* defined(HAVE_MSGPACK) && defined(HAVE_ZMQ) */ diff --git a/src/process_image.c b/src/process_image.c index d529a444..bca2a7f2 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -130,7 +130,7 @@ static struct image *file_wait_open_read(const char *filename, if ( !wait_message_done ) { STATUS("Waiting for '%s'\n", - image->filename); + filename); wait_message_done = 1; } @@ -194,6 +194,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, if ( pargs->msgpack_obj != NULL ) { set_last_task(last_task, "unpacking messagepack object"); image = unpack_msgpack_data(pargs->msgpack_obj, + iargs->dtempl, iargs->no_image_data); if ( image == NULL ) return; } else { |