diff --git a/src/bioxelnodes/bioxel/layer.py b/src/bioxelnodes/bioxel/layer.py index 5dd9a32..c212cc7 100644 --- a/src/bioxelnodes/bioxel/layer.py +++ b/src/bioxelnodes/bioxel/layer.py @@ -7,8 +7,7 @@ try: except ImportError: vdb = None -from . import scipy -from . import scipy as ndi +from scipy import ndimage as ndi # 3rd-party import transforms3d @@ -83,21 +82,17 @@ class Layer(): for f in range(self.frame_count): mask_frame = mask[f, :, :, :] if smooth > 0: - mask_frame = scipy.median_filter(mask_frame.astype(np.float32), - mode="nearest", - size=smooth) - # mask_frame = scipy.median_filter( - # mask_frame.astype(np.float32), size=2) + mask_frame = ndi.median_filter(mask_frame.astype(np.float32), + mode="nearest", + size=smooth) mask_frames += (mask_frame,) elif mask.ndim == 3: for f in range(self.frame_count): mask_frame = mask[:, :, :] if smooth > 0: - mask_frame = scipy.median_filter(mask_frame.astype(np.float32), - mode="nearest", - size=smooth) - # mask_frame = scipy.median_filter( - # mask_frame.astype(np.float32), size=2) + mask_frame = ndi.median_filter(mask_frame.astype(np.float32), + mode="nearest", + size=smooth) mask_frames += (mask_frame,) else: raise Exception("Mask shape order should be TXYZ or XYZ") @@ -127,9 +122,9 @@ class Layer(): frame = data[f, :, :, :, :] if smooth > 0: - frame = scipy.median_filter(frame.astype(np.float32), - mode="nearest", - size=smooth) + frame = ndi.median_filter(frame.astype(np.float32), + mode="nearest", + size=smooth) factors = np.divide(self.shape, shape) zoom_factors = [1 / f for f in factors] diff --git a/src/bioxelnodes/bioxel/scipy/__init__.py b/src/bioxelnodes/bioxel/scipy/__init__.py deleted file mode 100644 index c8b0ded..0000000 --- a/src/bioxelnodes/bioxel/scipy/__init__.py +++ /dev/null @@ -1,5 +0,0 @@ -from ._interpolation import zoom -from ._filters import (gaussian_filter, - median_filter, - maximum_filter, - minimum_filter) diff --git a/src/bioxelnodes/bioxel/scipy/_filters.py b/src/bioxelnodes/bioxel/scipy/_filters.py deleted file mode 100644 index 1e54ff4..0000000 --- a/src/bioxelnodes/bioxel/scipy/_filters.py +++ /dev/null @@ -1,1823 +0,0 @@ -# Copyright (C) 2003-2005 Peter J. Verveer -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# 1. Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# -# 2. Redistributions in binary form must reproduce the above -# copyright notice, this list of conditions and the following -# disclaimer in the documentation and/or other materials provided -# with the distribution. -# -# 3. The name of the author may not be used to endorse or promote -# products derived from this software without specific prior -# written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS -# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY -# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE -# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, -# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -from collections.abc import Iterable -import numbers -import warnings -import numpy as np -import operator - -from ._utils import normalize_axis_index -from . import _ni_support -from . import _nd_image - - -def _invalid_origin(origin, lenw): - return (origin < -(lenw // 2)) or (origin > (lenw - 1) // 2) - - -def _complex_via_real_components(func, input, weights, output, cval, **kwargs): - """Complex convolution via a linear combination of real convolutions.""" - complex_input = input.dtype.kind == 'c' - complex_weights = weights.dtype.kind == 'c' - if complex_input and complex_weights: - # real component of the output - func(input.real, weights.real, output=output.real, - cval=np.real(cval), **kwargs) - output.real -= func(input.imag, weights.imag, output=None, - cval=np.imag(cval), **kwargs) - # imaginary component of the output - func(input.real, weights.imag, output=output.imag, - cval=np.real(cval), **kwargs) - output.imag += func(input.imag, weights.real, output=None, - cval=np.imag(cval), **kwargs) - elif complex_input: - func(input.real, weights, output=output.real, cval=np.real(cval), - **kwargs) - func(input.imag, weights, output=output.imag, cval=np.imag(cval), - **kwargs) - else: - if np.iscomplexobj(cval): - raise ValueError("Cannot provide a complex-valued cval when the " - "input is real.") - func(input, weights.real, output=output.real, cval=cval, **kwargs) - func(input, weights.imag, output=output.imag, cval=cval, **kwargs) - return output - - -def correlate1d(input, weights, axis=-1, output=None, mode="reflect", - cval=0.0, origin=0): - """Calculate a 1-D correlation along the given axis. - - The lines of the array along the given axis are correlated with the - given weights. - - Parameters - ---------- - %(input)s - weights : array - 1-D sequence of numbers. - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - - Returns - ------- - result : ndarray - Correlation result. Has the same shape as `input`. - - Examples - -------- - >>> from scipy.ndimage import correlate1d - >>> correlate1d([2, 8, 0, 4, 1, 9, 9, 0], weights=[1, 3]) - array([ 8, 26, 8, 12, 7, 28, 36, 9]) - """ - input = np.asarray(input) - weights = np.asarray(weights) - complex_input = input.dtype.kind == 'c' - complex_weights = weights.dtype.kind == 'c' - if complex_input or complex_weights: - if complex_weights: - weights = weights.conj() - weights = weights.astype(np.complex128, copy=False) - kwargs = dict(axis=axis, mode=mode, origin=origin) - output = _ni_support._get_output(output, input, complex_output=True) - return _complex_via_real_components(correlate1d, input, weights, - output, cval, **kwargs) - - output = _ni_support._get_output(output, input) - weights = np.asarray(weights, dtype=np.float64) - if weights.ndim != 1 or weights.shape[0] < 1: - raise RuntimeError('no filter weights given') - if not weights.flags.contiguous: - weights = weights.copy() - axis = normalize_axis_index(axis, input.ndim) - if _invalid_origin(origin, len(weights)): - raise ValueError('Invalid origin; origin must satisfy ' - '-(len(weights) // 2) <= origin <= ' - '(len(weights)-1) // 2') - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.correlate1d(input, weights, axis, output, mode, cval, - origin) - return output - - -def convolve1d(input, weights, axis=-1, output=None, mode="reflect", - cval=0.0, origin=0): - """Calculate a 1-D convolution along the given axis. - - The lines of the array along the given axis are convolved with the - given weights. - - Parameters - ---------- - %(input)s - weights : ndarray - 1-D sequence of numbers. - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - - Returns - ------- - convolve1d : ndarray - Convolved array with same shape as input - - Examples - -------- - >>> from scipy.ndimage import convolve1d - >>> convolve1d([2, 8, 0, 4, 1, 9, 9, 0], weights=[1, 3]) - array([14, 24, 4, 13, 12, 36, 27, 0]) - """ - weights = weights[::-1] - origin = -origin - if not len(weights) & 1: - origin -= 1 - weights = np.asarray(weights) - if weights.dtype.kind == 'c': - # pre-conjugate here to counteract the conjugation in correlate1d - weights = weights.conj() - return correlate1d(input, weights, axis, output, mode, cval, origin) - - -def _gaussian_kernel1d(sigma, order, radius): - """ - Computes a 1-D Gaussian convolution kernel. - """ - if order < 0: - raise ValueError('order must be non-negative') - exponent_range = np.arange(order + 1) - sigma2 = sigma * sigma - x = np.arange(-radius, radius+1) - phi_x = np.exp(-0.5 / sigma2 * x ** 2) - phi_x = phi_x / phi_x.sum() - - if order == 0: - return phi_x - else: - # f(x) = q(x) * phi(x) = q(x) * exp(p(x)) - # f'(x) = (q'(x) + q(x) * p'(x)) * phi(x) - # p'(x) = -1 / sigma ** 2 - # Implement q'(x) + q(x) * p'(x) as a matrix operator and apply to the - # coefficients of q(x) - q = np.zeros(order + 1) - q[0] = 1 - D = np.diag(exponent_range[1:], 1) # D @ q(x) = q'(x) - P = np.diag(np.ones(order)/-sigma2, -1) # P @ q(x) = q(x) * p'(x) - Q_deriv = D + P - for _ in range(order): - q = Q_deriv.dot(q) - q = (x[:, None] ** exponent_range).dot(q) - return q * phi_x - - -def gaussian_filter1d(input, sigma, axis=-1, order=0, output=None, - mode="reflect", cval=0.0, truncate=4.0, *, radius=None): - """1-D Gaussian filter. - - Parameters - ---------- - %(input)s - sigma : scalar - standard deviation for Gaussian kernel - %(axis)s - order : int, optional - An order of 0 corresponds to convolution with a Gaussian - kernel. A positive order corresponds to convolution with - that derivative of a Gaussian. - %(output)s - %(mode_reflect)s - %(cval)s - truncate : float, optional - Truncate the filter at this many standard deviations. - Default is 4.0. - radius : None or int, optional - Radius of the Gaussian kernel. If specified, the size of - the kernel will be ``2*radius + 1``, and `truncate` is ignored. - Default is None. - - Returns - ------- - gaussian_filter1d : ndarray - - Notes - ----- - The Gaussian kernel will have size ``2*radius + 1`` along each axis. If - `radius` is None, a default ``radius = round(truncate * sigma)`` will be - used. - - Examples - -------- - >>> from scipy.ndimage import gaussian_filter1d - >>> import numpy as np - >>> gaussian_filter1d([1.0, 2.0, 3.0, 4.0, 5.0], 1) - array([ 1.42704095, 2.06782203, 3. , 3.93217797, 4.57295905]) - >>> gaussian_filter1d([1.0, 2.0, 3.0, 4.0, 5.0], 4) - array([ 2.91948343, 2.95023502, 3. , 3.04976498, 3.08051657]) - >>> import matplotlib.pyplot as plt - >>> rng = np.random.default_rng() - >>> x = rng.standard_normal(101).cumsum() - >>> y3 = gaussian_filter1d(x, 3) - >>> y6 = gaussian_filter1d(x, 6) - >>> plt.plot(x, 'k', label='original data') - >>> plt.plot(y3, '--', label='filtered, sigma=3') - >>> plt.plot(y6, ':', label='filtered, sigma=6') - >>> plt.legend() - >>> plt.grid() - >>> plt.show() - - """ - sd = float(sigma) - # make the radius of the filter equal to truncate standard deviations - lw = int(truncate * sd + 0.5) - if radius is not None: - lw = radius - if not isinstance(lw, numbers.Integral) or lw < 0: - raise ValueError('Radius must be a nonnegative integer.') - # Since we are calling correlate, not convolve, revert the kernel - weights = _gaussian_kernel1d(sigma, order, lw)[::-1] - return correlate1d(input, weights, axis, output, mode, cval, 0) - - -def gaussian_filter(input, sigma, order=0, output=None, - mode="reflect", cval=0.0, truncate=4.0, *, radius=None, - axes=None): - """Multidimensional Gaussian filter. - - Parameters - ---------- - %(input)s - sigma : scalar or sequence of scalars - Standard deviation for Gaussian kernel. The standard - deviations of the Gaussian filter are given for each axis as a - sequence, or as a single number, in which case it is equal for - all axes. - order : int or sequence of ints, optional - The order of the filter along each axis is given as a sequence - of integers, or as a single number. An order of 0 corresponds - to convolution with a Gaussian kernel. A positive order - corresponds to convolution with that derivative of a Gaussian. - %(output)s - %(mode_multiple)s - %(cval)s - truncate : float, optional - Truncate the filter at this many standard deviations. - Default is 4.0. - radius : None or int or sequence of ints, optional - Radius of the Gaussian kernel. The radius are given for each axis - as a sequence, or as a single number, in which case it is equal - for all axes. If specified, the size of the kernel along each axis - will be ``2*radius + 1``, and `truncate` is ignored. - Default is None. - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. When `axes` is - specified, any tuples used for `sigma`, `order`, `mode` and/or `radius` - must match the length of `axes`. The ith entry in any of these tuples - corresponds to the ith entry in `axes`. - - Returns - ------- - gaussian_filter : ndarray - Returned array of same shape as `input`. - - Notes - ----- - The multidimensional filter is implemented as a sequence of - 1-D convolution filters. The intermediate arrays are - stored in the same data type as the output. Therefore, for output - types with a limited precision, the results may be imprecise - because intermediate results may be stored with insufficient - precision. - - The Gaussian kernel will have size ``2*radius + 1`` along each axis. If - `radius` is None, the default ``radius = round(truncate * sigma)`` will be - used. - - Examples - -------- - >>> from scipy.ndimage import gaussian_filter - >>> import numpy as np - >>> a = np.arange(50, step=2).reshape((5,5)) - >>> a - array([[ 0, 2, 4, 6, 8], - [10, 12, 14, 16, 18], - [20, 22, 24, 26, 28], - [30, 32, 34, 36, 38], - [40, 42, 44, 46, 48]]) - >>> gaussian_filter(a, sigma=1) - array([[ 4, 6, 8, 9, 11], - [10, 12, 14, 15, 17], - [20, 22, 24, 25, 27], - [29, 31, 33, 34, 36], - [35, 37, 39, 40, 42]]) - - >>> from scipy import datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = gaussian_filter(ascent, sigma=5) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - input = np.asarray(input) - output = _ni_support._get_output(output, input) - - axes = _ni_support._check_axes(axes, input.ndim) - num_axes = len(axes) - orders = _ni_support._normalize_sequence(order, num_axes) - sigmas = _ni_support._normalize_sequence(sigma, num_axes) - modes = _ni_support._normalize_sequence(mode, num_axes) - radiuses = _ni_support._normalize_sequence(radius, num_axes) - axes = [(axes[ii], sigmas[ii], orders[ii], modes[ii], radiuses[ii]) - for ii in range(num_axes) if sigmas[ii] > 1e-15] - if len(axes) > 0: - for axis, sigma, order, mode, radius in axes: - gaussian_filter1d(input, sigma, axis, order, output, - mode, cval, truncate, radius=radius) - input = output - else: - output[...] = input[...] - return output - - -def prewitt(input, axis=-1, output=None, mode="reflect", cval=0.0): - """Calculate a Prewitt filter. - - Parameters - ---------- - %(input)s - %(axis)s - %(output)s - %(mode_multiple)s - %(cval)s - - Returns - ------- - prewitt : ndarray - Filtered array. Has the same shape as `input`. - - See Also - -------- - sobel: Sobel filter - - Notes - ----- - This function computes the one-dimensional Prewitt filter. - Horizontal edges are emphasised with the horizontal transform (axis=0), - vertical edges with the vertical transform (axis=1), and so on for higher - dimensions. These can be combined to give the magnitude. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> import numpy as np - >>> ascent = datasets.ascent() - >>> prewitt_h = ndimage.prewitt(ascent, axis=0) - >>> prewitt_v = ndimage.prewitt(ascent, axis=1) - >>> magnitude = np.sqrt(prewitt_h ** 2 + prewitt_v ** 2) - >>> magnitude *= 255 / np.max(magnitude) # Normalization - >>> fig, axes = plt.subplots(2, 2, figsize = (8, 8)) - >>> plt.gray() - >>> axes[0, 0].imshow(ascent) - >>> axes[0, 1].imshow(prewitt_h) - >>> axes[1, 0].imshow(prewitt_v) - >>> axes[1, 1].imshow(magnitude) - >>> titles = ["original", "horizontal", "vertical", "magnitude"] - >>> for i, ax in enumerate(axes.ravel()): - ... ax.set_title(titles[i]) - ... ax.axis("off") - >>> plt.show() - - """ - input = np.asarray(input) - axis = normalize_axis_index(axis, input.ndim) - output = _ni_support._get_output(output, input) - modes = _ni_support._normalize_sequence(mode, input.ndim) - correlate1d(input, [-1, 0, 1], axis, output, modes[axis], cval, 0) - axes = [ii for ii in range(input.ndim) if ii != axis] - for ii in axes: - correlate1d(output, [1, 1, 1], ii, output, modes[ii], cval, 0,) - return output - - -def sobel(input, axis=-1, output=None, mode="reflect", cval=0.0): - """Calculate a Sobel filter. - - Parameters - ---------- - %(input)s - %(axis)s - %(output)s - %(mode_multiple)s - %(cval)s - - Returns - ------- - sobel : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - This function computes the axis-specific Sobel gradient. - The horizontal edges can be emphasised with the horizontal transform (axis=0), - the vertical edges with the vertical transform (axis=1) and so on for higher - dimensions. These can be combined to give the magnitude. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> import numpy as np - >>> ascent = datasets.ascent().astype('int32') - >>> sobel_h = ndimage.sobel(ascent, 0) # horizontal gradient - >>> sobel_v = ndimage.sobel(ascent, 1) # vertical gradient - >>> magnitude = np.sqrt(sobel_h**2 + sobel_v**2) - >>> magnitude *= 255.0 / np.max(magnitude) # normalization - >>> fig, axs = plt.subplots(2, 2, figsize=(8, 8)) - >>> plt.gray() # show the filtered result in grayscale - >>> axs[0, 0].imshow(ascent) - >>> axs[0, 1].imshow(sobel_h) - >>> axs[1, 0].imshow(sobel_v) - >>> axs[1, 1].imshow(magnitude) - >>> titles = ["original", "horizontal", "vertical", "magnitude"] - >>> for i, ax in enumerate(axs.ravel()): - ... ax.set_title(titles[i]) - ... ax.axis("off") - >>> plt.show() - - """ - input = np.asarray(input) - axis = normalize_axis_index(axis, input.ndim) - output = _ni_support._get_output(output, input) - modes = _ni_support._normalize_sequence(mode, input.ndim) - correlate1d(input, [-1, 0, 1], axis, output, modes[axis], cval, 0) - axes = [ii for ii in range(input.ndim) if ii != axis] - for ii in axes: - correlate1d(output, [1, 2, 1], ii, output, modes[ii], cval, 0) - return output - - -def generic_laplace(input, derivative2, output=None, mode="reflect", - cval=0.0, - extra_arguments=(), - extra_keywords=None): - """ - N-D Laplace filter using a provided second derivative function. - - Parameters - ---------- - %(input)s - derivative2 : callable - Callable with the following signature:: - - derivative2(input, axis, output, mode, cval, - *extra_arguments, **extra_keywords) - - See `extra_arguments`, `extra_keywords` below. - %(output)s - %(mode_multiple)s - %(cval)s - %(extra_keywords)s - %(extra_arguments)s - - Returns - ------- - generic_laplace : ndarray - Filtered array. Has the same shape as `input`. - - """ - if extra_keywords is None: - extra_keywords = {} - input = np.asarray(input) - output = _ni_support._get_output(output, input) - axes = list(range(input.ndim)) - if len(axes) > 0: - modes = _ni_support._normalize_sequence(mode, len(axes)) - derivative2(input, axes[0], output, modes[0], cval, - *extra_arguments, **extra_keywords) - for ii in range(1, len(axes)): - tmp = derivative2(input, axes[ii], output.dtype, modes[ii], cval, - *extra_arguments, **extra_keywords) - output += tmp - else: - output[...] = input[...] - return output - - -def laplace(input, output=None, mode="reflect", cval=0.0): - """N-D Laplace filter based on approximate second derivatives. - - Parameters - ---------- - %(input)s - %(output)s - %(mode_multiple)s - %(cval)s - - Returns - ------- - laplace : ndarray - Filtered array. Has the same shape as `input`. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.laplace(ascent) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - def derivative2(input, axis, output, mode, cval): - return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0) - return generic_laplace(input, derivative2, output, mode, cval) - - -def gaussian_laplace(input, sigma, output=None, mode="reflect", - cval=0.0, **kwargs): - """Multidimensional Laplace filter using Gaussian second derivatives. - - Parameters - ---------- - %(input)s - sigma : scalar or sequence of scalars - The standard deviations of the Gaussian filter are given for - each axis as a sequence, or as a single number, in which case - it is equal for all axes. - %(output)s - %(mode_multiple)s - %(cval)s - Extra keyword arguments will be passed to gaussian_filter(). - - Returns - ------- - gaussian_laplace : ndarray - Filtered array. Has the same shape as `input`. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> ascent = datasets.ascent() - - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - - >>> result = ndimage.gaussian_laplace(ascent, sigma=1) - >>> ax1.imshow(result) - - >>> result = ndimage.gaussian_laplace(ascent, sigma=3) - >>> ax2.imshow(result) - >>> plt.show() - """ - input = np.asarray(input) - - def derivative2(input, axis, output, mode, cval, sigma, **kwargs): - order = [0] * input.ndim - order[axis] = 2 - return gaussian_filter(input, sigma, order, output, mode, cval, - **kwargs) - - return generic_laplace(input, derivative2, output, mode, cval, - extra_arguments=(sigma,), - extra_keywords=kwargs) - - -def generic_gradient_magnitude(input, derivative, output=None, - mode="reflect", cval=0.0, - extra_arguments=(), extra_keywords=None): - """Gradient magnitude using a provided gradient function. - - Parameters - ---------- - %(input)s - derivative : callable - Callable with the following signature:: - - derivative(input, axis, output, mode, cval, - *extra_arguments, **extra_keywords) - - See `extra_arguments`, `extra_keywords` below. - `derivative` can assume that `input` and `output` are ndarrays. - Note that the output from `derivative` is modified inplace; - be careful to copy important inputs before returning them. - %(output)s - %(mode_multiple)s - %(cval)s - %(extra_keywords)s - %(extra_arguments)s - - Returns - ------- - generic_gradient_matnitude : ndarray - Filtered array. Has the same shape as `input`. - - """ - if extra_keywords is None: - extra_keywords = {} - input = np.asarray(input) - output = _ni_support._get_output(output, input) - axes = list(range(input.ndim)) - if len(axes) > 0: - modes = _ni_support._normalize_sequence(mode, len(axes)) - derivative(input, axes[0], output, modes[0], cval, - *extra_arguments, **extra_keywords) - np.multiply(output, output, output) - for ii in range(1, len(axes)): - tmp = derivative(input, axes[ii], output.dtype, modes[ii], cval, - *extra_arguments, **extra_keywords) - np.multiply(tmp, tmp, tmp) - output += tmp - # This allows the sqrt to work with a different default casting - np.sqrt(output, output, casting='unsafe') - else: - output[...] = input[...] - return output - - -def gaussian_gradient_magnitude(input, sigma, output=None, - mode="reflect", cval=0.0, **kwargs): - """Multidimensional gradient magnitude using Gaussian derivatives. - - Parameters - ---------- - %(input)s - sigma : scalar or sequence of scalars - The standard deviations of the Gaussian filter are given for - each axis as a sequence, or as a single number, in which case - it is equal for all axes. - %(output)s - %(mode_multiple)s - %(cval)s - Extra keyword arguments will be passed to gaussian_filter(). - - Returns - ------- - gaussian_gradient_magnitude : ndarray - Filtered array. Has the same shape as `input`. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.gaussian_gradient_magnitude(ascent, sigma=5) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - input = np.asarray(input) - - def derivative(input, axis, output, mode, cval, sigma, **kwargs): - order = [0] * input.ndim - order[axis] = 1 - return gaussian_filter(input, sigma, order, output, mode, - cval, **kwargs) - - return generic_gradient_magnitude(input, derivative, output, mode, - cval, extra_arguments=(sigma,), - extra_keywords=kwargs) - - -def _correlate_or_convolve(input, weights, output, mode, cval, origin, - convolution): - input = np.asarray(input) - weights = np.asarray(weights) - complex_input = input.dtype.kind == 'c' - complex_weights = weights.dtype.kind == 'c' - if complex_input or complex_weights: - if complex_weights and not convolution: - # As for np.correlate, conjugate weights rather than input. - weights = weights.conj() - kwargs = dict( - mode=mode, origin=origin, convolution=convolution - ) - output = _ni_support._get_output(output, input, complex_output=True) - - return _complex_via_real_components(_correlate_or_convolve, input, - weights, output, cval, **kwargs) - - origins = _ni_support._normalize_sequence(origin, input.ndim) - weights = np.asarray(weights, dtype=np.float64) - wshape = [ii for ii in weights.shape if ii > 0] - if len(wshape) != input.ndim: - raise RuntimeError('filter weights array has incorrect shape.') - if convolution: - weights = weights[tuple([slice(None, None, -1)] * weights.ndim)] - for ii in range(len(origins)): - origins[ii] = -origins[ii] - if not weights.shape[ii] & 1: - origins[ii] -= 1 - for origin, lenw in zip(origins, wshape): - if _invalid_origin(origin, lenw): - raise ValueError('Invalid origin; origin must satisfy ' - '-(weights.shape[k] // 2) <= origin[k] <= ' - '(weights.shape[k]-1) // 2') - - if not weights.flags.contiguous: - weights = weights.copy() - output = _ni_support._get_output(output, input) - temp_needed = np.may_share_memory(input, output) - if temp_needed: - # input and output arrays cannot share memory - temp = output - output = _ni_support._get_output(output.dtype, input) - if not isinstance(mode, str) and isinstance(mode, Iterable): - raise RuntimeError("A sequence of modes is not supported") - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.correlate(input, weights, output, mode, cval, origins) - if temp_needed: - temp[...] = output - output = temp - return output - - -def correlate(input, weights, output=None, mode='reflect', cval=0.0, - origin=0): - """ - Multidimensional correlation. - - The array is correlated with the given kernel. - - Parameters - ---------- - %(input)s - weights : ndarray - array of weights, same number of dimensions as input - %(output)s - %(mode_reflect)s - %(cval)s - %(origin_multiple)s - - Returns - ------- - result : ndarray - The result of correlation of `input` with `weights`. - - See Also - -------- - convolve : Convolve an image with a kernel. - - Examples - -------- - Correlation is the process of moving a filter mask often referred to - as kernel over the image and computing the sum of products at each location. - - >>> from scipy.ndimage import correlate - >>> import numpy as np - >>> input_img = np.arange(25).reshape(5,5) - >>> print(input_img) - [[ 0 1 2 3 4] - [ 5 6 7 8 9] - [10 11 12 13 14] - [15 16 17 18 19] - [20 21 22 23 24]] - - Define a kernel (weights) for correlation. In this example, it is for sum of - center and up, down, left and right next elements. - - >>> weights = [[0, 1, 0], - ... [1, 1, 1], - ... [0, 1, 0]] - - We can calculate a correlation result: - For example, element ``[2,2]`` is ``7 + 11 + 12 + 13 + 17 = 60``. - - >>> correlate(input_img, weights) - array([[ 6, 10, 15, 20, 24], - [ 26, 30, 35, 40, 44], - [ 51, 55, 60, 65, 69], - [ 76, 80, 85, 90, 94], - [ 96, 100, 105, 110, 114]]) - - """ - return _correlate_or_convolve(input, weights, output, mode, cval, - origin, False) - - -def convolve(input, weights, output=None, mode='reflect', cval=0.0, - origin=0): - """ - Multidimensional convolution. - - The array is convolved with the given kernel. - - Parameters - ---------- - %(input)s - weights : array_like - Array of weights, same number of dimensions as input - %(output)s - %(mode_reflect)s - cval : scalar, optional - Value to fill past edges of input if `mode` is 'constant'. Default - is 0.0 - origin : int, optional - Controls the origin of the input signal, which is where the - filter is centered to produce the first element of the output. - Positive values shift the filter to the right, and negative values - shift the filter to the left. Default is 0. - - Returns - ------- - result : ndarray - The result of convolution of `input` with `weights`. - - See Also - -------- - correlate : Correlate an image with a kernel. - - Notes - ----- - Each value in result is :math:`C_i = \\sum_j{I_{i+k-j} W_j}`, where - W is the `weights` kernel, - j is the N-D spatial index over :math:`W`, - I is the `input` and k is the coordinate of the center of - W, specified by `origin` in the input parameters. - - Examples - -------- - Perhaps the simplest case to understand is ``mode='constant', cval=0.0``, - because in this case borders (i.e., where the `weights` kernel, centered - on any one value, extends beyond an edge of `input`) are treated as zeros. - - >>> import numpy as np - >>> a = np.array([[1, 2, 0, 0], - ... [5, 3, 0, 4], - ... [0, 0, 0, 7], - ... [9, 3, 0, 0]]) - >>> k = np.array([[1,1,1],[1,1,0],[1,0,0]]) - >>> from scipy import ndimage - >>> ndimage.convolve(a, k, mode='constant', cval=0.0) - array([[11, 10, 7, 4], - [10, 3, 11, 11], - [15, 12, 14, 7], - [12, 3, 7, 0]]) - - Setting ``cval=1.0`` is equivalent to padding the outer edge of `input` - with 1.0's (and then extracting only the original region of the result). - - >>> ndimage.convolve(a, k, mode='constant', cval=1.0) - array([[13, 11, 8, 7], - [11, 3, 11, 14], - [16, 12, 14, 10], - [15, 6, 10, 5]]) - - With ``mode='reflect'`` (the default), outer values are reflected at the - edge of `input` to fill in missing values. - - >>> b = np.array([[2, 0, 0], - ... [1, 0, 0], - ... [0, 0, 0]]) - >>> k = np.array([[0,1,0], [0,1,0], [0,1,0]]) - >>> ndimage.convolve(b, k, mode='reflect') - array([[5, 0, 0], - [3, 0, 0], - [1, 0, 0]]) - - This includes diagonally at the corners. - - >>> k = np.array([[1,0,0],[0,1,0],[0,0,1]]) - >>> ndimage.convolve(b, k) - array([[4, 2, 0], - [3, 2, 0], - [1, 1, 0]]) - - With ``mode='nearest'``, the single nearest value in to an edge in - `input` is repeated as many times as needed to match the overlapping - `weights`. - - >>> c = np.array([[2, 0, 1], - ... [1, 0, 0], - ... [0, 0, 0]]) - >>> k = np.array([[0, 1, 0], - ... [0, 1, 0], - ... [0, 1, 0], - ... [0, 1, 0], - ... [0, 1, 0]]) - >>> ndimage.convolve(c, k, mode='nearest') - array([[7, 0, 3], - [5, 0, 2], - [3, 0, 1]]) - - """ - return _correlate_or_convolve(input, weights, output, mode, cval, - origin, True) - - -def uniform_filter1d(input, size, axis=-1, output=None, - mode="reflect", cval=0.0, origin=0): - """Calculate a 1-D uniform filter along the given axis. - - The lines of the array along the given axis are filtered with a - uniform filter of given size. - - Parameters - ---------- - %(input)s - size : int - length of uniform filter - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - - Returns - ------- - result : ndarray - Filtered array. Has same shape as `input`. - - Examples - -------- - >>> from scipy.ndimage import uniform_filter1d - >>> uniform_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3) - array([4, 3, 4, 1, 4, 6, 6, 3]) - """ - input = np.asarray(input) - axis = normalize_axis_index(axis, input.ndim) - if size < 1: - raise RuntimeError('incorrect filter size') - complex_output = input.dtype.kind == 'c' - output = _ni_support._get_output(output, input, - complex_output=complex_output) - if (size // 2 + origin < 0) or (size // 2 + origin >= size): - raise ValueError('invalid origin') - mode = _ni_support._extend_mode_to_code(mode) - if not complex_output: - _nd_image.uniform_filter1d(input, size, axis, output, mode, cval, - origin) - else: - _nd_image.uniform_filter1d(input.real, size, axis, output.real, mode, - np.real(cval), origin) - _nd_image.uniform_filter1d(input.imag, size, axis, output.imag, mode, - np.imag(cval), origin) - return output - - -def uniform_filter(input, size=3, output=None, mode="reflect", - cval=0.0, origin=0, *, axes=None): - """Multidimensional uniform filter. - - Parameters - ---------- - %(input)s - size : int or sequence of ints, optional - The sizes of the uniform filter are given for each axis as a - sequence, or as a single number, in which case the size is - equal for all axes. - %(output)s - %(mode_multiple)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. When `axes` is - specified, any tuples used for `size`, `origin`, and/or `mode` - must match the length of `axes`. The ith entry in any of these tuples - corresponds to the ith entry in `axes`. - - Returns - ------- - uniform_filter : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - The multidimensional filter is implemented as a sequence of - 1-D uniform filters. The intermediate arrays are stored - in the same data type as the output. Therefore, for output types - with a limited precision, the results may be imprecise because - intermediate results may be stored with insufficient precision. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.uniform_filter(ascent, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - input = np.asarray(input) - output = _ni_support._get_output(output, input, - complex_output=input.dtype.kind == 'c') - axes = _ni_support._check_axes(axes, input.ndim) - num_axes = len(axes) - sizes = _ni_support._normalize_sequence(size, num_axes) - origins = _ni_support._normalize_sequence(origin, num_axes) - modes = _ni_support._normalize_sequence(mode, num_axes) - axes = [(axes[ii], sizes[ii], origins[ii], modes[ii]) - for ii in range(num_axes) if sizes[ii] > 1] - if len(axes) > 0: - for axis, size, origin, mode in axes: - uniform_filter1d(input, int(size), axis, output, mode, - cval, origin) - input = output - else: - output[...] = input[...] - return output - - -def minimum_filter1d(input, size, axis=-1, output=None, - mode="reflect", cval=0.0, origin=0): - """Calculate a 1-D minimum filter along the given axis. - - The lines of the array along the given axis are filtered with a - minimum filter of given size. - - Parameters - ---------- - %(input)s - size : int - length along which to calculate 1D minimum - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - - Returns - ------- - result : ndarray. - Filtered image. Has the same shape as `input`. - - 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 - - - Examples - -------- - >>> from scipy.ndimage import minimum_filter1d - >>> minimum_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3) - array([2, 0, 0, 0, 1, 1, 0, 0]) - """ - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError('Complex type not supported') - axis = normalize_axis_index(axis, input.ndim) - if size < 1: - raise RuntimeError('incorrect filter size') - output = _ni_support._get_output(output, input) - if (size // 2 + origin < 0) or (size // 2 + origin >= size): - raise ValueError('invalid origin') - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.min_or_max_filter1d(input, size, axis, output, mode, cval, - origin, 1) - return output - - -def maximum_filter1d(input, size, axis=-1, output=None, - mode="reflect", cval=0.0, origin=0): - """Calculate a 1-D maximum filter along the given axis. - - The lines of the array along the given axis are filtered with a - maximum filter of given size. - - Parameters - ---------- - %(input)s - size : int - Length along which to calculate the 1-D maximum. - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - - Returns - ------- - maximum1d : ndarray, 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 - - Examples - -------- - >>> from scipy.ndimage import maximum_filter1d - >>> maximum_filter1d([2, 8, 0, 4, 1, 9, 9, 0], size=3) - array([8, 8, 8, 4, 9, 9, 9, 9]) - """ - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError('Complex type not supported') - axis = normalize_axis_index(axis, input.ndim) - if size < 1: - raise RuntimeError('incorrect filter size') - output = _ni_support._get_output(output, input) - if (size // 2 + origin < 0) or (size // 2 + origin >= size): - raise ValueError('invalid origin') - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.min_or_max_filter1d(input, size, axis, output, mode, cval, - origin, 0) - return output - - -def _min_or_max_filter(input, size, footprint, structure, output, mode, - cval, origin, minimum, axes=None): - if (size is not None) and (footprint is not None): - warnings.warn("ignoring size because footprint is set", - UserWarning, stacklevel=3) - if structure is None: - if footprint is None: - if size is None: - raise RuntimeError("no footprint provided") - separable = True - else: - footprint = np.asarray(footprint, dtype=bool) - if not footprint.any(): - raise ValueError("All-zero footprint is not supported.") - if footprint.all(): - size = footprint.shape - footprint = None - separable = True - else: - separable = False - else: - structure = np.asarray(structure, dtype=np.float64) - separable = False - if footprint is None: - footprint = np.ones(structure.shape, bool) - else: - footprint = np.asarray(footprint, dtype=bool) - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError("Complex type not supported") - output = _ni_support._get_output(output, input) - temp_needed = np.may_share_memory(input, output) - if temp_needed: - # input and output arrays cannot share memory - temp = output - output = _ni_support._get_output(output.dtype, input) - axes = _ni_support._check_axes(axes, input.ndim) - num_axes = len(axes) - if separable: - origins = _ni_support._normalize_sequence(origin, num_axes) - sizes = _ni_support._normalize_sequence(size, num_axes) - modes = _ni_support._normalize_sequence(mode, num_axes) - axes = [(axes[ii], sizes[ii], origins[ii], modes[ii]) - for ii in range(len(axes)) if sizes[ii] > 1] - if minimum: - filter_ = minimum_filter1d - else: - filter_ = maximum_filter1d - if len(axes) > 0: - for axis, size, origin, mode in axes: - filter_(input, int(size), axis, output, mode, cval, origin) - input = output - else: - output[...] = input[...] - else: - origins = _ni_support._normalize_sequence(origin, num_axes) - if num_axes < input.ndim: - if footprint.ndim != num_axes: - raise RuntimeError("footprint array has incorrect shape") - footprint = np.expand_dims( - footprint, - tuple(ax for ax in range(input.ndim) if ax not in axes) - ) - # set origin = 0 for any axes not being filtered - origins_temp = [0,] * input.ndim - for o, ax in zip(origins, axes): - origins_temp[ax] = o - origins = origins_temp - - fshape = [ii for ii in footprint.shape if ii > 0] - if len(fshape) != input.ndim: - raise RuntimeError('footprint array has incorrect shape.') - for origin, lenf in zip(origins, fshape): - if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf): - raise ValueError("invalid origin") - if not footprint.flags.contiguous: - footprint = footprint.copy() - if structure is not None: - if len(structure.shape) != input.ndim: - raise RuntimeError("structure array has incorrect shape") - if num_axes != structure.ndim: - structure = np.expand_dims( - structure, - tuple(ax for ax in range(structure.ndim) if ax not in axes) - ) - if not structure.flags.contiguous: - structure = structure.copy() - if not isinstance(mode, str) and isinstance(mode, Iterable): - raise RuntimeError( - "A sequence of modes is not supported for non-separable " - "footprints") - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.min_or_max_filter(input, footprint, structure, output, - mode, cval, origins, minimum) - if temp_needed: - temp[...] = output - output = temp - return output - - -def minimum_filter(input, size=None, footprint=None, output=None, - mode="reflect", cval=0.0, origin=0, *, axes=None): - """Calculate a multidimensional minimum filter. - - Parameters - ---------- - %(input)s - %(size_foot)s - %(output)s - %(mode_multiple)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. When `axes` is - specified, any tuples used for `size`, `origin`, and/or `mode` - must match the length of `axes`. The ith entry in any of these tuples - corresponds to the ith entry in `axes`. - - Returns - ------- - minimum_filter : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - A sequence of modes (one per axis) is only supported when the footprint is - separable. Otherwise, a single mode string must be provided. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.minimum_filter(ascent, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - return _min_or_max_filter(input, size, footprint, None, output, mode, - cval, origin, 1, axes) - - -def maximum_filter(input, size=None, footprint=None, output=None, - mode="reflect", cval=0.0, origin=0, *, axes=None): - """Calculate a multidimensional maximum filter. - - Parameters - ---------- - %(input)s - %(size_foot)s - %(output)s - %(mode_multiple)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. When `axes` is - specified, any tuples used for `size`, `origin`, and/or `mode` - must match the length of `axes`. The ith entry in any of these tuples - corresponds to the ith entry in `axes`. - - Returns - ------- - maximum_filter : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - A sequence of modes (one per axis) is only supported when the footprint is - separable. Otherwise, a single mode string must be provided. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.maximum_filter(ascent, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - return _min_or_max_filter(input, size, footprint, None, output, mode, - cval, origin, 0, axes) - - -def _rank_filter(input, rank, size=None, footprint=None, output=None, - mode="reflect", cval=0.0, origin=0, operation='rank', - axes=None): - if (size is not None) and (footprint is not None): - warnings.warn("ignoring size because footprint is set", - UserWarning, stacklevel=3) - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError('Complex type not supported') - axes = _ni_support._check_axes(axes, input.ndim) - num_axes = len(axes) - origins = _ni_support._normalize_sequence(origin, num_axes) - if footprint is None: - if size is None: - raise RuntimeError("no footprint or filter size provided") - sizes = _ni_support._normalize_sequence(size, num_axes) - footprint = np.ones(sizes, dtype=bool) - else: - footprint = np.asarray(footprint, dtype=bool) - if num_axes < input.ndim: - # set origin = 0 for any axes not being filtered - origins_temp = [0,] * input.ndim - for o, ax in zip(origins, axes): - origins_temp[ax] = o - origins = origins_temp - - if not isinstance(mode, str) and isinstance(mode, Iterable): - # set mode = 'constant' for any axes not being filtered - modes = _ni_support._normalize_sequence(mode, num_axes) - modes_temp = ['constant'] * input.ndim - for m, ax in zip(modes, axes): - modes_temp[ax] = m - mode = modes_temp - - # insert singleton dimension along any non-filtered axes - if footprint.ndim != num_axes: - raise RuntimeError("footprint array has incorrect shape") - footprint = np.expand_dims( - footprint, - tuple(ax for ax in range(input.ndim) if ax not in axes) - ) - fshape = [ii for ii in footprint.shape if ii > 0] - if len(fshape) != input.ndim: - raise RuntimeError('footprint array has incorrect shape.') - for origin, lenf in zip(origins, fshape): - if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf): - raise ValueError('invalid origin') - if not footprint.flags.contiguous: - footprint = footprint.copy() - filter_size = np.where(footprint, 1, 0).sum() - if operation == 'median': - rank = filter_size // 2 - elif operation == 'percentile': - percentile = rank - if percentile < 0.0: - percentile += 100.0 - if percentile < 0 or percentile > 100: - raise RuntimeError('invalid percentile') - if percentile == 100.0: - rank = filter_size - 1 - else: - rank = int(float(filter_size) * percentile / 100.0) - if rank < 0: - rank += filter_size - if rank < 0 or rank >= filter_size: - raise RuntimeError('rank not within filter footprint size') - if rank == 0: - return minimum_filter(input, None, footprint, output, mode, cval, - origins, axes=None) - elif rank == filter_size - 1: - return maximum_filter(input, None, footprint, output, mode, cval, - origins, axes=None) - else: - output = _ni_support._get_output(output, input) - temp_needed = np.may_share_memory(input, output) - if temp_needed: - # input and output arrays cannot share memory - temp = output - output = _ni_support._get_output(output.dtype, input) - if not isinstance(mode, str) and isinstance(mode, Iterable): - raise RuntimeError( - "A sequence of modes is not supported by non-separable rank " - "filters") - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.rank_filter(input, rank, footprint, output, mode, cval, - origins) - if temp_needed: - temp[...] = output - output = temp - return output - - -def rank_filter(input, rank, size=None, footprint=None, output=None, - mode="reflect", cval=0.0, origin=0, *, axes=None): - """Calculate a multidimensional rank filter. - - Parameters - ---------- - %(input)s - rank : int - The rank parameter may be less than zero, i.e., rank = -1 - indicates the largest element. - %(size_foot)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. - - Returns - ------- - rank_filter : ndarray - Filtered array. Has the same shape as `input`. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.rank_filter(ascent, rank=42, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - rank = operator.index(rank) - return _rank_filter(input, rank, size, footprint, output, mode, cval, - origin, 'rank', axes=axes) - - -def median_filter(input, size=None, footprint=None, output=None, - mode="reflect", cval=0.0, origin=0, *, axes=None): - """ - Calculate a multidimensional median filter. - - Parameters - ---------- - %(input)s - %(size_foot)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. - - Returns - ------- - median_filter : ndarray - Filtered array. Has the same shape as `input`. - - See Also - -------- - scipy.signal.medfilt2d - - Notes - ----- - For 2-dimensional images with ``uint8``, ``float32`` or ``float64`` dtypes - the specialised function `scipy.signal.medfilt2d` may be faster. It is - however limited to constant mode with ``cval=0``. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.median_filter(ascent, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - return _rank_filter(input, 0, size, footprint, output, mode, cval, - origin, 'median', axes=axes) - - -def percentile_filter(input, percentile, size=None, footprint=None, - output=None, mode="reflect", cval=0.0, origin=0, *, - axes=None): - """Calculate a multidimensional percentile filter. - - Parameters - ---------- - %(input)s - percentile : scalar - The percentile parameter may be less than zero, i.e., - percentile = -20 equals percentile = 80 - %(size_foot)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin_multiple)s - axes : tuple of int or None, optional - If None, `input` is filtered along all axes. Otherwise, - `input` is filtered along the specified axes. - - Returns - ------- - percentile_filter : ndarray - Filtered array. Has the same shape as `input`. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - >>> fig = plt.figure() - >>> plt.gray() # show the filtered result in grayscale - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.percentile_filter(ascent, percentile=20, size=20) - >>> ax1.imshow(ascent) - >>> ax2.imshow(result) - >>> plt.show() - """ - return _rank_filter(input, percentile, size, footprint, output, mode, - cval, origin, 'percentile', axes=axes) - - -def generic_filter1d(input, function, filter_size, axis=-1, - output=None, mode="reflect", cval=0.0, origin=0, - extra_arguments=(), extra_keywords=None): - """Calculate a 1-D filter along the given axis. - - `generic_filter1d` iterates over the lines of the array, calling the - given function at each line. The arguments of the line are the - input line, and the output line. The input and output lines are 1-D - double arrays. The input line is extended appropriately according - to the filter size and origin. The output line must be modified - in-place with the result. - - Parameters - ---------- - %(input)s - function : {callable, scipy.LowLevelCallable} - Function to apply along given axis. - filter_size : scalar - Length of the filter. - %(axis)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin)s - %(extra_arguments)s - %(extra_keywords)s - - Returns - ------- - generic_filter1d : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - This function also accepts low-level callback functions with one of - the following signatures and wrapped in `scipy.LowLevelCallable`: - - .. code:: c - - int function(double *input_line, npy_intp input_length, - double *output_line, npy_intp output_length, - void *user_data) - int function(double *input_line, intptr_t input_length, - double *output_line, intptr_t output_length, - void *user_data) - - The calling function iterates over the lines of the input and output - arrays, calling the callback function at each line. The current line - is extended according to the border conditions set by the calling - function, and the result is copied into the array that is passed - through ``input_line``. The length of the input line (after extension) - is passed through ``input_length``. The callback function should apply - the filter and store the result in the array passed through - ``output_line``. The length of the output line is passed through - ``output_length``. ``user_data`` is the data pointer provided - to `scipy.LowLevelCallable` as-is. - - The callback function must return an integer error status that is zero - if something went wrong and one otherwise. If an error occurs, you should - normally set the python error status with an informative message - before returning, otherwise a default error message is set by the - calling function. - - In addition, some other low-level function pointer specifications - are accepted, but these are for backward compatibility only and should - not be used in new code. - - """ - if extra_keywords is None: - extra_keywords = {} - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError('Complex type not supported') - output = _ni_support._get_output(output, input) - if filter_size < 1: - raise RuntimeError('invalid filter size') - axis = normalize_axis_index(axis, input.ndim) - if (filter_size // 2 + origin < 0) or (filter_size // 2 + origin >= - filter_size): - raise ValueError('invalid origin') - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.generic_filter1d(input, function, filter_size, axis, output, - mode, cval, origin, extra_arguments, - extra_keywords) - return output - - -def generic_filter(input, function, size=None, footprint=None, - output=None, mode="reflect", cval=0.0, origin=0, - extra_arguments=(), extra_keywords=None): - """Calculate a multidimensional filter using the given function. - - At each element the provided function is called. The input values - within the filter footprint at that element are passed to the function - as a 1-D array of double values. - - Parameters - ---------- - %(input)s - function : {callable, scipy.LowLevelCallable} - Function to apply at each element. - %(size_foot)s - %(output)s - %(mode_reflect)s - %(cval)s - %(origin_multiple)s - %(extra_arguments)s - %(extra_keywords)s - - Returns - ------- - generic_filter : ndarray - Filtered array. Has the same shape as `input`. - - Notes - ----- - This function also accepts low-level callback functions with one of - the following signatures and wrapped in `scipy.LowLevelCallable`: - - .. code:: c - - int callback(double *buffer, npy_intp filter_size, - double *return_value, void *user_data) - int callback(double *buffer, intptr_t filter_size, - double *return_value, void *user_data) - - The calling function iterates over the elements of the input and - output arrays, calling the callback function at each element. The - elements within the footprint of the filter at the current element are - passed through the ``buffer`` parameter, and the number of elements - within the footprint through ``filter_size``. The calculated value is - returned in ``return_value``. ``user_data`` is the data pointer provided - to `scipy.LowLevelCallable` as-is. - - The callback function must return an integer error status that is zero - if something went wrong and one otherwise. If an error occurs, you should - normally set the python error status with an informative message - before returning, otherwise a default error message is set by the - calling function. - - In addition, some other low-level function pointer specifications - are accepted, but these are for backward compatibility only and should - not be used in new code. - - Examples - -------- - Import the necessary modules and load the example image used for - filtering. - - >>> import numpy as np - >>> from scipy import datasets - >>> from scipy.ndimage import zoom, generic_filter - >>> import matplotlib.pyplot as plt - >>> ascent = zoom(datasets.ascent(), 0.5) - - Compute a maximum filter with kernel size 5 by passing a simple NumPy - aggregation function as argument to `function`. - - >>> maximum_filter_result = generic_filter(ascent, np.amax, [5, 5]) - - While a maximmum filter could also directly be obtained using - `maximum_filter`, `generic_filter` allows generic Python function or - `scipy.LowLevelCallable` to be used as a filter. Here, we compute the - range between maximum and minimum value as an example for a kernel size - of 5. - - >>> def custom_filter(image): - ... return np.amax(image) - np.amin(image) - >>> custom_filter_result = generic_filter(ascent, custom_filter, [5, 5]) - - Plot the original and filtered images. - - >>> fig, axes = plt.subplots(3, 1, figsize=(3, 9)) - >>> plt.gray() # show the filtered result in grayscale - >>> top, middle, bottom = axes - >>> for ax in axes: - ... ax.set_axis_off() # remove coordinate system - >>> top.imshow(ascent) - >>> top.set_title("Original image") - >>> middle.imshow(maximum_filter_result) - >>> middle.set_title("Maximum filter, Kernel: 5x5") - >>> bottom.imshow(custom_filter_result) - >>> bottom.set_title("Custom filter, Kernel: 5x5") - >>> fig.tight_layout() - - """ - if (size is not None) and (footprint is not None): - warnings.warn("ignoring size because footprint is set", - UserWarning, stacklevel=2) - if extra_keywords is None: - extra_keywords = {} - input = np.asarray(input) - if np.iscomplexobj(input): - raise TypeError('Complex type not supported') - origins = _ni_support._normalize_sequence(origin, input.ndim) - if footprint is None: - if size is None: - raise RuntimeError("no footprint or filter size provided") - sizes = _ni_support._normalize_sequence(size, input.ndim) - footprint = np.ones(sizes, dtype=bool) - else: - footprint = np.asarray(footprint, dtype=bool) - fshape = [ii for ii in footprint.shape if ii > 0] - if len(fshape) != input.ndim: - raise RuntimeError('filter footprint array has incorrect shape.') - for origin, lenf in zip(origins, fshape): - if (lenf // 2 + origin < 0) or (lenf // 2 + origin >= lenf): - raise ValueError('invalid origin') - if not footprint.flags.contiguous: - footprint = footprint.copy() - output = _ni_support._get_output(output, input) - mode = _ni_support._extend_mode_to_code(mode) - _nd_image.generic_filter(input, function, footprint, output, mode, - cval, origins, extra_arguments, extra_keywords) - return output diff --git a/src/bioxelnodes/bioxel/scipy/_interpolation.py b/src/bioxelnodes/bioxel/scipy/_interpolation.py deleted file mode 100644 index 0479344..0000000 --- a/src/bioxelnodes/bioxel/scipy/_interpolation.py +++ /dev/null @@ -1,314 +0,0 @@ -import numpy as np -import warnings - -from . import _ni_support -from . import _nd_image -from ._utils import normalize_axis_index - - -def _prepad_for_spline_filter(input, mode, cval): - if mode in ['nearest', 'grid-constant']: - npad = 12 - if mode == 'grid-constant': - padded = np.pad(input, npad, mode='constant', - constant_values=cval) - elif mode == 'nearest': - padded = np.pad(input, npad, mode='edge') - else: - # other modes have exact boundary conditions implemented so - # no prepadding is needed - npad = 0 - padded = input - return padded, npad - - -def spline_filter1d(input, order=3, axis=-1, output=np.float64, - mode='mirror'): - """ - Calculate a 1-D spline filter along the given axis. - - The lines of the array along the given axis are filtered by a - spline filter. The order of the spline must be >= 2 and <= 5. - - Parameters - ---------- - %(input)s - order : int, optional - The order of the spline, default is 3. - axis : int, optional - The axis along which the spline filter is applied. Default is the last - axis. - output : ndarray or dtype, optional - The array in which to place the output, or the dtype of the returned - array. Default is ``numpy.float64``. - %(mode_interp_mirror)s - - Returns - ------- - spline_filter1d : ndarray - The filtered input. - - See Also - -------- - spline_filter : Multidimensional spline filter. - - Notes - ----- - All of the interpolation functions in `ndimage` do spline interpolation of - the input image. If using B-splines of `order > 1`, the input image - values have to be converted to B-spline coefficients first, which is - done by applying this 1-D filter sequentially along all - axes of the input. All functions that require B-spline coefficients - will automatically filter their inputs, a behavior controllable with - the `prefilter` keyword argument. For functions that accept a `mode` - parameter, the result will only be correct if it matches the `mode` - used when filtering. - - For complex-valued `input`, this function processes the real and imaginary - components independently. - - .. versionadded:: 1.6.0 - Complex-valued support added. - - Examples - -------- - We can filter an image using 1-D spline along the given axis: - - >>> from scipy.ndimage import spline_filter1d - >>> import numpy as np - >>> import matplotlib.pyplot as plt - >>> orig_img = np.eye(20) # create an image - >>> orig_img[10, :] = 1.0 - >>> sp_filter_axis_0 = spline_filter1d(orig_img, axis=0) - >>> sp_filter_axis_1 = spline_filter1d(orig_img, axis=1) - >>> f, ax = plt.subplots(1, 3, sharex=True) - >>> for ind, data in enumerate([[orig_img, "original image"], - ... [sp_filter_axis_0, "spline filter (axis=0)"], - ... [sp_filter_axis_1, "spline filter (axis=1)"]]): - ... ax[ind].imshow(data[0], cmap='gray_r') - ... ax[ind].set_title(data[1]) - >>> plt.tight_layout() - >>> plt.show() - - """ - if order < 0 or order > 5: - raise RuntimeError('spline order not supported') - input = np.asarray(input) - complex_output = np.iscomplexobj(input) - output = _ni_support._get_output(output, input, - complex_output=complex_output) - if complex_output: - spline_filter1d(input.real, order, axis, output.real, mode) - spline_filter1d(input.imag, order, axis, output.imag, mode) - return output - if order in [0, 1]: - output[...] = np.array(input) - else: - mode = _ni_support._extend_mode_to_code(mode) - axis = normalize_axis_index(axis, input.ndim) - _nd_image.spline_filter1d(input, order, axis, output, mode) - return output - - -def spline_filter(input, order=3, output=np.float64, mode='mirror'): - """ - Multidimensional spline filter. - - Parameters - ---------- - %(input)s - order : int, optional - The order of the spline, default is 3. - output : ndarray or dtype, optional - The array in which to place the output, or the dtype of the returned - array. Default is ``numpy.float64``. - %(mode_interp_mirror)s - - Returns - ------- - spline_filter : ndarray - Filtered array. Has the same shape as `input`. - - See Also - -------- - spline_filter1d : Calculate a 1-D spline filter along the given axis. - - Notes - ----- - The multidimensional filter is implemented as a sequence of - 1-D spline filters. The intermediate arrays are stored - in the same data type as the output. Therefore, for output types - with a limited precision, the results may be imprecise because - intermediate results may be stored with insufficient precision. - - For complex-valued `input`, this function processes the real and imaginary - components independently. - - .. versionadded:: 1.6.0 - Complex-valued support added. - - Examples - -------- - We can filter an image using multidimentional splines: - - >>> from scipy.ndimage import spline_filter - >>> import numpy as np - >>> import matplotlib.pyplot as plt - >>> orig_img = np.eye(20) # create an image - >>> orig_img[10, :] = 1.0 - >>> sp_filter = spline_filter(orig_img, order=3) - >>> f, ax = plt.subplots(1, 2, sharex=True) - >>> for ind, data in enumerate([[orig_img, "original image"], - ... [sp_filter, "spline filter"]]): - ... ax[ind].imshow(data[0], cmap='gray_r') - ... ax[ind].set_title(data[1]) - >>> plt.tight_layout() - >>> plt.show() - - """ - if order < 2 or order > 5: - raise RuntimeError('spline order not supported') - input = np.asarray(input) - complex_output = np.iscomplexobj(input) - output = _ni_support._get_output(output, input, - complex_output=complex_output) - if complex_output: - spline_filter(input.real, order, output.real, mode) - spline_filter(input.imag, order, output.imag, mode) - return output - if order not in [0, 1] and input.ndim > 0: - for axis in range(input.ndim): - spline_filter1d(input, order, axis, output=output, mode=mode) - input = output - else: - output[...] = input[...] - return output - - -def zoom(input, zoom, output=None, order=3, mode='constant', cval=0.0, - prefilter=True, *, grid_mode=False): - """ - Zoom an array. - - The array is zoomed using spline interpolation of the requested order. - - Parameters - ---------- - %(input)s - zoom : float or sequence - The zoom factor along the axes. If a float, `zoom` is the same for each - axis. If a sequence, `zoom` should contain one value for each axis. - %(output)s - order : int, optional - The order of the spline interpolation, default is 3. - The order has to be in the range 0-5. - %(mode_interp_constant)s - %(cval)s - %(prefilter)s - grid_mode : bool, optional - If False, the distance from the pixel centers is zoomed. Otherwise, the - distance including the full pixel extent is used. For example, a 1d - signal of length 5 is considered to have length 4 when `grid_mode` is - False, but length 5 when `grid_mode` is True. See the following - visual illustration: - - .. code-block:: text - - | pixel 1 | pixel 2 | pixel 3 | pixel 4 | pixel 5 | - |<-------------------------------------->| - vs. - |<----------------------------------------------->| - - The starting point of the arrow in the diagram above corresponds to - coordinate location 0 in each mode. - - Returns - ------- - zoom : ndarray - The zoomed input. - - Notes - ----- - For complex-valued `input`, this function zooms the real and imaginary - components independently. - - .. versionadded:: 1.6.0 - Complex-valued support added. - - Examples - -------- - >>> from scipy import ndimage, datasets - >>> import matplotlib.pyplot as plt - - >>> fig = plt.figure() - >>> ax1 = fig.add_subplot(121) # left side - >>> ax2 = fig.add_subplot(122) # right side - >>> ascent = datasets.ascent() - >>> result = ndimage.zoom(ascent, 3.0) - >>> ax1.imshow(ascent, vmin=0, vmax=255) - >>> ax2.imshow(result, vmin=0, vmax=255) - >>> plt.show() - - >>> print(ascent.shape) - (512, 512) - - >>> print(result.shape) - (1536, 1536) - """ - if order < 0 or order > 5: - raise RuntimeError('spline order not supported') - input = np.asarray(input) - if input.ndim < 1: - raise RuntimeError('input and output rank must be > 0') - zoom = _ni_support._normalize_sequence(zoom, input.ndim) - output_shape = tuple( - [int(round(ii * jj)) for ii, jj in zip(input.shape, zoom)]) - complex_output = np.iscomplexobj(input) - output = _ni_support._get_output(output, input, shape=output_shape, - complex_output=complex_output) - if complex_output: - # import under different name to avoid confusion with zoom parameter - from scipy.ndimage._interpolation import zoom as _zoom - - kwargs = dict(order=order, mode=mode, prefilter=prefilter) - _zoom(input.real, zoom, output=output.real, - cval=np.real(cval), **kwargs) - _zoom(input.imag, zoom, output=output.imag, - cval=np.imag(cval), **kwargs) - return output - if prefilter and order > 1: - padded, npad = _prepad_for_spline_filter(input, mode, cval) - filtered = spline_filter(padded, order, output=np.float64, mode=mode) - else: - npad = 0 - filtered = input - if grid_mode: - # warn about modes that may have surprising behavior - suggest_mode = None - if mode == 'constant': - suggest_mode = 'grid-constant' - elif mode == 'wrap': - suggest_mode = 'grid-wrap' - if suggest_mode is not None: - warnings.warn( - (f"It is recommended to use mode = {suggest_mode} instead of {mode} " - f"when grid_mode is True."), - stacklevel=2 - ) - mode = _ni_support._extend_mode_to_code(mode) - - zoom_div = np.array(output_shape) - zoom_nominator = np.array(input.shape) - if not grid_mode: - zoom_div -= 1 - zoom_nominator -= 1 - - # Zooming to infinite values is unpredictable, so just choose - # zoom factor 1 instead - zoom = np.divide(zoom_nominator, zoom_div, - out=np.ones_like(input.shape, dtype=np.float64), - where=zoom_div != 0) - zoom = np.ascontiguousarray(zoom) - _nd_image.zoom_shift(filtered, zoom, None, output, order, mode, cval, npad, - grid_mode) - return output diff --git a/src/bioxelnodes/bioxel/scipy/_ni_support.py b/src/bioxelnodes/bioxel/scipy/_ni_support.py deleted file mode 100644 index c312772..0000000 --- a/src/bioxelnodes/bioxel/scipy/_ni_support.py +++ /dev/null @@ -1,120 +0,0 @@ -# Copyright (C) 2003-2005 Peter J. Verveer -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# 1. Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# -# 2. Redistributions in binary form must reproduce the above -# copyright notice, this list of conditions and the following -# disclaimer in the documentation and/or other materials provided -# with the distribution. -# -# 3. The name of the author may not be used to endorse or promote -# products derived from this software without specific prior -# written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS -# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY -# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE -# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, -# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -from collections.abc import Iterable -import operator -import warnings -import numpy as np - - -def _extend_mode_to_code(mode): - """Convert an extension mode to the corresponding integer code. - """ - if mode == 'nearest': - return 0 - elif mode == 'wrap': - return 1 - elif mode in ['reflect', 'grid-mirror']: - return 2 - elif mode == 'mirror': - return 3 - elif mode == 'constant': - return 4 - elif mode == 'grid-wrap': - return 5 - elif mode == 'grid-constant': - return 6 - else: - raise RuntimeError('boundary mode not supported') - - -def _normalize_sequence(input, rank): - """If input is a scalar, create a sequence of length equal to the - rank by duplicating the input. If input is a sequence, - check if its length is equal to the length of array. - """ - is_str = isinstance(input, str) - if not is_str and isinstance(input, Iterable): - normalized = list(input) - if len(normalized) != rank: - err = "sequence argument must have length equal to input rank" - raise RuntimeError(err) - else: - normalized = [input] * rank - return normalized - - -def _get_output(output, input, shape=None, complex_output=False): - if shape is None: - shape = input.shape - if output is None: - if not complex_output: - output = np.zeros(shape, dtype=input.dtype.name) - else: - complex_type = np.promote_types(input.dtype, np.complex64) - output = np.zeros(shape, dtype=complex_type) - elif isinstance(output, (type, np.dtype)): - # Classes (like `np.float32`) and dtypes are interpreted as dtype - if complex_output and np.dtype(output).kind != 'c': - warnings.warn( - "promoting specified output dtype to complex", stacklevel=3) - output = np.promote_types(output, np.complex64) - output = np.zeros(shape, dtype=output) - elif isinstance(output, str): - output = np.dtype(output) - if complex_output and output.kind != 'c': - raise RuntimeError("output must have complex dtype") - elif not issubclass(output.type, np.number): - raise RuntimeError("output must have numeric dtype") - output = np.zeros(shape, dtype=output) - elif output.shape != shape: - raise RuntimeError("output shape not correct") - elif complex_output and output.dtype.kind != 'c': - raise RuntimeError("output must have complex dtype") - return output - - -def _check_axes(axes, ndim): - if axes is None: - return tuple(range(ndim)) - elif np.isscalar(axes): - axes = (operator.index(axes),) - elif isinstance(axes, Iterable): - for ax in axes: - axes = tuple(operator.index(ax) for ax in axes) - if ax < -ndim or ax > ndim - 1: - raise ValueError(f"specified axis: {ax} is out of range") - axes = tuple(ax % ndim if ax < 0 else ax for ax in axes) - else: - message = "axes must be an integer, iterable of integers, or None" - raise ValueError(message) - if len(tuple(set(axes))) != len(axes): - raise ValueError("axes must be unique") - return axes diff --git a/src/bioxelnodes/bioxel/scipy/_utils.py b/src/bioxelnodes/bioxel/scipy/_utils.py deleted file mode 100644 index 5a4b9be..0000000 --- a/src/bioxelnodes/bioxel/scipy/_utils.py +++ /dev/null @@ -1,10 +0,0 @@ - -def normalize_axis_index(axis, ndim): - # Check if `axis` is in the correct range and normalize it - if axis < -ndim or axis >= ndim: - msg = f"axis {axis} is out of bounds for array of dimension {ndim}" - raise Exception(msg) - - if axis < 0: - axis = axis + ndim - return axis diff --git a/src/bioxelnodes/bioxel/skimage/__init__.py b/src/bioxelnodes/bioxel/skimage/__init__.py deleted file mode 100644 index 79b5d74..0000000 --- a/src/bioxelnodes/bioxel/skimage/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from ._warps import resize diff --git a/src/bioxelnodes/bioxel/skimage/_utils.py b/src/bioxelnodes/bioxel/skimage/_utils.py deleted file mode 100644 index f07a428..0000000 --- a/src/bioxelnodes/bioxel/skimage/_utils.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np - - -def convert_to_float(image, preserve_range): - """Convert input image to float image with the appropriate range. - - Parameters - ---------- - image : ndarray - Input image. - preserve_range : bool - Determines if the range of the image should be kept or transformed - using img_as_float. Also see - https://scikit-image.org/docs/dev/user_guide/data_types.html - - Notes - ----- - * Input images with `float32` data type are not upcast. - - Returns - ------- - image : ndarray - Transformed version of the input. - - """ - if image.dtype == np.float16: - return image.astype(np.float32) - if preserve_range: - # Convert image to double only if it is not single or double - # precision float - if image.dtype.char not in 'df': - image = image.astype(float) - else: - from .dtype import img_as_float - - image = img_as_float(image) - return image diff --git a/src/bioxelnodes/bioxel/skimage/_warps.py b/src/bioxelnodes/bioxel/skimage/_warps.py deleted file mode 100644 index f8f4a4f..0000000 --- a/src/bioxelnodes/bioxel/skimage/_warps.py +++ /dev/null @@ -1,319 +0,0 @@ -import warnings -import numpy as np -from .. import scipy as ndi -from ._utils import convert_to_float - - -def _clip_warp_output(input_image, output_image, mode, cval, clip): - """Clip output image to range of values of input image. - - Note that this function modifies the values of `output_image` in-place - and it is only modified if ``clip=True``. - - Parameters - ---------- - input_image : ndarray - Input image. - output_image : ndarray - Output image, which is modified in-place. - - Other parameters - ---------------- - mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'} - Points outside the boundaries of the input are filled according - to the given mode. Modes match the behaviour of `numpy.pad`. - cval : float - Used in conjunction with mode 'constant', the value outside - the image boundaries. - clip : bool - Whether to clip the output to the range of values of the input image. - This is enabled by default, since higher order interpolation may - produce values outside the given input range. - - """ - if clip: - min_val = np.min(input_image) - if np.isnan(min_val): - # NaNs detected, use NaN-safe min/max - min_func = np.nanmin - max_func = np.nanmax - min_val = min_func(input_image) - else: - min_func = np.min - max_func = np.max - max_val = max_func(input_image) - - # Check if cval has been used such that it expands the effective input - # range - preserve_cval = ( - mode == 'constant' - and not min_val <= cval <= max_val - and min_func(output_image) <= cval <= max_func(output_image) - ) - - # expand min/max range to account for cval - if preserve_cval: - # cast cval to the same dtype as the input image - cval = input_image.dtype.type(cval) - min_val = min(min_val, cval) - max_val = max(max_val, cval) - - # Convert array-like types to ndarrays (gh-7159) - min_val, max_val = np.asarray(min_val), np.asarray(max_val) - np.clip(output_image, min_val, max_val, out=output_image) - - -def _validate_interpolation_order(image_dtype, order): - """Validate and return spline interpolation's order. - - Parameters - ---------- - image_dtype : dtype - Image dtype. - order : int, optional - The order of the spline interpolation. The order has to be in - the range 0-5. See `skimage.transform.warp` for detail. - - Returns - ------- - order : int - if input order is None, returns 0 if image_dtype is bool and 1 - otherwise. Otherwise, image_dtype is checked and input order - is validated accordingly (order > 0 is not supported for bool - image dtype) - - """ - - if order is None: - return 0 if image_dtype == bool else 1 - - if order < 0 or order > 5: - raise ValueError( - "Spline interpolation order has to be in the " "range 0-5.") - - if image_dtype == bool and order != 0: - raise ValueError( - "Input image dtype is bool. Interpolation is not defined " - "with bool data type. Please set order to 0 or explicitly " - "cast input image to another data type." - ) - - return order - - -def _preprocess_resize_output_shape(image, output_shape): - """Validate resize output shape according to input image. - - Parameters - ---------- - image: ndarray - Image to be resized. - output_shape: iterable - Size of the generated output image `(rows, cols[, ...][, dim])`. If - `dim` is not provided, the number of channels is preserved. - - Returns - ------- - image: ndarray - The input image, but with additional singleton dimensions appended in - the case where ``len(output_shape) > input.ndim``. - output_shape: tuple - The output image converted to tuple. - - Raises - ------ - ValueError: - If output_shape length is smaller than the image number of - dimensions - - Notes - ----- - The input image is reshaped if its number of dimensions is not - equal to output_shape_length. - - """ - output_shape = tuple(output_shape) - output_ndim = len(output_shape) - input_shape = image.shape - if output_ndim > image.ndim: - # append dimensions to input_shape - input_shape += (1,) * (output_ndim - image.ndim) - image = np.reshape(image, input_shape) - elif output_ndim == image.ndim - 1: - # multichannel case: append shape of last axis - output_shape = output_shape + (image.shape[-1],) - elif output_ndim < image.ndim: - raise ValueError( - "output_shape length cannot be smaller than the " - "image number of dimensions" - ) - - return image, output_shape - - -def _to_ndimage_mode(mode): - """Convert from `numpy.pad` mode name to the corresponding ndimage mode.""" - mode_translation_dict = dict( - constant='constant', - edge='nearest', - symmetric='reflect', - reflect='mirror', - wrap='wrap', - ) - if mode not in mode_translation_dict: - raise ValueError( - f"Unknown mode: '{mode}', or cannot translate mode. The " - f"mode should be one of 'constant', 'edge', 'symmetric', " - f"'reflect', or 'wrap'. See the documentation of numpy.pad for " - f"more info." - ) - return _fix_ndimage_mode(mode_translation_dict[mode]) - - -def _fix_ndimage_mode(mode): - # SciPy 1.6.0 introduced grid variants of constant and wrap which - # have less surprising behavior for images. Use these when available - grid_modes = {'constant': 'grid-constant', 'wrap': 'grid-wrap'} - return grid_modes.get(mode, mode) - - -def resize( - image, - output_shape, - order=None, - mode='reflect', - cval=0, - clip=True, - preserve_range=False, - anti_aliasing=None, - anti_aliasing_sigma=None, -): - """Resize image to match a certain size. - - Performs interpolation to up-size or down-size N-dimensional images. Note - that anti-aliasing should be enabled when down-sizing images to avoid - aliasing artifacts. For downsampling with an integer factor also see - `skimage.transform.downscale_local_mean`. - - Parameters - ---------- - image : ndarray - Input image. - output_shape : iterable - Size of the generated output image `(rows, cols[, ...][, dim])`. If - `dim` is not provided, the number of channels is preserved. In case the - number of input channels does not equal the number of output channels a - n-dimensional interpolation is applied. - - Returns - ------- - resized : ndarray - Resized version of the input. - - Other parameters - ---------------- - order : int, optional - The order of the spline interpolation, default is 0 if - image.dtype is bool and 1 otherwise. The order has to be in - the range 0-5. See `skimage.transform.warp` for detail. - mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional - Points outside the boundaries of the input are filled according - to the given mode. Modes match the behaviour of `numpy.pad`. - cval : float, optional - Used in conjunction with mode 'constant', the value outside - the image boundaries. - clip : bool, optional - Whether to clip the output to the range of values of the input image. - This is enabled by default, since higher order interpolation may - produce values outside the given input range. - preserve_range : bool, optional - Whether to keep the original range of values. Otherwise, the input - image is converted according to the conventions of `img_as_float`. - Also see https://scikit-image.org/docs/dev/user_guide/data_types.html - anti_aliasing : bool, optional - Whether to apply a Gaussian filter to smooth the image prior - to downsampling. It is crucial to filter when downsampling - the image to avoid aliasing artifacts. If not specified, it is set to - True when downsampling an image whose data type is not bool. - It is also set to False when using nearest neighbor interpolation - (``order`` == 0) with integer input data type. - anti_aliasing_sigma : {float, tuple of floats}, optional - Standard deviation for Gaussian filtering used when anti-aliasing. - By default, this value is chosen as (s - 1) / 2 where s is the - downsampling factor, where s > 1. For the up-size case, s < 1, no - anti-aliasing is performed prior to rescaling. - - Notes - ----- - Modes 'reflect' and 'symmetric' are similar, but differ in whether the edge - pixels are duplicated during the reflection. As an example, if an array - has values [0, 1, 2] and was padded to the right by four values using - symmetric, the result would be [0, 1, 2, 2, 1, 0, 0], while for reflect it - would be [0, 1, 2, 1, 0, 1, 2]. - - Examples - -------- - >>> from skimage import data - >>> from skimage.transform import resize - >>> image = data.camera() - >>> resize(image, (100, 100)).shape - (100, 100) - - """ - - image, output_shape = _preprocess_resize_output_shape(image, output_shape) - input_shape = image.shape - input_type = image.dtype - - if input_type == np.float16: - image = image.astype(np.float32) - - if anti_aliasing is None: - anti_aliasing = ( - not input_type == bool - and not (np.issubdtype(input_type, np.integer) and order == 0) - and any(x < y for x, y in zip(output_shape, input_shape)) - ) - - if input_type == bool and anti_aliasing: - raise ValueError("anti_aliasing must be False for boolean images") - - factors = np.divide(input_shape, output_shape) - order = _validate_interpolation_order(input_type, order) - if order > 0: - image = convert_to_float(image, preserve_range) - - # Translate modes used by np.pad to those used by scipy.ndimage - ndi_mode = _to_ndimage_mode(mode) - if anti_aliasing: - if anti_aliasing_sigma is None: - anti_aliasing_sigma = np.maximum(0, (factors - 1) / 2) - else: - anti_aliasing_sigma = np.atleast_1d(anti_aliasing_sigma) * np.ones_like( - factors - ) - if np.any(anti_aliasing_sigma < 0): - raise ValueError( - "Anti-aliasing standard deviation must be " - "greater than or equal to zero" - ) - elif np.any((anti_aliasing_sigma > 0) & (factors <= 1)): - warnings.warn( - "Anti-aliasing standard deviation greater than zero but " - "not down-sampling along all axes", - stacklevel=2 - ) - filtered = ndi.gaussian_filter( - image, anti_aliasing_sigma, cval=cval, mode=ndi_mode - ) - else: - filtered = image - - zoom_factors = [1 / f for f in factors] - out = ndi.zoom( - filtered, zoom_factors, order=order, mode=ndi_mode, cval=cval, grid_mode=True - ) - - _clip_warp_output(image, out, mode, cval, clip) - - return out diff --git a/src/bioxelnodes/bioxel/skimage/dtype.py b/src/bioxelnodes/bioxel/skimage/dtype.py deleted file mode 100644 index 07bc3b4..0000000 --- a/src/bioxelnodes/bioxel/skimage/dtype.py +++ /dev/null @@ -1,605 +0,0 @@ -import warnings -from warnings import warn - -import numpy as np - - -__all__ = [ - 'img_as_float32', - 'img_as_float64', - 'img_as_float', - 'img_as_int', - 'img_as_uint', - 'img_as_ubyte', - 'img_as_bool', - 'dtype_limits', -] - -# Some of these may or may not be aliases depending on architecture & platform -_integer_types = ( - np.int8, - np.byte, - np.int16, - np.short, - np.int32, - np.int64, - np.longlong, - np.int_, - np.intp, - np.intc, - int, - np.uint8, - np.ubyte, - np.uint16, - np.ushort, - np.uint32, - np.uint64, - np.ulonglong, - np.uint, - np.uintp, - np.uintc, -) -_integer_ranges = {t: (np.iinfo(t).min, np.iinfo(t).max) - for t in _integer_types} -dtype_range = { - bool: (False, True), - np.bool_: (False, True), - float: (-1, 1), - np.float16: (-1, 1), - np.float32: (-1, 1), - np.float64: (-1, 1), -} - -with warnings.catch_warnings(): - warnings.filterwarnings('ignore', category=DeprecationWarning) - - # np.bool8 is a deprecated alias of np.bool_ - if hasattr(np, 'bool8'): - dtype_range[np.bool8] = (False, True) - -dtype_range.update(_integer_ranges) - -_supported_types = list(dtype_range.keys()) - - -def dtype_limits(image, clip_negative=False): - """Return intensity limits, i.e. (min, max) tuple, of the image's dtype. - - Parameters - ---------- - image : ndarray - Input image. - clip_negative : bool, optional - If True, clip the negative range (i.e. return 0 for min intensity) - even if the image dtype allows negative values. - - Returns - ------- - imin, imax : tuple - Lower and upper intensity limits. - """ - imin, imax = dtype_range[image.dtype.type] - if clip_negative: - imin = 0 - return imin, imax - - -def _dtype_itemsize(itemsize, *dtypes): - """Return first of `dtypes` with itemsize greater than `itemsize` - - Parameters - ---------- - itemsize: int - The data type object element size. - - Other Parameters - ---------------- - *dtypes: - Any Object accepted by `np.dtype` to be converted to a data - type object - - Returns - ------- - dtype: data type object - First of `dtypes` with itemsize greater than `itemsize`. - - """ - return next(dt for dt in dtypes if np.dtype(dt).itemsize >= itemsize) - - -def _dtype_bits(kind, bits, itemsize=1): - """Return dtype of `kind` that can store a `bits` wide unsigned int - - Parameters: - kind: str - Data type kind. - bits: int - Desired number of bits. - itemsize: int - The data type object element size. - - Returns - ------- - dtype: data type object - Data type of `kind` that can store a `bits` wide unsigned int - - """ - - s = next( - i - for i in (itemsize,) + (2, 4, 8) - if bits < (i * 8) or (bits == (i * 8) and kind == 'u') - ) - - return np.dtype(kind + str(s)) - - -def _scale(a, n, m, copy=True): - """Scale an array of unsigned/positive integers from `n` to `m` bits. - - Numbers can be represented exactly only if `m` is a multiple of `n`. - - Parameters - ---------- - a : ndarray - Input image array. - n : int - Number of bits currently used to encode the values in `a`. - m : int - Desired number of bits to encode the values in `out`. - copy : bool, optional - If True, allocates and returns new array. Otherwise, modifies - `a` in place. - - Returns - ------- - out : array - Output image array. Has the same kind as `a`. - """ - kind = a.dtype.kind - if n > m and a.max() < 2**m: - mnew = int(np.ceil(m / 2) * 2) - if mnew > m: - dtype = f'int{mnew}' - else: - dtype = f'uint{mnew}' - n = int(np.ceil(n / 2) * 2) - warn( - f'Downcasting {a.dtype} to {dtype} without scaling because max ' - f'value {a.max()} fits in {dtype}', - stacklevel=3, - ) - return a.astype(_dtype_bits(kind, m)) - elif n == m: - return a.copy() if copy else a - elif n > m: - # downscale with precision loss - if copy: - b = np.empty(a.shape, _dtype_bits(kind, m)) - np.floor_divide(a, 2 ** (n - m), out=b, - dtype=a.dtype, casting='unsafe') - return b - else: - a //= 2 ** (n - m) - return a - elif m % n == 0: - # exact upscale to a multiple of `n` bits - if copy: - b = np.empty(a.shape, _dtype_bits(kind, m)) - np.multiply(a, (2**m - 1) // (2**n - 1), out=b, dtype=b.dtype) - return b - else: - a = a.astype(_dtype_bits(kind, m, a.dtype.itemsize), copy=False) - a *= (2**m - 1) // (2**n - 1) - return a - else: - # upscale to a multiple of `n` bits, - # then downscale with precision loss - o = (m // n + 1) * n - if copy: - b = np.empty(a.shape, _dtype_bits(kind, o)) - np.multiply(a, (2**o - 1) // (2**n - 1), out=b, dtype=b.dtype) - b //= 2 ** (o - m) - return b - else: - a = a.astype(_dtype_bits(kind, o, a.dtype.itemsize), copy=False) - a *= (2**o - 1) // (2**n - 1) - a //= 2 ** (o - m) - return a - - -def _convert(image, dtype, force_copy=False, uniform=False): - """ - Convert an image to the requested data-type. - - Warnings are issued in case of precision loss, or when negative values - are clipped during conversion to unsigned integer types (sign loss). - - Floating point values are expected to be normalized and will be clipped - to the range [0.0, 1.0] or [-1.0, 1.0] when converting to unsigned or - signed integers respectively. - - Numbers are not shifted to the negative side when converting from - unsigned to signed integer types. Negative values will be clipped when - converting to unsigned integers. - - Parameters - ---------- - image : ndarray - Input image. - dtype : dtype - Target data-type. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - uniform : bool, optional - Uniformly quantize the floating point range to the integer range. - By default (uniform=False) floating point values are scaled and - rounded to the nearest integers, which minimizes back and forth - conversion errors. - - .. versionchanged:: 0.15 - ``_convert`` no longer warns about possible precision or sign - information loss. See discussions on these warnings at: - https://github.com/scikit-image/scikit-image/issues/2602 - https://github.com/scikit-image/scikit-image/issues/543#issuecomment-208202228 - https://github.com/scikit-image/scikit-image/pull/3575 - - References - ---------- - .. [1] DirectX data conversion rules. - https://msdn.microsoft.com/en-us/library/windows/desktop/dd607323%28v=vs.85%29.aspx - .. [2] Data Conversions. In "OpenGL ES 2.0 Specification v2.0.25", - pp 7-8. Khronos Group, 2010. - .. [3] Proper treatment of pixels as integers. A.W. Paeth. - In "Graphics Gems I", pp 249-256. Morgan Kaufmann, 1990. - .. [4] Dirty Pixels. J. Blinn. In "Jim Blinn's corner: Dirty Pixels", - pp 47-57. Morgan Kaufmann, 1998. - - """ - image = np.asarray(image) - dtypeobj_in = image.dtype - if dtype is np.floating: - dtypeobj_out = np.dtype('float64') - else: - dtypeobj_out = np.dtype(dtype) - dtype_in = dtypeobj_in.type - dtype_out = dtypeobj_out.type - kind_in = dtypeobj_in.kind - kind_out = dtypeobj_out.kind - itemsize_in = dtypeobj_in.itemsize - itemsize_out = dtypeobj_out.itemsize - - # Below, we do an `issubdtype` check. Its purpose is to find out - # whether we can get away without doing any image conversion. This happens - # when: - # - # - the output and input dtypes are the same or - # - when the output is specified as a type, and the input dtype - # is a subclass of that type (e.g. `np.floating` will allow - # `float32` and `float64` arrays through) - - if np.issubdtype(dtype_in, dtype): - if force_copy: - image = image.copy() - return image - - if not (dtype_in in _supported_types and dtype_out in _supported_types): - raise ValueError( - f'Cannot convert from {dtypeobj_in} to ' f'{dtypeobj_out}.') - - if kind_in in 'ui': - imin_in = np.iinfo(dtype_in).min - imax_in = np.iinfo(dtype_in).max - if kind_out in 'ui': - imin_out = np.iinfo(dtype_out).min - imax_out = np.iinfo(dtype_out).max - - # any -> binary - if kind_out == 'b': - return image > dtype_in(dtype_range[dtype_in][1] / 2) - - # binary -> any - if kind_in == 'b': - result = image.astype(dtype_out) - if kind_out != 'f': - result *= dtype_out(dtype_range[dtype_out][1]) - return result - - # float -> any - if kind_in == 'f': - if kind_out == 'f': - # float -> float - return image.astype(dtype_out) - - if np.min(image) < -1.0 or np.max(image) > 1.0: - raise ValueError("Images of type float must be between -1 and 1.") - # floating point -> integer - # use float type that can represent output integer type - computation_type = _dtype_itemsize( - itemsize_out, dtype_in, np.float32, np.float64 - ) - - if not uniform: - if kind_out == 'u': - image_out = np.multiply( - image, imax_out, dtype=computation_type) - else: - image_out = np.multiply( - image, (imax_out - imin_out) / 2, dtype=computation_type - ) - image_out -= 1.0 / 2.0 - np.rint(image_out, out=image_out) - np.clip(image_out, imin_out, imax_out, out=image_out) - elif kind_out == 'u': - image_out = np.multiply( - image, imax_out + 1, dtype=computation_type) - np.clip(image_out, 0, imax_out, out=image_out) - else: - image_out = np.multiply( - image, (imax_out - imin_out + 1.0) / 2.0, dtype=computation_type - ) - np.floor(image_out, out=image_out) - np.clip(image_out, imin_out, imax_out, out=image_out) - return image_out.astype(dtype_out) - - # signed/unsigned int -> float - if kind_out == 'f': - # use float type that can exactly represent input integers - computation_type = _dtype_itemsize( - itemsize_in, dtype_out, np.float32, np.float64 - ) - - if kind_in == 'u': - # using np.divide or np.multiply doesn't copy the data - # until the computation time - image = np.multiply(image, 1.0 / imax_in, dtype=computation_type) - # DirectX uses this conversion also for signed ints - # if imin_in: - # np.maximum(image, -1.0, out=image) - elif kind_in == 'i': - # From DirectX conversions: - # The most negative value maps to -1.0f - # Every other value is converted to a float (call it c) - # and then result = c * (1.0f / (2⁽ⁿ⁻¹⁾-1)). - - image = np.multiply(image, 1.0 / imax_in, dtype=computation_type) - np.maximum(image, -1.0, out=image) - - else: - image = np.add(image, 0.5, dtype=computation_type) - image *= 2 / (imax_in - imin_in) - - return np.asarray(image, dtype_out) - - # unsigned int -> signed/unsigned int - if kind_in == 'u': - if kind_out == 'i': - # unsigned int -> signed int - image = _scale(image, 8 * itemsize_in, 8 * itemsize_out - 1) - return image.view(dtype_out) - else: - # unsigned int -> unsigned int - return _scale(image, 8 * itemsize_in, 8 * itemsize_out) - - # signed int -> unsigned int - if kind_out == 'u': - image = _scale(image, 8 * itemsize_in - 1, 8 * itemsize_out) - result = np.empty(image.shape, dtype_out) - np.maximum(image, 0, out=result, dtype=image.dtype, casting='unsafe') - return result - - # signed int -> signed int - if itemsize_in > itemsize_out: - return _scale(image, 8 * itemsize_in - 1, 8 * itemsize_out - 1) - - image = image.astype(_dtype_bits('i', itemsize_out * 8)) - image -= imin_in - image = _scale(image, 8 * itemsize_in, 8 * itemsize_out, copy=False) - image += imin_out - return image.astype(dtype_out) - - -def convert(image, dtype, force_copy=False, uniform=False): - warn( - "The use of this function is discouraged as its behavior may change " - "dramatically in scikit-image 1.0. This function will be removed " - "in scikit-image 1.0.", - FutureWarning, - stacklevel=2, - ) - return _convert(image=image, dtype=dtype, force_copy=force_copy, uniform=uniform) - - -if _convert.__doc__ is not None: - convert.__doc__ = ( - _convert.__doc__ - + """ - - Warns - ----- - FutureWarning: - .. versionadded:: 0.17 - - The use of this function is discouraged as its behavior may change - dramatically in scikit-image 1.0. This function will be removed - in scikit-image 1.0. - """ - ) - - -def img_as_float32(image, force_copy=False): - """Convert an image to single-precision (32-bit) floating point format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of float32 - Output image. - - Notes - ----- - The range of a floating point image is [0.0, 1.0] or [-1.0, 1.0] when - converting from unsigned or signed datatypes, respectively. - If the input image has a float type, intensity values are not modified - and can be outside the ranges [0.0, 1.0] or [-1.0, 1.0]. - - """ - return _convert(image, np.float32, force_copy) - - -def img_as_float64(image, force_copy=False): - """Convert an image to double-precision (64-bit) floating point format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of float64 - Output image. - - Notes - ----- - The range of a floating point image is [0.0, 1.0] or [-1.0, 1.0] when - converting from unsigned or signed datatypes, respectively. - If the input image has a float type, intensity values are not modified - and can be outside the ranges [0.0, 1.0] or [-1.0, 1.0]. - - """ - return _convert(image, np.float64, force_copy) - - -def img_as_float(image, force_copy=False): - """Convert an image to floating point format. - - This function is similar to `img_as_float64`, but will not convert - lower-precision floating point arrays to `float64`. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of float - Output image. - - Notes - ----- - The range of a floating point image is [0.0, 1.0] or [-1.0, 1.0] when - converting from unsigned or signed datatypes, respectively. - If the input image has a float type, intensity values are not modified - and can be outside the ranges [0.0, 1.0] or [-1.0, 1.0]. - - """ - return _convert(image, np.floating, force_copy) - - -def img_as_uint(image, force_copy=False): - """Convert an image to 16-bit unsigned integer format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of uint16 - Output image. - - Notes - ----- - Negative input values will be clipped. - Positive values are scaled between 0 and 65535. - - """ - return _convert(image, np.uint16, force_copy) - - -def img_as_int(image, force_copy=False): - """Convert an image to 16-bit signed integer format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of int16 - Output image. - - Notes - ----- - The values are scaled between -32768 and 32767. - If the input data-type is positive-only (e.g., uint8), then - the output image will still only have positive values. - - """ - return _convert(image, np.int16, force_copy) - - -def img_as_ubyte(image, force_copy=False): - """Convert an image to 8-bit unsigned integer format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of ubyte (uint8) - Output image. - - Notes - ----- - Negative input values will be clipped. - Positive values are scaled between 0 and 255. - - """ - return _convert(image, np.uint8, force_copy) - - -def img_as_bool(image, force_copy=False): - """Convert an image to boolean format. - - Parameters - ---------- - image : ndarray - Input image. - force_copy : bool, optional - Force a copy of the data, irrespective of its current dtype. - - Returns - ------- - out : ndarray of bool (`bool_`) - Output image. - - Notes - ----- - The upper half of the input dtype's positive range is True, and the lower - half is False. All negative values (if present) are False. - - """ - return _convert(image, bool, force_copy)