diff --git a/.travis.yml b/.travis.yml index d3c10f2c..7cc52075 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,8 +18,8 @@ addons: packages: - ccache - libfreeimage3 - - texlive - - texlive-latex-extra + - texlive + - texlive-latex-extra - dvipng - python-qt4 env: diff --git a/CONTRIBUTING.txt b/CONTRIBUTING.txt index 29dba55d..38bb784d 100644 --- a/CONTRIBUTING.txt +++ b/CONTRIBUTING.txt @@ -77,6 +77,8 @@ For a more detailed discussion, read these :doc:`detailed documents Travis fails, you can find out why by clicking on the "failed" icon (red cross) and inspecting the build and test log. + * A pull request must be approved by two core team members before merging. + 5. Document changes If your change introduces any API modifications, please update @@ -127,7 +129,8 @@ Guidelines * All code should be documented, to the same `standard `_ as NumPy and SciPy. * For new functionality, always add an example to the gallery. -* No changes are ever committed without review. Ask on the +* No changes are ever committed without review and approval by two core + team members. Ask on the `mailing list `_ if you get no response to your pull request. **Never merge your own pull request.** diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index c94755d8..9e028354 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -218,3 +218,6 @@ - Damian Eads Structuring elements in morphology module. + +- Egor Panfilov + Inpainting with biharmonic equation diff --git a/TODO.txt b/TODO.txt index bed2793a..fcb8952f 100644 --- a/TODO.txt +++ b/TODO.txt @@ -13,6 +13,7 @@ Version 0.14 parameters `(dist, theta)`, LineModelND has the more general parameters `(origin, direction)`. * Remove deprecated old syntax support for ``skimage.transform.integrate``. +* Remove deprecated ``skimage.data.lena`` and corresponding data files. Version 0.13 @@ -34,7 +35,6 @@ Version 0.13 _shared/interpolation.pyx, transform/_geometric.py, and transform/_warps.py - Version 0.12 ------------ * Change `label` to mark background as 0, not -1, which is consistent with diff --git a/doc/examples/filters/plot_inpaint.py b/doc/examples/filters/plot_inpaint.py new file mode 100644 index 00000000..104694ec --- /dev/null +++ b/doc/examples/filters/plot_inpaint.py @@ -0,0 +1,58 @@ +""" +=========== +Inpainting +=========== +Inpainting [1]_ is the process of reconstructing lost or deteriorated +parts of images and videos. + +The reconstruction is supposed to be performed in fully automatic way by +exploiting the information presented in non-damaged regions. + +In this example, we show how the masked pixels get inpainted by +inpainting algorithm based on 'biharmonic equation'-assumption [2]_ [3]_. + +.. [1] Wikipedia. Inpainting + https://en.wikipedia.org/wiki/Inpainting +.. [2] Wikipedia. Biharmonic equation + https://en.wikipedia.org/wiki/Biharmonic_equation +.. [3] N.S.Hoang, S.B.Damelin, "On surface completion and image + inpainting by biharmonic functions: numerical aspects", + http://www.ima.umn.edu/~damelin/biharmonic +""" + +import numpy as np +import matplotlib.pyplot as plt + +from skimage import data, color +from skimage.restoration import inpaint + +image_orig = data.astronaut() + +# Create mask with three defect regions: left, middle, right respectively +mask = np.zeros(image_orig.shape[:-1]) +mask[20:60, 0:20] = 1 +mask[200:300, 150:170] = 1 +mask[50:100, 400:430] = 1 + +# Defect image over the same region in each color channel +image_defect = image_orig.copy() +for layer in range(image_defect.shape[-1]): + image_defect[np.where(mask)] = 0 + +image_result = inpaint.inpaint_biharmonic(image_defect, mask, multichannel=True) + +fig, axes = plt.subplots(ncols=3, nrows=1) + +axes[0].set_title('Defected image') +axes[0].imshow(image_orig) +axes[0].set_xticks([]), axes[0].set_yticks([]) + +axes[1].set_title('Defect mask') +axes[1].imshow(mask, cmap=plt.cm.gray) +axes[1].set_xticks([]), axes[1].set_yticks([]) + +axes[2].set_title('Inpainted image') +axes[2].imshow(image_result) +axes[2].set_xticks([]), axes[2].set_yticks([]) + +plt.show() diff --git a/doc/examples/segmentation/plot_threshold_adaptive.py b/doc/examples/segmentation/plot_threshold_adaptive.py index 9bf7b0c6..6f473abb 100644 --- a/doc/examples/segmentation/plot_threshold_adaptive.py +++ b/doc/examples/segmentation/plot_threshold_adaptive.py @@ -26,7 +26,7 @@ image = data.page() global_thresh = threshold_otsu(image) binary_global = image > global_thresh -block_size = 40 +block_size = 35 binary_adaptive = threshold_adaptive(image, block_size, offset=10) fig, axes = plt.subplots(nrows=3, figsize=(7, 8)) diff --git a/doc/release/release_dev.txt b/doc/release/release_dev.txt index 3bc895d1..237b6b57 100644 --- a/doc/release/release_dev.txt +++ b/doc/release/release_dev.txt @@ -17,7 +17,7 @@ New Features - ``skimage.util.apply_parallel`` (#1493) - Plugin for ``imageio`` library (#1575) - +- Inpainting algorithm (#1804) Improvements ------------ diff --git a/doc/source/user_guide/numpy_images.txt b/doc/source/user_guide/numpy_images.txt index f4b3cf42..00f49bbb 100644 --- a/doc/source/user_guide/numpy_images.txt +++ b/doc/source/user_guide/numpy_images.txt @@ -77,7 +77,7 @@ disk: :: >>> nrows, ncols = camera.shape >>> row, col = np.ogrid[:nrows, :ncols] >>> cnt_row, cnt_col = nrows / 2, ncols / 2 - >>> outer_disk_mask = ((row - cnt_row)**2 + (col - cnt_col)**2 < + >>> outer_disk_mask = ((row - cnt_row)**2 + (col - cnt_col)**2 > ... (nrows / 2)**2) >>> camera[outer_disk_mask] = 0 diff --git a/doc/source/user_guide/plugins.txt b/doc/source/user_guide/plugins.txt index e17b3f59..5782a604 100644 --- a/doc/source/user_guide/plugins.txt +++ b/doc/source/user_guide/plugins.txt @@ -49,17 +49,18 @@ Any plugin in the ``_plugins`` directory is automatically examined by system:: >>> import skimage.io as io - >>> io.plugins() + >>> io.find_available_plugins() {'gtk': ['imshow'], - 'matplotlib': ['imshow', 'imsave'], - 'pil': ['imread'], - 'qt': ['imshow'], - 'test': ['imsave', 'imshow', 'imread']} + 'matplotlib': ['imshow', 'imread', 'imread_collection'], + 'pil': ['imread', 'imsave', 'imread_collection'], + 'qt': ['imshow', 'imsave', 'imread', 'imread_collection'], + 'test': ['imsave', 'imshow', 'imread', 'imread_collection'],} or only those already loaded:: - >>> io.plugins(loaded=True) - {'pil': ['imread']} + >>> io.find_available_plugins(loaded=True) + {'matplotlib': ['imshow', 'imread', 'imread_collection'], + 'pil': ['imread', 'imsave', 'imread_collection']} A plugin is loaded using the ``use_plugin`` command:: @@ -78,7 +79,7 @@ last plugin loaded is used. To query a plugin's capabilities, use ``plugin_info``:: >>> io.plugin_info('pil') - >>> + >>> {'description': 'Image reading via the Python Imaging Library', - 'provides': 'imread'} + 'provides': 'imread, imsave'} diff --git a/skimage/__init__.py b/skimage/__init__.py index 673562a7..833be389 100644 --- a/skimage/__init__.py +++ b/skimage/__init__.py @@ -158,7 +158,9 @@ else: if sys.version.startswith('2.6'): - warnings.warn("Python 2.6 is deprecated and will not be supported in scikit-image 0.13+") + msg = ("Python 2.6 is deprecated and will not be supported in " + "scikit-image 0.13+") + warnings.warn(msg, stacklevel=2) del warnings, functools, osp, imp, sys diff --git a/skimage/_shared/_warnings.py b/skimage/_shared/_warnings.py index 78d38fbd..d397d8e9 100644 --- a/skimage/_shared/_warnings.py +++ b/skimage/_shared/_warnings.py @@ -1,11 +1,20 @@ -__all__ = ['all_warnings', 'expected_warnings'] - from contextlib import contextmanager import sys import warnings import inspect import re +__all__ = ['all_warnings', 'expected_warnings', 'warn'] + + +def warn(message, category=None, stacklevel=2): + """A version of `warnings.warn` with a default stacklevel of 2. + """ + if category is not None: + warnings.warn(message, category=category, stacklevel=stacklevel) + else: + warnings.warn(message, stacklevel=stacklevel) + @contextmanager def all_warnings(): @@ -67,7 +76,7 @@ def all_warnings(): @contextmanager def expected_warnings(matching): """Context for use in testing to catch known warnings matching regexes - + Parameters ---------- matching : list of strings or compiled regexes @@ -84,15 +93,16 @@ def expected_warnings(matching): ----- Uses `all_warnings` to ensure all warnings are raised. Upon exiting, it checks the recorded warnings for the desired matching - pattern(s). + pattern(s). Raises a ValueError if any match was not found or an unexpected - warning was raised. - Allows for three types of behaviors: "and", "or", and "optional" matches. + warning was raised. + Allows for three types of behaviors: "and", "or", and "optional" matches. This is done to accomodate different build enviroments or loop conditions that may produce different warnings. The behaviors can be combined. If you pass multiple patterns, you get an orderless "and", where all of the warnings must be raised. - If you use the "|" operator in a pattern, you can catch one of several warnings. + If you use the "|" operator in a pattern, you can catch one of several + warnings. Finally, you can use "|\A\Z" in a pattern to signify it as optional. """ @@ -100,7 +110,7 @@ def expected_warnings(matching): # enter context yield w # exited user context, check the recorded warnings - remaining = [m for m in matching if not '\A\Z' in m.split('|')] + remaining = [m for m in matching if '\A\Z' not in m.split('|')] for warn in w: found = False for match in matching: diff --git a/skimage/_shared/utils.py b/skimage/_shared/utils.py index cc688fd7..96bc0fd7 100644 --- a/skimage/_shared/utils.py +++ b/skimage/_shared/utils.py @@ -6,10 +6,10 @@ import types import six -from ._warnings import all_warnings +from ._warnings import all_warnings, warn __all__ = ['deprecated', 'get_bound_method_class', 'all_warnings', - 'safe_as_int', 'assert_nD'] + 'safe_as_int', 'assert_nD', 'warn'] class skimage_deprecation(Warning): @@ -170,7 +170,7 @@ def _mode_deprecations(mode): """Used to update deprecated mode names in `skimage._shared.interpolation.pyx`.""" if mode.lower() == 'nearest': - warnings.warn(skimage_deprecation( + warn(skimage_deprecation( "Mode 'nearest' has been renamed to 'edge'. Mode 'nearest' will be " "removed in a future release.")) mode = 'edge' diff --git a/skimage/color/colorlabel.py b/skimage/color/colorlabel.py index 76cdd6f2..10d4e114 100644 --- a/skimage/color/colorlabel.py +++ b/skimage/color/colorlabel.py @@ -1,8 +1,8 @@ -import warnings import itertools import numpy as np +from .._shared.utils import warn from .. import img_as_float from . import rgb_colors from .colorconv import rgb2gray, gray2rgb @@ -148,7 +148,7 @@ def _label2rgb_overlay(label, image=None, colors=None, alpha=0.3, raise ValueError("`image` and `label` must be the same shape") if image.min() < 0: - warnings.warn("Negative intensities in `image` are not supported") + warn("Negative intensities in `image` are not supported") image = img_as_float(rgb2gray(image)) image = gray2rgb(image) * image_alpha + (1 - image_alpha) diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py index 0bf2ff91..ff8184b5 100644 --- a/skimage/data/__init__.py +++ b/skimage/data/__init__.py @@ -10,6 +10,7 @@ import os as _os from .. import data_dir from ..io import imread, use_plugin +from .._shared.utils import deprecated from ._binary_blobs import binary_blobs __all__ = ['load', @@ -56,6 +57,7 @@ def camera(): return load("camera.png") +@deprecated('skimage.data.astronaut') def lena(): """Colour "Lena" image. diff --git a/skimage/data/astronaut_GRAY_hog.npy b/skimage/data/astronaut_GRAY_hog.npy new file mode 100644 index 00000000..1e825849 Binary files /dev/null and b/skimage/data/astronaut_GRAY_hog.npy differ diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py index 62ce2ecd..b42f3169 100644 --- a/skimage/exposure/_adapthist.py +++ b/skimage/exposure/_adapthist.py @@ -19,7 +19,7 @@ import numpy as np from .. import img_as_float, img_as_uint from ..color.adapt_rgb import adapt_rgb, hsv_value from ..exposure import rescale_intensity -from .._shared.utils import skimage_deprecation, warnings +from .._shared.utils import skimage_deprecation, warn NR_OF_GREY = 2 ** 14 # number of grayscale levels to use in CLAHE algorithm @@ -77,9 +77,9 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, image = rescale_intensity(image, out_range=(0, NR_OF_GREY - 1)) if kernel_size is None: - warnings.warn('`ntiles_*` have been deprecated in favor of ' - '`kernel_size`. The `ntiles_*` keyword arguments ' - 'will be removed in v0.14', skimage_deprecation) + warn('`ntiles_*` have been deprecated in favor of ' + '`kernel_size`. The `ntiles_*` keyword arguments ' + 'will be removed in v0.14', skimage_deprecation) ntiles_x = ntiles_x or 8 ntiles_y = ntiles_y or 8 kernel_size = (np.round(image.shape[0] / ntiles_y), diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 49b59c2b..9cdb03d1 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -1,9 +1,9 @@ from __future__ import division -import warnings import numpy as np from ..color import rgb2gray from ..util.dtype import dtype_range, dtype_limits +from .._shared.utils import warn __all__ = ['histogram', 'cumulative_distribution', 'equalize_hist', @@ -60,9 +60,9 @@ def histogram(image, nbins=256): """ sh = image.shape if len(sh) == 3 and sh[-1] < 4: - warnings.warn("This might be a color image. The histogram will be " - "computed on the flattened image. You can instead " - "apply this function to each color channel.") + warn("This might be a color image. The histogram will be " + "computed on the flattened image. You can instead " + "apply this function to each color channel.") # For integer types, histogramming with bincount is more efficient. if np.issubdtype(image.dtype, np.integer): @@ -292,12 +292,12 @@ def rescale_intensity(image, in_range='image', out_range='dtype'): if in_range is None: in_range = 'image' msg = "`in_range` should not be set to None. Use {!r} instead." - warnings.warn(msg.format(in_range)) + warn(msg.format(in_range)) if out_range is None: out_range = 'dtype' msg = "`out_range` should not be set to None. Use {!r} instead." - warnings.warn(msg.format(out_range)) + warn(msg.format(out_range)) imin, imax = intensity_range(image, in_range) omin, omax = intensity_range(image, out_range, clip_negative=(imin >= 0)) diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py index 760102e7..c5d7ad90 100644 --- a/skimage/feature/_hog.py +++ b/skimage/feature/_hog.py @@ -2,11 +2,12 @@ from __future__ import division import numpy as np from .._shared.utils import assert_nD from . import _hoghistogram +import warnings def hog(image, orientations=9, pixels_per_cell=(8, 8), - cells_per_block=(3, 3), visualise=False, normalise=False, - feature_vector=True): + cells_per_block=(3, 3), visualise=False, transform_sqrt=False, + feature_vector=True, normalise=None): """Extract Histogram of Oriented Gradients (HOG) for a given image. Compute a Histogram of Oriented Gradients (HOG) by @@ -29,12 +30,16 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), Number of cells in each block. visualise : bool, optional Also return an image of the HOG. - normalise : bool, optional + transform_sqrt : bool, optional Apply power law compression to normalise the image before - processing. + processing. DO NOT use this if the image contains negative + values. Also see `notes` section below. feature_vector : bool, optional Return the data as a feature vector by calling .ravel() on the result just before returning. + normalise : bool, deprecated + The parameter is deprecated. Use `transform_sqrt` for power law + compression. `normalise` has been deprecated. Returns ------- @@ -51,6 +56,13 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), Human Detection, IEEE Computer Society Conference on Computer Vision and Pattern Recognition 2005 San Diego, CA, USA + Notes + ----- + Power law compression, also known as Gamma correction, is used to reduce + the effects of shadowing and illumination variations. The compression makes + the dark regions lighter. When the kwarg `transform_sqrt` is set to + ``True``, the function computes the square root of each color channel + and then applies the hog algorithm to the image. """ image = np.atleast_2d(image) @@ -66,7 +78,13 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), assert_nD(image, 2) - if normalise: + if normalise is not None: + raise ValueError("The normalise parameter was removed due to incorrect " + "behavior; it only applied a square root instead of a " + "true normalization. If you wish to duplicate the old " + "behavior, set ``transform_sqrt=True``.") + + if transform_sqrt: image = np.sqrt(image) """ @@ -173,7 +191,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), overlapping grid of blocks covering the detection window into a combined feature vector for use in the window classifier. """ - + if feature_vector: normalised_blocks = normalised_blocks.ravel() diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index a8fa7d78..fa680f95 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -136,7 +136,7 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0): -------- >>> from skimage.feature import hessian_matrix >>> square = np.zeros((5, 5)) - >>> square[2, 2] = 1 + >>> square[2, 2] = -1.0 / 1591.54943092 >>> Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) >>> Hxx array([[ 0., 0., 0., 0., 0.], @@ -149,22 +149,25 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0): image = _prepare_grayscale_input_2D(image) - # window extent to the left and right, which covers > 99% of the normal - # distribution + # Window extent which covers > 99% of the normal distribution. window_ext = max(1, np.ceil(3 * sigma)) ky, kx = np.mgrid[-window_ext:window_ext + 1, -window_ext:window_ext + 1] - # second derivative Gaussian kernels + # Second derivative Gaussian kernels. gaussian_exp = np.exp(-(kx ** 2 + ky ** 2) / (2 * sigma ** 2)) kernel_xx = 1 / (2 * np.pi * sigma ** 4) * (kx ** 2 / sigma ** 2 - 1) kernel_xx *= gaussian_exp - kernel_xx /= kernel_xx.sum() kernel_xy = 1 / (2 * np.pi * sigma ** 6) * (kx * ky) kernel_xy *= gaussian_exp - kernel_xy /= kernel_xx.sum() kernel_yy = kernel_xx.transpose() + # Remove small kernel values. + eps = np.finfo(kernel_xx.dtype).eps + kernel_xx[np.abs(kernel_xx) < eps * np.abs(kernel_xx).max()] = 0 + kernel_xy[np.abs(kernel_xy) < eps * np.abs(kernel_xy).max()] = 0 + kernel_yy[np.abs(kernel_yy) < eps * np.abs(kernel_yy).max()] = 0 + Hxx = ndi.convolve(image, kernel_xx, mode=mode, cval=cval) Hxy = ndi.convolve(image, kernel_xy, mode=mode, cval=cval) Hyy = ndi.convolve(image, kernel_yy, mode=mode, cval=cval) @@ -277,7 +280,7 @@ def hessian_matrix_eigvals(Hxx, Hxy, Hyy): -------- >>> from skimage.feature import hessian_matrix, hessian_matrix_eigvals >>> square = np.zeros((5, 5)) - >>> square[2, 2] = 1 + >>> square[2, 2] = -1 / 1591.54943092 >>> Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) >>> hessian_matrix_eigvals(Hxx, Hxy, Hyy)[0] array([[ 0., 0., 0., 0., 0.], @@ -796,7 +799,7 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): return corners_subpix -def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, +def corner_peaks(image, min_distance=1, threshold_abs=None, threshold_rel=0.1, exclude_border=True, indices=True, num_peaks=np.inf, footprint=None, labels=None): """Find corners in corner measure response image. @@ -820,18 +823,13 @@ def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, [ 0., 0., 1., 1., 0.], [ 0., 0., 1., 1., 0.], [ 0., 0., 0., 0., 0.]]) - >>> peak_local_max(response, exclude_border=False) + >>> peak_local_max(response) array([[2, 2], [2, 3], [3, 2], [3, 3]]) - >>> corner_peaks(response, exclude_border=False) + >>> corner_peaks(response) array([[2, 2]]) - >>> corner_peaks(response, exclude_border=False, min_distance=0) - array([[2, 2], - [2, 3], - [3, 2], - [3, 3]]) """ diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py index 421abeec..0cc1b26d 100644 --- a/skimage/feature/peak.py +++ b/skimage/feature/peak.py @@ -3,46 +3,49 @@ import scipy.ndimage as ndi from ..filters import rank_order -def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, - exclude_border=True, indices=True, num_peaks=np.inf, - footprint=None, labels=None): - """ - Find peaks in an image, and return them as coordinates or a boolean array. +def peak_local_max(image, min_distance=1, threshold_abs=None, + threshold_rel=None, exclude_border=True, indices=True, + num_peaks=np.inf, footprint=None, labels=None): + """Find peaks in an image as coordinate list or boolean mask. Peaks are the local maxima in a region of `2 * min_distance + 1` (i.e. peaks are separated by at least `min_distance`). - NOTE: If peaks are flat (i.e. multiple adjacent pixels have identical + If peaks are flat (i.e. multiple adjacent pixels have identical intensities), the coordinates of all such pixels are returned. + If both `threshold_abs` and `threshold_rel` are provided, the maximum + of the two is chosen as the minimum intensity threshold of peaks. + Parameters ---------- - image : ndarray of floats + image : ndarray Input image. - min_distance : int + min_distance : int, optional Minimum number of pixels separating peaks in a region of `2 * min_distance + 1` (i.e. peaks are separated by at least `min_distance`). If `exclude_border` is True, this value also excludes a border `min_distance` from the image boundary. To find the maximum number of peaks, use `min_distance=1`. - threshold_abs : float - Minimum intensity of peaks. - threshold_rel : float - Minimum intensity of peaks calculated as `max(image) * threshold_rel`. - exclude_border : bool + threshold_abs : float, optional + Minimum intensity of peaks. By default, the absolute threshold is + the minimum intensity of the image. + threshold_rel : float, optional + Minimum intensity of peaks, calculated as `max(image) * threshold_rel`. + exclude_border : bool, optional If True, `min_distance` excludes peaks from the border of the image as well as from each other. - indices : bool - If True, the output will be an array representing peak coordinates. - If False, the output will be a boolean array shaped as `image.shape` - with peaks present at True elements. - num_peaks : int + indices : bool, optional + If True, the output will be an array representing peak + coordinates. If False, the output will be a boolean array shaped as + `image.shape` with peaks present at True elements. + num_peaks : int, optional Maximum number of peaks. When the number of peaks exceeds `num_peaks`, return `num_peaks` peaks based on highest peak intensity. footprint : ndarray of bools, optional If provided, `footprint == 1` represents the local region within which to search for peaks at every point in `image`. Overrides - `min_distance`, except for border exclusion if `exclude_border=True`. + `min_distance` (also for `exclude_border`). labels : ndarray of ints, optional If provided, each unique region `labels == value` represents a unique region to search for peaks. Zero is reserved for background. @@ -58,10 +61,10 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, Notes ----- The peak local maximum function returns the coordinates of local peaks - (maxima) in a image. A maximum filter is used for finding local maxima. - This operation dilates the original image. After comparison between - dilated and original image, peak_local_max function returns the - coordinates of peaks where dilated image = original. + (maxima) in an image. A maximum filter is used for finding local maxima. + This operation dilates the original image. After comparison of the dilated + and original image, this function returns the coordinates or a mask of the + peaks where the dilated image equals the original image. Examples -------- @@ -90,7 +93,9 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, array([[10, 10, 10]]) """ + out = np.zeros_like(image, dtype=np.bool) + # In the case of labels, recursively build and return an output # operating on each label separately if labels is not None: @@ -123,7 +128,6 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, else: return out - image = image.copy() # Non maximum filter if footprint is not None: image_max = ndi.maximum_filter(image, footprint=footprint, @@ -131,25 +135,33 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, else: size = 2 * min_distance + 1 image_max = ndi.maximum_filter(image, size=size, mode='constant') - mask = (image == image_max) - image *= mask + mask = image == image_max - if exclude_border: + if exclude_border and (footprint is not None or min_distance > 0): # zero out the image borders - for i in range(image.ndim): - image = image.swapaxes(0, i) - image[:min_distance] = 0 - image[-min_distance:] = 0 - image = image.swapaxes(0, i) + for i in range(mask.ndim): + mask = mask.swapaxes(0, i) + remove = (footprint.shape[i] if footprint is not None + else 2 * min_distance) + mask[:remove // 2] = mask[-remove // 2:] = False + mask = mask.swapaxes(0, i) # find top peak candidates above a threshold - peak_threshold = max(np.max(image.ravel()) * threshold_rel, threshold_abs) + thresholds = [] + if threshold_abs is None: + threshold_abs = image.min() + thresholds.append(threshold_abs) + if threshold_rel is not None: + thresholds.append(threshold_rel * image.max()) + if thresholds: + mask &= image > max(thresholds) # get coordinates of peaks - coordinates = np.argwhere(image > peak_threshold) + coordinates = np.transpose(mask.nonzero()) if coordinates.shape[0] > num_peaks: - intensities = image.flat[np.ravel_multi_index(coordinates.transpose(),image.shape)] + intensities = image.flat[np.ravel_multi_index(coordinates.transpose(), + image.shape)] idx_maxsort = np.argsort(intensities)[::-1] coordinates = coordinates[idx_maxsort][:num_peaks] diff --git a/skimage/feature/tests/test_brief.py b/skimage/feature/tests/test_brief.py index e70a1148..7f6ae9b7 100644 --- a/skimage/feature/tests/test_brief.py +++ b/skimage/feature/tests/test_brief.py @@ -16,44 +16,46 @@ def test_color_image_unsupported_error(): def test_normal_mode(): """Verify the computed BRIEF descriptors with expected for normal mode.""" - img = rgb2gray(data.lena()) + img = data.coins() - keypoints = corner_peaks(corner_harris(img), min_distance=5) + keypoints = corner_peaks(corner_harris(img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor = BRIEF(descriptor_size=8, sigma=2) extractor.extract(img, keypoints[:8]) - expected = np.array([[ True, False, True, False, True, True, False, False], + expected = np.array([[False, True, False, False, True, False, True, False], + [ True, False, True, True, False, True, False, False], + [ True, False, False, True, False, True, False, True], + [ True, True, True, True, False, True, False, True], + [ True, True, True, False, False, True, True, True], [False, False, False, False, True, False, False, False], - [ True, True, True, True, True, True, True, True], - [ True, False, True, True, False, True, False, True], - [False, True, True, True, True, True, True, True], - [ True, False, False, False, False, True, False, True], - [False, True, True, True, False, False, True, False], - [False, False, False, False, True, False, False, False]], dtype=bool) + [False, True, False, False, True, False, True, False], + [False, False, False, False, False, False, False, False]], dtype=bool) assert_array_equal(extractor.descriptors, expected) def test_uniform_mode(): """Verify the computed BRIEF descriptors with expected for uniform mode.""" - img = rgb2gray(data.lena()) + img = data.coins() - keypoints = corner_peaks(corner_harris(img), min_distance=5) + keypoints = corner_peaks(corner_harris(img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor = BRIEF(descriptor_size=8, sigma=2, mode='uniform') extractor.extract(img, keypoints[:8]) - expected = np.array([[ True, False, True, False, False, True, False, False], - [False, True, False, False, True, True, True, True], - [ True, False, False, False, False, False, False, False], - [False, True, True, False, False, False, True, False], - [False, False, False, False, False, False, True, False], - [False, True, False, False, True, False, False, False], - [False, False, True, True, False, False, True, True], - [ True, True, False, False, False, False, False, False]], dtype=bool) + expected = np.array([[False, False, False, True, True, True, False, False], + [ True, True, True, False, True, False, False, True], + [ True, True, True, False, True, True, False, True], + [ True, True, True, True, False, True, False, True], + [ True, True, True, True, True, True, False, False], + [ True, True, True, True, True, True, True, True], + [False, False, False, True, True, True, True, True], + [False, True, False, True, False, True, True, True]], dtype=bool) assert_array_equal(extractor.descriptors, expected) diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index d6c3418f..c3919fe8 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -42,21 +42,21 @@ def test_hessian_matrix(): square = np.zeros((5, 5)) square[2, 2] = 1 Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) - assert_array_equal(Hxx, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(Hxy, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(Hyy, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) + assert_almost_equal(Hxx, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(Hxy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(Hyy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) def test_structure_tensor_eigvals(): @@ -81,16 +81,16 @@ def test_hessian_matrix_eigvals(): square[2, 2] = 1 Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) l1, l2 = hessian_matrix_eigvals(Hxx, Hxy, Hyy) - assert_array_equal(l1, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(l2, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) + assert_almost_equal(l1, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(l2, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) @test_parallel() @@ -107,21 +107,25 @@ def test_square_image(): im[:25, :25] = 1. # Moravec - results = peak_local_max(corner_moravec(im)) + results = peak_local_max(corner_moravec(im), + min_distance=10, threshold_rel=0) # interest points along edge assert len(results) == 57 # Harris - results = peak_local_max(corner_harris(im, method='k')) + results = peak_local_max(corner_harris(im, method='k'), + min_distance=10, threshold_rel=0) # interest at corner assert len(results) == 1 - results = peak_local_max(corner_harris(im, method='eps')) + results = peak_local_max(corner_harris(im, method='eps'), + min_distance=10, threshold_rel=0) # interest at corner assert len(results) == 1 # Shi-Tomasi - results = peak_local_max(corner_shi_tomasi(im)) + results = peak_local_max(corner_shi_tomasi(im), + min_distance=10, threshold_rel=0) # interest at corner assert len(results) == 1 @@ -133,18 +137,22 @@ def test_noisy_square_image(): im = im + np.random.uniform(size=im.shape) * .2 # Moravec - results = peak_local_max(corner_moravec(im)) + results = peak_local_max(corner_moravec(im), + min_distance=10, threshold_rel=0) # undefined number of interest points assert results.any() # Harris - results = peak_local_max(corner_harris(im, sigma=1.5, method='k')) + results = peak_local_max(corner_harris(im, method='k'), + min_distance=10, threshold_rel=0) assert len(results) == 1 - results = peak_local_max(corner_harris(im, sigma=1.5, method='eps')) + results = peak_local_max(corner_harris(im, method='eps'), + min_distance=10, threshold_rel=0) assert len(results) == 1 # Shi-Tomasi - results = peak_local_max(corner_shi_tomasi(im, sigma=1.5)) + results = peak_local_max(corner_shi_tomasi(im, sigma=1.5), + min_distance=10, threshold_rel=0) assert len(results) == 1 @@ -156,11 +164,13 @@ def test_squared_dot(): # Moravec fails # Harris - results = peak_local_max(corner_harris(im)) + results = peak_local_max(corner_harris(im), + min_distance=10, threshold_rel=0) assert (results == np.array([[6, 6]])).all() # Shi-Tomasi - results = peak_local_max(corner_shi_tomasi(im)) + results = peak_local_max(corner_shi_tomasi(im), + min_distance=10, threshold_rel=0) assert (results == np.array([[6, 6]])).all() @@ -173,20 +183,26 @@ def test_rotated_img(): im_rotated = im.T # Moravec - results = peak_local_max(corner_moravec(im)) - results_rotated = peak_local_max(corner_moravec(im_rotated)) + results = peak_local_max(corner_moravec(im), + min_distance=10, threshold_rel=0) + results_rotated = peak_local_max(corner_moravec(im_rotated), + min_distance=10, threshold_rel=0) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() # Harris - results = peak_local_max(corner_harris(im)) - results_rotated = peak_local_max(corner_harris(im_rotated)) + results = peak_local_max(corner_harris(im), + min_distance=10, threshold_rel=0) + results_rotated = peak_local_max(corner_harris(im_rotated), + min_distance=10, threshold_rel=0) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() # Shi-Tomasi - results = peak_local_max(corner_shi_tomasi(im)) - results_rotated = peak_local_max(corner_shi_tomasi(im_rotated)) + results = peak_local_max(corner_shi_tomasi(im), + min_distance=10, threshold_rel=0) + results_rotated = peak_local_max(corner_shi_tomasi(im_rotated), + min_distance=10, threshold_rel=0) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() @@ -195,7 +211,8 @@ def test_subpix_edge(): img = np.zeros((50, 50)) img[:25, :25] = 255 img[25:, 25:] = 255 - corner = peak_local_max(corner_harris(img), num_peaks=1) + corner = peak_local_max(corner_harris(img), + min_distance=10, threshold_rel=0, num_peaks=1) subpix = corner_subpix(img, corner) assert_array_equal(subpix[0], (24.5, 24.5)) @@ -203,7 +220,8 @@ def test_subpix_edge(): def test_subpix_dot(): img = np.zeros((50, 50)) img[25, 25] = 255 - corner = peak_local_max(corner_harris(img), num_peaks=1) + corner = peak_local_max(corner_harris(img), + min_distance=10, threshold_rel=0, num_peaks=1) subpix = corner_subpix(img, corner) assert_array_equal(subpix[0], (25, 25)) @@ -214,7 +232,8 @@ def test_subpix_no_class(): assert_array_equal(subpix[0], (np.nan, np.nan)) img[25, 25] = 1e-10 - corner = peak_local_max(corner_harris(img), num_peaks=1) + corner = peak_local_max(corner_harris(img), + min_distance=10, threshold_rel=0, num_peaks=1) subpix = corner_subpix(img, np.array([[25, 25]])) assert_array_equal(subpix[0], (np.nan, np.nan)) @@ -223,7 +242,7 @@ def test_subpix_border(): img = np.zeros((50, 50)) img[1:25,1:25] = 255 img[25:-1,25:-1] = 255 - corner = corner_peaks(corner_harris(img), min_distance=1) + corner = corner_peaks(corner_harris(img), threshold_rel=0) subpix = corner_subpix(img, corner, window_size=11) ref = np.array([[ 0.52040816, 0.52040816], [ 0.52040816, 24.47959184], @@ -244,21 +263,23 @@ def test_num_peaks(): for i in range(20): n = np.random.random_integers(20) - results = peak_local_max(img_corners, num_peaks=n) + results = peak_local_max(img_corners, + min_distance=10, threshold_rel=0, num_peaks=n) assert (results.shape[0] == n) def test_corner_peaks(): - response = np.zeros((5, 5)) - response[2:4, 2:4] = 1 + response = np.zeros((10, 10)) + response[2:5, 2:5] = 1 - corners = corner_peaks(response, exclude_border=False) + corners = corner_peaks(response, exclude_border=False, min_distance=10, + threshold_rel=0) assert len(corners) == 1 - corners = corner_peaks(response, exclude_border=False, min_distance=0) + corners = corner_peaks(response, exclude_border=False, min_distance=1) assert len(corners) == 4 - corners = corner_peaks(response, exclude_border=False, min_distance=0, + corners = corner_peaks(response, exclude_border=False, min_distance=1, indices=False) assert np.sum(corners) == 4 @@ -323,7 +344,8 @@ def test_corner_fast_lena(): [492, 139], [494, 169], [496, 266]]) - actual = corner_peaks(corner_fast(img, 12, 0.3)) + actual = corner_peaks(corner_fast(img, 12, 0.3), + min_distance=10, threshold_rel=0) assert_array_equal(actual, expected) @@ -340,11 +362,22 @@ def test_corner_orientations_even_shape_error(): @test_parallel() -def test_corner_orientations_lena(): - img = rgb2gray(data.lena()) - corners = corner_peaks(corner_fast(img, 11, 0.35)) - expected = np.array([-1.9195897 , -3.03159624, -1.05991162, -2.89573739, - -2.61607644, 2.98660159]) +def test_corner_orientations_astronaut(): + img = rgb2gray(data.astronaut()) + corners = corner_peaks(corner_fast(img, 11, 0.35), + min_distance=10, threshold_abs=0, threshold_rel=0.1) + expected = np.array([-1.75220190e+00, 2.01197383e+00, -2.01162417e+00, + -1.88247204e-01, 1.19134149e+00, -6.61151410e-01, + -2.99143370e+00, 2.17103132e+00, -7.52950306e-04, + 1.25854853e+00, 2.43573659e+00, -1.69230287e+00, + -9.88548213e-01, 1.47154532e+00, -1.65449964e+00, + 1.09650167e+00, 1.07812134e+00, -1.68885773e+00, + -1.64397304e+00, 3.09780364e+00, -3.49561988e-01, + -1.46554357e+00, -2.81524886e+00, 8.12701702e-01, + 2.47305654e+00, -1.63869275e+00, 5.46905279e-02, + -4.40598471e-01, 3.14918803e-01, -1.76069982e+00, + 3.05330950e+00, 2.39291733e+00, -1.22091334e-01, + -3.09279990e-01, 1.45931342e+00]) actual = corner_orientations(img, corners, octagon(3, 2)) assert_almost_equal(actual, expected) @@ -352,7 +385,8 @@ def test_corner_orientations_lena(): def test_corner_orientations_square(): square = np.zeros((12, 12)) square[3:9, 3:9] = 1 - corners = corner_peaks(corner_fast(square, 9), min_distance=1) + corners = corner_peaks(corner_fast(square, 9), + min_distance=1, threshold_rel=0) actual_orientations = corner_orientations(square, corners, octagon(3, 2)) actual_orientations_degrees = np.rad2deg(actual_orientations) expected_orientations_degree = np.array([ 45., 135., -45., -135.]) diff --git a/skimage/feature/tests/test_hog.py b/skimage/feature/tests/test_hog.py index 3afc0d89..0a31b910 100644 --- a/skimage/feature/tests/test_hog.py +++ b/skimage/feature/tests/test_hog.py @@ -2,6 +2,7 @@ import os import numpy as np from scipy import ndimage as ndi import skimage as si +from skimage import color from skimage import data from skimage import feature from skimage import img_as_float @@ -21,13 +22,13 @@ def test_histogram_of_oriented_gradients_output_size(): def test_histogram_of_oriented_gradients_output_correctness(): - img = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8.npy')) - correct_output = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8_hog.npy')) - - output = feature.hog(img, orientations=9, pixels_per_cell=(8, 8), + img = color.rgb2gray(data.astronaut()) + correct_output = np.load(os.path.join(si.data_dir, 'astronaut_GRAY_hog.npy')) + + output = feature.hog(img, orientations=9, pixels_per_cell=(8, 8), cells_per_block=(3, 3), feature_vector=True, - normalise=False, visualise=False) - + transform_sqrt=False, visualise=False) + assert_almost_equal(output, correct_output) @@ -48,7 +49,7 @@ def test_hog_basic_orientations_and_data_types(): # 1) create image (with float values) where upper half is filled by # zeros, bottom half by 100 # 2) create unsigned integer version of this image - # 3) calculate feature.hog() for both images, both with 'normalise' + # 3) calculate feature.hog() for both images, both with 'transform_sqrt' # option enabled and disabled # 4) verify that all results are equal where expected # 5) verify that computed feature vector is as expected @@ -69,16 +70,16 @@ def test_hog_basic_orientations_and_data_types(): (hog_float, hog_img_float) = feature.hog( image_float, orientations=4, pixels_per_cell=(8, 8), - cells_per_block=(1, 1), visualise=True, normalise=False) + cells_per_block=(1, 1), visualise=True, transform_sqrt=False) (hog_uint8, hog_img_uint8) = feature.hog( image_uint8, orientations=4, pixels_per_cell=(8, 8), - cells_per_block=(1, 1), visualise=True, normalise=False) + cells_per_block=(1, 1), visualise=True, transform_sqrt=False) (hog_float_norm, hog_img_float_norm) = feature.hog( image_float, orientations=4, pixels_per_cell=(8, 8), - cells_per_block=(1, 1), visualise=True, normalise=True) + cells_per_block=(1, 1), visualise=True, transform_sqrt=True) (hog_uint8_norm, hog_img_uint8_norm) = feature.hog( image_uint8, orientations=4, pixels_per_cell=(8, 8), - cells_per_block=(1, 1), visualise=True, normalise=True) + cells_per_block=(1, 1), visualise=True, transform_sqrt=True) # set to True to enable manual debugging with graphical output, # must be False for automatic testing @@ -100,11 +101,11 @@ def test_hog_basic_orientations_and_data_types(): plt.subplot(2, 3, 3) plt.imshow(hog_img_float_norm) plt.colorbar() - plt.title('HOG result (normalise) visualisation (float img)') + plt.title('HOG result (transform_sqrt) visualisation (float img)') plt.subplot(2, 3, 6) plt.imshow(hog_img_uint8_norm) plt.colorbar() - plt.title('HOG result (normalise) visualisation (uint8 img)') + plt.title('HOG result (transform_sqrt) visualisation (uint8 img)') plt.show() # results (features and visualisation) for float and uint8 images must @@ -112,7 +113,7 @@ def test_hog_basic_orientations_and_data_types(): assert_almost_equal(hog_float, hog_uint8) assert_almost_equal(hog_img_float, hog_img_uint8) - # resulting features should be almost equal when 'normalise' is enabled + # resulting features should be almost equal when 'transform_sqrt' is enabled # or disabled (for current simple testing image) assert_almost_equal(hog_float, hog_float_norm, decimal=4) assert_almost_equal(hog_float, hog_uint8_norm, decimal=4) @@ -156,7 +157,7 @@ def test_hog_orientations_circle(): (hog, hog_img) = feature.hog(image, orientations=orientations, pixels_per_cell=(8, 8), cells_per_block=(1, 1), visualise=True, - normalise=False) + transform_sqrt=False) # set to True to enable manual debugging with graphical output, # must be False for automatic testing @@ -187,5 +188,9 @@ def test_hog_orientations_circle(): assert_almost_equal(actual, desired, decimal=1) +def test_hog_normalise_none_error_raised(): + img = np.array([1, 2, 3]) + assert_raises(ValueError, feature.hog, img, normalise=True) + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/feature/tests/test_match.py b/skimage/feature/tests/test_match.py index 695116ed..770dc37b 100644 --- a/skimage/feature/tests/test_match.py +++ b/skimage/feature/tests/test_match.py @@ -29,18 +29,20 @@ def test_binary_descriptors_lena_rotation_crosscheck_false(): """Verify matched keypoints and their corresponding masks results between lena image and its rotated version with the expected keypoint pairs with cross_check disabled.""" - img = data.lena() + img = data.astronaut() img = rgb2gray(img) tform = tf.SimilarityTransform(scale=1, rotation=0.15, translation=(0, 0)) rotated_img = tf.warp(img, tform, clip=False) extractor = BRIEF(descriptor_size=512) - keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + keypoints1 = corner_peaks(corner_harris(img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor.extract(img, keypoints1) descriptors1 = extractor.descriptors - keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5) + keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor.extract(rotated_img, keypoints2) descriptors2 = extractor.descriptors @@ -50,10 +52,10 @@ def test_binary_descriptors_lena_rotation_crosscheck_false(): 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]) - exp_matches2 = np.array([33, 0, 35, 7, 1, 35, 3, 2, 3, 6, 4, 9, - 11, 10, 28, 7, 8, 5, 31, 14, 13, 15, 21, 16, - 16, 13, 17, 18, 19, 21, 22, 23, 0, 24, 1, 24, - 23, 0, 26, 27, 25, 34, 28, 14, 29, 30, 21]) + exp_matches2 = np.array([ 0, 31, 2, 3, 1, 4, 6, 4, 38, 5, 27, 7, + 13, 10, 9, 27, 7, 11, 15, 8, 23, 14, 12, 16, + 10, 25, 18, 19, 21, 20, 41, 24, 25, 26, 28, 27, + 22, 23, 29, 30, 31, 32, 35, 33, 34, 30, 36]) assert_equal(matches[:, 0], exp_matches1) assert_equal(matches[:, 1], exp_matches2) @@ -62,29 +64,31 @@ def test_binary_descriptors_lena_rotation_crosscheck_true(): """Verify matched keypoints and their corresponding masks results between lena image and its rotated version with the expected keypoint pairs with cross_check enabled.""" - img = data.lena() + img = data.astronaut() img = rgb2gray(img) tform = tf.SimilarityTransform(scale=1, rotation=0.15, translation=(0, 0)) rotated_img = tf.warp(img, tform, clip=False) extractor = BRIEF(descriptor_size=512) - keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + keypoints1 = corner_peaks(corner_harris(img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor.extract(img, keypoints1) descriptors1 = extractor.descriptors - keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5) + keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5, + threshold_abs=0, threshold_rel=0.1) extractor.extract(rotated_img, keypoints2) descriptors2 = extractor.descriptors matches = match_descriptors(descriptors1, descriptors2, cross_check=True) - exp_matches1 = np.array([ 0, 1, 2, 4, 6, 7, 9, 10, 11, 12, 13, 15, - 16, 17, 19, 20, 21, 24, 26, 27, 28, 29, 30, 35, - 36, 38, 39, 40, 42, 44, 45]) - exp_matches2 = np.array([33, 0, 35, 1, 3, 2, 6, 4, 9, 11, 10, 7, - 8, 5, 14, 13, 15, 16, 17, 18, 19, 21, 22, 24, - 23, 26, 27, 25, 28, 29, 30]) + exp_matches1 = np.array([ 0, 2, 3, 4, 5, 6, 9, 11, 12, 13, 14, 17, + 18, 19, 21, 22, 23, 26, 27, 28, 29, 31, 32, 33, + 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 46]) + exp_matches2 = np.array([ 0, 2, 3, 1, 4, 6, 5, 7, 13, 10, 9, 11, + 15, 8, 14, 12, 16, 18, 19, 21, 20, 24, 25, 26, + 28, 27, 22, 23, 29, 30, 31, 32, 35, 33, 34, 36]) assert_equal(matches[:, 0], exp_matches1) assert_equal(matches[:, 1], exp_matches2) diff --git a/skimage/feature/tests/test_orb.py b/skimage/feature/tests/test_orb.py index 8895d857..d95d2d4a 100644 --- a/skimage/feature/tests/test_orb.py +++ b/skimage/feature/tests/test_orb.py @@ -2,11 +2,10 @@ import numpy as np from numpy.testing import assert_equal, assert_almost_equal, run_module_suite from skimage.feature import ORB from skimage import data -from skimage.color import rgb2gray from skimage._shared.testing import test_parallel -img = rgb2gray(data.lena()) +img = data.coins() @test_parallel() @@ -14,22 +13,21 @@ def test_keypoints_orb_desired_no_of_keypoints(): detector_extractor = ORB(n_keypoints=10, fast_n=12, fast_threshold=0.20) detector_extractor.detect(img) - exp_rows = np.array([ 435. , 435.6 , 376. , 455. , 434.88, 269. , - 375.6 , 310.8 , 413. , 311.04]) - exp_cols = np.array([ 180. , 180. , 156. , 176. , 180. , 111. , - 156. , 172.8, 70. , 172.8]) + exp_rows = np.array([ 141. , 108. , 214.56 , 131. , 214.272, + 67. , 206. , 177. , 108. , 141. ]) + exp_cols = np.array([ 323. , 328. , 282.24 , 292. , 281.664, + 85. , 260. , 284. , 328.8 , 267. ]) - exp_scales = np.array([ 1. , 1.2 , 1. , 1. , 1.44 , 1. , - 1.2 , 1.2 , 1. , 1.728]) + exp_scales = np.array([ 323. , 328. , 282.24 , 292. , 281.664, + 85. , 260. , 284. , 328.8 , 267. ]) - exp_orientations = np.array([-175.64733392, -167.94842949, -148.98350192, - -142.03599837, -176.08535837, -53.08162354, - -150.89208271, 97.7693776 , -173.4479964 , - 38.66312042]) - exp_response = np.array([ 0.96770745, 0.81027306, 0.72376257, - 0.5626413 , 0.5097993 , 0.44351774, - 0.39154173, 0.39084861, 0.39063076, - 0.37602487]) + exp_orientations = np.array([ -53.97446153, 59.5055285 , -96.01885186, + -149.70789506, -94.70171899, -45.76429535, + -51.49752849, 113.57081195, 63.30428063, + -79.56091118]) + exp_response = np.array([ 1.01168357, 0.82934145, 0.67784179, 0.57176438, + 0.56637459, 0.52248355, 0.43696175, 0.42992376, + 0.37700486, 0.36126832]) assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) @@ -48,20 +46,16 @@ def test_keypoints_orb_less_than_desired_no_of_keypoints(): fast_threshold=0.33, downscale=2, n_scales=2) detector_extractor.detect(img) - exp_rows = np.array([ 67., 247., 269., 413., 435., 230., 264., - 330., 372.]) - exp_cols = np.array([ 157., 146., 111., 70., 180., 136., 336., - 148., 156.]) + exp_rows = np.array([ 58., 65., 108., 140., 203.]) + exp_cols = np.array([ 291., 130., 293., 202., 267.]) - exp_scales = np.array([ 1., 1., 1., 1., 1., 2., 2., 2., 2.]) + exp_scales = np.array([1., 1., 1., 1., 1.]) - exp_orientations = np.array([-105.76503839, -96.28973044, -53.08162354, - -173.4479964 , -175.64733392, -106.07927215, - -163.40016243, 75.80865813, -154.73195911]) + exp_orientations = np.array([-158.26941428, -59.42996346, 151.93905955, + -79.46341354, -56.90052451]) - exp_response = np.array([ 0.13197835, 0.24931321, 0.44351774, - 0.39063076, 0.96770745, 0.04935129, - 0.21431068, 0.15826555, 0.42403573]) + exp_response = np.array([ 0.2667641 , 0.04009017, -0.17641695, -0.03243431, + 0.26521259]) assert_almost_equal(exp_rows, detector_extractor.keypoints[:, 0]) assert_almost_equal(exp_cols, detector_extractor.keypoints[:, 1]) @@ -78,27 +72,26 @@ def test_keypoints_orb_less_than_desired_no_of_keypoints(): def test_descriptor_orb(): detector_extractor = ORB(fast_n=12, fast_threshold=0.20) - exp_descriptors = np.array([[ True, False, True, True, False, False, False, False, False, False], - [False, False, True, True, False, True, True, False, True, True], - [ True, False, False, False, True, False, True, True, True, False], - [ True, False, False, True, False, True, True, False, False, False], - [False, True, True, True, False, False, False, True, True, False], - [False, False, False, False, False, True, False, True, True, True], - [False, True, True, True, True, False, False, True, False, True], - [ True, True, True, False, True, True, True, True, False, False], - [ True, True, False, True, True, True, True, False, False, False], - [ True, False, False, False, False, True, False, False, True, True], - [ True, False, False, False, True, True, True, False, False, False], - [False, False, True, False, True, False, False, True, False, False], - [False, False, True, True, False, False, False, False, False, True], - [ True, True, False, False, False, True, True, True, True, True], - [ True, True, True, False, False, True, False, True, True, False], - [False, True, True, False, False, True, True, True, True, True], - [ True, True, True, False, False, False, False, True, True, True], - [False, False, False, False, True, False, False, True, True, False], - [False, True, False, False, True, False, False, False, True, True], - [ True, False, True, False, False, False, True, True, False, False]], dtype=bool) - + exp_descriptors = np.array([[0, 1, 1, 1, 0, 1, 0, 1, 0, 1], + [1, 1, 1, 0, 0, 1, 0, 0, 1, 1], + [1, 0, 1, 1, 0, 0, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 1, 0], + [1, 1, 0, 1, 1, 1, 0, 0, 1, 1], + [1, 1, 0, 1, 0, 0, 1, 0, 1, 1], + [0, 0, 1, 0, 1, 0, 0, 1, 1, 0], + [1, 0, 0, 0, 1, 0, 0, 0, 0, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 0, 1, 0, 1, 0, 0, 1, 1], + [1, 1, 1, 0, 0, 0, 1, 1, 1, 0], + [1, 1, 1, 1, 1, 1, 0, 0, 0, 0], + [1, 1, 1, 0, 1, 1, 1, 1, 0, 0], + [1, 1, 0, 0, 1, 0, 0, 1, 0, 1], + [1, 1, 0, 0, 0, 0, 1, 0, 0, 1], + [0, 0, 0, 0, 1, 1, 1, 0, 1, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 1], + [0, 0, 0, 0, 0, 1, 1, 0, 1, 1], + [0, 0, 0, 0, 1, 0, 1, 0, 1, 1]], dtype=bool) detector_extractor.detect(img) detector_extractor.extract(img, detector_extractor.keypoints, detector_extractor.scales, diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py index e6256208..40292cd9 100644 --- a/skimage/feature/tests/test_peak.py +++ b/skimage/feature/tests/test_peak.py @@ -1,6 +1,6 @@ import numpy as np from numpy.testing import (assert_array_almost_equal as assert_close, - assert_equal) + assert_equal, assert_raises) from scipy import ndimage as ndi from skimage.feature import peak @@ -70,12 +70,14 @@ def test_num_peaks(): image[1, 5] = 12 image[3, 5] = 8 image[5, 3] = 7 - assert len(peak.peak_local_max(image, min_distance=1)) == 5 - peaks_limited = peak.peak_local_max(image, min_distance=1, num_peaks=2) + assert len(peak.peak_local_max(image, min_distance=1, threshold_abs=0)) == 5 + peaks_limited = peak.peak_local_max( + image, min_distance=1, threshold_abs=0, num_peaks=2) assert len(peaks_limited) == 2 assert (1, 3) in peaks_limited assert (1, 5) in peaks_limited - peaks_limited = peak.peak_local_max(image, min_distance=1, num_peaks=4) + peaks_limited = peak.peak_local_max( + image, min_distance=1, threshold_abs=0, num_peaks=4) assert len(peaks_limited) == 4 assert (1, 3) in peaks_limited assert (1, 5) in peaks_limited @@ -270,9 +272,11 @@ def test_disk(): result = peak.peak_local_max(image, labels=np.ones((10, 20)), footprint=footprint, min_distance=1, threshold_rel=0, - indices=False, exclude_border=False) + threshold_abs=-1, indices=False, + exclude_border=False) assert np.all(result) - result = peak.peak_local_max(image, footprint=footprint) + result = peak.peak_local_max(image, footprint=footprint, threshold_abs=-1, + indices=False, exclude_border=False) assert np.all(result) @@ -280,11 +284,14 @@ def test_3D(): image = np.zeros((30, 30, 30)) image[15, 15, 15] = 1 image[5, 5, 5] = 1 - assert_equal(peak.peak_local_max(image), [[15, 15, 15]]) - assert_equal(peak.peak_local_max(image, min_distance=6), [[15, 15, 15]]) - assert_equal(peak.peak_local_max(image, exclude_border=False), + assert_equal(peak.peak_local_max(image, min_distance=10, threshold_rel=0), + [[15, 15, 15]]) + assert_equal(peak.peak_local_max(image, min_distance=6, threshold_rel=0), + [[15, 15, 15]]) + assert_equal(peak.peak_local_max(image, min_distance=10, threshold_rel=0, + exclude_border=False), [[5, 5, 5], [15, 15, 15]]) - assert_equal(peak.peak_local_max(image, min_distance=5), + assert_equal(peak.peak_local_max(image, min_distance=5, threshold_rel=0), [[5, 5, 5], [15, 15, 15]]) @@ -292,14 +299,30 @@ def test_4D(): image = np.zeros((30, 30, 30, 30)) image[15, 15, 15, 15] = 1 image[5, 5, 5, 5] = 1 - assert_equal(peak.peak_local_max(image), [[15, 15, 15, 15]]) - assert_equal(peak.peak_local_max(image, min_distance=6), [[15, 15, 15, 15]]) - assert_equal(peak.peak_local_max(image, exclude_border=False), + assert_equal(peak.peak_local_max(image, min_distance=10, threshold_rel=0), + [[15, 15, 15, 15]]) + assert_equal(peak.peak_local_max(image, min_distance=6, threshold_rel=0), + [[15, 15, 15, 15]]) + assert_equal(peak.peak_local_max(image, min_distance=10, threshold_rel=0, + exclude_border=False), [[5, 5, 5, 5], [15, 15, 15, 15]]) - assert_equal(peak.peak_local_max(image, min_distance=5), + assert_equal(peak.peak_local_max(image, min_distance=5, threshold_rel=0), [[5, 5, 5, 5], [15, 15, 15, 15]]) +def test_threshold_rel_default(): + image = np.ones((5, 5)) + + image[2, 2] = 1 + assert len(peak.peak_local_max(image)) == 0 + + image[2, 2] = 2 + assert_equal(peak.peak_local_max(image), [[2, 2]]) + + image[2, 2] = 0 + assert len(peak.peak_local_max(image, min_distance=0)) == image.size - 1 + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() diff --git a/skimage/filters/_gaussian.py b/skimage/filters/_gaussian.py index ad4e621d..2a455fe9 100644 --- a/skimage/filters/_gaussian.py +++ b/skimage/filters/_gaussian.py @@ -1,10 +1,10 @@ import collections as coll import numpy as np from scipy import ndimage as ndi -import warnings from ..util import img_as_float from ..color import guess_spatial_dimensions +from .._shared.utils import warn __all__ = ['gaussian'] @@ -91,7 +91,7 @@ def gaussian(image, sigma, output=None, mode='nearest', cval=0, msg = ("Images with dimensions (M, N, 3) are interpreted as 2D+RGB " "by default. Use `multichannel=False` to interpret as " "3D image with last dimension of length 3.") - warnings.warn(RuntimeWarning(msg)) + warn(RuntimeWarning(msg)) multichannel = True if np.any(np.asarray(sigma) < 0.0): raise ValueError("Sigma values less than zero are not valid") diff --git a/skimage/filters/rank/generic.py b/skimage/filters/rank/generic.py index 5c6a352b..ec369153 100644 --- a/skimage/filters/rank/generic.py +++ b/skimage/filters/rank/generic.py @@ -16,17 +16,16 @@ References """ -import warnings import numpy as np from ... import img_as_ubyte -from ..._shared.utils import assert_nD +from ..._shared.utils import assert_nD, warn from . import generic_cy __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', - 'geometric_mean', 'subtract_mean', 'median', 'minimum', 'modal', - 'enhance_contrast', 'pop', 'threshold', 'tophat', 'noise_filter', + 'geometric_mean', 'subtract_mean', 'median', 'minimum', 'modal', + 'enhance_contrast', 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] @@ -65,8 +64,8 @@ def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): bitdepth = int(np.log2(max_bin)) if bitdepth > 10: - warnings.warn("Bitdepth of %d may result in bad rank filter " - "performance due to large number of bins." % bitdepth) + warn("Bitdepth of %d may result in bad rank filter " + "performance due to large number of bins." % bitdepth) return image, selem, out, mask, max_bin @@ -377,7 +376,7 @@ def geometric_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=Fal References ---------- - .. [1] Gonzalez, R. C. and Wood, R. E. "Digital Image Processing (3rd Edition)." + .. [1] Gonzalez, R. C. and Wood, R. E. "Digital Image Processing (3rd Edition)." Prentice-Hall Inc, 2006. """ diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index a47b6f20..03557caf 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/skimage/filters/tests/test_thresholding.py @@ -1,8 +1,11 @@ import numpy as np -from numpy.testing import assert_equal, assert_almost_equal +from numpy.testing import (assert_equal, + assert_almost_equal, + assert_raises) import skimage from skimage import data +from skimage._shared._warnings import expected_warnings from skimage.filters.thresholding import (threshold_adaptive, threshold_otsu, threshold_li, @@ -156,13 +159,15 @@ def test_otsu_coins_image_as_float(): assert 0.41 < threshold_otsu(coins) < 0.42 -def test_otsu_lena_image(): - img = skimage.img_as_ubyte(data.lena()) - assert 140 < threshold_otsu(img) < 142 - def test_otsu_astro_image(): img = skimage.img_as_ubyte(data.astronaut()) - assert 109 < threshold_otsu(img) < 111 + with expected_warnings(['grayscale']): + assert 109 < threshold_otsu(img) < 111 + + +def test_otsu_one_color_image(): + img = np.ones((10, 10), dtype=np.uint8) + assert_raises(TypeError, threshold_otsu, img) def test_li_camera_image(): camera = skimage.img_as_ubyte(data.camera()) @@ -198,6 +203,11 @@ def test_yen_coins_image_as_float(): assert 0.43 < threshold_yen(coins) < 0.44 +def test_adaptive_even_block_size_error(): + img = data.camera() + assert_raises(ValueError, threshold_adaptive, img, block_size=4) + + def test_isodata_camera_image(): camera = skimage.img_as_ubyte(data.camera()) diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 272641ae..331ea3f4 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -7,7 +7,7 @@ __all__ = ['threshold_adaptive', import numpy as np from scipy import ndimage as ndi from ..exposure import histogram -from .._shared.utils import assert_nD +from .._shared.utils import assert_nD, warn def threshold_adaptive(image, block_size, method='gaussian', offset=0, @@ -24,7 +24,7 @@ def threshold_adaptive(image, block_size, method='gaussian', offset=0, image : (N, M) ndarray Input image. block_size : int - Uneven size of pixel neighborhood which is used to calculate the + Odd size of pixel neighborhood which is used to calculate the threshold value (e.g. 3, 5, 7, ..., 21, ...). method : {'generic', 'gaussian', 'mean', 'median'}, optional Method used to determine adaptive threshold for local neighbourhood in @@ -67,6 +67,9 @@ def threshold_adaptive(image, block_size, method='gaussian', offset=0, >>> func = lambda arr: arr.mean() >>> binary_image2 = threshold_adaptive(image, 15, 'generic', param=func) """ + if block_size % 2 == 0: + raise ValueError("The kwarg ``block_size`` must be odd! Given " + "``block_size`` {0} is even.".format(block_size)) assert_nD(image, 2) thresh_image = np.zeros(image.shape, 'double') if method == 'generic': @@ -97,7 +100,7 @@ def threshold_otsu(image, nbins=256): Parameters ---------- image : array - Input image. + Grayscale input image. nbins : int, optional Number of bins used to calculate histogram. This value is ignored for integer arrays. @@ -118,7 +121,22 @@ def threshold_otsu(image, nbins=256): >>> image = camera() >>> thresh = threshold_otsu(image) >>> binary = image <= thresh + + Notes + ----- + The input image must be grayscale. """ + if image.shape[-1] in (3, 4): + msg = "threshold_otsu is expected to work correctly only for " \ + "grayscale images; image shape {0} looks like an RGB image" + warn(msg.format(image.shape)) + + # Check if the image is multi-colored or not + if image.min() == image.max(): + raise TypeError("threshold_otsu is expected to work with images " \ + "having more than one color. The input image seems " \ + "to have just one color {0}.".format(image.min())) + hist, bin_centers = histogram(image.ravel(), nbins) hist = hist.astype(float) diff --git a/skimage/future/graph/_ncut.py b/skimage/future/graph/_ncut.py index acdade69..faa8591d 100644 --- a/skimage/future/graph/_ncut.py +++ b/skimage/future/graph/_ncut.py @@ -1,8 +1,8 @@ try: import networkx as nx except ImportError: - import warnings - warnings.warn('RAGs require networkx') + from ..._shared.utils import warn + warn('RAGs require networkx') import numpy as np from scipy import sparse from . import _ncut_cy diff --git a/skimage/future/graph/graph_cut.py b/skimage/future/graph/graph_cut.py index 3b112b7e..3279bb7e 100644 --- a/skimage/future/graph/graph_cut.py +++ b/skimage/future/graph/graph_cut.py @@ -1,8 +1,9 @@ + try: import networkx as nx except ImportError: - import warnings - warnings.warn('RAGs require networkx') + from ..._shared.utils import warn + warn('RAGs require networkx') import numpy as np from . import _ncut from . import _ncut_cy diff --git a/skimage/future/graph/rag.py b/skimage/future/graph/rag.py index dae5bbd1..96f00b6f 100644 --- a/skimage/future/graph/rag.py +++ b/skimage/future/graph/rag.py @@ -143,6 +143,13 @@ class RAG(nx.Graph): if label_image is not None: fp = ndi.generate_binary_structure(label_image.ndim, connectivity) + # In the next ``ndi.generic_filter`` function, the kwarg + # ``output`` is used to provide a strided array with a single + # 64-bit floating point number, to which the function repeatedly + # writes. This is done because even if we don't care about the + # output, without this, a float array of the same shape as the + # input image will be created and that could be expensive in + # memory consumption. ndi.generic_filter( label_image, function=_add_edge_filter, diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 6174f831..047f12b4 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -36,7 +36,7 @@ THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. import cython import numpy as np import heap -import warnings +from .._shared.utils import warn cimport numpy as cnp cimport heap @@ -304,7 +304,8 @@ cdef class MCP: self.flat_costs = costs.astype(FLOAT_D, copy=False).ravel('F') except TypeError: self.flat_costs = costs.astype(FLOAT_D).flatten('F') - warnings.warn('Upgrading NumPy should decrease memory usage and increase speed.', Warning) + warn('Upgrading NumPy should decrease memory usage and increase' + ' speed.') size = self.flat_costs.shape[0] self.flat_cumulative_costs = np.empty(size, dtype=FLOAT_D) self.dim = len(costs.shape) diff --git a/skimage/io/_io.py b/skimage/io/_io.py index 614b0e23..fb3e6fbc 100644 --- a/skimage/io/_io.py +++ b/skimage/io/_io.py @@ -1,5 +1,4 @@ from io import BytesIO -import warnings import numpy as np import six @@ -8,7 +7,7 @@ from ..io.manage_plugins import call_plugin from ..color import rgb2grey from .util import file_or_url_context from ..exposure import is_low_contrast -from .._shared._warnings import all_warnings +from .._shared.utils import all_warnings, warn __all__ = ['imread', 'imread_collection', 'imsave', 'imshow', 'show'] @@ -129,7 +128,7 @@ def imsave(fname, arr, plugin=None, **plugin_args): if fname.lower().endswith(('.tiff', '.tif')): plugin = 'tifffile' if is_low_contrast(arr): - warnings.warn('%s is a low contrast image' % fname) + warn('%s is a low contrast image' % fname) return call_plugin('imsave', fname, arr, plugin=plugin, **plugin_args) diff --git a/skimage/io/_plugins/matplotlib_plugin.py b/skimage/io/_plugins/matplotlib_plugin.py index a9754a72..34afb7ff 100644 --- a/skimage/io/_plugins/matplotlib_plugin.py +++ b/skimage/io/_plugins/matplotlib_plugin.py @@ -1,10 +1,11 @@ from collections import namedtuple import numpy as np -import warnings import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable from ...util import dtype as dtypes from ...exposure import is_low_contrast from ...util.colormap import viridis +from ..._shared.utils import warn _default_colormap = 'gray' _nonstandard_colormap = viridis @@ -67,14 +68,14 @@ def _raise_warnings(image_properties): """ ip = image_properties if ip.unsupported_dtype: - warnings.warn("Non-standard image type; displaying image with " - "stretched contrast.") + warn("Non-standard image type; displaying image with " + "stretched contrast.") if ip.low_dynamic_range: - warnings.warn("Low image dynamic range; displaying image with " - "stretched contrast.") + warn("Low image dynamic range; displaying image with " + "stretched contrast.") if ip.out_of_range_float: - warnings.warn("Float image out of standard range; displaying " - "image with stretched contrast.") + warn("Float image out of standard range; displaying " + "image with stretched contrast.") def _get_display_range(image): @@ -110,7 +111,7 @@ def _get_display_range(image): return lo, hi, cmap -def imshow(im, *args, **kwargs): +def imshow(im, ax=None, show_cbar=None, **kwargs): """Show the input image and return the current axes. By default, the image is displayed in greyscale, rather than @@ -131,8 +132,11 @@ def imshow(im, *args, **kwargs): ---------- im : array, shape (M, N[, 3]) The image to display. - - *args, **kwargs : positional and keyword arguments + ax: `matplotlib.axes.Axes`, optional + The axis to use for the image, defaults to plt.gca(). + show_cbar: boolean, optional. + Whether to show the colorbar (used to override default behavior). + **kwargs : Keyword arguments These are passed directly to `matplotlib.pyplot.imshow`. Returns @@ -147,9 +151,15 @@ def imshow(im, *args, **kwargs): kwargs.setdefault('cmap', cmap) kwargs.setdefault('vmin', lo) kwargs.setdefault('vmax', hi) - ax_im = plt.imshow(im, *args, **kwargs) - if cmap != _default_colormap: - plt.colorbar() + + ax = ax or plt.gca() + ax_im = ax.imshow(im, **kwargs) + if (cmap != _default_colormap and show_cbar is not False) or show_cbar: + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + plt.colorbar(ax_im, cax=cax) + ax.set_adjustable('box-forced') + ax.get_figure().tight_layout() return ax_im imread = plt.imread diff --git a/skimage/io/_plugins/qt_plugin.py b/skimage/io/_plugins/qt_plugin.py index 4d061933..ad869851 100644 --- a/skimage/io/_plugins/qt_plugin.py +++ b/skimage/io/_plugins/qt_plugin.py @@ -1,3 +1,4 @@ +from ..._shared import warn from .util import prepare_for_display, window_manager import numpy as np @@ -10,7 +11,6 @@ try: QLabel, QMainWindow, QPixmap, QWidget) from PyQt4 import QtCore, QtGui import sip - import warnings except ImportError: window_manager._release('qt') @@ -119,8 +119,7 @@ if sip.SIP_VERSION >= 0x040c00: # doesn't work with earlier versions imread = imread_qt else: - warnings.warn(RuntimeWarning( - "sip version too old. QT imread disabled")) + warn(RuntimeWarning("sip version too old. QT imread disabled")) def imshow(arr, fancy=False): diff --git a/skimage/io/tests/test_mpl_imshow.py b/skimage/io/tests/test_mpl_imshow.py index e369273a..76422580 100644 --- a/skimage/io/tests/test_mpl_imshow.py +++ b/skimage/io/tests/test_mpl_imshow.py @@ -78,7 +78,11 @@ def test_low_dynamic_range(): def test_outside_standard_range(): plt.figure() - with expected_warnings(["out of standard range"]): + # Warning raised by matplotlib on Windows: + # "The CObject type is marked Pending Deprecation in Python 2.7. + # Please use capsule objects instead." + # Ref: https://docs.python.org/2/c-api/cobject.html + with expected_warnings(["out of standard range|CObject type is marked"]): ax_im = io.imshow(im_hi) assert ax_im.get_clim() == (im_hi.min(), im_hi.max()) assert n_subplots(ax_im) == 2 @@ -87,7 +91,11 @@ def test_outside_standard_range(): def test_nonstandard_type(): plt.figure() - with expected_warnings(["Low image dynamic range"]): + # Warning raised by matplotlib on Windows: + # "The CObject type is marked Pending Deprecation in Python 2.7. + # Please use capsule objects instead." + # Ref: https://docs.python.org/2/c-api/cobject.html + with expected_warnings(["Low image dynamic range|CObject type is marked"]): ax_im = io.imshow(im64) assert ax_im.get_clim() == (im64.min(), im64.max()) assert n_subplots(ax_im) == 2 diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 9e8df953..be5f2bd0 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -2,6 +2,7 @@ from ._find_contours import find_contours from ._marching_cubes import (marching_cubes, mesh_surface_area, correct_mesh_orientation) from ._regionprops import regionprops, perimeter +from .simple_metrics import mean_squared_error, normalized_root_mse, psnr from ._structural_similarity import structural_similarity from ._polygon import approximate_polygon, subdivide_polygon from ._pnpoly import points_in_poly, grid_points_in_poly @@ -34,4 +35,7 @@ __all__ = ['find_contours', 'profile_line', 'label', 'points_in_poly', - 'grid_points_in_poly'] + 'grid_points_in_poly', + 'mean_squared_error', + 'normalized_root_mse', + 'psnr'] diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 90c921cc..e7088d24 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -4,7 +4,7 @@ #cython: wraparound=False import numpy as np -import warnings +from .._shared.utils import warn cimport numpy as cnp @@ -47,7 +47,7 @@ ctypedef struct bginfo: cdef void get_bginfo(background_val, bginfo *ret) except *: if background_val is None: - warnings.warn(DeprecationWarning( + warn(DeprecationWarning( 'The default value for `background` will change to 0 in v0.12' )) ret.background_val = -1 diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index cd0e1b60..46d50c7c 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -1,5 +1,6 @@ import numpy as np import scipy.ndimage as ndi +from .._shared.utils import warn from . import _marching_cubes_cy @@ -239,11 +240,9 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), skimage.measure.mesh_surface_area """ - import warnings - warnings.warn( - DeprecationWarning("`correct_mesh_orientation` is deprecated for " - "removal as `marching_cubes` now guarantess " - "correct mesh orientation.")) + warn(DeprecationWarning("`correct_mesh_orientation` is deprecated for " + "removal as `marching_cubes` now guarantess " + "correct mesh orientation.")) verts = verts.copy() verts[:, 0] /= spacing[0] diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index b227d59f..360a5d6d 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -1,8 +1,7 @@ import math -import warnings import numpy as np from scipy import optimize -from .._shared.utils import skimage_deprecation +from .._shared.utils import skimage_deprecation, warn def _check_data_dim(data, dim): @@ -27,8 +26,7 @@ class BaseModel(object): @property def _params(self): - warnings.warn('`_params` attribute is deprecated, ' - 'use `params` instead.') + warn('`_params` attribute is deprecated, use `params` instead.') return self.params @@ -61,8 +59,8 @@ class LineModel(BaseModel): def __init__(self): self.params = None - warnings.warn(skimage_deprecation('`LineModel` is deprecated, ' - 'use `LineModelND` instead.')) + warn(skimage_deprecation('`LineModel` is deprecated, ' + 'use `LineModelND` instead.')) def estimate(self, data): """Estimate line model from data using total least squares. diff --git a/skimage/measure/simple_metrics.py b/skimage/measure/simple_metrics.py new file mode 100644 index 00000000..6083010b --- /dev/null +++ b/skimage/measure/simple_metrics.py @@ -0,0 +1,132 @@ +from __future__ import division + +import numpy as np +from ..util.dtype import dtype_range + +__all__ = ['mean_squared_error', 'normalized_root_mse', 'psnr'] + + +def _assert_compatible(im1, im2): + """Raise an error if the shape and dtype do not match.""" + if not im1.dtype == im2.dtype: + raise ValueError('Input images must have the same dtype.') + if not im1.shape == im2.shape: + raise ValueError('Input images must have the same dimensions.') + return + + +def _as_floats(im1, im2): + """Promote im1, im2 to nearest appropriate floating point precision.""" + float_type = np.result_type(im1.dtype, im2.dtype, np.float32) + if im1.dtype != float_type: + im1 = im1.astype(float_type) + if im2.dtype != float_type: + im2 = im2.astype(float_type) + return im1, im2 + + +def mean_squared_error(im1, im2): + """Compute the mean-squared error between two images. + + Parameters + ---------- + im1, im2 : ndarray + Image. Any dimensionality. + + Returns + ------- + mse : float + The mean-squared error (MSE) metric. + + """ + _assert_compatible(im1, im2) + im1, im2 = _as_floats(im1, im2) + return np.mean(np.square(im1 - im2), dtype=np.float64) + + +def normalized_root_mse(im_true, im_test, norm_type='Euclidean'): + """Compute the normalized root mean-squared error (NRMSE) between two + images. + + Parameters + ---------- + im_true : ndarray + Ground-truth image. + im_test : ndarray + Test image. + norm_type : {'Euclidean', 'min-max', 'mean'} + Controls the normalization method to use in the denominator of the + NRMSE. There is no standard method of normalization across the + literature [1]_. The methods available here are as follows: + + - 'Euclidean' : normalize by the Euclidean norm of ``im_true``. + - 'min-max' : normalize by the intensity range of ``im_true``. + - 'mean' : normalize by the mean of ``im_true``. + + Returns + ------- + nrmse : float + The NRMSE metric. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Root-mean-square_deviation + + """ + _assert_compatible(im_true, im_test) + im_true, im_test = _as_floats(im_true, im_test) + + norm_type = norm_type.lower() + if norm_type == 'euclidean': + denom = np.sqrt(np.mean((im_true*im_true), dtype=np.float64)) + elif norm_type == 'min-max': + denom = im_true.max() - im_true.min() + elif norm_type == 'mean': + denom = im_true.mean() + else: + raise ValueError("Unsupported norm_type") + return np.sqrt(mean_squared_error(im_true, im_test)) / denom + + +def psnr(im_true, im_test, dynamic_range=None): + """ Compute the peak signal to noise ratio (PSNR) for an image. + + Parameters + ---------- + im_true : ndarray + Ground-truth image. + im_test : ndarray + Test image. + dynamic_range : int + The dynamic range of the input image (distance between minimum and + maximum possible values). By default, this is estimated from the image + data-type. + + Returns + ------- + psnr : float + The PSNR metric. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio + + """ + _assert_compatible(im_true, im_test) + if dynamic_range is None: + dmin, dmax = dtype_range[im_true.dtype.type] + true_min, true_max = np.min(im_true), np.max(im_true) + if true_max > dmax or true_min < dmin: + raise ValueError( + "im_true has intensity values outside the range expected for " + "its data type. Please manually specify the dynamic_range") + if true_min >= 0: + # most common case (255 for uint8, 1 for float) + dynamic_range = dmax + else: + dynamic_range = dmax - dmin + + im_true, im_test = _as_floats(im_true, im_test) + + err = mean_squared_error(im_true, im_test) + return 10 * np.log10((dynamic_range ** 2) / err) diff --git a/skimage/measure/tests/test_simple_metrics.py b/skimage/measure/tests/test_simple_metrics.py new file mode 100644 index 00000000..b59a11a8 --- /dev/null +++ b/skimage/measure/tests/test_simple_metrics.py @@ -0,0 +1,61 @@ +import numpy as np +from numpy.testing import (run_module_suite, assert_equal, assert_raises, + assert_almost_equal) + +from skimage.measure import psnr, normalized_root_mse, mean_squared_error +import skimage.data + +np.random.seed(5) +cam = skimage.data.camera() +sigma = 20.0 +cam_noisy = np.clip(cam + sigma * np.random.randn(*cam.shape), 0, 255) +cam_noisy = cam_noisy.astype(cam.dtype) + + +def test_PSNR_vs_IPOL(): + # Tests vs. imdiff result from the following IPOL article and code: + # http://www.ipol.im/pub/art/2011/g_lmii/ + p_IPOL = 22.4497 + p = psnr(cam, cam_noisy) + assert_almost_equal(p, p_IPOL, decimal=4) + + +def test_PSNR_float(): + p_uint8 = psnr(cam, cam_noisy) + p_float64 = psnr(cam/255., cam_noisy/255., dynamic_range=1) + assert_almost_equal(p_uint8, p_float64, decimal=5) + + +def test_PSNR_errors(): + assert_raises(ValueError, psnr, cam, cam.astype(np.float32)) + assert_raises(ValueError, psnr, cam, cam[:-1, :]) + + +def test_NRMSE(): + x = np.ones(4) + y = np.asarray([0., 2., 2., 2.]) + assert_equal(normalized_root_mse(y, x, 'mean'), 1/np.mean(y)) + assert_equal(normalized_root_mse(y, x, 'Euclidean'), 1/np.sqrt(3)) + assert_equal(normalized_root_mse(y, x, 'min-max'), 1/(y.max()-y.min())) + + +def test_NRMSE_no_int_overflow(): + camf = cam.astype(np.float32) + cam_noisyf = cam_noisy.astype(np.float32) + assert_almost_equal(mean_squared_error(cam, cam_noisy), + mean_squared_error(camf, cam_noisyf)) + assert_almost_equal(normalized_root_mse(cam, cam_noisy), + normalized_root_mse(camf, cam_noisyf)) + + +def test_NRMSE_errors(): + x = np.ones(4) + assert_raises(ValueError, normalized_root_mse, + x.astype(np.uint8), x.astype(np.float32)) + assert_raises(ValueError, normalized_root_mse, x[:-1], x) + # invalid normalization name + assert_raises(ValueError, normalized_root_mse, x, x, 'foo') + + +if __name__ == "__main__": + run_module_suite() diff --git a/skimage/morphology/misc.py b/skimage/morphology/misc.py index 1a024f38..653c240d 100644 --- a/skimage/morphology/misc.py +++ b/skimage/morphology/misc.py @@ -1,7 +1,7 @@ import numpy as np import functools -import warnings from scipy import ndimage as ndi +from .._shared.utils import warn from .selem import _default_selem # Our function names don't exactly correspond to ndimages. @@ -37,8 +37,8 @@ def default_selem(func): return func(image, selem=selem, *args, **kwargs) return func_out - -def _check_dtype_supported(ar): + +def _check_dtype_supported(ar): # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7 if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)): raise TypeError("Only bool or integer image types are supported. " @@ -119,8 +119,8 @@ def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False): "`skimage.morphology.label`.") if len(component_sizes) == 2: - warnings.warn("Only one label was provided to `remove_small_objects`. " - "Did you mean to use a boolean array?") + warn("Only one label was provided to `remove_small_objects`. " + "Did you mean to use a boolean array?") too_small = component_sizes < min_size too_small_mask = too_small[ccs] @@ -181,35 +181,35 @@ def remove_small_holes(ar, min_size=64, connectivity=1, in_place=False): Notes ----- - If the array type is int, it is assumed that it contains already-labeled - objects. The labels are not kept in the output image (this function always - outputs a bool image). It is suggested that labeling is completed after + If the array type is int, it is assumed that it contains already-labeled + objects. The labels are not kept in the output image (this function always + outputs a bool image). It is suggested that labeling is completed after using this function. """ _check_dtype_supported(ar) - + #Creates warning if image is an integer image if ar.dtype != bool: - warnings.warn("Any labeled images will be returned as a boolean array. " - "Did you mean to use a boolean array?", UserWarning) - + warn("Any labeled images will be returned as a boolean array. " + "Did you mean to use a boolean array?", UserWarning) + if in_place: out = ar else: out = ar.copy() - + #Creating the inverse of ar if in_place: out = np.logical_not(out,out) else: out = np.logical_not(out) - + #removing small objects from the inverse of ar out = remove_small_objects(out, min_size, connectivity, in_place) - + if in_place: out = np.logical_not(out,out) else: out = np.logical_not(out) - + return out diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 97fab246..62985cc7 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -116,13 +116,13 @@ def denoise_tv_bregman(image, weight, max_iter=100, eps=1e-3, isotropic=True): return _denoise_tv_bregman(image, weight, max_iter, eps, isotropic) -def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): - """Perform total-variation denoising on 3D images. +def _denoise_tv_chambolle_nd(im, weight=0.1, eps=2.e-4, n_iter_max=200): + """Perform total-variation denoising on n-dimensional images. Parameters ---------- im : ndarray - 3-D input data to be denoised. + n-D input data to be denoised. weight : float, optional Denoising weight. The greater `weight`, the more denoising (at the expense of fidelity to `input`). @@ -146,36 +146,45 @@ def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): """ - px = np.zeros_like(im) - py = np.zeros_like(im) - pz = np.zeros_like(im) - gx = np.zeros_like(im) - gy = np.zeros_like(im) - gz = np.zeros_like(im) + ndim = im.ndim + p = np.zeros((im.ndim, ) + im.shape, dtype=im.dtype) + g = np.zeros_like(p) d = np.zeros_like(im) i = 0 while i < n_iter_max: - d = - px - py - pz - d[1:] += px[:-1] - d[:, 1:] += py[:, :-1] - d[:, :, 1:] += pz[:, :, :-1] - - out = im + d + if i > 0: + # d will be the (negative) divergence of p + d = -p.sum(0) + slices_d = [slice(None), ] * ndim + slices_p = [slice(None), ] * (ndim + 1) + for ax in range(ndim): + slices_d[ax] = slice(1, None) + slices_p[ax+1] = slice(0, -1) + slices_p[0] = ax + d[slices_d] += p[slices_p] + slices_d[ax] = slice(None) + slices_p[ax+1] = slice(None) + out = im + d + else: + out = im E = (d ** 2).sum() - gx[:-1] = np.diff(out, axis=0) - gy[:, :-1] = np.diff(out, axis=1) - gz[:, :, :-1] = np.diff(out, axis=2) - norm = np.sqrt(gx ** 2 + gy ** 2 + gz ** 2) + # g stores the gradients of out along each axis + # e.g. g[0] is the first order finite difference along axis 0 + slices_g = [slice(None), ] * (ndim + 1) + for ax in range(ndim): + slices_g[ax+1] = slice(0, -1) + slices_g[0] = ax + g[slices_g] = np.diff(out, axis=ax) + slices_g[ax+1] = slice(None) + + norm = np.sqrt((g ** 2).sum(axis=0))[np.newaxis, ...] E += weight * norm.sum() - norm *= 0.5 / weight + tau = 1. / (2.*ndim) + norm *= tau / weight norm += 1. - px -= 1. / 6. * gx - px /= norm - py -= 1. / 6. * gy - py /= norm - pz -= 1 / 6. * gz - pz /= norm + p -= tau * g + p /= norm E /= float(im.size) if i == 0: E_init = E @@ -189,89 +198,13 @@ def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): return out -def _denoise_tv_chambolle_2d(im, weight=0.2, eps=2.e-4, n_iter_max=200): - """Perform total-variation denoising on 2D images. - - Parameters - ---------- - im : ndarray - Input data to be denoised. - weight : float, optional - Denoising weight. The greater `weight`, the more denoising (at - the expense of fidelity to `input`) - eps : float, optional - Relative difference of the value of the cost function that determines - the stop criterion. The algorithm stops when: - - (E_(n-1) - E_n) < eps * E_0 - - n_iter_max : int, optional - Maximal number of iterations used for the optimization. - - Returns - ------- - out : ndarray - Denoised array of floats. - - Notes - ----- - The principle of total variation denoising is explained in - http://en.wikipedia.org/wiki/Total_variation_denoising. - - This code is an implementation of the algorithm of Rudin, Fatemi and Osher - that was proposed by Chambolle in [1]_. - - References - ---------- - .. [1] A. Chambolle, An algorithm for total variation minimization and - applications, Journal of Mathematical Imaging and Vision, - Springer, 2004, 20, 89-97. - - """ - - px = np.zeros_like(im) - py = np.zeros_like(im) - gx = np.zeros_like(im) - gy = np.zeros_like(im) - d = np.zeros_like(im) - i = 0 - while i < n_iter_max: - d = -px - py - d[1:] += px[:-1] - d[:, 1:] += py[:, :-1] - - out = im + d - E = (d ** 2).sum() - gx[:-1] = np.diff(out, axis=0) - gy[:, :-1] = np.diff(out, axis=1) - norm = np.sqrt(gx ** 2 + gy ** 2) - E += weight * norm.sum() - norm *= 0.5 / weight - norm += 1 - px -= 0.25 * gx - px /= norm - py -= 0.25 * gy - py /= norm - E /= float(im.size) - if i == 0: - E_init = E - E_previous = E - else: - if np.abs(E_previous - E) < eps * E_init: - break - else: - E_previous = E - i += 1 - return out - - -def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, +def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200, multichannel=False): - """Perform total-variation denoising on 2D and 3D images. + """Perform total-variation denoising on n-dimensional images. Parameters ---------- - im : ndarray (2d or 3d) of ints, uints or floats + im : ndarray of ints, uints or floats Input data to be denoised. `im` can be of any numeric type, but it is cast into an ndarray of floats for the computation of the denoised image. @@ -289,7 +222,7 @@ def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, multichannel : bool, optional Apply total-variation denoising separately for each channel. This option should be true for color images, otherwise the denoising is - also applied in the 3rd dimension. + also applied in the channels dimension. Returns ------- @@ -341,17 +274,11 @@ def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, if not im_type.kind == 'f': im = img_as_float(im) - if im.ndim == 2: - out = _denoise_tv_chambolle_2d(im, weight, eps, n_iter_max) - elif im.ndim == 3: - if multichannel: - out = np.zeros_like(im) - for c in range(im.shape[2]): - out[..., c] = _denoise_tv_chambolle_2d(im[..., c], weight, eps, - n_iter_max) - else: - out = _denoise_tv_chambolle_3d(im, weight, eps, n_iter_max) + if multichannel: + out = np.zeros_like(im) + for c in range(im.shape[-1]): + out[..., c] = _denoise_tv_chambolle_nd(im[..., c], weight, eps, + n_iter_max) else: - raise ValueError('only 2-d and 3-d images may be denoised with this ' - 'function') + out = _denoise_tv_chambolle_nd(im, weight, eps, n_iter_max) return out diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py new file mode 100644 index 00000000..bf696e04 --- /dev/null +++ b/skimage/restoration/inpaint.py @@ -0,0 +1,139 @@ +from __future__ import division + +import numpy as np +import skimage +from scipy import sparse +from scipy.sparse.linalg import spsolve +from scipy.ndimage.filters import laplace + + +def _get_neighborhood(nd_idx, radius, nd_shape): + bounds_lo = (nd_idx - radius).clip(min=0) + bounds_hi = (nd_idx + radius + 1).clip(max=nd_shape) + return bounds_lo, bounds_hi + + +def _inpaint_biharmonic_single_channel(img, mask, out, limits): + # Initialize sparse matrices + matrix_unknown = sparse.lil_matrix((np.sum(mask), out.size)) + matrix_known = sparse.lil_matrix((np.sum(mask), out.size)) + + # Find indexes of masked points in flatten array + mask_i = np.ravel_multi_index(np.where(mask), mask.shape) + + # Find masked points and prepare them to be easily enumerate over + mask_pts = np.array(np.where(mask)).T + + # Iterate over masked points + for mask_pt_n, mask_pt_idx in enumerate(mask_pts): + # Get bounded neighborhood of selected radius + b_lo, b_hi = _get_neighborhood(mask_pt_idx, 2, out.shape) + + # Create biharmonic coefficients ndarray + neigh_coef = np.zeros(b_hi - b_lo) + neigh_coef[tuple(mask_pt_idx - b_lo)] = 1 + neigh_coef = laplace(laplace(neigh_coef)) + + # Iterate over masked point's neighborhood + it_inner = np.nditer(neigh_coef, flags=['multi_index']) + for coef in it_inner: + if coef == 0: + continue + tmp_pt_idx = np.add(b_lo, it_inner.multi_index) + tmp_pt_i = np.ravel_multi_index(tmp_pt_idx, mask.shape) + + if mask[tuple(tmp_pt_idx)]: + matrix_unknown[mask_pt_n, tmp_pt_i] = coef + else: + matrix_known[mask_pt_n, tmp_pt_i] = coef + + # Prepare diagonal matrix + flat_diag_image = sparse.dia_matrix((out.flatten(), np.array([0])), + shape=(out.size, out.size)) + + # Calculate right hand side as a sum of known matrix's columns + matrix_known = matrix_known.tocsr() + rhs = -(matrix_known * flat_diag_image).sum(axis=1) + + # Solve linear system for masked points + matrix_unknown = matrix_unknown[:, mask_i] + matrix_unknown = sparse.csr_matrix(matrix_unknown) + result = spsolve(matrix_unknown, rhs) + + # Handle enormous values + result = np.clip(result, *limits) + + result = result.ravel() + + # Substitute masked points with inpainted versions + for mask_pt_n, mask_pt_idx in enumerate(mask_pts): + out[tuple(mask_pt_idx)] = result[mask_pt_n] + + return out + + +def inpaint_biharmonic(img, mask, multichannel=False): + """Inpaint masked points in image with biharmonic equations. + + Parameters + ---------- + img : (M[, N[, ..., P]][, C]) ndarray + Input image. + mask : (M[, N[, ..., P]]) ndarray + Array of pixels to be inpainted. Have to be the same shape as one + of the 'img' channels. Unknown pixels have to be represented with 1, + known pixels - with 0. + multichannel : boolean, optional + If True, the last `img` dimension is considered as a color channel, + otherwise as spatial. + + Returns + ------- + out : (M[, N[, ..., P]][, C]) ndarray + Input image with masked pixels inpainted. + + Example + ------- + >>> img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1)) + >>> mask = np.zeros_like(img) + >>> mask[2, 2:] = 1 + >>> mask[1, 3:] = 1 + >>> mask[0, 4:] = 1 + >>> out = inpaint_biharmonic(img, mask) + + References + ---------- + Algorithm is based on: + .. [1] N.S.Hoang, S.B.Damelin, "On surface completion and image inpainting + by biharmonic functions: numerical aspects", + http://www.ima.umn.edu/~damelin/biharmonic + """ + + if img.ndim < 1: + raise ValueError('Input array has to be at least 1D') + + img_baseshape = img.shape[:-1] if multichannel else img.shape + if img_baseshape != mask.shape: + raise ValueError('Input arrays have to be the same shape') + + if np.ma.isMaskedArray(img): + raise TypeError('Masked arrays are not supported') + + img = skimage.img_as_float(img) + mask = mask.astype(np.bool) + + if not multichannel: + img = img[..., np.newaxis] + + out = np.copy(img) + + for i in range(img.shape[-1]): + known_points = img[..., i][~mask] + limits = (np.min(known_points), np.max(known_points)) + _inpaint_biharmonic_single_channel(img[..., i], mask, + out[..., i], limits) + + if not multichannel: + out = out[..., 0] + + return out diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index fb199afd..58f0a261 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -1,7 +1,7 @@ import numpy as np from numpy.testing import run_module_suite, assert_raises, assert_equal -from skimage import restoration, data, color, img_as_float +from skimage import restoration, data, color, img_as_float, measure np.random.seed(1234) @@ -21,7 +21,7 @@ def test_denoise_tv_chambolle_2d(): # clip noise so that it does not exceed allowed range for float images. img = np.clip(img, 0, 1) # denoise - denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.25) + denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.1) # which dtype? assert denoised_astro.dtype in [np.float, np.float32, np.float64] from scipy import ndimage as ndi @@ -34,8 +34,17 @@ def test_denoise_tv_chambolle_2d(): def test_denoise_tv_chambolle_multichannel(): - denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=0.25) - denoised = restoration.denoise_tv_chambolle(astro, weight=0.25, + denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=0.1) + denoised = restoration.denoise_tv_chambolle(astro, weight=0.1, + multichannel=True) + assert_equal(denoised[..., 0], denoised0) + + # tile astronaut subset to generate 3D+channels data + astro3 = np.tile(astro[:64, :64, np.newaxis, :], [1, 1, 2, 1]) + # modify along tiled dimension to give non-zero gradient on 3rd axis + astro3[:, :, 0, :] = 2*astro3[:, :, 0, :] + denoised0 = restoration.denoise_tv_chambolle(astro3[..., 0], weight=0.1) + denoised = restoration.denoise_tv_chambolle(astro3, weight=0.1, multichannel=True) assert_equal(denoised[..., 0], denoised0) @@ -46,7 +55,7 @@ def test_denoise_tv_chambolle_float_result_range(): int_astro = np.multiply(img, 255).astype(np.uint8) assert np.max(int_astro) > 1 denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, - weight=0.25) + weight=0.1) # test if the value range of output float data is within [0.0:1.0] assert denoised_int_astro.dtype == np.float assert np.max(denoised_int_astro) <= 1.0 @@ -62,13 +71,45 @@ def test_denoise_tv_chambolle_3d(): mask += 20 * np.random.rand(*mask.shape) mask[mask < 0] = 0 mask[mask > 255] = 255 - res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=0.4) + res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=0.1) assert res.dtype == np.float assert res.std() * 255 < mask.std() - # test wrong number of dimensions - assert_raises(ValueError, restoration.denoise_tv_chambolle, - np.random.rand(8, 8, 8, 8)) + +def test_denoise_tv_chambolle_1d(): + """Apply the TV denoising algorithm on a 1D sinusoid.""" + x = 125 + 100*np.sin(np.linspace(0, 8*np.pi, 1000)) + x += 20 * np.random.rand(x.size) + x = np.clip(x, 0, 255) + res = restoration.denoise_tv_chambolle(x.astype(np.uint8), weight=0.1) + assert res.dtype == np.float + assert res.std() * 255 < x.std() + + +def test_denoise_tv_chambolle_4d(): + """ TV denoising for a 4D input.""" + im = 255 * np.random.rand(8, 8, 8, 8) + res = restoration.denoise_tv_chambolle(im.astype(np.uint8), weight=0.1) + assert res.dtype == np.float + assert res.std() * 255 < im.std() + + +def test_denoise_tv_chambolle_weighting(): + # make sure a specified weight gives consistent results regardless of + # the number of input image dimensions + rstate = np.random.RandomState(1234) + img2d = astro_gray.copy() + img2d += 0.15 * rstate.standard_normal(img2d.shape) + img2d = np.clip(img2d, 0, 1) + + # generate 4D image by tiling + img4d = np.tile(img2d[..., None, None], (1, 1, 2, 2)) + + w = 0.2 + denoised_2d = restoration.denoise_tv_chambolle(img2d, weight=w) + denoised_4d = restoration.denoise_tv_chambolle(img4d, weight=w) + assert measure.structural_similarity(denoised_2d, + denoised_4d[:, :, 0, 0]) > 0.99 def test_denoise_tv_bregman_2d(): diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py new file mode 100644 index 00000000..04e8c8c2 --- /dev/null +++ b/skimage/restoration/tests/test_inpaint.py @@ -0,0 +1,66 @@ +from __future__ import print_function, division + +import numpy as np +from numpy.testing import (run_module_suite, assert_allclose, + assert_raises) +from skimage.restoration import inpaint + + +def test_inpaint_biharmonic_2d(): + img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1)) + mask = np.zeros_like(img) + mask[2, 2:] = 1 + mask[1, 3:] = 1 + mask[0, 4:] = 1 + img[np.where(mask)] = 0 + out = inpaint.inpaint_biharmonic(img, mask) + ref = np.array( + [[0., 0.0625, 0.25000000, 0.5625000, 0.73925058], + [0., 0.0625, 0.25000000, 0.5478048, 0.76557821], + [0., 0.0625, 0.25842878, 0.5623079, 0.85927796], + [0., 0.0625, 0.25000000, 0.5625000, 1.00000000], + [0., 0.0625, 0.25000000, 0.5625000, 1.00000000]] + ) + assert_allclose(ref, out) + + +def test_inpaint_biharmonic_3d(): + img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1)) + img = np.dstack((img, img.T)) + mask = np.zeros_like(img) + mask[2, 2:, :] = 1 + mask[1, 3:, :] = 1 + mask[0, 4:, :] = 1 + img[np.where(mask)] = 0 + out = inpaint.inpaint_biharmonic(img, mask) + ref = np.dstack(( + np.array( + [[0.0000, 0.0625, 0.25000000, 0.56250000, 0.53752796], + [0.0000, 0.0625, 0.25000000, 0.44443780, 0.53762210], + [0.0000, 0.0625, 0.23693666, 0.46621112, 0.68615592], + [0.0000, 0.0625, 0.25000000, 0.56250000, 1.00000000], + [0.0000, 0.0625, 0.25000000, 0.56250000, 1.00000000]]), + np.array( + [[0.0000, 0.0000, 0.00000000, 0.00000000, 0.19621902], + [0.0625, 0.0625, 0.06250000, 0.17470756, 0.30140091], + [0.2500, 0.2500, 0.27241289, 0.35155440, 0.43068654], + [0.5625, 0.5625, 0.56250000, 0.56250000, 0.56250000], + [1.0000, 1.0000, 1.00000000, 1.00000000, 1.00000000]]) + )) + assert_allclose(ref, out) + + +def test_invalid_input(): + img, mask = np.zeros([]), np.zeros([]) + assert_raises(ValueError, inpaint.inpaint_biharmonic, img, mask) + + img, mask = np.zeros((2, 2)), np.zeros((4, 1)) + assert_raises(ValueError, inpaint.inpaint_biharmonic, img, mask) + + img = np.ma.array(np.zeros((2, 2)), mask=[[0, 0], [0, 0]]) + mask = np.zeros((2, 2)) + assert_raises(TypeError, inpaint.inpaint_biharmonic, img, mask) + + +if __name__ == '__main__': + run_module_suite() diff --git a/skimage/restoration/unwrap.py b/skimage/restoration/unwrap.py index 14f1eafa..7c21edb7 100644 --- a/skimage/restoration/unwrap.py +++ b/skimage/restoration/unwrap.py @@ -1,7 +1,8 @@ import numpy as np -import warnings from six import string_types +from .._shared.utils import warn + from ._unwrap_1d import unwrap_1d from ._unwrap_2d import unwrap_2d from ._unwrap_3d import unwrap_3d @@ -83,9 +84,9 @@ def unwrap_phase(image, wrap_around=False, seed=None): if wrap_around[0]: raise ValueError('`wrap_around` is not supported for 1D images') if image.ndim in (2, 3) and 1 in image.shape: - warnings.warn('Image has a length 1 dimension. Consider using an ' - 'array of lower dimensionality to use a more efficient ' - 'algorithm') + warn('Image has a length 1 dimension. Consider using an ' + 'array of lower dimensionality to use a more efficient ' + 'algorithm') if np.ma.isMaskedArray(image): mask = np.require(np.ma.getmaskarray(image), np.uint8, ['C']) diff --git a/skimage/segmentation/_felzenszwalb.py b/skimage/segmentation/_felzenszwalb.py index d9b82c2b..7a3e7758 100644 --- a/skimage/segmentation/_felzenszwalb.py +++ b/skimage/segmentation/_felzenszwalb.py @@ -1,6 +1,6 @@ -import warnings import numpy as np +from .._shared.utils import warn from ._felzenszwalb_cy import _felzenszwalb_grey @@ -56,8 +56,8 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20): # assume we got 2d image with multiple channels n_channels = image.shape[2] if n_channels != 3: - warnings.warn("Got image with %d channels. Is that really what you" - " wanted?" % image.shape[2]) + warn("Got image with %d channels. Is that really what you" + " wanted?" % image.shape[2]) segmentations = [] # compute quickshift for each channel for c in range(n_channels): diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index d3a6a759..7aee8088 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -8,10 +8,12 @@ Installing pyamg and using the 'cg_mg' mode of random_walker improves significantly the performance. """ -import warnings import numpy as np from scipy import sparse, ndimage as ndi +from .._shared.utils import warn + + # executive summary for next code block: try to import umfpack from # scipy, but make sure not to raise a fuss if it fails since it's only # needed to speed up a few cases. @@ -345,17 +347,17 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, mode = 'bf' if UmfpackContext is None and mode == 'cg': - warnings.warn('"cg" mode will be used, but it may be slower than ' - '"bf" because SciPy was built without UMFPACK. Consider' - ' rebuilding SciPy with UMFPACK; this will greatly ' - 'accelerate the conjugate gradient ("cg") solver. ' - 'You may also install pyamg and run the random_walker ' - 'function in "cg_mg" mode (see docstring).') + warn('"cg" mode will be used, but it may be slower than ' + '"bf" because SciPy was built without UMFPACK. Consider' + ' rebuilding SciPy with UMFPACK; this will greatly ' + 'accelerate the conjugate gradient ("cg") solver. ' + 'You may also install pyamg and run the random_walker ' + 'function in "cg_mg" mode (see docstring).') if (labels != 0).all(): - warnings.warn('Random walker only segments unlabeled areas, where ' - 'labels == 0. No zero valued areas in labels were ' - 'found. Returning provided labels.') + warn('Random walker only segments unlabeled areas, where ' + 'labels == 0. No zero valued areas in labels were ' + 'found. Returning provided labels.') if return_full_prob: # Find and iterate over valid labels @@ -438,8 +440,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, return_full_prob=return_full_prob) if mode == 'cg_mg': if not amg_loaded: - warnings.warn( - """pyamg (http://pyamg.org/)) is needed to use + warn("""pyamg (http://pyamg.org/)) is needed to use this mode, but is not installed. The 'cg' mode will be used instead.""") X = _solve_cg(lap_sparse, B, tol=tol, diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py index adc21a65..b558f78a 100644 --- a/skimage/segmentation/slic_superpixels.py +++ b/skimage/segmentation/slic_superpixels.py @@ -3,8 +3,8 @@ import collections as coll import numpy as np from scipy import ndimage as ndi -import warnings +from .._shared.utils import warn from ..util import img_as_float, regular_grid from ..segmentation._slic import (_slic_cython, _enforce_label_connectivity_cython) @@ -111,8 +111,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0, """ if enforce_connectivity is None: - warnings.warn('Deprecation: enforce_connectivity will default to' - ' True in future versions.') + warn('Deprecation: enforce_connectivity will default to' + ' True in future versions.') enforce_connectivity = False image = img_as_float(image) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c9249444..1e06e8b5 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -1,12 +1,11 @@ import six import math -import warnings import numpy as np from scipy import spatial from scipy import ndimage as ndi from .._shared.utils import (get_bound_method_class, safe_as_int, - _mode_deprecations) + _mode_deprecations, warn) from ..util import img_as_float from ._warps_cy import _warp_fast @@ -184,8 +183,7 @@ class ProjectiveTransform(GeometricTransform): @property def _matrix(self): - warnings.warn('`_matrix` attribute is deprecated, ' - 'use `params` instead.') + warn('`_matrix` attribute is deprecated, use `params` instead.') return self.params @property @@ -782,8 +780,7 @@ class PolynomialTransform(GeometricTransform): @property def _params(self): - warnings.warn('`_params` attribute is deprecated, ' - 'use `params` instead.') + warn('`_params` attribute is deprecated, use `params` instead.') return self.params def estimate(self, src, dst, order=2): @@ -1320,13 +1317,13 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, if order == 2: # When fixing this issue, make sure to fix the branches further # below in this function - warnings.warn("Bi-quadratic interpolation behavior has changed due " - "to a bug in the implementation of scikit-image. " - "The new version now serves as a wrapper " - "around SciPy's interpolation functions, which itself " - "is not verified to be a correct implementation. Until " - "skimage's implementation is fixed, we recommend " - "to use bi-linear or bi-cubic interpolation instead.") + warn("Bi-quadratic interpolation behavior has changed due " + "to a bug in the implementation of scikit-image. " + "The new version now serves as a wrapper " + "around SciPy's interpolation functions, which itself " + "is not verified to be a correct implementation. Until " + "skimage's implementation is fixed, we recommend " + "to use bi-linear or bi-cubic interpolation instead.") if order in (0, 1, 3) and not map_args: # use fast Cython version for specific interpolation orders and input diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 2f2cf553..1979e8d4 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -251,16 +251,21 @@ def rotate(image, angle, resize=False, center=None, order=1, mode='constant', center = np.array((cols, rows)) / 2. - 0.5 else: center = np.asarray(center) - tform1 = SimilarityTransform(translation=-center) + tform1 = SimilarityTransform(translation=center) tform2 = SimilarityTransform(rotation=np.deg2rad(angle)) - tform3 = SimilarityTransform(translation=center) - tform = tform1 + tform2 + tform3 + tform3 = SimilarityTransform(translation=-center) + tform = tform3 + tform2 + tform1 output_shape = None if resize: # determine shape of output image - corners = np.array([[1, 1], [1, rows], [cols, rows], [cols, 1]]) - corners = tform(corners - 1) + corners = np.array([ + [0, 0], + [0, rows - 1], + [cols - 1, rows - 1], + [cols - 1, 0] + ]) + corners = tform.inverse(corners) minc = corners[:, 0].min() minr = corners[:, 1].min() maxc = corners[:, 0].max() @@ -270,7 +275,7 @@ def rotate(image, angle, resize=False, center=None, order=1, mode='constant', output_shape = np.ceil((out_rows, out_cols)) # fit output image in new shape - translation = ((cols - out_cols) / 2., (rows - out_rows) / 2.) + translation = (minc, minr) tform4 = SimilarityTransform(translation=translation) tform = tform4 + tform diff --git a/skimage/transform/integral.py b/skimage/transform/integral.py index 91b07ac7..2ec7ab65 100644 --- a/skimage/transform/integral.py +++ b/skimage/transform/integral.py @@ -1,6 +1,8 @@ import numpy as np import collections -import warnings + +from .._shared.utils import warn + def integral_image(img): """Integral image / summed area table. @@ -81,10 +83,10 @@ def integrate(ii, start, end, *args): rows = start.shape[0] # handle deprecated input format else: - warnings.warn("The syntax 'integrate(ii, r0, c0, r1, c1)' is " - "deprecated, and will be phased out in release 0.14. " - "The new syntax is " - "'integrate(ii, (r0, c0), (r1, c1))'.") + warn("The syntax 'integrate(ii, r0, c0, r1, c1)' is " + "deprecated, and will be phased out in release 0.14. " + "The new syntax is " + "'integrate(ii, (r0, c0), (r1, c1))'.") if isinstance(start, collections.Iterable): rows = len(start) args = (start, end) + args diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index a8dfcb87..95d1eb59 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -132,6 +132,20 @@ def test_rotate_center(): assert_almost_equal(x0, x) +def test_rotate_resize_center(): + x = np.zeros((10, 10), dtype=np.double) + x[0, 0] = 1 + + ref_x45 = np.zeros((14, 14), dtype=np.double) + ref_x45[6, 0] = 1 + ref_x45[7, 0] = 1 + + x45 = rotate(x, 45, resize=True, center=(3, 3), order=0) + # new dimension should be d = sqrt(2 * (10/2)^2) + assert x45.shape == (14, 14) + assert_equal(x45, ref_x45) + + def test_rescale(): # same scale factor x = np.zeros((5, 5), dtype=np.double) diff --git a/skimage/viewer/__init__.py b/skimage/viewer/__init__.py index 443e5659..f4b4c8ec 100644 --- a/skimage/viewer/__init__.py +++ b/skimage/viewer/__init__.py @@ -1,6 +1,6 @@ -import warnings +from .._shared.utils import warn from .viewers import ImageViewer, CollectionViewer from .qt import has_qt if not has_qt: - warnings.warn('Viewer requires Qt') + warn('Viewer requires Qt') diff --git a/skimage/viewer/utils/core.py b/skimage/viewer/utils/core.py index 7f632d43..c44f0991 100644 --- a/skimage/viewer/utils/core.py +++ b/skimage/viewer/utils/core.py @@ -1,15 +1,14 @@ -import warnings - import numpy as np from ..qt import QtWidgets, has_qt, FigureManagerQT, FigureCanvasQTAgg +from ..._shared.utils import warn import matplotlib as mpl from matplotlib.figure import Figure from matplotlib import _pylab_helpers from matplotlib.colors import LinearSegmentedColormap if has_qt and 'agg' not in mpl.get_backend().lower(): - warnings.warn("Recommended matplotlib backend is `Agg` for full " - "skimage.viewer functionality.") + warn("Recommended matplotlib backend is `Agg` for full " + "skimage.viewer functionality.") __all__ = ['init_qtapp', 'start_qtapp', 'RequiredAttr', 'figimage',