indexamajig: Add --median-filter

5 files changed, 137 insertions, 15 deletions

diff --git a/libcrystfel/src/filters.c b/libcrystfel/src/filters.c
index 0b84eb14..041b9686 100644
--- a/libcrystfel/src/filters.c
+++ b/libcrystfel/src/filters.c
@@ -96,7 +96,7 @@ void filter_cm(struct image *image)
 }
 
 
-void filter_noise(struct image *image, float *old)
+void filter_noise(struct image *image)
 {
 	int x, y;
 
@@ -106,8 +106,6 @@ void filter_noise(struct image *image, float *old)
 		int dx, dy;
 		int val = image->data[x+image->width*y];
 
-		if ( old != NULL ) old[x+image->width*y] = val;
-
 		/* FIXME: This isn't really the right thing to do
 		 * at the edges. */
 		if ( (x==0) || (x==image->width-1)
@@ -140,3 +138,113 @@ void filters_fudge_gslcblas()
 {
         STATUS("%p\n", cblas_sgemm);
 }
+
+
+#define SWAP(a,b) { float t=(a);(a)=(b);(b)=t; }
+static float kth_smallest(float *a, int n, int k)
+{
+	long i, j, l, m;
+	float x;
+
+	l = 0;
+	m = n-1;
+
+	while ( l < m ) {
+		x=a[k];
+		i=l;
+		j=m;
+		do {
+			while (a[i]<x) i++;
+			while (x<a[j]) j--;
+			if ( i<=j ) {
+				SWAP(a[i],a[j]);
+				i++;
+				j--;
+			}
+		} while (i<=j);
+		if ( j<k ) l = i;
+		if ( k<i ) m = j;
+	}
+	return a[k];
+}
+#undef SWAP
+
+
+void filter_median(struct image *image, int size)
+{
+	int counter;
+	int nn;
+	float *buffer;
+	float *localBg;
+	int pn;
+	int i;
+
+	if ( size <= 0 ) return;
+
+	nn = 2*size+1;
+	nn = nn*nn;
+
+	buffer = calloc(nn, sizeof(float));
+	localBg = calloc(image->width*image->height, sizeof(float));
+	if ( (buffer == NULL) || (localBg == NULL) ) {
+		ERROR("Failed to allocate LB buffers.\n");
+		return;
+	}
+
+	/* Determine local background
+	 * (median over window width either side of current pixel) */
+	for ( pn=0; pn<image->det->n_panels; pn++ ) {
+
+		int fs, ss;
+		struct panel *p;
+		int p_w, p_h;
+
+		p = &image->det->panels[pn];
+		p_w = (p->max_fs-p->min_fs)+1;
+		p_h = (p->max_ss-p->min_ss)+1;
+
+		for ( fs=0; fs<p_w; fs++ ) {
+		for ( ss=0; ss<p_h; ss++ ) {
+
+			int ifs, iss;
+			int e;
+
+			counter = 0;
+			e = fs+p->min_fs;
+			e += (ss+p->min_ss)*image->width;
+
+			// Loop over median window
+			for ( ifs=-size; ifs<=size; ifs++ ) {
+			for ( iss=-size; iss<=size; iss++ ) {
+
+				int idx;
+
+				if ( (fs+ifs) < 0 ) continue;
+				if ( (fs+ifs) >= p_w ) continue;
+				if ( (ss+iss) < 0 ) continue;
+				if ( (ss+iss) >= p_h ) continue;
+
+				idx = fs+ifs+p->min_fs;
+				idx += (ss+iss+p->min_ss)*image->width;
+
+				buffer[counter++] = image->data[idx];
+
+			}
+			}
+
+			// Find median value
+			localBg[e] = kth_smallest(buffer, counter, counter/2);
+
+		}
+		}
+	}
+
+
+	/* Do the background subtraction */
+	for ( i=0; i<image->width*image->height; i++ ) {
+		image->data[i] -= localBg[i];
+	}
+
+	free(localBg);
+	free(buffer);
+}
diff --git a/libcrystfel/src/filters.h b/libcrystfel/src/filters.h
index bda480c9..c8aeb52c 100644
--- a/libcrystfel/src/filters.h
+++ b/libcrystfel/src/filters.h
@@ -3,11 +3,11 @@
  *
  * Image filtering
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010,2012 Thomas White <taw@physics.org>
+ *   2010,2012-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -35,7 +35,8 @@
 
 
 extern void filter_cm(struct image *image);
-extern void filter_noise(struct image *image, float *old);
+extern void filter_noise(struct image *image);
+extern void filter_median(struct image *image, int size);
 
 
 #endif	/* FILTERS_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 9b868b9f..27ec31fc 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -111,6 +111,10 @@ static void show_help(const char *s)
 "                         pixels in each 3x3 region to zero if any of them\n"
 "                         have negative values.  Intensity measurement will\n"
 "                         be performed on the image as it was before this.\n"
+"    --median-filter=<n> Apply a median filter to the image data.  Intensity\n"
+"                         measurement will be performed on the image as it\n"
+"                         was before this.  The side length of the median\n"
+"                         filter box will be <n>.  Default: 0 (no filter).\n"
 "    --no-sat-corr       Don't correct values of saturated peaks using a\n"
 "                         table included in the HDF5 file.\n"
 "    --threshold=<n>     Only accept peaks above <n> ADU.  Default: 800.\n"
@@ -183,6 +187,7 @@ int main(int argc, char *argv[])
 	iargs.cell = NULL;
 	iargs.cmfilter = 0;
 	iargs.noisefilter = 0;
+	iargs.median_filter = 0;
 	iargs.satcorr = 1;
 	iargs.closer = 0;
 	iargs.bgsub = 1;
@@ -266,6 +271,7 @@ int main(int argc, char *argv[])
 		{"min-integration-snr",1, NULL,               12},
 		{"tolerance",          1, NULL,               13},
 		{"int-radius",         1, NULL,               14},
+		{"median-filter",      1, NULL,               15},
 
 		{0, 0, NULL, 0}
 	};
@@ -386,6 +392,10 @@ int main(int argc, char *argv[])
 			intrad = strdup(optarg);
 			break;
 
+			case 15 :
+			iargs.median_filter = atoi(optarg);
+			break;
+
 			case 0 :
 			break;
 
diff --git a/src/process_image.c b/src/process_image.c
index 010054c8..a8eeb663 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -136,23 +136,25 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	if ( iargs->cmfilter ) filter_cm(&image);
 
-	/* Take snapshot of image after CM subtraction but before
-	 * the aggressive noise filter. */
+	/* Take snapshot of image after CM subtraction but before applying
+	 * horrible noise filters to it */
 	data_size = image.width * image.height * sizeof(float);
 	data_for_measurement = malloc(data_size);
+	memcpy(data_for_measurement, image.data, data_size);
+
+	if ( iargs->median_filter > 0 ) {
+		filter_median(&image, iargs->median_filter);
+	}
 
 	if ( iargs->noisefilter ) {
-		filter_noise(&image, data_for_measurement);
-	} else {
-		memcpy(data_for_measurement, image.data, data_size);
+		filter_noise(&image);
 	}
 
 	switch ( iargs->peaks ) {
 
 		case PEAK_HDF5:
-		// Get peaks from HDF5
-		if (get_peaks(&image, hdfile,
-			iargs->hdf5_peak_path)) {
+		/* Get peaks from HDF5 */
+		if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) {
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
 		if ( !iargs->no_revalidate ) {
diff --git a/src/process_image.h b/src/process_image.h
index f2f034d9..5fe11e33 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -46,6 +46,7 @@ struct index_args
 	UnitCell *cell;
 	int cmfilter;
 	int noisefilter;
+	int median_filter;
 	int satcorr;
 	int closer;
 	int bgsub;