Merge branch 'tom/imagedata'

5 files changed, 194 insertions, 149 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);
 
 	}