Skip to content

Commit

Permalink
ENH: min(max)imum_filter1d using the MINLIST algorithm, fixes scipy#3898
Browse files Browse the repository at this point in the history 


This PR implements the minimum_filter1d and maximum_filter1d functions
using Richard Harter's implementation of the MINLIST algorithm. While
slower (~30%) than the MINLINE algorithm currently in place for random
inputs, its performance is consistent for all inputs, and is much
faster for long, ordered seqeunces.
  • Loading branch information
@jaimefrio
jaimefrio committed Sep 1, 2014
1 parent 6761b2c commit 598c922
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 90 deletions.
22 changes: 22 additions & 0 deletions scipy/ndimage/filters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -790,6 +790,17 @@ def minimum_filter1d(input, size, axis=-1, output=None,
%(mode)s
%(cval)s
%(origin)s
Notes
-----
This function implements the MINLIST algorithm [1]_, as described by
Richard Harter [2]_, and has a guaranteed O(n) performance, `n` being
the `input` length, regardless of filter size.
References
----------
.. [1] http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.2777
.. [2] http://www.richardhartersworld.com/cri/2001/slidingmin.html
"""
input = numpy.asarray(input)
if numpy.iscomplexobj(input):
Expand Down Expand Up @@ -831,6 +842,17 @@ def maximum_filter1d(input, size, axis=-1, output=None,
Maximum-filtered array with same shape as input.
None if `output` is not None
Notes
-----
This function implements the MAXLIST algorithm [1]_, as described by
Richard Harter [2]_, and has a guaranteed O(n) performance, `n` being
the `input` length, regardless of filter size.
References
----------
.. [1] http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.2777
.. [2] http://www.richardhartersworld.com/cri/2001/slidingmin.html
"""
input = numpy.asarray(input)
if numpy.iscomplexobj(input):
Expand Down
159 changes: 69 additions & 90 deletions scipy/ndimage/src/ni_filters.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,100 +335,57 @@ NI_UniformFilter1D(PyArrayObject *input, npy_intp filter_size,
return PyErr_Occurred() ? 0 : 1;
}

#define WINDOW_MIN(ptr, win_size, min, min_off) \
do { \
npy_intp j_; \
min = *ptr++; \
min_off = 0; \
for (j_ = 1; j_ < win_size; j_++) { \
double val = *ptr++; \
if (val <= min) { \
min = val; \
min_off = 0; \
} \
else { \
min_off++; \
} \
} \
} while (0)

#define UPDATE_WINDOW_MIN(ptr, win_size, min, min_off) \
do { \
double val = *ptr++; \
if (val <= min) { \
min = val; \
min_off = 0; \
} \
else { \
min_off++; \
if (min_off >= win_size) { \
double *tmp = ptr - win_size; \
WINDOW_MIN(tmp, win_size, min, min_off); \
} \
} \
} while (0)

#define WINDOW_MAX(ptr, win_size, max, max_off) \
do { \
npy_intp j_; \
max = *ptr++; \
max_off = 0; \
for (j_ = 1; j_ < win_size; j_++) { \
double val = *ptr++; \
if (val >= max) { \
max = val; \
max_off = 0; \
} \
else { \
max_off++; \
} \
} \
} while (0)

#define UPDATE_WINDOW_MAX(ptr, win_size, max, max_off) \
do { \
double val = *ptr++; \
if (val >= max) { \
max = val; \
max_off = 0; \
} \
else { \
max_off++; \
if (max_off >= win_size) { \
double *tmp = ptr - win_size; \
WINDOW_MAX(tmp, win_size, max, max_off); \
} \
} \
} while (0)
#define INCREASE_RING_PTR(ptr) \
(ptr)++;\
if ((ptr) >= end) { \
(ptr) = ring; \
}

#define DECREASE_RING_PTR(ptr) \
if ((ptr) == ring) { \
(ptr) = end; \
} \
(ptr)--;

int
NI_MinOrMaxFilter1D(PyArrayObject *input, npy_intp filter_size,
int axis, PyArrayObject *output, NI_ExtendMode mode,
int axis, PyArrayObject *output, NI_ExtendMode mode,
double cval, npy_intp origin, int minimum)
{
npy_intp lines, kk, jj, ll, length, size1, size2;
npy_intp lines, kk, ll, length, size1, size2;
int more;
double *ibuffer = NULL, *obuffer = NULL;
NI_LineBuffer iline_buffer, oline_buffer;
struct pairs {
double value;
npy_intp death;
} *ring = NULL, *minpair, *end, *last;

size1 = filter_size / 2;
size2 = filter_size - size1 - 1;
/* allocate and initialize the line buffers: */
lines = -1;
if (!NI_AllocateLineBuffer(input, axis, size1 + origin, size2 - origin,
&lines, BUFFER_SIZE, &ibuffer))
&lines, BUFFER_SIZE, &ibuffer))
goto exit;
if (!NI_AllocateLineBuffer(output, axis, 0, 0, &lines, BUFFER_SIZE,
&obuffer))
&obuffer))
goto exit;
if (!NI_InitLineBuffer(input, axis, size1 + origin, size2 - origin,
lines, ibuffer, mode, cval, &iline_buffer))
lines, ibuffer, mode, cval, &iline_buffer))
goto exit;
if (!NI_InitLineBuffer(output, axis, 0, 0, lines, obuffer, mode, 0.0,
&oline_buffer))
&oline_buffer))
goto exit;
length = input->nd > 0 ? input->dimensions[axis] : 1;

/* ring is a dequeue of pairs implemented as a circular array */
ring = malloc(filter_size * sizeof(struct pairs));
if (!ring) {
goto exit;
}
end = ring + filter_size;

/* iterate over all the array lines: */
do {
/* copy lines from array to buffer: */
Expand All @@ -439,21 +396,44 @@ NI_MinOrMaxFilter1D(PyArrayObject *input, npy_intp filter_size,
/* get lines: */
double *iline = NI_GET_LINE(iline_buffer, kk);
double *oline = NI_GET_LINE(oline_buffer, kk);
double minmax;
npy_intp minmax_off;

if (minimum) {
WINDOW_MIN(iline, filter_size - 1, minmax, minmax_off);
for (ll = 0; ll < length; ll++) {
UPDATE_WINDOW_MIN(iline, filter_size, minmax, minmax_off);
*oline++ = minmax;
}

/* This check could be moved out to the Python wrapper */
if (filter_size == 1) {
memcpy(oline, iline, sizeof(double) * length);
}
else {
WINDOW_MAX(iline, filter_size - 1, minmax, minmax_off);
for (ll = 0; ll < length; ll++) {
UPDATE_WINDOW_MAX(iline, filter_size, minmax, minmax_off);
*oline++ = minmax;
/*
* Original code by Richard Harter, adapted from:
* http://www.richardhartersworld.com/cri/2001/slidingmin.html
*/
minpair = ring;
minpair->value = *iline++;
minpair->death = filter_size;
last = ring;

for (ll = 1; ll < filter_size + length - 1; ll++) {
double val = *iline++;
if (minpair->death == ll) {
INCREASE_RING_PTR(minpair)
}
if ((minimum && val <= minpair->value) ||
(!minimum && val >= minpair->value)) {
minpair->value = val;
minpair->death = ll + filter_size;
last = minpair;
}
else {
while ((minimum && last->value >= val) ||
(!minimum && last->value <= val)) {
DECREASE_RING_PTR(last)
}
INCREASE_RING_PTR(last)
last->value = val;
last->death = ll + filter_size;
}
if (ll >= filter_size - 1) {
*oline++ = minpair->value;
}
}
}
}
Expand All @@ -463,15 +443,14 @@ NI_MinOrMaxFilter1D(PyArrayObject *input, npy_intp filter_size,
} while(more);

exit:
if (ibuffer) free(ibuffer);
if (obuffer) free(obuffer);
free(ring);
free(ibuffer);
free(obuffer);
return PyErr_Occurred() ? 0 : 1;
}

#undef UPDATE_WINDOW_MAX
#undef WINDOW_MAX
#undef UPDATE_WINDOW_MIN
#undef WINDOW_MIN
#undef DECREASE_RING_PTR
#undef INCREASE_RING_PTR


#define CASE_MIN_OR_MAX_POINT(_pi, _offsets, _filter_size, _cval, \
Expand Down
34 changes: 34 additions & 0 deletions scipy/ndimage/tests/test_filters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,3 +94,37 @@ def test_gaussian_truncate():
nonzero_indices = np.where(y != 0)[0]
n = nonzero_indices.ptp() + 1
assert_equal(n, 15)

def test_minmaximum_filter1d():
# Regression gh-3898
in_ = np.arange(10)
out = sndi.minimum_filter1d(in_, 1)
assert_equal(in_, out)
out = sndi.maximum_filter1d(in_, 1)
assert_equal(in_, out)
# Test reflect
out = sndi.minimum_filter1d(in_, 5, mode='reflect')
assert_equal([0, 0, 0, 1, 2, 3, 4, 5, 6, 7], out)
out = sndi.maximum_filter1d(in_, 5, mode='reflect')
assert_equal([2, 3, 4, 5, 6, 7, 8, 9, 9, 9], out)
#Test constant
out = sndi.minimum_filter1d(in_, 5, mode='constant', cval=-1)
assert_equal([-1, -1, 0, 1, 2, 3, 4, 5, -1, -1], out)
out = sndi.maximum_filter1d(in_, 5, mode='constant', cval=10)
assert_equal([10, 10, 4, 5, 6, 7, 8, 9, 10, 10], out)
# Test nearest
out = sndi.minimum_filter1d(in_, 5, mode='nearest')
assert_equal([0, 0, 0, 1, 2, 3, 4, 5, 6, 7], out)
out = sndi.maximum_filter1d(in_, 5, mode='nearest')
assert_equal([2, 3, 4, 5, 6, 7, 8, 9, 9, 9], out)
# Test wrap
out = sndi.minimum_filter1d(in_, 5, mode='wrap')
assert_equal([0, 0, 0, 1, 2, 3, 4, 5, 0, 0], out)
out = sndi.maximum_filter1d(in_, 5, mode='wrap')
assert_equal([9, 9, 4, 5, 6, 7, 8, 9, 9, 9], out)






0 comments on commit 598c922

Please sign in to comment.