Index: scipy/signal/medianfilter.c =================================================================== --- scipy/signal/medianfilter.c (revision 5274) +++ scipy/signal/medianfilter.c (working copy) @@ -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) Index: scipy/signal/tests/test_signaltools.py =================================================================== --- scipy/signal/tests/test_signaltools.py (revision 5274) +++ scipy/signal/tests/test_signaltools.py (working copy) @@ -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):