Source code for pydl.median
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
[docs]
def median(array, width=None, axis=None, even=False):
"""Replicate the IDL ``MEDIAN()`` function.
Parameters
----------
array : array-like
Compute the median of this array.
width : :class:`int`, optional
Size of the neighborhood in which to compute the median (*i.e.*,
perform median filtering). If omitted, the median of the whole
array is returned.
axis : :class:`int`, optional
Compute the median over this axis for a multi-dimensional array. If
ommitted, the median over the entire array will be returned. If
set, this function will behave as though `even` is ``True``.
even : :class:`bool`, optional
If set to ``True``, the median of arrays with an even number of elements
will be the average of the middle two values.
Returns
-------
array-like
The median of the array.
Raises
------
:exc:`ValueError`
If `width` is set, and the input `array` is not 1 or 2 dimensional.
Notes
-----
* For arrays with an even number of elements, the :func:`numpy.median`
function behaves like ``MEDIAN(array, /EVEN)``, so the absence of
the `even` keyword has to turn *off* that behavior.
* For median filtering, this uses :func:`scipy.signal.medfilt` and
:func:`scipy.signal.medfilt2d` under the hood, but patches up the
values on the array boundaries to match the return values of the
IDL ``MEDIAN()`` function.
References
----------
https://www.nv5geospatialsoftware.com/docs/median.html
"""
import numpy as np
from scipy.signal import medfilt, medfilt2d
if width is None:
if axis is None:
f = array.flatten()
if f.size % 2 == 1 or even:
return np.median(array)
else:
i = f.argsort()
return f[i[f.size//2]]
else:
return np.median(array, axis=axis)
else:
if array.ndim == 1:
medarray = medfilt(array, min(width, array.size))
istart = int((width - 1)/2)
iend = array.size - int((width + 1)/2)
i = np.arange(array.size)
w = (i < istart) | (i > iend)
medarray[w] = array[w]
return medarray
elif array.ndim == 2:
medarray = medfilt2d(array, min(width, array.size))
istart = int((width-1)/2)
iend = (array.shape[0] - int((width+1)/2), array.shape[1] - int((width+1)/2))
i = np.arange(array.shape[0])
j = np.arange(array.shape[1])
w = ((i < istart) | (i > iend[0]), (j < istart) | (j > iend[1]))
medarray[w[0], :] = array[w[0], :]
medarray[:, w[1]] = array[:, w[1]]
return medarray
else:
raise ValueError("Invalid number of dimensions for input array!")