Skip to content
This repository
Browse code

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).
  • Loading branch information...
commit abfc8c2ec3d8b011a840b542436e0a372a3282a8 1 parent c8ce689
Jarrod Millman jarrodmillman authored

Showing 2 changed files with 129 additions and 292 deletions. Show diff stats Hide diff stats

  1. +105 289 scipy/signal/medianfilter.c
  2. +24 3 scipy/signal/tests/test_signaltools.py
394 scipy/signal/medianfilter.c
... ... @@ -1,311 +1,127 @@
1 1
2 2 /*--------------------------------------------------------------------*/
3 3
4   -/*
5   - * This Quickselect routine is based on the algorithm described in
6   - * "Numerical recipies in C", Second Edition,
7   - * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
8   - */
9   -
10 4 #include "Python.h"
11 5 #define NO_IMPORT_ARRAY
12 6 #include "numpy/noprefix.h"
13 7
  8 +
  9 +/* defined below */
14 10 void f_medfilt2(float*,float*,intp*,intp*);
15 11 void d_medfilt2(double*,double*,intp*,intp*);
16 12 void b_medfilt2(unsigned char*,unsigned char*,intp*,intp*);
17 13 extern char *check_malloc (int);
18 14
19   -#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
20   -
21   -float f_quick_select(float arr[], int n)
22   -{
23   - int low, high ;
24   - int median;
25   - int middle, ll, hh;
26   -
27   - low = 0 ; high = n-1 ; median = (low + high) / 2;
28   - for (;;) {
29   - if (high <= low) /* One element only */
30   - return arr[median] ;
31   -
32   - if (high == low + 1) { /* Two elements only */
33   - if (arr[low] > arr[high])
34   - ELEM_SWAP(arr[low], arr[high]) ;
35   - return arr[median] ;
36   - }
37   -
38   - /* Find median of low, middle and high items; swap into position low */
39   - middle = (low + high) / 2;
40   - if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
41   - if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
42   - if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
43   -
44   - /* Swap low item (now in position middle) into position (low+1) */
45   - ELEM_SWAP(arr[middle], arr[low+1]) ;
46   -
47   - /* Nibble from each end towards middle, swapping items when stuck */
48   - ll = low + 1;
49   - hh = high;
50   - for (;;) {
51   - do ll++; while (arr[low] > arr[ll]) ;
52   - do hh--; while (arr[hh] > arr[low]) ;
53   -
54   - if (hh < ll)
55   - break;
56   -
57   - ELEM_SWAP(arr[ll], arr[hh]) ;
58   - }
59   -
60   - /* Swap middle item (in position low) back into correct position */
61   - ELEM_SWAP(arr[low], arr[hh]) ;
62   -
63   - /* Re-set active partition */
64   - if (hh <= median)
65   - low = ll;
66   - if (hh >= median)
67   - high = hh - 1;
68   - }
69   -}
70   -
71   -#undef ELEM_SWAP
72   -
73   -
74   -#define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
75   -
76   -double d_quick_select(double arr[], int n)
77   -{
78   - int low, high ;
79   - int median;
80   - int middle, ll, hh;
81   -
82   - low = 0 ; high = n-1 ; median = (low + high) / 2;
83   - for (;;) {
84   - if (high <= low) /* One element only */
85   - return arr[median] ;
86   -
87   - if (high == low + 1) { /* Two elements only */
88   - if (arr[low] > arr[high])
89   - ELEM_SWAP(arr[low], arr[high]) ;
90   - return arr[median] ;
91   - }
92   -
93   - /* Find median of low, middle and high items; swap into position low */
94   - middle = (low + high) / 2;
95   - if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
96   - if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
97   - if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
98   -
99   - /* Swap low item (now in position middle) into position (low+1) */
100   - ELEM_SWAP(arr[middle], arr[low+1]) ;
101   -
102   - /* Nibble from each end towards middle, swapping items when stuck */
103   - ll = low + 1;
104   - hh = high;
105   - for (;;) {
106   - do ll++; while (arr[low] > arr[ll]) ;
107   - do hh--; while (arr[hh] > arr[low]) ;
108   -
109   - if (hh < ll)
110   - break;
111   -
112   - ELEM_SWAP(arr[ll], arr[hh]) ;
113   - }
114   -
115   - /* Swap middle item (in position low) back into correct position */
116   - ELEM_SWAP(arr[low], arr[hh]) ;
117   -
118   - /* Re-set active partition */
119   - if (hh <= median)
120   - low = ll;
121   - if (hh >= median)
122   - high = hh - 1;
123   - }
124   -}
125   -
126   -#undef ELEM_SWAP
127   -
128   -#define ELEM_SWAP(a,b) { register unsigned char t=(a);(a)=(b);(b)=t; }
129   -
130   -unsigned char b_quick_select(unsigned char arr[], int n)
131   -{
132   - int low, high ;
133   - int median;
134   - int middle, ll, hh;
135   -
136   - low = 0 ; high = n-1 ; median = (low + high) / 2;
137   - for (;;) {
138   - if (high <= low) /* One element only */
139   - return arr[median] ;
140   -
141   - if (high == low + 1) { /* Two elements only */
142   - if (arr[low] > arr[high])
143   - ELEM_SWAP(arr[low], arr[high]) ;
144   - return arr[median] ;
145   - }
146   -
147   - /* Find median of low, middle and high items; swap into position low */
148   - middle = (low + high) / 2;
149   - if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
150   - if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
151   - if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
152   -
153   - /* Swap low item (now in position middle) into position (low+1) */
154   - ELEM_SWAP(arr[middle], arr[low+1]) ;
155   -
156   - /* Nibble from each end towards middle, swapping items when stuck */
157   - ll = low + 1;
158   - hh = high;
159   - for (;;) {
160   - do ll++; while (arr[low] > arr[ll]) ;
161   - do hh--; while (arr[hh] > arr[low]) ;
162   -
163   - if (hh < ll)
164   - break;
165   -
166   - ELEM_SWAP(arr[ll], arr[hh]) ;
167   - }
168   -
169   - /* Swap middle item (in position low) back into correct position */
170   - ELEM_SWAP(arr[low], arr[hh]) ;
171   -
172   - /* Re-set active partition */
173   - if (hh <= median)
174   - low = ll;
175   - if (hh >= median)
176   - high = hh - 1;
177   - }
178   -}
179   -
180   -#undef ELEM_SWAP
181 15
182   -/* 2-D median filter with zero-padding on edges. */
183   -void d_medfilt2(double* in, double* out, intp* Nwin, intp* Ns)
184   -{
185   - int nx, ny, hN[2];
186   - int pre_x, pre_y, pos_x, pos_y;
187   - int subx, suby, k, totN;
188   - double *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
189   -
190   - totN = Nwin[0] * Nwin[1];
191   - myvals = (double *) check_malloc( totN * sizeof(double));
192   -
193   - hN[0] = Nwin[0] >> 1;
194   - hN[1] = Nwin[1] >> 1;
195   - ptr1 = in;
196   - fptr1 = out;
197   - for (ny = 0; ny < Ns[0]; ny++)
198   - for (nx = 0; nx < Ns[1]; nx++) {
199   - pre_x = hN[1];
200   - pre_y = hN[0];
201   - pos_x = hN[1];
202   - pos_y = hN[0];
203   - if (nx < hN[1]) pre_x = nx;
204   - if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
205   - if (ny < hN[0]) pre_y = ny;
206   - if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
207   - fptr2 = myvals;
208   - ptr2 = ptr1 - pre_x - pre_y*Ns[1];
209   - for (suby = -pre_y; suby <= pos_y; suby++) {
210   - for (subx = -pre_x; subx <= pos_x; subx++)
211   - *fptr2++ = *ptr2++;
212   - ptr2 += Ns[1] - (pre_x + pos_x + 1);
213   - }
214   - ptr1++;
215   -
216   - /* Zero pad */
217   - for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
218   - *fptr2++ = 0.0;
  16 +/* The QUICK_SELECT routine is based on Hoare's Quickselect algorithm,
  17 + * with unrolled recursion.
  18 + * Author: Thouis R. Jones, 2008
  19 + */
219 20
220   - /* *fptr1++ = median(myvals,totN); */
221   - *fptr1++ = d_quick_select(myvals,totN);
222   - }
  21 +#define ELEM_SWAP(t, a, x, y) {register t temp = (a)[x]; (a)[x] = (a)[y]; (a)[y] = temp;}
  22 +#define FIRST_LOWEST(x, y, z) (((x) < (y)) && ((x) < (z)))
  23 +#define FIRST_HIGHEST(x, y, z) (((x) > (y)) && ((x) > (z)))
  24 +#define LOWEST_IDX(a, x, y) (((a)[x] < (a)[y]) ? (x) : (y))
  25 +#define HIGHEST_IDX(a, x, y) (((a)[x] > (a)[y]) ? (x) : (y))
  26 +
  27 +/* 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 */
  28 +#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)))
  29 +
  30 +#define QUICK_SELECT(NAME, TYPE) \
  31 +TYPE NAME(TYPE arr[], int n) \
  32 +{ \
  33 + int lo, hi, mid, md; \
  34 + int median_idx; \
  35 + int ll, hh; \
  36 + TYPE piv; \
  37 + \
  38 + lo = 0; hi = n-1; \
  39 + median_idx = (n - 1) / 2; /* lower of middle values for even-length arrays */ \
  40 + \
  41 + while (1) { \
  42 + if ((hi - lo) < 2) { \
  43 + if (arr[hi] < arr[lo]) ELEM_SWAP(TYPE, arr, lo, hi); \
  44 + return arr[median_idx]; \
  45 + } \
  46 + \
  47 + mid = (hi + lo) / 2; \
  48 + /* put the median of lo,mid,hi at position lo - this will be the pivot */ \
  49 + md = MEDIAN_IDX(arr, lo, mid, hi); \
  50 + ELEM_SWAP(TYPE, arr, lo, md); \
  51 + \
  52 + /* Nibble from each end towards middle, swapping misordered items */ \
  53 + piv = arr[lo]; \
  54 + for (ll = lo+1, hh = hi;; ll++, hh--) { \
  55 + while (arr[ll] < piv) ll++; \
  56 + while (arr[hh] > piv) hh--; \
  57 + if (hh < ll) break; \
  58 + ELEM_SWAP(TYPE, arr, ll, hh); \
  59 + } \
  60 + /* move pivot to top of lower partition */ \
  61 + ELEM_SWAP(TYPE, arr, hh, lo); \
  62 + /* set lo, hi for new range to search */ \
  63 + if (hh < median_idx) /* search upper partition */ \
  64 + lo = hh+1; \
  65 + else if (hh > median_idx) /* search lower partition */ \
  66 + hi = hh-1; \
  67 + else \
  68 + return piv; \
  69 + } \
223 70 }
224 71
225 72
226 73 /* 2-D median filter with zero-padding on edges. */
227   -void f_medfilt2(float* in, float* out, intp* Nwin, intp* Ns)
228   -{
229   - int nx, ny, hN[2];
230   - int pre_x, pre_y, pos_x, pos_y;
231   - int subx, suby, k, totN;
232   - float *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
233   -
234   - totN = Nwin[0] * Nwin[1];
235   - myvals = (float *) check_malloc( totN * sizeof(float));
236   -
237   - hN[0] = Nwin[0] >> 1;
238   - hN[1] = Nwin[1] >> 1;
239   - ptr1 = in;
240   - fptr1 = out;
241   - for (ny = 0; ny < Ns[0]; ny++)
242   - for (nx = 0; nx < Ns[1]; nx++) {
243   - pre_x = hN[1];
244   - pre_y = hN[0];
245   - pos_x = hN[1];
246   - pos_y = hN[0];
247   - if (nx < hN[1]) pre_x = nx;
248   - if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
249   - if (ny < hN[0]) pre_y = ny;
250   - if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
251   - fptr2 = myvals;
252   - ptr2 = ptr1 - pre_x - pre_y*Ns[1];
253   - for (suby = -pre_y; suby <= pos_y; suby++) {
254   - for (subx = -pre_x; subx <= pos_x; subx++)
255   - *fptr2++ = *ptr2++;
256   - ptr2 += Ns[1] - (pre_x + pos_x + 1);
257   - }
258   - ptr1++;
259   -
260   - /* Zero pad */
261   - for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
262   - *fptr2++ = 0.0;
263   -
264   - /* *fptr1++ = median(myvals,totN); */
265   - *fptr1++ = f_quick_select(myvals,totN);
266   - }
  74 +#define MEDIAN_FILTER_2D(NAME, TYPE, SELECT) \
  75 +void NAME(TYPE* in, TYPE* out, intp* Nwin, intp* Ns) \
  76 +{ \
  77 + int nx, ny, hN[2]; \
  78 + int pre_x, pre_y, pos_x, pos_y; \
  79 + int subx, suby, k, totN; \
  80 + TYPE *myvals, *fptr1, *fptr2, *ptr1, *ptr2; \
  81 + \
  82 + totN = Nwin[0] * Nwin[1]; \
  83 + myvals = (TYPE *) check_malloc( totN * sizeof(TYPE)); \
  84 + \
  85 + hN[0] = Nwin[0] >> 1; \
  86 + hN[1] = Nwin[1] >> 1; \
  87 + ptr1 = in; \
  88 + fptr1 = out; \
  89 + for (ny = 0; ny < Ns[0]; ny++) \
  90 + for (nx = 0; nx < Ns[1]; nx++) { \
  91 + pre_x = hN[1]; \
  92 + pre_y = hN[0]; \
  93 + pos_x = hN[1]; \
  94 + pos_y = hN[0]; \
  95 + if (nx < hN[1]) pre_x = nx; \
  96 + if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1; \
  97 + if (ny < hN[0]) pre_y = ny; \
  98 + if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1; \
  99 + fptr2 = myvals; \
  100 + ptr2 = ptr1 - pre_x - pre_y*Ns[1]; \
  101 + for (suby = -pre_y; suby <= pos_y; suby++) { \
  102 + for (subx = -pre_x; subx <= pos_x; subx++) \
  103 + *fptr2++ = *ptr2++; \
  104 + ptr2 += Ns[1] - (pre_x + pos_x + 1); \
  105 + } \
  106 + ptr1++; \
  107 + \
  108 + /* Zero pad */ \
  109 + for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++) \
  110 + *fptr2++ = 0.0; \
  111 + \
  112 + /* *fptr1++ = median(myvals,totN); */ \
  113 + *fptr1++ = SELECT(myvals,totN); \
  114 + } \
  115 + free(myvals); \
267 116 }
268 117
269 118
270   -/* 2-D median filter with zero-padding on edges. */
271   -void b_medfilt2(unsigned char *in, unsigned char *out, intp* Nwin, intp* Ns)
272   -{
273   - int nx, ny, hN[2];
274   - int pre_x, pre_y, pos_x, pos_y;
275   - int subx, suby, k, totN;
276   - unsigned char *myvals, *fptr1, *fptr2, *ptr1, *ptr2;
277   -
278   - totN = Nwin[0] * Nwin[1];
279   - myvals = (unsigned char *) check_malloc( totN * sizeof(unsigned char));
  119 +/* define quick_select for floats, doubles, and unsigned characters */
  120 +QUICK_SELECT(f_quick_select, float)
  121 +QUICK_SELECT(d_quick_select, double)
  122 +QUICK_SELECT(b_quick_select, unsigned char)
280 123
281   - hN[0] = Nwin[0] >> 1;
282   - hN[1] = Nwin[1] >> 1;
283   - ptr1 = in;
284   - fptr1 = out;
285   - for (ny = 0; ny < Ns[0]; ny++)
286   - for (nx = 0; nx < Ns[1]; nx++) {
287   - pre_x = hN[1];
288   - pre_y = hN[0];
289   - pos_x = hN[1];
290   - pos_y = hN[0];
291   - if (nx < hN[1]) pre_x = nx;
292   - if (nx >= Ns[1] - hN[1]) pos_x = Ns[1] - nx - 1;
293   - if (ny < hN[0]) pre_y = ny;
294   - if (ny >= Ns[0] - hN[0]) pos_y = Ns[0] - ny - 1;
295   - fptr2 = myvals;
296   - ptr2 = ptr1 - pre_x - pre_y*Ns[1];
297   - for (suby = -pre_y; suby <= pos_y; suby++) {
298   - for (subx = -pre_x; subx <= pos_x; subx++)
299   - *fptr2++ = *ptr2++;
300   - ptr2 += Ns[1] - (pre_x + pos_x + 1);
301   - }
302   - ptr1++;
303   -
304   - /* Zero pad */
305   - for (k = (pre_x + pos_x + 1)*(pre_y + pos_y + 1); k < totN; k++)
306   - *fptr2++ = 0.0;
307   -
308   - /* *fptr1++ = median(myvals,totN); */
309   - *fptr1++ = b_quick_select(myvals,totN);
310   - }
311   -}
  124 +/* define medfilt for floats, doubles, and unsigned characters */
  125 +MEDIAN_FILTER_2D(f_medfilt2, float, f_quick_select)
  126 +MEDIAN_FILTER_2D(d_medfilt2, double, d_quick_select)
  127 +MEDIAN_FILTER_2D(b_medfilt2, unsigned char, b_quick_select)
27 scipy/signal/tests/test_signaltools.py
@@ -28,9 +28,30 @@ def test_complex(self):
28 28
29 29 class TestMedFilt(TestCase):
30 30 def test_basic(self):
31   - f = [[3,4,5],[2,3,4],[1,2,5]]
32   - d = signal.medfilt(f)
33   - assert_array_equal(d, [[0,3,0],[2,3,3],[0,2,0]])
  31 + f = [[50, 50, 50, 50, 50, 92, 18, 27, 65, 46],
  32 + [50, 50, 50, 50, 50, 0, 72, 77, 68, 66],
  33 + [50, 50, 50, 50, 50, 46, 47, 19, 64, 77],
  34 + [50, 50, 50, 50, 50, 42, 15, 29, 95, 35],
  35 + [50, 50, 50, 50, 50, 46, 34, 9, 21, 66],
  36 + [70, 97, 28, 68, 78, 77, 61, 58, 71, 42],
  37 + [64, 53, 44, 29, 68, 32, 19, 68, 24, 84],
  38 + [ 3, 33, 53, 67, 1, 78, 74, 55, 12, 83],
  39 + [ 7, 11, 46, 70, 60, 47, 24, 43, 61, 26],
  40 + [32, 61, 88, 7, 39, 4, 92, 64, 45, 61]]
  41 +
  42 + d = signal.medfilt(f, [7, 3])
  43 + e = signal.medfilt2d(np.array(f, np.float), [7, 3])
  44 + assert_array_equal(d, [[ 0, 50, 50, 50, 42, 15, 15, 18, 27, 0],
  45 + [ 0, 50, 50, 50, 50, 42, 19, 21, 29, 0],
  46 + [50, 50, 50, 50, 50, 47, 34, 34, 46, 35],
  47 + [50, 50, 50, 50, 50, 50, 42, 47, 64, 42],
  48 + [50, 50, 50, 50, 50, 50, 46, 55, 64, 35],
  49 + [33, 50, 50, 50, 50, 47, 46, 43, 55, 26],
  50 + [32, 50, 50, 50, 50, 47, 46, 45, 55, 26],
  51 + [ 7, 46, 50, 50, 47, 46, 46, 43, 45, 21],
  52 + [ 0, 32, 33, 39, 32, 32, 43, 43, 43, 0],
  53 + [ 0, 7, 11, 7, 4, 4, 19, 19, 24, 0]])
  54 + assert_array_equal(d, e)
34 55
35 56 class TestWiener(TestCase):
36 57 def test_basic(self):

0 comments on commit abfc8c2

Please sign in to comment.
Something went wrong with that request. Please try again.