Attempt to make the peak integration suck (slightly) less

3 files changed, 89 insertions, 99 deletions

diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 66907298..e1a828e5 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -151,23 +151,20 @@ static int cull_peaks(struct image *image)
 /* 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 *pbg, double *pmax, double *sigma,
-                   int centroid, int bgsub)
+                   double *pfs, double *pss, double *intensity, double *sigma)
 {
 	signed int fs, ss;
 	double lim, out_lim, mid_lim;
 	double lim_sq, out_lim_sq, mid_lim_sq;
-	double total = 0.0;
-	double fsct = 0.0;
-	double ssct = 0.0;
-	double noise = 0.0;
-	int noise_counts = 0;
-	double max = 0.0;
-	struct panel *p = NULL;
-        int pixel_counts = 0;
-        double noise_mean = 0.0;
-	double noise_meansq = 0.0;
+	double pk_total;
+	int pk_counts;
+	double fsct, ssct;
+	double bg_tot = 0.0;
+	int bg_counts = 0;
+	struct panel *p;
+	double bg_mean, bg_var;
+	double bg_tot_sq = 0.0;
+	double var;
 	struct beam_params *beam;
 	double aduph;
 
@@ -189,6 +186,7 @@ int integrate_peak(struct image *image, int cfs, int css,
 	mid_lim_sq = pow(mid_lim, 2.0);
 	out_lim_sq = pow(out_lim, 2.0);
 
+	/* Estimate the background */
 	for ( fs=-out_lim; fs<+out_lim; fs++ ) {
 	for ( ss=-out_lim; ss<+out_lim; ss++ ) {
 
@@ -199,11 +197,9 @@ int integrate_peak(struct image *image, int cfs, int css,
 		struct panel *p2;
 		int idx;
 
-		/* Outer mask radius */
+		/* Restrict to annulus */
 		if ( fs*fs + ss*ss > out_lim_sq ) continue;
-
-		if ( ((fs+cfs)>=image->width) || ((fs+cfs)<0) ) continue;
-		if ( ((ss+css)>=image->height) || ((ss+css)<0) ) continue;
+		if ( fs*fs + ss*ss < mid_lim_sq ) continue;
 
 		/* Strayed off one panel? */
 		p2 = find_panel(image->det, fs+cfs, ss+css);
@@ -227,65 +223,80 @@ int integrate_peak(struct image *image, int cfs, int css,
 
 		val = image->data[idx];
 
-		if ( val > max ) max = val;
+		bg_tot += val;
+		bg_tot_sq += pow(val, 2.0);
+		bg_counts++;
 
-		/* If outside inner mask, estimate noise from this region */
-		if ( fs*fs + ss*ss > mid_lim_sq ) {
+	}
+	}
 
-			/* Noise
-			 * noise and noise_meansq are both in photons (^2) */
-			noise += val / aduph;
-			noise_counts++;
-			noise_meansq += pow(val/aduph, 2.0);
+	if ( bg_counts == 0 ) return 1;
+	bg_mean = bg_tot / bg_counts;
+	bg_var = (bg_tot_sq/bg_counts) - pow(bg_mean, 2.0);
 
-		} else if ( fs*fs + ss*ss < lim_sq ) {
+	/* Measure the peak */
+	pk_total = 0.0;
+	pk_counts = 0;
+	fsct = 0.0;  ssct = 0.0;
+	for ( fs=-lim; fs<+lim; fs++ ) {
+	for ( ss=-lim; ss<+lim; ss++ ) {
 
-			/* Peak */
-			pixel_counts++;
-			total += val;
-			fsct += val*(cfs+fs);
-			ssct += val*(css+ss);
+		double val;
+		double tt = 0.0;
+		double phi, pa, pb, pol;
+		uint16_t flags;
+		struct panel *p2;
+		int idx;
 
-		}
+		/* Inner mask radius */
+		if ( fs*fs + ss*ss > lim_sq ) continue;
 
-	}
-	}
+		/* Strayed off one panel? */
+		p2 = find_panel(image->det, fs+cfs, ss+css);
+		if ( p2 != p ) return 1;
 
-	noise_mean = noise / noise_counts;  /* photons */
+		idx = fs+cfs+image->width*(ss+css);
 
-	/* The centroid is excitingly undefined if there is no intensity */
-	if ( centroid && (total != 0) ) {
-		*pfs = ((double)fsct / total) + 0.5;
-		*pss = ((double)ssct / total) + 0.5;
-	} else {
-		*pfs = (double)cfs + 0.5;
-		*pss = (double)css + 0.5;
-	}
-	if ( bgsub ) {
-		*intensity = total - aduph * pixel_counts*noise_mean; /* ADU */
-	} else {
-		*intensity = total;  /* ADU */
-	}
+		/* Veto this peak if we tried to integrate in a bad region */
+		if ( image->flags != NULL ) {
 
-	if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
+			flags = image->flags[idx];
 
-	if ( sigma != NULL ) {
+			/* It must have all the "good" bits to be valid */
+			if ( !((flags & image->det->mask_good)
+			                   == image->det->mask_good) ) return 1;
 
-		/* First term is standard deviation of background per pixel
-		 * sqrt(pixel_counts) - increase of error for integrated value
-		 * sqrt(2) - increase of error for background subtraction  */
-		*sigma = sqrt(noise_meansq/noise_counts-(noise_mean*noise_mean))
-		          * sqrt(2.0*pixel_counts) * aduph;
+			/* If it has any of the "bad" bits, reject */
+			if ( flags & image->det->mask_bad ) return 1;
 
-	}
+		}
+
+		val = image->data[idx] - bg_mean;
+
+		/* At least half a photon? */
+		if ( val < aduph / 2.0 ) continue;
+
+		pk_counts++;
+		pk_total += val;
+
+		fsct += val*(cfs+fs);
+		ssct += val*(css+ss);
 
-	if ( pbg != NULL ) {
-		*pbg = aduph * (noise / noise_counts);
 	}
-	if ( pmax != NULL ) {
-		*pmax = max;
 	}
 
+	if ( pk_counts == 0 ) return 1;
+
+	*pfs = ((double)fsct / pk_total) + 0.5;
+	*pss = ((double)ssct / pk_total) + 0.5;
+
+	var = pk_counts * bg_var;
+	var += aduph * pk_total;
+	if ( var < 0.0 ) return 1;
+
+	if ( intensity != NULL ) *intensity = pk_total;
+	if ( sigma != NULL ) *sigma = sqrt(var);
+
 	return 0;
 }
 
@@ -396,12 +407,9 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		assert(mask_fs >= p->min_fs);
 		assert(mask_ss >= p->min_ss);
 
-		/* Centroid peak and get better coordinates.
-		 * Don't bother doing polarisation/SA correction, because the
-		 * intensity of this peak is only an estimate at this stage. */
+		/* Centroid peak and get better coordinates. */
 		r = integrate_peak(image, mask_fs, mask_ss,
-		                   &f_fs, &f_ss, &intensity,
-		                   &pbg, &pmax, &sigma, 1, 1);
+		                   &f_fs, &f_ss, &intensity, &sigma);
 
 		if ( r ) {
 			/* Bad region - don't detect peak */
@@ -654,7 +662,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 		}
 
 		r = integrate_peak(image, pfs, pss, &fs, &ss,
-		                   &intensity, &bg, &max, &sigma, 0, bgsub);
+		                   &intensity, &sigma);
 
 		/* Record intensity and set redundancy to 1 on success */
 		if ( r == 0 ) {
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 8b64a081..e5e31b91 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -45,10 +45,10 @@ extern double peak_lattice_agreement(struct image *image, UnitCell *cell,
 extern int peak_sanity_check(struct image *image);
 
 /* Exported so it can be poked by integration_check */
-extern int integrate_peak(struct image *image, int cfs, int css,
-                          double *pfs, double *pss, double *intensity,
-                          double *pbg, double *pmax, double *sigma,
-                          int centroid, int bgsub);
+extern int integrate_peak(struct image *image,
+                          int cfs, int css,
+                          double *pfs, double *pss,
+                          double *intensity, double *sigma);
 
 extern void estimate_resolution(RefList *list, UnitCell *cell,
                                 double *min, double *max);
diff --git a/tests/integration_check.c b/tests/integration_check.c
index aa9e2705..1e439611 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -59,8 +59,7 @@ static void third_integration_check(struct image *image, int n_trials,
 			image->data[fs+image->width*ss] = poisson_noise(10.0);
 		}
 		}
-		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
-			       &bg, &max, &sigma, 1, 1);
+		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
 
 		mean_intensity += intensity;
 		mean_bg += bg;
@@ -122,8 +121,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
 			}
 		}
 		}
-		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
-			       &bg, &max, &sigma, 1, 1);
+		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
 
 		mean_intensity += intensity;
 		mean_bg += bg;
@@ -182,8 +180,7 @@ static void fifth_integration_check(struct image *image, int n_trials,
 			}
 		}
 		}
-		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
-			       &bg, &max, &sigma, 1, 0);
+		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
 
 		mean_intensity += intensity;
 		mean_bg += bg;
@@ -244,8 +241,7 @@ static void sixth_integration_check(struct image *image, int n_trials,
 			}
 		}
 		}
-		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity,
-			       &bg, &max, &sigma, 1, 0);
+		integrate_peak(image, 64, 64, &fsp, &ssp, &intensity, &sigma);
 
 		mean_intensity += intensity;
 		mean_bg += bg;
@@ -281,7 +277,7 @@ static void sixth_integration_check(struct image *image, int n_trials,
 int main(int argc, char *argv[])
 {
 	struct image image;
-	double fsp, ssp, intensity, bg, max, sigma;
+	double fsp, ssp, intensity, sigma;
 	int fs, ss;
 	FILE *fh;
 	unsigned int seed;
@@ -326,18 +322,13 @@ int main(int argc, char *argv[])
 	image.beam->adu_per_photon = 100.0;
 
 	/* First check: no intensity -> zero intensity and bg */
-	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity,
-	               &bg, &max, &sigma, 1, 1);
-	STATUS("  First check: intensity = %.2f, bg = %.2f, max = %.2f,"
-	       " sigma = %.2f\n", intensity, bg, max, sigma);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma);
+	STATUS("  First check: intensity = %.2f, sigma = %.2f\n",
+	       intensity, sigma);
 	if ( intensity != 0.0 ) {
 		ERROR("Intensity should be zero.\n");
 		fail = 1;
 	}
-	if ( bg != 0.0 ) {
-		ERROR("Background should be zero.\n");
-		fail = 1;
-	}
 
 	/* Second check: uniform peak gives correct value and no bg */
 	for ( fs=0; fs<image.width; fs++ ) {
@@ -346,22 +337,13 @@ int main(int argc, char *argv[])
 		image.data[fs+image.width*ss] = 1000.0;
 	}
 	}
-	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity,
-	               &bg, &max, &sigma, 1, 1);
-	STATUS(" Second check: intensity = %.2f, bg = %.2f, max = %.2f,"
-	       " sigma = %.2f\n", intensity, bg, max, sigma);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma);
+	STATUS(" Second check: intensity = %.2f, sigma = %.2f\n",
+	       intensity, sigma);
 	if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
 		ERROR("Intensity should be close to 1000*pi*integr_r^2\n");
 		fail = 1;
 	}
-	if ( bg != 0.0 ) {
-		ERROR("Background should be zero.\n");
-		fail = 1;
-	}
-	if ( max != 1000.0 ) {
-		ERROR("Max should be 1000.\n");
-		fail = 1;
-	}
 	if ( sigma != 0.0 ) {
 		ERROR("Sigma should be zero.\n");
 		fail = 1;