Offset peak locations from HDF5 or CXI files by 0.5,0.5

CrystFEL considers all peak locations to be distances from the corner of
the detector panel, in pixel units, consistent with its description of
detector geometry.  In contrast, Cheetah considers the peak locations
to be pixel indices in the data array.  Therefore, a half-pixel offset
is needed when importing the peak lists.

For users who need the old behaviour, this commit adds a new option
indexamajig --no-half-pixel-shift to deactivate this offset.

7 files changed, 156 insertions, 25 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 8619fd92..5855c866 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -1,7 +1,7 @@
 .\"
 .\" indexamajig man page
 .\"
-.\" Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+.\" Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
 .\"                       a research centre of the Helmholtz Association.
 .\"
 .\" Part of CrystFEL - crystallography with a FEL
@@ -45,9 +45,11 @@ You can control the peak detection on the command line.  Firstly, you can choose
 
 \fB--peaks=cxi\fR works similarly to this, but expects four separate HDF5 datasets beneath \fIpath\fR, \fBnPeaks\fR, \fBpeakXPosRaw\fR, \fBpeakYPosRaw\fR and \fBpeakTotalIntensity\fR.  See the specification for the CXI file format at http://www.cxidb.org/ for more details.
 
+CrystFEL considers all peak locations to be distances from the corner of the detector panel, in pixel units, consistent with its description of detector geometry (see 'man crystfel_geometry').  The software which generates the HDF5 or CXI files, including Cheetah, may instead consider the peak locations to be pixel indices in the data array.  Therefore, the peak coordinates from \fB--peaks=cxi\fR or \fB--peaks=hdf5\fR will by default have 0.5 added to them.  Use \fB--no-half-pixel-shift\fR if this isn't what you want.
+
 If you use \fB--peaks=zaef\fR, indexamajig will use a simple gradient search after Zaefferer (2000).  You can control the overall threshold and minimum squared gradient for finding a peak using \fB--threshold\fR and \fB--min-gradient\fR.  The threshold has arbitrary units matching the pixel values in the data, and the minimum gradient has the equivalent squared units.  Peaks will be rejected if the 'foot point' is further away from the 'summit' of the peak by more than the inner integration radius (see below).  They will also be rejected if the peak is closer than twice the inner integration radius from another peak.
 
-If you instead use \fB--peaks=peakfinder8\fR, indexamajig will user the "peakfinder8" peak finding algorithm describerd in Barty et al. (2014). Pixels above a radius-dependent intensity threshold are considered as candidate peaks (although the user sets an absolute minimum threshold for candidate peaks). Peaks are then only accepted if their signal to noise level over the local background is sufficiently high. Peaks can include multiple pixels and the user can reject a peak if it includes too many or too few. The distance of a peak from the center of the detector can also be used as a filtering criterion. Please notice that the peakfinder8 will not report more than 2048 peaks for each panel: any additional peak is ignored.
+If you instead use \fB--peaks=peakfinder8\fR, indexamajig will use the "peakfinder8" peak finding algorithm describerd in Barty et al. (2014). Pixels above a radius-dependent intensity threshold are considered as candidate peaks (although the user sets an absolute minimum threshold for candidate peaks). Peaks are then only accepted if their signal to noise level over the local background is sufficiently high. Peaks can include multiple pixels and the user can reject a peak if it includes too many or too few. The distance of a peak from the center of the detector can also be used as a filtering criterion. Note that the peakfinder8 will not report more than 2048 peaks for each panel: any additional peak is ignored.
 
 You can suppress peak detection altogether for a panel in the geometry file by specifying the "no_index" value for the panel as non-zero.
 
@@ -335,6 +337,11 @@ Normally, peaks which contain one or more pixels above max_adu (defined in the d
 When using \fB--peaks=hdf5\fR or \fB--peaks=cxi\fR, the peaks will be put through some of 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), 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--no-half-pixel-shift\fR
+.PD
+CrystFEL considers all peak locations to be distances from the corner of the detector panel, in pixel units, consistent with its description of detector geometry (see 'man crystfel_geometry').  The software which generates the HDF5 or CXI files, including Cheetah, may instead consider the peak locations to be pixel indices in the data array.  Therefore, the peak coordinates from \fB--peaks=cxi\fR or \fB--peaks=hdf5\fR will by default have 0.5 added to them.  This option \fBdisables\fR this half-pixel offset.
+
+.PD 0
 .IP \fB--check-hdf5-snr\fR
 .PD
 With this option with \fB--peaks=hdf5\fR, the peaks will additionally be checked to see that they satisfy the minimum SNR specified with \fB--min-snr\fR.
diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt
index 4cd922d2..d58c5725 100644
--- a/doc/reference/libcrystfel/CrystFEL-sections.txt
+++ b/doc/reference/libcrystfel/CrystFEL-sections.txt
@@ -439,7 +439,9 @@ add_copy_hdf5_field
 new_copy_hdf5_field_list
 free_copy_hdf5_field_list
 get_peaks
+get_peaks_2
 get_peaks_cxi
+get_peaks_cxi_2
 hdfile_is_scalar
 check_path_existence
 fill_event_list
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index 9bcdb1b7..76617ee4 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -44,6 +44,19 @@
 #include "utils.h"
 
 
+/**
+ * SECTION:hdf5-file
+ * @short_description: HDF5 file handling
+ * @title: HDF5 file handling
+ * @section_id:
+ * @see_also:
+ * @include: "hdf5-file.h"
+ * @Image:
+ *
+ * Routines for accessing HDF5 files.
+ **/
+
+
 struct hdf5_write_location {
 
 	const char      *location;
@@ -320,9 +333,34 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line)
 }
 
 
-/* Get peaks from HDF5, in "CXI format" (as in "CXIDB") */
-int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
-                  struct filename_plus_event *fpe)
+/**
+ * get_peaks_cxi_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5, in "CXI format" (as in "CXIDB").  The data should be in
+ * a set of arrays under @p.  The number of peaks should be in a 1D array at
+ * @p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at
+ * @p/peakXPosRaw and @p/peakYPosRaw respectively (sorry about the naming).  The
+ * first dimension of these arrays should be the event number (as given by
+ * @fpe).  The intensities are expected to be at @p/peakTotalIntensity in a
+ * similar 2D array.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry').  The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array.  In this case, the peak coordinates should
+ * have 0.5 added to them.  This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+                    struct filename_plus_event *fpe, int half_pixel_shift)
 {
 	char path_n[1024];
 	char path_x[1024];
@@ -338,6 +376,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 	float *buf_y;
 	float *buf_i;
 
+	double peak_offset = half_pixel_shift ? 0.5 : 0.0;
+
 	if ( (fpe != NULL) && (fpe->ev != NULL)
 	  && (fpe->ev->dim_entries != NULL) )
 	{
@@ -375,8 +415,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 		float fs, ss, val;
 		struct panel *p;
 
-		fs = buf_x[pk];
-		ss = buf_y[pk];
+		fs = buf_x[pk] + peak_offset;
+		ss = buf_y[pk] + peak_offset;
 		val = buf_i[pk];
 
 		p = find_orig_panel(image->det, fs, ss);
@@ -395,7 +435,53 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
 }
 
 
-int get_peaks(struct image *image, struct hdfile *f, const char *p)
+/**
+ * get_peaks_cxi:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions.  Use get_peaks_cxi_2() instead.
+ *
+ * This function is equivalent to get_peaks_cxi_2(@image, @f, @p, @fpe, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
+                  struct filename_plus_event *fpe)
+{
+	return get_peaks_cxi_2(image, f, p, fpe, 1);
+}
+
+
+/**
+ * get_peaks_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5.  The peak list should be located at @p in the HDF5 file,
+ * a 2D array where the first dimension equals the number of peaks and second
+ * dimension is three.  The first two columns contain the fast scan and slow
+ * scan coordinates, respectively, of the peaks.  The third column contains the
+ * intensities.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry').  The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array.  In this case, the peak coordinates should
+ * have 0.5 added to them.  This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+                int half_pixel_shift)
 {
 	hid_t dh, sh;
 	hsize_t size[2];
@@ -405,6 +491,7 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 	herr_t r;
 	int tw;
 	char *np;
+	double peak_offset = half_pixel_shift ? 0.5 : 0.0;
 
 	if ( image->event != NULL ) {
 		np = retrieve_full_path(image->event, p);
@@ -473,8 +560,8 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 		float fs, ss, val;
 		struct panel *p;
 
-		fs = buf[tw*i+0];
-		ss = buf[tw*i+1];
+		fs = buf[tw*i+0] + peak_offset;
+		ss = buf[tw*i+1] + peak_offset;
 		val = buf[tw*i+2];
 
 		p = find_orig_panel(image->det, fs, ss);
@@ -499,6 +586,26 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
 }
 
 
+/**
+ * get_peaks:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions.  Use get_peaks_2() instead.
+ *
+ * This function is equivalent to get_peaks_2(@image, @f, @p, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks(struct image *image, struct hdfile *f, const char *p)
+{
+	return get_peaks_2(image, f, p, 1);
+}
+
+
 static void cleanup(hid_t fh)
 {
 	int n_ids, i;
diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h
index 8c89eb93..f22d4d54 100644
--- a/libcrystfel/src/hdf5-file.h
+++ b/libcrystfel/src/hdf5-file.h
@@ -3,11 +3,11 @@
  *
  * Read/write HDF5 data files
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2015 Thomas White <taw@physics.org>
+ *   2009-2017 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
 
  *
@@ -78,9 +78,16 @@ extern void hdfile_close(struct hdfile *f);
 
 extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
 
+extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+                       int half_pixel_shift);
+
 extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
                          struct filename_plus_event *fpe);
 
+extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+                           struct filename_plus_event *fpe,
+                           int half_pixel_shift);
+
 extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
 extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 28a7bf85..f48fac17 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -152,16 +152,17 @@ static void show_help(const char *s)
 " --temp-dir=<path>        Put the temporary folder under <path>.\n"
 "\n"
 "\nOptions you probably won't need:\n\n"
-"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
-"     --no-use-saturated   During the initial peak search, reject\n"
-"                           peaks which contain pixels above max_adu.\n"
-"     --no-revalidate      Don't re-integrate and check HDF5 peaks for\n"
-"                           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=<cond>    Show debugging information about reflections.\n"
-"     --no-refine          Skip the prediction refinement step.\n"
-"     --profile            Show timing data for performance monitoring.\n"
+"     --no-check-prefix     Don't attempt to correct the --prefix.\n"
+"     --no-use-saturated    During the initial peak search, reject\n"
+"                            peaks which contain pixels above max_adu.\n"
+"     --no-revalidate       Don't re-integrate and check HDF5 peaks for\n"
+"                            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=<cond>     Show debugging information about reflections.\n"
+"     --no-refine           Skip the prediction refinement step.\n"
+"     --profile             Show timing data for performance monitoring.\n"
+"     --no-half-pixel-shift Don't offset the HDF5 peak locations by 0.5 px.\n"
 "\nLow-level options for the felix indexer:\n\n"
 "     --felix-options      Change the default arguments passed to the indexer.\n"
 "                          Given as a list of comma separated list of \n"
@@ -243,6 +244,7 @@ int main(int argc, char *argv[])
 	iargs.peaks = PEAK_ZAEF;
 	iargs.beam = &beam;
 	iargs.hdf5_peak_path = NULL;
+	iargs.half_pixel_shift = 1;
 	iargs.copyme = NULL;
 	iargs.pk_inn = -1.0;
 	iargs.pk_mid = -1.0;
@@ -297,6 +299,7 @@ int main(int argc, char *argv[])
 		{"check-hdf5-snr",     0, &iargs.check_hdf5_snr,     1},
 		{"no-refine",          0, &no_refine,                1},
 		{"profile",            0, &iargs.profile,            1},
+		{"no-half-pixel-shift",0, &iargs.half_pixel_shift,   0},
 
 		/* Long-only options which don't actually do anything */
 		{"no-sat-corr",        0, &iargs.satcorr,            0},
diff --git a/src/process_image.c b/src/process_image.c
index e6e6d22a..cd5c0866 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -160,7 +160,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	switch ( iargs->peaks ) {
 
 		case PEAK_HDF5:
-		if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) {
+		if ( get_peaks_2(&image, hdfile, iargs->hdf5_peak_path,
+		                 iargs->half_pixel_shift) )
+		{
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
@@ -172,8 +174,10 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		break;
 
 		case PEAK_CXI:
-		if ( get_peaks_cxi(&image, hdfile, iargs->hdf5_peak_path,
-		                   pargs->filename_p_e) ) {
+		if ( get_peaks_cxi_2(&image, hdfile, iargs->hdf5_peak_path,
+		                     pargs->filename_p_e,
+		                     iargs->half_pixel_shift) )
+		{
 			ERROR("Failed to get peaks from CXI file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
diff --git a/src/process_image.h b/src/process_image.h
index bc8b31fd..d63679cf 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -67,6 +67,7 @@ struct index_args
 	float tols[4];
 	struct beam_params *beam;
 	char *hdf5_peak_path;
+	int half_pixel_shift;
 	float pk_inn;
 	float pk_mid;
 	float pk_out;