aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2009-11-27 16:10:57 +0100
committerThomas White <taw@physics.org>2009-11-27 16:10:57 +0100
commitef92cb3eebfb74c865cf0e10266ba8c46ffc8a9a (patch)
tree6e601f0a02cdf240ec065f71ca58a0ad31d3e76d
parent93cc3f1334294f8737dc7167f361af1209f3eaa8 (diff)
Poisson function returns integer count - do all downstream calculations as integers
-rw-r--r--src/detector.c51
-rw-r--r--src/image.h4
-rw-r--r--src/utils.c2
-rw-r--r--src/utils.h2
4 files changed, 28 insertions, 31 deletions
diff --git a/src/detector.c b/src/detector.c
index 52f24b1e..fd84a9c5 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -37,56 +37,56 @@
/* Bleed excess intensity into neighbouring pixels */
-static void bloom_values(double *tmp, int x, int y,
- int width, int height, double val)
+static void bloom_values(int *tmp, int x, int y,
+ int width, int height, int val)
{
- double overflow;
+ 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.0;
+ tmp[x-1 + width*y] += overflow / 6;
if ( y > 0 ) {
- tmp[x-1 + width*(y-1)] += overflow / 12.0;
+ tmp[x-1 + width*(y-1)] += overflow / 12;
}
if ( y < height-1 ) {
- tmp[x-1 + width*(y+1)] += overflow / 12.0;
+ tmp[x-1 + width*(y+1)] += overflow / 12;
}
}
if ( x < width-1 ) {
- tmp[x+1 + width*y] += overflow / 6.0;
+ tmp[x+1 + width*y] += overflow / 6;
if ( y > 0 ) {
- tmp[x+1 + width*(y-1)] += overflow / 12.0;
+ tmp[x+1 + width*(y-1)] += overflow / 12;
}
if ( y < height-1 ) {
- tmp[x+1 + width*(y+1)] += overflow / 12.0;
+ tmp[x+1 + width*(y+1)] += overflow / 12;
}
}
if ( y > 0 ) {
- tmp[x + width*(y-1)] += overflow / 6.0;
+ tmp[x + width*(y-1)] += overflow / 6;
}
if ( y < height-1 ) {
- tmp[x + width*(y+1)] += overflow / 6.0;
+ tmp[x + width*(y+1)] += overflow / 6;
}
}
-static uint16_t *bloom(double *hdr_in, int width, int height)
+static uint16_t *bloom(int *hdr_in, int width, int height)
{
int x, y;
uint16_t *data;
- double *tmp;
- double *hdr;
+ int *tmp;
+ int *hdr;
int did_something;
data = malloc(width * height * sizeof(uint16_t));
- tmp = malloc(width * height * sizeof(double));
- hdr = malloc(width * height * sizeof(double));
+ tmp = malloc(width * height * sizeof(int));
+ hdr = malloc(width * height * sizeof(int));
- memcpy(hdr, hdr_in, width*height*sizeof(double));
+ memcpy(hdr, hdr_in, width*height*sizeof(int));
/* Apply DQE (once only) */
for ( x=0; x<width; x++ ) {
@@ -97,13 +97,13 @@ static uint16_t *bloom(double *hdr_in, int width, int height)
do {
- memset(tmp, 0, width*height*sizeof(double));
+ memset(tmp, 0, width*height*sizeof(int));
did_something = 0;
for ( x=0; x<width; x++ ) {
for ( y=0; y<height; y++ ) {
- double hdval;
+ int hdval;
hdval = hdr[x + width*y];
@@ -122,7 +122,7 @@ static uint16_t *bloom(double *hdr_in, int width, int height)
/* Prepare new input if we're going round for another pass */
if ( did_something ) {
- memcpy(hdr, tmp, width*height*sizeof(double));
+ memcpy(hdr, tmp, width*height*sizeof(int));
}
} while ( did_something );
@@ -169,7 +169,8 @@ void record_image(struct image *image)
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
- double counts, intensity, sa, water;
+ int counts;
+ double intensity, sa, water;
double complex val;
double dsq, proj_area;
@@ -181,7 +182,6 @@ void record_image(struct image *image)
image->xray_energy,
BEAM_RADIUS, WATER_RADIUS);
- //printf("%e, %e, ", intensity, water);
intensity += water;
/* Area of pixel as seen from crystal (approximate) */
@@ -194,10 +194,7 @@ void record_image(struct image *image)
/* Projected area of pixel divided by distance squared */
sa = proj_area / (dsq + Lsq);
- counts = intensity * ph_per_e * sa;
- //printf("%e counts\n", counts);
-
- counts = poisson_noise(counts);
+ counts = poisson_noise(intensity * ph_per_e * sa * DQE);
image->hdr[x + image->width*y] = counts;
@@ -212,7 +209,7 @@ void record_image(struct image *image)
* sizeof(uint16_t));
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
- double val;
+ int val;
val = image->hdr[x + image->width*y];
if ( val > SATURATION ) val = SATURATION;
image->data[x + image->width*y] = (uint16_t)val;
diff --git a/src/image.h b/src/image.h
index e9db7cad..a7c6c355 100644
--- a/src/image.h
+++ b/src/image.h
@@ -62,8 +62,8 @@ struct threevec
/* Structure describing an image */
struct image {
- double *hdr; /* Actual counts */
- uint16_t *data; /* Integer counts after DQE/bloom */
+ int *hdr; /* Actual counts */
+ uint16_t *data; /* Integer counts after bloom */
double complex *sfacs;
struct threevec *qvecs;
double *twotheta;
diff --git a/src/utils.c b/src/utils.c
index 9aa235a2..36ae533e 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -74,7 +74,7 @@ void progress_bar(int val, int total, const char *text)
}
-double poisson_noise(double expected)
+int poisson_noise(double expected)
{
double L;
int k = 0;
diff --git a/src/utils.h b/src/utils.h
index bb0dd77e..31ecdf7b 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -69,7 +69,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, const char *text);
-extern double poisson_noise(double expected);
+extern int poisson_noise(double expected);
/* Keep these ones inline, to avoid function call overhead */
static inline struct quaternion invalid_quaternion(void)