indexamajig: Add --camera-length-estimate

Unfortunately, PinkIndexer needs the real camera length for its centre
refinement.  Giving a fake value and scaling the resulting shift does
not work - the indexing rate drops with even a small error.

Ideally, this would work in the same way as --wavelength-estimate, by
using a static value from the geometry file if it's given.  However,
this is rather complicated to implement because of the way all the units
stuff is implemented.  Therefore, this is left as an improvement for the
future.

8 files changed, 42 insertions, 52 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 6dffdffd..f94bec01 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -386,8 +386,9 @@ Check that most of the peaks can be accounted for by the indexing solution.  Thi
 
 .PD 0
 .IP \fB--wavelength-estimate=\fIm\fR
+.IP \fB--camera-length-estimate=\fIm\fR
 .PD
-Some indexing algorithms need to know the wavelength of the incident radiation in advance, e.g. to prepare an internal look-up table.  However, if the wavelength is taken from image headers, then the wavelength is not available at start-up.  In this case, you will be prompted to add this option to give an approximate wavelength, in metres.  A warning will be generated if the actual wavelength differs from this value by more than 10%.
+Some indexing algorithms need to know the camera length or the wavelength of the incident radiation in advance, e.g. to prepare an internal look-up table.  However, if these values are taken from image headers, then they are not available at start-up.  In this case, you will be prompted to add one of these options to give approximate values (in metres).  A warning will be generated if the actual value differs from this value by more than 10%.
 
 .PD 0
 .IP \fB--max-indexer-threads=\fIn\fR
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 1c78cc58..37ced03d 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -204,6 +204,7 @@ static char *friendly_indexer_name(IndexingMethod m)
 
 static void *prepare_method(IndexingMethod *m, UnitCell *cell,
                             double wavelength_estimate,
+                            double clen_estimate,
                             struct xgandalf_options *xgandalf_opts,
                             struct pinkIndexer_options* pinkIndexer_opts,
                             struct felix_options *felix_opts,
@@ -254,7 +255,8 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 
 		case INDEXING_PINKINDEXER :
 		priv = pinkIndexer_prepare(m, cell, pinkIndexer_opts,
-		                           wavelength_estimate);
+		                           wavelength_estimate,
+		                           clen_estimate);
 		break;
 
 		default :
@@ -331,6 +333,7 @@ IndexingPrivate *setup_indexing(const char *method_list,
                                 float *tols,
                                 IndexingFlags flags,
                                 double wavelength_estimate,
+                                double clen_estimate,
                                 int n_threads,
                                 struct taketwo_options *ttopts,
                                 struct xgandalf_options *xgandalf_opts,
@@ -402,6 +405,7 @@ IndexingPrivate *setup_indexing(const char *method_list,
 
 		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
 		                                          wavelength_estimate,
+		                                          clen_estimate,
 		                                          xgandalf_opts,
 		                                          pinkIndexer_opts,
 		                                          felix_opts,
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 88b9b67e..a4aef671 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -171,6 +171,7 @@ extern IndexingPrivate *setup_indexing(const char *methods,
                                        float *ltl,
                                        IndexingFlags flags,
                                        double wavelength_estimate,
+                                       double clen_estimate,
                                        int n_threads,
                                        struct taketwo_options *ttopts,
                                        struct xgandalf_options *xgandalf_opts,
diff --git a/libcrystfel/src/indexers/pinkindexer.c b/libcrystfel/src/indexers/pinkindexer.c
index 149bb5bd..5cb979a9 100644
--- a/libcrystfel/src/indexers/pinkindexer.c
+++ b/libcrystfel/src/indexers/pinkindexer.c
@@ -43,8 +43,6 @@
 #include "peaks.h"
 
 
-#define FAKE_CLEN (0.25)
-
 struct pinkIndexer_options {
 	unsigned int considered_peaks_count;
 	unsigned int angle_resolution;
@@ -81,39 +79,6 @@ static void restoreReciprocalCell(UnitCell *cell, LatticeTransform_t* appliedRed
 static void makeRightHanded(UnitCell* cell);
 
 
-/* Return mean clen in m */
-static double mean_clen(struct detgeom *dg)
-{
-	int i;
-	double total = 0.0;
-	for ( i=0; i<dg->n_panels; i++ ) {
-		total += dg->panels[i].cnz * dg->panels[i].pixel_pitch;
-	}
-	return total / dg->n_panels;
-}
-
-
-static void scale_detector_shift(double fake_clen,
-                                 struct detgeom *dg,
-                                 double inx, double iny,
-                                 double *pdx, double *pdy)
-{
-	int i;
-	double mean = mean_clen(dg);
-	for ( i=0; i<dg->n_panels; i++ ) {
-		if ( !within_tolerance(dg->panels[i].cnz*dg->panels[i].pixel_pitch, mean, 2.0) ) {
-			ERROR("WARNING: Detector is not flat enough to apply "
-			      "detector position offset\n");
-			*pdx = 0.0;
-			*pdy = 0.0;
-			return;
-		}
-	}
-	*pdx = (mean/fake_clen)*inx;
-	*pdy = (mean/fake_clen)*iny;
-}
-
-
 int run_pinkIndexer(struct image *image, void *ipriv, int n_threads)
 {
 	struct pinkIndexer_private_data *pinkIndexer_private_data = ipriv;
@@ -206,19 +171,14 @@ int run_pinkIndexer(struct image *image, void *ipriv, int n_threads)
 			ERROR("pinkIndexer: problem with returned cell!\n");
 		} else {
 
-			double dx, dy;
 			Crystal *cr = crystal_new();
 			if ( cr == NULL ) {
 				ERROR("Failed to allocate crystal.\n");
 				return 0;
 			}
 			crystal_set_cell(cr, new_cell_trans);
-			scale_detector_shift(FAKE_CLEN,
-			                     image->detgeom,
-			                     center_shift[0],
-			                     center_shift[1],
-			                     &dx, &dy);
-			crystal_set_det_shift(cr, dx, dy);
+			crystal_set_det_shift(cr, center_shift[0],
+			                          center_shift[1]);
 			image_add_crystal(image, cr);
 			indexed++;
 
@@ -233,7 +193,8 @@ int run_pinkIndexer(struct image *image, void *ipriv, int n_threads)
 void *pinkIndexer_prepare(IndexingMethod *indm,
                           UnitCell *cell,
                           struct pinkIndexer_options *pinkIndexer_opts,
-                          double wavelength_estimate)
+                          double wavelength_estimate,
+                          double clen_estimate)
 {
 	float beamEenergy_eV;
 
@@ -245,9 +206,14 @@ void *pinkIndexer_prepare(IndexingMethod *indm,
 		beamEenergy_eV = J_to_eV(ph_lambda_to_en(wavelength_estimate));
 	}
 
+	if ( isnan(clen_estimate) ) {
+		ERROR("PinkIndexer requires a camera length estimate.  "
+		      "Try again with --camera-length-estimate=xx\n");
+		return NULL;
+	}
+
 	if ( cell == NULL ) {
-		ERROR("Unit cell information is required for "
-		      "PinkIndexer.\n");
+		ERROR("Unit cell information is required for PinkIndexer.\n");
 		return NULL;
 	}
 
@@ -294,7 +260,7 @@ void *pinkIndexer_prepare(IndexingMethod *indm,
 	Lattice_t sampleReciprocalLattice_1_per_A = lattice;
 	float detectorRadius_m = 0.03; //fake, only for prediction
 	ExperimentSettings *experimentSettings = ExperimentSettings_new(beamEenergy_eV,
-	                                                                FAKE_CLEN,
+	                                                                clen_estimate,
 	                                                                detectorRadius_m,
 	                                                                divergenceAngle_deg,
 	                                                                nonMonochromaticity,
@@ -420,7 +386,8 @@ int run_pinkIndexer(struct image *image, void *ipriv, int n_threads)
 extern void *pinkIndexer_prepare(IndexingMethod *indm,
                                  UnitCell *cell,
                                  struct pinkIndexer_options *pinkIndexer_opts,
-                                 double wavelength_estimate)
+                                 double wavelength_estimate,
+                                 double clen_estimate)
 {
 	ERROR("This copy of CrystFEL was compiled without PINKINDEXER support.\n");
 	ERROR("To use PINKINDEXER indexing, recompile with PINKINDEXER.\n");
diff --git a/libcrystfel/src/indexers/pinkindexer.h b/libcrystfel/src/indexers/pinkindexer.h
index 358a8221..122bac87 100644
--- a/libcrystfel/src/indexers/pinkindexer.h
+++ b/libcrystfel/src/indexers/pinkindexer.h
@@ -42,7 +42,8 @@ extern int run_pinkIndexer(struct image *image, void *ipriv, int n_threads);
 extern void *pinkIndexer_prepare(IndexingMethod *indm,
                                  UnitCell *cell,
                                  struct pinkIndexer_options *pinkIndexer_opts,
-                                 double wavelength_estimate);
+                                 double wavelength_estimate,
+                                 double clen_estimate);
 
 extern void pinkIndexer_cleanup(void *pp);
 
diff --git a/src/gui_index.c b/src/gui_index.c
index 9b72ca10..f5f07f73 100644
--- a/src/gui_index.c
+++ b/src/gui_index.c
@@ -529,7 +529,9 @@ static void run_indexing_once(struct crystfelproject *proj)
 	ipriv = setup_indexing(methods, cell,
 	                       proj->indexing_params.tols,
 	                       indexing_flags(&proj->indexing_params),
-	                       proj->cur_image->lambda, 1,
+	                       proj->cur_image->lambda,
+	                       detgeom_mean_camera_length(proj->cur_image->detgeom),
+	                       1,
 	                       taketwoopts, xgandalf_opts,
 	                       pinkIndexer_opts, felix_opts,
 	                       NULL);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index d5d3d31a..eb39a32c 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -467,6 +467,15 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
 			return EINVAL;
 		}
 		break;
+
+		case 415 :
+		if (sscanf(arg, "%f", &args->iargs.clen_estimate) != 1)
+		{
+			ERROR("Invalid value for --camera-length-estimate\n");
+			return EINVAL;
+		}
+		break;
+
 		/* ---------- Integration ---------- */
 
 		case 501 :
@@ -682,6 +691,7 @@ int main(int argc, char *argv[])
 	args.iargs.no_image_data = 0;
 	args.iargs.no_mask_data = 0;
 	args.iargs.wavelength_estimate = NAN;
+	args.iargs.clen_estimate = NAN;
 	args.iargs.n_threads = 1;
 
 	argp_program_version_hook = show_version;
@@ -784,6 +794,8 @@ int main(int argc, char *argv[])
 		        "Estimate of the incident radiation wavelength, in metres."},
 		{"max-indexer-threads", 414, "n", 0,
 		        "Maximum number of threads allowed for indexing engines."},
+		{"camera-length-estimate", 415, "metres", 0,
+		        "Estimate of the camera length, in metres."},
 
 		{NULL, 0, 0, OPTION_DOC, "Integration options:", 5},
 		{"integration", 501, "method", OPTION_NO_USAGE, "Integration method"},
@@ -996,6 +1008,7 @@ int main(int argc, char *argv[])
 		                                  args.iargs.tols,
 		                                  flags,
 		                                  args.iargs.wavelength_estimate,
+		                                  args.iargs.clen_estimate,
 		                                  args.iargs.n_threads,
 		                                  taketwo_opts,
 		                                  xgandalf_opts,
diff --git a/src/process_image.h b/src/process_image.h
index 703f0cde..c3db8551 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -7,7 +7,7 @@
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
+ *   2010-2021 Thomas White <taw@physics.org>
  *   2014-2017 Valerio Mariani <valerio.mariani@desy.de>
  *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
  *
@@ -99,6 +99,7 @@ struct index_args
 	int no_image_data;
 	int no_mask_data;
 	float wavelength_estimate;
+	float clen_estimate;
 	int n_threads;
 };