diff options
| -rw-r--r-- | doc/man/indexamajig.1 | 30 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.c | 106 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.h | 6 | ||||
| -rw-r--r-- | src/indexamajig.c | 49 | ||||
| -rw-r--r-- | src/process_image.c | 13 | ||||
| -rw-r--r-- | src/process_image.h | 10 |
6 files changed, 210 insertions, 4 deletions
diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index c0fe8277..bdcb7aa3 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 user 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--window-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--sig-fac-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--sig-fac-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,33 @@ 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 for \fB--peaks=zaef\fR and \fB--peaks=peakfinder8\fR, \fB--min-snr=9\fR for \fB--peaks=peakfinder9\fR. + +.PD 0 +.IP \fB--window-radius=<n>\fR +.PD +(peakFinder9 only) defines the radius in which the local background is estimated. Should be > max radius of peak. Default is 3. + +.PD 0 +.IP \fB--sig-fac-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--sig-fac-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. Can be set to -INFINITY to turn of speedup and search with max precision! Default is -INFINITY. .PD 0 .IP \fB--min-pix-count=\fR\fIcnt\fR diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 28c79538..7f6e0ba7 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -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. * @@ -563,6 +564,111 @@ int search_peaks_peakfinder8(struct image *image, int max_n_peaks, } +//#define HAVE_FDIP +#ifdef HAVE_FDIP + +#include "fastDiffractionImageProcessing/adaptions/crystfel/peakFinder9.h" +#include "fastDiffractionImageProcessing/adaptions/crystfel/mask.h" +#include "fastDiffractionImageProcessing/peakList.h" + +#include "fastDiffractionImageProcessing/peakFinder9.h" //debug, only for eclipse +#include "fastDiffractionImageProcessing/detectorRawFormat.h" //debug, only for eclipse + +int search_peaks_peakfinder9(struct image *image, float sig_fac_biggest_pix, + float sig_fac_peak_pix, float sig_fac_whole_peak, float min_sig, + float min_peak_over_neighbour, int window_radius) +{ + 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; + + peakFinder9_accuracyConstants_t accuracy_consts; + accuracy_consts.sigmaFactorBiggestPixel = sig_fac_biggest_pix; + accuracy_consts.sigmaFactorPeakPixel = sig_fac_peak_pix; + accuracy_consts.sigmaFactorWholePeak = sig_fac_whole_peak; + accuracy_consts.minimumSigma = min_sig; + accuracy_consts.minimumPeakOversizeOverNeighbours = + min_peak_over_neighbour; + accuracy_consts.windowRadius = (unsigned int) window_radius; + + long NpeaksMax = 10000; //more peaks per panel should not appear + + float *data_copy = NULL, *data_copy_new; + peakList_t peakList; + if (allocatePeakList(&peakList, NpeaksMax)) { + return 1; + } + + for (int panel_number = 0; panel_number < image->det->n_panels; + panel_number++) { + + if (image->det->panels[panel_number].no_index) + continue; + + int w = image->det->panels[panel_number].w; + int h = image->det->panels[panel_number].h; + + detectorRawFormat_t det_size_one_panel; + 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 = (float*) 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 (int 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 sig_fac_biggest_pix, + float sig_fac_peak_pix, float sig_fac_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 0; +} +#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..e63287d4 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -9,6 +9,7 @@ * 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,11 @@ 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 sig_fac_biggest_pix, + float sig_fac_peak_pix, float sig_fac_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..20d9faa1 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -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. * @@ -48,6 +49,7 @@ #include <sys/types.h> #include <sys/stat.h> #include <fcntl.h> +#include <float.h> #include "utils.h" #include "hdf5-file.h" @@ -84,7 +86,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,7 +100,7 @@ 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,pekafinder8, pekafinder9 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" @@ -109,6 +111,15 @@ static void show_help(const char *s) " (peakfinder8 only) Default: 0\n" " --max-res=<n> Maximum resolution for peak search (in pixels)\n" " (peakfinder8 only) Default: 1200\n" +" --window-radius=<n> (peakFinder9 only) defines the radius in which the " +" local background is estimated\n" +" --sig-fac-biggest-pix=<n> (peakFinder9 only) min snr of the biggest pixel in " +" the peak\n" +" --sig-fac-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" @@ -212,6 +223,7 @@ static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copym int main(int argc, char *argv[]) { int c; + unsigned int tmp_enum; char *filename = NULL; char *outfile = NULL; FILE *fh; @@ -260,6 +272,12 @@ int main(int argc, char *argv[]) iargs.min_res = 0; iargs.max_res = 1200; iargs.local_bg_radius = 3; + iargs.sig_fac_biggest_pix = 7.0; /* peak finder 9 */ + iargs.sig_fac_peak_pix = 6.0; + iargs.sig_fac_whole_peak = 9.0; + iargs.min_sig = 11.0; + iargs.min_peak_over_neighbour = -INFINITY; + iargs.window_radius = 3; iargs.check_hdf5_snr = 0; iargs.det = NULL; iargs.peaks = PEAK_ZAEF; @@ -405,6 +423,11 @@ int main(int argc, char *argv[]) {"serial-start", 1, NULL, 44}, {"felix-domega", 1, NULL, 45}, {"felix-max-internal-angle", 1, NULL, 46}, + {"sig-fac-biggest-pix" ,1, NULL, 47}, + {"sig-fac-peak-pix" ,1, NULL, 48}, + {"min-sig" ,1, NULL, 49}, + {"min-peak-over-neighbour" ,1, NULL, 50}, + {"window-radius" ,1, NULL, 51}, {0, 0, NULL, 0} }; @@ -504,6 +527,7 @@ int main(int argc, char *argv[]) case 11 : iargs.min_snr = strtof(optarg, NULL); + iargs.sig_fac_whole_peak = iargs.min_snr; break; case 13 : @@ -731,6 +755,25 @@ int main(int argc, char *argv[]) } break; + case 47: + iargs.sig_fac_biggest_pix = strtof(optarg, NULL); + break; + + case 48: + iargs.sig_fac_peak_pix = strtof(optarg, NULL); + break; + + case 49: + iargs.min_sig = strtof(optarg, NULL); + break; + + case 50: + iargs.min_peak_over_neighbour = strtof(optarg, NULL); + break; + + case 51: + iargs.window_radius = atoi(optarg); + case 0 : break; @@ -789,6 +832,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..be4fa094 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -223,6 +223,19 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } break; + case PEAK_PEAKFINDER9: + if ( search_peaks_peakfinder9(&image, + iargs->sig_fac_biggest_pix, + iargs->sig_fac_peak_pix, + iargs->sig_fac_whole_peak, iargs->min_sig, + iargs->min_peak_over_neighbour, + iargs->window_radius) ) + { + ERROR("Couldn't get enough memory to perform " + "peakFinder9.\n"); + } + break; + } restore_image_data(image.dp, image.det, prefilter); diff --git a/src/process_image.h b/src/process_image.h index 6d170a39..2b5bd19e 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -9,6 +9,7 @@ * 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,12 @@ struct index_args int max_pix_count; int local_bg_radius; int min_peaks; + float sig_fac_biggest_pix; + float sig_fac_peak_pix; + float sig_fac_whole_peak; + float min_sig; + float min_peak_over_neighbour; + int window_radius; struct imagefile_field_list *copyme; int integrate_saturated; int use_saturated; |