Add Poisson noise to image

3 files changed, 25 insertions, 0 deletions

diff --git a/src/detector.c b/src/detector.c
index 1a7d36ee..fae37126 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -177,6 +177,8 @@ void record_image(struct image *image)
 		counts = intensity * ph_per_e * sa;
 		//printf("%e counts\n", counts);
 
+		counts += poisson_noise(counts);
+
 		image->hdr[x + image->width*y] = counts;
 
 	}
diff --git a/src/utils.c b/src/utils.c
index d5e1b874..faf78502 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -71,3 +71,25 @@ void progress_bar(int val, int total)
 
 	fflush(stdout);
 }
+
+
+double poisson_noise(double expected)
+{
+	double L;
+	int k;
+	double p = 1.0;
+
+	L = exp(-expected);
+
+	do {
+
+		double r;
+
+		k++;
+		r = (double)random()/RAND_MAX;
+		p *= r;
+
+	} while ( p > L );
+
+	return k - 1;
+}
diff --git a/src/utils.h b/src/utils.h
index 9e6458cf..54cfdcfb 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -55,6 +55,7 @@ extern double angle_between(double x1, double y1, double z1,
 extern size_t skipspace(const char *s);
 extern void chomp(char *s);
 extern void progress_bar(int val, int total);
+extern double poisson_noise(double expected);
 
 /* Keep these ones inline, to avoid function call overhead */
 static inline double distance(double x1, double y1, double x2, double y2)