indexamajig: Record number of saturated reflections in stream

5 files changed, 30 insertions, 10 deletions

diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index ac011d14..afa9e4a7 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -104,6 +104,7 @@ typedef struct _imagefeaturelist ImageFeatureList;
  *    int                     height;
  *
  *    RefList                 *reflections;
+ *    long long unsigned int  n_saturated;
  *
  *    ImageFeatureList        *features;
  * };
@@ -169,6 +170,7 @@ struct image {
 
 	/* Integrated (or about-to-be-integrated) reflections */
 	RefList                 *reflections;
+	long long int           n_saturated;  /* Number of overloads */
 
 	/* Detected peaks */
 	ImageFeatureList        *features;
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 4ae4bb30..f5fa6c70 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -216,7 +216,7 @@ static int integrate_peak(struct image *image, int cfs, int css,
                           double *pfs, double *pss,
                           double *intensity, double *sigma,
                           double ir_inn, double ir_mid, double ir_out,
-                          int use_max_adu, int *bgPkMask)
+                          int use_max_adu, int *bgPkMask, int *saturated)
 {
 	signed int dfs, dss;
 	double lim_sq, out_lim_sq, mid_lim_sq;
@@ -292,7 +292,10 @@ static int integrate_peak(struct image *image, int cfs, int css,
 		val = image->data[idx];
 
 		/* Veto peak if it contains saturation in bg region */
-		if ( use_max_adu && (val > p->max_adu) ) return 1;
+		if ( use_max_adu && (val > p->max_adu) ) {
+			if ( saturated != NULL ) *saturated = 1;
+			return 1;
+		}
 
 		bg_tot += val;
 		bg_tot_sq += pow(val, 2.0);
@@ -347,7 +350,10 @@ static int integrate_peak(struct image *image, int cfs, int css,
 		val = image->data[idx] - bg_mean;
 
 		/* Veto peak if it contains saturation */
-		if ( use_max_adu && (image->data[idx] > p->max_adu) ) return 1;
+		if ( use_max_adu && (val > p->max_adu) ) {
+			if ( saturated != NULL ) *saturated = 1;
+			return 1;
+		}
 
 		pk_counts++;
 		pk_total += val;
@@ -481,7 +487,7 @@ static void search_peaks_in_panel(struct image *image, float threshold,
 		/* Centroid peak and get better coordinates. */
 		r = integrate_peak(image, mask_fs, mask_ss,
 		                   &f_fs, &f_ss, &intensity, &sigma,
-		                   ir_inn, ir_mid, ir_out, 0, NULL);
+		                   ir_inn, ir_mid, ir_out, 0, NULL, NULL);
 
 		if ( r ) {
 			/* Bad region - don't detect peak */
@@ -730,6 +736,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 		signed int h, k, l;
 		struct panel *p;
 		int pnum, j, found;
+		int saturated;
 
 		refl = il[i].refl;
 
@@ -781,7 +788,9 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub,
 
 		r = integrate_peak(image, pfs, pss, &fs, &ss,
 		                   &intensity, &sigma, ir_inn, ir_mid, ir_out,
-		                   1, bgMasks[pnum]);
+		                   1, bgMasks[pnum], &saturated);
+
+		if ( saturated ) image->n_saturated++;
 
 		/* I/sigma(I) cutoff */
 		if ( intensity/sigma < min_snr ) r = 1;
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 6627e7ea..b8ff9238 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -261,6 +261,9 @@ void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f)
 		             i->diffracting_resolution/1e9,
 		             1e9 / i->diffracting_resolution);
 
+		fprintf(ofh, "num_saturated_reflections"
+		             " = %lli\n", i->n_saturated);
+
 	}
 
 	if ( i->det != NULL ) {
@@ -343,6 +346,7 @@ int read_chunk(FILE *fh, struct image *image)
 	image->features = NULL;
 	image->reflections = NULL;
 	image->indexed_cell = NULL;
+	image->n_saturated = 0;
 
 	do {
 
@@ -417,6 +421,10 @@ int read_chunk(FILE *fh, struct image *image)
 			have_ev = 1;
 		}
 
+		if ( strncmp(line, "num_saturated_reflections = ", 28) == 0 ) {
+			image->n_saturated = atoi(line+28);
+		}
+
 		if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) {
 			if ( read_peaks(fh, image) ) {
 				ERROR("Failed while reading peaks\n");
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 21b091a5..c313bda9 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -186,6 +186,7 @@ static void process_image(const struct index_args *iargs,
 	image.det = copy_geom(iargs->det);
 	image.copyme = iargs->copyme;
 	image.reflections = NULL;
+	image.n_saturated = 0;
 	image.id = cookie;
 	image.filename = pargs->filename;
 	image.beam = beam;
diff --git a/tests/integration_check.c b/tests/integration_check.c
index acadbdf8..7962e414 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -65,7 +65,7 @@ static void third_integration_check(struct image *image, int n_trials,
 
 		r = integrate_peak(image, 64, 64, &fsp, &ssp,
 		                   &intensity, &sigma, 10.0, 15.0, 17.0,
-		                   0, NULL);
+		                   0, NULL, NULL);
 
 		if ( r == 0 ) {
 			mean_intensity += intensity;
@@ -127,7 +127,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
 
 		r = integrate_peak(image, 64, 64, &fsp, &ssp,
 		                   &intensity, &sigma, 10.0, 15.0, 17.0,
-		                   0, NULL);
+		                   0, NULL, NULL);
 
 		if ( r == 0 ) {
 			mean_intensity += intensity;
@@ -209,7 +209,7 @@ int main(int argc, char *argv[])
 
 	/* First check: no intensity -> no peak, or very low intensity */
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0, 0, NULL);
+	                   10.0, 15.0, 17.0, 0, NULL, NULL);
 	STATUS("  First check: integrate_peak() returned %i", r);
 	if ( r == 0 ) {
 
@@ -235,7 +235,7 @@ int main(int argc, char *argv[])
 	}
 
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0, 0, NULL);
+	                   10.0, 15.0, 17.0, 0, NULL, NULL);
 	if ( r ) {
 		ERROR(" Second check: integrate_peak() returned %i (wrong).\n",
 		      r);
@@ -277,7 +277,7 @@ int main(int argc, char *argv[])
 	}
 
 	r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma,
-	                   10.0, 15.0, 17.0, 0, NULL);
+	                   10.0, 15.0, 17.0, 0, NULL, NULL);
 	if ( r ) {
 		ERROR("   Fifth check: integrate_peak() returned %i (wrong).\n",
 		      r);