aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/detector.c135
-rw-r--r--src/detector.h3
-rw-r--r--src/hdf5-file.c4
-rw-r--r--src/image.h2
-rw-r--r--src/indexamajig.c4
-rw-r--r--src/intensities.c8
-rw-r--r--src/pattern_sim.c17
-rw-r--r--src/peaks.c2
-rw-r--r--src/render.c23
9 files changed, 35 insertions, 163 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;
+ }
}
}
diff --git a/src/detector.h b/src/detector.h
index 4361d580..f71768c6 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -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);
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 3fc5de8c..d8ff7038 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -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");
diff --git a/src/image.h b/src/image.h
index 1c3d2a89..a8103b92 100644
--- a/src/image.h
+++ b/src/image.h
@@ -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;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index c7f0279c..5204b15e 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -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);
}
diff --git a/src/intensities.c b/src/intensities.c
index 4ec496c3..b0eef542 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -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);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index dc6ac354..a6a2156b 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -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);
}
diff --git a/src/peaks.c b/src/peaks.c
index 8e757c40..506263b2 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -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;
diff --git a/src/render.c b/src/render.c
index 9bdb9a75..b02469f6 100644
--- a/src/render.c
+++ b/src/render.c
@@ -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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 15:54:31 +0000