Add image_write()

This is a one-to-one transfer of the old hdf5_write_image()

5 files changed, 359 insertions, 0 deletions

diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h
index 15495cef..3afbcd0b 100644
--- a/libcrystfel/src/datatemplate_priv.h
+++ b/libcrystfel/src/datatemplate_priv.h
@@ -159,6 +159,10 @@ struct panel_template
 };
 
 
+#define PANEL_WIDTH(p) ((p)->orig_max_fs - (p)->orig_min_fs + 1)
+#define PANEL_HEIGHT(p) ((p)->orig_max_ss - (p)->orig_min_ss + 1)
+
+
 struct dt_badregion
 {
 	char name[1024];
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;
+}
diff --git a/libcrystfel/src/image-hdf5.h b/libcrystfel/src/image-hdf5.h
index 07cdf96e..45cbe7bd 100644
--- a/libcrystfel/src/image-hdf5.h
+++ b/libcrystfel/src/image-hdf5.h
@@ -64,4 +64,8 @@ extern char **image_hdf5_expand_frames(const DataTemplate *dtempl,
 
 extern int is_hdf5_file(const char *filename);
 
+extern int image_hdf5_write(const struct image *image,
+                            const DataTemplate *dtempl,
+                            const char *filename);
+
 #endif	/* IMAGE_HDF5_H */
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index d2fb30ba..c8d0f329 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -916,3 +916,16 @@ void mark_resolution_range_as_bad(struct image *image,
 
 	}
 }
+
+
+int image_write(const struct image *image,
+                const DataTemplate *dtempl,
+                const char *filename)
+{
+	if ( is_hdf5_file(filename) ) {
+		return image_hdf5_write(image, dtempl, filename);
+	}
+
+	ERROR("Can only write to HDF5 files.\n");
+	return 1;
+}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index a2f160af..3d42322a 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -208,6 +208,10 @@ extern int image_set_zero_data(struct image *image,
 extern int image_set_zero_mask(struct image *image,
                                const DataTemplate *dtempl);
 
+extern int image_write(const struct image *image,
+                       const DataTemplate *dtempl,
+                       const char *filename);
+
 /* Use within libcrystfel only */
 extern int create_detgeom(struct image *image,
                           const DataTemplate *dtempl);