[Scipy-svn] r5279 - in trunk/scipy/signal: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Dec 19 14:34:05 CST 2008


Author: jarrod.millman
Date: 2008-12-19 14:34:02 -0600 (Fri, 19 Dec 2008)
New Revision: 5279

Modified:
   trunk/scipy/signal/medianfilter.c
   trunk/scipy/signal/tests/test_signaltools.py
Log:
Applying patch from Ray Jones to medianfilter to remove code based
on Numerical Recipes.  It includes the following:
 * new quickselect (from scratch) in medianfilter.c
 * new macros for {f,d,b}_medfilt2 in medianfilter.c
 * modified test for medfilt and medfilt2 (larger array, non-square filter).



Modified: trunk/scipy/signal/medianfilter.c
===================================================================
--- trunk/scipy/signal/medianfilter.c	2008-12-19 17:55:32 UTC (rev 5278)
+++ trunk/scipy/signal/medianfilter.c	2008-12-19 20:34:02 UTC (rev 5279)
@@ -1,311 +1,127 @@
 
 /*--------------------------------------------------------------------*/
 
-/*
- *  This Quickselect routine is based on the algorithm described in
- *  "Numerical recipies in C", Second Edition,
- *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
- */
-
 #include "Python.h"
 #define NO_IMPORT_ARRAY
 #include "numpy/noprefix.h"
 
+
+/* defined below */
 void f_medfilt2(float*,float*,intp*,intp*);
 void d_medfilt2(double*,double*,intp*,intp*);
 void b_medfilt2(unsigned char*,unsigned char*,intp*,intp*);
 extern char *check_malloc (int);
 
-#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
 
-float f_quick_select(float arr[], int n) 
-{
-    int low, high ;
-    int median;
-    int middle, ll, hh;
+/* The QUICK_SELECT routine is based on Hoare's Quickselect algorithm,
+ * with unrolled recursion. 
+ * Author: Thouis R. Jones, 2008 
+ */
 
-    low = 0 ; high = n-1 ; median = (low + high) / 2;
-    for (;;) {
-        if (high <= low) /* One element only */
-            return arr[median] ;
+#define ELEM_SWAP(t, a, x, y) {register t temp = (a)[x]; (a)[x] = (a)[y]; (a)[y] = temp;}
+#define FIRST_LOWEST(x, y, z) (((x) < (y)) && ((x) < (z)))
+#define FIRST_HIGHEST(x, y, z) (((x) > (y)) && ((x) > (z)))
+#define LOWEST_IDX(a, x, y) (((a)[x] < (a)[y]) ? (x) : (y))
+#define HIGHEST_IDX(a, x, y) (((a)[x] > (a)[y]) ? (x) : (y))
 
-        if (high == low + 1) {  /* Two elements only */
-            if (arr[low] > arr[high])
-                ELEM_SWAP(arr[low], arr[high]) ;
-            return arr[median] ;
-        }
+/* if (l is index of lowest) {return lower of mid,hi} else if (l is index of highest) {return higher of mid,hi} else return l */
+#define MEDIAN_IDX(a, l, m, h) (FIRST_LOWEST((a)[l], (a)[m], (a)[h]) ? LOWEST_IDX(a, m, h) : (FIRST_HIGHEST((a)[l], (a)[m], (a)[h]) ? HIGHEST_IDX(a, m, h) : (l)))
 
-    /* Find median of low, middle and high items; swap into position low */
-    middle = (low + high) / 2;
-    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
-    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
-    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
-
-    /* Swap low item (now in position middle) into position (low+1) */
-    ELEM_SWAP(arr[middle], arr[low+1]) ;
-
-    /* Nibble from each end towards middle, swapping items when stuck */
-    ll = low + 1;
-    hh = high;
-    for (;;) {
-        do ll++; while (arr[low] > arr[ll]) ;
-        do hh--; while (arr[hh]  > arr[low]) ;
-
-        if (hh < ll)
-        break;
-
-        ELEM_SWAP(arr[ll], arr[hh]) ;
-    }
-
-    /* Swap middle item (in position low) back into correct position */
-    ELEM_SWAP(arr[low], arr[hh]) ;
-
-    /* Re-set active partition */
-    if (hh <= median)
-        low = ll;
-        if (hh >= median)
-        high = hh - 1;
-    }
+#define QUICK_SELECT(NAME, TYPE)                                        \
+TYPE NAME(TYPE arr[], int n)                                            \
+{                                                                       \
+    int lo, hi, mid, md;                                                \
+    int median_idx;                                                     \
+    int ll, hh;                                                         \
+    TYPE piv;                                                           \
+                                                                        \
+    lo = 0; hi = n-1;                                                   \
+    median_idx = (n - 1) / 2; /* lower of middle values for even-length arrays */ \
+                                                                        \
+    while (1) {                                                         \
+        if ((hi - lo) < 2) {                                            \
+            if (arr[hi] < arr[lo]) ELEM_SWAP(TYPE, arr, lo, hi);        \
+            return arr[median_idx];                                     \
+        }                                                               \
+                                                                        \
+        mid = (hi + lo) / 2;                                            \
+        /* put the median of lo,mid,hi at position lo - this will be the pivot */ \
+        md = MEDIAN_IDX(arr, lo, mid, hi);                              \
+        ELEM_SWAP(TYPE, arr, lo, md);                                   \
+                                                                        \
+        /* Nibble from each end towards middle, swapping misordered items */ \
+        piv = arr[lo];                                                  \
+        for (ll = lo+1, hh = hi;; ll++, hh--) {                         \
+	    while (arr[ll] < piv) ll++;					\
+	    while (arr[hh] > piv) hh--;					\
+	    if (hh < ll) break;						\
+	    ELEM_SWAP(TYPE, arr, ll, hh);				\
+        }                                                               \
+        /* move pivot to top of lower partition */                      \
+        ELEM_SWAP(TYPE, arr, hh, lo);                                   \
+        /* set lo, hi for new range to search */                        \
+        if (hh < median_idx) /* search upper partition */               \
+            lo = hh+1;                                                  \
+        else if (hh > median_idx) /* search lower partition */          \
+            hi = hh-1;                                                  \
+        else                                                            \
+            return piv;                                                 \
+    }                                                                   \
 }
 
-#undef ELEM_SWAP
 
-
-#define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
-
-double d_quick_select(double arr[], int n) 
-{
-    int low, high ;
-    int median;
-    int middle, ll, hh;
-
-    low = 0 ; high = n-1 ; median = (low + high) / 2;
-    for (;;) {
-        if (high <= low) /* One element only */
-            return arr[median] ;
-
-        if (high == low + 1) {  /* Two elements only */
-            if (arr[low] > arr[high])
-                ELEM_SWAP(arr[low], arr[high]) ;
-            return arr[median] ;
-        }
-
-    /* Find median of low, middle and high items; swap into position low */
-    middle = (low + high) / 2;
-    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
-    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
-    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
-
-    /* Swap low item (now in position middle) into position (low+1) */
-    ELEM_SWAP(arr[middle], arr[low+1]) ;
-
-    /* Nibble from each end towards middle, swapping items when stuck */
-    ll = low + 1;
-    hh = high;
-    for (;;) {
-        do ll++; while (arr[low] > arr[ll]) ;
-        do hh--; while (arr[hh]  > arr[low]) ;
-
-        if (hh < ll)
-        break;
-
-        ELEM_SWAP(arr[ll], arr[hh]) ;
-    }
-
-    /* Swap middle item (in position low) back into correct position */
-    ELEM_SWAP(arr[low], arr[hh]) ;
-
-    /* Re-set active partition */
-    if (hh <= median)
-        low = ll;
-        if (hh >= median)
-        high = hh - 1;
-    }
-}
-
-#undef ELEM_SWAP
-
-#define ELEM_SWAP(a,b) { register unsigned char t=(a);(a)=(b);(b)=t; }
-
-unsigned char b_quick_select(unsigned char arr[], int n) 
-{
-    int low, high ;
-    int median;
-    int middle, ll, hh;
-
-    low = 0 ; high = n-1 ; median = (low + high) / 2;
-    for (;;) {
-        if (high <= low) /* One element only */
-            return arr[median] ;
-
-        if (high == low + 1) {  /* Two elements only */
-            if (arr[low] > arr[high])
-                ELEM_SWAP(arr[low], arr[high]) ;
-            return arr[median] ;
-        }
-
-    /* Find median of low, middle and high items; swap into position low */
-    middle = (low + high) / 2;
-    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
-    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
-    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
-
-    /* Swap low item (now in position middle) into position (low+1) */
-    ELEM_SWAP(arr[middle], arr[low+1]) ;
-
-    /* Nibble from each end towards middle, swapping items when stuck */
-    ll = low + 1;
-    hh = high;
-    for (;;) {
-        do ll++; while (arr[low] > arr[ll]) ;
-        do hh--; while (arr[hh]  > arr[low]) ;
-
-        if (hh < ll)
-        break;
-
-        ELEM_SWAP(arr[ll], arr[hh]) ;
-    }
-
-    /* Swap middle item (in position low) back into correct position */
-    ELEM_SWAP(arr[low], arr[hh]) ;
-
-    /* Re-set active partition */
-    if (hh <= median)
-        low = ll;
-        if (hh >= median)
-        high = hh - 1;
-    }
-}
-
-#undef ELEM_SWAP
-
 /* 2-D median filter with zero-padding on edges. */
-void d_medfilt2(double* in, double* out, intp* Nwin, intp* Ns)
-{ 
-  int nx, ny, hN[2];
-  int pre_x, pre_y, pos_x, pos_y;
-  int subx, suby, k, totN;
-  double *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
-
-  totN = Nwin[0] * Nwin[1];
-  myvals = (double *) check_malloc( totN * sizeof(double));
-
-  hN[0] = Nwin[0] >> 1;
-  hN[1] = Nwin[1] >> 1;
-  ptr1 = in;
-  fptr1 = out;
-  for (ny = 0; ny < Ns[0]; ny++)
-    for (nx = 0; nx < Ns[1]; nx++) {
-      pre_x = hN[1];
-      pre_y = hN[0];
-      pos_x = hN[1];
-      pos_y = hN[0];
-      if (nx < hN[1]) pre_x = nx;
-      if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
-      if (ny < hN[0]) pre_y = ny;
-      if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
-      fptr2 = myvals;
-      ptr2 = ptr1 - pre_x - pre_y*Ns[1];
-      for (suby = -pre_y; suby <= pos_y; suby++) {
-	for (subx = -pre_x; subx <= pos_x; subx++) 	  
-	  *fptr2++ = *ptr2++;
-	ptr2 += Ns[1] - (pre_x + pos_x + 1);
-      }
-      ptr1++;
-
-      /* Zero pad */
-      for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
-	*fptr2++ = 0.0;
-
-      /*      *fptr1++ = median(myvals,totN); */
-      *fptr1++ = d_quick_select(myvals,totN);
-    }
+#define MEDIAN_FILTER_2D(NAME, TYPE, SELECT)                            \
+void NAME(TYPE* in, TYPE* out, intp* Nwin, intp* Ns)                    \
+{                                                                       \
+    int nx, ny, hN[2];                                                  \
+    int pre_x, pre_y, pos_x, pos_y;                                     \
+    int subx, suby, k, totN;                                            \
+    TYPE *myvals, *fptr1, *fptr2, *ptr1, *ptr2;                         \
+                                                                        \
+    totN = Nwin[0] * Nwin[1];                                           \
+    myvals = (TYPE *) check_malloc( totN * sizeof(TYPE));               \
+                                                                        \
+    hN[0] = Nwin[0] >> 1;                                               \
+    hN[1] = Nwin[1] >> 1;                                               \
+    ptr1 = in;                                                          \
+    fptr1 = out;                                                        \
+    for (ny = 0; ny < Ns[0]; ny++)                                      \
+        for (nx = 0; nx < Ns[1]; nx++) {                                \
+            pre_x = hN[1];                                              \
+            pre_y = hN[0];                                              \
+            pos_x = hN[1];                                              \
+            pos_y = hN[0];                                              \
+            if (nx < hN[1]) pre_x = nx;                                 \
+            if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;            \
+            if (ny < hN[0]) pre_y = ny;                                 \
+            if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;            \
+            fptr2 = myvals;                                             \
+            ptr2 = ptr1 - pre_x - pre_y*Ns[1];                          \
+            for (suby = -pre_y; suby <= pos_y; suby++) {                \
+                for (subx = -pre_x; subx <= pos_x; subx++)              \
+                    *fptr2++ = *ptr2++;                                 \
+                ptr2 += Ns[1] - (pre_x + pos_x + 1);                    \
+            }                                                           \
+            ptr1++;                                                     \
+                                                                        \
+            /* Zero pad */                                              \
+            for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++) \
+                *fptr2++ = 0.0;                                         \
+                                                                        \
+            /*      *fptr1++ = median(myvals,totN); */                  \
+            *fptr1++ = SELECT(myvals,totN);                             \
+        }                                                               \
+    free(myvals);                                                       \
 }
 
 
-/* 2-D median filter with zero-padding on edges. */
-void f_medfilt2(float* in, float* out, intp* Nwin, intp* Ns)
-{ 
-  int nx, ny, hN[2];
-  int pre_x, pre_y, pos_x, pos_y;
-  int subx, suby, k, totN;
-  float *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
+/* define quick_select for floats, doubles, and unsigned characters */
+QUICK_SELECT(f_quick_select, float)
+QUICK_SELECT(d_quick_select, double)
+QUICK_SELECT(b_quick_select, unsigned char)
 
-  totN = Nwin[0] * Nwin[1];
-  myvals = (float *) check_malloc( totN * sizeof(float));
-
-  hN[0] = Nwin[0] >> 1;
-  hN[1] = Nwin[1] >> 1;
-  ptr1 = in;
-  fptr1 = out;
-  for (ny = 0; ny < Ns[0]; ny++)
-    for (nx = 0; nx < Ns[1]; nx++) {
-      pre_x = hN[1];
-      pre_y = hN[0];
-      pos_x = hN[1];
-      pos_y = hN[0];
-      if (nx < hN[1]) pre_x = nx;
-      if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
-      if (ny < hN[0]) pre_y = ny;
-      if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
-      fptr2 = myvals;
-      ptr2 = ptr1 - pre_x - pre_y*Ns[1];
-      for (suby = -pre_y; suby <= pos_y; suby++) {
-	for (subx = -pre_x; subx <= pos_x; subx++) 	  
-	  *fptr2++ = *ptr2++;
-	ptr2 += Ns[1] - (pre_x + pos_x + 1);
-      }
-      ptr1++;
-
-      /* Zero pad */
-      for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
-	*fptr2++ = 0.0;
-
-      /*      *fptr1++ = median(myvals,totN); */
-      *fptr1++ = f_quick_select(myvals,totN);
-    }
-}
-
-
-/* 2-D median filter with zero-padding on edges. */
-void b_medfilt2(unsigned char *in, unsigned char *out, intp* Nwin, intp* Ns)
-{ 
-  int nx, ny, hN[2];
-  int pre_x, pre_y, pos_x, pos_y;
-  int subx, suby, k, totN;
-  unsigned char *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
-
-  totN = Nwin[0] * Nwin[1];
-  myvals = (unsigned char *) check_malloc( totN * sizeof(unsigned char));
-
-  hN[0] = Nwin[0] >> 1;
-  hN[1] = Nwin[1] >> 1;
-  ptr1 = in;
-  fptr1 = out;
-  for (ny = 0; ny < Ns[0]; ny++)
-    for (nx = 0; nx < Ns[1]; nx++) {
-      pre_x = hN[1];
-      pre_y = hN[0];
-      pos_x = hN[1];
-      pos_y = hN[0];
-      if (nx < hN[1]) pre_x = nx;
-      if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
-      if (ny < hN[0]) pre_y = ny;
-      if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
-      fptr2 = myvals;
-      ptr2 = ptr1 - pre_x - pre_y*Ns[1];
-      for (suby = -pre_y; suby <= pos_y; suby++) {
-	for (subx = -pre_x; subx <= pos_x; subx++) 	  
-	  *fptr2++ = *ptr2++;
-	ptr2 += Ns[1] - (pre_x + pos_x + 1);
-      }
-      ptr1++;
-
-      /* Zero pad */
-      for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
-	*fptr2++ = 0.0;
-
-      /*      *fptr1++ = median(myvals,totN); */
-      *fptr1++ = b_quick_select(myvals,totN);
-    }
-}
+/* define medfilt for floats, doubles, and unsigned characters */
+MEDIAN_FILTER_2D(f_medfilt2, float, f_quick_select)
+MEDIAN_FILTER_2D(d_medfilt2, double, d_quick_select)
+MEDIAN_FILTER_2D(b_medfilt2, unsigned char, b_quick_select)

Modified: trunk/scipy/signal/tests/test_signaltools.py
===================================================================
--- trunk/scipy/signal/tests/test_signaltools.py	2008-12-19 17:55:32 UTC (rev 5278)
+++ trunk/scipy/signal/tests/test_signaltools.py	2008-12-19 20:34:02 UTC (rev 5279)
@@ -28,9 +28,30 @@
 
 class TestMedFilt(TestCase):
     def test_basic(self):
-        f = [[3,4,5],[2,3,4],[1,2,5]]
-        d = signal.medfilt(f)
-        assert_array_equal(d, [[0,3,0],[2,3,3],[0,2,0]])
+        f = [[50, 50, 50, 50, 50, 92, 18, 27, 65, 46],
+             [50, 50, 50, 50, 50,  0, 72, 77, 68, 66],
+             [50, 50, 50, 50, 50, 46, 47, 19, 64, 77],
+             [50, 50, 50, 50, 50, 42, 15, 29, 95, 35],
+             [50, 50, 50, 50, 50, 46, 34,  9, 21, 66],
+             [70, 97, 28, 68, 78, 77, 61, 58, 71, 42],
+             [64, 53, 44, 29, 68, 32, 19, 68, 24, 84],
+             [ 3, 33, 53, 67,  1, 78, 74, 55, 12, 83],
+             [ 7, 11, 46, 70, 60, 47, 24, 43, 61, 26],
+             [32, 61, 88,  7, 39,  4, 92, 64, 45, 61]]
+        
+        d = signal.medfilt(f, [7, 3])
+        e = signal.medfilt2d(np.array(f, np.float), [7, 3])
+        assert_array_equal(d, [[ 0, 50, 50, 50, 42, 15, 15, 18, 27,  0],
+                               [ 0, 50, 50, 50, 50, 42, 19, 21, 29,  0],
+                               [50, 50, 50, 50, 50, 47, 34, 34, 46, 35],
+                               [50, 50, 50, 50, 50, 50, 42, 47, 64, 42],
+                               [50, 50, 50, 50, 50, 50, 46, 55, 64, 35],
+                               [33, 50, 50, 50, 50, 47, 46, 43, 55, 26],
+                               [32, 50, 50, 50, 50, 47, 46, 45, 55, 26],
+                               [ 7, 46, 50, 50, 47, 46, 46, 43, 45, 21],
+                               [ 0, 32, 33, 39, 32, 32, 43, 43, 43,  0],
+                               [ 0,  7, 11,  7,  4,  4, 19, 19, 24,  0]])
+        assert_array_equal(d, e)
 
 class TestWiener(TestCase):
     def test_basic(self):



More information about the Scipy-svn mailing list