Merge branch 'pf8_faster'

10 files changed, 672 insertions, 269 deletions

diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c
index 1ffa2777..bbcf1aa5 100644
--- a/libcrystfel/src/datatemplate.c
+++ b/libcrystfel/src/datatemplate.c
@@ -36,6 +36,7 @@
 
 #include "utils.h"
 #include "datatemplate.h"
+#include "image.h"
 
 #include "datatemplate_priv.h"
 
@@ -1734,3 +1735,364 @@ double data_template_get_wavelength_if_possible(const DataTemplate *dt)
 		return NAN;
 	}
 }
+
+
+static int separate_value_and_units(const char *from,
+                                    char **pvalue,
+                                    char **punits)
+{
+	char *sp;
+	char *fromcpy;
+	char *unitscpy;
+
+	if ( from == NULL ) return 1;
+
+	fromcpy = strdup(from);
+	if ( fromcpy == NULL ) return 1;
+
+	sp = strchr(fromcpy, ' ');
+	if ( sp == NULL ) {
+		unitscpy = NULL;
+	} else {
+		unitscpy = strdup(sp+1);
+		sp[0] = '\0';
+	}
+
+	*pvalue = fromcpy;
+	*punits = unitscpy;
+	return 0;
+}
+
+
+/* default_scale is a value to be used if both of the following
+ * conditions are met:
+ *
+ *  1. The value is a reference to image headers/metadata,
+ *      rather than a literal number.
+ *  2. No units are specified in the number.
+ *
+ * This is totally horrible.  Sorry.  Blame history.
+ */
+static int im_get_length(struct image *image, const char *from,
+                         double default_scale, double *pval)
+{
+	char *value_str;
+	char *units;
+
+	if ( from == NULL ) return 1;
+
+	if ( separate_value_and_units(from, &value_str, &units) ) return 1;
+
+	if ( units == NULL ) {
+
+		/* No units given */
+
+		if ( convert_float(value_str, pval) == 0 ) {
+
+			/* Literal value with no units */
+			free(value_str);
+			return 0;
+
+		} else {
+
+			int r;
+			r = image_read_header_float(image, value_str, pval);
+			free(value_str);
+
+			if ( r == 0 ) {
+				/* Value read from headers with no units */
+				*pval *= default_scale;
+				return 0;
+			} else {
+				/* Failed to read value from headers */
+				return 1;
+			}
+		}
+
+	} else {
+
+		/* Units are specified */
+
+		double scale;
+
+		if ( strcmp(units, "mm") == 0 ) {
+			scale = 1e-3;
+		} else if ( strcmp(units, "m") == 0 ) {
+			scale = 1.0;
+		} else {
+			ERROR("Invalid length unit '%s'\n", units);
+			free(value_str);
+			free(units);
+			return 1;
+		}
+
+		if ( convert_float(value_str, pval) == 0 ) {
+
+			/* Literal value, units specified */
+			free(value_str);
+			free(units);
+			*pval *= scale;
+			return 0;
+
+		} else {
+
+			int r;
+			r = image_read_header_float(image, value_str, pval);
+			free(value_str);
+
+			if ( r == 0 ) {
+				/* Value read from headers, units specified */
+				*pval *= scale;
+				return 0;
+			} else {
+				/* Failed to read value from headers */
+				return 1;
+			}
+		}
+	}
+}
+
+
+static int safe_strcmp(const char *a, const char *b)
+{
+	if ( (a==NULL) && (b==NULL) ) return 0;
+	if ( (a!=NULL) && (b!=NULL) ) return strcmp(a, b);
+	return 1;
+}
+
+
+static int all_panels_reference_same_clen(const DataTemplate *dtempl)
+{
+	int i;
+	char *first_val = NULL;
+	char *first_units = NULL;
+	int fail = 0;
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		struct panel_template *p = &dtempl->panels[i];
+		char *val;
+		char *units;
+		if ( separate_value_and_units(p->cnz_from, &val, &units) ) {
+			/* Parse error */
+			return 0;
+		}
+		if ( i == 0 ) {
+			first_val = val;
+			first_units = units;
+		} else {
+			if ( safe_strcmp(val, first_val) != 0 ) fail = 1;
+			if ( safe_strcmp(units, first_units) != 0 ) fail = 1;
+			free(val);
+			free(units);
+		}
+	}
+
+	free(first_val);
+	free(first_units);
+	return fail;
+}
+
+
+static int all_coffsets_small(const DataTemplate *dtempl)
+{
+	int i;
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		struct panel_template *p = &dtempl->panels[i];
+		if ( p->cnz_offset > 10.0*p->pixel_pitch ) return 0;
+	}
+
+	return 1;
+}
+
+
+static int all_panels_same_clen(const DataTemplate *dtempl)
+{
+	int i;
+	double *zvals;
+	double total = 0.0;
+	double mean;
+
+	zvals = malloc(sizeof(double)*dtempl->n_panels);
+	if ( zvals == NULL ) return 0;
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		struct panel_template *p = &dtempl->panels[i];
+		if ( im_get_length(NULL, p->cnz_from, 1e-3, &zvals[i]) ) {
+			/* Can't get length because it used a header reference */
+			free(zvals);
+			return 0;
+		}
+		total += zvals[i];
+	}
+
+	mean = total/dtempl->n_panels;
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		struct panel_template *p = &dtempl->panels[i];
+		if ( fabs(zvals[i] - mean) > 10.0*p->pixel_pitch ) return 0;
+	}
+
+	free(zvals);
+
+	return 1;
+}
+
+
+static int all_panels_perpendicular_to_beam(const DataTemplate *dtempl)
+{
+	int i;
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+		double z_diff;
+		struct panel_template *p = &dtempl->panels[i];
+		z_diff = p->fsz*PANEL_WIDTH(p) + p->ssz*PANEL_HEIGHT(p);
+		if ( z_diff > 10.0*p->pixel_pitch ) return 0;
+	}
+	return 1;
+}
+
+
+static int detector_flat(const DataTemplate *dtempl)
+{
+	return all_panels_perpendicular_to_beam(dtempl)
+	    && ( (all_panels_reference_same_clen(dtempl) && all_coffsets_small(dtempl))
+	          || all_panels_same_clen(dtempl) );
+}
+
+
+struct detgeom *create_detgeom(struct image *image,
+                               const DataTemplate *dtempl,
+                               int two_d_only)
+{
+	struct detgeom *detgeom;
+	int i;
+
+	if ( dtempl == NULL ) {
+		ERROR("NULL data template!\n");
+		return NULL;
+	}
+
+	detgeom = malloc(sizeof(struct detgeom));
+	if ( detgeom == NULL ) return NULL;
+
+	detgeom->panels = malloc(dtempl->n_panels*sizeof(struct detgeom_panel));
+	if ( detgeom->panels == NULL ) {
+		free(detgeom);
+		return NULL;
+	}
+
+	detgeom->n_panels = dtempl->n_panels;
+
+	if ( two_d_only ) {
+		if ( !detector_flat(dtempl) ) return NULL;
+		if ( dtempl->shift_x_from != NULL ) return NULL;
+		if ( dtempl->shift_y_from != NULL ) return NULL;
+	}
+
+	for ( i=0; i<dtempl->n_panels; i++ ) {
+
+		struct detgeom_panel *p = &detgeom->panels[i];
+		struct panel_template *tmpl = &dtempl->panels[i];
+		double shift_x, shift_y;
+
+		p->name = safe_strdup(tmpl->name);
+
+		p->pixel_pitch = tmpl->pixel_pitch;
+
+		/* NB cnx,cny are in pixels, cnz is in m */
+		p->cnx = tmpl->cnx;
+		p->cny = tmpl->cny;
+
+		if ( im_get_length(image, tmpl->cnz_from, 1e-3, &p->cnz) )
+		{
+			if ( two_d_only ) {
+				p->cnz = NAN;
+			} else {
+				ERROR("Failed to read length from '%s'\n", tmpl->cnz_from);
+				return NULL;
+			}
+		}
+
+		/* Apply offset (in m) and then convert cnz from
+		 * m to pixels */
+		p->cnz += tmpl->cnz_offset;
+		p->cnz /= p->pixel_pitch;
+
+		/* Apply overall shift (already in m) */
+		if ( dtempl->shift_x_from != NULL ) {
+			if ( im_get_length(image, dtempl->shift_x_from, 1.0, &shift_x) ) {
+				ERROR("Failed to read length from '%s'\n",
+				      dtempl->shift_x_from);
+				return NULL;
+			}
+			if ( im_get_length(image, dtempl->shift_y_from, 1.0, &shift_y) ) {
+				ERROR("Failed to read length from '%s'\n",
+				      dtempl->shift_y_from);
+				return NULL;
+			}
+		} else {
+			shift_x = 0.0;
+			shift_y = 0.0;
+		}
+
+		if ( !isnan(shift_x) ) {
+			p->cnx += shift_x / p->pixel_pitch;
+		}
+		if ( !isnan(shift_y) ) {
+			p->cny += shift_y / p->pixel_pitch;
+		}
+
+		p->max_adu = tmpl->max_adu;
+
+		switch ( tmpl->adu_scale_unit ) {
+
+			case ADU_PER_PHOTON:
+			p->adu_per_photon = tmpl->adu_scale;
+			break;
+
+			case ADU_PER_EV:
+			p->adu_per_photon = tmpl->adu_scale
+				* ph_lambda_to_eV(image->lambda);
+			break;
+
+			default:
+			p->adu_per_photon = 1.0;
+			ERROR("Invalid ADU/ph scale unit (%i)\n",
+			      tmpl->adu_scale_unit);
+			break;
+
+		}
+
+		p->w = tmpl->orig_max_fs - tmpl->orig_min_fs + 1;
+		p->h = tmpl->orig_max_ss - tmpl->orig_min_ss + 1;
+
+		p->fsx = tmpl->fsx;
+		p->fsy = tmpl->fsy;
+		p->fsz = tmpl->fsz;
+		p->ssx = tmpl->ssx;
+		p->ssy = tmpl->ssy;
+		p->ssz = tmpl->ssz;
+
+	}
+
+	return detgeom;
+}
+
+
+/**
+ * Create a detgeom structure from the DataTemplate, if possible, and ignoring
+ * 3D information.
+ *
+ * This procedure will create a detgeom structure provided that the detector
+ *  is close to lying in a single flat plane perpendicular to the beam
+ *  direction.  If certain things (e.g. panel z-positions) refer to headers,
+ *  it might not be possible to determine that the detector is really flat
+ *  until an image is loaded.  Therefore you must gracefully handle a NULL
+ *  return value from this routine.
+ *
+ * \returns the detgeom structure, or NULL if impossible.
+ */
+struct detgeom *data_template_get_2d_detgeom_if_possible(const DataTemplate *dt)
+{
+	return create_detgeom(NULL, dt, 1);
+}
diff --git a/libcrystfel/src/datatemplate.h b/libcrystfel/src/datatemplate.h
index 28cabcf5..ddc210e6 100644
--- a/libcrystfel/src/datatemplate.h
+++ b/libcrystfel/src/datatemplate.h
@@ -96,6 +96,8 @@ extern struct rg_collection *data_template_get_rigid_groups(const DataTemplate *
 
 extern double data_template_get_wavelength_if_possible(const DataTemplate *dt);
 
+extern struct detgeom *data_template_get_2d_detgeom_if_possible(const DataTemplate *dt);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h
index ff383705..26b90d6b 100644
--- a/libcrystfel/src/datatemplate_priv.h
+++ b/libcrystfel/src/datatemplate_priv.h
@@ -32,6 +32,8 @@
 #ifndef DATATEMPLATE_PRIV_H
 #define DATATEMPLATE_PRIV_H
 
+#include "detgeom.h"
+
 /* Maximum number of dimensions expected in data files */
 #define MAX_DIMS (16)
 
@@ -234,5 +236,8 @@ struct _datatemplate
 };
 
 extern double convert_to_m(double val, int units);
+extern struct detgeom *create_detgeom(struct image *image,
+                                      const DataTemplate *dtempl,
+                                      int two_d_only);
 
 #endif	/* DATATEMPLATE_PRIV_H */
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 3448421e..84088ead 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -412,11 +412,6 @@ static struct header_cache_entry *cached_header(struct image *image, const char
 {
 	struct header_cache_entry *ce;
 
-	if ( image == NULL ) {
-		ERROR("Attempt to retrieve a header value without an image\n");
-		return NULL;
-	}
-
 	ce = find_cache_entry(image, from);
 	if ( ce != NULL ) return ce;
 
@@ -434,6 +429,8 @@ int image_read_header_float(struct image *image, const char *from, double *val)
 {
 	struct header_cache_entry *ce;
 
+	if ( image == NULL ) return 1;
+
 	ce = cached_header(image, from);
 	if ( ce == NULL ) return 1;
 
@@ -486,231 +483,6 @@ static DataSourceType file_type(const char *filename)
 }
 
 
-static int separate_value_and_units(const char *from,
-                                    char **pvalue,
-                                    char **punits)
-{
-	char *sp;
-	char *fromcpy;
-	char *unitscpy;
-
-	if ( from == NULL ) return 1;
-
-	fromcpy = strdup(from);
-	if ( fromcpy == NULL ) return 1;
-
-	sp = strchr(fromcpy, ' ');
-	if ( sp == NULL ) {
-		unitscpy = NULL;
-	} else {
-		unitscpy = strdup(sp+1);
-		sp[0] = '\0';
-	}
-
-	*pvalue = fromcpy;
-	*punits = unitscpy;
-	return 0;
-}
-
-
-/* default_scale is a value to be used if both of the following
- * conditions are met:
- *
- *  1. The value is a reference to image headers/metadata,
- *      rather than a literal number.
- *  2. No units are specified in the number.
- *
- * This is totally horrible.  Sorry.  Blame history.
- */
-static int im_get_length(struct image *image, const char *from,
-                         double default_scale, double *pval)
-{
-	char *value_str;
-	char *units;
-
-	if ( from == NULL ) return 1;
-
-	if ( separate_value_and_units(from, &value_str, &units) ) return 1;
-
-	if ( units == NULL ) {
-
-		/* No units given */
-
-		if ( convert_float(value_str, pval) == 0 ) {
-
-			/* Literal value with no units */
-			free(value_str);
-			return 0;
-
-		} else {
-
-			int r;
-			r = image_read_header_float(image, value_str, pval);
-			free(value_str);
-
-			if ( r == 0 ) {
-				/* Value read from headers with no units */
-				*pval *= default_scale;
-				return 0;
-			} else {
-				/* Failed to read value from headers */
-				return 1;
-			}
-		}
-
-	} else {
-
-		/* Units are specified */
-
-		double scale;
-
-		if ( strcmp(units, "mm") == 0 ) {
-			scale = 1e-3;
-		} else if ( strcmp(units, "m") == 0 ) {
-			scale = 1.0;
-		} else {
-			ERROR("Invalid length unit '%s'\n", units);
-			free(value_str);
-			free(units);
-			return 1;
-		}
-
-		if ( convert_float(value_str, pval) == 0 ) {
-
-			/* Literal value, units specified */
-			free(value_str);
-			free(units);
-			*pval *= scale;
-			return 0;
-
-		} else {
-
-			int r;
-			r = image_read_header_float(image, value_str, pval);
-			free(value_str);
-
-			if ( r == 0 ) {
-				/* Value read from headers, units specified */
-				*pval *= scale;
-				return 0;
-			} else {
-				/* Failed to read value from headers */
-				return 1;
-			}
-		}
-	}
-}
-
-
-int create_detgeom(struct image *image, const DataTemplate *dtempl)
-{
-	struct detgeom *detgeom;
-	int i;
-
-	if ( dtempl == NULL ) {
-		ERROR("NULL data template!\n");
-		return 1;
-	}
-
-	detgeom = malloc(sizeof(struct detgeom));
-	if ( detgeom == NULL ) return 1;
-
-	detgeom->panels = malloc(dtempl->n_panels*sizeof(struct detgeom_panel));
-	if ( detgeom->panels == NULL ) {
-		free(detgeom);
-		return 1;
-	}
-
-	detgeom->n_panels = dtempl->n_panels;
-
-	for ( i=0; i<dtempl->n_panels; i++ ) {
-
-		struct detgeom_panel *p = &detgeom->panels[i];
-		struct panel_template *tmpl = &dtempl->panels[i];
-		double shift_x, shift_y;
-
-		p->name = safe_strdup(tmpl->name);
-
-		p->pixel_pitch = tmpl->pixel_pitch;
-
-		/* NB cnx,cny are in pixels, cnz is in m */
-		p->cnx = tmpl->cnx;
-		p->cny = tmpl->cny;
-		if ( im_get_length(image, tmpl->cnz_from,
-		                   1e-3, &p->cnz) )
-		{
-			ERROR("Failed to read length from '%s'\n", tmpl->cnz_from);
-			return 1;
-		}
-
-		/* Apply offset (in m) and then convert cnz from
-		 * m to pixels */
-		p->cnz += tmpl->cnz_offset;
-		p->cnz /= p->pixel_pitch;
-
-		/* Apply overall shift (already in m) */
-		if ( dtempl->shift_x_from != NULL ) {
-			if ( im_get_length(image, dtempl->shift_x_from, 1.0, &shift_x) ) {
-				ERROR("Failed to read length from '%s'\n",
-				      dtempl->shift_x_from);
-				return 1;
-			}
-			if ( im_get_length(image, dtempl->shift_y_from, 1.0, &shift_y) ) {
-				ERROR("Failed to read length from '%s'\n",
-				      dtempl->shift_y_from);
-				return 1;
-			}
-		} else {
-			shift_x = 0.0;
-			shift_y = 0.0;
-		}
-
-		if ( !isnan(shift_x) ) {
-			p->cnx += shift_x / p->pixel_pitch;
-		}
-		if ( !isnan(shift_y) ) {
-			p->cny += shift_y / p->pixel_pitch;
-		}
-
-		p->max_adu = tmpl->max_adu;
-
-		switch ( tmpl->adu_scale_unit ) {
-
-			case ADU_PER_PHOTON:
-			p->adu_per_photon = tmpl->adu_scale;
-			break;
-
-			case ADU_PER_EV:
-			p->adu_per_photon = tmpl->adu_scale
-				* ph_lambda_to_eV(image->lambda);
-			break;
-
-			default:
-			p->adu_per_photon = 1.0;
-			ERROR("Invalid ADU/ph scale unit (%i)\n",
-			      tmpl->adu_scale_unit);
-			break;
-
-		}
-
-		p->w = tmpl->orig_max_fs - tmpl->orig_min_fs + 1;
-		p->h = tmpl->orig_max_ss - tmpl->orig_min_ss + 1;
-
-		p->fsx = tmpl->fsx;
-		p->fsy = tmpl->fsy;
-		p->fsz = tmpl->fsz;
-		p->ssx = tmpl->ssx;
-		p->ssy = tmpl->ssy;
-		p->ssz = tmpl->ssz;
-
-	}
-
-	image->detgeom = detgeom;
-
-	return 0;
-}
-
-
 int image_set_zero_data(struct image *image,
                         const DataTemplate *dtempl)
 {
@@ -1299,7 +1071,8 @@ struct image *image_create_for_simulation(const DataTemplate *dtempl)
 		return NULL;
 	}
 
-	if ( create_detgeom(image, dtempl) ) {
+	image->detgeom = create_detgeom(image, dtempl, 0);
+	if ( image->detgeom == NULL ) {
 		image_free(image);
 		return NULL;
 	}
@@ -1350,9 +1123,9 @@ static int do_image_read(struct image *image, const DataTemplate *dtempl,
 	}
 
 	profile_start("create-detgeom");
-	r = create_detgeom(image, dtempl);
+	image->detgeom = create_detgeom(image, dtempl, 0);
 	profile_end("create-detgeom");
-	if ( r ) {
+	if ( image->detgeom == NULL ) {
 		ERROR("Failed to read geometry information\n");
 		return 1;
 	}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 3f1ed219..10ae0ba7 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -270,9 +270,6 @@ 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);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c
index 537c2365..3420a7a7 100644
--- a/libcrystfel/src/peakfinder8.c
+++ b/libcrystfel/src/peakfinder8.c
@@ -34,6 +34,7 @@
 #include <float.h>
 #include <math.h>
 #include <stdlib.h>
+#include <profile.h>
 
 #include "peakfinder8.h"
 #include "detgeom.h"
@@ -43,13 +44,24 @@
 /** \file peakfinder8.h */
 
 // CrystFEL-only block 1
+
 struct radius_maps
 {
 	float **r_maps;
+	int *n_pixels;
 	int n_rmaps;
 };
 
 
+struct radial_stats_pixels
+{
+	int n_panels;
+	int *n_pixels; // n_pixels[panel]
+	int **pidx;	   // pixel_index[panel][0..n_pixels]
+	int **radius;  // pixel_radius[panel][0..n_pixels]
+};
+
+
 struct peakfinder_mask
 {
 	char **masks;
@@ -101,7 +113,138 @@ struct peakfinder_peak_data
 };
 
 
-// CrystFEL-only block 2
+static struct radial_stats_pixels *compute_rstats_pixels(struct radius_maps *rmaps)
+{
+	int p;
+	int i;
+
+	struct radial_stats_pixels *rsp = NULL;
+	rsp = (struct radial_stats_pixels *)malloc(sizeof(struct radial_stats_pixels));
+	if ( rsp == NULL ) {
+		return NULL;
+	}
+	rsp->n_pixels = (int *)malloc(rmaps->n_rmaps * sizeof(int));
+	if ( rsp->n_pixels == NULL ) {
+		free(rsp);
+		return NULL;
+	}
+	rsp->pidx = (int **)malloc(rmaps->n_rmaps * sizeof(int *));
+	if ( rsp->pidx == NULL ) {
+		free(rsp->n_pixels);
+		free(rsp);
+		return NULL;
+	}
+	rsp->radius = (int **)malloc(rmaps->n_rmaps * sizeof(int *));
+	if ( rsp->radius == NULL ) {
+		free(rsp->n_pixels);
+		free(rsp->pidx);
+		free(rsp);
+		return NULL;
+	}
+	srand(0);
+
+	int n_pixels_per_bin = 100; // Can make this a parameter
+
+	// Assuming 5000 is the maximum possible radius
+	int n_pixels[5000] = {0}; // selected pixels per bin
+	int n_tot_pixels[5000] = {0}; // total pixels per bin
+	int panel[5000][n_pixels_per_bin]; // panel ID of selected pixels
+	int idx[5000][n_pixels_per_bin]; // index of selected pixels
+	int radius;
+
+	for ( p = 0; p < rmaps->n_rmaps; p++ ) {
+		rsp->n_pixels[p] = 0;
+		for ( i = 0; i < rmaps->n_pixels[p]; i++ ) {
+			// Reservoir sampling:
+			radius = (int)rint(rmaps->r_maps[p][i]);
+			n_tot_pixels[radius] += 1;
+
+			if ( n_pixels[radius] < n_pixels_per_bin ) {
+				panel[radius][n_pixels[radius]] = p;
+				idx[radius][n_pixels[radius]] = i;
+
+				n_pixels[radius] += 1;
+				rsp->n_pixels[p] += 1;
+			} else {
+				int rand_i = rand() % n_tot_pixels[radius];
+				if ( rand_i < n_pixels_per_bin ) {
+					rsp->n_pixels[panel[radius][rand_i]] -= 1;
+					rsp->n_pixels[p] += 1;
+
+					panel[radius][rand_i] = p;
+					idx[radius][rand_i] = i;
+				}
+			}
+		}
+	}
+
+	int *sidx = (int *)malloc(rmaps->n_rmaps * sizeof(int));
+	if ( sidx == NULL ) {
+		free(rsp->n_pixels);
+		free(rsp->pidx);
+		free(rsp->radius);
+		free(rsp);
+		return NULL;
+	}
+	for ( p = 0; p < rmaps->n_rmaps; p++ ) {
+		rsp->pidx[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int));
+		if ( rsp->pidx[p] == NULL ) {
+			for ( i = 0; i < p; i++ ) {
+				free(rsp->pidx[i]);
+				free(rsp->radius[i]);
+			}
+			free(rsp->pidx);
+			free(rsp->radius);
+			free(rsp->n_pixels);
+			free(rsp);
+			free(sidx);
+			return NULL;
+		}
+		rsp->radius[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int));
+		if ( rsp->radius[p] == NULL ) {
+			for ( i = 0; i < p; i++ ) {
+				free(rsp->pidx[i]);
+				free(rsp->radius[i]);
+			}
+			free(rsp->pidx[p]);
+			free(rsp->pidx);
+			free(rsp->radius);
+			free(rsp->n_pixels);
+			free(rsp);
+			free(sidx);
+			return NULL;
+		}
+		sidx[p] = 0;
+	}
+
+	for ( radius = 0; radius < 5000; radius++ ) {
+		for ( i = 0; i < n_pixels[radius]; i++ ) {
+			p = panel[radius][i];
+			rsp->pidx[p][sidx[p]] = idx[radius][i];
+			rsp->radius[p][sidx[p]] = radius;
+			sidx[p] += 1;
+		}
+	}
+	free(sidx);
+
+	rsp->n_panels = rmaps->n_rmaps;
+	return rsp;
+}
+
+static void free_rstats_pixels(struct radial_stats_pixels *rsp)
+{
+	int i;
+	for ( i = 0; i < rsp->n_panels; i++ ) {
+		free(rsp->pidx[i]);
+		free(rsp->radius[i]);
+	}
+	free(rsp->pidx);
+	free(rsp->radius);
+	free(rsp->n_pixels);
+	free(rsp);
+}
+
+
 static struct radius_maps *compute_radius_maps(struct detgeom *det)
 {
 	int i, u, iss, ifs;
@@ -119,6 +262,12 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det)
 		return NULL;
 	}
 
+	rm->n_pixels = (int *)malloc(det->n_panels*sizeof(int*));
+	if ( rm->r_maps == NULL ) {
+		free(rm);
+		return NULL;
+	}
+
 	rm->n_rmaps = det->n_panels;
 
 	for( i=0 ; i<det->n_panels ; i++ ) {
@@ -133,7 +282,7 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det)
 			free(rm);
 			return NULL;
 		}
-
+		rm->n_pixels[i] = p.h * p.w;
 		for ( iss=0 ; iss<p.h ; iss++ ) {
 			for ( ifs=0; ifs<p.w; ifs++ ) {
 
@@ -161,10 +310,53 @@ static void free_radius_maps(struct radius_maps *r_maps)
 		free(r_maps->r_maps[i]);
 	}
 	free(r_maps->r_maps);
+	free(r_maps->n_pixels);
 	free(r_maps);
 }
 
 
+// CrystFEL-only block 2
+struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode)
+{
+	struct pf8_private_data *data = NULL;
+	if ( det == NULL ) {
+		return NULL;
+	}
+
+	data = (struct pf8_private_data *)malloc(sizeof(struct pf8_private_data));
+	if ( data == NULL ) {
+		return NULL;
+	}
+	data->rmaps = compute_radius_maps(det);
+	if ( data->rmaps == NULL ) {
+		free(data);
+		return NULL;
+	}
+	if ( fast_mode ) {
+		data->rpixels = compute_rstats_pixels(data->rmaps);
+		if ( data->rpixels == NULL ) {
+			free_radius_maps(data->rmaps);
+			free(data);
+			return NULL;
+		}
+	} else {
+		data->rpixels = NULL;
+	}
+	data->fast_mode = fast_mode;
+	return data;
+}
+
+
+void free_pf8_private_data(struct pf8_private_data *data)
+{
+	free_radius_maps(data->rmaps);
+	if ( data->fast_mode ) {
+		free_rstats_pixels(data->rpixels);
+	}
+	free(data);
+}
+
+
 static struct peakfinder_mask *create_peakfinder_mask(struct image *img,
                                                       struct radius_maps *rmps,
                                                       int min_res,
@@ -385,6 +577,31 @@ static void fill_radial_bins(float *data,
 	}
 }
 
+static void fill_radial_bins_fast(float *data, int w, int h, int n_pixels,
+				  int *pidx, int *radius, char *mask,
+				  float *rthreshold, float *lthreshold,
+				  float *roffset, float *rsigma, int *rcount)
+{
+	int i;
+
+	int curr_r;
+	float value;
+
+	for (i = 0; i < n_pixels; i++)
+	{
+		if (mask[pidx[i]] != 0)
+		{
+			curr_r = radius[i];
+			value = data[pidx[i]];
+			if (value < rthreshold[curr_r] && value > lthreshold[curr_r])
+			{
+				roffset[curr_r] += value;
+				rsigma[curr_r] += (value * value);
+				rcount[curr_r] += 1;
+			}
+		}
+	}
+}
 
 static void compute_radial_stats(float *rthreshold,
                                  float *lthreshold,
@@ -1021,9 +1238,13 @@ int peakfinder8(struct image *img, int max_n_peaks,
                 float threshold, float min_snr,
                 int min_pix_count, int max_pix_count,
                 int local_bg_radius, int min_res,
-                int max_res, int use_saturated)
+                int max_res, int use_saturated,
+                int fast_mode, struct pf8_private_data *private_data)
 {
+	struct pf8_private_data *geomdata;
 	struct radius_maps *rmaps;
+	struct radial_stats_pixels *rspixels;
+
 	struct peakfinder_mask *pfmask;
 	struct peakfinder_panel_data *pfdata;
 	struct radial_stats *rstats;
@@ -1040,18 +1261,28 @@ int peakfinder8(struct image *img, int max_n_peaks,
 
 	if ( img->detgeom == NULL) return 1;
 
-	rmaps = compute_radius_maps(img->detgeom);
-	if ( rmaps == NULL ) return 1;
+	profile_start("pf8-rmaps");
+	if ( private_data == NULL ) {
+		geomdata = prepare_peakfinder8(img->detgeom, fast_mode);
+	} else {
+		geomdata = private_data;
+	}
+	rmaps = geomdata->rmaps;
+	rspixels = geomdata->rpixels;
+	profile_end("pf8-rmaps");
+	if (geomdata == NULL) return 1;
 
+	profile_start("pf8-mask");
 	pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res);
+	profile_end("pf8-mask");
 	if ( pfmask == NULL ) {
-		free_radius_maps(rmaps);
+		if ( private_data == NULL ) free_pf8_private_data(geomdata);
 		return 1;
 	}
 
 	pfdata = allocate_panel_data(img->detgeom->n_panels);
 	if ( pfdata == NULL) {
-		free_radius_maps(rmaps);
+		if ( private_data == NULL ) free_pf8_private_data(geomdata);
 		free_peakfinder_mask(pfmask);
 		return 1;
 	}
@@ -1077,7 +1308,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
 
 	rstats = allocate_radial_stats(num_rad_bins);
 	if ( rstats == NULL ) {
-		free_radius_maps(rmaps);
+		if ( private_data == NULL ) free_pf8_private_data(geomdata);
 		free_peakfinder_mask(pfmask);
 		free_panel_data(pfdata);
 		return 1;
@@ -1087,7 +1318,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
 		rstats->rthreshold[i] = 1e9;
 		rstats->lthreshold[i] = -1e9;
 	}
-
+	profile_start("pf8-rstats");
 	for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) {
 
 		for ( i=0; i<num_rad_bins; i++ ) {
@@ -1097,18 +1328,31 @@ int peakfinder8(struct image *img, int max_n_peaks,
 		}
 
 		for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
-
-			fill_radial_bins(pfdata->panel_data[pi],
-			                 pfdata->panel_w[pi],
-			                 pfdata->panel_h[pi],
-			                 rmaps->r_maps[pi],
-			                 pfmask->masks[pi],
-			                 rstats->rthreshold,
-			                 rstats->lthreshold,
-			                 rstats->roffset,
-			                 rstats->rsigma,
-			                 rstats->rcount);
-
+			if ( fast_mode ) {
+				fill_radial_bins_fast(pfdata->panel_data[pi],
+						      pfdata->panel_w[pi],
+						      pfdata->panel_h[pi],
+						      rspixels->n_pixels[pi],
+						      rspixels->pidx[pi],
+						      rspixels->radius[pi],
+						      pfmask->masks[pi],
+						      rstats->rthreshold,
+						      rstats->lthreshold,
+						      rstats->roffset,
+						      rstats->rsigma,
+						      rstats->rcount);
+			} else {
+				fill_radial_bins(pfdata->panel_data[pi],
+						 pfdata->panel_w[pi],
+						 pfdata->panel_h[pi],
+						 rmaps->r_maps[pi],
+						 pfmask->masks[pi],
+						 rstats->rthreshold,
+						 rstats->lthreshold,
+						 rstats->roffset,
+						 rstats->rsigma,
+						 rstats->rcount);
+			}
 		}
 
 		compute_radial_stats(rstats->rthreshold,
@@ -1121,10 +1365,11 @@ int peakfinder8(struct image *img, int max_n_peaks,
 		                     threshold);
 
 	}
+	profile_end("pf8-rstats");
 
 	pkdata = allocate_peak_data(max_n_peaks);
 	if ( pkdata == NULL ) {
-		free_radius_maps(rmaps);
+		if ( private_data == NULL ) free_pf8_private_data(geomdata);
 		free_peakfinder_mask(pfmask);
 		free_panel_data(pfdata);
 		free_radial_stats(rstats);
@@ -1132,7 +1377,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
 	}
 
 	remaining_max_num_peaks = max_n_peaks;
-
+	profile_start("pf8-search");
 	for ( pi=0 ; pi<img->detgeom->n_panels ; pi++) {
 
 		int peaks_to_add;
@@ -1165,10 +1410,11 @@ int peakfinder8(struct image *img, int max_n_peaks,
 		                       NULL);
 
 		if ( ret != 0 ) {
-			free_radius_maps(rmaps);
+			if ( private_data == NULL ) free_pf8_private_data(geomdata);
 			free_peakfinder_mask(pfmask);
 			free_panel_data(pfdata);
 			free_radial_stats(rstats);
+			profile_end("pf8-search");
 			return 1;
 		}
 
@@ -1199,8 +1445,9 @@ int peakfinder8(struct image *img, int max_n_peaks,
 			                  NULL);
 		}
 	}
+	profile_end("pf8-search");
 
-	free_radius_maps(rmaps);
+	if ( private_data == NULL ) free_pf8_private_data(geomdata);
 	free_peakfinder_mask(pfmask);
 	free_panel_data(pfdata);
 	free_radial_stats(rstats);
diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h
index 28212e49..e5e86038 100644
--- a/libcrystfel/src/peakfinder8.h
+++ b/libcrystfel/src/peakfinder8.h
@@ -43,11 +43,23 @@ extern "C" {
  * The "peakfinder8" peak search algorithm, originally implemented in Cheetah
  */
 
+struct pf8_private_data
+{
+    int fast_mode;
+    struct radius_maps *rmaps;
+    struct radial_stats_pixels *rpixels;
+};
+
+struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode);
+
+void free_pf8_private_data(struct pf8_private_data *data);
+
 extern int peakfinder8(struct image *img, int max_n_peaks,
                        float threshold, float min_snr,
                        int mix_pix_count, int max_pix_count,
                        int local_bg_radius, int min_res,
-                       int max_res, int use_saturated);
+                       int max_res, int use_saturated,
+                       int fast_mode, struct pf8_private_data *private_data);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index d260b66d..57d2c92e 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -454,10 +454,11 @@ void search_peaks(struct image *image, float threshold, float min_sq_gradient,
  * the actual \ref peakfinder8 function, found in \ref peakfinder8.h.
  */
 int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
-                              float threshold, float min_snr,
-                              int min_pix_count, int max_pix_count,
-                              int local_bg_radius, int min_res,
-                              int max_res, int use_saturated)
+                             float threshold, float min_snr,
+                             int min_pix_count, int max_pix_count,
+                             int local_bg_radius, int min_res,
+                             int max_res, int use_saturated,
+                             int fast_mode, void *private_data)
 {
 	if ( image->features != NULL ) {
 		image_feature_list_free(image->features);
@@ -467,7 +468,8 @@ int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
 	return peakfinder8(image, max_n_peaks, threshold, min_snr,
 	                   min_pix_count, max_pix_count,
 	                   local_bg_radius, min_res,
-	                   max_res, use_saturated);
+	                   max_res, use_saturated,
+			   fast_mode, private_data);
 }
 
 
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 2f100db2..d6d5c1f9 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -73,7 +73,8 @@ extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
                                     float threshold, float min_snr,
                                     int mix_pix_count, int max_pix_count,
                                     int local_bg_radius, int min_res,
-                                    int max_res, int use_saturated);
+                                    int max_res, int use_saturated,
+                                    int fast_mode, void *private_data);
 
 extern int search_peaks_peakfinder9(struct image *image,
                                     float min_snr_biggest_pix,
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 95345f40..0195e0ff 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -49,6 +49,7 @@
 #include "reflist.h"
 #include "reflist-utils.h"
 #include "datatemplate.h"
+#include "datatemplate_priv.h"
 #include "detgeom.h"
 #include "libcrystfel-version.h"
 
@@ -896,7 +897,8 @@ struct image *stream_read_chunk(Stream *st, StreamFlags srf)
 			if ( have_filename && have_ev ) {
 				/* Success */
 				if ( srf & STREAM_DATA_DETGEOM ) {
-					if ( create_detgeom(image, st->dtempl_read) ) {
+					image->detgeom = create_detgeom(image, st->dtempl_read, 0);
+					if ( image->detgeom == NULL ) {
 						image_free(image);
 						return NULL;
 					}