indexamajig: Add --integrate-found

6 files changed, 86 insertions, 4 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 3519dc5f..9933627c 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -334,6 +334,11 @@ Normally, reflections which contain one or more pixels above max_adu (defined in
 .PD
 When using \fB--peaks=hdf5\fR, the peaks will be put through the same checks as if you were using \fB--peaks=zaef\fR.  These checks reject peaks which are too close to panel edges, are saturated (unless you use \fB--use-saturated\fR), fall short of the minimum SNR value given by \fB--min-snr\fR, have other nearby peaks (closer than twice the inner integration radius, see \fB--int-radius\fR), or have any part in a bad region.  Using this option skips this validation step, and uses the peaks directly.
 
+.PD 0
+.IP \fB--integrate-found\fR
+.PD
+Without this option, peaks will be predicted in each pattern based on the crystal orientation from autoindexing and the parameters in the beam description file.  With this option, indices will be assigned to the peaks found by the peak search, and the positions of those peaks used for the final integration stage.
+
 .SH BUGS
 Don't run more than one indexamajig jobs simultaneously in the same working
 directory - they'll overwrite each other's DirAx or MOSFLM files, causing subtle
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 60fab488..79a438ec 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -33,6 +33,7 @@
 
 #include <stdlib.h>
 #include <assert.h>
+#include <fenv.h>
 
 #include "utils.h"
 #include "cell.h"
@@ -276,6 +277,65 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 }
 
 
+RefList *select_intersections(struct image *image, UnitCell *cell)
+{
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	const double min_dist = 0.25;
+	RefList *list;
+	int i;
+
+	/* Round towards nearest */
+	fesetround(1);
+
+	/* Cell basis vectors for this image */
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	list = reflist_new();
+	if ( list == NULL ) return NULL;
+
+	/* Loop over peaks, checking proximity to nearest reflection */
+	for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+		struct imagefeature *f;
+		struct rvec q;
+		double h, k, l, hd, kd, ld;
+		double dsq;
+
+		f = image_get_feature(image->features, i);
+		if ( f == NULL ) continue;
+
+		/* Reciprocal space position of found peak */
+		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
+
+		/* Decimal and fractional Miller indices of nearest
+		 * reciprocal lattice point */
+		hd = q.u * ax + q.v * ay + q.w * az;
+		kd = q.u * bx + q.v * by + q.w * bz;
+		ld = q.u * cx + q.v * cy + q.w * cz;
+		h = lrint(hd);
+		k = lrint(kd);
+		l = lrint(ld);
+
+		/* Check distance */
+		dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
+
+		if ( sqrt(dsq) < min_dist ) {
+
+			Reflection *refl;
+
+			refl = add_refl(list, h, k, l);
+			set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
+
+		}
+
+	}
+
+	return list;
+}
+
+
 /* Calculate partialities and apply them to the image's reflections */
 void update_partialities(struct image *image)
 {
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 94b496e6..14a3a7d6 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -43,9 +43,10 @@
 extern "C" {
 #endif
 
-RefList *find_intersections(struct image *image, UnitCell *cell);
+extern RefList *find_intersections(struct image *image, UnitCell *cell);
+extern RefList *select_intersections(struct image *image, UnitCell *cell);
 
-void update_partialities(struct image *image);
+extern void update_partialities(struct image *image);
 
 #ifdef __cplusplus
 }
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 32bcba15..c8cceb2b 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -301,9 +301,17 @@ static void process_image(const struct index_args *iargs,
 	              iargs->config_insane, iargs->tols);
 
 	if ( image.indexed_cell != NULL ) {
+
 		pargs->indexable = 1;
-		image.reflections = find_intersections(&image,
-				image.indexed_cell);
+
+		if ( iargs->integrate_found ) {
+			image.reflections = select_intersections(&image,
+			                                    image.indexed_cell);
+		} else {
+			image.reflections = find_intersections(&image,
+			                                    image.indexed_cell);
+		}
+
 		if (image.reflections != NULL) {
 			integrate_reflections(&image,
 			                      iargs->config_closer,
@@ -314,6 +322,7 @@ static void process_image(const struct index_args *iargs,
 			                      iargs->ir_out,
 			                      iargs->integrate_saturated);
 		}
+
 	} else {
 		image.reflections = NULL;
 	}
diff --git a/src/im-sandbox.h b/src/im-sandbox.h
index 3851fb40..872f84ad 100644
--- a/src/im-sandbox.h
+++ b/src/im-sandbox.h
@@ -69,6 +69,7 @@ struct index_args
 	int integrate_saturated;
 	int use_saturated;
 	int no_revalidate;
+	int integrate_found;
 };
 
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index ae03986c..b132331c 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -167,6 +167,9 @@ static void show_help(const char *s)
 "                           peaks which contain pixels above max_adu.\n"
 "     --no-revalidate      Don't re-integrate and check HDF5 peaks for\n"
 "                           validity.\n"
+"     --integrate-found    Skip the spot prediction step, and just integrate\n"
+"                           the intensities of the spots found by the initial\n"
+"                           peak search.\n"
 );
 }
 
@@ -251,6 +254,7 @@ int main(int argc, char *argv[])
 	int integrate_saturated = 0;
 	int use_saturated = 0;
 	int no_revalidate = 0;
+	int integrate_found = 0;
 
 	copyme = new_copy_hdf5_field_list();
 	if ( copyme == NULL ) {
@@ -300,6 +304,7 @@ int main(int argc, char *argv[])
 		{"integrate-saturated",0, &integrate_saturated,1},
 		{"use-saturated",0, &use_saturated,            1},
 		{"no-revalidate",      0, &no_revalidate,      1},
+		{"integrate-found",    0, &integrate_found,    1},
 		{0, 0, NULL, 0}
 	};
 
@@ -649,6 +654,7 @@ int main(int argc, char *argv[])
 	iargs.use_saturated = use_saturated;
 	iargs.integrate_saturated = integrate_saturated;
 	iargs.no_revalidate = no_revalidate;
+	iargs.integrate_found = integrate_found;
 
 	create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
 	                    use_this_one_instead, ofh);