Added background subtraction and sigma

5 files changed, 68 insertions, 30 deletions

diff --git a/src/geometry.c b/src/geometry.c
index f58f21e1..1e61c26d 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -285,7 +285,7 @@ double integrate_all(struct image *image, RefList *reflections)
 		get_detector_pos(refl, &xp, &yp);
 
 		if ( integrate_peak(image, xp, yp, &x, &y,
-                                    &intensity, NULL, NULL, 0, 0) ) continue;
+                                    &intensity, NULL, NULL, NULL, 0, 0) ) continue;
 
 		itot += intensity;
 	}
diff --git a/src/peaks.c b/src/peaks.c
index 0af195e6..aa56da8b 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -4,6 +4,7 @@
  * Peak search and other image analysis
  *
  * (c) 2006-2011 Thomas White <taw@physics.org>
+ *	    2011 Andrew Martin
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -142,7 +143,7 @@ static int cull_peaks(struct image *image)
  * 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 *pbg, double *pmax, double *sigma,
                    int do_polar, int centroid)
 {
 	signed int fs, ss;
@@ -155,6 +156,10 @@ int integrate_peak(struct image *image, int cfs, int css,
 	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;
+
 
 	p = find_panel(image->det, cfs, css);
 	if ( p == NULL ) return 1;
@@ -203,16 +208,6 @@ int integrate_peak(struct image *image, int cfs, int css,
 
 		val = image->data[idx];
 
-		/* Inner mask */
-		if ( fs*fs + ss*ss > lim_sq ) {
-			/* Estimate noise from this region */
-			noise += fabs(val);
-			noise_counts++;
-			continue;
-		}
-
-		if ( val > max ) max = val;
-
 		if ( do_polar ) {
 
 			tt = get_tt(image, fs+cfs, ss+css);
@@ -226,27 +221,54 @@ int integrate_peak(struct image *image, int cfs, int css,
 
 		}
 
-		total += val;
+		if ( val > max ) max = val;
 
+		/* Inner mask */
+		if ( fs*fs + ss*ss > lim_sq ) {
+			/* Estimate noise from this region */
+			noise += val; //fabs(val);
+			noise_counts++;
+			/* Estimate the standard deviation of the noise */
+			noise_meansq += fabs(val)*fabs(val) ;
+			continue;
+		}
+
+		pixel_counts ++;
+		total += val;
 		fsct += val*(cfs+fs);
 		ssct += val*(css+ss);
 
 	}
 	}
 
+        noise_mean = noise / noise_counts;
+
 	/* The centroid is excitingly undefined if there is no intensity */
 	if ( centroid && (total != 0) ) {
 		*pfs = (double)fsct / total;
 		*pss = (double)ssct / total;
-		*intensity = total;
+		*intensity = total - pixel_counts*noise_mean;
 	} else {
 		*pfs = (double)cfs;
 		*pss = (double)css;
-		*intensity = total;
+		*intensity = total - pixel_counts*noise_mean;
 	}
 
 	if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
 
+        if ( sigma != NULL ) {
+		/*
+		 * 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);
+		/* printf(" counts %d %d \n", noise_counts, pixel_counts);
+		   printf(" intensity, bg, diff, %f, %f, %f \n", total, pixel_counts*noise_mean, total - pixel_counts*noise_mean);
+		   printf(" sigma = %f \n", *sigma);	*/
+	}
+
 	if ( pbg != NULL ) {
 		*pbg = (noise / noise_counts);
 	}
@@ -362,7 +384,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		 * 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_fs, mask_ss,
-		                   &f_fs, &f_ss, &intensity, NULL, NULL, 0, 1);
+		                   &f_fs, &f_ss, &intensity, NULL, NULL, NULL, 0, 1);
 		if ( r ) {
 			/* Bad region - don't detect peak */
 			nrej_bad++;
@@ -587,6 +609,7 @@ void integrate_reflections(struct image *image, int polar, int use_closer)
 		double d;
 		int idx;
 		double bg, max;
+		double sigma;
 		struct panel *p;
 		double pfs, pss;
 		int r;
@@ -616,11 +639,12 @@ void integrate_reflections(struct image *image, int polar, int use_closer)
 		}
 
 		r = integrate_peak(image, pfs, pss, &fs, &ss,
-		                   &intensity, &bg, &max, polar, 0);
+		                   &intensity, &bg, &max, &sigma, polar, 0);
 
 		/* Record intensity and set redundancy to 1 on success */
 		if ( r == 0 ) {
 			set_int(refl, intensity);
+			set_esd_intensity(refl,sigma);
 			set_redundancy(refl, 1);
 		}
 
diff --git a/src/peaks.h b/src/peaks.h
index 7a6557f1..2faa2b5e 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -35,7 +35,7 @@ extern RefList *find_projected_peaks(struct image *image, UnitCell *cell,
 
 extern int integrate_peak(struct image *image, int xp, int yp,
                           double *xc, double *yc, double *intensity,
-                          double *pbg, double *pmax,
+                          double *pbg, double *pmax, double *sigma,
                           int do_polar, int centroid);
 
 #endif	/* PEAKS_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index e5a1f4d0..10597337 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -4,7 +4,7 @@
  * Assemble and process FEL Bragg intensities
  *
  * (c) 2006-2011 Thomas White <taw@physics.org>
- *
+ *  	    2011 Andrew Martin
  * Part of CrystFEL - crystallography with a FEL
  *
  */
@@ -171,8 +171,12 @@ static void merge_pattern(RefList *model, RefList *new, int max_only,
 		} else if ( pass == 2 ) {
 
 			double dev = get_sum_squared_dev(model_version);
-			set_sum_squared_dev(model_version,
-			                   dev + pow(intensity-model_int, 2.0));
+			double refl_esd = get_esd_intensity(refl);
+			
+			set_sum_squared_dev(model_version, 
+					      	  	// dev + pow(refl_esd,2.0) ); 
+						        dev    + pow(intensity, 2.0) ); 
+				                  // dev + pow(intensity-model_int, 2.0));			
 
 			if ( hist_vals != NULL ) {
 				int p = *hist_n;
@@ -188,6 +192,8 @@ static void merge_pattern(RefList *model, RefList *new, int max_only,
 		}
 
 	}
+
+
 }
 
 
@@ -365,11 +371,19 @@ static void merge_all(FILE *fh, RefList *model,
 		      refl = next_refl(refl, iter) ) {
 
 			double sum_squared_dev = get_sum_squared_dev(refl);
+			double intensity = get_intensity(refl);
 			int red = get_redundancy(refl);
-
-			set_esd_intensity(refl,
-			                  sqrt(sum_squared_dev)/(double)red);
-
+			int h, k, l;
+			get_indices(refl,&h,&k,&l);
+	 
+			/*
+			 set_esd_intensity(refl,
+			         sqrt(sum_squared_dev)/(double)red); 
+			 */
+								
+			 set_esd_intensity(refl,
+			        sqrt((sum_squared_dev/(double)red) - pow(intensity,2.0) )/(double)red);
+			
 		}
 	}
 
diff --git a/tests/integration_check.c b/tests/integration_check.c
index f80181ee..1e6e64db 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -27,7 +27,7 @@
 int main(int argc, char *argv[])
 {
 	struct image image;
-	double fsp, ssp, intensity, bg, max;
+	double fsp, ssp, intensity, bg, max, sigma;
 	int fs, ss;
 
 	image.data = malloc(128*128*sizeof(float));
@@ -63,7 +63,7 @@ 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, 0, 1);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
 	STATUS(" First check: intensity = %.2f bg = %.2f max = %.2f\n",
 	       intensity, bg, max);
 	if ( intensity != 0.0 ) {
@@ -82,7 +82,7 @@ int main(int argc, char *argv[])
 		image.data[fs+image.width*ss] = 1000.0;
 	}
 	}
-	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
 	STATUS("Second check: intensity = %.2f, bg = %.2f, max = %.2f\n",
 	       intensity, bg, max);
 	if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {
@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
 		image.data[fs+image.width*ss] = poisson_noise(10.0) - 10.0;
 	}
 	}
-	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
 	STATUS(" Third check: intensity = %.2f, bg = %.2f, max = %.2f\n",
 	       intensity, bg, max);
 	if ( fabs(intensity) > 100.0 ) {
@@ -122,7 +122,7 @@ int main(int argc, char *argv[])
 		image.data[fs+image.width*ss] = 1000.0;
 	}
 	}
-	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, 0, 1);
+	integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &bg, &max, &sigma, 0, 1);
 	STATUS("Fourth check: intensity = %.2f, bg = %.2f, max = %.2f\n",
 	       intensity, bg, max);
 	if ( fabs(intensity - M_PI*9.0*9.0*1000.0) > 4000.0 ) {