Completely revamped implementation of peakfinder8

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);
-
 		}
 	}