Coverage for /usr/lib/python3/dist-packages/scipy/ndimage/_filters.py: 13%
494 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1# Copyright (C) 2003-2005 Peter J. Verveer
2#
3# Redistribution and use in source and binary forms, with or without
4# modification, are permitted provided that the following conditions
5# are met:
6#
7# 1. Redistributions of source code must retain the above copyright
8# notice, this list of conditions and the following disclaimer.
9#
10# 2. Redistributions in binary form must reproduce the above
11# copyright notice, this list of conditions and the following
12# disclaimer in the documentation and/or other materials provided
13# with the distribution.
14#
15# 3. The name of the author may not be used to endorse or promote
16# products derived from this software without specific prior
17# written permission.
18#
19# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
20# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
23# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
25# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
27# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31from collections.abc import Iterable
32import numbers
33import warnings
34import numpy
35import operator
36from numpy.core.multiarray import normalize_axis_index
37from . import _ni_support
38from . import _nd_image
39from . import _ni_docstrings
41__all__ = ['correlate1d', 'convolve1d', 'gaussian_filter1d', 'gaussian_filter',
42 'prewitt', 'sobel', 'generic_laplace', 'laplace',
43 'gaussian_laplace', 'generic_gradient_magnitude',
44 'gaussian_gradient_magnitude', 'correlate', 'convolve',
45 'uniform_filter1d', 'uniform_filter', 'minimum_filter1d',
46 'maximum_filter1d', 'minimum_filter', 'maximum_filter',
47 'rank_filter', 'median_filter', 'percentile_filter',
48 'generic_filter1d', 'generic_filter']
51def _invalid_origin(origin, lenw):
52 return (origin < -(lenw // 2)) or (origin > (lenw - 1) // 2)
55def _complex_via_real_components(func, input, weights, output, cval, **kwargs):
56 """Complex convolution via a linear combination of real convolutions."""
57 complex_input = input.dtype.kind == 'c'
58 complex_weights = weights.dtype.kind == 'c'
59 if complex_input and complex_weights:
60 # real component of the output
61 func(input.real, weights.real, output=output.real,
62 cval=numpy.real(cval), **kwargs)
63 output.real -= func(input.imag, weights.imag, output=None,
64 cval=numpy.imag(cval), **kwargs)
65 # imaginary component of the output
66 func(input.real, weights.imag, output=output.imag,
67 cval=numpy.real(cval), **kwargs)
68 output.imag += func(input.imag, weights.real, output=None,
69 cval=numpy.imag(cval), **kwargs)
70 elif complex_input:
71 func(input.real, weights, output=output.real, cval=numpy.real(cval),
72 **kwargs)
73 func(input.imag, weights, output=output.imag, cval=numpy.imag(cval),
74 **kwargs)
75 else:
76 if numpy.iscomplexobj(cval):
77 raise ValueError("Cannot provide a complex-valued cval when the "
78 "input is real.")
79 func(input, weights.real, output=output.real, cval=cval, **kwargs)
80 func(input, weights.imag, output=output.imag, cval=cval, **kwargs)
81 return output
84@_ni_docstrings.docfiller
85def correlate1d(input, weights, axis=-1, output=None, mode="reflect",
86 cval=0.0, origin=0):
87 """Calculate a 1-D correlation along the given axis.
89 The lines of the array along the given axis are correlated with the
90 given weights.
92 Parameters
93 ----------
94 %(input)s
95 weights : array
96 1-D sequence of numbers.
97 %(axis)s
98 %(output)s
99 %(mode_reflect)s
100 %(cval)s
101 %(origin)s
103 Examples
104 --------
105 >>> from scipy.ndimage import correlate1d
106 >>> correlate1d([2, 8, 0, 4, 1, 9, 9, 0], weights=[1, 3])
107 array([ 8, 26, 8, 12, 7, 28, 36, 9])
108 """
109 input = numpy.asarray(input)
110 weights = numpy.asarray(weights)
111 complex_input = input.dtype.kind == 'c'
112 complex_weights = weights.dtype.kind == 'c'
113 if complex_input or complex_weights:
114 if complex_weights:
115 weights = weights.conj()
116 weights = weights.astype(numpy.complex128, copy=False)
117 kwargs = dict(axis=axis, mode=mode, origin=origin)
118 output = _ni_support._get_output(output, input, complex_output=True)
119 return _complex_via_real_components(correlate1d, input, weights,
120 output, cval, **kwargs)
122 output = _ni_support._get_output(output, input)
123 weights = numpy.asarray(weights, dtype=numpy.float64)
124 if weights.ndim != 1 or weights.shape[0] < 1:
125 raise RuntimeError('no filter weights given')
126 if not weights.flags.contiguous:
127 weights = weights.copy()
128 axis = normalize_axis_index(axis, input.ndim)
129 if _invalid_origin(origin, len(weights)):
130 raise ValueError('Invalid origin; origin must satisfy '
131 '-(len(weights) // 2) <= origin <= '
132 '(len(weights)-1) // 2')
133 mode = _ni_support._extend_mode_to_code(mode)
134 _nd_image.correlate1d(input, weights, axis, output, mode, cval,
135 origin)
136 return output
139@_ni_docstrings.docfiller
140def convolve1d(input, weights, axis=-1, output=None, mode="reflect",
141 cval=0.0, origin=0):
142 """Calculate a 1-D convolution along the given axis.
144 The lines of the array along the given axis are convolved with the
145 given weights.
147 Parameters
148 ----------
149 %(input)s
150 weights : ndarray
151 1-D sequence of numbers.
152 %(axis)s
153 %(output)s
154 %(mode_reflect)s
155 %(cval)s
156 %(origin)s
158 Returns
159 -------
160 convolve1d : ndarray
161 Convolved array with same shape as input
163 Examples
164 --------
165 >>> from scipy.ndimage import convolve1d
166 >>> convolve1d([2, 8, 0, 4, 1, 9, 9, 0], weights=[1, 3])
167 array([14, 24, 4, 13, 12, 36, 27, 0])
168 """
169 weights = weights[::-1]
170 origin = -origin
171 if not len(weights) & 1:
172 origin -= 1
173 weights = numpy.asarray(weights)
174 if weights.dtype.kind == 'c':
175 # pre-conjugate here to counteract the conjugation in correlate1d
176 weights = weights.conj()
177 return correlate1d(input, weights, axis, output, mode, cval, origin)
180def _gaussian_kernel1d(sigma, order, radius):
181 """
182 Computes a 1-D Gaussian convolution kernel.
183 """
184 if order < 0:
185 raise ValueError('order must be non-negative')
186 exponent_range = numpy.arange(order + 1)
187 sigma2 = sigma * sigma
188 x = numpy.arange(-radius, radius+1)
189 phi_x = numpy.exp(-0.5 / sigma2 * x ** 2)
190 phi_x = phi_x / phi_x.sum()
192 if order == 0:
193 return phi_x
194 else:
195 # f(x) = q(x) * phi(x) = q(x) * exp(p(x))
196 # f'(x) = (q'(x) + q(x) * p'(x)) * phi(x)
197 # p'(x) = -1 / sigma ** 2
198 # Implement q'(x) + q(x) * p'(x) as a matrix operator and apply to the
199 # coefficients of q(x)
200 q = numpy.zeros(order + 1)
201 q[0] = 1
202 D = numpy.diag(exponent_range[1:], 1) # D @ q(x) = q'(x)
203 P = numpy.diag(numpy.ones(order)/-sigma2, -1) # P @ q(x) = q(x) * p'(x)
204 Q_deriv = D + P
205 for _ in range(order):
206 q = Q_deriv.dot(q)
207 q = (x[:, None] ** exponent_range).dot(q)
208 return q * phi_x
211@_ni_docstrings.docfiller
212def gaussian_filter1d(input, sigma, axis=-1, order=0, output=None,
213 mode="reflect", cval=0.0, truncate=4.0, *, radius=None):
214 """1-D Gaussian filter.
216 Parameters
217 ----------
218 %(input)s
219 sigma : scalar
220 standard deviation for Gaussian kernel
221 %(axis)s
222 order : int, optional
223 An order of 0 corresponds to convolution with a Gaussian
224 kernel. A positive order corresponds to convolution with
225 that derivative of a Gaussian.
226 %(output)s
227 %(mode_reflect)s
228 %(cval)s
229 truncate : float, optional
230 Truncate the filter at this many standard deviations.
231 Default is 4.0.
232 radius : None or int, optional
233 Radius of the Gaussian kernel. If specified, the size of
234 the kernel will be ``2*radius + 1``, and `truncate` is ignored.
235 Default is None.
237 Returns
238 -------
239 gaussian_filter1d : ndarray
241 Notes
242 -----
243 The Gaussian kernel will have size ``2*radius + 1`` along each axis. If
244 `radius` is None, a default ``radius = round(truncate * sigma)`` will be
245 used.
247 Examples
248 --------
249 >>> from scipy.ndimage import gaussian_filter1d
250 >>> import numpy as np
251 >>> gaussian_filter1d([1.0, 2.0, 3.0, 4.0, 5.0], 1)
252 array([ 1.42704095, 2.06782203, 3. , 3.93217797, 4.57295905])
253 >>> gaussian_filter1d([1.0, 2.0, 3.0, 4.0, 5.0], 4)
254 array([ 2.91948343, 2.95023502, 3. , 3.04976498, 3.08051657])
255 >>> import matplotlib.pyplot as plt
256 >>> rng = np.random.default_rng()
257 >>> x = rng.standard_normal(101).cumsum()
258 >>> y3 = gaussian_filter1d(x, 3)
259 >>> y6 = gaussian_filter1d(x, 6)
260 >>> plt.plot(x, 'k', label='original data')
261 >>> plt.plot(y3, '--', label='filtered, sigma=3')
262 >>> plt.plot(y6, ':', label='filtered, sigma=6')
263 >>> plt.legend()
264 >>> plt.grid()
265 >>> plt.show()
267 """
268 sd = float(sigma)
269 # make the radius of the filter equal to truncate standard deviations
270 lw = int(truncate * sd + 0.5)
271 if radius is not None:
272 lw = radius
273 if not isinstance(lw, numbers.Integral) or lw < 0:
274 raise ValueError('Radius must be a nonnegative integer.')
275 # Since we are calling correlate, not convolve, revert the kernel
276 weights = _gaussian_kernel1d(sigma, order, lw)[::-1]
277 return correlate1d(input, weights, axis, output, mode, cval, 0)
280@_ni_docstrings.docfiller
281def gaussian_filter(input, sigma, order=0, output=None,
282 mode="reflect", cval=0.0, truncate=4.0, *, radius=None,
283 axes=None):
284 """Multidimensional Gaussian filter.
286 Parameters
287 ----------
288 %(input)s
289 sigma : scalar or sequence of scalars
290 Standard deviation for Gaussian kernel. The standard
291 deviations of the Gaussian filter are given for each axis as a
292 sequence, or as a single number, in which case it is equal for
293 all axes.
294 order : int or sequence of ints, optional
295 The order of the filter along each axis is given as a sequence
296 of integers, or as a single number. An order of 0 corresponds
297 to convolution with a Gaussian kernel. A positive order
298 corresponds to convolution with that derivative of a Gaussian.
299 %(output)s
300 %(mode_multiple)s
301 %(cval)s
302 truncate : float, optional
303 Truncate the filter at this many standard deviations.
304 Default is 4.0.
305 radius : None or int or sequence of ints, optional
306 Radius of the Gaussian kernel. The radius are given for each axis
307 as a sequence, or as a single number, in which case it is equal
308 for all axes. If specified, the size of the kernel along each axis
309 will be ``2*radius + 1``, and `truncate` is ignored.
310 Default is None.
311 axes : tuple of int or None, optional
312 If None, `input` is filtered along all axes. Otherwise,
313 `input` is filtered along the specified axes. When `axes` is
314 specified, any tuples used for `sigma`, `order`, `mode` and/or `radius`
315 must match the length of `axes`. The ith entry in any of these tuples
316 corresponds to the ith entry in `axes`.
318 Returns
319 -------
320 gaussian_filter : ndarray
321 Returned array of same shape as `input`.
323 Notes
324 -----
325 The multidimensional filter is implemented as a sequence of
326 1-D convolution filters. The intermediate arrays are
327 stored in the same data type as the output. Therefore, for output
328 types with a limited precision, the results may be imprecise
329 because intermediate results may be stored with insufficient
330 precision.
332 The Gaussian kernel will have size ``2*radius + 1`` along each axis. If
333 `radius` is None, the default ``radius = round(truncate * sigma)`` will be
334 used.
336 Examples
337 --------
338 >>> from scipy.ndimage import gaussian_filter
339 >>> import numpy as np
340 >>> a = np.arange(50, step=2).reshape((5,5))
341 >>> a
342 array([[ 0, 2, 4, 6, 8],
343 [10, 12, 14, 16, 18],
344 [20, 22, 24, 26, 28],
345 [30, 32, 34, 36, 38],
346 [40, 42, 44, 46, 48]])
347 >>> gaussian_filter(a, sigma=1)
348 array([[ 4, 6, 8, 9, 11],
349 [10, 12, 14, 15, 17],
350 [20, 22, 24, 25, 27],
351 [29, 31, 33, 34, 36],
352 [35, 37, 39, 40, 42]])
354 >>> from scipy import datasets
355 >>> import matplotlib.pyplot as plt
356 >>> fig = plt.figure()
357 >>> plt.gray() # show the filtered result in grayscale
358 >>> ax1 = fig.add_subplot(121) # left side
359 >>> ax2 = fig.add_subplot(122) # right side
360 >>> ascent = datasets.ascent()
361 >>> result = gaussian_filter(ascent, sigma=5)
362 >>> ax1.imshow(ascent)
363 >>> ax2.imshow(result)
364 >>> plt.show()
365 """
366 input = numpy.asarray(input)
367 output = _ni_support._get_output(output, input)
369 axes = _ni_support._check_axes(axes, input.ndim)
370 num_axes = len(axes)
371 orders = _ni_support._normalize_sequence(order, num_axes)
372 sigmas = _ni_support._normalize_sequence(sigma, num_axes)
373 modes = _ni_support._normalize_sequence(mode, num_axes)
374 radiuses = _ni_support._normalize_sequence(radius, num_axes)
375 axes = [(axes[ii], sigmas[ii], orders[ii], modes[ii], radiuses[ii])
376 for ii in range(num_axes) if sigmas[ii] > 1e-15]
377 if len(axes) > 0:
378 for axis, sigma, order, mode, radius in axes:
379 gaussian_filter1d(input, sigma, axis, order, output,
380 mode, cval, truncate, radius=radius)
381 input = output
382 else:
383 output[...] = input[...]
384 return output
387@_ni_docstrings.docfiller
388def prewitt(input, axis=-1, output=None, mode="reflect", cval=0.0):
389 """Calculate a Prewitt filter.
391 Parameters
392 ----------
393 %(input)s
394 %(axis)s
395 %(output)s
396 %(mode_multiple)s
397 %(cval)s
399 Examples
400 --------
401 >>> from scipy import ndimage, datasets
402 >>> import matplotlib.pyplot as plt
403 >>> fig = plt.figure()
404 >>> plt.gray() # show the filtered result in grayscale
405 >>> ax1 = fig.add_subplot(121) # left side
406 >>> ax2 = fig.add_subplot(122) # right side
407 >>> ascent = datasets.ascent()
408 >>> result = ndimage.prewitt(ascent)
409 >>> ax1.imshow(ascent)
410 >>> ax2.imshow(result)
411 >>> plt.show()
412 """
413 input = numpy.asarray(input)
414 axis = normalize_axis_index(axis, input.ndim)
415 output = _ni_support._get_output(output, input)
416 modes = _ni_support._normalize_sequence(mode, input.ndim)
417 correlate1d(input, [-1, 0, 1], axis, output, modes[axis], cval, 0)
418 axes = [ii for ii in range(input.ndim) if ii != axis]
419 for ii in axes:
420 correlate1d(output, [1, 1, 1], ii, output, modes[ii], cval, 0,)
421 return output
424@_ni_docstrings.docfiller
425def sobel(input, axis=-1, output=None, mode="reflect", cval=0.0):
426 """Calculate a Sobel filter.
428 Parameters
429 ----------
430 %(input)s
431 %(axis)s
432 %(output)s
433 %(mode_multiple)s
434 %(cval)s
436 Notes
437 -----
438 This function computes the axis-specific Sobel gradient.
439 The horizontal edges can emphasised with the horizontal trasform (axis=0),
440 the vertical edges with the vertical transform (axis=1) and so on for higher
441 dimensions. These can be combined to give the magnitude.
443 Examples
444 --------
445 >>> from scipy import ndimage, datasets
446 >>> import matplotlib.pyplot as plt
447 >>> import numpy as np
448 >>> ascent = datasets.ascent().astype('int32')
449 >>> sobel_h = ndimage.sobel(ascent, 0) # horizontal gradient
450 >>> sobel_v = ndimage.sobel(ascent, 1) # vertical gradient
451 >>> magnitude = np.sqrt(sobel_h**2 + sobel_v**2)
452 >>> magnitude *= 255.0 / np.max(magnitude) # normalization
453 >>> fig, axs = plt.subplots(2, 2, figsize=(8, 8))
454 >>> plt.gray() # show the filtered result in grayscale
455 >>> axs[0, 0].imshow(ascent)
456 >>> axs[0, 1].imshow(sobel_h)
457 >>> axs[1, 0].imshow(sobel_v)
458 >>> axs[1, 1].imshow(magnitude)
459 >>> titles = ["original", "horizontal", "vertical", "magnitude"]
460 >>> for i, ax in enumerate(axs.ravel()):
461 ... ax.set_title(titles[i])
462 ... ax.axis("off")
463 >>> plt.show()
465 """
466 input = numpy.asarray(input)
467 axis = normalize_axis_index(axis, input.ndim)
468 output = _ni_support._get_output(output, input)
469 modes = _ni_support._normalize_sequence(mode, input.ndim)
470 correlate1d(input, [-1, 0, 1], axis, output, modes[axis], cval, 0)
471 axes = [ii for ii in range(input.ndim) if ii != axis]
472 for ii in axes:
473 correlate1d(output, [1, 2, 1], ii, output, modes[ii], cval, 0)
474 return output
477@_ni_docstrings.docfiller
478def generic_laplace(input, derivative2, output=None, mode="reflect",
479 cval=0.0,
480 extra_arguments=(),
481 extra_keywords=None):
482 """
483 N-D Laplace filter using a provided second derivative function.
485 Parameters
486 ----------
487 %(input)s
488 derivative2 : callable
489 Callable with the following signature::
491 derivative2(input, axis, output, mode, cval,
492 *extra_arguments, **extra_keywords)
494 See `extra_arguments`, `extra_keywords` below.
495 %(output)s
496 %(mode_multiple)s
497 %(cval)s
498 %(extra_keywords)s
499 %(extra_arguments)s
500 """
501 if extra_keywords is None:
502 extra_keywords = {}
503 input = numpy.asarray(input)
504 output = _ni_support._get_output(output, input)
505 axes = list(range(input.ndim))
506 if len(axes) > 0:
507 modes = _ni_support._normalize_sequence(mode, len(axes))
508 derivative2(input, axes[0], output, modes[0], cval,
509 *extra_arguments, **extra_keywords)
510 for ii in range(1, len(axes)):
511 tmp = derivative2(input, axes[ii], output.dtype, modes[ii], cval,
512 *extra_arguments, **extra_keywords)
513 output += tmp
514 else:
515 output[...] = input[...]
516 return output
519@_ni_docstrings.docfiller
520def laplace(input, output=None, mode="reflect", cval=0.0):
521 """N-D Laplace filter based on approximate second derivatives.
523 Parameters
524 ----------
525 %(input)s
526 %(output)s
527 %(mode_multiple)s
528 %(cval)s
530 Examples
531 --------
532 >>> from scipy import ndimage, datasets
533 >>> import matplotlib.pyplot as plt
534 >>> fig = plt.figure()
535 >>> plt.gray() # show the filtered result in grayscale
536 >>> ax1 = fig.add_subplot(121) # left side
537 >>> ax2 = fig.add_subplot(122) # right side
538 >>> ascent = datasets.ascent()
539 >>> result = ndimage.laplace(ascent)
540 >>> ax1.imshow(ascent)
541 >>> ax2.imshow(result)
542 >>> plt.show()
543 """
544 def derivative2(input, axis, output, mode, cval):
545 return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
546 return generic_laplace(input, derivative2, output, mode, cval)
549@_ni_docstrings.docfiller
550def gaussian_laplace(input, sigma, output=None, mode="reflect",
551 cval=0.0, **kwargs):
552 """Multidimensional Laplace filter using Gaussian second derivatives.
554 Parameters
555 ----------
556 %(input)s
557 sigma : scalar or sequence of scalars
558 The standard deviations of the Gaussian filter are given for
559 each axis as a sequence, or as a single number, in which case
560 it is equal for all axes.
561 %(output)s
562 %(mode_multiple)s
563 %(cval)s
564 Extra keyword arguments will be passed to gaussian_filter().
566 Examples
567 --------
568 >>> from scipy import ndimage, datasets
569 >>> import matplotlib.pyplot as plt
570 >>> ascent = datasets.ascent()
572 >>> fig = plt.figure()
573 >>> plt.gray() # show the filtered result in grayscale
574 >>> ax1 = fig.add_subplot(121) # left side
575 >>> ax2 = fig.add_subplot(122) # right side
577 >>> result = ndimage.gaussian_laplace(ascent, sigma=1)
578 >>> ax1.imshow(result)
580 >>> result = ndimage.gaussian_laplace(ascent, sigma=3)
581 >>> ax2.imshow(result)
582 >>> plt.show()
583 """
584 input = numpy.asarray(input)
586 def derivative2(input, axis, output, mode, cval, sigma, **kwargs):
587 order = [0] * input.ndim
588 order[axis] = 2
589 return gaussian_filter(input, sigma, order, output, mode, cval,
590 **kwargs)
592 return generic_laplace(input, derivative2, output, mode, cval,
593 extra_arguments=(sigma,),
594 extra_keywords=kwargs)
597@_ni_docstrings.docfiller
598def generic_gradient_magnitude(input, derivative, output=None,
599 mode="reflect", cval=0.0,
600 extra_arguments=(), extra_keywords=None):
601 """Gradient magnitude using a provided gradient function.
603 Parameters
604 ----------
605 %(input)s
606 derivative : callable
607 Callable with the following signature::
609 derivative(input, axis, output, mode, cval,
610 *extra_arguments, **extra_keywords)
612 See `extra_arguments`, `extra_keywords` below.
613 `derivative` can assume that `input` and `output` are ndarrays.
614 Note that the output from `derivative` is modified inplace;
615 be careful to copy important inputs before returning them.
616 %(output)s
617 %(mode_multiple)s
618 %(cval)s
619 %(extra_keywords)s
620 %(extra_arguments)s
621 """
622 if extra_keywords is None:
623 extra_keywords = {}
624 input = numpy.asarray(input)
625 output = _ni_support._get_output(output, input)
626 axes = list(range(input.ndim))
627 if len(axes) > 0:
628 modes = _ni_support._normalize_sequence(mode, len(axes))
629 derivative(input, axes[0], output, modes[0], cval,
630 *extra_arguments, **extra_keywords)
631 numpy.multiply(output, output, output)
632 for ii in range(1, len(axes)):
633 tmp = derivative(input, axes[ii], output.dtype, modes[ii], cval,
634 *extra_arguments, **extra_keywords)
635 numpy.multiply(tmp, tmp, tmp)
636 output += tmp
637 # This allows the sqrt to work with a different default casting
638 numpy.sqrt(output, output, casting='unsafe')
639 else:
640 output[...] = input[...]
641 return output
644@_ni_docstrings.docfiller
645def gaussian_gradient_magnitude(input, sigma, output=None,
646 mode="reflect", cval=0.0, **kwargs):
647 """Multidimensional gradient magnitude using Gaussian derivatives.
649 Parameters
650 ----------
651 %(input)s
652 sigma : scalar or sequence of scalars
653 The standard deviations of the Gaussian filter are given for
654 each axis as a sequence, or as a single number, in which case
655 it is equal for all axes.
656 %(output)s
657 %(mode_multiple)s
658 %(cval)s
659 Extra keyword arguments will be passed to gaussian_filter().
661 Returns
662 -------
663 gaussian_gradient_magnitude : ndarray
664 Filtered array. Has the same shape as `input`.
666 Examples
667 --------
668 >>> from scipy import ndimage, datasets
669 >>> import matplotlib.pyplot as plt
670 >>> fig = plt.figure()
671 >>> plt.gray() # show the filtered result in grayscale
672 >>> ax1 = fig.add_subplot(121) # left side
673 >>> ax2 = fig.add_subplot(122) # right side
674 >>> ascent = datasets.ascent()
675 >>> result = ndimage.gaussian_gradient_magnitude(ascent, sigma=5)
676 >>> ax1.imshow(ascent)
677 >>> ax2.imshow(result)
678 >>> plt.show()
679 """
680 input = numpy.asarray(input)
682 def derivative(input, axis, output, mode, cval, sigma, **kwargs):
683 order = [0] * input.ndim
684 order[axis] = 1
685 return gaussian_filter(input, sigma, order, output, mode,
686 cval, **kwargs)
688 return generic_gradient_magnitude(input, derivative, output, mode,
689 cval, extra_arguments=(sigma,),
690 extra_keywords=kwargs)
693def _correlate_or_convolve(input, weights, output, mode, cval, origin,
694 convolution):
695 input = numpy.asarray(input)
696 weights = numpy.asarray(weights)
697 complex_input = input.dtype.kind == 'c'
698 complex_weights = weights.dtype.kind == 'c'
699 if complex_input or complex_weights:
700 if complex_weights and not convolution:
701 # As for numpy.correlate, conjugate weights rather than input.
702 weights = weights.conj()
703 kwargs = dict(
704 mode=mode, origin=origin, convolution=convolution
705 )
706 output = _ni_support._get_output(output, input, complex_output=True)
708 return _complex_via_real_components(_correlate_or_convolve, input,
709 weights, output, cval, **kwargs)
711 origins = _ni_support._normalize_sequence(origin, input.ndim)
712 weights = numpy.asarray(weights, dtype=numpy.float64)
713 wshape = [ii for ii in weights.shape if ii > 0]
714 if len(wshape) != input.ndim:
715 raise RuntimeError('filter weights array has incorrect shape.')
716 if convolution:
717 weights = weights[tuple([slice(None, None, -1)] * weights.ndim)]
718 for ii in range(len(origins)):
719 origins[ii] = -origins[ii]
720 if not weights.shape[ii] & 1:
721 origins[ii] -= 1
722 for origin, lenw in zip(origins, wshape):
723 if _invalid_origin(origin, lenw):
724 raise ValueError('Invalid origin; origin must satisfy '
725 '-(weights.shape[k] // 2) <= origin[k] <= '
726 '(weights.shape[k]-1) // 2')
728 if not weights.flags.contiguous:
729 weights = weights.copy()
730 output = _ni_support._get_output(output, input)
731 temp_needed = numpy.may_share_memory(input, output)
732 if temp_needed:
733 # input and output arrays cannot share memory
734 temp = output
735 output = _ni_support._get_output(output.dtype, input)
736 if not isinstance(mode, str) and isinstance(mode, Iterable):
737 raise RuntimeError("A sequence of modes is not supported")
738 mode = _ni_support._extend_mode_to_code(mode)
739 _nd_image.correlate(input, weights, output, mode, cval, origins)
740 if temp_needed:
741 temp[...] = output
742 output = temp
743 return output
746@_ni_docstrings.docfiller
747def correlate(input, weights, output=None, mode='reflect', cval=0.0,
748 origin=0):
749 """
750 Multidimensional correlation.
752 The array is correlated with the given kernel.
754 Parameters
755 ----------
756 %(input)s
757 weights : ndarray
758 array of weights, same number of dimensions as input
759 %(output)s
760 %(mode_reflect)s
761 %(cval)s
762 %(origin_multiple)s
764 Returns
765 -------
766 result : ndarray
767 The result of correlation of `input` with `weights`.
769 See Also
770 --------
771 convolve : Convolve an image with a kernel.
773 Examples
774 --------
775 Correlation is the process of moving a filter mask often referred to
776 as kernel over the image and computing the sum of products at each location.
778 >>> from scipy.ndimage import correlate
779 >>> import numpy as np
780 >>> input_img = np.arange(25).reshape(5,5)
781 >>> print(input_img)
782 [[ 0 1 2 3 4]
783 [ 5 6 7 8 9]
784 [10 11 12 13 14]
785 [15 16 17 18 19]
786 [20 21 22 23 24]]
788 Define a kernel (weights) for correlation. In this example, it is for sum of
789 center and up, down, left and right next elements.
791 >>> weights = [[0, 1, 0],
792 ... [1, 1, 1],
793 ... [0, 1, 0]]
795 We can calculate a correlation result:
796 For example, element ``[2,2]`` is ``7 + 11 + 12 + 13 + 17 = 60``.
798 >>> correlate(input_img, weights)
799 array([[ 6, 10, 15, 20, 24],
800 [ 26, 30, 35, 40, 44],
801 [ 51, 55, 60, 65, 69],
802 [ 76, 80, 85, 90, 94],
803 [ 96, 100, 105, 110, 114]])
805 """
806 return _correlate_or_convolve(input, weights, output, mode, cval,
807 origin, False)
810@_ni_docstrings.docfiller
811def convolve(input, weights, output=None, mode='reflect', cval=0.0,
812 origin=0):
813 """
814 Multidimensional convolution.
816 The array is convolved with the given kernel.
818 Parameters
819 ----------
820 %(input)s
821 weights : array_like
822 Array of weights, same number of dimensions as input
823 %(output)s
824 %(mode_reflect)s
825 cval : scalar, optional
826 Value to fill past edges of input if `mode` is 'constant'. Default
827 is 0.0
828 origin : int, optional
829 Controls the origin of the input signal, which is where the
830 filter is centered to produce the first element of the output.
831 Positive values shift the filter to the right, and negative values
832 shift the filter to the left. Default is 0.
834 Returns
835 -------
836 result : ndarray
837 The result of convolution of `input` with `weights`.
839 See Also
840 --------
841 correlate : Correlate an image with a kernel.
843 Notes
844 -----
845 Each value in result is :math:`C_i = \\sum_j{I_{i+k-j} W_j}`, where
846 W is the `weights` kernel,
847 j is the N-D spatial index over :math:`W`,
848 I is the `input` and k is the coordinate of the center of
849 W, specified by `origin` in the input parameters.
851 Examples
852 --------
853 Perhaps the simplest case to understand is ``mode='constant', cval=0.0``,
854 because in this case borders (i.e., where the `weights` kernel, centered
855 on any one value, extends beyond an edge of `input`) are treated as zeros.
857 >>> import numpy as np
858 >>> a = np.array([[1, 2, 0, 0],
859 ... [5, 3, 0, 4],
860 ... [0, 0, 0, 7],
861 ... [9, 3, 0, 0]])
862 >>> k = np.array([[1,1,1],[1,1,0],[1,0,0]])
863 >>> from scipy import ndimage
864 >>> ndimage.convolve(a, k, mode='constant', cval=0.0)
865 array([[11, 10, 7, 4],
866 [10, 3, 11, 11],
867 [15, 12, 14, 7],
868 [12, 3, 7, 0]])
870 Setting ``cval=1.0`` is equivalent to padding the outer edge of `input`
871 with 1.0's (and then extracting only the original region of the result).
873 >>> ndimage.convolve(a, k, mode='constant', cval=1.0)
874 array([[13, 11, 8, 7],
875 [11, 3, 11, 14],
876 [16, 12, 14, 10],
877 [15, 6, 10, 5]])
879 With ``mode='reflect'`` (the default), outer values are reflected at the
880 edge of `input` to fill in missing values.
882 >>> b = np.array([[2, 0, 0],
883 ... [1, 0, 0],
884 ... [0, 0, 0]])
885 >>> k = np.array([[0,1,0], [0,1,0], [0,1,0]])
886 >>> ndimage.convolve(b, k, mode='reflect')
887 array([[5, 0, 0],
888 [3, 0, 0],
889 [1, 0, 0]])
891 This includes diagonally at the corners.
893 >>> k = np.array([[1,0,0],[0,1,0],[0,0,1]])
894 >>> ndimage.convolve(b, k)
895 array([[4, 2, 0],
896 [3, 2, 0],
897 [1, 1, 0]])
899 With ``mode='nearest'``, the single nearest value in to an edge in
900 `input` is repeated as many times as needed to match the overlapping
901 `weights`.
903 >>> c = np.array([[2, 0, 1],
904 ... [1, 0, 0],
905 ... [0, 0, 0]])
906 >>> k = np.array([[0, 1, 0],
907 ... [0, 1, 0],
908 ... [0, 1, 0],
909 ... [0, 1, 0],
910 ... [0, 1, 0]])
911 >>> ndimage.convolve(c, k, mode='nearest')
912 array([[7, 0, 3],
913 [5, 0, 2],
914 [3, 0, 1]])
916 """
917 return _correlate_or_convolve(input, weights, output, mode, cval,
918 origin, True)
921@_ni_docstrings.docfiller
922def uniform_filter1d(input, size, axis=-1, output=None,
923 mode="reflect", cval=0.0, origin=0):
924 """Calculate a 1-D uniform filter along the given axis.
926 The lines of the array along the given axis are filtered with a
927 uniform filter of given size.
929 Parameters
930 ----------
931 %(input)s
932 size : int
933 length of uniform filter
934 %(axis)s
935 %(output)s
936 %(mode_reflect)s
937 %(cval)s
938 %(origin)s
940 Examples
941 --------
942 >>> from scipy.ndimage import uniform_filter1d
943 >>> uniform_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3)
944 array([4, 3, 4, 1, 4, 6, 6, 3])
945 """
946 input = numpy.asarray(input)
947 axis = normalize_axis_index(axis, input.ndim)
948 if size < 1:
949 raise RuntimeError('incorrect filter size')
950 complex_output = input.dtype.kind == 'c'
951 output = _ni_support._get_output(output, input,
952 complex_output=complex_output)
953 if (size // 2 + origin < 0) or (size // 2 + origin >= size):
954 raise ValueError('invalid origin')
955 mode = _ni_support._extend_mode_to_code(mode)
956 if not complex_output:
957 _nd_image.uniform_filter1d(input, size, axis, output, mode, cval,
958 origin)
959 else:
960 _nd_image.uniform_filter1d(input.real, size, axis, output.real, mode,
961 numpy.real(cval), origin)
962 _nd_image.uniform_filter1d(input.imag, size, axis, output.imag, mode,
963 numpy.imag(cval), origin)
964 return output
967@_ni_docstrings.docfiller
968def uniform_filter(input, size=3, output=None, mode="reflect",
969 cval=0.0, origin=0, *, axes=None):
970 """Multidimensional uniform filter.
972 Parameters
973 ----------
974 %(input)s
975 size : int or sequence of ints, optional
976 The sizes of the uniform filter are given for each axis as a
977 sequence, or as a single number, in which case the size is
978 equal for all axes.
979 %(output)s
980 %(mode_multiple)s
981 %(cval)s
982 %(origin_multiple)s
983 axes : tuple of int or None, optional
984 If None, `input` is filtered along all axes. Otherwise,
985 `input` is filtered along the specified axes. When `axes` is
986 specified, any tuples used for `size`, `origin`, and/or `mode`
987 must match the length of `axes`. The ith entry in any of these tuples
988 corresponds to the ith entry in `axes`.
990 Returns
991 -------
992 uniform_filter : ndarray
993 Filtered array. Has the same shape as `input`.
995 Notes
996 -----
997 The multidimensional filter is implemented as a sequence of
998 1-D uniform filters. The intermediate arrays are stored
999 in the same data type as the output. Therefore, for output types
1000 with a limited precision, the results may be imprecise because
1001 intermediate results may be stored with insufficient precision.
1003 Examples
1004 --------
1005 >>> from scipy import ndimage, datasets
1006 >>> import matplotlib.pyplot as plt
1007 >>> fig = plt.figure()
1008 >>> plt.gray() # show the filtered result in grayscale
1009 >>> ax1 = fig.add_subplot(121) # left side
1010 >>> ax2 = fig.add_subplot(122) # right side
1011 >>> ascent = datasets.ascent()
1012 >>> result = ndimage.uniform_filter(ascent, size=20)
1013 >>> ax1.imshow(ascent)
1014 >>> ax2.imshow(result)
1015 >>> plt.show()
1016 """
1017 input = numpy.asarray(input)
1018 output = _ni_support._get_output(output, input,
1019 complex_output=input.dtype.kind == 'c')
1020 axes = _ni_support._check_axes(axes, input.ndim)
1021 num_axes = len(axes)
1022 sizes = _ni_support._normalize_sequence(size, num_axes)
1023 origins = _ni_support._normalize_sequence(origin, num_axes)
1024 modes = _ni_support._normalize_sequence(mode, num_axes)
1025 axes = [(axes[ii], sizes[ii], origins[ii], modes[ii])
1026 for ii in range(num_axes) if sizes[ii] > 1]
1027 if len(axes) > 0:
1028 for axis, size, origin, mode in axes:
1029 uniform_filter1d(input, int(size), axis, output, mode,
1030 cval, origin)
1031 input = output
1032 else:
1033 output[...] = input[...]
1034 return output
1037@_ni_docstrings.docfiller
1038def minimum_filter1d(input, size, axis=-1, output=None,
1039 mode="reflect", cval=0.0, origin=0):
1040 """Calculate a 1-D minimum filter along the given axis.
1042 The lines of the array along the given axis are filtered with a
1043 minimum filter of given size.
1045 Parameters
1046 ----------
1047 %(input)s
1048 size : int
1049 length along which to calculate 1D minimum
1050 %(axis)s
1051 %(output)s
1052 %(mode_reflect)s
1053 %(cval)s
1054 %(origin)s
1056 Notes
1057 -----
1058 This function implements the MINLIST algorithm [1]_, as described by
1059 Richard Harter [2]_, and has a guaranteed O(n) performance, `n` being
1060 the `input` length, regardless of filter size.
1062 References
1063 ----------
1064 .. [1] http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.2777
1065 .. [2] http://www.richardhartersworld.com/cri/2001/slidingmin.html
1068 Examples
1069 --------
1070 >>> from scipy.ndimage import minimum_filter1d
1071 >>> minimum_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3)
1072 array([2, 0, 0, 0, 1, 1, 0, 0])
1073 """
1074 input = numpy.asarray(input)
1075 if numpy.iscomplexobj(input):
1076 raise TypeError('Complex type not supported')
1077 axis = normalize_axis_index(axis, input.ndim)
1078 if size < 1:
1079 raise RuntimeError('incorrect filter size')
1080 output = _ni_support._get_output(output, input)
1081 if (size // 2 + origin < 0) or (size // 2 + origin >= size):
1082 raise ValueError('invalid origin')
1083 mode = _ni_support._extend_mode_to_code(mode)
1084 _nd_image.min_or_max_filter1d(input, size, axis, output, mode, cval,
1085 origin, 1)
1086 return output
1089@_ni_docstrings.docfiller
1090def maximum_filter1d(input, size, axis=-1, output=None,
1091 mode="reflect", cval=0.0, origin=0):
1092 """Calculate a 1-D maximum filter along the given axis.
1094 The lines of the array along the given axis are filtered with a
1095 maximum filter of given size.
1097 Parameters
1098 ----------
1099 %(input)s
1100 size : int
1101 Length along which to calculate the 1-D maximum.
1102 %(axis)s
1103 %(output)s
1104 %(mode_reflect)s
1105 %(cval)s
1106 %(origin)s
1108 Returns
1109 -------
1110 maximum1d : ndarray, None
1111 Maximum-filtered array with same shape as input.
1112 None if `output` is not None
1114 Notes
1115 -----
1116 This function implements the MAXLIST algorithm [1]_, as described by
1117 Richard Harter [2]_, and has a guaranteed O(n) performance, `n` being
1118 the `input` length, regardless of filter size.
1120 References
1121 ----------
1122 .. [1] http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.2777
1123 .. [2] http://www.richardhartersworld.com/cri/2001/slidingmin.html
1125 Examples
1126 --------
1127 >>> from scipy.ndimage import maximum_filter1d
1128 >>> maximum_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3)
1129 array([8, 8, 8, 4, 9, 9, 9, 9])
1130 """
1131 input = numpy.asarray(input)
1132 if numpy.iscomplexobj(input):
1133 raise TypeError('Complex type not supported')
1134 axis = normalize_axis_index(axis, input.ndim)
1135 if size < 1:
1136 raise RuntimeError('incorrect filter size')
1137 output = _ni_support._get_output(output, input)
1138 if (size // 2 + origin < 0) or (size // 2 + origin >= size):
1139 raise ValueError('invalid origin')
1140 mode = _ni_support._extend_mode_to_code(mode)
1141 _nd_image.min_or_max_filter1d(input, size, axis, output, mode, cval,
1142 origin, 0)
1143 return output
1146def _min_or_max_filter(input, size, footprint, structure, output, mode,
1147 cval, origin, minimum, axes=None):
1148 if (size is not None) and (footprint is not None):
1149 warnings.warn("ignoring size because footprint is set", UserWarning, stacklevel=3)
1150 if structure is None:
1151 if footprint is None:
1152 if size is None:
1153 raise RuntimeError("no footprint provided")
1154 separable = True
1155 else:
1156 footprint = numpy.asarray(footprint, dtype=bool)
1157 if not footprint.any():
1158 raise ValueError("All-zero footprint is not supported.")
1159 if footprint.all():
1160 size = footprint.shape
1161 footprint = None
1162 separable = True
1163 else:
1164 separable = False
1165 else:
1166 structure = numpy.asarray(structure, dtype=numpy.float64)
1167 separable = False
1168 if footprint is None:
1169 footprint = numpy.ones(structure.shape, bool)
1170 else:
1171 footprint = numpy.asarray(footprint, dtype=bool)
1172 input = numpy.asarray(input)
1173 if numpy.iscomplexobj(input):
1174 raise TypeError('Complex type not supported')
1175 output = _ni_support._get_output(output, input)
1176 temp_needed = numpy.may_share_memory(input, output)
1177 if temp_needed:
1178 # input and output arrays cannot share memory
1179 temp = output
1180 output = _ni_support._get_output(output.dtype, input)
1181 axes = _ni_support._check_axes(axes, input.ndim)
1182 num_axes = len(axes)
1183 if separable:
1184 origins = _ni_support._normalize_sequence(origin, num_axes)
1185 sizes = _ni_support._normalize_sequence(size, num_axes)
1186 modes = _ni_support._normalize_sequence(mode, num_axes)
1187 axes = [(axes[ii], sizes[ii], origins[ii], modes[ii])
1188 for ii in range(len(axes)) if sizes[ii] > 1]
1189 if minimum:
1190 filter_ = minimum_filter1d
1191 else:
1192 filter_ = maximum_filter1d
1193 if len(axes) > 0:
1194 for axis, size, origin, mode in axes:
1195 filter_(input, int(size), axis, output, mode, cval, origin)
1196 input = output
1197 else:
1198 output[...] = input[...]
1199 else:
1200 origins = _ni_support._normalize_sequence(origin, input.ndim)
1201 if num_axes < input.ndim:
1202 if footprint.ndim != num_axes:
1203 raise RuntimeError("footprint array has incorrect shape")
1204 footprint = numpy.expand_dims(
1205 footprint,
1206 tuple(ax for ax in range(input.ndim) if ax not in axes)
1207 )
1208 fshape = [ii for ii in footprint.shape if ii > 0]
1209 if len(fshape) != input.ndim:
1210 raise RuntimeError('footprint array has incorrect shape.')
1211 for origin, lenf in zip(origins, fshape):
1212 if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf):
1213 raise ValueError('invalid origin')
1214 if not footprint.flags.contiguous:
1215 footprint = footprint.copy()
1216 if structure is not None:
1217 if len(structure.shape) != input.ndim:
1218 raise RuntimeError('structure array has incorrect shape')
1219 if num_axes != structure.ndim:
1220 structure = numpy.expand_dims(
1221 structure,
1222 tuple(ax for ax in range(structure.ndim) if ax not in axes)
1223 )
1224 if not structure.flags.contiguous:
1225 structure = structure.copy()
1226 if not isinstance(mode, str) and isinstance(mode, Iterable):
1227 raise RuntimeError(
1228 "A sequence of modes is not supported for non-separable "
1229 "footprints")
1230 mode = _ni_support._extend_mode_to_code(mode)
1231 _nd_image.min_or_max_filter(input, footprint, structure, output,
1232 mode, cval, origins, minimum)
1233 if temp_needed:
1234 temp[...] = output
1235 output = temp
1236 return output
1239@_ni_docstrings.docfiller
1240def minimum_filter(input, size=None, footprint=None, output=None,
1241 mode="reflect", cval=0.0, origin=0, *, axes=None):
1242 """Calculate a multidimensional minimum filter.
1244 Parameters
1245 ----------
1246 %(input)s
1247 %(size_foot)s
1248 %(output)s
1249 %(mode_multiple)s
1250 %(cval)s
1251 %(origin_multiple)s
1252 axes : tuple of int or None, optional
1253 If None, `input` is filtered along all axes. Otherwise,
1254 `input` is filtered along the specified axes. When `axes` is
1255 specified, any tuples used for `size`, `origin`, and/or `mode`
1256 must match the length of `axes`. The ith entry in any of these tuples
1257 corresponds to the ith entry in `axes`.
1259 Returns
1260 -------
1261 minimum_filter : ndarray
1262 Filtered array. Has the same shape as `input`.
1264 Notes
1265 -----
1266 A sequence of modes (one per axis) is only supported when the footprint is
1267 separable. Otherwise, a single mode string must be provided.
1269 Examples
1270 --------
1271 >>> from scipy import ndimage, datasets
1272 >>> import matplotlib.pyplot as plt
1273 >>> fig = plt.figure()
1274 >>> plt.gray() # show the filtered result in grayscale
1275 >>> ax1 = fig.add_subplot(121) # left side
1276 >>> ax2 = fig.add_subplot(122) # right side
1277 >>> ascent = datasets.ascent()
1278 >>> result = ndimage.minimum_filter(ascent, size=20)
1279 >>> ax1.imshow(ascent)
1280 >>> ax2.imshow(result)
1281 >>> plt.show()
1282 """
1283 return _min_or_max_filter(input, size, footprint, None, output, mode,
1284 cval, origin, 1, axes)
1287@_ni_docstrings.docfiller
1288def maximum_filter(input, size=None, footprint=None, output=None,
1289 mode="reflect", cval=0.0, origin=0, *, axes=None):
1290 """Calculate a multidimensional maximum filter.
1292 Parameters
1293 ----------
1294 %(input)s
1295 %(size_foot)s
1296 %(output)s
1297 %(mode_multiple)s
1298 %(cval)s
1299 %(origin_multiple)s
1300 axes : tuple of int or None, optional
1301 If None, `input` is filtered along all axes. Otherwise,
1302 `input` is filtered along the specified axes. When `axes` is
1303 specified, any tuples used for `size`, `origin`, and/or `mode`
1304 must match the length of `axes`. The ith entry in any of these tuples
1305 corresponds to the ith entry in `axes`.
1307 Returns
1308 -------
1309 maximum_filter : ndarray
1310 Filtered array. Has the same shape as `input`.
1312 Notes
1313 -----
1314 A sequence of modes (one per axis) is only supported when the footprint is
1315 separable. Otherwise, a single mode string must be provided.
1317 Examples
1318 --------
1319 >>> from scipy import ndimage, datasets
1320 >>> import matplotlib.pyplot as plt
1321 >>> fig = plt.figure()
1322 >>> plt.gray() # show the filtered result in grayscale
1323 >>> ax1 = fig.add_subplot(121) # left side
1324 >>> ax2 = fig.add_subplot(122) # right side
1325 >>> ascent = datasets.ascent()
1326 >>> result = ndimage.maximum_filter(ascent, size=20)
1327 >>> ax1.imshow(ascent)
1328 >>> ax2.imshow(result)
1329 >>> plt.show()
1330 """
1331 return _min_or_max_filter(input, size, footprint, None, output, mode,
1332 cval, origin, 0, axes)
1335@_ni_docstrings.docfiller
1336def _rank_filter(input, rank, size=None, footprint=None, output=None,
1337 mode="reflect", cval=0.0, origin=0, operation='rank',
1338 axes=None):
1339 if (size is not None) and (footprint is not None):
1340 warnings.warn("ignoring size because footprint is set", UserWarning, stacklevel=3)
1341 input = numpy.asarray(input)
1342 if numpy.iscomplexobj(input):
1343 raise TypeError('Complex type not supported')
1344 axes = _ni_support._check_axes(axes, input.ndim)
1345 num_axes = len(axes)
1346 origins = _ni_support._normalize_sequence(origin, num_axes)
1347 if footprint is None:
1348 if size is None:
1349 raise RuntimeError("no footprint or filter size provided")
1350 sizes = _ni_support._normalize_sequence(size, num_axes)
1351 footprint = numpy.ones(sizes, dtype=bool)
1352 else:
1353 footprint = numpy.asarray(footprint, dtype=bool)
1354 if num_axes < input.ndim:
1355 # set origin = 0 for any axes not being filtered
1356 origins_temp = [0,] * input.ndim
1357 for o, ax in zip(origins, axes):
1358 origins_temp[ax] = o
1359 origins = origins_temp
1361 if not isinstance(mode, str) and isinstance(mode, Iterable):
1362 # set mode = 'constant' for any axes not being filtered
1363 modes = _ni_support._normalize_sequence(mode, num_axes)
1364 modes_temp = ['constant'] * input.ndim
1365 for m, ax in zip(modes, axes):
1366 modes_temp[ax] = m
1367 mode = modes_temp
1369 # insert singleton dimension along any non-filtered axes
1370 if footprint.ndim != num_axes:
1371 raise RuntimeError("footprint array has incorrect shape")
1372 footprint = numpy.expand_dims(
1373 footprint,
1374 tuple(ax for ax in range(input.ndim) if ax not in axes)
1375 )
1376 fshape = [ii for ii in footprint.shape if ii > 0]
1377 if len(fshape) != input.ndim:
1378 raise RuntimeError('footprint array has incorrect shape.')
1379 for origin, lenf in zip(origins, fshape):
1380 if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf):
1381 raise ValueError('invalid origin')
1382 if not footprint.flags.contiguous:
1383 footprint = footprint.copy()
1384 filter_size = numpy.where(footprint, 1, 0).sum()
1385 if operation == 'median':
1386 rank = filter_size // 2
1387 elif operation == 'percentile':
1388 percentile = rank
1389 if percentile < 0.0:
1390 percentile += 100.0
1391 if percentile < 0 or percentile > 100:
1392 raise RuntimeError('invalid percentile')
1393 if percentile == 100.0:
1394 rank = filter_size - 1
1395 else:
1396 rank = int(float(filter_size) * percentile / 100.0)
1397 if rank < 0:
1398 rank += filter_size
1399 if rank < 0 or rank >= filter_size:
1400 raise RuntimeError('rank not within filter footprint size')
1401 if rank == 0:
1402 return minimum_filter(input, None, footprint, output, mode, cval,
1403 origins, axes=None)
1404 elif rank == filter_size - 1:
1405 return maximum_filter(input, None, footprint, output, mode, cval,
1406 origins, axes=None)
1407 else:
1408 output = _ni_support._get_output(output, input)
1409 temp_needed = numpy.may_share_memory(input, output)
1410 if temp_needed:
1411 # input and output arrays cannot share memory
1412 temp = output
1413 output = _ni_support._get_output(output.dtype, input)
1414 if not isinstance(mode, str) and isinstance(mode, Iterable):
1415 raise RuntimeError(
1416 "A sequence of modes is not supported by non-separable rank "
1417 "filters")
1418 mode = _ni_support._extend_mode_to_code(mode)
1419 _nd_image.rank_filter(input, rank, footprint, output, mode, cval,
1420 origins)
1421 if temp_needed:
1422 temp[...] = output
1423 output = temp
1424 return output
1427@_ni_docstrings.docfiller
1428def rank_filter(input, rank, size=None, footprint=None, output=None,
1429 mode="reflect", cval=0.0, origin=0, *, axes=None):
1430 """Calculate a multidimensional rank filter.
1432 Parameters
1433 ----------
1434 %(input)s
1435 rank : int
1436 The rank parameter may be less than zero, i.e., rank = -1
1437 indicates the largest element.
1438 %(size_foot)s
1439 %(output)s
1440 %(mode_reflect)s
1441 %(cval)s
1442 %(origin_multiple)s
1443 axes : tuple of int or None, optional
1444 If None, `input` is filtered along all axes. Otherwise,
1445 `input` is filtered along the specified axes.
1447 Returns
1448 -------
1449 rank_filter : ndarray
1450 Filtered array. Has the same shape as `input`.
1452 Examples
1453 --------
1454 >>> from scipy import ndimage, datasets
1455 >>> import matplotlib.pyplot as plt
1456 >>> fig = plt.figure()
1457 >>> plt.gray() # show the filtered result in grayscale
1458 >>> ax1 = fig.add_subplot(121) # left side
1459 >>> ax2 = fig.add_subplot(122) # right side
1460 >>> ascent = datasets.ascent()
1461 >>> result = ndimage.rank_filter(ascent, rank=42, size=20)
1462 >>> ax1.imshow(ascent)
1463 >>> ax2.imshow(result)
1464 >>> plt.show()
1465 """
1466 rank = operator.index(rank)
1467 return _rank_filter(input, rank, size, footprint, output, mode, cval,
1468 origin, 'rank', axes=axes)
1471@_ni_docstrings.docfiller
1472def median_filter(input, size=None, footprint=None, output=None,
1473 mode="reflect", cval=0.0, origin=0, *, axes=None):
1474 """
1475 Calculate a multidimensional median filter.
1477 Parameters
1478 ----------
1479 %(input)s
1480 %(size_foot)s
1481 %(output)s
1482 %(mode_reflect)s
1483 %(cval)s
1484 %(origin_multiple)s
1485 axes : tuple of int or None, optional
1486 If None, `input` is filtered along all axes. Otherwise,
1487 `input` is filtered along the specified axes.
1489 Returns
1490 -------
1491 median_filter : ndarray
1492 Filtered array. Has the same shape as `input`.
1494 See Also
1495 --------
1496 scipy.signal.medfilt2d
1498 Notes
1499 -----
1500 For 2-dimensional images with ``uint8``, ``float32`` or ``float64`` dtypes
1501 the specialised function `scipy.signal.medfilt2d` may be faster. It is
1502 however limited to constant mode with ``cval=0``.
1504 Examples
1505 --------
1506 >>> from scipy import ndimage, datasets
1507 >>> import matplotlib.pyplot as plt
1508 >>> fig = plt.figure()
1509 >>> plt.gray() # show the filtered result in grayscale
1510 >>> ax1 = fig.add_subplot(121) # left side
1511 >>> ax2 = fig.add_subplot(122) # right side
1512 >>> ascent = datasets.ascent()
1513 >>> result = ndimage.median_filter(ascent, size=20)
1514 >>> ax1.imshow(ascent)
1515 >>> ax2.imshow(result)
1516 >>> plt.show()
1517 """
1518 return _rank_filter(input, 0, size, footprint, output, mode, cval,
1519 origin, 'median', axes=axes)
1522@_ni_docstrings.docfiller
1523def percentile_filter(input, percentile, size=None, footprint=None,
1524 output=None, mode="reflect", cval=0.0, origin=0, *,
1525 axes=None):
1526 """Calculate a multidimensional percentile filter.
1528 Parameters
1529 ----------
1530 %(input)s
1531 percentile : scalar
1532 The percentile parameter may be less than zero, i.e.,
1533 percentile = -20 equals percentile = 80
1534 %(size_foot)s
1535 %(output)s
1536 %(mode_reflect)s
1537 %(cval)s
1538 %(origin_multiple)s
1539 axes : tuple of int or None, optional
1540 If None, `input` is filtered along all axes. Otherwise,
1541 `input` is filtered along the specified axes.
1543 Returns
1544 -------
1545 percentile_filter : ndarray
1546 Filtered array. Has the same shape as `input`.
1548 Examples
1549 --------
1550 >>> from scipy import ndimage, datasets
1551 >>> import matplotlib.pyplot as plt
1552 >>> fig = plt.figure()
1553 >>> plt.gray() # show the filtered result in grayscale
1554 >>> ax1 = fig.add_subplot(121) # left side
1555 >>> ax2 = fig.add_subplot(122) # right side
1556 >>> ascent = datasets.ascent()
1557 >>> result = ndimage.percentile_filter(ascent, percentile=20, size=20)
1558 >>> ax1.imshow(ascent)
1559 >>> ax2.imshow(result)
1560 >>> plt.show()
1561 """
1562 return _rank_filter(input, percentile, size, footprint, output, mode,
1563 cval, origin, 'percentile', axes=axes)
1566@_ni_docstrings.docfiller
1567def generic_filter1d(input, function, filter_size, axis=-1,
1568 output=None, mode="reflect", cval=0.0, origin=0,
1569 extra_arguments=(), extra_keywords=None):
1570 """Calculate a 1-D filter along the given axis.
1572 `generic_filter1d` iterates over the lines of the array, calling the
1573 given function at each line. The arguments of the line are the
1574 input line, and the output line. The input and output lines are 1-D
1575 double arrays. The input line is extended appropriately according
1576 to the filter size and origin. The output line must be modified
1577 in-place with the result.
1579 Parameters
1580 ----------
1581 %(input)s
1582 function : {callable, scipy.LowLevelCallable}
1583 Function to apply along given axis.
1584 filter_size : scalar
1585 Length of the filter.
1586 %(axis)s
1587 %(output)s
1588 %(mode_reflect)s
1589 %(cval)s
1590 %(origin)s
1591 %(extra_arguments)s
1592 %(extra_keywords)s
1594 Notes
1595 -----
1596 This function also accepts low-level callback functions with one of
1597 the following signatures and wrapped in `scipy.LowLevelCallable`:
1599 .. code:: c
1601 int function(double *input_line, npy_intp input_length,
1602 double *output_line, npy_intp output_length,
1603 void *user_data)
1604 int function(double *input_line, intptr_t input_length,
1605 double *output_line, intptr_t output_length,
1606 void *user_data)
1608 The calling function iterates over the lines of the input and output
1609 arrays, calling the callback function at each line. The current line
1610 is extended according to the border conditions set by the calling
1611 function, and the result is copied into the array that is passed
1612 through ``input_line``. The length of the input line (after extension)
1613 is passed through ``input_length``. The callback function should apply
1614 the filter and store the result in the array passed through
1615 ``output_line``. The length of the output line is passed through
1616 ``output_length``. ``user_data`` is the data pointer provided
1617 to `scipy.LowLevelCallable` as-is.
1619 The callback function must return an integer error status that is zero
1620 if something went wrong and one otherwise. If an error occurs, you should
1621 normally set the python error status with an informative message
1622 before returning, otherwise a default error message is set by the
1623 calling function.
1625 In addition, some other low-level function pointer specifications
1626 are accepted, but these are for backward compatibility only and should
1627 not be used in new code.
1629 """
1630 if extra_keywords is None:
1631 extra_keywords = {}
1632 input = numpy.asarray(input)
1633 if numpy.iscomplexobj(input):
1634 raise TypeError('Complex type not supported')
1635 output = _ni_support._get_output(output, input)
1636 if filter_size < 1:
1637 raise RuntimeError('invalid filter size')
1638 axis = normalize_axis_index(axis, input.ndim)
1639 if (filter_size // 2 + origin < 0) or (filter_size // 2 + origin >=
1640 filter_size):
1641 raise ValueError('invalid origin')
1642 mode = _ni_support._extend_mode_to_code(mode)
1643 _nd_image.generic_filter1d(input, function, filter_size, axis, output,
1644 mode, cval, origin, extra_arguments,
1645 extra_keywords)
1646 return output
1649@_ni_docstrings.docfiller
1650def generic_filter(input, function, size=None, footprint=None,
1651 output=None, mode="reflect", cval=0.0, origin=0,
1652 extra_arguments=(), extra_keywords=None):
1653 """Calculate a multidimensional filter using the given function.
1655 At each element the provided function is called. The input values
1656 within the filter footprint at that element are passed to the function
1657 as a 1-D array of double values.
1659 Parameters
1660 ----------
1661 %(input)s
1662 function : {callable, scipy.LowLevelCallable}
1663 Function to apply at each element.
1664 %(size_foot)s
1665 %(output)s
1666 %(mode_reflect)s
1667 %(cval)s
1668 %(origin_multiple)s
1669 %(extra_arguments)s
1670 %(extra_keywords)s
1672 Notes
1673 -----
1674 This function also accepts low-level callback functions with one of
1675 the following signatures and wrapped in `scipy.LowLevelCallable`:
1677 .. code:: c
1679 int callback(double *buffer, npy_intp filter_size,
1680 double *return_value, void *user_data)
1681 int callback(double *buffer, intptr_t filter_size,
1682 double *return_value, void *user_data)
1684 The calling function iterates over the elements of the input and
1685 output arrays, calling the callback function at each element. The
1686 elements within the footprint of the filter at the current element are
1687 passed through the ``buffer`` parameter, and the number of elements
1688 within the footprint through ``filter_size``. The calculated value is
1689 returned in ``return_value``. ``user_data`` is the data pointer provided
1690 to `scipy.LowLevelCallable` as-is.
1692 The callback function must return an integer error status that is zero
1693 if something went wrong and one otherwise. If an error occurs, you should
1694 normally set the python error status with an informative message
1695 before returning, otherwise a default error message is set by the
1696 calling function.
1698 In addition, some other low-level function pointer specifications
1699 are accepted, but these are for backward compatibility only and should
1700 not be used in new code.
1702 Examples
1703 --------
1704 Import the necessary modules and load the example image used for
1705 filtering.
1707 >>> import numpy as np
1708 >>> from scipy import datasets
1709 >>> from scipy.ndimage import generic_filter
1710 >>> import matplotlib.pyplot as plt
1711 >>> ascent = datasets.ascent()
1713 Compute a maximum filter with kernel size 10 by passing a simple NumPy
1714 aggregation function as argument to `function`.
1716 >>> maximum_filter_result = generic_filter(ascent, np.amax, [10, 10])
1718 While a maximmum filter could also directly be obtained using
1719 `maximum_filter`, `generic_filter` allows generic Python function or
1720 `scipy.LowLevelCallable` to be used as a filter. Here, we compute the
1721 range between maximum and minimum value as an example for a kernel size
1722 of 5.
1724 >>> def custom_filter(image):
1725 ... return np.amax(image) - np.amin(image)
1726 >>> custom_filter_result = generic_filter(ascent, custom_filter, [5, 5])
1728 Plot the original and filtered images.
1730 >>> fig, axes = plt.subplots(3, 1, figsize=(4, 12))
1731 >>> plt.gray() # show the filtered result in grayscale
1732 >>> top, middle, bottom = axes
1733 >>> for ax in axes:
1734 ... ax.set_axis_off() # remove coordinate system
1735 >>> top.imshow(ascent)
1736 >>> top.set_title("Original image")
1737 >>> middle.imshow(maximum_filter_result)
1738 >>> middle.set_title("Maximum filter, Kernel: 10x10")
1739 >>> bottom.imshow(custom_filter_result)
1740 >>> bottom.set_title("Custom filter, Kernel: 5x5")
1741 >>> fig.tight_layout()
1743 """
1744 if (size is not None) and (footprint is not None):
1745 warnings.warn("ignoring size because footprint is set", UserWarning, stacklevel=2)
1746 if extra_keywords is None:
1747 extra_keywords = {}
1748 input = numpy.asarray(input)
1749 if numpy.iscomplexobj(input):
1750 raise TypeError('Complex type not supported')
1751 origins = _ni_support._normalize_sequence(origin, input.ndim)
1752 if footprint is None:
1753 if size is None:
1754 raise RuntimeError("no footprint or filter size provided")
1755 sizes = _ni_support._normalize_sequence(size, input.ndim)
1756 footprint = numpy.ones(sizes, dtype=bool)
1757 else:
1758 footprint = numpy.asarray(footprint, dtype=bool)
1759 fshape = [ii for ii in footprint.shape if ii > 0]
1760 if len(fshape) != input.ndim:
1761 raise RuntimeError('filter footprint array has incorrect shape.')
1762 for origin, lenf in zip(origins, fshape):
1763 if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf):
1764 raise ValueError('invalid origin')
1765 if not footprint.flags.contiguous:
1766 footprint = footprint.copy()
1767 output = _ni_support._get_output(output, input)
1768 mode = _ni_support._extend_mode_to_code(mode)
1769 _nd_image.generic_filter(input, function, footprint, output, mode,
1770 cval, origins, extra_arguments, extra_keywords)
1771 return output