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..)

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);