Implement loading of saturation map

3 files changed, 100 insertions, 1 deletions

diff --git a/libcrystfel/src/image-hdf5.c b/libcrystfel/src/image-hdf5.c
index 25e9d3f7..7af0e3a6 100644
--- a/libcrystfel/src/image-hdf5.c
+++ b/libcrystfel/src/image-hdf5.c
@@ -578,6 +578,27 @@ int image_hdf5_read(struct image *image,
 }
 
 
+float *image_hdf5_read_satmap(struct panel_template *p,
+                              const char *filename,
+                              const char *event,
+                              const char *map_location)
+{
+	float *map = NULL;
+
+	if ( load_hdf5_hyperslab(p, filename, event,
+	                         (void *)&map, H5T_NATIVE_FLOAT,
+	                         sizeof(float), 1, map_location) )
+	{
+		ERROR("Failed to load saturation map data\n");
+		free(map);
+		return NULL;
+	}
+
+	return map;
+}
+
+
+
 int image_hdf5_read_mask(struct panel_template *p,
                          const char *filename, const char *event,
                          int *bad, const char *mask_location,
diff --git a/libcrystfel/src/image-hdf5.h b/libcrystfel/src/image-hdf5.h
index 96b21558..4c36cf89 100644
--- a/libcrystfel/src/image-hdf5.h
+++ b/libcrystfel/src/image-hdf5.h
@@ -50,6 +50,11 @@ extern int image_hdf5_read_mask(struct panel_template *p,
                                 const char *mask_location,
                                 int mask_good, int mask_bad);
 
+extern float *image_hdf5_read_satmap(struct panel_template *p,
+                                     const char *filename,
+                                     const char *event,
+                                     const char *map_location);
+
 extern ImageFeatureList *image_hdf5_read_peaks_cxi(const DataTemplate *dtempl,
                                                    const char *filename,
                                                    const char *event,
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 8d9f1640..1423c4d3 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -967,7 +967,80 @@ static int create_badmap(struct image *image,
 static int create_satmap(struct image *image,
                          const DataTemplate *dtempl)
 {
-	/* FIXME: Implementation */
+	int i;
+	int any;
+
+	/* The panels will be treated separately, but we'll only bother at all
+	 * if at least one of them has a saturation map. */
+	any = 0;
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		if ( dtempl->panels[i].satmap != NULL ) {
+			any = 1;
+			break;
+		}
+	}
+
+	if ( !any ) return 0;
+
+	image->sat = malloc(dtempl->n_panels * sizeof(float *));
+	if ( image->sat == NULL ) {
+		ERROR("Failed to allocate saturation map\n");
+		return 1;
+	}
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+
+		struct panel_template *p = &dtempl->panels[i];
+
+		if ( p->satmap == NULL ) {
+
+			/* At least one other panel has a saturation map,
+			 * but it isn't this one.  Therefore make a fake
+			 * saturation map */
+
+			long int j;
+			int p_w, p_h;
+
+			p_w = p->orig_max_fs - p->orig_min_fs + 1;
+			p_h = p->orig_max_ss - p->orig_min_ss + 1;
+
+			image->sat[i] = malloc(p_w*p_h*sizeof(float));
+
+			if ( image->sat[i] != NULL ) {
+				for ( j=0; j<p_w*p_h; j++ ) {
+					image->sat[i][j] = INFINITY;
+				}
+			}
+
+		} else {
+
+			const char *map_fn;
+
+			if ( p->satmap_file == NULL ) {
+				map_fn = image->filename;
+			} else {
+				map_fn = p->satmap_file;
+			}
+
+			if ( is_hdf5_file(map_fn) ) {
+				image->sat[i] = image_hdf5_read_satmap(p, map_fn,
+				                                       image->ev,
+				                                       p->satmap);
+
+			} else {
+				ERROR("Saturation map must be in HDF5 format\n");
+				return 1;
+			}
+		}
+
+		if ( image->sat[i] == NULL ) {
+			ERROR("Failed to allocate saturation map (panel %s)\n",
+			      p->name);
+			return 1;
+		}
+
+	}
+
 	return 0;
 }