这篇主要是想调用threshold_multiotsu，但是该函数只有在scikit-image0.16.0版本以上才有，但是0.16.0版本以上仅支持python3.6及以上版本。但是我现在所有的环境都是python3.5，为了能使用该函数，从scikit-image中将这部分内容抠取出来。
新建文件test.py
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from skimage import data
from multi_level.multi_thresh import threshold_multiotsu
# Setting the font size for all plots.
matplotlib.rcParams['font.size'] = 9
# The input image.
image = data.camera()
# Applying multi-Otsu threshold for the default value, generating
# three classes.
thresholds = threshold_multiotsu(image)
# Using the threshold values, we generate the three regions.
regions = np.digitize(image, bins=thresholds)
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(10, 3.5))
# Plotting the original image.
ax[0].imshow(image, cmap='gray')
ax[0].set_title('Original')
ax[0].axis('off')
# Plotting the histogram and the two thresholds obtained from
# multi-Otsu.
ax[1].hist(image.ravel(), bins=255)
ax[1].set_title('Histogram')
for thresh in thresholds:
ax[1].axvline(thresh, color='r')
# Plotting the Multi Otsu result.
ax[2].imshow(regions, cmap='Accent')
ax[2].set_title('Multi-Otsu result')
ax[2].axis('off')
plt.subplots_adjust()
plt.show()新建文件multi_thresh.py
import numpy as np
import pyximport
pyximport.install()
from multi_level._multiotsu import (_get_multiotsu_thresh_indices_lut, _get_multiotsu_thresh_indices)
_integer_types = (np.byte, np.ubyte, # 8 bits
np.short, np.ushort, # 16 bits
np.intc, np.uintc, # 16 or 32 or 64 bits
np.int_, np.uint, # 32 or 64 bits
np.longlong, np.ulonglong) # 64 bits
_integer_ranges = {t: (np.iinfo(t).min, np.iinfo(t).max)
for t in _integer_types}
dtype_range = {np.bool_: (False, True),
np.bool8: (False, True),
np.float16: (-1, 1),
np.float32: (-1, 1),
np.float64: (-1, 1)}
dtype_range.update(_integer_ranges)
_supported_types = list(dtype_range.keys())
def dtype_limits(image, clip_negative=False):
"""Return intensity limits, i.e. (min, max) tuple, of the image's dtype.
Parameters
----------
image : ndarray
Input image.
clip_negative : bool, optional
If True, clip the negative range (i.e. return 0 for min intensity)
even if the image dtype allows negative values.
Returns
-------
imin, imax : tuple
Lower and upper intensity limits.
"""
imin, imax = dtype_range[image.dtype.type]
if clip_negative:
imin = 0
return imin, imax
def _offset_array(arr, low_boundary, high_boundary):
"""Offset the array to get the lowest value at 0 if negative."""
if low_boundary < 0:
offset = low_boundary
dyn_range = high_boundary - low_boundary
# get smallest dtype that can hold both minimum and offset maximum
offset_dtype = np.promote_types(np.min_scalar_type(dyn_range),
np.min_scalar_type(low_boundary))
if arr.dtype != offset_dtype:
# prevent overflow errors when offsetting
arr = arr.astype(offset_dtype)
arr = arr - offset
else:
offset = 0
return arr, offset
def _bincount_histogram(image, source_range):
"""
Efficient histogram calculation for an image of integers.
This function is significantly more efficient than np.histogram but
works only on images of integers. It is based on np.bincount.
Parameters
----------
image : array
Input image.
source_range : string
'image' determines the range from the input image.
'dtype' determines the range from the expected range of the images
of that data type.
Returns
-------
hist : array
The values of the histogram.
bin_centers : array
The values at the center of the bins.
"""
if source_range not in ['image', 'dtype']:
raise ValueError('Incorrect value for `source_range` argument: {}'.format(source_range))
if source_range == 'image':
image_min = int(image.min().astype(np.int64))
image_max = int(image.max().astype(np.int64))
elif source_range == 'dtype':
image_min, image_max = dtype_limits(image, clip_negative=False)
image, offset = _offset_array(image, image_min, image_max)
hist = np.bincount(image.ravel(), minlength=image_max - image_min + 1)
bin_centers = np.arange(image_min, image_max + 1)
if source_range == 'image':
idx = max(image_min, 0)
hist = hist[idx:]
return hist, bin_centers
def histogram(image, nbins=256, source_range='image', normalize=False):
"""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.
Parameters
----------
image : array
Input image.
nbins : int, optional
Number of bins used to calculate histogram. This value is ignored for
integer arrays.
source_range : string, optional
'image' (default) determines the range from the input image.
'dtype' determines the range from the expected range of the images
of that data type.
normalize : bool, optional
If True, normalize the histogram by the sum of its values.
Returns
-------
hist : array
The values of the histogram.
bin_centers : array
The values at the center of the bins.
See Also
--------
cumulative_distribution
Examples
--------
>>> from skimage import data, exposure, img_as_float
>>> image = img_as_float(data.camera())
>>> np.histogram(image, bins=2)
(array([107432, 154712]), array([ 0. , 0.5, 1. ]))
>>> exposure.histogram(image, nbins=2)
(array([107432, 154712]), array([ 0.25, 0.75]))
"""
sh = image.shape
if len(sh) == 3 and sh[-1] < 4:
raise("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.")
image = image.flatten()
# For integer types, histogramming with bincount is more efficient.
if np.issubdtype(image.dtype, np.integer):
hist, bin_centers = _bincount_histogram(image, source_range)
else:
if source_range == 'image':
hist_range = None
elif source_range == 'dtype':
hist_range = dtype_limits(image, clip_negative=False)
else:
ValueError('Wrong value for the `source_range` argument')
hist, bin_edges = np.histogram(image, bins=nbins, range=hist_range)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.
if normalize:
hist = hist / np.sum(hist)
return hist, bin_centers
def threshold_multiotsu(image, classes=3, nbins=256):
r"""Generate `classes`-1 threshold values to divide gray levels in `image`.
The threshold values are chosen to maximize the total sum of pairwise
variances between the thresholded graylevel classes. See Notes and [1]_
for more details.
Parameters
----------
image : (N, M) ndarray
Grayscale input image.
classes : int, optional
Number of classes to be thresholded, i.e. the number of resulting
regions.
nbins : int, optional
Number of bins used to calculate the histogram. This value is ignored
for integer arrays.
Returns
-------
thresh : array
Array containing the threshold values for the desired classes.
Raises
------
ValueError
If ``image`` contains less grayscale value then the desired
number of classes.
Notes
-----
This implementation relies on a Cython function whose complexity
is :math:`O\left(\frac{Ch^{C-1}}{(C-1)!}\right)`, where :math:`h`
is the number of histogram bins and :math:`C` is the number of
classes desired.
The input image must be grayscale.
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
.. [2] Tosa, Y., "Multi-Otsu Threshold", a java plugin for ImageJ.
Available at:
<http://imagej.net/plugins/download/Multi_OtsuThreshold.java>
Examples
--------
>>> from skimage.color import label2rgb
>>> from skimage import data
>>> image = data.camera()
>>> thresholds = threshold_multiotsu(image)
>>> regions = np.digitize(image, bins=thresholds)
>>> regions_colorized = label2rgb(regions)
"""
if len(image.shape) > 2 and image.shape[-1] in (3, 4):
msg = ("threshold_multiotsu is expected to work correctly only for "
"grayscale images; image shape {0} looks like an RGB image")
raise (msg.format(image.shape))
# calculating the histogram and the probability of each gray level.
prob, bin_centers = histogram(image.ravel(),
nbins=nbins,
source_range='image',
normalize=True)
prob = prob.astype('float32')
nvalues = np.count_nonzero(prob)
if nvalues < classes:
msg = ("The input image has only {} different values. "
"It can not be thresholded in {} classes")
raise ValueError(msg.format(nvalues, classes))
elif nvalues == classes:
thresh_idx = np.where(prob > 0)[0][:-1]
else:
# Get threshold indices
try:
thresh_idx = _get_multiotsu_thresh_indices_lut(prob, classes - 1)
except MemoryError:
# Don't use LUT if the number of bins is too large (if the
# image is uint16 for example): in this case, the
# allocated memory is too large.
thresh_idx = _get_multiotsu_thresh_indices(prob, classes - 1)
thresh = bin_centers[thresh_idx]
return thresh新建文件_multiotsu.pyx
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
cimport cython
def _get_multiotsu_thresh_indices_lut(float [::1] prob,
Py_ssize_t thresh_count):
"""Finds the indices of Otsu thresholds according to the values
occurence probabilities.
This implementation uses a LUT to reduce the number of floating
point operations (see [1]_). The use of the LUT reduces the
computation time at the price of more memory consumption.
Parameters
----------
prob : array
Value occurence probabilities.
thresh_count : int
The desired number of thresholds (classes-1).
Returns
-------
py_thresh_indices : ndarray
The indices of the desired thresholds.
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
"""
cdef Py_ssize_t nbins = prob.shape[0]
py_thresh_indices = np.empty(thresh_count, dtype=np.intp)
cdef Py_ssize_t[::1] thresh_indices = py_thresh_indices
cdef Py_ssize_t[::1] current_indices = np.empty(thresh_count, dtype=np.intp)
cdef float [::1] var_btwcls = np.zeros((nbins * (nbins + 1)) / 2,
dtype=np.float32)
cdef float [::1] zeroth_moment = np.empty(nbins, dtype=np.float32)
cdef float [::1] first_moment = np.empty(nbins, dtype=np.float32)
with nogil:
_set_var_btwcls_lut(prob, nbins, var_btwcls, zeroth_moment,
first_moment)
_set_thresh_indices_lut(var_btwcls, hist_idx=0,
thresh_idx=0, nbins=nbins,
thresh_count=thresh_count, sigma_max=0,
current_indices=current_indices,
thresh_indices=thresh_indices)
return py_thresh_indices
cdef void _set_var_btwcls_lut(float [::1] prob,
Py_ssize_t nbins,
float [::1] var_btwcls,
float [::1] zeroth_moment,
float [::1] first_moment) nogil:
"""Builds the lookup table containing the variance between classes.
The variance between classes are stored in
``var_btwcls``. ``zeroth_moment`` and ``first_moment`` are buffers
for storing the first row of respectively the zeroth and first
order moments lookup table (respectively H, P and S in [1]_).
Parameters
----------
prob : array
Value occurence probabilities.
nbins : int
The number of intensity values.
var_btwcls : array
The upper triangular part of the lookup table containing the
variance between classes (referred to as H in [1]_). Its size
is equal to nbins*(nbins + 1) / 2.
zeroth_moment : array
First row of the zeroth order moments LUT (referred to as P in
[1]_).
first_moment : array
First row of the first order moments LUT (referred to as S in
[1]_).
Notes
-----
Only the first rows of the moments lookup tables are necessary to
build the lookup table containing the variance between
classes. ``var_btwcls`` is stored in the compressed upper
triangular matrix form (i.e. the seros of the lower triangular
part are not stored).
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
"""
cdef cnp.intp_t i, j, idx
cdef float zeroth_moment_ij, first_moment_ij
zeroth_moment[0] = prob[0]
first_moment[0] = prob[0]
for i in range(1, nbins):
zeroth_moment[i] = zeroth_moment[i - 1] + prob[i]
first_moment[i] = first_moment[i - 1] + i * prob[i]
if zeroth_moment[i] > 0:
var_btwcls[i] = (first_moment[i]**2) / zeroth_moment[i]
idx = nbins
for i in range(1, nbins):
for j in range(i, nbins):
zeroth_moment_ij = zeroth_moment[j] - zeroth_moment[i - 1]
if zeroth_moment_ij > 0:
first_moment_ij = first_moment[j] - first_moment[i - 1]
var_btwcls[idx] = (first_moment_ij**2) / zeroth_moment_ij
idx += 1
cdef float _get_var_btwclas_lut(float [::1] var_btwcls, Py_ssize_t i,
Py_ssize_t j, Py_ssize_t nbins) nogil:
"""Returns the variance between classes stored in compressed upper
triangular matrix form at the desired 2D indices.
Parameters
----------
var_btwcls : array
The lookup table containing the variance between classes in
compressed upper triangular matrix form.
i, j : int
2D indices in the uncompressed lookup table.
nbins : int
The number of columns in the lookup table.
Returns
-------
value : float
The value of the lookup table corresponding to index (i, j).
"""
cdef cnp.intp_t idx = (i * (2 * nbins - i + 1)) / 2 + j - i
return var_btwcls[idx]
cdef float _set_thresh_indices_lut(float[::1] var_btwcls, Py_ssize_t hist_idx,
Py_ssize_t thresh_idx, Py_ssize_t nbins,
Py_ssize_t thresh_count, float sigma_max,
Py_ssize_t[::1] current_indices,
Py_ssize_t[::1] thresh_indices) nogil:
"""Recursive function for finding the indices of the thresholds
maximizing the variance between classes sigma.
This implementation use a lookup table of variance between classes
to perform a brute force evaluation of sigma over all the
combinations of threshold to find the indices maximizing sigma
(see [1]_).
Parameters
----------
var_btwcls : array
The upper triangular part of the lookup table containing the
variance between classes (referred to as H in [1]_). Its size
is equal to nbins*(nbins + 1) / 2.
hist_idx : int
Current index in the histogram.
thresh_idx : int
Current index in thresh_indices.
nbins : int
number of bins used in the histogram
thresh_count : int
The desired number of thresholds (classes-1).
sigma_max : float
Current maximum variance between classes.
current_indices : array
Current evalueted threshold indices.
thresh_indices : array
The indices of thresholds maximizing the variance between
classes.
Returns
-------
max_sigma : float
Maximum variance between classes.
Notes
-----
For any candidate current_indices {t_0, ..., t_n}, sigma equals
var_btwcls[0, t_0] + var_btwcls[t_0+1, t_1] + ...
+ var_btwcls[t_(i-1) + 1, t_i] + ... + var_btwcls[t_n, nbins-1]
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
"""
cdef cnp.intp_t idx
cdef float sigma
if thresh_idx < thresh_count:
for idx in range(hist_idx, nbins - thresh_count + thresh_idx):
current_indices[thresh_idx] = idx
sigma_max = _set_thresh_indices_lut(var_btwcls, hist_idx=idx + 1,
thresh_idx=thresh_idx + 1,
nbins=nbins,
thresh_count=thresh_count,
sigma_max=sigma_max,
current_indices=current_indices,
thresh_indices=thresh_indices)
else:
sigma = (_get_var_btwclas_lut(var_btwcls, 0, current_indices[0], nbins)
+ _get_var_btwclas_lut(var_btwcls,
current_indices[thresh_count - 1] + 1,
nbins - 1, nbins))
for idx in range(thresh_count - 1):
sigma += _get_var_btwclas_lut(var_btwcls, current_indices[idx] + 1,
current_indices[idx + 1], nbins)
if sigma > sigma_max:
sigma_max = sigma
thresh_indices[:] = current_indices[:]
return sigma_max
def _get_multiotsu_thresh_indices(float [::1] prob, Py_ssize_t thresh_count):
"""Finds the indices of Otsu thresholds according to the values
occurence probabilities.
This implementation, as opposed to `_get_multiotsu_thresh_indices_lut`,
does not use LUT. It is therefore slower.
Parameters
----------
prob : array
Value occurence probabilities.
thresh_count : int
The desired number of threshold.
Returns
-------
py_thresh_indices : array
The indices of the desired thresholds.
"""
cdef Py_ssize_t nbins = prob.shape[0]
py_thresh_indices = np.empty(thresh_count, dtype=np.intp)
cdef Py_ssize_t[::1] thresh_indices = py_thresh_indices
cdef Py_ssize_t[::1] current_indices = np.empty(thresh_count, dtype=np.intp)
cdef float [::1] zeroth_moment = np.empty(nbins, dtype=np.float32)
cdef float [::1] first_moment = np.empty(nbins, dtype=np.float32)
with nogil:
_set_moments_lut_first_row(prob, nbins, zeroth_moment, first_moment)
_set_thresh_indices(zeroth_moment, first_moment, hist_idx=0,
thresh_idx=0, nbins=nbins,
thresh_count=thresh_count, sigma_max=0,
current_indices=current_indices,
thresh_indices=thresh_indices)
return py_thresh_indices
cdef void _set_moments_lut_first_row(float [::1] prob,
Py_ssize_t nbins,
float [::1] zeroth_moment,
float [::1] first_moment) nogil:
"""Builds the first rows of the zeroth and first moments lookup table
necessary to the computation of the variance between class.
Parameters
----------
prob : array
Value occurence probabilities.
nbins : int
The number of intensity values.
zeroth_moment : array
First row of the zeroth order moments LUT (referred to as P in
[1]_).
first_moment : array
First row of the first order moments LUT (referred to as S in
[1]_).
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
"""
cdef cnp.intp_t i
zeroth_moment[0] = prob[0]
first_moment[0] = prob[0]
for i in range(1, nbins):
zeroth_moment[i] = zeroth_moment[i - 1] + prob[i]
first_moment[i] = first_moment[i - 1] + i * prob[i]
cdef float _get_var_btwclas(float [::1] zeroth_moment,
float [::1] first_moment,
Py_ssize_t i, Py_ssize_t j) nogil:
"""Computes the variance between two classes.
Parameters
----------
zeroth_moment : array
First row of the zeroth order moments LUT (referred to as P in
[1]_).
first_moment : array
First row of the first order moments LUT (referred to as S in
[1]_).
i, j : int
The indices of the two considred classes.
Returns
-------
value : float
The variance between the classes i and j.
"""
cdef float zeroth_moment_ij, first_moment_ij
if i == 0:
if zeroth_moment[i] > 0:
return (first_moment[j]**2) / zeroth_moment[j]
else:
zeroth_moment_ij = zeroth_moment[j] - zeroth_moment[i - 1]
if zeroth_moment_ij > 0:
first_moment_ij = first_moment[j] - first_moment[i - 1]
return (first_moment_ij**2) / zeroth_moment_ij
return 0
cdef float _set_thresh_indices(float[::1] zeroth_moment,
float[::1] first_moment,
Py_ssize_t hist_idx,
Py_ssize_t thresh_idx, Py_ssize_t nbins,
Py_ssize_t thresh_count, float sigma_max,
Py_ssize_t[::1] current_indices,
Py_ssize_t[::1] thresh_indices) nogil:
"""Recursive function for finding the indices of the thresholds
maximizing the variance between classes sigma.
This implementation uses the first rows of the zeroth and first
moments lookup table to compute the variance between class and
performs a brute force evaluation of sigma over all the
combinations of threshold to find the indices maximizing sigma
(see [1]_)..
Parameters
----------
zeroth_moment : array
First row of the zeroth order moments LUT (referred to as P in
[1]_).
first_moment : array
First row of the first order moments LUT (referred to as S in
[1]_).
hist_idx : int
Current index in the histogram.
thresh_idx : int
Current index in thresh_indices.
nbins : int
number of bins used in the histogram
thresh_count : int
The desired number of thresholds (classes-1).
sigma_max : float
Current maximum variance between classes.
current_indices : array
Current evalueted threshold indices.
thresh_indices : array
The indices of thresholds maximizing the variance between
classes.
Returns
-------
max_sigma : float
Maximum variance between classes.
References
----------
.. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
multilevel thresholding", Journal of Information Science and
Engineering 17 (5): 713-727, 2001. Available at:
<http://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
:DOI:`10.6688/JISE.2001.17.5.1`
"""
cdef cnp.intp_t idx
cdef float sigma
if thresh_idx < thresh_count:
for idx in range(hist_idx, nbins - thresh_count + thresh_idx):
current_indices[thresh_idx] = idx
sigma_max = _set_thresh_indices(zeroth_moment,
first_moment,
hist_idx=idx + 1,
thresh_idx=thresh_idx + 1,
nbins=nbins,
thresh_count=thresh_count,
sigma_max=sigma_max,
current_indices=current_indices,
thresh_indices=thresh_indices)
else:
sigma = (_get_var_btwclas(zeroth_moment, first_moment, 0,
current_indices[0])
+ _get_var_btwclas(zeroth_moment, first_moment,
current_indices[thresh_count - 1] + 1,
nbins - 1))
for idx in range(thresh_count - 1):
sigma += _get_var_btwclas(zeroth_moment, first_moment,
current_indices[idx] + 1,
current_indices[idx + 1])
if sigma > sigma_max:
sigma_max = sigma
thresh_indices[:] = current_indices[:]
return sigma_max
新建文件setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
from Cython.Distutils import build_ext
import numpy as np
# ext_module = cythonize("TestOMP.pyx")
ext_module = Extension(
"_multiotsu",
["_multiotsu.pyx"],
)
setup(
cmdclass={'build_ext': build_ext},
ext_modules=[ext_module],
include_dirs=[np.get_include()]
) 编译_multiotsu.pyx，首先cd到该文件所在目录，然后使用python setup.py build_ext --inplace，然后会出现以下文件夹
7869
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
