Fix integrate_peak() problems which caused incorrect peak positions

5 files changed, 30 insertions, 27 deletions

diff --git a/src/facetron.c b/src/facetron.c
index 38d701e6..f16de2f7 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -207,7 +207,8 @@ static void integrate_image(struct image *image)
 		 * pattern? */
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[j].x, spots[j].y,
-		                    &xc, &yc, &i_partial, NULL, NULL, 1, 1) ) {
+		                    &xc, &yc, &i_partial, NULL, NULL,
+		                    1, 1, 0) ) {
 			continue;
 		}
 
diff --git a/src/geometry.c b/src/geometry.c
index 263be4af..f3146fb1 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -284,7 +284,7 @@ double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
 		float x, y, intensity;
 
 		if ( integrate_peak(image, cpeaks[i].x, cpeaks[i].y, &x, &y,
-                                    &intensity, NULL, NULL, 0, 0) ) continue;
+                                    &intensity, NULL, NULL, 0, 0, 0) ) continue;
 
 		itot += intensity;
 	}
diff --git a/src/peaks.c b/src/peaks.c
index 1d548aa0..263040dd 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -49,12 +49,6 @@
 /* Window size for Zaefferer peak detection */
 #define PEAK_WINDOW_SIZE (10)
 
-/* Integration radius for peaks */
-#define INTEGRATION_RADIUS (10)
-
-/* Integration radius for background estimation */
-#define OUT_INTEGRATION_RADIUS (15)
-
 
 static int in_streak(int x, int y)
 {
@@ -188,26 +182,33 @@ static int cull_peaks(struct image *image)
 int integrate_peak(struct image *image, int xp, int yp,
                    float *xc, float *yc, float *intensity,
                    double *pbg, double *pmax,
-                   int do_polar, int do_sa)
+                   int do_polar, int do_sa, int centroid)
 {
 	signed int x, y;
-	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
-	const int out_lim = OUT_INTEGRATION_RADIUS * OUT_INTEGRATION_RADIUS;
+	int lim, out_lim;
 	double total = 0.0;
 	int xct = 0;
 	int yct = 0;
 	double noise = 0.0;
 	int noise_counts = 0;
 	double max = 0.0;
+	struct panel *p = NULL;
+
+	p = find_panel(image->det, xp, yp);
+	if ( p == NULL ) return 1;
+	if ( p->no_index ) return 1;
+
+	lim = p->peak_sep/2;
+	out_lim = lim + 1;
 
-	for ( x=-OUT_INTEGRATION_RADIUS; x<+OUT_INTEGRATION_RADIUS; x++ ) {
-	for ( y=-OUT_INTEGRATION_RADIUS; y<+OUT_INTEGRATION_RADIUS; y++ ) {
+	for ( x=-out_lim; x<+out_lim; x++ ) {
+	for ( y=-out_lim; y<+out_lim; y++ ) {
 
-		struct panel *p = NULL;
 		double val, sa, pix_area, Lsq, dsq, proj_area;
 		float tt = 0.0;
 		double phi, pa, pb, pol;
 		uint16_t flags;
+		struct panel *p2;
 
 		/* Outer mask radius */
 		if ( x*x + y*y > out_lim ) continue;
@@ -215,6 +216,10 @@ int integrate_peak(struct image *image, int xp, int yp,
 		if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue;
 		if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue;
 
+		/* Strayed off one panel? */
+		p2 = find_panel(image->det, x+xp, y+yp);
+		if ( p2 != p ) return 1;
+
 		/* Veto this peak if we tried to integrate in a bad region */
 		if ( image->flags != NULL ) {
 			flags = image->flags[(x+xp)+image->width*(y+yp)];
@@ -233,11 +238,6 @@ int integrate_peak(struct image *image, int xp, int yp,
 
 		if ( val > max ) max = val;
 
-		p = find_panel(image->det, x+xp, y+yp);
-		if ( p == NULL ) return 1;
-
-		if ( p->no_index ) return 1;
-
 		if ( do_sa ) {
 
 			/* Area of one pixel */
@@ -281,7 +281,7 @@ int integrate_peak(struct image *image, int xp, int yp,
 	}
 
 	/* The centroid is excitingly undefined if there is no intensity */
-	if ( total != 0 ) {
+	if ( centroid && (total != 0) ) {
 		*xc = (float)xct / total;
 		*yc = (float)yct / total;
 		*intensity = total;
@@ -420,7 +420,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
 		 * Don't bother doing polarisation/SA correction, because the
 		 * intensity of this peak is only an estimate at this stage. */
 		r = integrate_peak(image, mask_x, mask_y,
-		                   &fx, &fy, &intensity, NULL, NULL, 0, 0);
+		                   &fx, &fy, &intensity, NULL, NULL, 0, 0, 1);
 		if ( r ) {
 			/* Bad region - don't detect peak */
 			nrej_bad++;
@@ -766,7 +766,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 				 * revised coordinates. */
 				r = integrate_peak(image, f->x, f->y, &x, &y,
 					           &intensity, &bg, &max,
-					           polar, sa);
+					           polar, sa, 1);
 				if ( r ) {
 					/* The original peak (which also went
 					 * through integrate_peak(), but with
@@ -788,7 +788,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 					           image->cpeaks[i].x,
 					           image->cpeaks[i].y,
 					           &x, &y, &intensity, &bg, &max,
-					           polar, sa);
+					           polar, sa, 1);
 				if ( r ) {
 					/* Plain old ordinary peak veto */
 					n_veto++;
@@ -809,7 +809,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 				           image->cpeaks[i].x,
 				           image->cpeaks[i].y,
 				           &x, &y, &intensity, &bg, &max, polar,
-				           sa);
+				           sa, 0);
 			if ( r ) {
 				/* Plain old ordinary peak veto */
 				n_veto++;
diff --git a/src/peaks.h b/src/peaks.h
index cfc3d322..ba1eea30 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -39,6 +39,6 @@ extern int find_projected_peaks(struct image *image, UnitCell *cell,
 extern int integrate_peak(struct image *image, int xp, int yp,
                           float *xc, float *yc, float *intensity,
                           double *pbg, double *pmax,
-                          int do_polar, int do_sa);
+                          int do_polar, int do_sa, int centroid);
 
 #endif	/* PEAKS_H */
diff --git a/src/post-refinement.c b/src/post-refinement.c
index a9382773..d835a2ff 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -245,7 +245,8 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
 		 * pattern? */
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[h].x, spots[h].y,
-		                    &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
+		                    &xc, &yc, &I_partial, NULL, NULL,
+		                    1, 1, 0) ) {
 			continue;
 		}
 
@@ -316,7 +317,8 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 		 * pattern? */
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[h].x, spots[h].y,
-		                    &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
+		                    &xc, &yc, &I_partial, NULL, NULL,
+		                    1, 1, 0) ) {
 			continue;
 		}