aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2009-11-18 15:13:09 +0100
committerThomas White <taw@bitwiz.org.uk>2009-11-18 15:13:09 +0100
commitc77d160603298560b89243dda34fd81ee52d3d10 (patch)
tree1708513688c4dad90fec576480f420a6f77917dd
parent239dfda559a4982079b3e7e824eeec2e3ebf3696 (diff)
Bloom simulation
-rw-r--r--src/detector.c103
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;
}