diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 5900487c..ebed896a 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -132,9 +132,10 @@ Dense DAISY feature description, circle perimeter drawing. - François Boulogne - Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, + Drawing: Andres Method for circle perimeter, ellipse perimeter, Bezier curve, anti-aliasing. Circular and elliptical Hough Transforms + Thresholding Various fixes - Thouis Jones diff --git a/doc/examples/filters/plot_threshold_minimum.py b/doc/examples/filters/plot_threshold_minimum.py deleted file mode 100644 index 950d34d2..00000000 --- a/doc/examples/filters/plot_threshold_minimum.py +++ /dev/null @@ -1,36 +0,0 @@ -""" -================================== -Minimum Algorithm For Thresholding -================================== - -The minimum algorithm takes a histogram of the image and smooths it -repeatedly until there are only two peaks in the histogram. Then it -finds the minimum value between the two peaks. After smoothing the -histogram, there can be multiple pixel values with the minimum histogram -count, so you can pick the 'min', 'mid', or 'max' of these values. - -""" -import matplotlib.pyplot as plt - -from skimage import data -from skimage.filters.thresholding import threshold_minimum - -image = data.camera() - -threshold = threshold_minimum(image, bias='min') -binarized = image > threshold - -fig, axes = plt.subplots(nrows=2, figsize=(7, 8)) -ax0, ax1 = axes -plt.gray() - -ax0.imshow(image) -ax0.set_title('Original image') - -ax1.imshow(binarized) -ax1.set_title('Result') - -for ax in axes: - ax.axis('off') - -plt.show() diff --git a/doc/examples/segmentation/plot_local_otsu.py b/doc/examples/segmentation/plot_local_otsu.py deleted file mode 100644 index e51ef5e9..00000000 --- a/doc/examples/segmentation/plot_local_otsu.py +++ /dev/null @@ -1,60 +0,0 @@ -""" -==================== -Local Otsu Threshold -==================== - -This example shows how Otsu's threshold [1]_ method can be applied locally. For -each pixel, an "optimal" threshold is determined by maximizing the variance -between two classes of pixels of the local neighborhood defined by a -structuring element. - -The example compares the local threshold with the global threshold. - -.. note: local is much slower than global thresholding - -.. [1] http://en.wikipedia.org/wiki/Otsu's_method - -""" - -from skimage import data -from skimage.morphology import disk -from skimage.filters import threshold_otsu, rank -from skimage.util import img_as_ubyte - -import matplotlib -import matplotlib.pyplot as plt - -matplotlib.rcParams['font.size'] = 9 -img = img_as_ubyte(data.page()) - -radius = 15 -selem = disk(radius) - -local_otsu = rank.otsu(img, selem) -threshold_global_otsu = threshold_otsu(img) -global_otsu = img >= threshold_global_otsu - -fig, ax = plt.subplots(2, 2, figsize=(8, 5), sharex=True, sharey=True, - subplot_kw={'adjustable': 'box-forced'}) -ax0, ax1, ax2, ax3 = ax.ravel() -plt.tight_layout() - -fig.colorbar(ax0.imshow(img, cmap=plt.cm.gray), - ax=ax0, orientation='horizontal') -ax0.set_title('Original') -ax0.axis('off') - -fig.colorbar(ax1.imshow(local_otsu, cmap=plt.cm.gray), - ax=ax1, orientation='horizontal') -ax1.set_title('Local Otsu (radius=%d)' % radius) -ax1.axis('off') - -ax2.imshow(img >= local_otsu, cmap=plt.cm.gray) -ax2.set_title('Original >= Local Otsu' % threshold_global_otsu) -ax2.axis('off') - -ax3.imshow(global_otsu, cmap=plt.cm.gray) -ax3.set_title('Global Otsu (threshold = %d)' % threshold_global_otsu) -ax3.axis('off') - -plt.show() diff --git a/doc/examples/segmentation/plot_otsu.py b/doc/examples/segmentation/plot_otsu.py deleted file mode 100644 index c279a65b..00000000 --- a/doc/examples/segmentation/plot_otsu.py +++ /dev/null @@ -1,49 +0,0 @@ -""" -============ -Thresholding -============ - -Thresholding is used to create a binary image. This example uses Otsu's method -to calculate the threshold value. - -Otsu's method calculates an "optimal" threshold (marked by a red line in the -histogram below) by maximizing the variance between two classes of pixels, -which are separated by the threshold. Equivalently, this threshold minimizes -the intra-class variance. - -.. [1] http://en.wikipedia.org/wiki/Otsu's_method - -""" -import matplotlib -import matplotlib.pyplot as plt - -from skimage.data import camera -from skimage.filters import threshold_otsu - - -matplotlib.rcParams['font.size'] = 9 - - -image = camera() -thresh = threshold_otsu(image) -binary = image > thresh - -#fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5)) -fig = plt.figure(figsize=(8, 2.5)) -ax1 = plt.subplot(1, 3, 1, adjustable='box-forced') -ax2 = plt.subplot(1, 3, 2) -ax3 = plt.subplot(1, 3, 3, sharex=ax1, sharey=ax1, adjustable='box-forced') - -ax1.imshow(image, cmap=plt.cm.gray) -ax1.set_title('Original') -ax1.axis('off') - -ax2.hist(image) -ax2.set_title('Histogram') -ax2.axvline(thresh, color='r') - -ax3.imshow(binary, cmap=plt.cm.gray) -ax3.set_title('Thresholded') -ax3.axis('off') - -plt.show() diff --git a/doc/examples/segmentation/plot_threshold_adaptive.py b/doc/examples/segmentation/plot_threshold_adaptive.py deleted file mode 100644 index 6f473abb..00000000 --- a/doc/examples/segmentation/plot_threshold_adaptive.py +++ /dev/null @@ -1,48 +0,0 @@ -""" -===================== -Adaptive Thresholding -===================== - -Thresholding is the simplest way to segment objects from a background. If that -background is relatively uniform, then you can use a global threshold value to -binarize the image by pixel-intensity. If there's large variation in the -background intensity, however, adaptive thresholding (a.k.a. local or dynamic -thresholding) may produce better results. - -Here, we binarize an image using the `threshold_adaptive` function, which -calculates thresholds in regions of size `block_size` surrounding each pixel -(i.e. local neighborhoods). Each threshold value is the weighted mean of the -local neighborhood minus an offset value. - -""" -import matplotlib.pyplot as plt - -from skimage import data -from skimage.filters import threshold_otsu, threshold_adaptive - - -image = data.page() - -global_thresh = threshold_otsu(image) -binary_global = image > global_thresh - -block_size = 35 -binary_adaptive = threshold_adaptive(image, block_size, offset=10) - -fig, axes = plt.subplots(nrows=3, figsize=(7, 8)) -ax0, ax1, ax2 = axes -plt.gray() - -ax0.imshow(image) -ax0.set_title('Image') - -ax1.imshow(binary_global) -ax1.set_title('Global thresholding') - -ax2.imshow(binary_adaptive) -ax2.set_title('Adaptive thresholding') - -for ax in axes: - ax.axis('off') - -plt.show() diff --git a/doc/examples/segmentation/plot_thresholding.py b/doc/examples/segmentation/plot_thresholding.py new file mode 100644 index 00000000..9a3ccb0d --- /dev/null +++ b/doc/examples/segmentation/plot_thresholding.py @@ -0,0 +1,73 @@ +""" +============ +Thresholding +============ + +Thresholding is used to create a binary image from a grayscale image [1]_. + +.. [1] https://en.wikipedia.org/wiki/Thresholding_%28image_processing%29 + +.. seealso:: + A more comprehensive presentation on + :ref:`sphx_glr_auto_examples_xx_applications_plot_thresholding.py` + +""" + +###################################################################### +# We illustrate how to apply one of these thresholding algorithms. +# Otsu's method [2]_ calculates an "optimal" threshold (marked by a red line in the +# histogram below) by maximizing the variance between two classes of pixels, +# which are separated by the threshold. Equivalently, this threshold minimizes +# the intra-class variance. +# +# .. [2] http://en.wikipedia.org/wiki/Otsu's_method +# + +import matplotlib.pyplot as plt +from skimage import data +from skimage.filters import threshold_otsu + + +image = data.camera() +thresh = threshold_otsu(image) +binary = image > thresh + +fig, axes = plt.subplots(ncols=3, figsize=(8, 2.5)) +ax = axes.ravel() +ax[0] = plt.subplot(1, 3, 1, adjustable='box-forced') +ax[1] = plt.subplot(1, 3, 2) +ax[2] = plt.subplot(1, 3, 3, sharex=ax[0], sharey=ax[0], adjustable='box-forced') + +ax[0].imshow(image, cmap=plt.cm.gray) +ax[0].set_title('Original') +ax[0].axis('off') + +ax[1].hist(image.ravel(), bins=256) +ax[1].set_title('Histogram') +ax[1].axvline(thresh, color='r') + +ax[2].imshow(binary, cmap=plt.cm.gray) +ax[2].set_title('Thresholded') +ax[2].axis('off') + +plt.show() + + +###################################################################### +# If you are not familiar with the details of the different algorithms and the +# underlying assumptions, it is often difficult to know which algorithm will give +# the best results. Therefore, Scikit-image includes a function to evaluate +# thresholding algorithms provided by the library. At a glance, you can select +# the best algorithm for you data without a deep understanding of their +# mechanisms. +# + +from skimage.filters import try_all_threshold + +img = data.page() + +# Here, we specify a radius for local thresholding algorithms. +# If it is not specified, only global algorithms are called. +fig, ax = try_all_threshold(img, radius=20, + figsize=(10, 8), verbose=False) +plt.show() diff --git a/doc/examples/xx_applications/plot_thresholding.py b/doc/examples/xx_applications/plot_thresholding.py new file mode 100644 index 00000000..92e4556b --- /dev/null +++ b/doc/examples/xx_applications/plot_thresholding.py @@ -0,0 +1,256 @@ +""" +============ +Thresholding +============ + +Thresholding is used to create a binary image from a grayscale image [1]_. +It is the simplest way to segment objects from a background. + +Thresholding algorithms implemented in scikit-image can be separated in two +categories: + +- Histogram-based. The histogram of the pixels' intensity is used and + certain assumptions are made on the properties of this histogram (e.g. bimodal). +- Local. To process a pixel, only the neighboring pixels are used. + These algorithms often require more computation time. + + +If you are not familiar with the details of the different algorithms and the +underlying assumptions, it is often difficult to know which algorithm will give +the best results. Therefore, Scikit-image includes a function to evaluate +thresholding algorithms provided by the library. At a glance, you can select +the best algorithm for you data without a deep understanding of their +mechanisms. + +.. [1] https://en.wikipedia.org/wiki/Thresholding_%28image_processing%29 + +.. seealso:: + Presentation on + :ref:`sphx_glr_auto_examples_xx_applications_plot_rank_filters.py`. +""" +import matplotlib +import matplotlib.pyplot as plt + +from skimage import data +from skimage.filters import try_all_threshold + +img = data.page() + +# Here, we specify a radius for local thresholding algorithms. +# If it is not specified, only global algorithms are called. +fig, ax = try_all_threshold(img, radius=20, + figsize=(10, 8), verbose=False) +plt.show() + +###################################################################### +# How to apply a threshold? +# ========================= +# +# Now, we illustrate how to apply one of these thresholding algorithms. +# This example uses the mean value of pixel intensities. It is a simple +# and naive threshold value, which is sometimes used as a guess value. +# + +from skimage.filters import threshold_mean + + +image = data.camera() +thresh = threshold_mean(image) +binary = image > thresh + +fig, axes = plt.subplots(ncols=2, figsize=(8, 3)) +ax = axes.ravel() + +ax[0].imshow(image, cmap=plt.cm.gray) +ax[0].set_title('Original image') + +ax[1].imshow(binary, cmap=plt.cm.gray) +ax[1].set_title('Result') + +for a in ax: + a.axis('off') + +plt.show() + +###################################################################### +# Bimodal histogram +# ================= +# +# For pictures with a bimodal histogram, more specific algorithms can be used. +# For instance, the minimum algorithm takes a histogram of the image and smooths it +# repeatedly until there are only two peaks in the histogram. Then it +# finds the minimum value between the two peaks. After smoothing the +# histogram, there can be multiple pixel values with the minimum histogram +# count, so you can pick the 'min', 'mid', or 'max' of these values. +# + +from skimage.filters import threshold_minimum + + +image = data.camera() + +thresh_min = threshold_minimum(image, bias='min') +binary_min = image > thresh_min +thresh_mid = threshold_minimum(image, bias='mid') +binary_mid = image > thresh_mid +thresh_max = threshold_minimum(image, bias='max') +binary_max = image > thresh_max + +fig, axes = plt.subplots(4, 2, figsize=(10, 10)) +ax = axes.ravel() + +ax[0].imshow(image, cmap=plt.cm.gray) +ax[0].set_title('Original') +ax[0].axis('off') + +ax[1].hist(image.ravel(), bins=256) +ax[1].set_title('Histogram') + +ax[2].imshow(binary_min, cmap=plt.cm.gray) +ax[2].set_title('Thresholded (min)') + +ax[3].hist(image.ravel(), bins=256) +ax[3].axvline(thresh_min, color='r') + +ax[4].imshow(binary_mid, cmap=plt.cm.gray) +ax[4].set_title('Thresholded (mid)') +ax[5].hist(image.ravel(), bins=256) +ax[5].axvline(thresh_mid, color='r') + +ax[6].imshow(binary_max, cmap=plt.cm.gray) +ax[6].set_title('Thresholded (max)') +ax[7].hist(image.ravel(), bins=256) +ax[7].axvline(thresh_max, color='r') + +for a in ax[::2]: + a.axis('off') +plt.show() + +###################################################################### +# Otsu's method [2]_ calculates an "optimal" threshold (marked by a red line in the +# histogram below) by maximizing the variance between two classes of pixels, +# which are separated by the threshold. Equivalently, this threshold minimizes +# the intra-class variance. +# +# .. [2] http://en.wikipedia.org/wiki/Otsu's_method +# + +from skimage.filters import threshold_otsu + + +image = data.camera() +thresh = threshold_otsu(image) +binary = image > thresh + +fig, axes = plt.subplots(ncols=3, figsize=(8, 2.5)) +ax = axes.ravel() +ax[0] = plt.subplot(1, 3, 1, adjustable='box-forced') +ax[1] = plt.subplot(1, 3, 2) +ax[2] = plt.subplot(1, 3, 3, sharex=ax[0], sharey=ax[0], adjustable='box-forced') + +ax[0].imshow(image, cmap=plt.cm.gray) +ax[0].set_title('Original') +ax[0].axis('off') + +ax[1].hist(image.ravel(), bins=256) +ax[1].set_title('Histogram') +ax[1].axvline(thresh, color='r') + +ax[2].imshow(binary, cmap=plt.cm.gray) +ax[2].set_title('Thresholded') +ax[2].axis('off') + +plt.show() + +###################################################################### +# Local thresholding +# ================== +# +# If the image background is relatively uniform, then you can use a global +# threshold value as presented above. However, if there is large variation in the +# background intensity, adaptive thresholding (a.k.a. local or dynamic +# thresholding) may produce better results. Note that local is much slower than +# global thresholding. +# +# Here, we binarize an image using the `threshold_adaptive` function, which +# calculates thresholds in regions with a characteristic size `block_size` surrounding +# each pixel (i.e. local neighborhoods). Each threshold value is the weighted mean +# of the local neighborhood minus an offset value. +# + +from skimage.filters import threshold_otsu, threshold_adaptive + + +image = data.page() + +global_thresh = threshold_otsu(image) +binary_global = image > global_thresh + +block_size = 35 +binary_adaptive = threshold_adaptive(image, block_size, offset=10) + +fig, axes = plt.subplots(nrows=3, figsize=(7, 8)) +ax = axes.ravel() +plt.gray() + +ax[0].imshow(image) +ax[0].set_title('Original') + +ax[1].imshow(binary_global) +ax[1].set_title('Global thresholding') + +ax[2].imshow(binary_adaptive) +ax[2].set_title('Adaptive thresholding') + +for a in ax: + a.axis('off') + +plt.show() + +###################################################################### +# Now, we show how Otsu's threshold [2]_ method can be applied locally. For +# each pixel, an "optimal" threshold is determined by maximizing the variance +# between two classes of pixels of the local neighborhood defined by a +# structuring element. +# +# The example compares the local threshold with the global threshold. +# + +from skimage.morphology import disk +from skimage.filters import threshold_otsu, rank +from skimage.util import img_as_ubyte + + +img = img_as_ubyte(data.page()) + +radius = 15 +selem = disk(radius) + +local_otsu = rank.otsu(img, selem) +threshold_global_otsu = threshold_otsu(img) +global_otsu = img >= threshold_global_otsu + +fig, axes = plt.subplots(2, 2, figsize=(8, 5), sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) +ax = axes.ravel() +plt.tight_layout() + +fig.colorbar(ax[0].imshow(img, cmap=plt.cm.gray), + ax=ax[0], orientation='horizontal') +ax[0].set_title('Original') +ax[0].axis('off') + +fig.colorbar(ax[1].imshow(local_otsu, cmap=plt.cm.gray), + ax=ax[1], orientation='horizontal') +ax[1].set_title('Local Otsu (radius=%d)' % radius) +ax[1].axis('off') + +ax[2].imshow(img >= local_otsu, cmap=plt.cm.gray) +ax[2].set_title('Original >= Local Otsu' % threshold_global_otsu) +ax[2].axis('off') + +ax[3].imshow(global_otsu, cmap=plt.cm.gray) +ax[3].set_title('Global Otsu (threshold = %d)' % threshold_global_otsu) +ax[3].axis('off') + +plt.show() diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index 9f722318..61fc850c 100644 --- a/skimage/filters/__init__.py +++ b/skimage/filters/__init__.py @@ -8,7 +8,8 @@ from .edges import (sobel, sobel_h, sobel_v, from ._rank_order import rank_order from ._gabor import gabor_kernel, gabor from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen, - threshold_isodata, threshold_li, threshold_minimum) + threshold_isodata, threshold_li, threshold_minimum, + threshold_mean, threshold_triangle, try_all_threshold) from . import rank from .rank import median @@ -44,6 +45,7 @@ __all__ = ['inverse', 'rank_order', 'gabor_kernel', 'gabor', + 'try_all_threshold', 'threshold_adaptive', 'threshold_otsu', 'threshold_yen', diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 2896c982..06493880 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -1,10 +1,15 @@ +import math import numpy as np from scipy import ndimage as ndi from scipy.ndimage import filters as ndif +from collections import OrderedDict from ..exposure import histogram from .._shared.utils import assert_nD, warn +from ..morphology import disk +from ..filters.rank import otsu -__all__ = ['threshold_adaptive', +__all__ = ['try_all_threshold', + 'threshold_adaptive', 'threshold_otsu', 'threshold_yen', 'threshold_isodata', @@ -14,6 +19,144 @@ __all__ = ['threshold_adaptive', 'threshold_triangle'] +def _try_all(image, methods=None, figsize=None, num_cols=2, verbose=True): + """Returns a figure comparing the outputs of different methods. + + Parameters + ---------- + image : (N, M) ndarray + Input image. + methods : dict, optional + Names and associated functions. + Functions must take and return an image. + figsize : tuple, optional + Figure size (in inches). + num_cols : int, optional + Number of columns. + verbose : bool, optional + Print function name for each method. + + Returns + ------- + fig, ax : tuple + Matplotlib figure and axes. + """ + from matplotlib import pyplot as plt + + num_rows = math.ceil((len(methods) + 1.) / num_cols) + num_rows = int(num_rows) # Python 2.7 support + fig, ax = plt.subplots(num_rows, num_cols, figsize=figsize, + sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) + ax = ax.ravel() + + ax[0].imshow(image, cmap=plt.cm.gray) + ax[0].set_title('Original') + + i = 1 + for name, func in methods.items(): + ax[i].imshow(func(image), cmap=plt.cm.gray) + ax[i].set_title(name) + i += 1 + if verbose: + print(func.__orifunc__) + + for a in ax: + a.axis('off') + + fig.tight_layout() + return fig, ax + + +def try_all_threshold(image, radius=None, figsize=(8, 5), verbose=True): + """Returns a figure comparing the outputs of different thresholding methods. + + Parameters + ---------- + image : (N, M) ndarray + Input image. + radius : int, optional + Lengthscale used for local methods. + If None, local methods are ignored. + figsize : tuple, optional + Figure size (in inches). + verbose : bool, optional + Print function name for each method. + + Returns + ------- + fig, ax : tuple + Matplotlib figure and axes. + + Notes + ----- + The following algorithms are used: + + * isodata + * li + * mean + * minimum + * otsu + * triangle + * yen + * adaptive threshold (local) + * rank otsu (local) + + Examples + -------- + >>> from skimage.data import text + >>> fig, ax = try_all_threshold(text(), radius=20, + ... figsize=(10, 6), verbose=False) + """ + + def include_selem(func, *args, **kwargs): + """ + A wrapper function to embed a threshold range for local algorithms. + """ + def wrapper(im): + return func(im, *args, **kwargs) + try: + wrapper.__orifunc__ = func.__orifunc__ + except AttributeError: + wrapper.__orifunc__ = func.__module__ + '.' + func.__name__ + return wrapper + + def thresh(func): + """ + A wrapper function to return a thresholded image. + """ + def wrapper(im): + return im > func(im) + try: + wrapper.__orifunc__ = func.__orifunc__ + except AttributeError: + wrapper.__orifunc__ = func.__module__ + '.' + func.__name__ + return wrapper + + # Global algorithms. + methods = OrderedDict({'Isodata': thresh(threshold_isodata), + 'Li': thresh(threshold_li), + 'Mean': thresh(threshold_mean), + 'Minimum': thresh(threshold_minimum), + 'Otsu': thresh(threshold_otsu), + 'Triangle': thresh(threshold_triangle), + 'Yen': thresh(threshold_yen)}) + + # Local algorithms. + if radius is not None: + selem = disk(radius) + local_otsu = include_selem(otsu, selem) + methods['Local Otsu'] = thresh(local_otsu) + + block_size = 2 * int(radius) + 1 + adaptive_threshold = include_selem(threshold_adaptive, block_size, + offset=10) + methods['Adaptive threshold'] = adaptive_threshold + + return _try_all(image, figsize=figsize, + methods=methods, verbose=verbose) + + def threshold_adaptive(image, block_size, method='gaussian', offset=0, mode='reflect', param=None): """Applies an adaptive threshold to an array. @@ -142,8 +285,8 @@ def threshold_otsu(image, nbins=256): # Check if the image is multi-colored or not if image.min() == image.max(): - raise ValueError("threshold_otsu is expected to work with images " \ - "having more than one color. The input image seems " \ + raise ValueError("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)