segmentation.py

import numpy as np
from skimage import io, restoration
from skimage import img_as_float
from scipy import ndimage


def hysteresis_thresholding(im, high_threshold=0.5, low_threshold=0.2):
    """
    Hysteresis thresholding, using two thresholds.

    Regions of values larger than low threshold are kept only if they have
    some pixels with values larger than high threshold.

    Parameters
    ----------

    im: ndarray
        grayscale image (can be n-D)

    high_threshold: float
        Value of the high threshold

    low_threshold: float
        Value of the low threshold
    """
    if high_threshold < low_threshold:
        raise ValueError('high threshold must be higher than low threshold')
    im = img_as_float(im)
    high_mask = im >= high_threshold
    low_mask = im >= low_threshold
    #
    # Segment the low-mask, then only keep low-segments that have
    # some high_mask component in them
    #
    strel = np.ones((3,) * im.ndim, bool)
    labels, count = ndimage.label(low_mask, strel)
    if count == 0:
        return low_mask

    sums = (np.array(ndimage.sum(high_mask, labels,
                             np.arange(count, dtype=np.int32) + 1),
                             copy=False, ndmin=1))
    good_label = np.zeros((count + 1,), bool)
    good_label[1:] = sums > 0
    output_mask = good_label[labels]
    return output_mask


im = io.imread('blobs.jpg', as_grey=True)
# Tune amount of denoising here
tv = restoration.denoise_tv_bregman(im, 8)
# Tune high and low threshold here
mask = hysteresis_thresholding(tv, 0.55, 0.38)

import matplotlib.pyplot as plt
plt.imshow(im, cmap='gray')
plt.contour(mask, [0.5], linewidths=[2])
plt.axis('off')
plt.show()