From 2970de33ee015c908dda02839f1e388476d6d02e Mon Sep 17 00:00:00 2001 From: Jeremy Metz Date: Tue, 17 Feb 2015 14:36:58 +0000 Subject: [PATCH 1/5] Initial draft of Li thresholding (comments left in) --- skimage/filters/__init__.py | 3 +- skimage/filters/thresholding.py | 97 ++++++++++++++++++++++++++++++++- 2 files changed, 98 insertions(+), 2 deletions(-) diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index 795819c6..7cffb140 100644 --- a/skimage/filters/__init__.py +++ b/skimage/filters/__init__.py @@ -9,7 +9,7 @@ from .edges import (sobel, hsobel, vsobel, sobel_h, sobel_v, from ._rank_order import rank_order from ._gabor import gabor_kernel, gabor_filter from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen, - threshold_isodata) + threshold_isodata, threshold_li) from . import rank from .rank import median @@ -66,4 +66,5 @@ __all__ = ['inverse', 'threshold_otsu', 'threshold_yen', 'threshold_isodata', + 'threshold_li', 'rank'] diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index e10bca5c..5bd1712d 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -1,7 +1,8 @@ __all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen', - 'threshold_isodata'] + 'threshold_isodata', + 'threshold_li',] import numpy as np import scipy.ndimage @@ -302,3 +303,97 @@ def threshold_isodata(image, nbins=256, return_all=False): return thresholds else: return thresholds[0] + + +def threshold_li(image, nbins=256): + """Return threshold value based on Li's Minimum Cross Entropy method. + + Parameters + ---------- + image : array + Input image. + nbins : int, optional + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + + Returns + ------- + threshold : float + Upper threshold value. All pixels intensities that less or equal of + this value assumed as foreground. + + References + ---------- + .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" + Pattern Recognition, 26(4): 617-625 + .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum + Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776 + .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding + Techniques and Quantitative Performance Evaluation" Journal of + Electronic Imaging, 13(1): 146-165 + http://citeseer.ist.psu.edu/sezgin04survey.html + + Ported to skimage by J. Metz from ImageJ plugin by G.Landini + + Examples + -------- + >>> from skimage.data import camera + >>> image = camera() + >>> thresh = threshold_li(image) + >>> binary = image <= thresh + """ + + tolerance=0.5 + num_pixels = image.size + + # Calculate the mean gray-level + mean = image.mean() + #for (int ih = 0 + 1; ih < 256; ih++ ) //0 + 1? + # mean += (double)ih * data[ih]; + #mean /= num_pixels; + # Initial estimate + new_thresh = mean + old_thresh = new_thresh + 2*tolerance + + # Stop the iterations when the difference between the + # new and old threshold values is less than the tolerance + while abs( new_thresh - old_thresh ) > tolerance: + old_thresh = new_thresh + threshold = int(old_thresh + 0.5) # range + # Calculate the means of background and object pixels + # Background + #sum_back = 0 + #num_back = 0; + #for (int ih = 0; ih <= threshold; ih++ ) { + # sum_back += (double)ih * data[ih]; + # num_back += data[ih]; + #} + #mean_back = ( num_back == 0 ? 0.0 : ( sum_back / ( double ) num_back ) ); + mean_back = image[ image <= threshold ].mean() + # Object + #sum_obj = 0; + #num_obj = 0; + #for (int ih = threshold + 1; ih < 256; ih++ ) { + # sum_obj += (double)ih * data[ih]; + # num_obj += data[ih]; + #} + #mean_obj = ( num_obj == 0 ? 0.0 : ( sum_obj / ( double ) num_obj ) ); + mean_obj = image[ image > threshold ].mean() + + # Calculate the new threshold: Equation (7) in Ref. 2 + # //new_thresh = simple_round ( ( mean_back - mean_obj ) / ( Math.log ( mean_back ) - Math.log ( mean_obj ) ) ); + # //simple_round ( double x ) { + # // return ( int ) ( IS_NEG ( x ) ? x - .5 : x + .5 ); + # //} + # // + # //#define IS_NEG( x ) ( ( x ) < -DBL_EPSILON ) + # //DBL_EPSILON = 2.220446049250313E-16 + temp = (mean_back - mean_obj) / (np.log(mean_back) - np.log(mean_obj)) + + if temp < 0: # (temp < -2.220446049250313E-16) + new_thresh = int(temp - 0.5) + else: + new_thresh = int(temp + 0.5) + + return threshold #; + From 979d7515f0a34f2a32bee475c79ff842dbd9be0e Mon Sep 17 00:00:00 2001 From: Jeremy Metz Date: Tue, 17 Feb 2015 14:45:31 +0000 Subject: [PATCH 2/5] Cleaned up comments, removed unnecessary vars --- skimage/filters/thresholding.py | 30 ++---------------------------- 1 file changed, 2 insertions(+), 28 deletions(-) diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 5bd1712d..77e1a00d 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -344,13 +344,9 @@ def threshold_li(image, nbins=256): """ tolerance=0.5 - num_pixels = image.size - # Calculate the mean gray-level mean = image.mean() - #for (int ih = 0 + 1; ih < 256; ih++ ) //0 + 1? - # mean += (double)ih * data[ih]; - #mean /= num_pixels; + # Initial estimate new_thresh = mean old_thresh = new_thresh + 2*tolerance @@ -362,32 +358,10 @@ def threshold_li(image, nbins=256): threshold = int(old_thresh + 0.5) # range # Calculate the means of background and object pixels # Background - #sum_back = 0 - #num_back = 0; - #for (int ih = 0; ih <= threshold; ih++ ) { - # sum_back += (double)ih * data[ih]; - # num_back += data[ih]; - #} - #mean_back = ( num_back == 0 ? 0.0 : ( sum_back / ( double ) num_back ) ); mean_back = image[ image <= threshold ].mean() # Object - #sum_obj = 0; - #num_obj = 0; - #for (int ih = threshold + 1; ih < 256; ih++ ) { - # sum_obj += (double)ih * data[ih]; - # num_obj += data[ih]; - #} - #mean_obj = ( num_obj == 0 ? 0.0 : ( sum_obj / ( double ) num_obj ) ); mean_obj = image[ image > threshold ].mean() - # Calculate the new threshold: Equation (7) in Ref. 2 - # //new_thresh = simple_round ( ( mean_back - mean_obj ) / ( Math.log ( mean_back ) - Math.log ( mean_obj ) ) ); - # //simple_round ( double x ) { - # // return ( int ) ( IS_NEG ( x ) ? x - .5 : x + .5 ); - # //} - # // - # //#define IS_NEG( x ) ( ( x ) < -DBL_EPSILON ) - # //DBL_EPSILON = 2.220446049250313E-16 temp = (mean_back - mean_obj) / (np.log(mean_back) - np.log(mean_obj)) if temp < 0: # (temp < -2.220446049250313E-16) @@ -395,5 +369,5 @@ def threshold_li(image, nbins=256): else: new_thresh = int(temp + 0.5) - return threshold #; + return threshold From 5eef38bba5c6021103c4bdea3a63ee22a2ffac29 Mon Sep 17 00:00:00 2001 From: Jeremy Metz Date: Tue, 17 Feb 2015 14:55:27 +0000 Subject: [PATCH 3/5] Removed nbins input (from use with histogram input) --- skimage/filters/thresholding.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 77e1a00d..6487f4c3 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -305,16 +305,13 @@ def threshold_isodata(image, nbins=256, return_all=False): return thresholds[0] -def threshold_li(image, nbins=256): +def threshold_li(image): """Return threshold value based on Li's Minimum Cross Entropy method. Parameters ---------- image : array Input image. - nbins : int, optional - Number of bins used to calculate histogram. This value is ignored for - integer arrays. Returns ------- From 3cbae843602532fb23b0ee9886a96dd122807694 Mon Sep 17 00:00:00 2001 From: Jeremy Metz Date: Mon, 23 Feb 2015 11:22:12 +0000 Subject: [PATCH 4/5] PEP8 fixes, License file added, tests added * Added 7 tests (3 in TestSimpleImage, 4 as separate functions) * Added ImageJ license file * Updated Li code * threshold set from image instead of fixed at 0.5 * subtract min to handle -ve images (as log is integral part of alg.) --- skimage/filters/LICENSE.txt | 28 ++++++++++ skimage/filters/tests/test_thresholding.py | 35 ++++++++++++ skimage/filters/thresholding.py | 62 ++++++++++++---------- 3 files changed, 97 insertions(+), 28 deletions(-) create mode 100644 skimage/filters/LICENSE.txt diff --git a/skimage/filters/LICENSE.txt b/skimage/filters/LICENSE.txt new file mode 100644 index 00000000..3507299a --- /dev/null +++ b/skimage/filters/LICENSE.txt @@ -0,0 +1,28 @@ +threshold_li based on imagej code: + +Copyright (c) 2009 - 2015, Board of Regents of the University of +Wisconsin-Madison, Broad Institute of MIT and Harvard, and Max Planck +Institute of Molecular Cell Biology and Genetics. +All rights reserved. + +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. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 COPYRIGHT HOLDER OR CONTRIBUTORS 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. diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index 5de758e8..c360a030 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/skimage/filters/tests/test_thresholding.py @@ -5,6 +5,7 @@ import skimage from skimage import data from skimage.filters.thresholding import (threshold_adaptive, threshold_otsu, + threshold_li, threshold_yen, threshold_isodata) @@ -28,6 +29,17 @@ class TestSimpleImage(): image = np.float64(self.image) assert 2 <= threshold_otsu(image) < 3 + def test_li(self): + assert int(threshold_li(self.image)) == 2 + + def test_li_negative_int(self): + image = self.image - 2 + assert int(threshold_li(image)) == 0 + + def test_li_float_image(self): + image = np.float64(self.image) + assert 2 <= threshold_li(image) < 3 + def test_yen(self): assert threshold_yen(self.image) == 2 @@ -152,6 +164,29 @@ def test_otsu_astro_image(): img = skimage.img_as_ubyte(data.astronaut()) assert 109 < threshold_otsu(img) < 111 +def test_li_camera_image(): + camera = skimage.img_as_ubyte(data.camera()) + assert 63 < threshold_li(camera) < 65 + + +def test_li_coins_image(): + coins = skimage.img_as_ubyte(data.coins()) + assert 95 < threshold_li(coins) < 97 + + +def test_li_coins_image_as_float(): + coins = skimage.img_as_float(data.coins()) + assert 0.37 < threshold_li(coins) < 0.38 + + +def test_li_lena_image(): + img = skimage.img_as_ubyte(data.lena()) + assert 127 < threshold_li(img) < 129 + +def test_li_astro_image(): + img = skimage.img_as_ubyte(data.astronaut()) + assert 66 < threshold_li(img) < 68 + def test_yen_camera_image(): camera = skimage.img_as_ubyte(data.camera()) assert 197 < threshold_yen(camera) < 199 diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 6487f4c3..d5907f8e 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -1,8 +1,8 @@ __all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen', - 'threshold_isodata', - 'threshold_li',] + 'threshold_isodata', + 'threshold_li', ] import numpy as np import scipy.ndimage @@ -306,7 +306,7 @@ def threshold_isodata(image, nbins=256, return_all=False): def threshold_li(image): - """Return threshold value based on Li's Minimum Cross Entropy method. + """Return threshold value based on adaptation of Li's Minimum Cross Entropy method. Parameters ---------- @@ -321,16 +321,17 @@ def threshold_li(image): References ---------- - .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" + .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" Pattern Recognition, 26(4): 617-625 - .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum + .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776 - .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding - Techniques and Quantitative Performance Evaluation" Journal of - Electronic Imaging, 13(1): 146-165 + .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding + Techniques and Quantitative Performance Evaluation" Journal of + Electronic Imaging, 13(1): 146-165 http://citeseer.ist.psu.edu/sezgin04survey.html - - Ported to skimage by J. Metz from ImageJ plugin by G.Landini + .. [4] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold + + Adapted for skimage by J. Metz from ImageJ plugin by G.Landini Examples -------- @@ -339,32 +340,37 @@ def threshold_li(image): >>> thresh = threshold_li(image) >>> binary = image <= thresh """ + # Requires positive image ( log(mean)) + offset = image.min() + # Can't use fixed tolerance for float image + imrange = image.max()-offset + image -= offset - tolerance=0.5 - # Calculate the mean gray-level + tolerance = 0.5 * imrange / 256.0 + # Calculate the mean gray-level mean = image.mean() - - # Initial estimate + + # Initial estimate new_thresh = mean - old_thresh = new_thresh + 2*tolerance + old_thresh = new_thresh + 2 * tolerance # Stop the iterations when the difference between the - # new and old threshold values is less than the tolerance - while abs( new_thresh - old_thresh ) > tolerance: + # new and old threshold values is less than the tolerance + while abs(new_thresh - old_thresh) > tolerance: old_thresh = new_thresh - threshold = int(old_thresh + 0.5) # range - # Calculate the means of background and object pixels - # Background - mean_back = image[ image <= threshold ].mean() - # Object - mean_obj = image[ image > threshold ].mean() + threshold = old_thresh + tolerance # range + # Calculate the means of background and object pixels + # Background + mean_back = image[image <= threshold].mean() + # Object + mean_obj = image[image > threshold].mean() temp = (mean_back - mean_obj) / (np.log(mean_back) - np.log(mean_obj)) - if temp < 0: # (temp < -2.220446049250313E-16) - new_thresh = int(temp - 0.5) + if temp < 0: + new_thresh = temp - tolerance else: - new_thresh = int(temp + 0.5) - - return threshold + new_thresh = temp + tolerance + + return threshold + offset From 2bc9e82465d1d9f1323cf3f1c93b953fe1cc36ba Mon Sep 17 00:00:00 2001 From: Jeremy Metz Date: Mon, 23 Feb 2015 16:40:38 +0000 Subject: [PATCH 5/5] Minor stylistic changes, removed lena test --- skimage/filters/tests/test_thresholding.py | 4 ---- skimage/filters/thresholding.py | 11 ++++------- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index c360a030..a47b6f20 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/skimage/filters/tests/test_thresholding.py @@ -179,10 +179,6 @@ def test_li_coins_image_as_float(): assert 0.37 < threshold_li(coins) < 0.38 -def test_li_lena_image(): - img = skimage.img_as_ubyte(data.lena()) - assert 127 < threshold_li(img) < 129 - def test_li_astro_image(): img = skimage.img_as_ubyte(data.astronaut()) assert 66 < threshold_li(img) < 68 diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index d5907f8e..2921017b 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -324,7 +324,7 @@ def threshold_li(image): .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" Pattern Recognition, 26(4): 617-625 .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum - Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776 + Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776 .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding Techniques and Quantitative Performance Evaluation" Journal of Electronic Imaging, 13(1): 146-165 @@ -340,10 +340,10 @@ def threshold_li(image): >>> thresh = threshold_li(image) >>> binary = image <= thresh """ - # Requires positive image ( log(mean)) + # Requires positive image (because of log(mean)) offset = image.min() - # Can't use fixed tolerance for float image - imrange = image.max()-offset + # Can not use fixed tolerance for float image + imrange = image.max() - offset image -= offset tolerance = 0.5 * imrange / 256.0 @@ -360,9 +360,7 @@ def threshold_li(image): old_thresh = new_thresh threshold = old_thresh + tolerance # range # Calculate the means of background and object pixels - # Background mean_back = image[image <= threshold].mean() - # Object mean_obj = image[image > threshold].mean() temp = (mean_back - mean_obj) / (np.log(mean_back) - np.log(mean_obj)) @@ -373,4 +371,3 @@ def threshold_li(image): new_thresh = temp + tolerance return threshold + offset -