Take account of bad regions of detector

3 files changed, 74 insertions, 13 deletions

diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 5573267b..d635923d 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -318,6 +318,11 @@ int hdf5_read(struct hdfile *f, struct image *image)
 {
 	herr_t r;
 	float *buf;
+	uint8_t *flags;
+	hid_t mask_dh;
+
+	image->height = f->nx;
+	image->width = f->ny;
 
 	buf = malloc(sizeof(float)*f->nx*f->ny);
 
@@ -327,10 +332,22 @@ int hdf5_read(struct hdfile *f, struct image *image)
 		ERROR("Couldn't read data\n");
 		return 1;
 	}
-
 	image->data = buf;
-	image->height = f->nx;
-	image->width = f->ny;
+
+	mask_dh = H5Dopen(f->fh, "/processing/hitfinder/masks", H5P_DEFAULT);
+	if ( mask_dh <= 0 ) {
+		ERROR("Couldn't open flags\n");
+	} else {
+		flags = malloc(sizeof(uint8_t)*f->nx*f->ny);
+		r = H5Dread(mask_dh, H5T_NATIVE_B8, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, flags);
+		if ( r < 0 ) {
+			ERROR("Couldn't read flags\n");
+			image->flags = NULL;
+		} else {
+			image->flags = flags;
+		}
+	}
 
 	/* Read wavelength from file */
 	image->lambda = get_wavelength(f);
diff --git a/src/image.h b/src/image.h
index db1b491e..e3a520ca 100644
--- a/src/image.h
+++ b/src/image.h
@@ -81,6 +81,7 @@ struct reflhit {
 struct image {
 
 	float			*data;
+	uint8_t			*flags;
 	double			*twotheta;
 	UnitCell		*indexed_cell;
 	UnitCell		*candidate_cells[MAX_CELL_CANDIDATES];
diff --git a/src/peaks.c b/src/peaks.c
index 22e12949..0d4292d3 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -133,7 +133,9 @@ static void cull_peaks(struct image *image)
 }
 
 
-static void integrate_peak(struct image *image, int xp, int yp,
+/* Returns non-zero if peak has been vetoed.
+ * i.e. don't use result if return value is not zero. */
+static int integrate_peak(struct image *image, int xp, int yp,
                            float *xc, float *yc, float *intensity,
                            int do_polar)
 {
@@ -150,6 +152,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		double val, sa, pix_area, Lsq, dsq, proj_area;
 		float tt;
 		double phi, pa, pb, pol;
+		uint8_t flags;
 
 		/* Circular mask */
 		if ( x*x + y*y > lim ) continue;
@@ -157,6 +160,12 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue;
 		if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue;
 
+		/* Veto this peak if we tried to integrate in a bad region */
+		if ( image->flags != NULL ) {
+			flags = image->flags[(x+xp)+image->width*(y+yp)];
+			if ( flags & 0x03 ) return 1;
+		}
+
 		p = find_panel(&image->det, x+xp, y+yp);
 
 		/* Area of one pixel */
@@ -203,6 +212,8 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		*yc = (float)yp;
 		*intensity = 0;
 	}
+
+	return 0;
 }
 
 
@@ -219,6 +230,7 @@ void search_peaks(struct image *image)
 	int nrej_hot = 0;
 	int nrej_pro = 0;
 	int nrej_fra = 0;
+	int nrej_bad = 0;
 	int nacc = 0;
 
 	data = image->data;
@@ -240,6 +252,7 @@ void search_peaks(struct image *image)
 		int sx, sy;
 		double max;
 		unsigned int did_something;
+		int r;
 
 		/* Overall threshold */
 		if ( data[x+width*y] < 800 ) continue;
@@ -313,7 +326,13 @@ void search_peaks(struct image *image)
 		/* Centroid peak and get better coordinates.
 		 * Don't bother doing polarisation correction, because the
 		 * intensity of this peak is only an estimate at this stage. */
-		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0);
+		r = integrate_peak(image, mask_x, mask_y,
+		                   &fx, &fy, &intensity, 0);
+		if ( r ) {
+			/* Bad region - don't detect peak */
+			nrej_bad++;
+			continue;
+		}
 
 		/* It is possible for the centroid to fall outside the image */
 		if ( (fx < 0.0) || (fx > image->width)
@@ -337,8 +356,9 @@ void search_peaks(struct image *image)
 	}
 	}
 
-	STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame\n",
-	       nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra);
+	STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame,"
+	       "%i in bad regions.\n",
+	       nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra, nrej_bad);
 
 	cull_peaks(image);
 }
@@ -492,6 +512,8 @@ void output_intensities(struct image *image, UnitCell *cell,
 	int n_found;
 	int n_indclose = 0;
 	int n_foundclose = 0;
+	int n_veto = 0;
+	int n_veto_second = 0;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
@@ -543,17 +565,35 @@ void output_intensities(struct image *image, UnitCell *cell,
 		                          &d, &idx);
 		if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
 
+			int r;
+
 			/* f->intensity was measured on the filtered pattern,
 			 * so instead re-integrate using old coordinates.
 			 * This will produce further revised coordinates. */
-			integrate_peak(image, f->x, f->y, &x, &y, &intensity,
-			               !unpolar);
+			r = integrate_peak(image, f->x, f->y, &x, &y,
+			                   &intensity, !unpolar);
+			if ( r ) {
+				/* The original peak (which also went through
+				 * integrate_peak(), but with the mangled
+				 * image data) would have been rejected if it
+				 * was in a bad region.  Integration of the same
+				 * peak included a bad region this time. */
+				n_veto_second++;
+				continue;
+			}
 			intensity = f->intensity;
 
 		} else {
 
-			integrate_peak(image, hits[i].x, hits[i].y,
-			               &x, &y, &intensity, !unpolar);
+			int r;
+
+			r = integrate_peak(image, hits[i].x, hits[i].y,
+			                   &x, &y, &intensity, !unpolar);
+			if ( r ) {
+				/* Plain old ordinary peak veto */
+				n_veto++;
+				continue;
+			}
 
 		}
 
@@ -592,9 +632,12 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 	printf("Peak statistics: %i peaks found by the peak search out of "
 	       "%i were close to indexed positions. "
-	       "%i indexed positions out of %i were close to detected peaks\n",
+	       "%i indexed positions out of %i were close to detected peaks.\n",
 	       n_foundclose, n_found, n_indclose, n_hits);
-
+	printf("%i integrations using indexed locations were aborted because "
+	       "they hit one or more bad pixels.\n", n_veto);
+	printf("%i integrations using peak search locations were aborted "
+	       "because they hit one or more bad pixels.\n", n_veto_second);
 	/* Blank line at end */
 	printf("\n");