WIP NEW histogram backprojection

doc/examples/segmentation/plot_backprojection.py

@@ -0,0 +1,74 @@
+"""
+========================
+Histogram Backprojection
+========================
+
+Histogram Backprojection is a tool for image segmentation based on object's
+color or intensity distribution.
+
+The histogram of an object, which you want to detect or track, is used as a
+feature to find the object in other images. So it takes two images as input,
+one is your object image and the next is the image where you want to find the
+object. The function then computes the histograms of these two images and
+return a **backprojected image** [1]_. A backprojected image is a grayscale image
+where each pixel denotes the probability of that pixel being part of the
+object. By thresholding this backprojected image with a suitable value gives
+the objects' mask.
+
+In below example, the brown region in the image is needed to be segmented.
+So, first 200x200 block of image is selected as object image. Then
+backprojection is applied to it. The result is thresholded with Otsu's
+thresholding to get the mask. Then a bitwise_and operation with input image
+and mask gives the segmented object.
+
+
+References
+----------
+
+.. [1] Swain, Michael J., and Dana H. Ballard. "Indexing via color histograms."
+       Active Perception and Robot Vision. Springer Berlin Heidelberg,
+       1992. 261-273.  :DOI:`10.1109/ICCV.1990.139558`
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage.segmentation import histogram_backprojection
+from skimage import data
+from skimage.filters import threshold_otsu
+from skimage import img_as_ubyte
+
+
+image = data.immunohistochemistry()
+image = img_as_ubyte(image)
+template = image[50:200, :150]
+
+# apply histogram backprojection
+backprojection = histogram_backprojection(image, template)
+
+# threshold the image with otsu's thresholding
+thresh = threshold_otsu(backprojection)
+mask = backprojection >= thresh
+# Set zero values where the mask is False
+res = np.where(mask[..., None], image, 0)
+
+fig, ax = plt.subplots(nrows=2, ncols=2)
+ax = ax.ravel()
+
+ax[0].imshow(image)
+ax[0].axis('off')
+ax[0].set_title('Input image')
+
+ax[1].imshow(template)
+ax[1].axis('off')
+ax[1].set_title('Template')
+
+ax[2].imshow(backprojection, cmap=plt.cm.gray)
+ax[2].axis('off')
+ax[2].set_title('Backprojected image')
+
+ax[3].imshow(res)
+ax[3].axis('off')
+ax[3].set_title('Segmented image')
+
+plt.show()

skimage/segmentation/__init__.py

@@ -12,6 +12,7 @@
                           morphological_chan_vese, inverse_gaussian_gradient,
                           circle_level_set, checkerboard_level_set)
 from ..morphology import flood, flood_fill
+from .backproject import histogram_backprojection
 
 
 __all__ = ['random_walker',
@@ -30,6 +31,7 @@
            'morphological_geodesic_active_contour',
            'morphological_chan_vese',
            'inverse_gaussian_gradient',
+           'histogram_backprojection',
            'circle_level_set',
            'checkerboard_level_set',
            'flood',

skimage/segmentation/backproject.py

@@ -0,0 +1,127 @@
+import numpy as np
+from scipy.ndimage import convolve
+
+from ..color import rgb2hsv
+from ..exposure import rescale_intensity, histogram
+from ..util import img_as_ubyte
+
+# kernel for final convolution
+disc = np.array([[0, 0, 1, 0, 0],
+                 [1, 1, 1, 1, 1],
+                 [1, 1, 1, 1, 1],
+                 [1, 1, 1, 1, 1],
+                 [0, 0, 1, 0, 0]], dtype=np.uint8)
+
+
+def _normalize_hsv(hsv):
+    """It changes the range of Hue and Saturation to [0, 179] and [0, 255]
+    respectively.
+
+    Parameters
+    ----------
+    hsv : HSV image array
+        HSV image array obtained as result of color.rgb2hsv()
+
+    Returns
+    -------
+    out : HSV image array
+        It is an uint8 image array in HSV format
+
+    """
+
+    return ([179, 255, 1] * hsv).astype(np.uint8)
+
+
+def histogram_backprojection(image, template, multichannel=True):
+    """Project the histogram of one image onto another
+
+    Parameters
+    ----------
+    image : (M, N, 3) ndarray or (M, N) ndarray
+        input image.
+    template : (M, N, 3) ndarray or (M, N) ndarray
+        Template for the histogram reference.
+    multichannel : bool, optional
+        Whether the last axis of the image is to be interpreted as multiple
+        channels or another spatial dimension.
+
+    Returns
+    -------
+    out : (M, N) ndarray
+        Backprojected image. Greyscale image that indicates the
+        level of correspondence to the histogram of the object image.
+
+    References
+    ----------
+    .. [1] Swain, Michael J., and Dana H. Ballard. "Indexing via color
+           histograms." Active Perception and Robot Vision. Springer Berlin
+           Heidelberg, 1992. 261-273.  DOI:`10.1109/ICCV.1990.139558`
+
+    """
+
+    # Both image should be of uint8 dtype
+    if image.dtype != np.uint8:
+        image = img_as_ubyte(image)
+    if template.dtype != np.uint8:
+        template = img_as_ubyte(template)
+
+    isGray = False
+    isColor = False
+
+    if image.ndim == 2 and template.ndim == 2:
+        isGray = True
+    elif multichannel and image.ndim == 3 and template.ndim == 3:
+        isColor = True
+    elif not multichannel and image.ndim >= 2 and template.ndim >= 2:
+        image = image.reshape(image.shape[:2])
+        template = template.reshape(template.shape[:2])
+        isGray = True
+    else:
+        raise ValueError("Both images should be 1-channel or 3-channel \
+            and multichannel should be True for 3-channel images")
+
+    # for grayscale image take 1D histogram of intensity values
+    if isGray:
+
+        # find histograms
+        hist1 = histogram(image, nbins=256, source_range='dtype')
+        hist2 = histogram(template, nbins=256, source_range='dtype')
+
+        # find their ratio hist2/hist1
+        R = np.float64(hist2) / (hist1 + 1)
+
+        # Now apply this ratio as the palette to original image, image
+        B = R[image]
+        B = np.minimum(B, 1)
+        B = rescale_intensity(B, out_range=(0, 255))
+        B = np.uint8(B)
+        B = convolve(B, disc)
+        return B
+
+    # if color image, take 2D histogram
+    elif isColor:
+        # convert images to hsv plane
+        hsv_image = rgb2hsv(image)
+        hsv_template = rgb2hsv(template)
+        hsv_image = _normalize_hsv(hsv_image)
+        hsv_template = _normalize_hsv(hsv_template)
+
+        # find their color 2D histograms
+        h1, s1, v1 = np.dsplit(hsv_image, (1, 2))
+        hist1, _, _ = np.histogram2d(h1.ravel(), s1.ravel(),
+                                     bins=[180, 256], range=[[0, 180], [0, 256]])
+        h2, s2, v2 = np.dsplit(hsv_template, (1, 2))
+        hist2, _, _ = np.histogram2d(h2.ravel(), s2.ravel(),
+                                     bins=[180, 256], range=[[0, 180], [0, 256]])
+
+        # find their ratio hist2/hist1
+        R = hist2 / (hist1 + 1)
+
+        # backproject
+        B = R[h1.ravel(), s1.ravel()]
+        B = np.minimum(B, 1)
+        B = B.reshape(image.shape[:2])
+        B = rescale_intensity(B, out_range=(0, 255))
+        B = convolve(B, disc)
+        B = np.clip(B, 0, 255).astype('uint8')
+        return B