Commit 86dd71e8 authored by Thomas White's avatar Thomas White
Browse files

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..)
parent b4664bc2
......@@ -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;
}
}
}
......
......@@ -38,8 +38,7 @@ struct detector
int n_panels;
};
extern void record_image(struct image *image, int do_water, int do_poisson,
int do_bloom);
extern void record_image(struct image *image, int do_water, int do_poisson);
extern struct panel *find_panel(struct detector *det, int x, int y);
......
......@@ -176,11 +176,11 @@ static double get_wavelength(struct hdfile *f)
int hdf5_read(struct hdfile *f, struct image *image)
{
herr_t r;
int16_t *buf;
float *buf;
buf = malloc(sizeof(float)*f->nx*f->ny);
r = H5Dread(f->dh, H5T_NATIVE_INT16, H5S_ALL, H5S_ALL,
r = H5Dread(f->dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
H5P_DEFAULT, buf);
if ( r < 0 ) {
ERROR("Couldn't read data\n");
......
......@@ -70,7 +70,7 @@ struct rvec
struct image {
int *hdr; /* Actual counts */
int16_t *data; /* Integer counts after bloom */
float *data; /* Integer counts after bloom */
double complex *sfacs;
double *twotheta;
struct molecule *molecule;
......
......@@ -261,11 +261,11 @@ int main(int argc, char *argv[])
ERROR("Couldn't open molecule.pdb\n");
return 1;
}
record_image(&image, 0, 0, 0);
record_image(&image, 0, 0);
hdf5_write("simulated.h5", image.data,
image.width, image.height,
H5T_NATIVE_INT16);
H5T_NATIVE_FLOAT);
}
......
......@@ -35,10 +35,10 @@ struct reflhit {
};
static int sum_nearby_points(int16_t *data, int width, int x, int y)
static double sum_nearby_points(float *data, int width, int x, int y)
{
int dx, dy;
int intensity = 0;
double intensity = 0;
for ( dx=-3; dx<=3; dx++ ) {
for ( dy=-3; dy<=3; dy++ ) {
......@@ -127,7 +127,7 @@ void output_intensities(struct image *image, UnitCell *cell)
image->orientation.y, image->orientation.z);
for ( i=0; i<n_hits; i++ ) {
int intensity;
double intensity;
/* Bounds check */
if ( hits[i].x + 3 >= image->width ) continue;
......@@ -138,7 +138,7 @@ void output_intensities(struct image *image, UnitCell *cell)
intensity = sum_nearby_points(image->data, image->width,
hits[i].x, hits[i].y);
printf("%3i %3i %3i %6i (at %i,%i)\n",
printf("%3i %3i %3i %6f (at %i,%i)\n",
hits[i].h, hits[i].k, hits[i].l, intensity,
hits[i].x, hits[i].y);
......
......@@ -57,9 +57,6 @@ static void show_help(const char *s)
"\n"
" --no-water Do not simulate water background.\n"
" --no-noise Do not calculate Poisson noise.\n"
" --no-bloom Do not calculate CCD bloom (intensities which are\n"
" above the recordable range will be clamped to\n"
" the maximum allowable value).\n"
" --no-sfac Pretend that all structure factors are 1.\n"
);
}
......@@ -109,12 +106,7 @@ static void show_details()
"algorithm. When the intensity is sufficiently high that Knuth's algorithm\n"
"would result in machine precision problems, a normal distribution with\n"
"standard deviation sqrt(I) is used instead.\n"
"\n"
"Bloom of the CCD is included. Any excess intensity in a particular pixel\n"
"is divided between the neighbouring pixels. Diagonal neighbours receive\n"
"half the contribution of adjacent pixels. This process is repeated for\n"
"every pixel until all pixels are below the saturation value. Note that this\n"
"process is slow for very saturated images.\n");
);
}
......@@ -164,7 +156,6 @@ int main(int argc, char *argv[])
int config_noimages = 0;
int config_nowater = 0;
int config_nonoise = 0;
int config_nobloom = 0;
int config_nosfac = 0;
int config_gpu = 0;
int config_powder = 0;
......@@ -184,7 +175,6 @@ int main(int argc, char *argv[])
{"no-images", 0, &config_noimages, 1},
{"no-water", 0, &config_nowater, 1},
{"no-noise", 0, &config_nonoise, 1},
{"no-bloom", 0, &config_nobloom, 1},
{"no-sfac", 0, &config_nosfac, 1},
{"powder", 0, &config_powder, 1},
{0, 0, NULL, 0}
......@@ -294,8 +284,7 @@ int main(int argc, char *argv[])
goto skip;
}
record_image(&image, !config_nowater, !config_nonoise,
!config_nobloom);
record_image(&image, !config_nowater, !config_nonoise);
if ( config_nearbragg ) {
output_intensities(&image, image.molecule->cell);
......@@ -329,7 +318,7 @@ int main(int argc, char *argv[])
/* Write the output file */
hdf5_write(filename, image.data,
image.width, image.height, H5T_NATIVE_INT16);
image.width, image.height, H5T_NATIVE_FLOAT);
}
......
......@@ -307,7 +307,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
void search_peaks(struct image *image)
{
int x, y, width, height;
int16_t *data;
float *data;
data = image->data;
width = image->width;
......
......@@ -24,25 +24,23 @@
#include "filters.h"
static void *render_bin(int16_t *in, int inw, int inh,
int binning, int16_t *maxp)
static void *render_bin(float *in, int inw, int inh, int binning, float *maxp)
{
int16_t *data;
float *data;
int x, y;
int w, h;
int16_t max;
float max;
w = inw / binning;
h = inh / binning; /* Some pixels might get discarded */
data = malloc(w*h*sizeof(int16_t));
max = 0;
data = malloc(w*h*sizeof(float));
max = 0.0;
for ( x=0; x<w; x++ ) {
for ( y=0; y<h; y++ ) {
/* Big enough to hold large values */
unsigned long long int total;
double total;
size_t xb, yb;
total = 0;
......@@ -65,10 +63,10 @@ static void *render_bin(int16_t *in, int inw, int inh,
}
int16_t *render_get_image_binned(DisplayWindow *dw, int binning, int16_t *max)
float *render_get_image_binned(DisplayWindow *dw, int binning, float *max)
{
struct image *image;
int16_t *data;
float *data;
if ( (dw->image == NULL) || (dw->image_dirty) ) {
......@@ -85,6 +83,7 @@ int16_t *render_get_image_binned(DisplayWindow *dw, int binning, int16_t *max)
if ( dw->image != NULL ) {
image->features = dw->image->features;
if ( dw->image->data != NULL ) free(dw->image->data);
free(dw->image);
}
dw->image = image;
......@@ -209,9 +208,9 @@ GdkPixbuf *render_get_image(DisplayWindow *dw)
{
int mw, mh, w, h;
guchar *data;
int16_t *hdr;
float *hdr;
size_t x, y;
int16_t max;
float max;
mw = hdfile_get_width(dw->hdfile);
mh = hdfile_get_height(dw->hdfile);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment