From c87b1ad90ec1c24afdb85b019535c29218a8132e Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 19:48:58 -0500 Subject: [PATCH 1/8] Add exposure module with histogram equalization function --- skimage/exposure/__init__.py | 1 + skimage/exposure/exposure.py | 65 ++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+) create mode 100644 skimage/exposure/__init__.py create mode 100644 skimage/exposure/exposure.py diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py new file mode 100644 index 00000000..51dcc5cf --- /dev/null +++ b/skimage/exposure/__init__.py @@ -0,0 +1 @@ +from exposure import equalize_hist, cumulative_distribution diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py new file mode 100644 index 00000000..1d328e3e --- /dev/null +++ b/skimage/exposure/exposure.py @@ -0,0 +1,65 @@ +import numpy as np + + +__all__ = ['cumulative_distribution', 'equalize_hist'] + + +def cumulative_distribution(img, nbins=256): + """Return cumulative distribution function (cdf) for the given image. + + Parameters + ---------- + img : array + Image array. + nbins : int + Number of bins for image histogram. + + Returns + ------- + img_cdf : array + Values of cumulative distribution function. + bin_edges : array + Bin edges for cdf. Length is ``len(img_cdf) + 1``. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Cumulative_distribution_function + + """ + hist, bin_edges = np.histogram(img.flat, nbins, density=True) + img_cdf = hist.cumsum() + return img_cdf, bin_edges + + +def equalize_hist(img, nbins=256, max_intensity=255): + """Return image after histogram equalization. + + Parameters + ---------- + img : array + Image array. + nbins : int + Number of bins for image histogram. + max_intensity : int + Maximum intensity of the returned image. + + Returns + ------- + out : array + Image array after histogram equalization. + + Notes + ----- + This function is adapted from [1] with the author's permission. + + References + ---------- + .. [1] http://www.janeriksolem.net/2009/06/histogram-equalization-with-python-and.html + .. [2] http://en.wikipedia.org/wiki/Histogram_equalization + + """ + cdf, bin_edges = cumulative_distribution(img, nbins) + cdf = max_intensity * cdf / cdf[-1] + out = np.interp(img.flat, bin_edges[:-1], cdf) + return out.reshape(img.shape) + From 87c2353845923241592ccd26c0fd553f5996a42c Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Tue, 13 Dec 2011 22:57:15 -0500 Subject: [PATCH 2/8] Change API for equalize_hist and cumulative_distribution. `equalize_hist`: "max_intensity" parameter no longer exists---img_as_float normalizes intensity range `cumulative_distribution`: Return centers of bins instead of the edges. Move `histogram` function from filter subpackage. Add test of equalize_hist Add example of histogram equalization. --- doc/examples/plot_equalize_hist.py | 55 +++++++++++++++++ skimage/exposure/__init__.py | 2 +- skimage/exposure/exposure.py | 81 +++++++++++++++++++------ skimage/exposure/tests/test_exposure.py | 37 +++++++++++ skimage/filter/thresholding.py | 39 +----------- 5 files changed, 159 insertions(+), 55 deletions(-) create mode 100644 doc/examples/plot_equalize_hist.py create mode 100644 skimage/exposure/tests/test_exposure.py diff --git a/doc/examples/plot_equalize_hist.py b/doc/examples/plot_equalize_hist.py new file mode 100644 index 00000000..59b8defc --- /dev/null +++ b/doc/examples/plot_equalize_hist.py @@ -0,0 +1,55 @@ +""" +====================== +Histogram Equalization +====================== + +This examples takes an image with low contrast and enhances its contrast using +histogram equalization. Histogram equalization enhances contrast by "spreading +out the most frequent intensity values" in an image [1]. The equalized image +has a roughly linear cumulative distribution function, as shown in this example. + +.. [1] http://en.wikipedia.org/wiki/Histogram_equalization + +""" +import matplotlib.pyplot as plt + +from skimage import data +from skimage.util.dtype import dtype_range +from skimage import exposure + + +def plot_hist(img, bins=256, ax=None): + """Plot histogram and cumulative histogram for image""" + ax = ax if ax is not None else plt.gca() + img_cdf, bins = exposure.cumulative_distribution(img, bins) + ax.hist(img.ravel(), bins=bins) + ax_right = ax.twinx() + ax_right.plot(bins, img_cdf, 'r') + xmin, xmax = dtype_range[img.dtype.type] + ax.set_xlim(xmin, xmax) + + ax.set_ylabel('# pixels') + ax.set_xlabel('pixel intensiy') + ax_right.set_ylabel('fraction of total intensity') + + +img_orig = data.camera() +# squeeze image intensities to lower image contrast +img = img_orig / 5 + 100 +img_eq = exposure.equalize_hist(img) + +plt.subplot(2, 2, 1) +plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.axis('off') +plt.subplot(2, 2, 2) +plot_hist(img) + +plt.subplot(2, 2, 3) +plt.imshow(img_eq, cmap=plt.cm.gray, vmin=0, vmax=1) +plt.axis('off') +plt.subplot(2, 2, 4) +plot_hist(img_eq) + +plt.subplots_adjust(left=0.05, hspace=0.25, wspace=0.3, top=0.95, bottom=0.1) +plt.show() + diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py index 51dcc5cf..1c64aad4 100644 --- a/skimage/exposure/__init__.py +++ b/skimage/exposure/__init__.py @@ -1 +1 @@ -from exposure import equalize_hist, cumulative_distribution +from exposure import histogram, equalize_hist, cumulative_distribution diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 1d328e3e..06172597 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -1,15 +1,63 @@ import numpy as np - -__all__ = ['cumulative_distribution', 'equalize_hist'] +import skimage -def cumulative_distribution(img, nbins=256): +__all__ = ['histogram', 'cumulative_distribution', 'equalize_hist'] + + +def histogram(image, nbins, density=True): + """Return histogram of image. + + Unlike `numpy.histogram`, this function returns the centers of bins and + does not rebin integer arrays. + + Parameters + ---------- + image : array + Input image. + nbins : int + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + density : {True | False} + If True, values represent the probability density function. + If False, values represent the number of pixels in bins. + See numpy.histogram for details. + + Returns + ------- + hist : array + The values of the histogram. + bin_centers : array + The values at the center of the bins. + """ + 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: + + if np.version.version >= '1.6': + hist, bin_edges = np.histogram(image.flat, nbins, density=density) + else: + hist, bin_edges = np.histogram(image.flat, nbins, normed=density) + + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2. + return hist, bin_centers + + +def cumulative_distribution(image, nbins=256): """Return cumulative distribution function (cdf) for the given image. Parameters ---------- - img : array + image : array Image array. nbins : int Number of bins for image histogram. @@ -18,34 +66,33 @@ def cumulative_distribution(img, nbins=256): ------- img_cdf : array Values of cumulative distribution function. - bin_edges : array - Bin edges for cdf. Length is ``len(img_cdf) + 1``. + bin_centers : array + Centers of bins. References ---------- .. [1] http://en.wikipedia.org/wiki/Cumulative_distribution_function """ - hist, bin_edges = np.histogram(img.flat, nbins, density=True) + hist, bin_centers = histogram(image, nbins) img_cdf = hist.cumsum() - return img_cdf, bin_edges + img_cdf = img_cdf / float(img_cdf[-1]) + return img_cdf, bin_centers -def equalize_hist(img, nbins=256, max_intensity=255): +def equalize_hist(image, nbins=256): """Return image after histogram equalization. Parameters ---------- - img : array + image : array Image array. nbins : int Number of bins for image histogram. - max_intensity : int - Maximum intensity of the returned image. Returns ------- - out : array + out : float array Image array after histogram equalization. Notes @@ -58,8 +105,8 @@ def equalize_hist(img, nbins=256, max_intensity=255): .. [2] http://en.wikipedia.org/wiki/Histogram_equalization """ - cdf, bin_edges = cumulative_distribution(img, nbins) - cdf = max_intensity * cdf / cdf[-1] - out = np.interp(img.flat, bin_edges[:-1], cdf) - return out.reshape(img.shape) + image = skimage.img_as_float(image) + cdf, bin_centers = cumulative_distribution(image, nbins) + out = np.interp(image.flat, bin_centers, cdf) + return out.reshape(image.shape) diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py new file mode 100644 index 00000000..b2d19b8a --- /dev/null +++ b/skimage/exposure/tests/test_exposure.py @@ -0,0 +1,37 @@ +import numpy as np + +import skimage +from skimage import data +from skimage import exposure + + +# squeeze image intensities to lower image contrast +test_img = data.camera() / 5 + 100 + + +def test_equalize_hist_ubyte(): + img_eq = exposure.equalize_hist(test_img) + + cdf, bin_edges = exposure.cumulative_distribution(img_eq) + check_cdf_slope(cdf) + + +def test_equalize_hist_float(): + img = skimage.img_as_float(test_img) + img_eq = exposure.equalize_hist(img) + + cdf, bin_edges = exposure.cumulative_distribution(img_eq) + check_cdf_slope(cdf) + + +def check_cdf_slope(cdf): + """Slope of cdf which should equal 1 for an equalized histogram.""" + norm_intensity = np.linspace(0, 1, len(cdf)) + slope, intercept = np.polyfit(norm_intensity, cdf, 1) + assert 0.9 < slope < 1.1 + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() + diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 800e18b8..856bd78d 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,5 +1,7 @@ import numpy as np +from skimage.exposure import histogram + __all__ = ['threshold_otsu'] @@ -50,40 +52,3 @@ def threshold_otsu(image, nbins=256): threshold = bin_centers[:-1][idx] return threshold - -def histogram(image, nbins): - """Return histogram of image. - - Unlike `numpy.histogram`, this function returns the centers of bins and - does not rebin integer arrays. - - Parameters - ---------- - image : array - Input image. - nbins : int - Number of bins used to calculate histogram. This value is ignored for - integer arrays. - - Returns - ------- - hist : array - The values of the histogram. - bin_centers : array - The values at the center of the bins. - """ - 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 return only non-zero bins - idx = np.nonzero(hist)[0][0] - return hist[idx:], bin_centers[idx:] - else: - hist, bin_edges = np.histogram(image, bins=nbins) - bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2. - return hist, bin_centers - From d4ca519ca5d95a000b03c10e9ef11dabfd831388 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 22 Dec 2011 10:48:28 -0800 Subject: [PATCH 3/8] Rename equalize_hist to equalize and minor cleanup --- doc/examples/plot_equalize_hist.py | 21 ++++++++++----------- skimage/exposure/__init__.py | 2 +- skimage/exposure/exposure.py | 9 +++++---- skimage/exposure/tests/test_exposure.py | 8 ++++---- 4 files changed, 20 insertions(+), 20 deletions(-) diff --git a/doc/examples/plot_equalize_hist.py b/doc/examples/plot_equalize_hist.py index 59b8defc..e2ff895e 100644 --- a/doc/examples/plot_equalize_hist.py +++ b/doc/examples/plot_equalize_hist.py @@ -5,7 +5,7 @@ Histogram Equalization This examples takes an image with low contrast and enhances its contrast using histogram equalization. Histogram equalization enhances contrast by "spreading -out the most frequent intensity values" in an image [1]. The equalized image +out the most frequent intensity values" in an image [1]_. The equalized image has a roughly linear cumulative distribution function, as shown in this example. .. [1] http://en.wikipedia.org/wiki/Histogram_equalization @@ -18,25 +18,24 @@ from skimage.util.dtype import dtype_range from skimage import exposure -def plot_hist(img, bins=256, ax=None): +def plot_hist(img, bins=256): """Plot histogram and cumulative histogram for image""" - ax = ax if ax is not None else plt.gca() img_cdf, bins = exposure.cumulative_distribution(img, bins) - ax.hist(img.ravel(), bins=bins) - ax_right = ax.twinx() - ax_right.plot(bins, img_cdf, 'r') + plt.hist(img.ravel(), bins=bins) + ax_cdf = plt.twinx() + ax_cdf.plot(bins, img_cdf, 'r') xmin, xmax = dtype_range[img.dtype.type] - ax.set_xlim(xmin, xmax) + plt.xlim(xmin, xmax) - ax.set_ylabel('# pixels') - ax.set_xlabel('pixel intensiy') - ax_right.set_ylabel('fraction of total intensity') + plt.ylabel('# pixels') + plt.xlabel('pixel intensiy') + ax_cdf.set_ylabel('fraction of total intensity') img_orig = data.camera() # squeeze image intensities to lower image contrast img = img_orig / 5 + 100 -img_eq = exposure.equalize_hist(img) +img_eq = exposure.equalize(img) plt.subplot(2, 2, 1) plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255) diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py index 1c64aad4..ee5969ab 100644 --- a/skimage/exposure/__init__.py +++ b/skimage/exposure/__init__.py @@ -1 +1 @@ -from exposure import histogram, equalize_hist, cumulative_distribution +from exposure import histogram, equalize, cumulative_distribution diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 06172597..ca16c975 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -3,14 +3,15 @@ import numpy as np import skimage -__all__ = ['histogram', 'cumulative_distribution', 'equalize_hist'] +__all__ = ['histogram', 'cumulative_distribution', 'equalize'] def histogram(image, nbins, density=True): """Return histogram of image. Unlike `numpy.histogram`, this function returns the centers of bins and - does not rebin integer arrays. + does not rebin integer arrays. For integer arrays, each integer value has + its own bin, which improves speed and intensity-resolution. Parameters ---------- @@ -80,7 +81,7 @@ def cumulative_distribution(image, nbins=256): return img_cdf, bin_centers -def equalize_hist(image, nbins=256): +def equalize(image, nbins=256): """Return image after histogram equalization. Parameters @@ -97,7 +98,7 @@ def equalize_hist(image, nbins=256): Notes ----- - This function is adapted from [1] with the author's permission. + This function is adapted from [1]_ with the author's permission. References ---------- diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py index b2d19b8a..22339de7 100644 --- a/skimage/exposure/tests/test_exposure.py +++ b/skimage/exposure/tests/test_exposure.py @@ -9,16 +9,16 @@ from skimage import exposure test_img = data.camera() / 5 + 100 -def test_equalize_hist_ubyte(): - img_eq = exposure.equalize_hist(test_img) +def test_equalize_ubyte(): + img_eq = exposure.equalize(test_img) cdf, bin_edges = exposure.cumulative_distribution(img_eq) check_cdf_slope(cdf) -def test_equalize_hist_float(): +def test_equalize_float(): img = skimage.img_as_float(test_img) - img_eq = exposure.equalize_hist(img) + img_eq = exposure.equalize(img) cdf, bin_edges = exposure.cumulative_distribution(img_eq) check_cdf_slope(cdf) From 62d3b5590b0e5d45d1128bc60daa1a9eb39ee1eb Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 22 Dec 2011 10:58:02 -0800 Subject: [PATCH 4/8] Remove `density` parameter in `histogram` The `density` (or `normed`) parameter was set in the original implementation of `equalize`, but it is unnecessary since `cumulative_distribution` renormalizes values. The only other function in scikits-image that calls `histogram` (`filter.threshold_otsu`) is not affected by renormalization. --- skimage/exposure/exposure.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index ca16c975..d120b386 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -6,7 +6,7 @@ import skimage __all__ = ['histogram', 'cumulative_distribution', 'equalize'] -def histogram(image, nbins, density=True): +def histogram(image, nbins=256): """Return histogram of image. Unlike `numpy.histogram`, this function returns the centers of bins and @@ -20,10 +20,6 @@ def histogram(image, nbins, density=True): nbins : int Number of bins used to calculate histogram. This value is ignored for integer arrays. - density : {True | False} - If True, values represent the probability density function. - If False, values represent the number of pixels in bins. - See numpy.histogram for details. Returns ------- @@ -43,12 +39,7 @@ def histogram(image, nbins, density=True): idx = np.nonzero(hist)[0][0] return hist[idx:], bin_centers[idx:] else: - - if np.version.version >= '1.6': - hist, bin_edges = np.histogram(image.flat, nbins, density=density) - else: - hist, bin_edges = np.histogram(image.flat, nbins, normed=density) - + hist, bin_edges = np.histogram(image.flat, nbins) bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2. return hist, bin_centers From 60f146c780a65df0ba1429f09be8cba7c22c27d9 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Dec 2011 01:53:37 -0800 Subject: [PATCH 5/8] DOC: Update contributors. --- CONTRIBUTORS.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 007f267a..bb17033a 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -26,9 +26,11 @@ - Tony Yu Reading of paletted images; build, bug and doc fixes. Code to generate skimage logo. + Otsu thresholding, histogram equalisation, and more. - Zachary Pincus - Tracing of low cost paths, FreeImage I/O plugin, iso-contours + Tracing of low cost paths, FreeImage I/O plugin, iso-contours, + and more. - Almar Klein Binary heap class for graph algorithms @@ -45,7 +47,7 @@ - Emmanuelle Guillart Total variation noise filtering, integration of CellProfiler's - mathematical morphology tools. + mathematical morphology tools, tutorials, and more. - Maƫl Primet Total variation noise filtering From d74295c7f1cf98af31abaa2fc4033b5269923ae1 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Dec 2011 01:53:59 -0800 Subject: [PATCH 6/8] DOC: Explain why bincount is used to histogram integer arrays. --- skimage/exposure/exposure.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index d120b386..42a60dbf 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -28,6 +28,8 @@ def histogram(image, nbins=256): bin_centers : array The values at the center of the bins. """ + + # For integer types, histogramming with bincount is more efficient. if np.issubdtype(image.dtype, np.integer): offset = 0 if np.min(image) < 0: From 19fa34f269c71c571f82363c286b376f7b106dfd Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Dec 2011 01:55:00 -0800 Subject: [PATCH 7/8] DOC: Rename equalize example to match function name. --- doc/examples/{plot_equalize_hist.py => plot_equalize.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename doc/examples/{plot_equalize_hist.py => plot_equalize.py} (100%) diff --git a/doc/examples/plot_equalize_hist.py b/doc/examples/plot_equalize.py similarity index 100% rename from doc/examples/plot_equalize_hist.py rename to doc/examples/plot_equalize.py From 4e054fce7c5c6e2f4aa7df847cd0541e64897eb9 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Dec 2011 02:01:08 -0800 Subject: [PATCH 8/8] DOC: Move around plot statements so that all axes are labeled. --- doc/examples/plot_equalize.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py index e2ff895e..9cb30dd7 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/plot_equalize.py @@ -22,15 +22,15 @@ def plot_hist(img, bins=256): """Plot histogram and cumulative histogram for image""" img_cdf, bins = exposure.cumulative_distribution(img, bins) plt.hist(img.ravel(), bins=bins) + plt.ylabel('Number of pixels') + plt.xlabel('Pixel intensiy') + ax_cdf = plt.twinx() ax_cdf.plot(bins, img_cdf, 'r') xmin, xmax = dtype_range[img.dtype.type] plt.xlim(xmin, xmax) - plt.ylabel('# pixels') - plt.xlabel('pixel intensiy') - ax_cdf.set_ylabel('fraction of total intensity') - + ax_cdf.set_ylabel('Fraction of total intensity') img_orig = data.camera() # squeeze image intensities to lower image contrast @@ -39,12 +39,14 @@ img_eq = exposure.equalize(img) plt.subplot(2, 2, 1) plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.title('Low contrast input image') plt.axis('off') plt.subplot(2, 2, 2) plot_hist(img) plt.subplot(2, 2, 3) plt.imshow(img_eq, cmap=plt.cm.gray, vmin=0, vmax=1) +plt.title('After\nhistogram equalization') plt.axis('off') plt.subplot(2, 2, 4) plot_hist(img_eq)