diff options
author | Valerio Mariani <valerio.mariani@desy.de> | 2017-03-22 15:13:16 +0100 |
---|---|---|
committer | Valerio Mariani <valerio.mariani@desy.de> | 2017-03-22 15:13:16 +0100 |
commit | 6aea5300ee1c70ac2e99a7f3e6de4d26fb7e243a (patch) | |
tree | a1e43cd63f19cdc812690e90a9431bed2cbc6084 /libcrystfel/src | |
parent | 984379c6c48a2325134b5f1edfff78934604e1ae (diff) |
Completely revamped implementation of peakfinder8
Diffstat (limited to 'libcrystfel/src')
-rw-r--r-- | libcrystfel/src/peakfinder8.c | 813 |
1 files changed, 484 insertions, 329 deletions
diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 9ad0e4bd..0e916a0c 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -1,7 +1,7 @@ /* * peakfinder8.c * - * The processing pipeline for one image + * The peakfinder8 algorithm * * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. @@ -50,6 +50,8 @@ struct radial_stats { float *roffset; float *rthreshold; + float *rsigma; + int *rcount; int n_rad_bins; }; @@ -63,12 +65,22 @@ struct peakfinder_panel_data }; +struct peakfinder_intern_data +{ + char *pix_in_peak_map; + int *infs; + int *inss; + int *peak_pixels; +}; + + struct peakfinder_peak_data { int num_found_peaks; int *npix; float *com_fs; float *com_ss; + int *com_index; float *tot_i; float *max_i; float *sigma; @@ -80,16 +92,16 @@ static double compute_r_sigma(float *rsigma, int *rcount, float* roffset, int i) { return sqrt(rsigma[i] / rcount[i] - - ((roffset[i] / rcount[i]) * - (roffset[i] / rcount[i]))); + ((roffset[i] / rcount[i]) * + (roffset[i] / rcount[i]))); } static void update_radial_stats (float *roffset, float *rsigma, int *rcount, - float data, int curr_radius) + float value, int curr_radius) { - roffset[curr_radius] += data; - rsigma[curr_radius] += (data * data); + roffset[curr_radius] += value; + rsigma[curr_radius] += (value * value); rcount[curr_radius] += 1; } @@ -112,12 +124,12 @@ static struct radius_maps* compute_radius_maps(struct detector *det) struct panel p; struct radius_maps *rm = NULL; - rm = malloc(sizeof(struct radius_maps)); + rm = (struct radius_maps *)malloc(sizeof(struct radius_maps)); if ( rm == NULL ) { return NULL; } - rm->r_maps = malloc(det->n_panels*sizeof(float*)); + rm->r_maps = (float **)malloc(det->n_panels*sizeof(float*)); if ( rm->r_maps == NULL ) { free(rm); return NULL; @@ -128,7 +140,7 @@ static struct radius_maps* compute_radius_maps(struct detector *det) for( i=0 ; i<det->n_panels ; i++ ) { p = det->panels[i]; - rm->r_maps[i] = malloc(p.h*p.w*sizeof(float)); + rm->r_maps[i] = (float *)malloc(p.h*p.w*sizeof(float)); if ( rm->r_maps[i] == NULL ) { for ( u = 0; u<i; u++ ) { @@ -175,8 +187,8 @@ static struct peakfinder_mask *create_peakfinder_mask(struct image *img, int i; struct peakfinder_mask *msk; - msk = malloc(sizeof(struct peakfinder_mask)); - msk->masks = malloc(img->det->n_panels*sizeof(char*)); + msk = (struct peakfinder_mask *)malloc(sizeof(struct peakfinder_mask)); + msk->masks =(char **) malloc(img->det->n_panels*sizeof(char*)); msk->n_masks = img->det->n_panels; for ( i=0; i<img->det->n_panels; i++) { @@ -185,7 +197,7 @@ static struct peakfinder_mask *create_peakfinder_mask(struct image *img, p = img->det->panels[i]; - msk->masks[i] = calloc(p.w*p.h,sizeof(char)); + msk->masks[i] = (char *)calloc(p.w*p.h,sizeof(char)); for ( iss=0 ; iss<p.h ; iss++ ) { for ( ifs=0 ; ifs<p.w ; ifs++ ) { @@ -203,7 +215,6 @@ static struct peakfinder_mask *create_peakfinder_mask(struct image *img, msk->masks[i][idx] = 1; } - } } @@ -225,22 +236,27 @@ static void free_peakfinder_mask(struct peakfinder_mask * pfmask) } -static int compute_num_radial_bins( int num_panels, int *w, int *h, - float **r_maps ) +static int compute_num_radial_bins(int num_panels, int *w, int *h, + float **r_maps ) { - int m; float max_r; int max_r_int; + int m; max_r = -1e9; + for ( m=0 ; m<num_panels ; m++ ) { - int i; + int ifs, iss; + int pidx; - for ( i=0 ; i<w[m]*h[m] ; i++ ) { - if ( r_maps[m][i] > max_r ) { - max_r = r_maps[m][i]; + for ( iss=0 ; iss<h[m] ; iss++ ) { + for ( ifs=0 ; ifs<w[m] ; ifs++ ) { + pidx = iss * w[m] + ifs; + if ( r_maps[m][pidx] > max_r ) { + max_r = r_maps[m][pidx]; + } } } } @@ -251,39 +267,40 @@ static int compute_num_radial_bins( int num_panels, int *w, int *h, } -static struct peakfinder_panel_data* allocate_panel_data(struct image *img) +static struct peakfinder_panel_data* allocate_panel_data(int num_panels) { struct peakfinder_panel_data *pfdata; - pfdata = malloc(sizeof(struct peakfinder_panel_data)); + pfdata = (struct peakfinder_panel_data *)malloc( + sizeof(struct peakfinder_panel_data)); if ( pfdata == NULL ) { return NULL; } - pfdata->panel_h = malloc(img->det->n_panels*sizeof(int)); + pfdata->panel_h = (int *)malloc(num_panels*sizeof(int)); if ( pfdata->panel_h == NULL ) { free(pfdata); return NULL; } - pfdata->panel_w = malloc(img->det->n_panels*sizeof(int)); + pfdata->panel_w = (int *)malloc(num_panels*sizeof(int)); if ( pfdata->panel_w == NULL ) { + free(pfdata->panel_h); free(pfdata); - free(pfdata->panel_w); return NULL; } - pfdata->panel_data = malloc(img->det->n_panels*sizeof(float*)); + pfdata->panel_data = (float **)malloc(num_panels*sizeof(float*)); if ( pfdata->panel_data == NULL ) { - free(pfdata); free(pfdata->panel_w); free(pfdata->panel_h); + free(pfdata); return NULL; } - pfdata->num_panels = img->det->n_panels; + pfdata->num_panels = num_panels; return pfdata; } @@ -302,16 +319,38 @@ static struct radial_stats* allocate_radial_stats(int num_rad_bins) { struct radial_stats* rstats; - rstats = malloc(sizeof(struct radial_stats)); + rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats)); + if ( rstats == NULL ) { + return NULL; + } - rstats->roffset = malloc(num_rad_bins*sizeof(float)); + rstats->roffset = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->roffset == NULL ) { + free(rstats); return NULL; } - rstats->rthreshold = malloc(num_rad_bins*sizeof(float)); + rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float)); if ( rstats->rthreshold == NULL ) { + free(rstats); + free(rstats->roffset); + return NULL; + } + + rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float)); + if ( rstats->rsigma == NULL ) { + free(rstats); free(rstats->roffset); + free(rstats->rthreshold); + return NULL; + } + + rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int)); + if ( rstats->rcount == NULL ) { + free(rstats); + free(rstats->roffset); + free(rstats->rthreshold); + free(rstats->rsigma); return NULL; } @@ -325,110 +364,92 @@ static void free_radial_stats(struct radial_stats *rstats) { free(rstats->roffset); free(rstats->rthreshold); + free(rstats->rsigma); + free(rstats->rcount); free(rstats); } -static int compute_radial_stats(float ** panel_data, - int num_panels, - int *w, - int *h, - float **r_maps, - char **masks, - float *rthreshold, - float *roffset, - int num_rad_bins, - float min_snr, - float acd_threshold, - int iterations) +static void fill_radial_bins(float *data, + int w, + int h, + float *r_map, + char *mask, + float *rthreshold, + float *roffset, + float *rsigma, + int *rcount) { + int iss, ifs; + int pidx; - int i, ri, pi; - int it_counter; int curr_r; - float data; - float this_offset, this_sigma; - float *rsigma; - int *rcount; + float value; + for ( iss = 0; iss<h ; iss++ ) { - for ( i=0 ; i<num_rad_bins ; i++) { - rthreshold[i] = 1e9; - } + for ( ifs = 0; ifs<w ; ifs++ ) { + pidx = iss * w + ifs; - rsigma = malloc(num_rad_bins*sizeof(float)); - if ( rsigma == NULL ) { - return 1; - } + if ( mask[pidx] != 0 ) { - rcount = malloc(num_rad_bins*sizeof(int)); - if ( rcount == NULL ) { - free(rsigma); - return 1; - } + curr_r = (int)rint(r_map[pidx]); - for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) { + value = data[pidx]; - for ( i=0; i<num_rad_bins; i++ ) { - roffset[i] = 0; - rsigma[i] = 0; - rcount[i] = 0; + if ( value < rthreshold[curr_r ] ) { + + update_radial_stats(roffset, + rsigma, + rcount, + value, + curr_r); + } + } } + } +} - for ( pi=0 ; pi<num_panels ; pi++ ) { - int i; +static void compute_radial_stats(float *rthreshold, + float *roffset, + float *rsigma, + int *rcount, + int num_rad_bins, + float min_snr, + float acd_threshold) +{ - for ( i = 0; i < w[pi]*h[pi] ; i++ ) { - if ( masks[pi][i] != 0 ) { + int ri; + float this_offset, this_sigma; - curr_r = (int)rint(r_maps[pi][i]); + for ( ri=0 ; ri<num_rad_bins ; ri++ ) { - data = panel_data[pi][i]; + if ( rcount[ri] == 0 ) { + roffset[ri] = 0; + rsigma[ri] = 0; + rthreshold[ri] = 1e9; + } else { + this_offset = roffset[ri]/rcount[ri]; - if ( data < rthreshold[curr_r ] ) { + this_sigma = compute_r_sigma(rsigma, + rcount, + roffset, + ri); - update_radial_stats(roffset, - rsigma, - rcount, - data, - curr_r); - } - } - } - } - - for ( ri=0 ; ri<num_rad_bins ; ri++ ) { - if ( rcount[ri] == 0 ) { - roffset[ri] = 0; - rsigma[ri] = 0; - rthreshold[ri] = 1e9; - } else { - this_offset = roffset[ri]/rcount[ri]; - - this_sigma = compute_r_sigma(rsigma, - rcount, - roffset, - ri); - - roffset[ri] = this_offset; - rsigma[ri] = this_sigma; + roffset[ri] = this_offset; + rsigma[ri] = this_sigma; - rthreshold[ri] = roffset[ri] + - min_snr*rsigma[ri]; + rthreshold[ri] = roffset[ri] + + min_snr*rsigma[ri]; - if ( rthreshold[ri] < acd_threshold ) { - rthreshold[ri] = acd_threshold; - } + if ( rthreshold[ri] < acd_threshold ) { + rthreshold[ri] = acd_threshold; } } } - free(rsigma); - free(rcount); - - return 0; } @@ -436,26 +457,27 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) { struct peakfinder_peak_data *pkdata; - pkdata = malloc(sizeof(struct peakfinder_peak_data)); + pkdata = (struct peakfinder_peak_data*)malloc( + sizeof(struct peakfinder_peak_data)); if ( pkdata == NULL ) { return NULL; } - pkdata->npix = malloc(max_num_peaks*sizeof(int)); + pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int)); if ( pkdata->npix == NULL ) { free(pkdata); free(pkdata->npix); return NULL; } - pkdata->com_fs = malloc(max_num_peaks*sizeof(float)); + pkdata->com_fs = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->com_fs == NULL ) { free(pkdata->npix); free(pkdata); return NULL; } - pkdata->com_ss = malloc(max_num_peaks*sizeof(float)); + pkdata->com_ss = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->com_ss == NULL ) { free(pkdata->npix); free(pkdata->com_fs); @@ -463,41 +485,55 @@ struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks) return NULL; } - pkdata->tot_i = malloc(max_num_peaks*sizeof(float)); + pkdata->com_index = (int *)malloc(max_num_peaks*sizeof(int)); + if ( pkdata->com_ss == NULL ) { + free(pkdata->npix); + free(pkdata->com_fs); + free(pkdata->com_ss); + free(pkdata); + return NULL; + } + + + pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->tot_i == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata); return NULL; } - pkdata->max_i = malloc(max_num_peaks*sizeof(float)); + pkdata->max_i = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->max_i == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata); return NULL; } - pkdata->sigma = malloc(max_num_peaks*sizeof(float)); + pkdata->sigma = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->sigma == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata); return NULL; } - pkdata->snr = malloc(max_num_peaks*sizeof(float)); + pkdata->snr = (float *)malloc(max_num_peaks*sizeof(float)); if ( pkdata->snr == NULL ) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata->sigma); @@ -513,6 +549,7 @@ static void free_peak_data(struct peakfinder_peak_data *pkdata) { free(pkdata->npix); free(pkdata->com_fs); free(pkdata->com_ss); + free(pkdata->com_index); free(pkdata->tot_i); free(pkdata->max_i); free(pkdata->sigma); @@ -521,11 +558,70 @@ static void free_peak_data(struct peakfinder_peak_data *pkdata) { } -static void peak_search(int p, int *inss, int *infs, int h, int w, - int multip_w, float *copy, char *mask, float *r_map, +static struct peakfinder_intern_data *allocate_peakfinder_intern_data( + int data_size, int max_pix_count) +{ + + struct peakfinder_intern_data *intern_data; + + intern_data = (struct peakfinder_intern_data *)malloc( + sizeof(struct peakfinder_intern_data)); + if ( intern_data == NULL ) { + return NULL; + } + + intern_data->pix_in_peak_map =(char *)calloc(data_size, sizeof(char)); + if ( intern_data->pix_in_peak_map == NULL ) { + free(intern_data); + return NULL; + } + + intern_data->infs =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->infs == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data); + return NULL; + } + + intern_data->inss =(int *)calloc(data_size, sizeof(int)); + if ( intern_data->inss == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data); + return NULL; + } + + intern_data->peak_pixels =(int *)calloc(max_pix_count, sizeof(int)); + if ( intern_data->peak_pixels == NULL ) { + free(intern_data->pix_in_peak_map); + free(intern_data->infs); + free(intern_data->inss); + free(intern_data); + return NULL; + } + + return intern_data; +} + + +static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid) +{ + free(pfid->peak_pixels); + free(pfid->pix_in_peak_map); + free(pfid->infs); + free(pfid->inss); + free(pfid); +} + + + +static void peak_search(int p, + struct peakfinder_intern_data *pfinter, + float *copy, char *mask, float *r_map, float *rthreshold, float *roffset, - char *pix_in_peak_map, int *peak_pixels, - int *num_pix_in_peak, float *sum_com_fs, + int *num_pix_in_peak, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs, float *sum_com_fs, float *sum_com_ss, float *sum_i, int max_pix_count) { @@ -542,24 +638,24 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, for ( k=0; k<search_n; k++ ) { - if ( (infs[p] + search_fs[k]) < 0 ) + if ( (pfinter->infs[p] + search_fs[k]) < 0 ) continue; - if ( (infs[p] + search_fs[k]) >= w ) + if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue; - if ( (inss[p] + search_ss[k]) < 0 ) + if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue; - if ( (inss[p] + search_ss[k]) >= h ) + if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue; - curr_fs = infs[p] + search_fs[k]; - curr_ss = inss[p] + search_ss[k]; - pi = curr_fs + curr_ss * multip_w; + curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs; + curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; curr_radius = (int)rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; if ( copy[pi] > curr_threshold - && pix_in_peak_map[pi] == 0 + && pfinter->pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) { @@ -569,12 +665,15 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, *sum_com_fs += curr_i * ((float)curr_fs); *sum_com_ss += curr_i * ((float)curr_ss); - inss[*num_pix_in_peak] = inss[p] + search_ss[k]; - infs[*num_pix_in_peak] = infs[p] + search_fs[k]; + pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + + search_ss[k]; + pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + + search_fs[k]; - pix_in_peak_map[pi] = 1; - if ( *num_pix_in_peak < max_pix_count ) - peak_pixels[*num_pix_in_peak] = pi; + pfinter->pix_in_peak_map[pi] = 1; + if ( *num_pix_in_peak < max_pix_count ) { + pfinter->peak_pixels[*num_pix_in_peak] = pi; + } *num_pix_in_peak = *num_pix_in_peak+1; } } @@ -582,16 +681,17 @@ static void peak_search(int p, int *inss, int *infs, int h, int w, static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, - int h, int w, int multip_w, float *copy, - float *r_map, float *rthreshold, - float *roffset, char *pix_in_peak_map, - char *mask, float *local_sigma, + float *copy, float *r_map, + float *rthreshold, float *roffset, + char *pix_in_peak_map, char *mask, int asic_size_fs, + int asic_size_ss, int aifs, int aiss, + int num_pix_fs,float *local_sigma, float *local_offset, float *background_max_i, - int com_idx, - int local_bg_radius) + int com_idx, int local_bg_radius) { int ssj, fsi; float pix_radius; + int curr_fs, curr_ss; int pi; int curr_radius; float curr_threshold; @@ -617,26 +717,28 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, if ( (com_fs_int + fsi) < 0 ) continue; - if ( (com_fs_int + fsi) >= w ) + if ( (com_fs_int + fsi) >= asic_size_fs ) continue; if ( (com_ss_int + ssj) < 0 ) continue; - if ( (com_ss_int + ssj) >= h ) + if ( (com_ss_int + ssj) >= asic_size_ss ) continue; pix_radius = sqrt(fsi * fsi + ssj * ssj); if ( pix_radius>ring_width ) continue; - pi = (com_fs_int +fsi) + (com_ss_int +ssj) * multip_w; + curr_fs = com_fs_int +fsi + aifs * asic_size_fs; + curr_ss = com_ss_int +ssj + aiss * asic_size_ss; + pi = curr_fs + curr_ss * num_pix_fs; curr_radius = rint(r_map[pi]); curr_threshold = rthreshold[curr_radius]; curr_i = copy[pi]; if ( copy[pi] < curr_threshold - && pix_in_peak_map[pi] == 0 - && mask[pi] != 0 ) { + && pix_in_peak_map[pi] == 0 + && mask[pi] != 0 ) { np_sigma++; sum_i += curr_i; @@ -653,7 +755,7 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, if ( np_sigma != 0 ) { *local_offset = sum_i / np_sigma; *local_sigma = sqrt(sum_i_squared / np_sigma - - ((sum_i / np_sigma) * (sum_i / np_sigma))); + ((sum_i / np_sigma) * (sum_i / np_sigma))); } else { local_radius = rint(r_map[(int)rint(com_idx)]); @@ -663,121 +765,75 @@ static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int, } -static int peakfinder8_multipanel(float *roffset, float *rthreshold, - float *data, char *mask, float *r_map, int w, - int h, int max_n_peaks, int *num_found_peaks, - int *npix, float *com_fs, float *com_ss, - float *tot_i, float *max_i, float *sigma, - float *snr, int min_pix_count, - int max_pix_count, int local_bg_radius, - float min_snr, int multip_w) +static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs, + int aiss, int aifs, float *rthreshold, + float *roffset, int *peak_count, + float *copy, struct peakfinder_intern_data *pfinter, + float *r_map, char *mask, int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr, int max_n_peaks) { - int i, p, iss, ifs; - int peak_count; - int peak_idx; - int num_pix_in_peak, lt_num_pix_in_pk; - int idx, com_idx; - int curr_radius; - int ring_width; - int curr_idx; - int curr_fs, curr_ss; - int *peak_pixels; - int *infs; - int *inss; - float peak_com_fs, peak_com_ss; - float sum_com_fs, sum_com_ss, sum_i; - float peak_com_fs_int, peak_com_ss_int; - float curr_threshold; - float peak_tot_i, pk_tot_i_raw; - float peak_max_i, pk_max_i_raw; - float local_sigma; - float local_offset; - float background_max_i; - float f_background_thresh; - float peak_snr; - float curr_i, curr_i_raw; - float *copy; - char *pix_in_peak_map; - - copy = (float*) malloc(w*h*sizeof(float)); - if ( copy == NULL ) { - return 1; - } - - memcpy(copy, data, w*h*sizeof(float)); - pix_in_peak_map = calloc(w*h, sizeof(char)); - if ( pix_in_peak_map == NULL ) { - free(copy); - return 1; - } - - infs = (int *) calloc(w*h, sizeof(int)); - if ( infs == NULL ) { - free(pix_in_peak_map); - free(copy); - return 1; - } - - inss = (int *) calloc(w*h, sizeof(int)); - if ( inss == NULL ) { - free(infs); - free(pix_in_peak_map); - free(copy); - return 1; - } - - // TODO Possible bug: should 0,0 be included? - peak_pixels = calloc(max_pix_count, sizeof(int)); - if ( peak_pixels == NULL ) { - free(inss); - free(infs); - free(pix_in_peak_map); - free(copy); - return 1; - } - - for ( i = 0; i < w*h; i++ ) { - copy[i] *= mask[i]; - } - - peak_count = 0; - num_pix_in_peak = 0; - - for ( iss=0 ; iss<h ; iss++ ) { - for ( ifs=0 ; ifs<w ; ifs++ ) { - - idx = ifs+multip_w*iss; - - curr_radius = (int)rint(r_map[idx]); - curr_threshold = rthreshold[curr_radius]; - - if ( copy[idx] > curr_threshold && - pix_in_peak_map[idx] == 0 ) { - - infs[0] = ifs; - inss[0] = iss; - - peak_pixels[0] = idx; + int pxss, pxfs; + int num_pix_in_peak; + + for ( pxss=1 ; pxss<asic_size_ss-1 ; pxss++ ) { + for ( pxfs=1 ; pxfs<asic_size_fs-1 ; pxfs++ ) { + + float curr_thresh; + int pxidx; + int curr_rad; + + pxidx = (pxss + aiss * asic_size_ss) * + num_pix_fs + pxfs + + aifs * asic_size_fs; + + curr_rad = (int)rint(r_map[pxidx]); + curr_thresh = rthreshold[curr_rad]; + + if ( copy[pxidx] > curr_thresh + && pfinter->pix_in_peak_map[pxidx] == 0 ) { + + float sum_com_fs, sum_com_ss; + float sum_i; + float peak_com_fs, peak_com_ss; + float peak_com_fs_int, peak_com_ss_int; + float peak_tot_i, pk_tot_i_raw; + float peak_max_i, pk_max_i_raw; + float peak_snr; + float local_sigma, local_offset; + float background_max_i; + float f_background_thresh; + int lt_num_pix_in_pk; + int ring_width; + int peak_idx; + int com_idx; + int p; + + pfinter->infs[0] = pxfs; + pfinter->inss[0] = pxss; + pfinter->peak_pixels[0] = pxidx; num_pix_in_peak = 1; + sum_i = 0; sum_com_fs = 0; sum_com_ss = 0; - do { - + do { lt_num_pix_in_pk = num_pix_in_peak; for ( p=0; p<num_pix_in_peak; p++ ) { - - peak_search(p, inss, infs, - h, w, multip_w, - copy, mask, + peak_search(p, + pfinter, copy, mask, r_map, rthreshold, roffset, - pix_in_peak_map, - peak_pixels, &num_pix_in_peak, + asic_size_fs, + asic_size_ss, + aifs, aiss, + num_pix_fs, &sum_com_fs, &sum_com_ss, &sum_i, @@ -795,20 +851,29 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, peak_com_ss = sum_com_ss / fabs(sum_i); com_idx = rint(peak_com_fs) + - rint(peak_com_ss) * multip_w; + rint(peak_com_ss) * num_pix_fs; - peak_com_fs_int = (int)rint(peak_com_fs); - peak_com_ss_int = (int)rint(peak_com_ss); + peak_com_fs_int = (int)rint(peak_com_fs) - + aifs * asic_size_fs; + peak_com_ss_int = (int)rint(peak_com_ss) - + aiss * asic_size_ss; local_sigma = 0.0; local_offset = 0.0; background_max_i = 0.0; + ring_width = 2 * local_bg_radius; + search_in_ring(ring_width, peak_com_fs_int, - peak_com_ss_int, h, w, multip_w, + peak_com_ss_int, copy, r_map, rthreshold, - roffset, pix_in_peak_map, - mask, &local_sigma, + roffset, + pfinter->pix_in_peak_map, + mask, asic_size_fs, + asic_size_ss, + aifs, aiss, + num_pix_fs, + &local_sigma, &local_offset, &background_max_i, com_idx, local_bg_radius); @@ -821,11 +886,22 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, sum_com_ss = 0; for ( peak_idx = 1 ; - peak_idx < num_pix_in_peak && - peak_idx <= max_pix_count ; + peak_idx < num_pix_in_peak + && peak_idx <= max_pix_count ; peak_idx++ ) { - curr_idx = peak_pixels[peak_idx]; + int curr_idx; + float curr_i; + float curr_i_raw; + int curr_fs, curr_ss; + + // BUG HERE, I THINK. PEAK_PIXELS + // HAS SIZE MAX_PIX_COUNT, BUT + // IN THE FOLLOWING LINE PEAK_IDX + // CAN HAVE VALUE EXACTLY MAX_PEAK_COUNT + // (SEE THE FOR LOOP ABOVE) + curr_idx = + pfinter->peak_pixels[peak_idx]; curr_i_raw = copy[curr_idx]; curr_i = curr_i_raw - local_offset; @@ -833,8 +909,8 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, peak_tot_i += curr_i; pk_tot_i_raw += curr_i_raw; - curr_fs = curr_idx % multip_w; - curr_ss = curr_idx / multip_w; + curr_fs = curr_idx % asic_size_fs; + curr_ss = curr_idx / asic_size_fs; sum_com_fs += curr_i * ((float)curr_fs); sum_com_ss += curr_i * ((float)curr_ss); @@ -844,10 +920,11 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, peak_max_i = curr_i; } + peak_com_fs = sum_com_fs / fabs(peak_tot_i); peak_com_ss = sum_com_ss / fabs(peak_tot_i); - peak_snr = peak_tot_i / local_sigma; + peak_snr = peak_tot_i / local_sigma; if (peak_snr < min_snr) { continue; @@ -856,70 +933,113 @@ static int peakfinder8_multipanel(float *roffset, float *rthreshold, f_background_thresh = 1; f_background_thresh *= background_max_i - local_offset; - if (peak_max_i < f_background_thresh) + if (peak_max_i < f_background_thresh) { continue; + } if ( num_pix_in_peak >= min_pix_count && num_pix_in_peak <= max_pix_count ) { - if ( peak_tot_i == 0 ) + int peak_com_idx; + + if ( peak_tot_i == 0 ) { continue; + } + + peak_com_idx = rint(peak_com_fs) + + rint(peak_com_ss) * + asic_size_fs; - if ( peak_count < max_n_peaks ) { + if ( *peak_count < max_n_peaks ) { int pidx; - pidx = peak_count; + pidx = *peak_count; npix[pidx] = num_pix_in_peak; com_fs[pidx] = peak_com_fs; com_ss[pidx] = peak_com_ss; + com_index[pidx] = peak_com_idx; tot_i[pidx] = peak_tot_i; max_i[pidx] = peak_max_i; sigma[pidx] = local_sigma; snr[pidx] = peak_snr; - peak_count++; - } - else { - peak_count++; + *peak_count = *peak_count + 1; + } else { + *peak_count = *peak_count + 1; } } } } } +} + + +static int peakfinder8_base(float *roffset, float *rthreshold, + float *data, char *mask, float *r_map, + int asic_size_fs, int num_asics_fs, + int asic_size_ss, int num_asics_ss, + int max_n_peaks, int *num_found_peaks, + int *npix, float *com_fs, + float *com_ss, int *com_index, float *tot_i, + float *max_i, float *sigma, float *snr, + int min_pix_count, int max_pix_count, + int local_bg_radius, float min_snr) +{ + int num_pix_fs, num_pix_ss, num_pix_tot; + int i, aifs, aiss; + int peak_count; + float *copy; + struct peakfinder_intern_data *pfinter; + + num_pix_fs = asic_size_fs * num_asics_fs; + num_pix_ss = asic_size_ss * num_asics_ss; + num_pix_tot = num_pix_fs * num_pix_ss; + + copy = (float *)calloc(num_pix_tot, sizeof(float)); + if ( copy == NULL ) { + return 1; + } + + memcpy(copy, data, num_pix_tot*sizeof(float)); + + for (i = 0; i < num_pix_tot; i++) { + copy[i] *= mask[i]; + } + + pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count); + if ( pfinter == NULL ) { + free(copy); + return 1; + } + + peak_count = 0; + + for ( aiss=0 ; aiss<num_asics_ss ; aiss++ ) { + for ( aifs=0 ; aifs<num_asics_fs ; aifs++ ) { + + process_panel(asic_size_fs, asic_size_ss, num_pix_fs, + aiss, aifs, rthreshold, roffset, + &peak_count, copy, pfinter, r_map, mask, + npix, com_fs, com_ss, com_index, tot_i, + max_i, sigma, snr, min_pix_count, + max_pix_count, local_bg_radius, min_snr, + max_n_peaks); + } + } *num_found_peaks = peak_count; - free(peak_pixels); - free(inss); - free(infs); - free(pix_in_peak_map); + free_peakfinder_intern_data(pfinter); free(copy); return 0; } -static int peakfinder8_singlepanel(float *roffset, float *rthreshold, - float *data, char *mask, float *r_map, int w, - int h, int max_n_peaks, int *num_found_peaks, - int *npix, float *com_fs, float *com_ss, - float *tot_i, float *max_i, float *sigma, - float *snr, int min_pix_count, - int max_pix_count, int local_bg_radius, - float min_snr) -{ - return peakfinder8_multipanel(roffset, rthreshold, data, mask, r_map, - w, h, max_n_peaks, num_found_peaks, npix, - com_fs, com_ss, tot_i, max_i, sigma, snr, - min_pix_count, max_pix_count, - local_bg_radius, min_snr, w); - -} - 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) + 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) { struct radius_maps *rmaps; struct peakfinder_mask *pfmask; @@ -928,8 +1048,9 @@ int peakfinder8(struct image *img, int max_n_peaks, struct peakfinder_peak_data *pkdata; int num_rad_bins; int pi; - int ret; + int i, it_counter; int num_found_peaks; + int remaining_max_num_peaks; int iterations; iterations = 5; @@ -949,7 +1070,7 @@ int peakfinder8(struct image *img, int max_n_peaks, return 1; } - pfdata = allocate_panel_data(img); + pfdata = allocate_panel_data(img->det->n_panels); if ( pfdata == NULL) { free_radius_maps(rmaps); free_peakfinder_mask(pfmask); @@ -960,6 +1081,7 @@ int peakfinder8(struct image *img, int max_n_peaks, pfdata->panel_h[pi] = img->det->panels[pi].h; pfdata->panel_w[pi] = img->det->panels[pi].w; pfdata->panel_data[pi] = img->dp[pi]; + pfdata->num_panels = img->det->n_panels; } num_rad_bins = compute_num_radial_bins(img->det->n_panels, @@ -975,18 +1097,40 @@ int peakfinder8(struct image *img, int max_n_peaks, return 1; } - ret = compute_radial_stats(pfdata->panel_data, pfdata->num_panels, - pfdata->panel_w, pfdata->panel_h, - rmaps->r_maps, - pfmask->masks, rstats->rthreshold, - rstats->roffset, num_rad_bins, - min_snr, threshold, iterations); - if ( ret != 0 ) { - free_radius_maps(rmaps); - free_peakfinder_mask(pfmask); - free_panel_data(pfdata); - free_radial_stats(rstats); - return 1; + for ( i=0 ; i<rstats->n_rad_bins ; i++) { + rstats->rthreshold[i] = 1e9; + } + + for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) { + + for ( i=0; i<num_rad_bins; i++ ) { + rstats->roffset[i] = 0; + rstats->rsigma[i] = 0; + rstats->rcount[i] = 0; + } + + 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->roffset, + rstats->rsigma, + rstats->rcount); + + } + + compute_radial_stats(rstats->rthreshold, + rstats->roffset, + rstats->rsigma, + rstats->rcount, + num_rad_bins, + min_snr, + threshold); + } pkdata = allocate_peak_data(max_n_peaks); @@ -998,10 +1142,13 @@ int peakfinder8(struct image *img, int max_n_peaks, return 1; } + remaining_max_num_peaks = max_n_peaks; + for ( pi=0 ; pi<img->det->n_panels ; pi++) { int peaks_to_add; int pki; + int ret; num_found_peaks = 0; @@ -1009,34 +1156,43 @@ int peakfinder8(struct image *img, int max_n_peaks, continue; } - ret = peakfinder8_singlepanel(rstats->roffset, - rstats->rthreshold, - pfdata->panel_data[pi], - pfmask->masks[pi], - rmaps->r_maps[pi], - pfdata->panel_w[pi], - pfdata->panel_h[pi], - max_n_peaks, - &num_found_peaks, - pkdata->npix, - pkdata->com_fs, - pkdata->com_ss, - pkdata->tot_i, - pkdata->max_i, - pkdata->sigma, - pkdata->snr, - min_pix_count, - max_pix_count, - local_bg_radius, - min_snr); + ret = peakfinder8_base(rstats->roffset, + rstats->rthreshold, + pfdata->panel_data[pi], + pfmask->masks[pi], + rmaps->r_maps[pi], + pfdata->panel_w[pi], 1, + pfdata->panel_h[pi], 1, + max_n_peaks, + &num_found_peaks, + pkdata->npix, + pkdata->com_fs, + pkdata->com_ss, + pkdata->com_index, + pkdata->tot_i, + pkdata->max_i, + pkdata->sigma, + pkdata->snr, + min_pix_count, + max_pix_count, + local_bg_radius, + min_snr); + + if ( ret != 0 ) { + free_radius_maps(rmaps); + free_peakfinder_mask(pfmask); + free_panel_data(pfdata); + free_radial_stats(rstats); + return 1; + } peaks_to_add = num_found_peaks; - if ( num_found_peaks > max_n_peaks ) { - peaks_to_add = max_n_peaks; + if ( num_found_peaks > remaining_max_num_peaks ) { + peaks_to_add = remaining_max_num_peaks; } - max_n_peaks -= peaks_to_add; + remaining_max_num_peaks -= peaks_to_add; for ( pki=0 ; pki<peaks_to_add ; pki++ ) { @@ -1059,7 +1215,6 @@ int peakfinder8(struct image *img, int max_n_peaks, img, pkdata->tot_i[pki], NULL); - } } |