indexamajig: Add --median-filter

2 files changed, 116 insertions, 7 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 */