indexamajig: Add --int-diag

7 files changed, 120 insertions, 25 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 57a2636e..351a2e45 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -327,6 +327,11 @@ Do not record peak search results in the stream.  You won't be able to check tha
 .PD
 Do not record integrated reflections in the stream.  The resulting output won't be usable for merging, but will be a lot smaller.  This option might be useful if you're only interested in things like unit cell parameters and orientations.
 
+.PD 0
+.IP \fB--int-diag=\fIcondition\fR\fR
+.PD
+Show detailed information about reflection integration when \fIcondition\fR is met.  The \fIcondition\fR can be \fBall\fR, \fBnone\fR, a set of Miller indices separated by commas, or \fBrandom\fR.  The default is \fB--int-diag=none\fR, and \fBrandom\fR means to show information about a random 1% of the peaks.
+
 .SH BUGS
 ReAx indexing is experimental.  It works very nicely for some people, and crashes for others.  In a future version, it will be improved and fully supported.
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index dec8b39f..c3d2730e 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -53,10 +53,6 @@
 #include "integration.h"
 
 
-#define VERBOSITY (0)
-// ((h==-6) && (k==0) && (l==-8))
-
-
 static void check_eigen(gsl_vector *e_val)
 {
 	int i;
@@ -182,6 +178,11 @@ struct intcontext
 	int n_saturated;
 	int n_implausible;
 	double limit;
+
+	IntDiag int_diag;
+	signed int int_diag_h;
+	signed int int_diag_k;
+	signed int int_diag_l;
 };
 
 
@@ -955,8 +956,6 @@ static int center_and_check_box(struct intcontext *ic, struct peak_box *bx,
 	if ( check_box(ic, bx, sat) ) return 1;
 	fit_bg(ic, bx);
 
-	if ( bx->verbose ) show_peak_box(ic, bx);
-
 	for ( i=0; i<10; i++ ) {
 
 		int p, q;
@@ -1013,7 +1012,6 @@ static int center_and_check_box(struct intcontext *ic, struct peak_box *bx,
 		}
 
 		fit_bg(ic, bx);
-		if ( bx->verbose ) show_peak_box(ic, bx);
 
 		if ( (ifs==0) && (iss==0) ) break;
 
@@ -1320,8 +1318,32 @@ static void refine_rigid_groups(struct intcontext *ic)
 }
 
 
+static int get_int_diag(struct intcontext *ic, Reflection *refl)
+{
+	if ( ic->int_diag == INTDIAG_NONE ) return 0;
+
+	if ( ic->int_diag == INTDIAG_ALL ) return 1;
+
+	if ( ic->int_diag == INTDIAG_RANDOM ) {
+		return random() < RAND_MAX/100;
+	}
+
+	if ( ic->int_diag == INTDIAG_INDICES ) {
+		signed int h, k, l;
+		get_indices(refl, &h, &k, &l);
+		if ( ic->int_diag_h != h ) return 0;
+		if ( ic->int_diag_k != k ) return 0;
+		if ( ic->int_diag_l != l ) return 0;
+		return 1;
+	}
+
+	return 0;
+}
+
+
 static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
-                             struct image *image,
+                             struct image *image, IntDiag int_diag,
+                             signed int idh, signed int idk, signed int idl,
                              double ir_inn, double ir_mid, double ir_out)
 {
 	RefList *list;
@@ -1344,6 +1366,10 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 	ic.n_saturated = 0;
 	ic.n_implausible = 0;
 	ic.cell = cell;
+	ic.int_diag = int_diag;
+	ic.int_diag_h = idh;
+	ic.int_diag_k = idk;
+	ic.int_diag_l = idl;
 	if ( init_intcontext(&ic) ) {
 		ERROR("Failed to initialise integration.\n");
 		return;
@@ -1392,9 +1418,7 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 		bx->pn = pn;
 
 		get_indices(refl, &h, &k, &l);
-		if ( VERBOSITY ) {
-			bx->verbose = 1;
-		}
+		bx->verbose = get_int_diag(&ic, refl);
 
 		/* Which reference profile? */
 		bx->rp = 0;//bx->pn;
@@ -1449,7 +1473,6 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 			get_indices(bx->refl, &h, &k, &l);
 			STATUS("NaN intensity for %i %i %i !\n", h, k, l);
 			STATUS("panel %s\n", image->det->panels[bx->pn].name);
-			show_peak_box(&ic, bx);
 			show_reference_profile(&ic, bx->rp);
 		}
 		if ( bx->intensity < 0.0 ) {
@@ -1458,7 +1481,6 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 			STATUS("Negative intensity (%f) for %i %i %i !\n",
 			       bx->intensity, h, k, l);
 			STATUS("panel %s\n", image->det->panels[bx->pn].name);
-			show_peak_box(&ic, bx);
 			show_reference_profile(&ic, bx->rp);
 		}
 #endif
@@ -1499,7 +1521,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	int pn;
 	struct panel *p;
 	int fid_fs, fid_ss;  /* Center coordinates, rounded,
-				* in overall data block */
+	                      * in overall data block */
 	int cfs, css;  /* Corner coordinates */
 	double intensity;
 	double sigma;
@@ -1514,11 +1536,11 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	get_detector_pos(refl, &pfs, &pss);
 
 	/* Explicit truncation of digits after the decimal point.
-		* This is actually the correct thing to do here, not
-		* e.g. lrint().  pfs/pss is the position of the spot, measured
-		* in numbers of pixels, from the panel corner (not the center
-		* of the first pixel).  So any coordinate from 2.0 to 2.9999
-		* belongs to pixel index 2. */
+	 * This is actually the correct thing to do here, not
+	 * e.g. lrint().  pfs/pss is the position of the spot, measured
+	 * in numbers of pixels, from the panel corner (not the center
+	 * of the first pixel).  So any coordinate from 2.0 to 2.9999
+	 * belongs to pixel index 2. */
 	fid_fs = pfs;
 	fid_ss = pss;
 	pn = find_panel_number(image->det, fid_fs, fid_ss);
@@ -1534,7 +1556,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	bx->p = p;
 	bx->pn = pn;
 	get_indices(refl, &h, &k, &l);
-	bx->verbose = VERBOSITY;
+	bx->verbose = get_int_diag(ic, refl);
 
 	if ( ic->meth & INTEGRATION_CENTER ) {
 		r = center_and_check_box(ic, bx, &saturated);
@@ -1542,7 +1564,6 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 		r = check_box(ic, bx, &saturated);
 		if ( !r ) {
 			fit_bg(ic, bx);
-			if ( bx->verbose ) show_peak_box(ic, bx);
 		}
 		bx->offs_fs = 0.0;
 		bx->offs_ss = 0.0;
@@ -1614,11 +1635,14 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	pfs += bx->offs_fs;
 	pss += bx->offs_ss;
 	set_detector_pos(refl, 0.0, pfs, pss);
+
+	if ( bx->verbose ) show_peak_box(ic, bx);
 }
 
 
 static void integrate_rings(IntegrationMethod meth, Crystal *cr,
-                            struct image *image,
+                            struct image *image, IntDiag int_diag,
+                            signed int idh, signed int idk, signed int idl,
                             double ir_inn, double ir_mid, double ir_out)
 {
 	RefList *list;
@@ -1643,6 +1667,10 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
 	ic.ir_inn = ir_inn;
 	ic.ir_mid = ir_mid;
 	ic.ir_out = ir_out;
+	ic.int_diag = int_diag;
+	ic.int_diag_h = idh;
+	ic.int_diag_k = idk;
+	ic.int_diag_l = idl;
 	ic.limit = 0.0;
 	if ( init_intcontext(&ic) ) {
 		ERROR("Failed to initialise integration.\n");
@@ -1783,7 +1811,9 @@ static void resolution_cutoff(RefList *list, UnitCell *cell)
 
 
 void integrate_all(struct image *image, IntegrationMethod meth,
-                   double ir_inn, double ir_mid, double ir_out)
+                   double ir_inn, double ir_mid, double ir_out,
+                   IntDiag int_diag,
+                   signed int idh, signed int idk, signed int idl)
 {
 	int i;
 
@@ -1796,11 +1826,13 @@ void integrate_all(struct image *image, IntegrationMethod meth,
 
 			case INTEGRATION_RINGS :
 			integrate_rings(meth, image->crystals[i], image,
+			                int_diag, idh, idk, idl,
 			                ir_inn, ir_mid, ir_out);
 			break;
 
 			case INTEGRATION_PROF2D :
 			integrate_prof2d(meth, image->crystals[i], image,
+			                 int_diag, idh, idk, idl,
 			                 ir_inn, ir_mid, ir_out);
 			break;
 
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
index 1273348d..68ce91ec 100644
--- a/libcrystfel/src/integration.h
+++ b/libcrystfel/src/integration.h
@@ -33,6 +33,16 @@
 #include <config.h>
 #endif
 
+
+typedef enum {
+
+	INTDIAG_NONE,
+	INTDIAG_RANDOM,
+	INTDIAG_ALL,
+	INTDIAG_INDICES
+
+} IntDiag;
+
 #define INTEGRATION_DEFAULTS_RINGS (INTEGRATION_RINGS)
 #define INTEGRATION_DEFAULTS_PROF2D (INTEGRATION_PROF2D | INTEGRATION_CENTER)
 
@@ -68,6 +78,9 @@ typedef enum {
 extern IntegrationMethod integration_method(const char *t, int *err);
 
 extern void integrate_all(struct image *image, IntegrationMethod meth,
-	                  double ir_inn, double ir_mid, double ir_out);
+	                  double ir_inn, double ir_mid, double ir_out,
+                          IntDiag int_diag,
+                          signed int idh, signed int idk, signed int idl);
+
 
 #endif	/* INTEGRATION_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 7d3747ba..823ebb79 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -144,6 +144,7 @@ static void show_help(const char *s)
 "                           validity.\n"
 "     --no-peaks-in-stream Do not record peak search results in the stream.\n"
 "     --no-refls-in-stream Do not record integrated reflections in the stream.\n"
+"     --int-diag=<r>       Show debugging information about this reflection.\n"
 );
 }
 
@@ -173,6 +174,7 @@ int main(int argc, char *argv[])
 	char *intrad = NULL;
 	char *int_str = NULL;
 	char *tempdir = NULL;
+	char *int_diag = NULL;
 
 	/* Defaults */
 	iargs.cell = NULL;
@@ -201,6 +203,7 @@ int main(int argc, char *argv[])
 	iargs.no_revalidate = 0;
 	iargs.stream_peaks = 1;
 	iargs.stream_refls = 1;
+	iargs.int_diag = INTDIAG_NONE;
 	iargs.copyme = new_copy_hdf5_field_list();
 	if ( iargs.copyme == NULL ) {
 		ERROR("Couldn't allocate HDF5 field list.\n");
@@ -256,6 +259,7 @@ int main(int argc, char *argv[])
 		{"median-filter",      1, NULL,               15},
 		{"integration",        1, NULL,               16},
 		{"temp-dir",           1, NULL,               17},
+		{"int-diag",           1, NULL,               18},
 
 		{0, 0, NULL, 0}
 	};
@@ -384,6 +388,10 @@ int main(int argc, char *argv[])
 			tempdir = strdup(optarg);
 			break;
 
+			case 18 :
+			int_diag = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -527,6 +535,36 @@ int main(int argc, char *argv[])
 		iargs.cell = NULL;
 	}
 
+	if ( int_diag != NULL ) {
+
+		int r;
+		signed int h, k, l;
+
+		if ( strcmp(int_diag, "random") == 0 ) {
+			iargs.int_diag = INTDIAG_RANDOM;
+		}
+
+		if ( strcmp(int_diag, "all") == 0 ) {
+			iargs.int_diag = INTDIAG_ALL;
+		}
+
+		r = sscanf(int_diag, "%i,%i,%i", &h, &k, &l);
+		if ( r == 3 ) {
+			iargs.int_diag = INTDIAG_INDICES;
+			iargs.int_diag_h = h;
+			iargs.int_diag_k = k;
+			iargs.int_diag_l = l;
+		}
+
+		if ( iargs.int_diag == INTDIAG_NONE ) {
+			ERROR("Invalid value for --int-diag.\n");
+			return 1;
+		}
+
+		free(int_diag);
+
+	}
+
 	ofd = open(outfile, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR);
 	if ( ofd == -1 ) {
 		ERROR("Failed to open stream '%s'\n", outfile);
diff --git a/src/process_image.c b/src/process_image.c
index b07e16b6..d1e3a7de 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -184,7 +184,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	/* Integrate all the crystals at once - need all the crystals so that
 	 * overlaps can be detected. */
 	integrate_all(&image, iargs->int_meth,
-	                      iargs->ir_inn, iargs->ir_mid, iargs->ir_out);
+	                      iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+	                      iargs->int_diag, iargs->int_diag_h,
+	                      iargs->int_diag_k, iargs->int_diag_l);
 
 	write_chunk(st, &image, hdfile,
 	            iargs->stream_peaks, iargs->stream_refls);
diff --git a/src/process_image.h b/src/process_image.h
index b9339cd1..90113ff9 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -73,6 +73,10 @@ struct index_args
 	int stream_peaks;
 	int stream_refls;
 	IntegrationMethod int_meth;
+	IntDiag int_diag;
+	signed int int_diag_h;
+	signed int int_diag_k;
+	signed int int_diag_l;
 };
 
 
diff --git a/tests/integration_check.c b/tests/integration_check.c
index a63e27f6..c469edca 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -138,6 +138,7 @@ int main(int argc, char *argv[])
 		ic.ir_out = ir_out;
 		ic.limit = 0.0;
 		ic.meth = INTEGRATION_RINGS;
+		ic.int_diag = INTDIAG_NONE;
 		if ( init_intcontext(&ic) ) {
 			ERROR("Failed to initialise integration.\n");
 			return 1;