Bloom simulation

1 files changed, 91 insertions, 12 deletions

diff --git a/src/detector.c b/src/detector.c
index 13d5a20b..9de99e6c 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -13,6 +13,7 @@
 #include <stdlib.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 
 #include "image.h"
 #include "utils.h"
@@ -23,28 +24,106 @@
 
 
 /* Detector's quantum efficiency */
-#define DQE (0.8)
+#define DQE (0.9)
 
 
-static uint16_t *bloom(double *hdr, int width, int height)
+/* Detector's saturation value */
+#define SATURATION (60000)
+
+
+/* Bleed excess intensity into neighbouring pixels */
+static void bloom_values(double *tmp, int x, int y,
+                          int width, int height, double val)
+{
+	double overflow;
+
+	overflow = val - SATURATION;
+
+	/* Intensity which bleeds off the edge of the detector is lost */
+	if ( x > 0 ) {
+		tmp[x-1 + width*y] += overflow / 6.0;
+		if ( y > 0 ) {
+			tmp[x-1 + width*(y-1)] += overflow / 12.0;
+		}
+		if ( y < height-1 ) {
+			tmp[x-1 + width*(y+1)] += overflow / 12.0;
+		}
+	}
+
+	if ( x < width-1 ) {
+		tmp[x+1 + width*y] += overflow / 6.0;
+		if ( y > 0 ) {
+			tmp[x+1 + width*(y-1)] += overflow / 12.0;
+		}
+		if ( y < height-1 ) {
+			tmp[x+1 + width*(y+1)] += overflow / 12.0;
+		}
+	}
+
+	if ( y > 0 ) {
+		tmp[x + width*(y-1)] += overflow / 6.0;
+	}
+	if ( y < height-1 ) {
+		tmp[x + width*(y+1)] += overflow / 6.0;
+	}
+}
+
+
+static uint16_t *bloom(double *hdr_in, int width, int height)
 {
 	int x, y;
 	uint16_t *data;
-	
+	double *tmp;
+	double *hdr;
+	int did_something;
+
 	data = malloc(width * height * sizeof(uint16_t));
-	
+	tmp = malloc(width * height * sizeof(double));
+	hdr = malloc(width * height * sizeof(double));
+
+	memcpy(hdr, hdr_in, width*height*sizeof(double));
+	do {
+
+		memset(tmp, 0, width*height*sizeof(double));
+		did_something = 0;
+
+		for ( x=0; x<width; x++ ) {
+		for ( y=0; y<height; y++ ) {
+
+			double hdval;
+
+			hdval = hdr[x + width*y];
+
+			/* Pixel not saturated? */
+			if ( hdval <= SATURATION ) {
+				tmp[x + width*y] += hdval;
+				continue;
+			}
+
+			bloom_values(tmp, x, y, width, height, hdval);
+			tmp[x + width*y] += SATURATION;
+			did_something = 1;
+
+		}
+		}
+
+		/* Prepare new input if we're going round for another pass */
+		if ( did_something ) {
+			memcpy(hdr, tmp, width*height*sizeof(double));
+		}
+
+	} while ( did_something );
+
+	/* Turn into integer array of counts */
 	for ( x=0; x<width; x++ ) {
 	for ( y=0; y<height; y++ ) {
-	
-		double hdval;
-		
-		hdval = hdr[x + width*y] * DQE;
-		
-		
-	
+		data[x + width*y] = tmp[x + width*y];
 	}
 	}
-	
+
+	free(tmp);
+	free(hdr);
+
 	return data;
 }