aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/peakfinder8.c297
-rw-r--r--libcrystfel/src/peakfinder8.h14
-rw-r--r--libcrystfel/src/peaks.c12
-rw-r--r--libcrystfel/src/peaks.h3
-rw-r--r--src/gui_peaksearch.c2
-rw-r--r--src/indexamajig.c20
-rw-r--r--src/process_image.c4
-rw-r--r--src/process_image.h2
8 files changed, 320 insertions, 34 deletions
diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c
index 537c2365..3420a7a7 100644
--- a/libcrystfel/src/peakfinder8.c
+++ b/libcrystfel/src/peakfinder8.c
@@ -34,6 +34,7 @@
#include <float.h>
#include <math.h>
#include <stdlib.h>
+#include <profile.h>
#include "peakfinder8.h"
#include "detgeom.h"
@@ -43,13 +44,24 @@
/** \file peakfinder8.h */
// CrystFEL-only block 1
+
struct radius_maps
{
float **r_maps;
+ int *n_pixels;
int n_rmaps;
};
+struct radial_stats_pixels
+{
+ int n_panels;
+ int *n_pixels; // n_pixels[panel]
+ int **pidx; // pixel_index[panel][0..n_pixels]
+ int **radius; // pixel_radius[panel][0..n_pixels]
+};
+
+
struct peakfinder_mask
{
char **masks;
@@ -101,7 +113,138 @@ struct peakfinder_peak_data
};
-// CrystFEL-only block 2
+static struct radial_stats_pixels *compute_rstats_pixels(struct radius_maps *rmaps)
+{
+ int p;
+ int i;
+
+ struct radial_stats_pixels *rsp = NULL;
+ rsp = (struct radial_stats_pixels *)malloc(sizeof(struct radial_stats_pixels));
+ if ( rsp == NULL ) {
+ return NULL;
+ }
+ rsp->n_pixels = (int *)malloc(rmaps->n_rmaps * sizeof(int));
+ if ( rsp->n_pixels == NULL ) {
+ free(rsp);
+ return NULL;
+ }
+ rsp->pidx = (int **)malloc(rmaps->n_rmaps * sizeof(int *));
+ if ( rsp->pidx == NULL ) {
+ free(rsp->n_pixels);
+ free(rsp);
+ return NULL;
+ }
+ rsp->radius = (int **)malloc(rmaps->n_rmaps * sizeof(int *));
+ if ( rsp->radius == NULL ) {
+ free(rsp->n_pixels);
+ free(rsp->pidx);
+ free(rsp);
+ return NULL;
+ }
+ srand(0);
+
+ int n_pixels_per_bin = 100; // Can make this a parameter
+
+ // Assuming 5000 is the maximum possible radius
+ int n_pixels[5000] = {0}; // selected pixels per bin
+ int n_tot_pixels[5000] = {0}; // total pixels per bin
+ int panel[5000][n_pixels_per_bin]; // panel ID of selected pixels
+ int idx[5000][n_pixels_per_bin]; // index of selected pixels
+ int radius;
+
+ for ( p = 0; p < rmaps->n_rmaps; p++ ) {
+ rsp->n_pixels[p] = 0;
+ for ( i = 0; i < rmaps->n_pixels[p]; i++ ) {
+ // Reservoir sampling:
+ radius = (int)rint(rmaps->r_maps[p][i]);
+ n_tot_pixels[radius] += 1;
+
+ if ( n_pixels[radius] < n_pixels_per_bin ) {
+ panel[radius][n_pixels[radius]] = p;
+ idx[radius][n_pixels[radius]] = i;
+
+ n_pixels[radius] += 1;
+ rsp->n_pixels[p] += 1;
+ } else {
+ int rand_i = rand() % n_tot_pixels[radius];
+ if ( rand_i < n_pixels_per_bin ) {
+ rsp->n_pixels[panel[radius][rand_i]] -= 1;
+ rsp->n_pixels[p] += 1;
+
+ panel[radius][rand_i] = p;
+ idx[radius][rand_i] = i;
+ }
+ }
+ }
+ }
+
+ int *sidx = (int *)malloc(rmaps->n_rmaps * sizeof(int));
+ if ( sidx == NULL ) {
+ free(rsp->n_pixels);
+ free(rsp->pidx);
+ free(rsp->radius);
+ free(rsp);
+ return NULL;
+ }
+ for ( p = 0; p < rmaps->n_rmaps; p++ ) {
+ rsp->pidx[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int));
+ if ( rsp->pidx[p] == NULL ) {
+ for ( i = 0; i < p; i++ ) {
+ free(rsp->pidx[i]);
+ free(rsp->radius[i]);
+ }
+ free(rsp->pidx);
+ free(rsp->radius);
+ free(rsp->n_pixels);
+ free(rsp);
+ free(sidx);
+ return NULL;
+ }
+ rsp->radius[p] = (int *)malloc(rsp->n_pixels[p] * sizeof(int));
+ if ( rsp->radius[p] == NULL ) {
+ for ( i = 0; i < p; i++ ) {
+ free(rsp->pidx[i]);
+ free(rsp->radius[i]);
+ }
+ free(rsp->pidx[p]);
+ free(rsp->pidx);
+ free(rsp->radius);
+ free(rsp->n_pixels);
+ free(rsp);
+ free(sidx);
+ return NULL;
+ }
+ sidx[p] = 0;
+ }
+
+ for ( radius = 0; radius < 5000; radius++ ) {
+ for ( i = 0; i < n_pixels[radius]; i++ ) {
+ p = panel[radius][i];
+ rsp->pidx[p][sidx[p]] = idx[radius][i];
+ rsp->radius[p][sidx[p]] = radius;
+ sidx[p] += 1;
+ }
+ }
+ free(sidx);
+
+ rsp->n_panels = rmaps->n_rmaps;
+ return rsp;
+}
+
+static void free_rstats_pixels(struct radial_stats_pixels *rsp)
+{
+ int i;
+ for ( i = 0; i < rsp->n_panels; i++ ) {
+ free(rsp->pidx[i]);
+ free(rsp->radius[i]);
+ }
+ free(rsp->pidx);
+ free(rsp->radius);
+ free(rsp->n_pixels);
+ free(rsp);
+}
+
+
static struct radius_maps *compute_radius_maps(struct detgeom *det)
{
int i, u, iss, ifs;
@@ -119,6 +262,12 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det)
return NULL;
}
+ rm->n_pixels = (int *)malloc(det->n_panels*sizeof(int*));
+ if ( rm->r_maps == NULL ) {
+ free(rm);
+ return NULL;
+ }
+
rm->n_rmaps = det->n_panels;
for( i=0 ; i<det->n_panels ; i++ ) {
@@ -133,7 +282,7 @@ static struct radius_maps *compute_radius_maps(struct detgeom *det)
free(rm);
return NULL;
}
-
+ rm->n_pixels[i] = p.h * p.w;
for ( iss=0 ; iss<p.h ; iss++ ) {
for ( ifs=0; ifs<p.w; ifs++ ) {
@@ -161,10 +310,53 @@ static void free_radius_maps(struct radius_maps *r_maps)
free(r_maps->r_maps[i]);
}
free(r_maps->r_maps);
+ free(r_maps->n_pixels);
free(r_maps);
}
+// CrystFEL-only block 2
+struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode)
+{
+ struct pf8_private_data *data = NULL;
+ if ( det == NULL ) {
+ return NULL;
+ }
+
+ data = (struct pf8_private_data *)malloc(sizeof(struct pf8_private_data));
+ if ( data == NULL ) {
+ return NULL;
+ }
+ data->rmaps = compute_radius_maps(det);
+ if ( data->rmaps == NULL ) {
+ free(data);
+ return NULL;
+ }
+ if ( fast_mode ) {
+ data->rpixels = compute_rstats_pixels(data->rmaps);
+ if ( data->rpixels == NULL ) {
+ free_radius_maps(data->rmaps);
+ free(data);
+ return NULL;
+ }
+ } else {
+ data->rpixels = NULL;
+ }
+ data->fast_mode = fast_mode;
+ return data;
+}
+
+
+void free_pf8_private_data(struct pf8_private_data *data)
+{
+ free_radius_maps(data->rmaps);
+ if ( data->fast_mode ) {
+ free_rstats_pixels(data->rpixels);
+ }
+ free(data);
+}
+
+
static struct peakfinder_mask *create_peakfinder_mask(struct image *img,
struct radius_maps *rmps,
int min_res,
@@ -385,6 +577,31 @@ static void fill_radial_bins(float *data,
}
}
+static void fill_radial_bins_fast(float *data, int w, int h, int n_pixels,
+ int *pidx, int *radius, char *mask,
+ float *rthreshold, float *lthreshold,
+ float *roffset, float *rsigma, int *rcount)
+{
+ int i;
+
+ int curr_r;
+ float value;
+
+ for (i = 0; i < n_pixels; i++)
+ {
+ if (mask[pidx[i]] != 0)
+ {
+ curr_r = radius[i];
+ value = data[pidx[i]];
+ if (value < rthreshold[curr_r] && value > lthreshold[curr_r])
+ {
+ roffset[curr_r] += value;
+ rsigma[curr_r] += (value * value);
+ rcount[curr_r] += 1;
+ }
+ }
+ }
+}
static void compute_radial_stats(float *rthreshold,
float *lthreshold,
@@ -1021,9 +1238,13 @@ int peakfinder8(struct image *img, int max_n_peaks,
float threshold, float min_snr,
int min_pix_count, int max_pix_count,
int local_bg_radius, int min_res,
- int max_res, int use_saturated)
+ int max_res, int use_saturated,
+ int fast_mode, struct pf8_private_data *private_data)
{
+ struct pf8_private_data *geomdata;
struct radius_maps *rmaps;
+ struct radial_stats_pixels *rspixels;
+
struct peakfinder_mask *pfmask;
struct peakfinder_panel_data *pfdata;
struct radial_stats *rstats;
@@ -1040,18 +1261,28 @@ int peakfinder8(struct image *img, int max_n_peaks,
if ( img->detgeom == NULL) return 1;
- rmaps = compute_radius_maps(img->detgeom);
- if ( rmaps == NULL ) return 1;
+ profile_start("pf8-rmaps");
+ if ( private_data == NULL ) {
+ geomdata = prepare_peakfinder8(img->detgeom, fast_mode);
+ } else {
+ geomdata = private_data;
+ }
+ rmaps = geomdata->rmaps;
+ rspixels = geomdata->rpixels;
+ profile_end("pf8-rmaps");
+ if (geomdata == NULL) return 1;
+ profile_start("pf8-mask");
pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res);
+ profile_end("pf8-mask");
if ( pfmask == NULL ) {
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
return 1;
}
pfdata = allocate_panel_data(img->detgeom->n_panels);
if ( pfdata == NULL) {
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
free_peakfinder_mask(pfmask);
return 1;
}
@@ -1077,7 +1308,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
rstats = allocate_radial_stats(num_rad_bins);
if ( rstats == NULL ) {
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
free_peakfinder_mask(pfmask);
free_panel_data(pfdata);
return 1;
@@ -1087,7 +1318,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
rstats->rthreshold[i] = 1e9;
rstats->lthreshold[i] = -1e9;
}
-
+ profile_start("pf8-rstats");
for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) {
for ( i=0; i<num_rad_bins; i++ ) {
@@ -1097,18 +1328,31 @@ int peakfinder8(struct image *img, int max_n_peaks,
}
for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
-
- fill_radial_bins(pfdata->panel_data[pi],
- pfdata->panel_w[pi],
- pfdata->panel_h[pi],
- rmaps->r_maps[pi],
- pfmask->masks[pi],
- rstats->rthreshold,
- rstats->lthreshold,
- rstats->roffset,
- rstats->rsigma,
- rstats->rcount);
-
+ if ( fast_mode ) {
+ fill_radial_bins_fast(pfdata->panel_data[pi],
+ pfdata->panel_w[pi],
+ pfdata->panel_h[pi],
+ rspixels->n_pixels[pi],
+ rspixels->pidx[pi],
+ rspixels->radius[pi],
+ pfmask->masks[pi],
+ rstats->rthreshold,
+ rstats->lthreshold,
+ rstats->roffset,
+ rstats->rsigma,
+ rstats->rcount);
+ } else {
+ fill_radial_bins(pfdata->panel_data[pi],
+ pfdata->panel_w[pi],
+ pfdata->panel_h[pi],
+ rmaps->r_maps[pi],
+ pfmask->masks[pi],
+ rstats->rthreshold,
+ rstats->lthreshold,
+ rstats->roffset,
+ rstats->rsigma,
+ rstats->rcount);
+ }
}
compute_radial_stats(rstats->rthreshold,
@@ -1121,10 +1365,11 @@ int peakfinder8(struct image *img, int max_n_peaks,
threshold);
}
+ profile_end("pf8-rstats");
pkdata = allocate_peak_data(max_n_peaks);
if ( pkdata == NULL ) {
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
free_peakfinder_mask(pfmask);
free_panel_data(pfdata);
free_radial_stats(rstats);
@@ -1132,7 +1377,7 @@ int peakfinder8(struct image *img, int max_n_peaks,
}
remaining_max_num_peaks = max_n_peaks;
-
+ profile_start("pf8-search");
for ( pi=0 ; pi<img->detgeom->n_panels ; pi++) {
int peaks_to_add;
@@ -1165,10 +1410,11 @@ int peakfinder8(struct image *img, int max_n_peaks,
NULL);
if ( ret != 0 ) {
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
free_peakfinder_mask(pfmask);
free_panel_data(pfdata);
free_radial_stats(rstats);
+ profile_end("pf8-search");
return 1;
}
@@ -1199,8 +1445,9 @@ int peakfinder8(struct image *img, int max_n_peaks,
NULL);
}
}
+ profile_end("pf8-search");
- free_radius_maps(rmaps);
+ if ( private_data == NULL ) free_pf8_private_data(geomdata);
free_peakfinder_mask(pfmask);
free_panel_data(pfdata);
free_radial_stats(rstats);
diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h
index 28212e49..e5e86038 100644
--- a/libcrystfel/src/peakfinder8.h
+++ b/libcrystfel/src/peakfinder8.h
@@ -43,11 +43,23 @@ extern "C" {
* The "peakfinder8" peak search algorithm, originally implemented in Cheetah
*/
+struct pf8_private_data
+{
+ int fast_mode;
+ struct radius_maps *rmaps;
+ struct radial_stats_pixels *rpixels;
+};
+
+struct pf8_private_data *prepare_peakfinder8(struct detgeom *det, int fast_mode);
+
+void free_pf8_private_data(struct pf8_private_data *data);
+
extern int peakfinder8(struct image *img, int max_n_peaks,
float threshold, float min_snr,
int mix_pix_count, int max_pix_count,
int local_bg_radius, int min_res,
- int max_res, int use_saturated);
+ int max_res, int use_saturated,
+ int fast_mode, struct pf8_private_data *private_data);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index d260b66d..57d2c92e 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -454,10 +454,11 @@ void search_peaks(struct image *image, float threshold, float min_sq_gradient,
* the actual \ref peakfinder8 function, found in \ref peakfinder8.h.
*/
int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
- float threshold, float min_snr,
- int min_pix_count, int max_pix_count,
- int local_bg_radius, int min_res,
- int max_res, int use_saturated)
+ float threshold, float min_snr,
+ int min_pix_count, int max_pix_count,
+ int local_bg_radius, int min_res,
+ int max_res, int use_saturated,
+ int fast_mode, void *private_data)
{
if ( image->features != NULL ) {
image_feature_list_free(image->features);
@@ -467,7 +468,8 @@ int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
return peakfinder8(image, max_n_peaks, threshold, min_snr,
min_pix_count, max_pix_count,
local_bg_radius, min_res,
- max_res, use_saturated);
+ max_res, use_saturated,
+ fast_mode, private_data);
}
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 2f100db2..d6d5c1f9 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -73,7 +73,8 @@ extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
float threshold, float min_snr,
int mix_pix_count, int max_pix_count,
int local_bg_radius, int min_res,
- int max_res, int use_saturated);
+ int max_res, int use_saturated,
+ int fast_mode, void *private_data);
extern int search_peaks_peakfinder9(struct image *image,
float min_snr_biggest_pix,
diff --git a/src/gui_peaksearch.c b/src/gui_peaksearch.c
index 6af0b3a6..976e08ec 100644
--- a/src/gui_peaksearch.c
+++ b/src/gui_peaksearch.c
@@ -80,7 +80,7 @@ void update_peaks(struct crystfelproject *proj)
proj->peak_search_params.local_bg_radius,
proj->peak_search_params.min_res,
proj->peak_search_params.max_res,
- 1);
+ 1, 0, NULL);
break;
case PEAK_PEAKFINDER9:
diff --git a/src/indexamajig.c b/src/indexamajig.c
index add455b4..95420438 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -61,6 +61,8 @@
#include <integration.h>
#include <image.h>
#include <datatemplate.h>
+#include <detgeom.h>
+#include <peakfinder8.h>
#include "im-sandbox.h"
#include "im-zmq.h"
@@ -210,6 +212,7 @@ static void write_harvest_file(struct index_args *args,
write_float(fh, 1, "min_squared_gradient_adu2", args->min_sq_gradient);
write_float(fh, 1, "min_snr", args->min_snr);
write_bool(fh, 1, "check_hdf5_snr", args->check_hdf5_snr);
+ write_bool(fh, 1, "peakfinder8_fast", args->peakfinder8_fast);
write_bool(fh, 1, "half_pixel_shift", args->half_pixel_shift);
write_int(fh, 1, "min_res_px", args->min_res);
write_int(fh, 1, "max_res_px", args->max_res);
@@ -591,6 +594,10 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
args->iargs.check_hdf5_snr = 1;
break;
+ case 322:
+ args->iargs.peakfinder8_fast = 1;
+ break;
+
/* ---------- Indexing ---------- */
case 400 :
@@ -902,6 +909,8 @@ int main(int argc, char *argv[])
args.iargs.min_sig = 11.0;
args.iargs.min_peak_over_neighbour = -INFINITY;
args.iargs.check_hdf5_snr = 0;
+ args.iargs.peakfinder8_fast = 0;
+ args.iargs.pf_private = NULL;
args.iargs.dtempl = NULL;
args.iargs.peaks = PEAK_ZAEF;
args.iargs.half_pixel_shift = 1;
@@ -1020,6 +1029,7 @@ int main(int argc, char *argv[])
"locations by 0.5 pixels"},
{"check-hdf5-snr", 321, NULL, OPTION_NO_USAGE, "Check SNR for peaks from HDF5, "
"CXI or MsgPack (see --min-snr)"},
+ {"peakfinder8-fast", 322, NULL, OPTION_NO_USAGE, "peakfinder8 fast execution"},
{NULL, 0, 0, OPTION_DOC, "Indexing options:", 4},
{"indexing", 400, "method", 0, "List of indexing methods"},
@@ -1349,11 +1359,21 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
+ struct pf8_private_data *pf8_data = NULL;
+ struct detgeom *detgeom = NULL;
+ if ( args.iargs.peaks == PEAK_PEAKFINDER8 ) {
+ detgeom = data_template_get_2d_detgeom_if_possible(args.iargs.dtempl);
+ pf8_data = prepare_peakfinder8(detgeom, args.iargs.peakfinder8_fast);
+ args.iargs.pf_private = pf8_data;
+ }
+
r = create_sandbox(&args.iargs, args.n_proc, args.prefix, args.basename,
fh, st, tmpdir, args.serial_start,
&args.zmq_params, &args.asapo_params,
timeout, args.profile);
+ if ( pf8_data != NULL ) free_pf8_private_data(pf8_data);
+ if ( detgeom != NULL) detgeom_free(detgeom);
cell_free(args.iargs.cell);
free(args.prefix);
free(args.temp_location);
diff --git a/src/process_image.c b/src/process_image.c
index de2d8792..f1ddb7f0 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -322,7 +322,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
iargs->local_bg_radius,
iargs->min_res,
iargs->max_res,
- iargs->use_saturated) ) {
+ iargs->use_saturated,
+ iargs->peakfinder8_fast,
+ iargs->pf_private) ) {
ERROR("Failed to find peaks in image %s"
"(event %s).\n",
image->filename, image->ev);
diff --git a/src/process_image.h b/src/process_image.h
index e2f792a5..166f6cd5 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -81,6 +81,8 @@ struct index_args
float min_snr_peak_pix;
float min_sig;
float min_peak_over_neighbour;
+ int peakfinder8_fast;
+ void *pf_private;
/* Hit finding */
int min_peaks;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 13:21:59 +0000