3D Rank Filters (scikit-image#4435)

* made core_3d * fixed reshaping error * Removed extraneous functions * PEP8 compliance * added tests * documentation, benchmarks * Added to release_dev.rst * added otsu * add otsu * added test for otsu * example * changes to documentation * changed contrast in brain image, html updates * added html demo for otsu 3D * autolevel implementation, same as new_rank_dev * fix bug and start adding test * modify loading 3d rank test and added file_hash * revert previous file_hash change and import through os instead of fetch * added data to registry and deleted * implemented and created basic tests for all generic rank filters except windowed_histogram since it is the only filter that applies vector per pixel; will work on percentiles and bilaterals after fixing the reviewed code * Apply suggestions from code review suggestions up to __init__.py Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu> * Simple code review fixes up to __init__.py * Apply suggestions from code review The ones after __init__.py Co-authored-by: Riadh Fezzani <rfezzani@gmail.com> Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu> * fixed bug and code reviews other than core_cy.pyx * modifying core_cy.pyx with error * Apply suggestions from code review Co-authored-by: Riadh Fezzani <rfezzani@gmail.com> * one line change to core_cy.pyx * change histo type from pointer to memoryview for all 2D and 3D in core, generic, bilateral, percentile * add check for selem/image size and fix bug * test selem of size 1 * test: removing 3D plot_local_equalize * style change and separate core_cy_3d.pyx * BUG: fix indexing order in is_in_mask_3D rank filter function * add plt.show() to equalization example * Update benchmarks/benchmark_rank.py Co-authored-by: Gregory R. Lee <grlee77@gmail.com> * updated benchmark * non-isotropic 3D tests Co-authored-by: Stanley_Wang <wangxi4062@gmail.com> Co-authored-by: Jaewon Chung <jaewonc78@gmail.com> Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu> Co-authored-by: Riadh Fezzani <rfezzani@gmail.com> Co-authored-by: Gregory Lee <grlee77@gmail.com> Co-authored-by: Emmanuelle Gouillart <emma@plot.ly>
cclauss · Sep 27, 2020 · cd5c1a7 · cd5c1a7
1 parent 853c3f8
commit cd5c1a7
Show file tree

Hide file tree

Showing 21 changed files with 1,426 additions and 422 deletions.
diff --git a/.gitignore b/.gitignore
@@ -13,6 +13,7 @@ __pycache__
 *.new
 *.md5
 *.old
+.DS_Store
 .pytest_cache
 .mypy_cache/
 temp.tif

diff --git a/benchmarks/benchmark_rank.py b/benchmarks/benchmark_rank.py
@@ -1,18 +1,31 @@
 import numpy as np
-from skimage import img_as_ubyte
 from skimage.filters import rank
 from skimage.filters.rank import __all__ as all_rank_filters
-from skimage.morphology import grey, disk
+from skimage.filters.rank import __3Dfilters as all_3d_rank_filters
+from skimage.morphology import disk, ball
 
 
 class RankSuite(object):
 
-    param_names = ["filter", "shape"]
+    param_names = ["filter_func", "shape"]
     params = [sorted(all_rank_filters), [(32, 32), (256, 256)]]
 
-    def setup(self, filter, shape):
+    def setup(self, filter_func, shape):
         self.image = np.random.randint(0, 255, size=shape, dtype=np.uint8)
         self.selem = disk(1)
 
-    def time_filter(self, filter, shape):
-        getattr(rank, filter)(self.image, self.selem)
+    def time_filter(self, filter_func, shape):
+        getattr(rank, filter_func)(self.image, self.selem)
+
+
+class Rank3DSuite(object):
+
+    param_names = ["filter3d", "shape3d"]
+    params = [sorted(all_3d_rank_filters), [(32, 32, 32), (128, 128, 128)]]
+
+    def setup(self, filter3d, shape3d):
+        self.volume = np.random.randint(0, 255, size=shape3d, dtype=np.uint8)
+        self.selem_3d = ball(1)
+
+    def time_3d_filters(self, filter3d, shape3d):
+        getattr(rank, filter3d)(self.volume, self.selem_3d)
diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py
@@ -64,7 +64,7 @@
 # noise.
 
 from skimage.filters.rank import median
-from skimage.morphology import disk
+from skimage.morphology import disk, ball
 
 noise = np.random.random(noisy_image.shape)
 noisy_image = img_as_ubyte(data.camera())
@@ -353,6 +353,8 @@
 # threshold is determined by maximizing the variance between two classes of
 # pixels of the local neighborhood defined by a structuring element.
 #
+# These algorithms can be used on both 2D and 3D images.
+#
 # The example compares the local threshold with the global threshold
 # :py:func:`skimage.filters.threshold_otsu`.
 #
@@ -366,6 +368,7 @@
 
 from skimage.filters.rank import otsu
 from skimage.filters import threshold_otsu
+from skimage import exposure
 
 p8 = data.page()
 
@@ -401,6 +404,45 @@
 
 plt.tight_layout()
 
+######################################################################
+# The example compares the local threshold with the global threshold in 3D
+
+brain = exposure.rescale_intensity(data.brain().astype(float))
+
+radius = 5
+neighborhood = ball(radius)
+
+# t_loc_otsu is an image
+t_loc_otsu = rank.otsu(brain, neighborhood)
+loc_otsu = brain >= t_loc_otsu
+
+# t_glob_otsu is a scalar
+t_glob_otsu = threshold_otsu(brain)
+glob_otsu = brain >= t_glob_otsu
+
+fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12),
+                         sharex=True, sharey=True)
+ax = axes.ravel()
+
+slice_index = 3
+
+fig.colorbar(ax[0].imshow(brain[slice_index], cmap=plt.cm.gray), ax=ax[0])
+ax[0].set_title('Original')
+
+fig.colorbar(ax[1].imshow(t_loc_otsu[slice_index], cmap=plt.cm.gray), ax=ax[1])
+ax[1].set_title('Local Otsu ($r=%d$)' % radius)
+
+ax[2].imshow(brain[slice_index] >= t_loc_otsu[slice_index], cmap=plt.cm.gray)
+ax[2].set_title('Original >= local Otsu' % t_glob_otsu)
+
+ax[3].imshow(glob_otsu[slice_index], cmap=plt.cm.gray)
+ax[3].set_title('Global Otsu ($t=%d$)' % t_glob_otsu)
+
+for a in ax:
+    a.axis('off')
+
+fig.tight_layout()
+
 ######################################################################
 # The following example shows how local Otsu thresholding handles a global
 # level shift applied to a synthetic image.

diff --git a/doc/examples/color_exposure/plot_local_equalize.py b/doc/examples/color_exposure/plot_local_equalize.py
@@ -13,6 +13,8 @@
 The local version [2]_ of the histogram equalization emphasized every local
 graylevel variations.
 
+These algorithms can be used on both 2D and 3D images.
+
 References
 ----------
 .. [1] https://en.wikipedia.org/wiki/Histogram_equalization
@@ -28,6 +30,7 @@
 from skimage.util import img_as_ubyte
 from skimage import exposure
 from skimage.morphology import disk
+from skimage.morphology import ball
 from skimage.filters import rank
 
 
@@ -93,6 +96,76 @@ def plot_img_and_hist(image, axes, bins=256):
 ax_cdf.set_ylabel('Fraction of total intensity')
 
 
+# prevent overlap of y-axis labels
+fig.tight_layout()
+
+
+######################################################################
+#
+# 3D Equalization
+# ===============
+#
+# 3D Volumes can also be equalized in a similar fashion.
+# Here the histograms are collected from the entire 3D image, but
+# only a single slice is shown for visual inspection.
+
+
+matplotlib.rcParams['font.size'] = 9
+
+
+def plot_img_and_hist(image, axes, bins=256):
+    """Plot an image along with its histogram and cumulative histogram.
+
+    """
+    ax_img, ax_hist = axes
+    ax_cdf = ax_hist.twinx()
+
+    # Display Slice of Image
+    ax_img.imshow(image[0], cmap=plt.cm.gray)
+    ax_img.set_axis_off()
+
+    # Display histogram
+    ax_hist.hist(image.ravel(), bins=bins)
+    ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
+    ax_hist.set_xlabel('Pixel intensity')
+
+    xmin, xmax = dtype_range[image.dtype.type]
+    ax_hist.set_xlim(xmin, xmax)
+
+    # Display cumulative distribution
+    img_cdf, bins = exposure.cumulative_distribution(image, bins)
+    ax_cdf.plot(bins, img_cdf, 'r')
+
+    return ax_img, ax_hist, ax_cdf
+
+
+# Load an example image
+img = img_as_ubyte(data.brain())
+
+# Global equalization
+img_rescale = exposure.equalize_hist(img)
+
+# Local equalization
+neighborhood = ball(3)
+img_eq = rank.equalize(img, selem=neighborhood)
+
+# Display results
+fig, axes = plt.subplots(2, 3, figsize=(8, 5))
+axes[0, 1] = plt.subplot(2, 3, 2, sharex=axes[0, 0], sharey=axes[0, 0])
+axes[0, 2] = plt.subplot(2, 3, 3, sharex=axes[0, 0], sharey=axes[0, 0])
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
+ax_img.set_title('Low contrast image')
+ax_hist.set_ylabel('Number of pixels')
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1])
+ax_img.set_title('Global equalize')
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
+ax_img.set_title('Local equalize')
+ax_cdf.set_ylabel('Fraction of total intensity')
+
+
 # prevent overlap of y-axis labels
 fig.tight_layout()
 plt.show()
diff --git a/doc/release/release_dev.rst b/doc/release/release_dev.rst
@@ -15,6 +15,8 @@ https://scikit-image.org
 New Features
 ------------
 
+- Added 3D support for many filters in skimage.filters.rank.
+
 - A new doc tutorial presenting a cell biology example has been added to the
   gallery (#4648). The scientific content benefited from a much appreciated
   review by Pierre Poulain and Fred Bernard, both assistant professors at

diff --git a/requirements/build.txt b/requirements/build.txt
@@ -1,4 +1,4 @@
-Cython>=0.29.13,!=0.29.18
+Cython>=0.29.21,!=0.29.18
 wheel
 # numpy 1.18.0 breaks builds on MacOSX
 # https://github.com/numpy/numpy/pull/15194

diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py
@@ -22,6 +22,7 @@
            'download_all',
            'astronaut',
            'binary_blobs',
+           'brain',
            'brick',
            'camera',
            'cat',
@@ -1054,3 +1055,25 @@ def lfw_subset():
 
     """
     return np.load(_fetch('data/lfw_subset.npy'))
+
+
+def brain():
+    """Subset of data from the University of North Carolina Volume Rendering
+    Test Data Set.
+
+    The full dataset is available at [1]_.
+
+    Returns
+    -------
+    image : (10, 256, 256) uint16 ndarray
+
+    Notes
+    -----
+    The 3D volume consists of 10 layers from the larger volume.
+
+    References
+    ----------
+    .. [1] https://graphics.stanford.edu/data/voldata/
+
+    """
+    return _load("data/brain.tiff")
diff --git a/skimage/data/_registry.py b/skimage/data/_registry.py
@@ -99,6 +99,7 @@
     "data/astronaut_GRAY_hog_L1.npy": "5d8ab22b166d1dd49c12caeff9d178ed28132efea3852b952e9d75f7f7f94954",
     "data/astronaut_GRAY_hog_L2-Hys.npy": "c4dd6e50d1129aada358311cf8880ce8c775f31e0e550fc322c16e43a96d56fe",
     "data/rank_filter_tests.npz": "efaf5699630f4a53255e91681dc72a965acd4a8aa1f84671c686fb93e7df046d",
+    "data/rank_filters_tests_3d.npz": "1741c2b978424e93558a07d345b2a0d9bfbb33c095c123da147fca066714ab16",
     "data/multi.fits": "5c71a83436762a52b1925f2f0d83881af7765ed50aede155af2800e54bbd5040",
     "data/simple.fits": "cd36087fdbb909b6ba506bbff6bcd4c5f4da3a41862608fbac5e8555ef53d40f",
     "data/palette_color.png": "c4e817035fb9f7730fe95cff1da3866dea01728efc72b6e703d78f7ab9717bdd",
@@ -130,13 +131,16 @@
     "registration/tests/data/TransformedX130Y130.png": "bb10c6ae3f91a313b0ac543efdb7ca69c4b95e55674c65a88472a6c4f4692a25",
     "registration/tests/data/TransformedX75Y75.png": "a1e9ead5f8e4a0f604271e1f9c50e89baf53f068f1d19fab2876af4938e695ea",
     "data/cells.tif": "2120cfe08e0396324793a10a905c9bbcb64b117215eb63b2c24b643e1600c8c9",
+    "data/brain.tiff": "bcdbaf424fbad7b1fb0f855f608c68e5a838f35affc323ff04ea17f678eef5c6",
     "data/kidney.tif": "80c0799bc58b08cf6eaa53ecd202305eb42fd7bc73746cb6c5064dbeae7e8476",
     "data/lily.tif": "395c2f0194c25b9824a8cd79266920362a0816bc9e906dd392adce2d8309af03",
     "data/mitosis.tif": "2751ba667c4067c5d30817cff004aa06f6f6287f1cdbb5b8c9c6a500308cb456",
 }
 
 registry_urls = {
     "data/cells.tif": "https://github.com/scikit-image/skimage-tutorials/raw/master/images/cells.tif",
+    "data/rank_filters_tests_3d.npz": "https://gitlab.com/scikit-image/data/-/raw/master/Tests_besides_Equalize_Otsu/add18_entropy/rank_filters_tests_3d.npz",
+    "data/brain.tiff": "https://gitlab.com/scikit-image/data/-/raw/master/brain.tiff",
     "data/eagle.png": "https://gitlab.com/scikit-image/data/-/raw/master/eagle.png",
     "data/kidney.tif": "https://gitlab.com/scikit-image/data/-/raw/master/kidney-tissue-fluorescence.tif",
     "data/lily.tif": "https://gitlab.com/scikit-image/data/-/raw/master/lily-of-the-valley-fluorescence.tif",

diff --git a/skimage/data/tests/test_data.py b/skimage/data/tests/test_data.py
@@ -153,6 +153,13 @@ def test_cells_3d():
     assert image.shape == (60, 256, 256)
 
 
+def test_brain_3d():
+    """Needs internet connection."""
+    path = fetch('data/brain.tiff')
+    image = io.imread(path)
+    assert image.shape == (10, 256, 256)
+
+
 def test_kidney_3d_multichannel():
     """Test that 3D multichannel image of kidney tissue can be loaded.
 

diff --git a/skimage/filters/rank/__init__.py b/skimage/filters/rank/__init__.py
@@ -1,8 +1,8 @@
-from .generic import (autolevel, bottomhat, equalize, gradient,
+from .generic import (autolevel, bottomhat, equalize, gradient, majority,
                       maximum, mean, geometric_mean, subtract_mean,
                       median, minimum, modal, enhance_contrast, pop,
                       threshold, tophat, noise_filter, entropy, otsu,
-                      sum, windowed_histogram, majority)
+                      sum, windowed_histogram)
 from ._percentile import (autolevel_percentile, gradient_percentile,
                           mean_percentile, subtract_mean_percentile,
                           enhance_contrast_percentile, percentile,
@@ -16,6 +16,7 @@
            'equalize',
            'gradient',
            'gradient_percentile',
+           'majority',
            'maximum',
            'mean',
            'geometric_mean',
@@ -41,5 +42,23 @@
            'entropy',
            'otsu',
            'percentile',
-           'windowed_histogram',
-           'majority']
+           'windowed_histogram']
+
+__3Dfilters = ['autolevel',
+               'equalize',
+               'gradient',
+               'majority',
+               'maximum',
+               'mean',
+               'geometric_mean',
+               'subtract_mean',
+               'median',
+               'minimum',
+               'modal',
+               'enhance_contrast',
+               'pop',
+               'sum',
+               'threshold',
+               'noise_filter',
+               'entropy',
+               'otsu']
diff --git a/skimage/filters/rank/bilateral_cy.pyx b/skimage/filters/rank/bilateral_cy.pyx
@@ -11,7 +11,7 @@ from .core_cy cimport dtype_t, dtype_t_out, _core
 cnp.import_array()
 
 cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
-                              Py_ssize_t* histo,
+                              Py_ssize_t[::1] histo,
                               double pop, dtype_t g,
                               Py_ssize_t n_bins, Py_ssize_t mid_bin,
                               double p0, double p1,
@@ -35,7 +35,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
 
 
 cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
-                             Py_ssize_t* histo,
+                             Py_ssize_t[::1] histo,
                              double pop, dtype_t g,
                              Py_ssize_t n_bins, Py_ssize_t mid_bin,
                              double p0, double p1,
@@ -54,7 +54,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
 
 
 cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
-                             Py_ssize_t* histo,
+                             Py_ssize_t[::1] histo,
                              double pop, dtype_t g,
                              Py_ssize_t n_bins, Py_ssize_t mid_bin,
                              double p0, double p1,

diff --git a/skimage/filters/rank/core_cy.pxd b/skimage/filters/rank/core_cy.pxd
@@ -15,7 +15,7 @@ cdef dtype_t _max(dtype_t a, dtype_t b) nogil
 cdef dtype_t _min(dtype_t a, dtype_t b) nogil
 
 
-cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
+cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t[::1], double,
                             dtype_t, Py_ssize_t, Py_ssize_t, double,
                             double, Py_ssize_t, Py_ssize_t) nogil,
                 dtype_t[:, ::1] image,