aboutsummaryrefslogtreecommitdiff
path: root/src/detector.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-02-26 10:58:13 +0100
committerThomas White <taw@physics.org>2010-02-26 13:40:11 +0100
commit86dd71e8640394f4e4f5aa71b2e5f51f5fea4a11 (patch)
treee7039c5d4e9b5477bc1e350e0330d157125d5c9f /src/detector.c
parentb4664bc2a399463bd4e01f331a153749b3947b8f (diff)
Handle images as floats rather than int16_t
Also, remove bloom - it's not a useful model for what really happens and takes too long (this isn't a detector simulation..)
Diffstat (limited to 'src/detector.c')
-rw-r--r--src/detector.c135
1 files changed, 10 insertions, 125 deletions
diff --git a/src/detector.c b/src/detector.c
index 80a88a1a..d605822b 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -24,11 +24,8 @@
/* Number of photons in pulse */
#define FLUENCE (1.0e13)
-/* Detector's quantum efficiency */
-#define DQE (0.9)
-
-/* Detector's saturation value */
-#define SATURATION (5000)
+/* Detector's quantum efficiency (ADU per photon, front detector) */
+#define DQE (167.0)
/* Radius of the water column */
#define WATER_RADIUS (3.0e-6 / 2.0)
@@ -37,113 +34,7 @@
#define BEAM_RADIUS (3.0e-6 / 2.0)
-/* Bleed excess intensity into neighbouring pixels */
-static void bloom_values(int *tmp, int x, int y,
- int width, int height, int val)
-{
- int 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;
- if ( y > 0 ) {
- tmp[x-1 + width*(y-1)] += overflow / 12;
- }
- if ( y < height-1 ) {
- tmp[x-1 + width*(y+1)] += overflow / 12;
- }
- }
-
- if ( x < width-1 ) {
- tmp[x+1 + width*y] += overflow / 6;
- if ( y > 0 ) {
- tmp[x+1 + width*(y-1)] += overflow / 12;
- }
- if ( y < height-1 ) {
- tmp[x+1 + width*(y+1)] += overflow / 12;
- }
- }
-
- if ( y > 0 ) {
- tmp[x + width*(y-1)] += overflow / 6;
- }
- if ( y < height-1 ) {
- tmp[x + width*(y+1)] += overflow / 6;
- }
-}
-
-
-static int16_t *bloom(int *hdr_in, int width, int height)
-{
- int x, y;
- int16_t *data;
- int *tmp;
- int *hdr;
- int did_something;
-
- data = malloc(width * height * sizeof(int16_t));
- tmp = malloc(width * height * sizeof(int));
- hdr = malloc(width * height * sizeof(int));
-
- memcpy(hdr, hdr_in, width*height*sizeof(int));
-
- /* Apply DQE (once only) */
- for ( x=0; x<width; x++ ) {
- for ( y=0; y<height; y++ ) {
- hdr[x + width*y] *= DQE;
- }
- }
-
- do {
-
- memset(tmp, 0, width*height*sizeof(int));
- did_something = 0;
-
- for ( x=0; x<width; x++ ) {
- for ( y=0; y<height; y++ ) {
-
- int 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(int));
- }
-
- } while ( did_something );
-
- /* Turn into integer array of counts */
- for ( x=0; x<width; x++ ) {
- for ( y=0; y<height; y++ ) {
- data[x + width*y] = (int16_t)tmp[x + width*y];
- }
- }
-
- free(tmp);
- free(hdr);
-
- return data;
-}
-
-
-void record_image(struct image *image, int do_water, int do_poisson,
- int do_bloom)
+void record_image(struct image *image, int do_water, int do_poisson)
{
int x, y;
double total_energy, energy_density;
@@ -222,19 +113,13 @@ void record_image(struct image *image, int do_water, int do_poisson,
progress_bar(x, image->width-1, "Post-processing");
}
- if ( do_bloom ) {
- image->data = bloom(image->hdr, image->width, image->height);
- } else {
- image->data = malloc(image->width * image->height
- * sizeof(int16_t));
- for ( x=0; x<image->width; x++ ) {
- for ( y=0; y<image->height; y++ ) {
- int val;
- val = image->hdr[x + image->width*y];
- if ( val > SATURATION ) val = SATURATION;
- image->data[x + image->width*y] = (int16_t)val;
- }
- }
+ image->data = malloc(image->width * image->height * sizeof(float));
+ for ( x=0; x<image->width; x++ ) {
+ for ( y=0; y<image->height; y++ ) {
+ int val;
+ val = image->hdr[x + image->width*y];
+ image->data[x + image->width*y] = val;
+ }
}
}