Merge branch 'tom/peakfinder9'

8 files changed, 313 insertions, 111 deletions

diff --git a/CMake/FindFDIP.cmake b/CMake/FindFDIP.cmake
index 6db3b993..3306b483 100644
--- a/CMake/FindFDIP.cmake
+++ b/CMake/FindFDIP.cmake
@@ -10,16 +10,16 @@ if (FDIP_INCLUDES)
   set (FDIP_FIND_QUIETLY TRUE)
 endif (FDIP_INCLUDES)
 
-find_path (FDIP_INCLUDES streakFinder.h
+find_path (FDIP_INCLUDES fastDiffractionImageProcessing/streakFinder.h
            PATHS
-		   ${CMAKE_INSTALL_PREFIX}/include
-           PATH_SUFFIXES fastDiffractionImageProcessing)
+           ${CMAKE_INSTALL_PREFIX}/include)
 
 find_library (FDIP_LIBRARIES fastDiffractionImageProcessing
               PATHS
-		      ${CMAKE_INSTALL_PREFIX}/lib)
+              ${CMAKE_INSTALL_PREFIX}/lib)
 
 # handle the QUIETLY and REQUIRED arguments and set FDIP_FOUND to TRUE if
 # all listed variables are TRUE
 include (FindPackageHandleStandardArgs)
+set(FDIP_INCLUDES ${FDIP_INCLUDES};${FDIP_INCLUDES}/fastDiffractionImageProcessing)
 find_package_handle_standard_args (FDIP DEFAULT_MSG FDIP_LIBRARIES FDIP_INCLUDES)
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index c0fe8277..37e80801 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -51,6 +51,8 @@ If you use \fB--peaks=zaef\fR, indexamajig will use a simple gradient search aft
 
 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.
 
+If you instead use \fB--peaks=peakfinder9\fR, indexamajig will use the "peakfinder9" peak finding algorithm described in the master thesis "Real-time image analysis and data compression in high throughput X-ray diffraction experiments" by Gevorkov. Other than peakFinder8, peakFinder9 uses local background estimation based on border pixels in a specified radius (\fB--local-bg-radius\fR). For being fast and precise, a hierarchy of conditions is used. First condition is only useful for speed consideration, it demands that a pixel that is the biggest pixel in a peak must be larger than every border pixel by a constant value (\fB--min-peak-over-neighbour\fR). Second condition ensures, that the pixel passing the previous condition is the highest pixel in the peak. It assumes, that peaks rise monotonically towards the biggest pixel. Third condition ensures, that the biggest pixel in the peak is significantly over the noise level (\fB--min-snr-biggest-pix\fR) by computing the local statistics from the border pixels in a specified radius. Fourth condition sums up all pixels belonging to the peak (\fB--min-snr-peak-pix\fR) and demands that the whole peak must be significantly over the noise level (\fB--min-snr\fR). Only if all conditions are passed, the peak is accepted.
+
 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.
 
 
@@ -254,7 +256,27 @@ Set the square of the gradient threshold for peak detection using \fB--peaks=zae
 .PD 0
 .IP \fB--min-snr=\fR\fIsnr\fR
 .PD
-Set the minimum I/sigma(I) for peak detection when using \fB--peaks=zaef\fR or \fB--peaks=peakfinder8\fR.  The default is \fB--min-snr=5\fR.
+Set the minimum I/sigma(I) for peak detection when using \fB--peaks=zaef\fR, \fB--peaks=peakfinder8\fR or \fB--peaks=peakfinder9\fR.  The default is \fB--min-snr=5\fR.
+
+.PD 0
+.IP \fB--min-snr-biggest-pix=<n>\fR
+.PD
+(peakFinder9 only) min snr of the biggest pixel in the peak, given as a factor of the standard deviation. Default is 7.0.
+
+.PD 0
+.IP \fB--min-snr-peak-pix=<n>\fR
+.PD
+(peakFinder9 only) min snr of a peak pixel, given as a factor of the standard deviation. Should be smaller or equal to sig_fac_biggest_pix. Default is 6.0.
+
+.PD 0
+.IP \fB--min-sig=<n>\fR
+.PD
+(peakFinder9 only) minimum standard deviation of the background. Prevents finding of peaks in erroneous or highly shadowed unmasked regions. Default is 11.0.
+
+.PD 0
+.IP \fB--min-peak-over-neighbour=<n>\fR
+.PD
+(peakFinder9 only) just for speed. Biggest pixel must be n higher than the pixels in window_radius distance to be a candidate for the biggest pixel in a peak. Should be chosen as a small positive number, a few times smaller than the weakest expected peak. The default is -INFINITY, which turns off the speedup and searches with maximum precision.
 
 .PD 0
 .IP \fB--min-pix-count=\fR\fIcnt\fR
@@ -269,7 +291,7 @@ Accepts peaks only if they include less than \fR\fIcnt\fR pixels, when using \fB
 .PD 0
 .IP \fB--local-bg-radius=\fR\fIr\fR
 .PD
-Radius (in pixels) used for the estimation of the local background when using \fB--peaks=peakfinder8\fR.  The default is \fB--local-bg-radius=3\fR.
+Radius (in pixels) used for the estimation of the local background when using \fB--peaks=peakfinder8 or --peaks=peakfinder9\fR.  The default is \fB--local-bg-radius=3\fR.
 
 .PD 0
 .IP \fB--min-res=\fR\fIpx\fR
@@ -495,7 +517,7 @@ For a full explanation of how the internal layout of the data file can be  descr
 You can use \fBlist_events\fR to prepare a list of each event in one or more input files.  Note that you only need to do this if you need to perform some sorting or filtering on this list.  If you want to process every event in a file, simply specify the filename in the input file.
 
 .SH AUTHOR
-This page was written by Thomas White.
+This page was written by Thomas White, Yaroslav Gevorkov and Valerio Mariani.
 
 .SH REPORTING BUGS
 Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index 28c79538..6c48bd3a 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -3,7 +3,7 @@
  *
  * Peak search and other image analysis
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
@@ -13,6 +13,7 @@
  *   2011      Andrew Martin <andrew.martin@desy.de>
  *   2011      Richard Kirian
  *   2017      Valerio Mariani <valerio.mariani@desy.de>
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -45,6 +46,12 @@
 #include <pthread.h>
 #include <fenv.h>
 
+#ifdef HAVE_FDIP
+#include "fastDiffractionImageProcessing/adaptions/crystfel/peakFinder9.h"
+#include "fastDiffractionImageProcessing/adaptions/crystfel/mask.h"
+#include "fastDiffractionImageProcessing/peakList.h"
+#endif
+
 #include "image.h"
 #include "utils.h"
 #include "peaks.h"
@@ -563,6 +570,107 @@ int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
 }
 
 
+#ifdef HAVE_FDIP
+
+int search_peaks_peakfinder9(struct image *image, float min_snr_biggest_pix,
+                             float min_snr_peak_pix, float min_snr_whole_peak,
+                             float min_sig, float min_peak_over_neighbour,
+                             int window_radius)
+{
+	peakFinder9_accuracyConstants_t accuracy_consts;
+	peakList_t peakList;
+	long NpeaksMax = 10000; //more peaks per panel should not appear
+	float *data_copy = NULL;
+	float *data_copy_new;
+	int panel_number;
+
+	if ( image->features != NULL ) {
+		image_feature_list_free(image->features);
+	}
+	image->features = image_feature_list_new();
+	image->num_peaks = 0;
+	image->num_saturated_peaks = 0;
+
+	accuracy_consts.minSNR_biggestPixel = min_snr_biggest_pix;
+	accuracy_consts.minSNR_peakPixel = min_snr_peak_pix;
+	accuracy_consts.minSNR_wholePeak = min_snr_whole_peak;
+	accuracy_consts.minimumSigma = min_sig;
+	accuracy_consts.minimumPeakOversizeOverNeighbours = min_peak_over_neighbour;
+	accuracy_consts.windowRadius = window_radius;
+
+	if ( allocatePeakList(&peakList, NpeaksMax) ) return 1;
+
+	for ( panel_number=0; panel_number<image->det->n_panels; panel_number++ ) {
+
+		int w, h;
+		int peak_number;
+		detectorRawFormat_t det_size_one_panel;
+
+		if ( image->det->panels[panel_number].no_index ) continue;
+
+		w = image->det->panels[panel_number].w;
+		h = image->det->panels[panel_number].h;
+
+		det_size_one_panel.asic_nx = w;
+		det_size_one_panel.asic_ny = h;
+		det_size_one_panel.nasics_x = 1;
+		det_size_one_panel.nasics_y = 1;
+		det_size_one_panel.pix_nx = w;
+		det_size_one_panel.pix_ny = h;
+		det_size_one_panel.pix_nn = w * h;
+
+		data_copy_new = realloc(data_copy, w*h*sizeof(*data_copy));
+		if ( data_copy_new == NULL ) {
+			if ( data_copy != NULL ) {
+				free(data_copy);
+			}
+			freePeakList(peakList);
+			return 1;
+		} else {
+			data_copy = data_copy_new;
+		}
+
+		mergeMaskAndDataIntoDataCopy(image->dp[panel_number], data_copy,
+		                             image->bad[panel_number],
+		                             &det_size_one_panel);
+
+		peakList.peakCount = 0;
+		image->num_peaks += peakFinder9_onePanel_noSlab(data_copy,
+		                                                &accuracy_consts,
+		                                                &det_size_one_panel,
+		                                                &peakList);
+
+		for ( peak_number=0; peak_number<peakList.peakCount; peak_number++) {
+			image_add_feature(image->features,
+			                  peakList.centerOfMass_rawX[peak_number],
+			                  peakList.centerOfMass_rawY[peak_number],
+			                  &image->det->panels[panel_number],
+			                  image,
+			                  peakList.totalIntensity[peak_number],
+			                  NULL);
+		}
+
+	}
+
+	freePeakList(peakList);
+	free(data_copy);
+	return 0;
+}
+
+#else
+
+int search_peaks_peakfinder9(struct image *image, float min_snr_biggest_pix,
+                             float min_snr_peak_pix, float min_snr_whole_peak,
+                             float min_sig, float min_peak_over_neighbour,
+                             int window_radius)
+{
+	ERROR("This copy of CrystFEL was compiled without peakfinder9 support.\n");
+	return 1;
+}
+
+#endif // HAVE_FDIP
+
+
 int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst)
 {
 	int n_feat = 0;
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index a5095127..3abef729 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -3,12 +3,13 @@
  *
  * Peak search and other image analysis
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2010-2015 Thomas White <taw@physics.org>
  *   2017      Valerio Mariani <valerio.mariani@desy.de>
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -55,6 +56,13 @@ extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
                                     int local_bg_radius, int min_res,
                                     int max_res, int use_saturated);
 
+extern int search_peaks_peakfinder9(struct image *image,
+                                    float min_snr_biggest_pix,
+                                    float min_snr_peak_pix,
+                                    float min_snr_whole_peak, float min_sig,
+                                    float min_peak_over_neighbour,
+                                    int window_radius);
+
 extern int peak_sanity_check(struct image *image, Crystal **crystals,
                              int n_cryst);
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 38d83d58..4a790a17 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,7 +3,7 @@
  *
  * Index patterns, output hkl+intensity etc.
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  * Copyright © 2012 Lorenzo Galli
@@ -14,6 +14,7 @@
  *   2012      Lorenzo Galli
  *   2012      Chunhong Yoon
  *   2017      Valerio Mariani <valerio.mariani@desy.de>
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -84,7 +85,7 @@ static void show_help(const char *s)
 "     --profile             Show timing data for performance monitoring\n"
 "     --temp-dir=<path>     Put the temporary folder under <path>\n"
 "\nPeak search options:\n\n"
-"     --peaks=<method>      Peak search method (zaef,peakfinder8,hdf5,cxi)\n"
+"     --peaks=<method>      Peak search method (zaef,peakfinder8,peakfinder9,hdf5,cxi)\n"
 "                            Default: zaef\n"
 "     --peak-radius=<r>     Integration radii for peak search\n"
 "     --min-peaks=<n>       Minimum number of peaks for indexing\n"
@@ -98,17 +99,24 @@ static void show_help(const char *s)
 "     --min-gradient=<n>    Minimum squared gradient\n"
 "                            (zaef only) Default: 100,000\n"
 "     --min-snr=<n>         Minimum signal/noise ratio for peaks\n"
-"                            (zaef,pekafinder8 only) Default: 5\n"
+"                            (zaef,peakfinder8, peakfinder9 only) Default: 5\n"
 "     --min-pix-count=<n>   Minimum number of pixels per peak\n"
 "                            (peakfinder8 only) Default: 2\n"
 "     --max-pix-count=<n>   Maximum number of pixels per peak\n"
 "                            (peakfinder8 only) Default: 200\n"
 "     --local-bg-radius=<n> Radius (pixels) for local background estimation\n"
-"                            (peakfinder8 only) Default: 3\n"
+"                            (peakfinder8, peakfinder9 only) Default: 3\n"
 "     --min-res=<n>         Minimum resolution for peak search (in pixels)\n"
 "                            (peakfinder8 only) Default: 0\n"
 "     --max-res=<n>         Maximum resolution for peak search (in pixels)\n"
 "                            (peakfinder8 only) Default: 1200\n"
+"     --min-snr-biggest-pix=<n>    (peakFinder9 only) min snr of the biggest pixel in "
+"                                  the peak\n"
+"     --min-snr-peak-pix=<n>       (peakFinder9 only) min snr of a peak pixel\n"
+"     --min-sig=<n>                (peakFinder9 only) minimum standard deviation of "
+"                                 the background\n"
+"     --min-peak-over-neighbour=<n>    (peakFinder9 only) just for speed. Biggest pixel"
+"                                     in peak must be n higher than this.\n"
 "     --no-use-saturated    Reject saturated peaks\n"
 "     --no-revalidate       Don't re-integrate and check HDF5 peaks\n"
 "     --no-half-pixel-shift\n"
@@ -260,6 +268,10 @@ int main(int argc, char *argv[])
 	iargs.min_res = 0;
 	iargs.max_res = 1200;
 	iargs.local_bg_radius = 3;
+	iargs.min_snr_biggest_pix = 7.0;    /* peak finder 9  */
+	iargs.min_snr_peak_pix = 6.0;
+	iargs.min_sig = 11.0;
+	iargs.min_peak_over_neighbour = -INFINITY;
 	iargs.check_hdf5_snr = 0;
 	iargs.det = NULL;
 	iargs.peaks = PEAK_ZAEF;
@@ -356,55 +368,59 @@ int main(int argc, char *argv[])
 		{"use-saturated",      0, &iargs.use_saturated,      1},
 
 		/* Long-only options with arguments */
-		{"peaks",              1, NULL,                2},
-		{"cell-reduction",     1, NULL,                3},
-		{"min-gradient",       1, NULL,                4},
-		{"record",             1, NULL,                5},
-		{"cpus",               1, NULL,                6},
-		{"cpugroup",           1, NULL,                7},
-		{"cpuoffset",          1, NULL,                8},
-		{"hdf5-peaks",         1, NULL,                9},
-		{"copy-hdf5-field",    1, NULL,               10},
-		{"min-snr",            1, NULL,               11},
-		{"tolerance",          1, NULL,               13},
-		{"int-radius",         1, NULL,               14},
-		{"median-filter",      1, NULL,               15},
-		{"integration",        1, NULL,               16},
-		{"temp-dir",           1, NULL,               17},
-		{"int-diag",           1, NULL,               18},
-		{"push-res",           1, NULL,               19},
-		{"res-push",           1, NULL,               19}, /* compat */
-		{"peak-radius",        1, NULL,               20},
-		{"highres",            1, NULL,               21},
-		{"fix-profile-radius", 1, NULL,               22},
-		{"fix-bandwidth",      1, NULL,               23},
-		{"fix-divergence",     1, NULL,               24},
-		{"felix-options",      1, NULL,               25},
-		{"min-pix-count",      1, NULL,               26},
-		{"max-pix-count",      1, NULL,               27},
-		{"local-bg-radius",    1, NULL,               28},
-		{"min-res",            1, NULL,               29},
-		{"max-res",            1, NULL,               30},
-		{"min-peaks",          1, NULL,               31},
-		{"taketwo-member-threshold", 1, NULL,         32},
-		{"taketwo-member-thresh",    1, NULL,         32}, /* compat */
-		{"taketwo-len-tolerance",    1, NULL,         33},
-		{"taketwo-len-tol",          1, NULL,         33}, /* compat */
-		{"taketwo-angle-tolerance",  1, NULL,         34},
-		{"taketwo-angle-tol",        1, NULL,         34}, /* compat */
-		{"taketwo-trace-tolerance",  1, NULL,         35},
-		{"taketwo-trace-tol",        1, NULL,         35}, /* compat */
-		{"felix-tthrange-min",       1, NULL,         36},
-		{"felix-tthrange-max",       1, NULL,         37},
-		{"felix-min-visits",         1, NULL,         38},
-		{"felix-min-completeness",   1, NULL,         39},
-		{"felix-max-uniqueness",     1, NULL,         40},
-		{"felix-num-voxels",         1, NULL,         41},
-		{"felix-fraction-max-visits",1, NULL,         42},
-		{"felix-sigma",              1, NULL,         43},
-		{"serial-start",             1, NULL,         44},
-		{"felix-domega",             1, NULL,         45},
-		{"felix-max-internal-angle", 1, NULL,         46},
+		{"peaks",              1, NULL,              302},
+		{"cell-reduction",     1, NULL,              303},
+		{"min-gradient",       1, NULL,              304},
+		{"record",             1, NULL,              305},
+		{"cpus",               1, NULL,              306},
+		{"cpugroup",           1, NULL,              307},
+		{"cpuoffset",          1, NULL,              308},
+		{"hdf5-peaks",         1, NULL,              309},
+		{"copy-hdf5-field",    1, NULL,              310},
+		{"min-snr",            1, NULL,              311},
+		{"tolerance",          1, NULL,              313},
+		{"int-radius",         1, NULL,              314},
+		{"median-filter",      1, NULL,              315},
+		{"integration",        1, NULL,              316},
+		{"temp-dir",           1, NULL,              317},
+		{"int-diag",           1, NULL,              318},
+		{"push-res",           1, NULL,              319},
+		{"res-push",           1, NULL,              319}, /* compat */
+		{"peak-radius",        1, NULL,              320},
+		{"highres",            1, NULL,              321},
+		{"fix-profile-radius", 1, NULL,              322},
+		{"fix-bandwidth",      1, NULL,              323},
+		{"fix-divergence",     1, NULL,              324},
+		{"felix-options",      1, NULL,              325},
+		{"min-pix-count",      1, NULL,              326},
+		{"max-pix-count",      1, NULL,              327},
+		{"local-bg-radius",    1, NULL,              328},
+		{"min-res",            1, NULL,              329},
+		{"max-res",            1, NULL,              330},
+		{"min-peaks",          1, NULL,              331},
+		{"taketwo-member-threshold", 1, NULL,        332},
+		{"taketwo-member-thresh",    1, NULL,        332}, /* compat */
+		{"taketwo-len-tolerance",    1, NULL,        333},
+		{"taketwo-len-tol",          1, NULL,        333}, /* compat */
+		{"taketwo-angle-tolerance",  1, NULL,        334},
+		{"taketwo-angle-tol",        1, NULL,        334}, /* compat */
+		{"taketwo-trace-tolerance",  1, NULL,        335},
+		{"taketwo-trace-tol",        1, NULL,        335}, /* compat */
+		{"felix-tthrange-min",       1, NULL,        336},
+		{"felix-tthrange-max",       1, NULL,        337},
+		{"felix-min-visits",         1, NULL,        338},
+		{"felix-min-completeness",   1, NULL,        339},
+		{"felix-max-uniqueness",     1, NULL,        340},
+		{"felix-num-voxels",         1, NULL,        341},
+		{"felix-fraction-max-visits",1, NULL,        342},
+		{"felix-sigma",              1, NULL,        343},
+		{"serial-start",             1, NULL,        344},
+		{"felix-domega",             1, NULL,        345},
+		{"felix-max-internal-angle", 1, NULL,        346},
+		{"min-snr-biggest-pix"              ,1, NULL,347},
+		{"min-snr-peak-pix"                 ,1, NULL,348},
+		{"min-sig"                          ,1, NULL,349},
+		{"min-peak-over-neighbour"          ,1, NULL,350},
 
 		{0, 0, NULL, 0}
 	};
@@ -462,11 +478,11 @@ int main(int argc, char *argv[])
 			iargs.threshold = strtof(optarg, NULL);
 			break;
 
-			case 2 :
+			case 302 :
 			speaks = strdup(optarg);
 			break;
 
-			case 3 :
+			case 303 :
 			ERROR("The option '--cell-reduction' is no longer "
 			      "used.\n"
 			      "The complete indexing behaviour is now "
@@ -475,62 +491,62 @@ int main(int argc, char *argv[])
 			      "available methods.\n");
 			return 1;
 
-			case 4 :
+			case 304 :
 			iargs.min_gradient = strtof(optarg, NULL);
 			break;
 
-			case 5 :
+			case 305 :
 			ERROR("The option '--record' is no longer used.\n"
 			      "Use '--no-peaks-in-stream' and"
 			      "'--no-refls-in-stream' if you need to control"
 			      "the contents of the stream.\n");
 			return 1;
 
-			case 6 :
-			case 7 :
-			case 8 :
+			case 306 :
+			case 307 :
+			case 308 :
 			ERROR("The options --cpus, --cpugroup and --cpuoffset"
 			      " are no longer used by indexamajig.\n");
 			break;
 
-			case 9 :
+			case 309 :
 			free(command_line_peak_path);
 			command_line_peak_path = strdup(optarg);
 			break;
 
-			case 10 :
+			case 310 :
 			add_imagefile_field(iargs.copyme, optarg);
 			break;
 
-			case 11 :
+			case 311 :
 			iargs.min_snr = strtof(optarg, NULL);
 			break;
 
-			case 13 :
+			case 313 :
 			toler = strdup(optarg);
 			break;
 
-			case 14 :
+			case 314 :
 			intrad = strdup(optarg);
 			break;
 
-			case 15 :
+			case 315 :
 			iargs.median_filter = atoi(optarg);
 			break;
 
-			case 16 :
+			case 316 :
 			int_str = strdup(optarg);
 			break;
 
-			case 17 :
+			case 317 :
 			tempdir = strdup(optarg);
 			break;
 
-			case 18 :
+			case 318 :
 			int_diag = strdup(optarg);
 			break;
 
-			case 19 :
+			case 319 :
 			if ( sscanf(optarg, "%f", &iargs.push_res) != 1 ) {
 				ERROR("Invalid value for --push-res\n");
 				return 1;
@@ -539,11 +555,11 @@ int main(int argc, char *argv[])
 			have_push_res = 1;
 			break;
 
-			case 20 :
+			case 320 :
 			pkrad = strdup(optarg);
 			break;
 
-			case 21 :
+			case 321 :
 			if ( sscanf(optarg, "%f", &iargs.highres) != 1 ) {
 				ERROR("Invalid value for --highres\n");
 				return 1;
@@ -552,7 +568,7 @@ int main(int argc, char *argv[])
 			iargs.highres = 1.0 / (iargs.highres/1e10);
 			break;
 
-			case 22 :
+			case 322 :
 			if ( sscanf(optarg, "%f", &iargs.fix_profile_r) != 1 ) {
 				ERROR("Invalid value for "
 				      "--fix-profile-radius\n");
@@ -560,50 +576,50 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 23 :
+			case 323 :
 			if ( sscanf(optarg, "%f", &iargs.fix_bandwidth) != 1 ) {
 				ERROR("Invalid value for --fix-bandwidth\n");
 				return 1;
 			}
 			break;
 
-			case 24 :
+			case 324 :
 			if ( sscanf(optarg, "%f", &iargs.fix_divergence) != 1 ) {
 				ERROR("Invalid value for --fix-divergence\n");
 				return 1;
 			}
 			break;
 
-			case 25 :
+			case 325 :
 			ERROR("--felix-options is no longer used.\n");
 			ERROR("See --help for the new Felix options.\n");
 			return 1;
 
-		        case 26:
+			case 326:
 			iargs.min_pix_count = atoi(optarg);
 			break;
 
-		        case 27:
+			case 327:
 			iargs.max_pix_count = atoi(optarg);
 			break;
 
-		        case 28:
+			case 328:
 			iargs.local_bg_radius = atoi(optarg);
 			break;
 
-		        case 29:
+			case 329:
 			iargs.min_res = atoi(optarg);
 			break;
 
-		        case 30:
+			case 330:
 			iargs.max_res = atoi(optarg);
 			break;
 
-			case 31:
+			case 331:
 			iargs.min_peaks = atoi(optarg);
 			break;
 
-			case 32:
+			case 332:
 			if ( sscanf(optarg, "%i", &iargs.taketwo_opts.member_thresh) != 1 )
 			{
 				ERROR("Invalid value for --taketwo-member-threshold\n");
@@ -611,7 +627,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 33:
+			case 333:
 			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.len_tol) != 1 )
 			{
 				ERROR("Invalid value for --taketwo-len-tolerance\n");
@@ -621,7 +637,7 @@ int main(int argc, char *argv[])
 			iargs.taketwo_opts.len_tol *= 1e10;
 			break;
 
-			case 34:
+			case 334:
 			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.angle_tol) != 1 )
 			{
 				ERROR("Invalid value for --taketwo-angle-tolerance\n");
@@ -631,7 +647,7 @@ int main(int argc, char *argv[])
 			iargs.taketwo_opts.angle_tol = deg2rad(iargs.taketwo_opts.angle_tol);
 			break;
 
-			case 35:
+			case 335:
 			if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.trace_tol) != 1 )
 			{
 				ERROR("Invalid value for --taketwo-trace-tolerance\n");
@@ -641,7 +657,7 @@ int main(int argc, char *argv[])
 			iargs.taketwo_opts.trace_tol = deg2rad(iargs.taketwo_opts.trace_tol);
 			break;
 
-			case 36:
+			case 336:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.ttmin) != 1 )
 			{
 				ERROR("Invalid value for --felix-tthrange-min\n");
@@ -650,7 +666,7 @@ int main(int argc, char *argv[])
 			iargs.felix_opts.ttmin = deg2rad(iargs.felix_opts.ttmin);
 			break;
 
-			case 37:
+			case 337:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.ttmax) != 1 )
 			{
 				ERROR("Invalid value for --felix-tthrange-max\n");
@@ -659,7 +675,7 @@ int main(int argc, char *argv[])
 			iargs.felix_opts.ttmax = deg2rad(iargs.felix_opts.ttmax);
 			break;
 
-			case 38:
+			case 338:
 			if ( sscanf(optarg, "%i", &iargs.felix_opts.min_visits) != 1 )
 			{
 				ERROR("Invalid value for --felix-min-visits\n");
@@ -667,7 +683,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 39:
+			case 339:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.min_completeness) != 1 )
 			{
 				ERROR("Invalid value for --felix-min-completeness\n");
@@ -675,7 +691,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 40:
+			case 340:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.max_uniqueness) != 1 )
 			{
 				ERROR("Invalid value for --felix-max-uniqueness\n");
@@ -683,7 +699,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 41:
+			case 341:
 			if ( sscanf(optarg, "%i", &iargs.felix_opts.n_voxels) != 1 )
 			{
 				ERROR("Invalid value for --felix-num-voxels\n");
@@ -691,7 +707,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 42:
+			case 342:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.fraction_max_visits) != 1 )
 			{
 				ERROR("Invalid value for --felix-fraction-max-visits\n");
@@ -699,7 +715,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 43:
+			case 343:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.sigma) != 1 )
 			{
 				ERROR("Invalid value for --felix-sigma\n");
@@ -707,7 +723,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 44:
+			case 344:
 			if ( sscanf(optarg, "%i", &serial_start) != 1 )
 			{
 				ERROR("Invalid value for --serial-start\n");
@@ -715,7 +731,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 45:
+			case 345:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.domega) != 1 )
 			{
 				ERROR("Invalid value for --felix-domega\n");
@@ -723,7 +739,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 
-			case 46:
+			case 346:
 			if ( sscanf(optarg, "%lf", &iargs.felix_opts.max_internal_angle) != 1 )
 			{
 				ERROR("Invalid value for --felix-max-internal-angle\n");
@@ -731,6 +747,22 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 347:
+			iargs.min_snr_biggest_pix = strtof(optarg, NULL);
+			break;
+
+			case 348:
+			iargs.min_snr_peak_pix = strtof(optarg, NULL);
+			break;
+
+			case 349:
+			iargs.min_sig = strtof(optarg, NULL);
+			break;
+
+			case 350:
+			iargs.min_peak_over_neighbour = strtof(optarg, NULL);
+			break;
+
 			case 0 :
 			break;
 
@@ -789,6 +821,8 @@ int main(int argc, char *argv[])
 		iargs.peaks = PEAK_HDF5;
 	} else if ( strcmp(speaks, "cxi") == 0 ) {
 		iargs.peaks = PEAK_CXI;
+	} else if ( strcmp(speaks, "peakfinder9") == 0 ) {
+		iargs.peaks = PEAK_PEAKFINDER9;
 	} else {
 		ERROR("Unrecognised peak detection method '%s'\n", speaks);
 		return 1;
diff --git a/src/process_image.c b/src/process_image.c
index 4b02e694..1d41ae6e 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -3,7 +3,7 @@
  *
  * The processing pipeline for one image
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
@@ -201,7 +201,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		             iargs->use_saturated);
 		break;
 
-	        case PEAK_PEAKFINDER8:
+		case PEAK_PEAKFINDER8:
 		if ( search_peaks_peakfinder8(&image, 2048,
 		                               iargs->threshold,
 		                               iargs->min_snr,
@@ -223,6 +223,26 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		}
 		break;
 
+		case PEAK_PEAKFINDER9:
+		if ( search_peaks_peakfinder9(&image,
+		                              iargs->min_snr_biggest_pix,
+		                              iargs->min_snr_peak_pix,
+		                              iargs->min_snr,
+		                              iargs->min_sig,
+		                              iargs->min_peak_over_neighbour,
+		                              iargs->local_bg_radius) )
+		{
+			if ( image.event != NULL ) {
+				ERROR("Failed to find peaks in image %s"
+				      "(event %s).\n", image.filename,
+				      get_event_string(image.event));
+			} else {
+				ERROR("Failed to find peaks in image %s.",
+				      image.filename);
+			}
+		}
+		break;
+
 	}
 
 	restore_image_data(image.dp, image.det, prefilter);
diff --git a/src/process_image.h b/src/process_image.h
index 6d170a39..52f063fd 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -3,12 +3,13 @@
  *
  * The processing pipeline for one image
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2010-2016 Thomas White <taw@physics.org>
  *   2014-2017 Valerio Mariani <valerio.mariani@desy.de>
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -44,7 +45,8 @@ struct index_args;
 
 
 enum {
-        PEAK_PEAKFINDER8,
+	PEAK_PEAKFINDER9,
+	PEAK_PEAKFINDER8,
 	PEAK_ZAEF,
 	PEAK_HDF5,
 	PEAK_CXI,
@@ -83,6 +85,10 @@ struct index_args
 	int max_pix_count;
 	int local_bg_radius;
 	int min_peaks;
+	float min_snr_biggest_pix;
+	float min_snr_peak_pix;
+	float min_sig;
+	float min_peak_over_neighbour;
 	struct imagefile_field_list *copyme;
 	int integrate_saturated;
 	int use_saturated;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index df98b6a1..1fa8012a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -51,7 +51,11 @@ add_test(prof2d_check prof2d_check)
 
 add_executable(ring_check ring_check.c)
 target_include_directories(ring_check PRIVATE ${COMMON_INCLUDES})
-target_link_libraries(ring_check ${COMMON_LIBRARIES})
+target_link_libraries(ring_check PRIVATE ${COMMON_LIBRARIES})
+if (FDIP_FOUND)
+  target_include_directories(ring_check PRIVATE ${FDIP_INCLUDES})
+  target_link_libraries(ring_check PRIVATE ${FDIP_LIBRARIES})
+endif (FDIP_FOUND)
 add_test(ring_check ring_check)
 
 add_executable(symmetry_check symmetry_check.c)