# Source code for molyso.generic.otsu

# -*- coding: utf-8 -*-
"""
otsu.py contains an implementation of Otsu's thresholding method, taken verbatim from scikit-image.
"""
# Taken from scikit-image, to remove dependency on scikit-image (for merely using one function)
# This file
# Copyright (C) 2011, the scikit-image team
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#  1. Redistributions of source code must retain the above copyright
#     notice, this list of conditions and the following disclaimer.
#  2. Redistributions in binary form must reproduce the above copyright
#     notice, this list of conditions and the following disclaimer in
#     the documentation and/or other materials provided with the
#     distribution.
#  3. Neither the name of skimage nor the names of its contributors may be
#     used to endorse or promote products derived from this software without
#     specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AS IS'' AND ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

import warnings
import numpy as np

[docs]def histogram(image, nbins=256):
"""
Return histogram of image.

Unlike numpy.histogram, this function returns the centers of bins and
does not rebin integer arrays. For integer arrays, each integer value has
its own bin, which improves speed and intensity-resolution.

The histogram is computed on the flattened image: for color images, the
function should be used separately on each channel to obtain a histogram
for each color channel.

:parameter image: Input image.
:parameter nbins: Number of bins used to calculate histogram. This value is ignored for integer arrays.
:type image: numpy.ndarray
:type nbins: int, optional
:return: The values of the histogram, the values at the center of the bins.
:rtype: tuple(numpy.ndarray, numpy.ndarray)
"""
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.")

# For integer types, histogramming with bincount is more efficient.
if np.issubdtype(image.dtype, np.integer):
offset = 0
if np.min(image) < 0:
offset = np.min(image)
hist = np.bincount(image.ravel() - offset)
bin_centers = np.arange(len(hist)) + offset

# clip histogram to start with a non-zero bin
idx = np.nonzero(hist)[0][0]
return hist[idx:], bin_centers[idx:]
else:
hist, bin_edges = np.histogram(image.flat, nbins)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.
return hist, bin_centers

[docs]def threshold_otsu(image, nbins=256):
"""
Return threshold value based on Otsu's method.

.. [1] Wikipedia, http://en.wikipedia.org/wiki/Otsu's_Method

:param image: Input image
:param nbins: Number of bins used to calculate histogram. This value is ignored for integer arrays.
:type image: numpy.ndarray
:type nbins: int, optional
:return: Upper threshold value. All pixels intensities that less or equal of this value assumed as foreground.
:rtype: float
"""
hist, bin_centers = histogram(image, nbins)
hist = hist.astype(float)

# class probabilities for all possible thresholds
weight1 = np.cumsum(hist)
weight2 = np.cumsum(hist[::-1])[::-1]
# class means for all possible thresholds
mean1 = np.cumsum(hist * bin_centers) / weight1
mean2 = (np.cumsum((hist * bin_centers)[::-1]) / weight2[::-1])[::-1]

# Clip ends to align class 1 and class 2 variables:
# The last value of weight1/mean1 should pair with zero values in
# weight2/mean2, which do not exist.
variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2

idx = np.argmax(variance12)
threshold = bin_centers[:-1][idx]
return threshold