Merge branch 'master' into tom/speed

Conflicts:
	libcrystfel/src/peaks.c
	src/indexamajig.c

14 files changed, 290 insertions, 152 deletions

diff --git a/AUTHORS b/AUTHORS
index 54a29e7b..5ce6a9c1 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,6 +1,9 @@
 * Thomas White <taw@physics.org>
   Lead author, architecture, "vision" and so on.
 
+* Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+  Peak integration improvements
+
 * Richard Kirian <rkirian@asu.edu>
   MOSFLM auto-indexing, bug fixing.
 
diff --git a/README b/README
index 700e646c..23237849 100644
--- a/README
+++ b/README
@@ -7,6 +7,7 @@ Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 Authors:
   Thomas White <taw@physics.org>
   Richard Kirian <rkirian@asu.edu>
+  Kenneth Beyerlein <kenneth.beyerlein@desy.de>
   Andrew Aquila <andrew.aquila@cfel.de>
   Andrew Martin <andrew.martin@desy.de>
   Lorenzo Galli <lorenzo.galli@desy.de>
diff --git a/data/diffraction.cl b/data/diffraction.cl
index f6318378..d128446e 100644
--- a/data/diffraction.cl
+++ b/data/diffraction.cl
@@ -33,7 +33,7 @@ const sampler_t sampler_c = CLK_NORMALIZED_COORDS_TRUE
 
 
 float4 get_q(float fs, float ss, float res, float clen, float k,
-             float *ttp, float corner_x, float corner_y,
+             float corner_x, float corner_y,
              float fsx, float fsy, float ssx, float ssy)
 {
 	float rx, ry, r;
@@ -51,7 +51,6 @@ float4 get_q(float fs, float ss, float res, float clen, float k,
 	r = sqrt(pow(rx, 2.0f) + pow(ry, 2.0f));
 
 	tt = atan2(r, clen);
-	*ttp = tt;
 
 	az = atan2(ry, rx);
 
@@ -211,7 +210,7 @@ float molecule_factor(global float *intensities, global float *flags,
 }
 
 
-kernel void diffraction(global float *diff, global float *ttp, float k,
+kernel void diffraction(global float *diff, float k,
                         int w, float corner_x, float corner_y,
                         float res, float clen, float16 cell,
                         global float *intensities,
@@ -219,7 +218,8 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
                         read_only image2d_t func_b,
                         read_only image2d_t func_c,
                         global float *flags,
-                        float fsx, float fsy, float ssx, float ssy)
+                        float fsx, float fsy, float ssx, float ssy,
+                        float xo, float yo)
 {
 	float tt;
 	float fs, ss;
@@ -230,11 +230,11 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
 	int idx;
 
 	/* Calculate fractional coordinates in fs/ss */
-	fs = convert_float(get_global_id(0));
-	ss = convert_float(get_global_id(1));
+	fs = convert_float(get_global_id(0)) + xo;
+	ss = convert_float(get_global_id(1)) + yo;
 
 	/* Get the scattering vector */
-	q = get_q(fs, ss, res, clen, k, &tt,
+	q = get_q(fs, ss, res, clen, k,
 	          corner_x, corner_y, fsx, fsy, ssx, ssy);
 
 	/* Calculate the diffraction */
@@ -245,5 +245,4 @@ kernel void diffraction(global float *diff, global float *ttp, float k,
 	/* Write the value to memory */
 	idx = convert_int_rtz(fs) + w*convert_int_rtz(ss);
 	diff[idx] = I_molecule * I_lattice;
-	ttp[idx] = tt;
 }
diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5
index 006af722..68abb341 100644
--- a/doc/man/crystfel_geometry.5
+++ b/doc/man/crystfel_geometry.5
@@ -47,7 +47,7 @@ of pixel values in the HDF5 file, defined in terms only of the "fast scan" and
 +x completes the right-handed coordinate system.
 
 .PP
-Naively speaking, this means that CrystFEL at the images from the "into the
+Naively speaking, this means that CrystFEL looks at the images from the "into the
 beam" perspective, but please avoid thinking of things in this way.  It's much
 better to consider the precise way in which the coordinates are mapped.
 
@@ -116,6 +116,8 @@ panel0/ss = -x
 .br
 ; the HDF5 file, is now given a position in the lab coordinate system.
 .br
+; Units are pixels.
+.br
 ; Note that "first point in the panel" is a conceptual simplification.  We refer
 .br
 ; to that corner, and to the very corner of the pixel - NOT, for example, to the
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index b7d809df..f212821d 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -608,8 +608,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
 		panel->res = atof(val);
 	} else if ( strcmp(key, "max_adu") == 0 ) {
 		panel->max_adu = atof(val);
-	} else if ( strcmp(key, "peak_sep") == 0 ) {
-		panel->peak_sep = atof(val);
 	} else if ( strcmp(key, "badrow_direction") == 0 ) {
 		panel->badrow = val[0]; /* First character only */
 		if ( (panel->badrow != 'x') && (panel->badrow != 'y')
@@ -687,8 +685,6 @@ static void parse_toplevel(struct detector *det, const char *key,
 			det->mask_good = v;
 		}
 
-	} else if ( strcmp(key, "peak_sep") == 0 ) {
-		det->defaults.peak_sep = atof(val);
 	} else if ( strcmp(key, "coffset") == 0 ) {
 		det->defaults.coffset = atof(val);
 	} else if ( parse_field_for_panel(&det->defaults, key, val, det) ) {
@@ -738,7 +734,6 @@ struct detector *get_detector_geometry(const char *filename)
 	det->defaults.res = -1.0;
 	det->defaults.badrow = '-';
 	det->defaults.no_index = 0;
-	det->defaults.peak_sep = 50.0;
 	det->defaults.fsx = 1.0;
 	det->defaults.fsy = 0.0;
 	det->defaults.ssx = 0.0;
@@ -887,7 +882,6 @@ struct detector *get_detector_geometry(const char *filename)
 		}
 		/* It's OK if the badrow direction is '0' */
 		/* It's not a problem if "no_index" is still zero */
-		/* The default peak_sep is OK (maybe) */
 		/* The default transformation matrix is at least valid */
 
 		if ( det->panels[i].max_fs > max_fs ) {
@@ -1250,7 +1244,6 @@ int write_detector_geometry(const char *filename, struct detector *det)
 		fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
 		fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
 		fprintf(fh, "%s/res = %g\n", p->name, p->res);
-		fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
 		fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
 		fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy);
 		fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index 2b4ac5c6..b740965c 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -62,7 +62,6 @@ struct panel
 	double   res;      /* Resolution in pixels per metre */
 	char     badrow;   /* 'x' or 'y' */
 	int      no_index; /* Don't index peaks in this panel if non-zero */
-	double   peak_sep; /* Characteristic peak separation */
 	char    *rigid_group;  /* Rigid group, or -1 for none */
 	double   adu_per_eV;   /* Number of ADU per eV */
 	double   max_adu;  /* Treat pixel as unreliable if higher than this */
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 583dff4c..ffe77290 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -9,6 +9,7 @@
  *
  * Authors:
  *   2010-2012 Thomas White <taw@physics.org>
+ *   2012      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *   2011      Andrew Martin <andrew.martin@desy.de>
  *   2011      Richard Kirian
  *
@@ -98,7 +99,6 @@ static int cull_peaks_in_panel(struct image *image, struct panel *p)
 		if ( ncol <= 3 ) continue;
 
 		/* Yes?  Delete them all... */
-		nelim = 0;
 		for ( j=0; j<n; j++ ) {
 			struct imagefeature *g;
 			g = image_get_feature(image->features, j);
@@ -146,13 +146,76 @@ static int cull_peaks(struct image *image)
 }
 
 
+/* cfs, css relative to panel origin */
+static int *make_BgMask(struct image *image, struct panel *p,
+                        double ir_out, int cfs, int css, double ir_in)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	int *mask;
+	int w, h;
+
+	w = p->max_fs - p->min_fs + 1;
+	h = p->max_ss - p->min_ss + 1;
+	mask = calloc(w*h, sizeof(int));
+	if ( mask == NULL ) return NULL;
+
+	/* Loop over all reflections */
+	for ( refl = first_refl(image->reflections, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		struct panel *p2;
+		double pk2_fs, pk2_ss;
+		signed int dfs, dss;
+		double pk2_cfs, pk2_css;
+
+		get_detector_pos(refl, &pk2_fs, &pk2_ss);
+
+		/* Determine if reflection is in the same panel */
+		p2 = find_panel(image->det, pk2_fs, pk2_ss);
+		if ( p2 != p ) continue;
+
+		pk2_cfs = pk2_fs - p->min_fs;
+		pk2_css = pk2_ss - p->min_ss;
+
+		for ( dfs=-ir_in; dfs<=ir_in; dfs++ ) {
+		for ( dss=-ir_in; dss<=ir_in; dss++ ) {
+
+			signed int fs, ss;
+
+			/* In peak region for this peak? */
+			if ( dfs*dfs + dss*dss > ir_in*ir_in ) continue;
+
+			fs = pk2_cfs + dfs;
+			ss = pk2_css + dss;
+
+			/* On panel? */
+			if ( fs >= w ) continue;
+			if ( ss >= h ) continue;
+			if ( fs < 0 ) continue;
+			if ( ss < 0 ) continue;
+
+			mask[fs + ss*w] = 1;
+
+		}
+		}
+
+	}
+
+	return mask;
+}
+
+
 /* Returns non-zero if peak has been vetoed.
  * i.e. don't use result if return value is not zero. */
-int integrate_peak(struct image *image, int cfs, int css,
-                   double *pfs, double *pss, double *intensity, double *sigma,
-                   double ir_inn, double ir_mid, double ir_out)
+static int integrate_peak(struct image *image, int cfs, int css,
+                          double *pfs, double *pss,
+                          double *intensity, double *sigma,
+                          double ir_inn, double ir_mid, double ir_out,
+                          int use_max_adu)
 {
-	signed int fs, ss;
+	signed int dfs, dss;
 	double lim_sq, out_lim_sq, mid_lim_sq;
 	double pk_total;
 	int pk_counts;
@@ -164,11 +227,21 @@ int integrate_peak(struct image *image, int cfs, int css,
 	double bg_tot_sq = 0.0;
 	double var;
 	double aduph;
+	int *bgPkMask;
+	int p_cfs, p_css, p_w, p_h;
 
 	p = find_panel(image->det, cfs, css);
 	if ( p == NULL ) return 1;
 	if ( p->no_index ) return 1;
 
+	/* 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;
+	bgPkMask = make_BgMask(image, p, ir_out, p_cfs, p_css, ir_inn);
+	if ( bgPkMask == NULL ) return 1;
+
 	aduph = p->adu_per_eV * ph_lambda_to_eV(image->lambda);
 
 	lim_sq = pow(ir_inn, 2.0);
@@ -176,23 +249,27 @@ int integrate_peak(struct image *image, int cfs, int css,
 	out_lim_sq = pow(ir_out, 2.0);
 
 	/* Estimate the background */
-	for ( fs=-ir_out; fs<=+ir_out; fs++ ) {
-	for ( ss=-ir_out; ss<=+ir_out; ss++ ) {
+	for ( dfs=-ir_out; dfs<=+ir_out; dfs++ ) {
+	for ( dss=-ir_out; dss<=+ir_out; dss++ ) {
 
 		double val;
 		uint16_t flags;
-		struct panel *p2;
 		int idx;
 
 		/* Restrict to annulus */
-		if ( fs*fs + ss*ss > out_lim_sq ) continue;
-		if ( fs*fs + ss*ss < mid_lim_sq ) continue;
+		if ( dfs*dfs + dss*dss > out_lim_sq ) continue;
+		if ( dfs*dfs + dss*dss < mid_lim_sq ) continue;
 
 		/* Strayed off one panel? */
-		p2 = find_panel(image->det, fs+cfs, ss+css);
-		if ( p2 != p ) return 1;
+		if ( p_cfs+dfs >= p_w ) continue;
+		if ( p_css+dss >= p_h ) continue;
+		if ( p_cfs+dfs < 0 ) continue;
+		if ( p_css+dss < 0 ) continue;
+
+		/* Check if there is a peak in the background region */
+		if ( bgPkMask[(p_cfs+dfs) + p_w*(p_css+dss)] ) continue;
 
-		idx = fs+cfs+image->width*(ss+css);
+		idx = dfs+cfs+image->width*(dss+css);
 
 		/* Veto this peak if we tried to integrate in a bad region */
 		if ( image->flags != NULL ) {
@@ -211,7 +288,7 @@ int integrate_peak(struct image *image, int cfs, int css,
 		val = image->data[idx];
 
 		/* Veto peak if it contains saturation in bg region */
-		if ( val > p->max_adu ) return 1;
+		if ( use_max_adu && (val > p->max_adu) ) return 1;
 
 		bg_tot += val;
 		bg_tot_sq += pow(val, 2.0);
@@ -228,22 +305,23 @@ int integrate_peak(struct image *image, int cfs, int css,
 	pk_total = 0.0;
 	pk_counts = 0;
 	fsct = 0.0;  ssct = 0.0;
-	for ( fs=-ir_inn; fs<=+ir_inn; fs++ ) {
-	for ( ss=-ir_inn; ss<=+ir_inn; ss++ ) {
+	for ( dfs=-ir_inn; dfs<=+ir_inn; dfs++ ) {
+	for ( dss=-ir_inn; dss<=+ir_inn; dss++ ) {
 
 		double val;
 		uint16_t flags;
-		struct panel *p2;
 		int idx;
 
 		/* Inner mask radius */
-		if ( fs*fs + ss*ss > lim_sq ) continue;
+		if ( dfs*dfs + dss*dss > lim_sq ) continue;
 
 		/* Strayed off one panel? */
-		p2 = find_panel(image->det, fs+cfs, ss+css);
-		if ( p2 != p ) return 1;
+		if ( p_cfs+dfs >= p_w ) continue;
+		if ( p_css+dss >= p_h ) continue;
+		if ( p_cfs+dfs < 0 ) continue;
+		if ( p_css+dss < 0 ) continue;
 
-		idx = fs+cfs+image->width*(ss+css);
+		idx = dfs+cfs+image->width*(dss+css);
 
 		/* Veto this peak if we tried to integrate in a bad region */
 		if ( image->flags != NULL ) {
@@ -262,13 +340,13 @@ int integrate_peak(struct image *image, int cfs, int css,
 		val = image->data[idx] - bg_mean;
 
 		/* Veto peak if it contains saturation */
-		if ( image->data[idx] > p->max_adu ) return 1;
+		if ( use_max_adu && (image->data[idx] > p->max_adu) ) return 1;
 
 		pk_counts++;
 		pk_total += val;
 
-		fsct += val*(cfs+fs);
-		ssct += val*(css+ss);
+		fsct += val*(cfs+dfs);
+		ssct += val*(css+dss);
 
 	}
 	}
@@ -285,6 +363,8 @@ int integrate_peak(struct image *image, int cfs, int css,
 	if ( intensity != NULL ) *intensity = pk_total;
 	if ( sigma != NULL ) *sigma = sqrt(var);
 
+	free(bgPkMask);
+
 	return 0;
 }
 
@@ -309,7 +389,6 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 	int nrej_snr = 0;
 	int nacc = 0;
 	int ncull;
-	const int pws = p->peak_sep/2;
 
 	data = image->data;
 	stride = image->width;
@@ -352,16 +431,14 @@ 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-pws/2,
-			                   p->min_ss);
-			      s_ss<=smallest(mask_ss+pws/2,
-			                     p->max_ss);
-			      s_ss++ ) {
-			for ( s_fs=biggest(mask_fs-pws/2,
-			                   p->min_fs);
-			      s_fs<=smallest(mask_fs+pws/2,
-			                     p->max_fs);
-			      s_fs++ ) {
+			for ( s_ss=biggest(mask_ss-ir_inn, p->min_ss);
+			      s_ss<=smallest(mask_ss+ir_inn, p->max_ss);
+			      s_ss++ )
+			{
+			for ( s_fs=biggest(mask_fs-ir_inn, p->min_fs);
+			      s_fs<=smallest(mask_fs+ir_inn, p->max_fs);
+			      s_fs++ )
+			{
 
 				if ( data[s_fs+stride*s_ss] > max ) {
 					max = data[s_fs+stride*s_ss];
@@ -374,8 +451,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 			}
 
 			/* Abort if drifted too far from the foot point */
-			if ( distance(mask_fs, mask_ss, fs, ss) >
-			     p->peak_sep/2.0 )
+			if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn )
 			{
 				break;
 			}
@@ -383,7 +459,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		} while ( did_something );
 
 		/* Too far from foot point? */
-		if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) {
+		if ( distance(mask_fs, mask_ss, fs, ss) > ir_inn ) {
 			nrej_dis++;
 			continue;
 		}
@@ -397,7 +473,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		/* Centroid peak and get better coordinates. */
 		r = integrate_peak(image, mask_fs, mask_ss,
 		                   &f_fs, &f_ss, &intensity, &sigma,
-		                   ir_inn, ir_mid, ir_out);
+		                   ir_inn, ir_mid, ir_out, 0);
 
 		if ( r ) {
 			/* Bad region - don't detect peak */
@@ -419,7 +495,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 
 		/* Check for a nearby feature */
 		image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
-		if ( d < p->peak_sep/2.0 ) {
+		if ( d < 2.0*ir_inn ) {
 			nrej_pro++;
 			continue;
 		}
@@ -445,6 +521,12 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 //	       "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
 //	       nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
 //	       nrej_snr, min_snr, ncull);
+
+	if ( ncull != 0 ) {
+		STATUS("WARNING: %i peaks were badrow culled.  This feature"
+		       " should not usually be used.\nConsider setting"
+		       " badrow=- in the geometry file.\n", ncull);
+	}
 }
 
 
@@ -655,7 +737,11 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 		}
 
 		r = integrate_peak(image, pfs, pss, &fs, &ss,
-		                   &intensity, &sigma, ir_inn, ir_mid, ir_out);
+		                   &intensity, &sigma, ir_inn, ir_mid, ir_out,
+		                   1);
+
+		/* I/sigma(I) cutoff */
+		if ( intensity/sigma < min_snr ) r = 1;
 
 		/* Record intensity and set redundancy to 1 on success */
 		if ( r == 0 ) {
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index cea25bf3..64154774 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -900,6 +900,7 @@ struct _reflistiterator {
  **/
 Reflection *first_refl(RefList *list, RefListIterator **piter)
 {
+	Reflection *refl;
 	RefListIterator *iter;
 
 	iter = malloc(sizeof(struct _reflistiterator));
@@ -908,7 +909,9 @@ Reflection *first_refl(RefList *list, RefListIterator **piter)
 	iter->stack_ptr = 0;
 	*piter = iter;
 
-	Reflection *refl = list->head;
+	if ( list == NULL ) return NULL;
+
+	refl = list->head;
 
 	do {
 
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 81f4d722..ce0ee33e 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -474,6 +474,7 @@ int skip_some_files(FILE *fh, int n)
 
 		rval = fgets(line, 1023, fh);
 		if ( rval == NULL ) continue;
+		chomp(line);
 		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
 
 	} while ( rval != NULL );
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 2b8e97e5..15ff39cc 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -167,57 +167,13 @@ static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
 }
 
 
-void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                         int na, int nb, int nc, UnitCell *ucell)
+static void do_panels(struct gpu_context *gctx, struct image *image,
+                      int *n_inf, int *n_neg, int *n_nan, double k, double frac,
+                      double xo, double yo)
 {
-	cl_int err;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
 	int i;
-	cl_float16 cell;
-	cl_int4 ncells;
-	int n_inf = 0;
-	int n_neg = 0;
-	int n_nan = 0;
 
-	if ( gctx == NULL ) {
-		ERROR("GPU setup failed.\n");
-		return;
-	}
-
-	cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	cell.s[0] = ax;  cell.s[1] = ay;  cell.s[2] = az;
-	cell.s[3] = bx;  cell.s[4] = by;  cell.s[5] = bz;
-	cell.s[6] = cx;  cell.s[7] = cy;  cell.s[8] = cz;
-
-	ncells.s[0] = na;
-	ncells.s[1] = nb;
-	ncells.s[2] = nc;
-	ncells.s[3] = 0;  /* unused */
-
-	/* Ensure all required LUTs are available */
-	check_sinc_lut(gctx, na);
-	check_sinc_lut(gctx, nb);
-	check_sinc_lut(gctx, nc);
-
-	if ( set_arg_float(gctx, 2, 1.0/image->lambda) ) return;
-	if ( set_arg_mem(gctx, 9, gctx->intensities) ) return;
-	if ( set_arg_mem(gctx, 10, gctx->sinc_luts[na-1]) ) return;
-	if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nb-1]) ) return;
-	if ( set_arg_mem(gctx, 12, gctx->sinc_luts[nc-1]) ) return;
-	if ( set_arg_mem(gctx, 13, gctx->flags) ) return;
-
-	/* Unit cell */
-	err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set unit cell: %s\n", clError(err));
-		return;
-	}
-
-	/* Allocate memory for the result */
-	image->data = calloc(image->width * image->height, sizeof(float));
-	image->twotheta = calloc(image->width * image->height, sizeof(double));
+	if ( set_arg_float(gctx, 1, k) ) return;
 
 	/* Iterate over panels */
 	for ( i=0; i<image->det->n_panels; i++ ) {
@@ -225,14 +181,12 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 		size_t dims[2];
 		size_t ldims[2] = {1, 1};
 		struct panel *p;
-		cl_mem tt;
-		size_t tt_size;
 		cl_mem diff;
 		size_t diff_size;
 		float *diff_ptr;
-		float *tt_ptr;
 		int pan_width, pan_height;
 		int fs, ss;
+		cl_int err;
 
 		p = &image->det->panels[i];
 
@@ -247,25 +201,19 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 			ERROR("Couldn't allocate diffraction memory\n");
 			return;
 		}
-		tt_size = pan_width * pan_height * sizeof(cl_float);
-		tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size,
-			            NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't allocate twotheta memory\n");
-			return;
-		}
 
 		if ( set_arg_mem(gctx, 0, diff) ) return;
-		if ( set_arg_mem(gctx, 1, tt) ) return;
-		if ( set_arg_int(gctx, 3, pan_width) ) return;
-		if ( set_arg_float(gctx, 4, p->cnx) ) return;
-		if ( set_arg_float(gctx, 5, p->cny) ) return;
-		if ( set_arg_float(gctx, 6, p->res) ) return;
-		if ( set_arg_float(gctx, 7, p->clen) ) return;
-		if ( set_arg_float(gctx, 14, p->fsx) ) return;
-		if ( set_arg_float(gctx, 15, p->fsy) ) return;
-		if ( set_arg_float(gctx, 16, p->ssx) ) return;
-		if ( set_arg_float(gctx, 17, p->ssy) ) return;
+		if ( set_arg_int(gctx, 2, pan_width) ) return;
+		if ( set_arg_float(gctx, 3, p->cnx) ) return;
+		if ( set_arg_float(gctx, 4, p->cny) ) return;
+		if ( set_arg_float(gctx, 5, p->res) ) return;
+		if ( set_arg_float(gctx, 6, p->clen) ) return;
+		if ( set_arg_float(gctx, 13, p->fsx) ) return;
+		if ( set_arg_float(gctx, 14, p->fsy) ) return;
+		if ( set_arg_float(gctx, 15, p->ssx) ) return;
+		if ( set_arg_float(gctx, 16, p->ssy) ) return;
+		if ( set_arg_float(gctx, 17, xo) ) return;
+		if ( set_arg_float(gctx, 18, yo) ) return;
 
 		dims[0] = pan_width;
 		dims[1] = pan_height;
@@ -288,43 +236,124 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 			      clError(err));
 			return;
 		}
-		tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ,
-		                            0, tt_size, 0, NULL, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't map tt buffer\n");
-			return;
-		}
 
 		for ( fs=0; fs<pan_width; fs++ ) {
 		for ( ss=0; ss<pan_height; ss++ ) {
 
-			float val, tt;
+			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++;
-			tt = tt_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->twotheta[tfs + image->width*tss] = tt;
+			image->data[tfs + image->width*tss] += frac*val;
 
 		}
 		}
 
 		clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr,
 		                        0, NULL, NULL);
-		clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr,
-		                        0, NULL, NULL);
 
 		clReleaseMemObject(diff);
-		clReleaseMemObject(tt);
 
 	}
+}
+
+
+void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
+                         int na, int nb, int nc, UnitCell *ucell)
+{
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	cl_float16 cell;
+	cl_int4 ncells;
+	cl_int err;
+	int n_inf = 0;
+	int n_neg = 0;
+	int n_nan = 0;
+	int fs, ss;
+	const int nkstep = 10;
+	const int nxs = 4;
+	const int nys = 4;
+	int kstep;
+	double klow, khigh, kinc, w;
+	double frac;
+	int xs, ys;
+	int n, ntot;
+
+	if ( gctx == NULL ) {
+		ERROR("GPU setup failed.\n");
+		return;
+	}
 
+	cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	cell.s[0] = ax;  cell.s[1] = ay;  cell.s[2] = az;
+	cell.s[3] = bx;  cell.s[4] = by;  cell.s[5] = bz;
+	cell.s[6] = cx;  cell.s[7] = cy;  cell.s[8] = cz;
+
+	ncells.s[0] = na;
+	ncells.s[1] = nb;
+	ncells.s[2] = nc;
+	ncells.s[3] = 0;  /* unused */
+
+	/* Ensure all required LUTs are available */
+	check_sinc_lut(gctx, na);
+	check_sinc_lut(gctx, nb);
+	check_sinc_lut(gctx, nc);
+
+	if ( set_arg_mem(gctx, 8, gctx->intensities) ) return;
+	if ( set_arg_mem(gctx, 9, gctx->sinc_luts[na-1]) ) return;
+	if ( set_arg_mem(gctx, 10, gctx->sinc_luts[nb-1]) ) return;
+	if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nc-1]) ) return;
+	if ( set_arg_mem(gctx, 12, gctx->flags) ) return;
+
+	/* Unit cell */
+	err = clSetKernelArg(gctx->kern, 7, sizeof(cl_float16), &cell);
+	if ( err != CL_SUCCESS ) {
+		ERROR("Couldn't set unit cell: %s\n", clError(err));
+		return;
+	}
+
+	/* Allocate memory for the result */
+	image->data = calloc(image->width * image->height, sizeof(float));
+
+	w = image->lambda * image->bw;
+	klow  = 1.0 / (image->lambda + (w/2.0));  /* Smallest k */
+	khigh = 1.0 / (image->lambda - (w/2.0));  /* Largest k */
+	kinc = (khigh-klow) / (nkstep+1);
+
+	ntot = nkstep * nxs * nys;
+	frac = 1.0 / ntot;
+
+	n = 0;
+	for ( xs=0; xs<nxs; xs++ ) {
+	for ( ys=0; ys<nys; ys++ ) {
+
+		double xo, yo;
+
+		xo = (1.0/nxs) * xs;
+		yo = (1.0/nys) * ys;
+
+		for ( kstep=0; kstep<nkstep; kstep++ ) {
+
+			double k;
+
+			k = klow + kstep*kinc;
+
+			do_panels(gctx, image, &n_inf, &n_neg, &n_nan, k, frac,
+			          xo, yo);
+
+			n++;
+			progress_bar(n, ntot, "Simulating");
+
+		}
+	}
+	}
 
 	if ( n_neg + n_inf + n_nan ) {
 		ERROR("WARNING: The GPU calculation produced %i negative"
@@ -332,6 +361,25 @@ void 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++ ) {
+
+		struct rvec q;
+		double twotheta, k;
+		int idx;
+
+		/* Calculate k this time round */
+		k = 1.0/image->lambda;
+
+		q = get_q(image, fs, ss, &twotheta, k);
+
+		idx = fs + image->width*ss;
+		image->twotheta[idx] = twotheta;
+
+	}
+	}
 }
 
 
diff --git a/src/diffraction.c b/src/diffraction.c
index 870569b7..b76154c2 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -420,7 +420,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 		image->data[idx] = I_lattice * I_molecule;
 		image->twotheta[idx] = twotheta;
 
-
 	}
 	progress_bar(fs, image->width-1, "Calculating diffraction");
 	}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 8061b1ae..2c1af284 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -291,9 +291,10 @@ static void process_image(const struct index_args *iargs,
 	image.indexed_cell = NULL;
 	image.det = copy_geom(iargs->det);
 	image.copyme = iargs->copyme;
-	image.beam = beam;
+	image.reflections = NULL;
 	image.id = cookie;
 	image.filename = pargs->filename;
+	image.beam = beam;
 
 	hdfile = hdfile_open(image.filename);
 	if ( hdfile == NULL ) return;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 03b009b4..e05dcdec 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -497,6 +497,8 @@ int main(int argc, char *argv[])
 	image.width = image.det->max_fs + 1;
 	image.height = image.det->max_ss + 1;
 	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+	image.bw = image.beam->bandwidth;
+	image.div = image.beam->divergence;
 
 	/* Load unit cell */
 	input_cell = load_cell_from_pdb(filename);
@@ -617,8 +619,7 @@ int main(int argc, char *argv[])
 			number++;
 
 			/* Write the output file */
-			hdf5_write(filename, image.data,
-			           image.width, image.height, H5T_NATIVE_FLOAT);
+			hdf5_write_image(filename, &image);
 
 		}
 
diff --git a/tests/integration_check.c b/tests/integration_check.c
index 2d18ac5a..d43e5fd2 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -64,7 +64,7 @@ static void third_integration_check(struct image *image, int n_trials,
 		}
 
 		r = integrate_peak(image, 64, 64, &fsp, &ssp,
-		                   &intensity, &sigma, 10.0, 15.0, 17.0);
+		                   &intensity, &sigma, 10.0, 15.0, 17.0, 0);
 
 		if ( r == 0 ) {
 			mean_intensity += intensity;
@@ -125,7 +125,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
 		}
 
 		r = integrate_peak(image, 64, 64, &fsp, &ssp,
-		                   &intensity, &sigma, 10.0, 15.0, 17.0);
+		                   &intensity, &sigma, 10.0, 15.0, 17.0, 0);
 
 		if ( r == 0 ) {
 			mean_intensity += intensity;
@@ -203,9 +203,11 @@ int main(int argc, char *argv[])
 	image.height = 128;
 	memset(image.data, 0, 128*128*sizeof(float));
 
+	image.reflections=reflist_new();
+
 	/* First check: no intensity -> no peak, or very low intensity */
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0);
+	                   10.0, 15.0, 17.0, 0);
 	STATUS("  First check: integrate_peak() returned %i", r);
 	if ( r == 0 ) {
 
@@ -231,7 +233,7 @@ int main(int argc, char *argv[])
 	}
 
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0);
+	                   10.0, 15.0, 17.0, 0);
 	if ( r ) {
 		ERROR(" Second check: integrate_peak() returned %i (wrong).\n",
 		      r);
@@ -273,7 +275,7 @@ int main(int argc, char *argv[])
 	}
 
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0);
+	                   10.0, 15.0, 17.0, 0);
 	if ( r ) {
 		ERROR("   Fifth check: integrate_peak() returned %i (wrong).\n",
 		      r);