import functools
import numpy as np
import scipy.fft as fft
from .._shared.utils import _supported_float_type
def _get_nd_butterworth_filter(
shape, factor, order, high_pass, real, dtype=np.float64, squared_butterworth=True
):
"""Create a N-dimensional Butterworth mask for an FFT
Parameters
----------
shape : tuple of int
Shape of the n-dimensional FFT and mask.
factor : float
Fraction of mask dimensions where the cutoff should be.
order : float
Controls the slope in the cutoff region.
high_pass : bool
Whether the filter is high pass (low frequencies attenuated) or
low pass (high frequencies are attenuated).
real : bool
Whether the FFT is of a real (True) or complex (False) image
squared_butterworth : bool, optional
When True, the square of the Butterworth filter is used.
Returns
-------
wfilt : ndarray
The FFT mask.
"""
ranges = []
for i, d in enumerate(shape):
# start and stop ensures center of mask aligns with center of FFT
axis = np.arange(-(d - 1) // 2, (d - 1) // 2 + 1) / (d * factor)
ranges.append(fft.ifftshift(axis**2))
# for real image FFT, halve the last axis
if real:
limit = d // 2 + 1
ranges[-1] = ranges[-1][:limit]
# q2 = squared Euclidean distance grid
q2 = functools.reduce(np.add, np.meshgrid(*ranges, indexing="ij", sparse=True))
q2 = q2.astype(dtype)
q2 = np.power(q2, order)
wfilt = 1 / (1 + q2)
if high_pass:
wfilt *= q2
if not squared_butterworth:
np.sqrt(wfilt, out=wfilt)
return wfilt
[文档]
def butterworth(
image,
cutoff_frequency_ratio=0.005,
high_pass=True,
order=2.0,
channel_axis=None,
*,
squared_butterworth=True,
npad=0,
):
"""Apply a Butterworth filter to enhance high or low frequency features.
This filter is defined in the Fourier domain.
Parameters
----------
image : (M[, N[, ..., P]][, C]) ndarray
Input image.
cutoff_frequency_ratio : float, optional
Determines the position of the cut-off relative to the shape of the
FFT. Receives a value between [0, 0.5].
high_pass : bool, optional
Whether to perform a high pass filter. If False, a low pass filter is
performed.
order : float, optional
Order of the filter which affects the slope near the cut-off. Higher
order means steeper slope in frequency space.
channel_axis : int, optional
If there is a channel dimension, provide the index here. If None
(default) then all axes are assumed to be spatial dimensions.
squared_butterworth : bool, optional
When True, the square of a Butterworth filter is used. See notes below
for more details.
npad : int, optional
Pad each edge of the image by `npad` pixels using `numpy.pad`'s
``mode='edge'`` extension.
Returns
-------
result : ndarray
The Butterworth-filtered image.
Notes
-----
A band-pass filter can be achieved by combining a high-pass and low-pass
filter. The user can increase `npad` if boundary artifacts are apparent.
The "Butterworth filter" used in image processing textbooks (e.g. [1]_,
[2]_) is often the square of the traditional Butterworth filters as
described by [3]_, [4]_. The squared version will be used here if
`squared_butterworth` is set to ``True``. The lowpass, squared Butterworth
filter is given by the following expression for the lowpass case:
.. math::
H_{low}(f) = \\frac{1}{1 + \\left(\\frac{f}{c f_s}\\right)^{2n}}
with the highpass case given by
.. math::
H_{hi}(f) = 1 - H_{low}(f)
where :math:`f=\\sqrt{\\sum_{d=0}^{\\mathrm{ndim}} f_{d}^{2}}` is the
absolute value of the spatial frequency, :math:`f_s` is the sampling
frequency, :math:`c` the ``cutoff_frequency_ratio``, and :math:`n` is the
filter `order` [1]_. When ``squared_butterworth=False``, the square root of
the above expressions are used instead.
Note that ``cutoff_frequency_ratio`` is defined in terms of the sampling
frequency, :math:`f_s`. The FFT spectrum covers the Nyquist range
(:math:`[-f_s/2, f_s/2]`) so ``cutoff_frequency_ratio`` should have a value
between 0 and 0.5. The frequency response (gain) at the cutoff is 0.5 when
``squared_butterworth`` is true and :math:`1/\\sqrt{2}` when it is false.
Examples
--------
Apply a high-pass and low-pass Butterworth filter to a grayscale and
color image respectively:
>>> from skimage.data import camera, astronaut
>>> from skimage.filters import butterworth
>>> high_pass = butterworth(camera(), 0.07, True, 8)
>>> low_pass = butterworth(astronaut(), 0.01, False, 4, channel_axis=-1)
References
----------
.. [1] Russ, John C., et al. The Image Processing Handbook, 3rd. Ed.
1999, CRC Press, LLC.
.. [2] Birchfield, Stan. Image Processing and Analysis. 2018. Cengage
Learning.
.. [3] Butterworth, Stephen. "On the theory of filter amplifiers."
Wireless Engineer 7.6 (1930): 536-541.
.. [4] https://en.wikipedia.org/wiki/Butterworth_filter
"""
if npad < 0:
raise ValueError("npad must be >= 0")
elif npad > 0:
center_slice = tuple(slice(npad, s + npad) for s in image.shape)
image = np.pad(image, npad, mode='edge')
fft_shape = (
image.shape if channel_axis is None else np.delete(image.shape, channel_axis)
)
is_real = np.isrealobj(image)
float_dtype = _supported_float_type(image.dtype, allow_complex=True)
if cutoff_frequency_ratio < 0 or cutoff_frequency_ratio > 0.5:
raise ValueError("cutoff_frequency_ratio should be in the range [0, 0.5]")
wfilt = _get_nd_butterworth_filter(
fft_shape,
cutoff_frequency_ratio,
order,
high_pass,
is_real,
float_dtype,
squared_butterworth,
)
axes = np.arange(image.ndim)
if channel_axis is not None:
axes = np.delete(axes, channel_axis)
abs_channel = channel_axis % image.ndim
post = image.ndim - abs_channel - 1
sl = (slice(None),) * abs_channel + (np.newaxis,) + (slice(None),) * post
wfilt = wfilt[sl]
if is_real:
butterfilt = fft.irfftn(
wfilt * fft.rfftn(image, axes=axes), s=fft_shape, axes=axes
)
else:
butterfilt = fft.ifftn(
wfilt * fft.fftn(image, axes=axes), s=fft_shape, axes=axes
)
if npad > 0:
butterfilt = butterfilt[center_slice]
return butterfilt
