Honour "-sat"

2 files changed, 59 insertions, 57 deletions

diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 5decd69a..9c5c8825 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -673,7 +673,7 @@ static void calculate_reference_profiles(struct intcontext *ic)
 }
 
 
-static int check_box(struct intcontext *ic, struct peak_box *bx)
+static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat)
 {
 	int p, q;
 	int n_pk = 0;
@@ -719,6 +719,11 @@ static int check_box(struct intcontext *ic, struct peak_box *bx)
 			bx->bm[p+ic->w*q] = BM_IG;
 		}
 
+		if ( (bx->bm[p+ic->w*q] != BM_IG)
+		  && (boxi(ic, bx, p, q) > bx->p->max_adu) ) {
+			*sat = 1;
+		}
+
 		/* Ignore if this pixel is closer to the next reciprocal lattice
 		 * point */
 		dv = get_q_for_panel(bx->p, fs, ss, NULL, ic->k);
@@ -910,13 +915,14 @@ static int suitable_reference(struct intcontext *ic, struct peak_box *bx)
 }
 
 
-static void measure_all_intensities(RefList *list, struct image *image,
-                                    UnitCell *cell)
+static void measure_all_intensities(IntegrationMethod meth, RefList *list,
+                                    struct image *image, UnitCell *cell)
 {
 	Reflection *refl;
 	RefListIterator *iter;
 	struct intcontext ic;
 	int i;
+	int n_saturated = 0;
 
 	ic.halfw = 4;
 	ic.image = image;
@@ -940,6 +946,7 @@ static void measure_all_intensities(RefList *list, struct image *image,
 		int fid_fs, fid_ss;  /* Center coordinates, rounded,
 		                      * in overall data block */
 		int cfs, css;  /* Corner coordinates */
+		int saturated;
 
 		set_redundancy(refl, 0);
 
@@ -967,11 +974,19 @@ static void measure_all_intensities(RefList *list, struct image *image,
 		/* Which reference profile? */
 		bx->rp = 0;//bx->pn;
 
-		if ( check_box(&ic, bx) ) {
+		if ( check_box(&ic, bx, &saturated) ) {
 			delete_box(&ic, bx);
 			continue;
 		}
 
+		if ( saturated ) {
+			n_saturated++;
+			if ( !(meth & INTEGRATION_SATURATED) ) {
+				delete_box(&ic, bx);
+				continue;
+			}
+		}
+
 		get_indices(refl, &h, &k, &l);
 		if ( (h==-24) && (k==6) && (l==-12) ) {
 			bx->verbose = 1;
@@ -1034,10 +1049,13 @@ static void measure_all_intensities(RefList *list, struct image *image,
 			set_redundancy(bx->refl, 1);
 		}
 	}
+
+	image->num_saturated_peaks = n_saturated;
 }
 
 
-static void estimate_mosaicity(Crystal *cr, struct image *image)
+static void estimate_mosaicity(IntegrationMethod meth, Crystal *cr,
+                               struct image *image)
 {
 	int msteps = 50;
 	int i;
@@ -1050,7 +1068,7 @@ static void estimate_mosaicity(Crystal *cr, struct image *image)
 	crystal_set_mosaicity(cr, mmax);
 	list = find_intersections(image, cr);
 	crystal_set_reflections(cr, list);
-	measure_all_intensities(list, image, crystal_get_cell(cr));
+	measure_all_intensities(meth, list, image, crystal_get_cell(cr));
 
 	for ( i=1; i<=msteps; i++ ) {
 
@@ -1208,10 +1226,9 @@ static void estimate_resolution(RefList *reflections, Crystal *cr,
 }
 
 
-static void integrate_prof2d(Crystal *cr, struct image *image, int use_closer,
-                             double min_snr,
-                             double ir_inn, double ir_mid, double ir_out,
-                             int integrate_saturated)
+static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
+                             struct image *image,  double min_snr,
+                             double ir_inn, double ir_mid, double ir_out)
 {
 	RefList *reflections;
 	UnitCell *cell;
@@ -1220,7 +1237,7 @@ static void integrate_prof2d(Crystal *cr, struct image *image, int use_closer,
 
 	/* Create initial list of reflections with nominal parameters */
 	reflections = find_intersections(image, cr);
-	measure_all_intensities(reflections, image, cell);
+	measure_all_intensities(meth, reflections, image, cell);
 
 	/* Find resolution limit of pattern using this list */
 	//estimate_resolution(reflections, cr, image);
@@ -1244,7 +1261,7 @@ static void integrate_prof2d(Crystal *cr, struct image *image, int use_closer,
 
 
 static void integrate_box(struct intcontext *ic, struct peak_box *bx,
-                          double *intensity, double *sigma, int *saturated)
+                          double *intensity, double *sigma)
 {
 	double pk_total;
 	int pk_counts;
@@ -1257,8 +1274,6 @@ static void integrate_box(struct intcontext *ic, struct peak_box *bx,
 	double aduph;
 	int p, q;
 
-	if ( saturated != NULL ) *saturated = 0;
-
 	aduph = bx->p->adu_per_eV * ph_lambda_to_eV(ic->image->lambda);
 
 	/* Measure the background */
@@ -1271,10 +1286,6 @@ static void integrate_box(struct intcontext *ic, struct peak_box *bx,
 		double bi = boxi(ic, bx, p, q);
 		if ( bx->bm[p + ic->w*q] != BM_BG ) continue;
 
-		if ( (saturated != NULL) && (bi > bx->p->max_adu) ) {
-			*saturated = 1;
-		}
-
 		bg_tot += bi;
 		bg_tot_sq += pow(bi, 2.0);
 		bg_counts++;
@@ -1296,10 +1307,6 @@ static void integrate_box(struct intcontext *ic, struct peak_box *bx,
 		double bi = boxi(ic, bx, p, q);
 		if ( bx->bm[p + ic->w*q] != BM_PK ) continue;
 
-		if ( (saturated != NULL) && (bi > bx->p->max_adu) ) {
-			*saturated = 1;
-		}
-
 		pk_counts++;
 		pk_total += (bi - bg_mean);
 		fsct += (bi-bg_mean)*p;
@@ -1319,9 +1326,9 @@ static void integrate_box(struct intcontext *ic, struct peak_box *bx,
 }
 
 
-static void integrate_rings(Crystal *cr, struct image *image, double min_snr,
-                            double ir_inn, double ir_mid, double ir_out,
-                            int integrate_saturated)
+static void integrate_rings(IntegrationMethod meth, Crystal *cr,
+                            struct image *image, double min_snr,
+                            double ir_inn, double ir_mid, double ir_out)
 {
 	RefList *list;
 	Reflection *refl;
@@ -1367,7 +1374,6 @@ static void integrate_rings(Crystal *cr, struct image *image, double min_snr,
 		double sigma, snr;
 		int saturated;
 		double one_over_d;
-		int r;
 
 		set_redundancy(refl, 0);
 
@@ -1392,43 +1398,39 @@ static void integrate_rings(Crystal *cr, struct image *image, double min_snr,
 		bx->p = p;
 		bx->pn = pn;
 
-		if ( check_box(&ic, bx) ) {
+		if ( check_box(&ic, bx, &saturated) ) {
 			delete_box(&ic, bx);
 			continue;
 		}
 
+		if ( saturated ) {
+			n_saturated++;
+			if ( !(meth & INTEGRATION_SATURATED) ) {
+				delete_box(&ic, bx);
+				continue;
+			}
+		}
+
 		get_indices(refl, &h, &k, &l);
-//		if ( (h==-1) && (k==14) && (l==13) ) {
-//			bx->verbose = 1;
-//		}
+		if ( (h==0) && (k==13) && (l==13) ) {
+			bx->verbose = 1;
+		}
 
 		if ( bx->verbose ) show_peak_box(&ic, bx);
 
-		integrate_box(&ic, bx, &intensity, &sigma, &saturated);
-
-		r = 0;
-		if ( saturated ) {
-			n_saturated++;
-			if ( !integrate_saturated ) r = 1;
-		}
+		integrate_box(&ic, bx, &intensity, &sigma);
 
 		/* I/sigma(I) cutoff
 		 * Rejects reflections below --min-integration-snr, or if the
 		 * SNR is clearly silly.  Silly indicates that the intensity
 		 * was zero. */
 		snr = fabs(intensity)/sigma;
-		if ( !r && (isnan(snr) || (snr < min_snr)) ) {
-			r = 1;
-		}
+		if ( (isnan(snr) || (snr < min_snr)) ) continue;
 
 		/* Record intensity and set redundancy to 1 on success */
-		if ( !r ) {
-			set_intensity(refl, intensity);
-			set_esd_intensity(refl, sigma);
-			set_redundancy(refl, 1);
-		} else {
-			set_redundancy(refl, 0);
-		}
+		set_intensity(refl, intensity);
+		set_esd_intensity(refl, sigma);
+		set_redundancy(refl, 1);
 
 		one_over_d = resolution(cell, h, k, l);
 		if ( one_over_d > limit ) limit = one_over_d;
@@ -1456,15 +1458,13 @@ void integrate_all(struct image *image, IntegrationMethod meth,
 			return;
 
 			case INTEGRATION_RINGS :
-			integrate_rings(image->crystals[i], image, min_snr,
-			                ir_inn, ir_mid, ir_out,
-				        integrate_saturated);
+			integrate_rings(meth, image->crystals[i], image,
+			                min_snr, ir_inn, ir_mid, ir_out);
 			return;
 
 			case INTEGRATION_PROF2D :
-			integrate_prof2d(image->crystals[i], image, use_closer,
-			                 min_snr, ir_inn, ir_mid, ir_out,
-				         integrate_saturated);
+			integrate_prof2d(meth, image->crystals[i], image,
+			                 min_snr, ir_inn, ir_mid, ir_out);
 			return;
 
 			default :
@@ -1497,6 +1497,9 @@ IntegrationMethod integration_method(const char *str, int *err)
 		} else if ( strcmp(methods[i], "none") == 0) {
 			return INTEGRATION_NONE;
 
+		} else if ( strcmp(methods[i], "sat") == 0) {
+			meth |= INTEGRATION_SATURATED;
+
 		} else {
 			ERROR("Bad integration method: '%s'\n", str);
 			if ( err != NULL ) *err = 1;
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
index 36f1d5bd..fe939e21 100644
--- a/libcrystfel/src/integration.h
+++ b/libcrystfel/src/integration.h
@@ -33,8 +33,8 @@
 #include <config.h>
 #endif
 
-#define INTEGRATION_DEFAULTS_RINGS (INTEGRATION_RINGS | INTEGRATION_SATURATED)
-#define INTEGRATION_DEFAULTS_PROF2D (INTEGRATION_PROF2D | INTEGRATION_SATURATED)
+#define INTEGRATION_DEFAULTS_RINGS (INTEGRATION_RINGS)
+#define INTEGRATION_DEFAULTS_PROF2D (INTEGRATION_PROF2D)
 
 /**
  * IntegrationMethod:
@@ -56,8 +56,7 @@ typedef enum {
 
 	/* Bits at the top of the IntegrationMethod are flags which modify the
 	 * behaviour of the integration. */
-	INTEGRATION_SATURATED = 256,
-	INTEGRATION_CLOSER = 512,
+	INTEGRATION_SATURATED = 256
 
 } IntegrationMethod;