Merge branch 'tom/imagedata'

17 files changed, 473 insertions, 361 deletions

diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 99526ff9..53510d69 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -412,33 +412,17 @@ int detector_has_clen_references(struct detector *det)
 }
 
 
-void record_image(struct image *image, int do_poisson, double background,
-                  gsl_rng *rng, double beam_radius, double nphotons)
+static void record_panel(struct panel *p, float *dp, int do_poisson,
+                         gsl_rng *rng, double ph_per_e, double background,
+			 double lambda,
+                         int *n_neg1, int *n_inf1, int *n_nan1,
+                         int *n_neg2, int *n_inf2, int *n_nan2,
+			 double *max_tt)
 {
-	int x, y;
-	double total_energy, energy_density;
-	double ph_per_e;
-	double area;
-	double max_tt = 0.0;
-	int n_inf1 = 0;
-	int n_neg1 = 0;
-	int n_nan1 = 0;
-	int n_inf2 = 0;
-	int n_neg2 = 0;
-	int n_nan2 = 0;
-
-	/* How many photons are scattered per electron? */
-	area = M_PI*pow(beam_radius, 2.0);
-	total_energy = nphotons * ph_lambda_to_en(image->lambda);
-	energy_density = total_energy / area;
-	ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0);
-	STATUS("Fluence = %8.2e photons, "
-	       "Energy density = %5.3f kJ/cm^2, "
-	       "Total energy = %5.3f microJ\n",
-	       nphotons, energy_density/1e7, total_energy*1e6);
+	int fs, ss;
 
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
+	for ( ss=0; ss>p->h; ss++ ) {
+	for ( fs=0; fs>p->w; fs++ ) {
 
 		double counts;
 		double cf;
@@ -447,28 +431,27 @@ void record_image(struct image *image, int do_poisson, double background,
 		double xs, ys, rx, ry;
 		double dsq, proj_area;
 		float dval;
-		struct panel *p;
-
-		intensity = (double)image->data[x + image->width*y];
-		if ( isinf(intensity) ) n_inf1++;
-		if ( intensity < 0.0 ) n_neg1++;
-		if ( isnan(intensity) ) n_nan1++;
+		double twotheta;
 
-		p = find_panel(image->det, x, y);
+		intensity = (double)dp[fs + p->w*ss];
+		if ( isinf(intensity) ) (*n_inf1)++;
+		if ( intensity < 0.0 ) (*n_neg1)++;
+		if ( isnan(intensity) ) (*n_nan1)++;
 
 		/* Area of one pixel */
 		pix_area = pow(1.0/p->res, 2.0);
 		Lsq = pow(p->clen, 2.0);
 
-		/* Area of pixel as seen from crystal (approximate) */
-		proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
-
 		/* Calculate distance from crystal to pixel */
-		xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
-		ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
+		xs = fs*p->fsx + ss*p->ssx;
+		ys = ss*p->fsy + ss*p->ssy;
 		rx = (xs + p->cnx) / p->res;
 		ry = (ys + p->cny) / p->res;
 		dsq = pow(rx, 2.0) + pow(ry, 2.0);
+		twotheta = atan2(sqrt(dsq), p->clen);
+
+		/* Area of pixel as seen from crystal (approximate) */
+		proj_area = pix_area * cos(twotheta);
 
 		/* Projected area of pixel divided by distance squared */
 		sa = proj_area / (dsq + Lsq);
@@ -484,36 +467,60 @@ void record_image(struct image *image, int do_poisson, double background,
 		dval = counts + poisson_noise(rng, background);
 
 		/* Convert to ADU */
-		dval *= p->adu_per_eV * ph_lambda_to_eV(image->lambda);
+		dval *= p->adu_per_eV * ph_lambda_to_eV(lambda);
 
 		/* Saturation */
 		if ( dval > p->max_adu ) dval = p->max_adu;
 
-		image->data[x + image->width*y] = dval;
+		dp[fs + p->w*ss] = dval;
 
 		/* Sanity checks */
-		if ( isinf(image->data[x+image->width*y]) ) n_inf2++;
-		if ( isnan(image->data[x+image->width*y]) ) n_nan2++;
-		if ( image->data[x+image->width*y] < 0.0 ) n_neg2++;
+		if ( isinf(dp[fs + p->w*ss]) ) n_inf2++;
+		if ( isnan(dp[fs + p->w*ss]) ) n_nan2++;
+		if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++;
 
-		if ( image->twotheta[x + image->width*y] > max_tt ) {
-			max_tt = image->twotheta[x + image->width*y];
-		}
+		if ( twotheta > *max_tt ) *max_tt = twotheta;
 
 	}
-	progress_bar(x, image->width-1, "Post-processing");
 	}
+}
 
-	STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
-	        rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
 
-	double tt_side = image->twotheta[(image->width/2)+image->width*0];
-	STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n",
-	        rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
+void record_image(struct image *image, int do_poisson, double background,
+                  gsl_rng *rng, double beam_radius, double nphotons)
+{
+	double total_energy, energy_density;
+	double ph_per_e;
+	double area;
+	double max_tt = 0.0;
+	int n_inf1 = 0;
+	int n_neg1 = 0;
+	int n_nan1 = 0;
+	int n_inf2 = 0;
+	int n_neg2 = 0;
+	int n_nan2 = 0;
+	int pn;
 
-	tt_side = image->twotheta[0+image->width*(image->height/2)];
-	STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n",
-	        rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
+	/* How many photons are scattered per electron? */
+	area = M_PI*pow(beam_radius, 2.0);
+	total_energy = nphotons * ph_lambda_to_en(image->lambda);
+	energy_density = total_energy / area;
+	ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0);
+	STATUS("Fluence = %8.2e photons, "
+	       "Energy density = %5.3f kJ/cm^2, "
+	       "Total energy = %5.3f microJ\n",
+	       nphotons, energy_density/1e7, total_energy*1e6);
+
+	for ( pn=0; pn<image->det->n_panels; pn++ ) {
+		record_panel(&image->det->panels[pn], image->dp[pn],
+		             do_poisson, rng, ph_per_e, background,
+			     image->lambda,
+		             &n_neg1, &n_inf1, &n_nan1,
+		             &n_neg2, &n_inf2, &n_nan2, &max_tt);
+	}
+
+	STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
+	        rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
 
 	STATUS("Halve the d values to get the voxel size for a synthesis.\n");
 
@@ -554,9 +561,7 @@ struct panel *find_panel(struct detector *det, double fs, double ss)
 }
 
 
-/* Like find_panel(), but uses the original panel bounds, i.e. referring to
- * what's in the HDF5 file */
-struct panel *find_orig_panel(struct detector *det, double fs, double ss)
+signed int find_orig_panel_number(struct detector *det, double fs, double ss)
 {
 	int p;
 
@@ -564,13 +569,20 @@ struct panel *find_orig_panel(struct detector *det, double fs, double ss)
 		if ( (fs >= det->panels[p].orig_min_fs)
 		  && (fs < det->panels[p].orig_max_fs+1)
 		  && (ss >= det->panels[p].orig_min_ss)
-		  && (ss < det->panels[p].orig_max_ss+1) )
-		{
-			return &det->panels[p];
-		}
+		  && (ss < det->panels[p].orig_max_ss+1) ) return p;
 	}
 
-	return NULL;
+	return -1;
+}
+
+
+/* Like find_panel(), but uses the original panel bounds, i.e. referring to
+ * what's in the HDF5 file */
+struct panel *find_orig_panel(struct detector *det, double fs, double ss)
+{
+	signed int pn = find_orig_panel_number(det, fs, ss);
+	if ( pn == -1 ) return NULL;
+	return &det->panels[pn];
 }
 
 
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index b3474b7a..5e544b75 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -205,6 +205,8 @@ extern struct panel *find_panel(struct detector *det, double fs, double ss);
 extern signed int find_panel_number(struct detector *det, double fs, double ss);
 extern struct panel *find_orig_panel(struct detector *det,
                                      double fs, double ss);
+extern signed int find_orig_panel_number(struct detector *det,
+                                         double fs, double ss);
 
 extern struct detector *get_detector_geometry(const char *filename,
                                               struct beam_params *beam);
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index e7b6eed4..b159327e 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -748,7 +748,7 @@ static struct hdf5_write_location *make_location_list(struct detector *det,
 }
 
 
-static void write_location(hid_t fh, const struct image *image,
+static void write_location(hid_t fh, struct detector *det, float *data,
                            struct hdf5_write_location *loc)
 {
 	hid_t sh, dh, ph;
@@ -784,7 +784,7 @@ static void write_location(hid_t fh, const struct image *image,
 		struct panel p;
 		int r;
 
-		p = image->det->panels[loc->panel_idxs[pi]];
+		p = det->panels[loc->panel_idxs[pi]];
 
 		f_offset[0] = p.orig_min_ss;
 		f_offset[1] = p.orig_min_fs;
@@ -810,14 +810,14 @@ static void write_location(hid_t fh, const struct image *image,
 		m_count[0] = p.max_ss - p.min_ss +1;
 		m_count[1] = p.max_fs - p.min_fs +1;
 
-		dimsm[0] = image->height;
-		dimsm[1] = image->width;
+		dimsm[0] = det->max_fs + 1;
+		dimsm[1] = det->max_ss + 1;
 		memspace = H5Screate_simple(2, dimsm, NULL);
 
 		r = H5Sselect_hyperslab(memspace, H5S_SELECT_SET,
 		                        m_offset, NULL, m_count, NULL);
 		r = H5Dwrite(dh, H5T_NATIVE_FLOAT, memspace,
-		             dh_dataspace, H5P_DEFAULT, image->data);
+		             dh_dataspace, H5P_DEFAULT, data);
 		if ( r < 0 ) {
 			ERROR("Couldn't write data\n");
 			H5Pclose(ph);
@@ -948,6 +948,34 @@ static void write_spectrum(hid_t fh, struct sample *spectrum, int spectrum_size,
 }
 
 
+static float *make_array_from_dp(const struct image *image)
+{
+	int i;
+	float *data;
+
+	data = malloc(image->width * image->height * sizeof(float));
+	if ( data == NULL ) {
+		ERROR("Failed to allocate data\n");
+		return NULL;
+	}
+
+	for ( i=0; i<image->det->n_panels; i++ ) {
+
+		int fs, ss;
+		struct panel *p = &image->det->panels[i];
+
+		for ( ss=0; ss<p->h; ss++ ) {
+		for ( fs=0; fs<p->w; fs++ ) {
+			int idx = p->min_fs+fs + image->width*(p->min_ss+ss);
+			data[idx] = image->dp[i][fs + p->w*ss];
+		}
+		}
+	}
+
+	return data;
+}
+
+
 int hdf5_write_image(const char *filename, const struct image *image,
                      char *element)
 {
@@ -957,12 +985,16 @@ int hdf5_write_image(const char *filename, const struct image *image,
 	struct hdf5_write_location *locations;
 	int num_locations;
 	const char *ph_en_loc;
+	float *data;
 
 	if ( image->det == NULL ) {
 		ERROR("Geometry not available\n");
 		return 1;
 	}
 
+	data = make_array_from_dp(image);
+	if ( data == NULL ) return 1;
+
 	fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
 	if ( fh < 0 ) {
 		ERROR("Couldn't create file: %s\n", filename);
@@ -979,7 +1011,7 @@ int hdf5_write_image(const char *filename, const struct image *image,
 	                               &num_locations);
 
 	for ( li=0; li<num_locations; li++ ) {
-		write_location(fh, image, &locations[li]);
+		write_location(fh, image->det, data, &locations[li]);
 	}
 
 	if ( image->beam == NULL || image->beam->photon_energy_from == NULL ) {
@@ -1001,7 +1033,7 @@ int hdf5_write_image(const char *filename, const struct image *image,
 	for ( li=0; li<num_locations; li ++ ) {
 		free(locations[li].panel_idxs);
 	}
-
+	free(data);
 	free(locations);
 	return 0;
 }
@@ -1062,18 +1094,28 @@ static void debodge_saturation(struct hdfile *f, struct image *image)
 
 	for ( i=0; i<size[0]; i++ ) {
 
-		unsigned int x, y;
+		unsigned int fs, ss;
 		float val;
+		struct panel *p;
+		signed int pn;
 
-		x = buf[3*i+0];
-		y = buf[3*i+1];
+		fs = buf[3*i+0];
+		ss = buf[3*i+1];
 		val = buf[3*i+2];
 
-		image->data[x+image->width*y] = val / 5.0;
-		image->data[x+1+image->width*y] = val / 5.0;
-		image->data[x-1+image->width*y] = val / 5.0;
-		image->data[x+image->width*(y+1)] = val / 5.0;
-		image->data[x+image->width*(y-1)] = val / 5.0;
+		/* 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;
 
 	}
 
@@ -1083,7 +1125,8 @@ static void debodge_saturation(struct hdfile *f, struct image *image)
 }
 
 
-static int unpack_panels(struct image *image, struct detector *det)
+static int unpack_panels(struct image *image, struct detector *det,
+                         float *data, uint16_t *flags)
 {
 	int pi;
 
@@ -1118,7 +1161,7 @@ static int unpack_panels(struct image *image, struct detector *det)
 			css = ss+p->min_ss;
 			idx = cfs + css*image->width;
 
-			image->dp[pi][fs+p->w*ss] = image->data[idx];
+			image->dp[pi][fs+p->w*ss] = data[idx];
 
 			if ( p->no_index ) bad = 1;
 
@@ -1126,18 +1169,18 @@ static int unpack_panels(struct image *image, struct detector *det)
 				bad = 1;
 			}
 
-			if ( image->flags != NULL ) {
+			if ( flags != NULL ) {
 
-				int flags;
+				int f;
 
-				flags = image->flags[idx];
+				f = flags[idx];
 
 				/* Bad if it's missing any of the "good" bits */
-				if ( (flags & image->det->mask_good)
+				if ( (f & image->det->mask_good)
 				       != image->det->mask_good ) bad = 1;
 
 				/* Bad if it has any of the "bad" bits. */
-				if ( flags & image->det->mask_bad ) bad = 1;
+				if ( f & image->det->mask_bad ) bad = 1;
 
 			}
 			image->bad[pi][fs+p->w*ss] = bad;
@@ -1403,17 +1446,15 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
 		free(buf);
 		return 1;
 	}
-	image->data = buf;
 
 	if ( image->det != NULL ) {
 		ERROR("WARNING: hdf5_read() called with geometry structure.\n");
 	}
 	image->det = simple_geometry(image);
 
+	unpack_panels(image, image->det, buf, NULL);
 	if ( satcorr ) debodge_saturation(f, image);
 
-	unpack_panels(image, image->det);
-
 	if ( image->beam != NULL ) {
 
 		fill_in_beam_parameters(image->beam, f, NULL, image);
@@ -1432,12 +1473,12 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
 }
 
 
-static void load_mask(struct hdfile *f, struct event *ev, char *mask,
-                      const char *mask_file,
-                      const char *pname, struct image *image,
-                      size_t p_w, size_t sum_p_h,
-                      hsize_t *f_offset, hsize_t *f_count,
-                      hsize_t *m_offset, hsize_t *m_count)
+static int load_mask(struct hdfile *f, struct event *ev, char *mask,
+                     const char *mask_file,
+                     const char *pname, struct image *image,
+                     size_t p_w, size_t sum_p_h, uint16_t *flags,
+                     hsize_t *f_offset, hsize_t *f_count,
+                     hsize_t *m_offset, hsize_t *m_count)
 {
 	hid_t mask_dataspace, mask_dh;
 	int exists;
@@ -1450,7 +1491,7 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask,
 		fh = H5Fopen(mask_file, H5F_ACC_RDONLY, H5P_DEFAULT);
 		if ( fh < 0 ) {
 			ERROR("Couldn't open mask file '%s'\n", mask_file);
-			return;
+			return 1;
 		}
 	} else {
 		fh = f->fh;
@@ -1491,7 +1532,7 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask,
 	}
 
 	r = H5Dread(mask_dh, H5T_NATIVE_UINT16, memspace,
-	            mask_dataspace, H5P_DEFAULT, image->flags);
+	            mask_dataspace, H5P_DEFAULT, flags);
 	if ( r < 0 ) {
 		ERROR("Couldn't read flags for panel %s\n", pname);
 		goto err;
@@ -1501,14 +1542,12 @@ static void load_mask(struct hdfile *f, struct event *ev, char *mask,
 	H5Dclose(mask_dh);
 	if ( ev != NULL ) free(mask);
 
-	return;
+	return 0;
 
 err:
 	if ( mask_file != NULL ) H5Fclose(fh);
 	if ( ev != NULL ) free(mask);
-	free(image->flags);
-	image->flags = NULL;
-	return;
+	return 1;
 }
 
 
@@ -1517,6 +1556,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 {
 	herr_t r;
 	float *buf;
+	uint16_t *flags;
 	int sum_p_h;
 	int p_w;
 	int pi;
@@ -1548,8 +1588,8 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 	image->width = p_w;
 	image->height = sum_p_h;
 
-	image->flags = calloc(p_w*sum_p_h,sizeof(uint16_t));
-	if ( image->flags == NULL ) {
+	flags = calloc(p_w*sum_p_h,sizeof(uint16_t));
+	if ( flags == NULL ) {
 		ERROR("Failed to allocate memory for flags\n");
 		return 1;
 	}
@@ -1695,9 +1735,11 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 		H5Sclose(memspace);
 
 		if ( p->mask != NULL ) {
-			load_mask(f, ev, p->mask, p->mask_file, p->name,
-			          image, p_w, sum_p_h,
-			          f_offset, f_count, m_offset, m_count);
+			if ( load_mask(f, ev, p->mask, p->mask_file, p->name,
+			               image, p_w, sum_p_h, flags,
+			               f_offset, f_count, m_offset, m_count) ) {
+				ERROR("Error loading bad pixel mask!\n");
+			}
 		}
 
 		free(f_offset);
@@ -1705,13 +1747,10 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 
 	}
 
-	image->data = buf;
-
-	if ( satcorr ) debodge_saturation(f, image);
-
 	fill_in_values(image->det, f, ev);
 
-	unpack_panels(image, image->det);
+	unpack_panels(image, image->det, buf, flags);
+	if ( satcorr ) debodge_saturation(f, image);
 
 	if ( image->beam != NULL ) {
 
@@ -1727,6 +1766,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
 
 	}
 
+	free(buf);
+	free(flags);
+
 	return 0;
 }
 
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 10e83905..52614d06 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -113,10 +113,6 @@ struct beam_params
  * <programlisting>
  * struct image
  * {
- *    float                   *data;
- *    uint16_t                *flags;
- *    double                  *twotheta;
- *
  *    Crystal                 **crystals;
  *    int                     n_crystals;
  *    IndexingMethod          indexed_by;
@@ -163,11 +159,6 @@ struct image;
 
 struct image {
 
-	/* The following three fields will be going away in the future */
-	float                   *data;
-	uint16_t                *flags;
-	double                  *twotheta;
-
 	float                   **dp;    /* Data in panel */
 	int                     **bad;   /* Bad pixels by panel */
 
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 9e3c5e77..69aad4c0 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -242,7 +242,7 @@ static int integrate_peak(struct image *image, int cfs, int css,
 	double bg_tot_sq = 0.0;
 	double var;
 	double aduph;
-	int p_cfs, p_css, p_w, p_h;
+	int p_cfs, p_css;
 	signed int pn;
 
 	pn = find_panel_number(image->det, cfs, css);
@@ -254,8 +254,6 @@ static int integrate_peak(struct image *image, int cfs, int css,
 	/* Determine regions where there is expected to be a peak */
 	p_cfs = cfs - p->min_fs;
 	p_css = css - p->min_ss;  /* Panel-relative coordinates */
-	p_w = p->max_fs - p->min_fs + 1;
-	p_h = p->max_ss - p->min_ss + 1;
 
 	aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda);
 
@@ -264,8 +262,8 @@ static int integrate_peak(struct image *image, int cfs, int css,
 	out_lim_sq = pow(ir_out, 2.0);
 
 	/* Estimate the background */
-	for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) {
 	for ( dss=-ir_out; dss<=+ir_out; dss++ ) {
+	for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) {
 
 		double val;
 		int idx;
@@ -275,7 +273,7 @@ static int integrate_peak(struct image *image, int cfs, int css,
 		if ( dfs*dfs + dss*dss < mid_lim_sq ) continue;
 
 		/* Strayed off one panel? */
-		if ( (p_cfs+dfs >= p_w) || (p_css+dss >= p_h)
+		if ( (p_cfs+dfs >= p->w) || (p_css+dss >= p->h)
 		  || (p_cfs+dfs < 0 ) || (p_css+dss < 0) ) return 4;
 
 		/* Wandered into a bad region? */
@@ -283,8 +281,8 @@ static int integrate_peak(struct image *image, int cfs, int css,
 			return 14;
 		}
 
-		idx = dfs+cfs+image->width*(dss+css);
-		val = image->data[idx];
+		idx = dfs+p_cfs+p->w*(dss+p_css);
+		val = image->dp[pn][idx];
 
 		/* Check if peak contains saturation in bg region */
 		if ( (saturated != NULL) && (val > p->max_adu) ) *saturated = 1;
@@ -304,8 +302,8 @@ static int integrate_peak(struct image *image, int cfs, int css,
 	pk_total = 0.0;
 	pk_counts = 0;
 	fsct = 0.0;  ssct = 0.0;
-	for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) {
 	for ( dss=-ir_inn; dss<=+ir_inn; dss++ ) {
+	for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) {
 
 		double val;
 		int idx;
@@ -314,7 +312,7 @@ static int integrate_peak(struct image *image, int cfs, int css,
 		if ( dfs*dfs + dss*dss > lim_sq ) continue;
 
 		/* Strayed off one panel? */
-		if ( (p_cfs+dfs >= p_w) || (p_css+dss >= p_h)
+		if ( (p_cfs+dfs >= p->w) || (p_css+dss >= p->h)
 		  || (p_cfs+dfs < 0 ) || (p_css+dss < 0) ) return 8;
 
 		/* Wandered into a bad region? */
@@ -322,8 +320,8 @@ static int integrate_peak(struct image *image, int cfs, int css,
 			return 15;
 		}
 
-		idx = dfs+cfs+image->width*(dss+css);
-		val = image->data[idx];
+		idx = dfs+p_cfs+p->w*(dss+p_css);
+		val = image->dp[pn][idx];
 
 		/* Check if peak contains saturation */
 		if ( (saturated != NULL) && (val > p->max_adu) ) *saturated = 1;
@@ -356,13 +354,13 @@ static int integrate_peak(struct image *image, int cfs, int css,
 
 
 static void search_peaks_in_panel(struct image *image, float threshold,
-                                  float min_gradient, float min_snr,
-                                  struct panel *p,
+                                  float min_gradient, float min_snr, int pn,
                                   double ir_inn, double ir_mid, double ir_out,
                                   int use_saturated)
 {
 	int fs, ss, stride;
 	float *data;
+	struct panel *p;
 	double d;
 	int idx;
 	double f_fs = 0.0;
@@ -378,11 +376,12 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 	int nacc = 0;
 	int ncull;
 
-	data = image->data;
-	stride = image->width;
+	p = &image->det->panels[pn];
+	data = image->dp[pn];
+	stride = p->w;
 
-	for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) {
-	for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) {
+	for ( ss=0; ss<=p->h; ss++ ) {
+	for ( fs=0; fs<=p->w; fs++ ) {
 
 		double dx1, dx2, dy1, dy2;
 		double dxs, dys;
@@ -425,12 +424,12 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 			max = data[mask_fs+stride*mask_ss];
 			did_something = 0;
 
-			for ( s_ss=biggest(mask_ss-ir_inn, p->min_ss);
-			      s_ss<=smallest(mask_ss+ir_inn, p->max_ss);
+			for ( s_ss=biggest(mask_ss-ir_inn, 0);
+			      s_ss<=smallest(mask_ss+ir_inn, p->h-1);
 			      s_ss++ )
 			{
-			for ( s_fs=biggest(mask_fs-ir_inn, p->min_fs);
-			      s_fs<=smallest(mask_fs+ir_inn, p->max_fs);
+			for ( s_fs=biggest(mask_fs-ir_inn, 0);
+			      s_fs<=smallest(mask_fs+ir_inn, p->w-1);
 			      s_fs++ )
 			{
 
@@ -459,13 +458,13 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		}
 
 		/* Should be enforced by bounds used above.  Muppet check. */
-		assert(mask_fs <= p->max_fs);
-		assert(mask_ss <= p->max_ss);
-		assert(mask_fs >= p->min_fs);
-		assert(mask_ss >= p->min_ss);
+		assert(mask_fs <= p->w);
+		assert(mask_ss <= p->h);
+		assert(mask_fs >= 0);
+		assert(mask_ss >= 0);
 
 		/* Centroid peak and get better coordinates. */
-		r = integrate_peak(image, mask_fs, mask_ss,
+		r = integrate_peak(image, mask_fs+p->min_fs, mask_ss+p->min_ss,
 		                   &f_fs, &f_ss, &intensity, &sigma,
 		                   ir_inn, ir_mid, ir_out, &saturated);
 
@@ -504,8 +503,8 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		}
 
 		/* Add using "better" coordinates */
-		image_add_feature(image->features, f_fs, f_ss, image, intensity,
-		                  NULL);
+		image_add_feature(image->features, f_fs, f_ss,
+		                  image, intensity, NULL);
 		nacc++;
 
 	}
@@ -551,11 +550,10 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
 
 	for ( i=0; i<image->det->n_panels; i++ ) {
 
-		struct panel *p = &image->det->panels[i];
+		if ( image->det->panels[i].no_index ) continue;
 
-		if ( p->no_index ) continue;
 		search_peaks_in_panel(image, threshold, min_gradient,
-		                      min_snr, p, ir_inn, ir_mid, ir_out,
+		                      min_snr, i, ir_inn, ir_mid, ir_out,
 		                      use_saturated);
 
 	}
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index bdfdb13b..1f937bf8 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -3,11 +3,11 @@
  *
  * Calculate diffraction patterns by Fourier methods (GPU version)
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2014 Thomas White <taw@physics.org>
+ *   2009-2015 Thomas White <taw@physics.org>
  *   2013      Alexandra Tolstikova
  *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
  *
@@ -253,20 +253,17 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
 			return 1;
 		}
 
-		for ( fs=0; fs<pan_width; fs++ ) {
 		for ( ss=0; ss<pan_height; ss++ ) {
+		for ( fs=0; fs<pan_width; fs++ ) {
 
 			float val;
-			int tfs, tss;
 
 			val = diff_ptr[fs + pan_width*ss];
 			if ( isinf(val) ) (*n_inf)++;
 			if ( val < 0.0 ) (*n_neg)++;
 			if ( isnan(val) ) (*n_nan)++;
 
-			tfs = p->min_fs + fs;
-			tss = p->min_ss + ss;
-			image->data[tfs + image->width*tss] += val;
+			image->dp[i][fs + pan_width*ss] += val;
 
 		}
 		}
@@ -294,7 +291,6 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	int n_inf = 0;
 	int n_neg = 0;
 	int n_nan = 0;
-	int fs, ss;
 	int i;
 
 	if ( gctx == NULL ) {
@@ -326,11 +322,19 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return 1;
 
 	/* Allocate memory for the result */
-	image->data = calloc(image->width * image->height, sizeof(float));
-	if ( image->data == NULL ) {
+	image->dp = malloc(image->det->n_panels * sizeof(float *));
+	if ( image->dp == NULL ) {
 		ERROR("Couldn't allocate memory for result.\n");
 		return 1;
 	}
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		struct panel *p = &image->det->panels[i];
+		image->dp[i] = calloc(p->w * p->h, sizeof(float));
+		if ( image->dp[i] == NULL ) {
+			ERROR("Couldn't allocate memory for panel %i\n", i);
+			return 1;
+		}
+	}
 
 	double tot = 0.0;
 	for ( i=0; i<image->nsamples; i++ ) {
@@ -356,18 +360,6 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 		      n_neg, n_inf, n_nan);
 	}
 
-	/* Calculate "2theta" values for detector geometry */
-	image->twotheta = calloc(image->width * image->height, sizeof(double));
-	for ( fs=0; fs<image->width; fs++ ) {
-	for ( ss=0; ss<image->height; ss++ ) {
-
-		double twotheta;
-		get_q(image, fs, ss, &twotheta, 1.0/image->lambda);
-		image->twotheta[fs + image->width*ss] = twotheta;
-
-	}
-	}
-
 	return 0;
 }
 
diff --git a/src/diffraction.c b/src/diffraction.c
index 5a936809..93674f52 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -361,31 +361,30 @@ static double molecule_factor(const double *intensities, const double *phases,
 }
 
 
-
-static void diffraction_at_k(struct image *image, const double *intensities,
-                             const double *phases, const unsigned char *flags,
-                             UnitCell *cell, GradientMethod m,
-                             const SymOpList *sym, double k,
-                             double ax, double ay, double az,
-                             double bx, double by, double bz,
-                             double cx, double cy, double cz,
-                             double *lut_a, double *lut_b, double *lut_c,
-                             double weight)
+static void diffraction_panel(struct image *image, const double *intensities,
+                              const double *phases, const unsigned char *flags,
+                              UnitCell *cell, GradientMethod m,
+                              const SymOpList *sym, double k,
+                              double ax, double ay, double az,
+                              double bx, double by, double bz,
+                              double cx, double cy, double cz,
+                              double *lut_a, double *lut_b, double *lut_c,
+                              int pn, double weight)
 {
-	unsigned int fs, ss;
+	int fs, ss;
 	const int nxs = 4;
 	const int nys = 4;
+	struct panel *p = &image->det->panels[pn];
 
 	weight /= nxs*nys;
 
-	for ( fs=0; fs<image->width; fs++ ) {
-	for ( ss=0; ss<image->height; ss++ ) {
+	for ( ss=0; ss<p->h; ss++ ) {
+	for ( fs=0; fs<p->w; fs++ ) {
 
 		int idx;
 		double f_lattice, I_lattice;
 		double I_molecule;
 		struct rvec q;
-		double twotheta;
 		int xs, ys;
 		float xo, yo;
 
@@ -395,7 +394,7 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 			xo = (1.0/nxs) * xs;
 			yo = (1.0/nys) * ys;
 
-			q = get_q(image, fs+xo, ss+yo, &twotheta, k);
+			q = get_q_for_panel(p, fs+xo, ss+yo, NULL, k);
 
 			f_lattice = lattice_factor(q, ax, ay, az,
 					           bx, by, bz,
@@ -411,14 +410,33 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 
 			I_lattice = pow(f_lattice, 2.0);
 
-			idx = fs + image->width*ss;
-			image->data[idx] += I_lattice * I_molecule * weight;
-			image->twotheta[idx] = twotheta;
+			idx = fs + p->w*ss;
+			image->dp[pn][idx] += I_lattice * I_molecule * weight;
 
 		}
 		}
 	}
-	progress_bar(fs, image->width-1, "Calculating diffraction");
+	progress_bar(ss, p->h-1, "Calculating diffraction");
+	}
+}
+
+
+static void diffraction_at_k(struct image *image, const double *intensities,
+                             const double *phases, const unsigned char *flags,
+                             UnitCell *cell, GradientMethod m,
+                             const SymOpList *sym, double k,
+                             double ax, double ay, double az,
+                             double bx, double by, double bz,
+                             double cx, double cy, double cz,
+                             double *lut_a, double *lut_b, double *lut_c,
+                             double weight)
+{
+	int i;
+
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		diffraction_panel(image, intensities, phases, flags, cell, m,
+		                  sym, k, ax, ay, az, bx, by, bz, cx, cy, cz,
+		                  lut_a, lut_b, lut_c, i, weight);
 	}
 }
 
@@ -710,12 +728,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
-	/* Allocate (and zero) the "diffraction array" */
-	image->data = calloc(image->width * image->height, sizeof(float));
-
-	/* Needed later for Lorentz calculation */
-	image->twotheta = malloc(image->width * image->height * sizeof(double));
-
 	lut_a = get_sinc_lut(na, no_fringes);
 	lut_b = get_sinc_lut(nb, no_fringes);
 	lut_c = get_sinc_lut(nc, no_fringes);
diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 004c92ae..e260e4a0 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -94,8 +94,6 @@ static gint displaywindow_closed(GtkWidget *window, DisplayWindow *dw)
 	if ( dw->image != NULL ) {
 
 		free(dw->image->filename);
-		free(dw->image->data);
-		free(dw->image->flags);
 		free(dw->image);
 	}
 
@@ -207,7 +205,7 @@ static int render_adsc_uint16(DisplayWindow *dw, const char *filename)
 
 		fs = dfs;
 		ss = dss;
-		val = image->data[fs + image->width * ss];
+		val = 0.0;//image->data[fs + image->width * ss]; FIXME!
 		if ( val < 0 ) {
 			out = 0;
 		} else if ( val > 65535 ) {
@@ -925,8 +923,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event)
 
 	if ( dw->not_ready_yet ) return 0;
 
-	float *old_data = dw->image->data;
-	uint16_t *old_flags = dw->image->flags;
 	float **old_dp = dw->image->dp;
 	int **old_bad = dw->image->bad;
 
@@ -934,8 +930,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event)
                           dw->ev_list->events[new_event], 0);
 	if ( fail ) {
 		ERROR("Couldn't load image");
-		dw->image->data = old_data;
-		dw->image->flags = old_flags;
 		dw->image->dp = old_dp;
 		dw->image->bad = old_bad;
 		return 1;
@@ -956,8 +950,6 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event)
 		free(old_dp[i]);
 		free(old_bad[i]);
 	}
-	free(old_data);
-	free(old_flags);
 	free(old_dp);
 	free(old_bad);
 	return 0;
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
index 0ac59ddf..4144bba2 100644
--- a/src/geoptimiser.c
+++ b/src/geoptimiser.c
@@ -2193,51 +2193,6 @@ struct rectangle
 };
 
 
-static int unpack_slab(struct image *image)
-{
-	struct detector *det = image->det;
-	int pi;
-
-	image->dp = malloc(det->n_panels * sizeof(float *));
-	image->bad = malloc(det->n_panels * sizeof(int *));
-	if ( (image->dp == NULL) || (image->bad == NULL) ) {
-		ERROR("Failed to allocate panels.\n");
-		return 1;
-	}
-
-	for ( pi=0; pi<det->n_panels; pi++ ) {
-
-		struct panel *p;
-		int fs, ss;
-
-		p = &det->panels[pi];
-		image->dp[pi] = malloc(p->w*p->h*sizeof(float));
-		image->bad[pi] = calloc(p->w*p->h, sizeof(int));
-		if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) ) {
-			ERROR("Failed to allocate panel\n");
-			return 1;
-		}
-
-		for ( ss=0; ss<p->h; ss++ ) {
-			for ( fs=0; fs<p->w; fs++ ) {
-
-				int idx;
-				int cfs, css;
-
-				cfs = fs+p->min_fs;
-				css = ss+p->min_ss;
-				idx = cfs + css*image->width;
-
-				image->dp[pi][fs+p->w*ss] = image->data[idx];
-				image->bad[pi][fs+p->w*ss] = 0;
-			}
-		}
-	}
-
-	return 0;
-}
-
-
 static int draw_detector(cairo_surface_t *surf, struct image *image,
                          struct rectangle rect)
 {
@@ -2249,7 +2204,6 @@ static int draw_detector(cairo_surface_t *surf, struct image *image,
 
 	cr = cairo_create(surf);
 
-	unpack_slab(image);
 	pixbufs = render_panels(image, 1, SCALE_GEOPTIMISER, 1, &n_pixbufs);
 
 	/* Blank grey background */
@@ -2304,25 +2258,51 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet,
 	cairo_status_t r;
 	cairo_surface_t *surf;
 
-	im.data = malloc((edet->width)*(edet->height)*sizeof(float));
-	if ( im.data == NULL ) {
-		ERROR("Failed to allocate memory to save data.\n");
-		return 1;
-	}
 	im.det = edet->det;
 	im.width = edet->width;
 	im.height = edet->height;
-	im.flags = NULL;
+	im.dp = malloc(edet->det->n_panels*sizeof(float *));
+	if ( im.dp == NULL ) {
+		ERROR("Failed to allocate data\n");
+		return 1;
+	}
+	for ( i=0; i<edet->det->n_panels; i++ ) {
 
-	for ( i=0; i<(edet->width)*(edet->height); i++) {
-		if ( data[i] == -10000.0) {
-			im.data[i] = 0.0;
-		} else if ( data[i] > 1.0) {
-			im.data[i] = 1.0;
-		} else {
-			im.data[i] = (float)data[i];
+		int fs, ss;
+		struct panel *p;
+
+		p = &edet->det->panels[i];
+
+		im.dp[i] = calloc(p->w * p->h, sizeof(float));
+		if ( im.dp[i] == NULL ) {
+			ERROR("Failed to allocate data\n");
+			return 1;
+		}
+
+		for ( ss=0; ss<p->h; ss++ ) {
+		for ( fs=0; fs<p->w; fs++ ) {
+
+			int idx;
+			int cfs, css;
+			float val;
+
+			cfs = fs+p->min_fs;
+			css = ss+p->min_ss;
+			idx = cfs + css*edet->width;
+
+			if ( data[idx] == -10000.0) {
+				val = 0.0;
+			} else if ( data[idx] > 1.0) {
+				val = 1.0;
+			} else {
+				val = (float)data[idx];
+			}
+			val *= 10.0; /* render_panels sets this as max */
+
+			im.dp[i][fs+p->w*ss] = val;
+
+		}
 		}
-		im.data[i] *= 10.0; /* render_panels sets this as max */
 	}
 
 	get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x,
@@ -2354,13 +2334,13 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet,
 	cairo_fill(cr);
 	cairo_destroy(cr);
 
-	r = cairo_surface_write_to_png(surf, filename);
-	if (r != CAIRO_STATUS_SUCCESS) {
-		free(im.data);
-		return 1;
+	for ( i=0; i<edet->det->n_panels; i++ ) {
+		free(im.dp[i]);
 	}
+	free(im.dp);
 
-	free(im.data);
+	r = cairo_surface_write_to_png(surf, filename);
+	if ( r != CAIRO_STATUS_SUCCESS ) return 1;
 
 	return 0;
 }
diff --git a/src/hdfsee-render.c b/src/hdfsee-render.c
index 71e2e1e8..4c6c3114 100644
--- a/src/hdfsee-render.c
+++ b/src/hdfsee-render.c
@@ -283,8 +283,9 @@ int render_tiff_fp(struct image *image, const char *filename)
 
 	line = _TIFFmalloc(TIFFScanlineSize(th));
 	for ( y=0; y<image->height; y++ ) {
-		memcpy(line, &image->data[(image->height-1-y)*image->width],
-		       image->width*4);
+		//memcpy(line, &image->data[(image->height-1-y)*image->width],
+		//       image->width*4);
+		//       FIXME!
 		TIFFWriteScanline(th, line, y, 0);
 	}
 	_TIFFfree(line);
@@ -325,7 +326,8 @@ int render_tiff_int16(struct image *image, const char *filename, double boost)
 	for ( y=0; y<image->height; y++ ) {
 	for ( x=0;x<image->width; x++ ) {
 		double val;
-		val = image->data[x+image->height*y];
+		//val = image->data[x+image->height*y];
+		val = 0.0;  // FIXME!
 		if ( val > max ) max = val;
 	}
 	}
@@ -336,7 +338,8 @@ int render_tiff_int16(struct image *image, const char *filename, double boost)
 
 			double val;
 
-			val = image->data[x+(image->height-1-y)*image->width];
+			//val = image->data[x+(image->height-1-y)*image->width];
+			val = 0.0;  // FIXME!
 			val *= ((double)boost/max);
 
 			/* Clamp to 16-bit range,
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 867183da..16e9cb88 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -170,6 +170,18 @@ static void calculate_partials(Crystal *cr,
 }
 
 
+static signed int panel_number(struct detector *det, struct panel *p)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+		if ( &det->panels[i] == p ) return i;
+	}
+
+	return -1;
+}
+
+
 static void draw_and_write_image(struct image *image, RefList *reflections,
                                  gsl_rng *rng, double background)
 {
@@ -177,7 +189,26 @@ static void draw_and_write_image(struct image *image, RefList *reflections,
 	RefListIterator *iter;
 	int i;
 
-	image->data = calloc(image->width*image->height, sizeof(float));
+	image->dp = malloc(image->det->n_panels*sizeof(float *));
+	if ( image->dp == NULL ) {
+		ERROR("Failed to allocate data\n");
+		return;
+	}
+	for ( i=0; i<image->det->n_panels; i++ ) {
+
+		int j;
+		struct panel *p = &image->det->panels[i];
+
+		image->dp[i] = calloc(p->w * p->h, sizeof(float));
+		if ( image->dp[i] == NULL ) {
+			ERROR("Failed to allocate data\n");
+			return;
+		}
+		for ( j=0; j<p->w*p->h; j++ ) {
+			image->dp[i][j] = poisson_noise(rng, background);
+		}
+
+	}
 
 	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
@@ -186,26 +217,32 @@ static void draw_and_write_image(struct image *image, RefList *reflections,
 		double Ip;
 		double dfs, dss;
 		int fs, ss;
+		struct panel *p;
+		signed int pn;
 
 		Ip = get_intensity(refl);
 
 		get_detector_pos(refl, &dfs, &dss);
-		fs = nearbyint(dfs);
-		ss = nearbyint(dss);
+		p = get_panel(refl);
+		pn = panel_number(image->det, p);  /* Yeurgh */
+		assert(pn != -1);
+
+		fs = nearbyint(dfs) - p->min_fs;
+		ss = nearbyint(dss) - p->min_ss;
 		assert(fs >= 0);
 		assert(ss >= 0);
-		assert(fs < image->width);
-		assert(ss < image->height);
-		image->data[fs + image->width*ss] += Ip;
-
-	}
+		assert(fs < p->w);
+		assert(ss < p->h);
+		image->dp[pn][fs + p->w*ss] += Ip;
 
-	for ( i=0; i<image->width*image->height; i++ ) {
-		image->data[i] += poisson_noise(rng, background);
 	}
 
 	hdf5_write_image(image->filename, image, NULL);
-	free(image->data);
+
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		free(image->dp[i]);
+	}
+	free(image->dp);
 }
 
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index bd97ead7..f724057e 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -249,8 +249,7 @@ int main(int argc, char *argv[])
 	int c;
 	struct image image;
 	struct gpu_context *gctx = NULL;
-	struct image *powder;
-	float *powder_data;
+	struct image powder;
 	char *intfile = NULL;
 	double *intensities;
 	char *rval;
@@ -293,6 +292,7 @@ int main(int argc, char *argv[])
 	double bandwidth = 0.01;
 	double photon_energy = 9000.0;
 	struct beam_params beam;
+	int i;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -667,7 +667,21 @@ int main(int argc, char *argv[])
 	/* Initialise stuff */
 	image.filename = NULL;
 	image.features = NULL;
-	image.flags = NULL;
+	image.bad = NULL;
+	image.dp = malloc(image.det->n_panels*sizeof(float *));
+	if ( image.dp == NULL ) {
+		ERROR("Failed to allocate data\n");
+		return 1;
+	}
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		image.dp[i] = calloc(image.det->panels[i].w *
+		                     image.det->panels[i].h,
+		                     sizeof(float));
+		if ( image.dp[i] == NULL ) {
+			ERROR("Failed to allocate data\n");
+			return 1;
+		}
+	}
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 	if ( config_random ) {
@@ -679,12 +693,23 @@ int main(int argc, char *argv[])
 		gsl_rng_set(rng, seed);
 	}
 
-	powder = calloc(1, sizeof(struct image));
-	powder->width = image.width;
-	powder->height = image.height;
-	powder->det = image.det;
-	powder_data = calloc(image.width*image.height, sizeof(float));
-	powder->data = powder_data;
+	powder.width = image.width;
+	powder.height = image.height;
+	powder.det = image.det;
+	powder.dp = malloc(image.det->n_panels*sizeof(float *));
+	if ( powder.dp == NULL ) {
+		ERROR("Failed to allocate powder data\n");
+		return 1;
+	}
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		powder.dp[i] = calloc(image.det->panels[i].w *
+		                      image.det->panels[i].h,
+		                      sizeof(float));
+		if ( powder.dp[i] == NULL ) {
+			ERROR("Failed to allocate powder data\n");
+			return 1;
+		}
+	}
 
 	/* Splurge a few useful numbers */
 	STATUS("Simulation parameters:\n");
@@ -732,6 +757,7 @@ int main(int argc, char *argv[])
 		int na, nb, nc;
 		double a, b, c, d;
 		UnitCell *cell;
+		int err = 0;
 
 		if ( random_size ) {
 
@@ -825,10 +851,6 @@ int main(int argc, char *argv[])
 
 		}
 
-		/* Ensure no residual information */
-		image.data = NULL;
-		image.twotheta = NULL;
-
 		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
 		STATUS("Particle size = %i x %i x %i"
 		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
@@ -836,8 +858,6 @@ int main(int argc, char *argv[])
 
 		if ( config_gpu ) {
 
-			int err;
-
 			if ( gctx == NULL ) {
 				gctx = setup_gpu(config_nosfac,
 				                 intensities, flags, sym_str,
@@ -845,13 +865,12 @@ int main(int argc, char *argv[])
 			}
 			err = get_diffraction_gpu(gctx, &image, na, nb, nc,
 			                          cell, no_fringes);
-			if ( err ) image.data = NULL;
 
 		} else {
 			get_diffraction(&image, na, nb, nc, intensities, phases,
 			                flags, cell, grad, sym, no_fringes);
 		}
-		if ( image.data == NULL ) {
+		if ( err ) {
 			ERROR("Diffraction calculation failed.\n");
 			done = 1;
 			goto skip;
@@ -862,18 +881,23 @@ int main(int argc, char *argv[])
 
 		if ( powder_fn != NULL ) {
 
-			int x, y, w;
+			int pn;
 
-			w = image.width;
+			for ( pn=0; pn<image.det->n_panels; pn++ ) {
+
+				size_t w, i;
+
+				w = image.det->panels[pn].w
+				  * image.det->panels[pn].h;
+
+				for ( i=0; i<w; i++ ) {
+					powder.dp[pn][i] += (double)image.dp[pn][i];
+				}
 
-			for ( x=0; x<image.width; x++ ) {
-			for ( y=0; y<image.height; y++ ) {
-				powder->data[x+w*y] += (double)image.data[x+w*y];
-			}
 			}
 
 			if ( !(ndone % 10) ) {
-				hdf5_write_image(powder_fn, powder, NULL);
+				hdf5_write_image(powder_fn, &powder, NULL);
 			}
 		}
 
@@ -895,10 +919,6 @@ int main(int argc, char *argv[])
 
 		}
 
-		/* Clean up */
-		free(image.data);
-		free(image.twotheta);
-
 		cell_free(cell);
 
 skip:
@@ -909,18 +929,18 @@ skip:
 	} while ( !done );
 
 	if ( powder_fn != NULL ) {
-		hdf5_write_image(powder_fn, powder, NULL);
+		hdf5_write_image(powder_fn, &powder, NULL);
 	}
 
 	if ( gctx != NULL ) {
 		cleanup_gpu(gctx);
 	}
 
-	free(image.det->panels);
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		free(powder.dp[i]);
+	}
 	free(intfile);
 	free(image.det);
-	free(powder->data);
-	free(powder);
 	cell_free(input_cell);
 	free(intensities);
 	free(outfile);
diff --git a/src/process_image.c b/src/process_image.c
index e7e3aa78..4fd19176 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -135,8 +135,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	int any_crystals;
 
 	image.features = NULL;
-	image.data = NULL;
-	image.flags = NULL;
 	image.copyme = iargs->copyme;
 	image.id = cookie;
 	image.filename = pargs->filename_p_e->filename;
@@ -343,8 +341,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	free(image.dp);
 	free(image.bad);
 
-	free(image.data);
-	if ( image.flags != NULL ) free(image.flags);
 	image_feature_list_free(image.features);
 	free_detector_geometry(image.det);
 	hdfile_close(hdfile);
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
index dbe59ce2..7c7c8e0c 100644
--- a/tests/gpu_sim_check.c
+++ b/tests/gpu_sim_check.c
@@ -3,11 +3,11 @@
  *
  * Check that GPU simulation agrees with CPU version
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2012-2014 Thomas White <taw@physics.org>
+ *   2012-2015 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -110,11 +110,19 @@ int main(int argc, char *argv[])
 	det = calloc(1, sizeof(struct detector));
 	det->n_panels = 2;
 	det->panels = calloc(2, sizeof(struct panel));
+	det->max_fs = 1024;
+	det->max_ss = 1024;
 
 	det->panels[0].min_fs = 0;
 	det->panels[0].max_fs = 1023;
 	det->panels[0].min_ss = 0;
 	det->panels[0].max_ss = 511;
+	det->panels[0].orig_min_fs = 0;
+	det->panels[0].orig_max_fs = 1023;
+	det->panels[0].orig_min_ss = 0;
+	det->panels[0].orig_max_ss = 511;
+	det->panels[0].w = 1024;
+	det->panels[0].h = 512;
 	det->panels[0].fsx = 1;
 	det->panels[0].fsy = 0;
 	det->panels[0].ssx = 0;
@@ -128,11 +136,18 @@ int main(int argc, char *argv[])
 	det->panels[0].clen = 100.0e-3;
 	det->panels[0].res = 9090.91;
 	det->panels[0].adu_per_eV = 1.0;
+	det->panels[0].data = NULL;
 
 	det->panels[1].min_fs = 0;
 	det->panels[1].max_fs = 1023;
 	det->panels[1].min_ss = 512;
 	det->panels[1].max_ss = 1023;
+	det->panels[1].orig_min_fs = 0;
+	det->panels[1].orig_max_fs = 1023;
+	det->panels[1].orig_min_ss = 512;
+	det->panels[1].orig_max_ss = 1023;
+	det->panels[1].w = 1024;
+	det->panels[1].h = 512;
 	det->panels[1].fsx = 1;
 	det->panels[1].fsy = 0;
 	det->panels[1].ssx = 0;
@@ -145,12 +160,15 @@ int main(int argc, char *argv[])
 	det->panels[1].cny = sep;
 	det->panels[1].clen = 100.0e-3;
 	det->panels[1].res = 9090.91;
-	det->panels[0].adu_per_eV = 1.0;
+	det->panels[1].adu_per_eV = 1.0;
+	det->panels[1].data = NULL;
 
 	cpu_image.det = det;
 	gpu_image.det = det;
 	cpu_image.beam = NULL;
 	gpu_image.beam = NULL;
+	cpu_image.spectrum = NULL;
+	gpu_image.spectrum = NULL;
 
 	cpu_image.lambda = ph_en_to_lambda(eV_to_J(6000));
 	gpu_image.lambda = ph_en_to_lambda(eV_to_J(6000));
@@ -171,6 +189,20 @@ int main(int argc, char *argv[])
 
 	sym = get_pointgroup("1");
 
+	cpu_image.dp = malloc(det->n_panels * sizeof(float *));
+	if ( cpu_image.dp == NULL ) {
+		ERROR("Couldn't allocate memory for result.\n");
+		return 1;
+	}
+	for ( i=0; i<det->n_panels; i++ ) {
+		struct panel *p = &det->panels[i];
+		cpu_image.dp[i] = calloc(p->w * p->h, sizeof(float));
+		if ( cpu_image.dp[i] == NULL ) {
+			ERROR("Couldn't allocate memory for panel %i\n", i);
+			return 1;
+		}
+	}
+
 	start = get_hires_seconds();
 	get_diffraction(&cpu_image, 8, 8, 8, NULL, NULL, NULL, cell,
 	                GRADIENT_MOSAIC, sym, 1);
@@ -184,18 +216,25 @@ int main(int argc, char *argv[])
 	gpu_min = +INFINITY;  gpu_max = -INFINITY;  gpu_tot = 0.0;
 	cpu_min = +INFINITY;  cpu_max = -INFINITY;  cpu_tot = 0.0;
 	dev = 0.0;
-	for ( i=0; i<1024*1024; i++ ) {
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		int j;
+		struct panel *p = &det->panels[i];
+
+		for ( j=0; j<p->w*p->h; j++ ) {
+
+			const double cpu = cpu_image.dp[i][j];
+			const double gpu = gpu_image.dp[i][j];
 
-		const double cpu = cpu_image.data[i];
-		const double gpu = gpu_image.data[i];
+			if ( cpu > cpu_max ) cpu_max = cpu;
+			if ( cpu < cpu_min ) cpu_min = cpu;
+			if ( gpu > gpu_max ) gpu_max = gpu;
+			if ( gpu < gpu_min ) gpu_min = gpu;
+			gpu_tot += gpu;
+			cpu_tot += cpu;
+			dev += fabs(gpu - cpu);
 
-		if ( cpu > cpu_max ) cpu_max = cpu;
-		if ( cpu < cpu_min ) cpu_min = cpu;
-		if ( gpu > gpu_max ) gpu_max = gpu;
-		if ( gpu < gpu_min ) gpu_min = gpu;
-		gpu_tot += gpu;
-		cpu_tot += cpu;
-		dev += fabs(gpu - cpu);
+		}
 
 	}
 	perc = 100.0*dev/cpu_tot;
@@ -204,25 +243,22 @@ int main(int argc, char *argv[])
 	STATUS("CPU: min=%8e, max=%8e, total=%8e\n", cpu_min, cpu_max, cpu_tot);
 	STATUS("dev = %8e (%5.2f%% of CPU total)\n", dev, perc);
 
-	cell_free(cell);
-	free_detector_geometry(det);
-
-	if ( perc > 1.0 ) {
+	if ( 1 ) {
 
 		STATUS("Test failed!  I'm writing cpu-sim.h5 and gpu-sim.h5"
 		       " for you to inspect.\n");
 
-		hdf5_write("cpu-sim.h5", cpu_image.data, cpu_image.width,
-		            cpu_image.height, H5T_NATIVE_FLOAT);
-
-		hdf5_write("gpu-sim.h5", gpu_image.data, gpu_image.width,
-		            gpu_image.height, H5T_NATIVE_FLOAT);
+		hdf5_write_image("cpu-sim.h5", &cpu_image, NULL);
+		hdf5_write_image("gpu-sim.h5", &gpu_image, NULL);
 
 		return 1;
 
 	}
 
 	gsl_rng_free(rng);
+	cell_free(cell);
+	free_detector_geometry(det);
+
 
 	return 0;
 }
diff --git a/tests/integration_check.c b/tests/integration_check.c
index ccb613c3..2acd0f34 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -65,7 +65,6 @@ int main(int argc, char *argv[])
 	fclose(fh);
 	gsl_rng_set(rng, seed);
 
-	image.flags = NULL;
 	image.beam = NULL;
 	image.lambda = ph_eV_to_lambda(9000.0);
 
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
index c7ebd8b2..f66b6513 100644
--- a/tests/prof2d_check.c
+++ b/tests/prof2d_check.c
@@ -77,7 +77,6 @@ int main(int argc, char *argv[])
 	fclose(fh);
 	gsl_rng_set(rng, seed);
 
-	image.flags = NULL;
 	image.beam = NULL;
 	image.lambda = ph_eV_to_lambda(9000.0);
 	image.bw = 0.000001;
diff --git a/tests/ring_check.c b/tests/ring_check.c
index 8ced12d5..99c433c2 100644
--- a/tests/ring_check.c
+++ b/tests/ring_check.c
@@ -62,7 +62,7 @@ static void third_integration_check(struct image *image, int n_trials,
 
 		for ( fs=0; fs<image->width; fs++ ) {
 		for ( ss=0; ss<image->height; ss++ ) {
-			image->data[fs+image->width*ss]
+			image->dp[0][fs+image->width*ss]
 			                           = poisson_noise(rng, 1000.0);
 		}
 		}
@@ -121,9 +121,9 @@ static void fourth_integration_check(struct image *image, int n_trials,
 		for ( fs=0; fs<image->width; fs++ ) {
 		for ( ss=0; ss<image->height; ss++ ) {
 			int idx = fs+image->width*ss;
-			image->data[idx] = poisson_noise(rng, 1000.0);
+			image->dp[0][idx] = poisson_noise(rng, 1000.0);
 			if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
-			image->data[idx] += 1000.0;
+			image->dp[0][idx] += 1000.0;
 			pcount++;
 		}
 		}
@@ -170,7 +170,6 @@ int main(int argc, char *argv[])
 	int r, npx;
 	double ex;
 	gsl_rng *rng;
-	int *bad;
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
@@ -179,10 +178,10 @@ int main(int argc, char *argv[])
 	fclose(fh);
 	gsl_rng_set(rng, seed);
 
-	image.data = malloc(128*128*sizeof(float));
-	bad = calloc(128*128, sizeof(int));
-	image.bad = &bad;
-	image.flags = NULL;
+	image.dp = malloc(sizeof(float *));
+	image.dp[0] = malloc(128*128*sizeof(float));
+	image.bad = malloc(sizeof(uint16_t *));
+	image.bad[0] = calloc(128*128, sizeof(uint16_t));
 	image.beam = NULL;
 	image.lambda = ph_eV_to_lambda(1000.0);
 
@@ -191,9 +190,9 @@ int main(int argc, char *argv[])
 	image.det->panels = calloc(1, sizeof(struct panel));
 
 	image.det->panels[0].min_fs = 0;
-	image.det->panels[0].max_fs = 128;
+	image.det->panels[0].max_fs = 127;
 	image.det->panels[0].min_ss = 0;
-	image.det->panels[0].max_ss = 128;
+	image.det->panels[0].max_ss = 127;
 	image.det->panels[0].fsx = 1.0;
 	image.det->panels[0].fsy = 0.0;
 	image.det->panels[0].ssx = 0.0;
@@ -206,12 +205,14 @@ int main(int argc, char *argv[])
 	image.det->panels[0].cny = -64.0;
 	image.det->panels[0].clen = 1.0;
 	image.det->panels[0].res = 1.0;
+	image.det->panels[0].w = 128;
+	image.det->panels[0].h = 128;
 	image.det->panels[0].adu_per_eV = 1.0/1000.0;  /* -> 1 adu per photon */
 	image.det->panels[0].max_adu = +INFINITY;  /* No cutoff */
 
 	image.width = 128;
 	image.height = 128;
-	memset(image.data, 0, 128*128*sizeof(float));
+	memset(image.dp[0], 0, 128*128*sizeof(float));
 
 	image.n_crystals = 0;
 	image.crystals = NULL;
@@ -238,7 +239,7 @@ int main(int argc, char *argv[])
 	for ( fs=0; fs<image.width; fs++ ) {
 	for ( ss=0; ss<image.height; ss++ ) {
 		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
-		image.data[fs+image.width*ss] = 1000.0;
+		image.dp[0][fs+image.width*ss] = 1000.0;
 		npx++;
 	}
 	}
@@ -278,9 +279,9 @@ int main(int argc, char *argv[])
 	npx = 0;
 	for ( fs=0; fs<image.width; fs++ ) {
 	for ( ss=0; ss<image.height; ss++ ) {
-		image.data[fs+image.width*ss] = 1000.0;
+		image.dp[0][fs+image.width*ss] = 1000.0;
 		if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
-		image.data[fs+image.width*ss] += 1000.0;
+		image.dp[0][fs+image.width*ss] += 1000.0;
 		npx++;
 	}
 	}
@@ -310,11 +311,11 @@ int main(int argc, char *argv[])
 
 	}
 
-
 	free(image.beam);
 	free(image.det->panels);
 	free(image.det);
-	free(image.data);
+	free(image.dp[0]);
+	free(image.dp);
 	gsl_rng_free(rng);
 
 	if ( fail ) return 1;