{ "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": [ "%matplotlib inline" ], "metadata": {} }, { "source": "
\n

Rank filters

\n

Rank filters are non-linear filters using the local gray-level ordering to\ncompute the filtered value. This ensemble of filters share a common base: the\nlocal gray-level histogram is computed on the neighborhood of a pixel (defined\nby a 2-D structuring element). If the filtered value is taken as the middle\nvalue of the histogram, we get the classical median filter.

\n

Rank filters can be used for several purposes such as:

\n\n

Some well known filters are specific cases of rank filters [1] e.g.\nmorphological dilation, morphological erosion, median filters.

\n

In this example, we will see how to filter a gray-level image using some of the\nlinear and non-linear filters available in skimage. We use the camera image\nfrom skimage.data for all comparisons.

\n\n\n\n\n\n
[1]Pierre Soille, On morphological operators based on rank filters, Pattern\nRecognition 35 (2002) 527-535.
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nfrom skimage import img_as_ubyte\nfrom skimage import data\n\nnoisy_image = img_as_ubyte(data.camera())\nhist = np.histogram(noisy_image, bins=np.arange(0, 256))\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))\nax1.imshow(noisy_image, interpolation='nearest', cmap=plt.cm.gray)\nax1.axis('off')\nax2.plot(hist[1][:-1], hist[0], lw=2)\nax2.set_title('Histogram of grey values')", "metadata": {} }, { "source": "
\n

Noise removal

\n

Some noise is added to the image, 1% of pixels are randomly set to 255, 1% are\nrandomly set to 0. The median filter is applied to remove the noise.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import median\nfrom skimage.morphology import disk\n\nnoise = np.random.random(noisy_image.shape)\nnoisy_image = img_as_ubyte(data.camera())\nnoisy_image[noise > 0.99] = 255\nnoisy_image[noise < 0.01] = 0\n\nfig, ax = plt.subplots(2, 2, figsize=(10, 7))\nax1, ax2, ax3, ax4 = ax.ravel()\n\nax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)\nax1.set_title('Noisy image')\nax1.axis('off')\n\nax2.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255, cmap=plt.cm.gray)\nax2.set_title('Median $r=1$')\nax2.axis('off')\n\nax3.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255, cmap=plt.cm.gray)\nax3.set_title('Median $r=5$')\nax3.axis('off')\n\nax4.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255, cmap=plt.cm.gray)\nax4.set_title('Median $r=20$')\nax4.axis('off')", "metadata": {} }, { "source": "
\n

The added noise is efficiently removed, as the image defaults are small (1\npixel wide), a small filter radius is sufficient. As the radius is increasing,\nobjects with bigger sizes are filtered as well, such as the camera tripod. The\nmedian filter is often used for noise removal because borders are preserved and\ne.g. salt and pepper noise typically does not distort the gray-level.

\n
\n

Image smoothing

\n

The example hereunder shows how a local mean filter smooths the camera man\nimage.

\n
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import mean\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7])\n\nloc_mean = mean(noisy_image, disk(10))\n\nax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(loc_mean, vmin=0, vmax=255, cmap=plt.cm.gray)\nax2.set_title('Local mean $r=10$')\nax2.axis('off')", "metadata": {} }, { "source": "
\n

One may be interested in smoothing an image while preserving important borders\n(median filters already achieved this), here we use the bilateral filter\nthat restricts the local neighborhood to pixel having a gray-level similar to\nthe central one.

\n
\n

Note

\n

A different implementation is available for color images in\nskimage.filter.denoise_bilateral.

\n
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import mean_bilateral\n\nnoisy_image = img_as_ubyte(data.camera())\n\nbilat = mean_bilateral(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10)\n\nfig, ax = plt.subplots(2, 2, figsize=(10, 7))\nax1, ax2, ax3, ax4 = ax.ravel()\n\nax1.imshow(noisy_image, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(bilat, cmap=plt.cm.gray)\nax2.set_title('Bilateral mean')\nax2.axis('off')\n\nax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)\nax3.axis('off')\n\nax4.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)\nax4.axis('off')", "metadata": {} }, { "source": "
\n

One can see that the large continuous part of the image (e.g. sky) is smoothed\nwhereas other details are preserved.

\n
\n

Contrast enhancement

\n

We compare here how the global histogram equalization is applied locally.

\n

The equalized image [2] has a roughly linear cumulative distribution function\nfor each pixel neighborhood. The local version [3] of the histogram\nequalization emphasizes every local gray-level variations.

\n\n\n\n\n\n
[2]http://en.wikipedia.org/wiki/Histogram_equalization
\n\n\n\n\n\n
[3]http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
\n
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage import exposure\nfrom skimage.filter import rank\n\nnoisy_image = img_as_ubyte(data.camera())\n\n# equalize globally and locally\nglob = exposure.equalize_hist(noisy_image) * 255\nloc = rank.equalize(noisy_image, disk(20))\n\n# extract histogram for each image\nhist = np.histogram(noisy_image, bins=np.arange(0, 256))\nglob_hist = np.histogram(glob, bins=np.arange(0, 256))\nloc_hist = np.histogram(loc, bins=np.arange(0, 256))\n\nfig, ax = plt.subplots(3, 2, figsize=(10, 10))\nax1, ax2, ax3, ax4, ax5, ax6 = ax.ravel()\n\nax1.imshow(noisy_image, interpolation='nearest', cmap=plt.cm.gray)\nax1.axis('off')\n\nax2.plot(hist[1][:-1], hist[0], lw=2)\nax2.set_title('Histogram of gray values')\n\nax3.imshow(glob, interpolation='nearest', cmap=plt.cm.gray)\nax3.axis('off')\n\nax4.plot(glob_hist[1][:-1], glob_hist[0], lw=2)\nax4.set_title('Histogram of gray values')\n\nax5.imshow(loc, interpolation='nearest', cmap=plt.cm.gray)\nax5.axis('off')\n\nax6.plot(loc_hist[1][:-1], loc_hist[0], lw=2)\nax6.set_title('Histogram of gray values')", "metadata": {} }, { "source": "
\n

Another way to maximize the number of gray-levels used for an image is to apply\na local auto-leveling, i.e. the gray-value of a pixel is proportionally\nremapped between local minimum and local maximum.

\n

The following example shows how local auto-level enhances the camara man\npicture.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import autolevel\n\nnoisy_image = img_as_ubyte(data.camera())\n\nauto = autolevel(noisy_image.astype(np.uint16), disk(20))\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7])\n\nax1.imshow(noisy_image, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(auto, cmap=plt.cm.gray)\nax2.set_title('Local autolevel')\nax2.axis('off')", "metadata": {} }, { "source": "
\n

This filter is very sensitive to local outliers, see the little white spot in\nthe left part of the sky. This is due to a local maximum which is very high\ncomparing to the rest of the neighborhood. One can moderate this using the\npercentile version of the auto-level filter which uses given percentiles (one\ninferior, one superior) in place of local minimum and maximum. The example\nbelow illustrates how the percentile parameters influence the local auto-level\nresult.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import autolevel_percentile\n\nimage = data.camera()\n\nselem = disk(20)\nloc_autolevel = autolevel(image, selem=selem)\nloc_perc_autolevel0 = autolevel_percentile(image, selem=selem, p0=.00, p1=1.0)\nloc_perc_autolevel1 = autolevel_percentile(image, selem=selem, p0=.01, p1=.99)\nloc_perc_autolevel2 = autolevel_percentile(image, selem=selem, p0=.05, p1=.95)\nloc_perc_autolevel3 = autolevel_percentile(image, selem=selem, p0=.1, p1=.9)\n\nfig, axes = plt.subplots(nrows=3, figsize=(7, 8))\nax0, ax1, ax2 = axes\nplt.gray()\n\nax0.imshow(np.hstack((image, loc_autolevel)), cmap=plt.cm.gray)\nax0.set_title('Original / auto-level')\n\nax1.imshow(\n np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255)\nax1.set_title('Percentile auto-level 0%,1%')\nax2.imshow(\n np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255)\nax2.set_title('Percentile auto-level 5% and 10%')\n\nfor ax in axes:\n ax.axis('off')", "metadata": {} }, { "source": "
\n

The morphological contrast enhancement filter replaces the central pixel by the\nlocal maximum if the original pixel value is closest to local maximum,\notherwise by the minimum local.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import enhance_contrast\n\nnoisy_image = img_as_ubyte(data.camera())\n\nenh = enhance_contrast(noisy_image, disk(5))\n\nfig, ax = plt.subplots(2, 2, figsize=[10, 7])\nax1, ax2, ax3, ax4 = ax.ravel()\n\nax1.imshow(noisy_image, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(enh, cmap=plt.cm.gray)\nax2.set_title('Local morphological contrast enhancement')\nax2.axis('off')\n\nax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)\nax3.axis('off')\n\nax4.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)\nax4.axis('off')", "metadata": {} }, { "source": "
\n

The percentile version of the local morphological contrast enhancement uses\npercentile p0 and p1 instead of the local minimum and maximum.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import enhance_contrast_percentile\n\nnoisy_image = img_as_ubyte(data.camera())\n\npenh = enhance_contrast_percentile(noisy_image, disk(5), p0=.1, p1=.9)\n\nfig, ax = plt.subplots(2, 2, figsize=[10, 7])\nax1, ax2, ax3, ax4 = ax.ravel()\n\nax1.imshow(noisy_image, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(penh, cmap=plt.cm.gray)\nax2.set_title('Local percentile morphological\\n contrast enhancement')\nax2.axis('off')\n\nax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)\nax3.axis('off')\n\nax4.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)\nax4.axis('off')", "metadata": {} }, { "source": "
\n

Image threshold

\n

The Otsu threshold [1]_ method can be applied locally using the local gray-\nlevel distribution. In the example below, for each pixel, an "optimal"\nthreshold is determined by maximizing the variance between two classes of\npixels of the local neighborhood defined by a structuring element.

\n

The example compares the local threshold with the global threshold\nskimage.filter.threshold_otsu.

\n
\n

Note

\n

Local is much slower than global thresholding. A function for global Otsu\nthresholding can be found in : skimage.filter.threshold_otsu.

\n
\n\n\n\n\n\n
[4]http://en.wikipedia.org/wiki/Otsu's_method
\n
\n

Docutils System Messages

\n
\n

System Message: ERROR/3 (<string>, line 7); backlink

\nUnknown target name: "1".
\n
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import otsu\nfrom skimage.filter import threshold_otsu\n\np8 = data.page()\n\nradius = 10\nselem = disk(radius)\n\n# t_loc_otsu is an image\nt_loc_otsu = otsu(p8, selem)\nloc_otsu = p8 >= t_loc_otsu\n\n# t_glob_otsu is a scalar\nt_glob_otsu = threshold_otsu(p8)\nglob_otsu = p8 >= t_glob_otsu\n\nfig, ax = plt.subplots(2, 2)\nax1, ax2, ax3, ax4 = ax.ravel()\n\nfig.colorbar(ax1.imshow(p8, cmap=plt.cm.gray), ax=ax1)\nax1.set_title('Original')\nax1.axis('off')\n\nfig.colorbar(ax2.imshow(t_loc_otsu, cmap=plt.cm.gray), ax=ax2)\nax2.set_title('Local Otsu ($r=%d$)' % radius)\nax2.axis('off')\n\nax3.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)\nax3.set_title('Original >= local Otsu' % t_glob_otsu)\nax3.axis('off')\n\nax4.imshow(glob_otsu, cmap=plt.cm.gray)\nax4.set_title('Global Otsu ($t=%d$)' % t_glob_otsu)\nax4.axis('off')", "metadata": {} }, { "source": "
\n

The following example shows how local Otsu thresholding handles a global level\nshift applied to a synthetic image.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nn = 100\ntheta = np.linspace(0, 10 * np.pi, n)\nx = np.sin(theta)\nm = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8)\n\nradius = 10\nt = rank.otsu(m, disk(radius))\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\n\nax1.imshow(m)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(m >= t, interpolation='nearest')\nax2.set_title('Local Otsu ($r=%d$)' % radius)\nax2.axis('off')", "metadata": {} }, { "source": "
\n

Image morphology

\n

Local maximum and local minimum are the base operators for gray-level\nmorphology.

\n
\n

Note

\n

skimage.dilate and skimage.erode are equivalent filters (see below for\ncomparison).

\n
\n

Here is an example of the classical morphological gray-level filters: opening,\nclosing and morphological gradient.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage.filter.rank import maximum, minimum, gradient\n\nnoisy_image = img_as_ubyte(data.camera())\n\nclosing = maximum(minimum(noisy_image, disk(5)), disk(5))\nopening = minimum(maximum(noisy_image, disk(5)), disk(5))\ngrad = gradient(noisy_image, disk(5))\n\n# display results\nfig, ax = plt.subplots(2, 2, figsize=[10, 7])\nax1, ax2, ax3, ax4 = ax.ravel()\n\nax1.imshow(noisy_image, cmap=plt.cm.gray)\nax1.set_title('Original')\nax1.axis('off')\n\nax2.imshow(closing, cmap=plt.cm.gray)\nax2.set_title('Gray-level closing')\nax2.axis('off')\n\nax3.imshow(opening, cmap=plt.cm.gray)\nax3.set_title('Gray-level opening')\nax3.axis('off')\n\nax4.imshow(grad, cmap=plt.cm.gray)\nax4.set_title('Morphological gradient')\nax4.axis('off')", "metadata": {} }, { "source": "
\n

Feature extraction

\n

Local histograms can be exploited to compute local entropy, which is related to\nthe local image complexity. Entropy is computed using base 2 logarithm i.e. the\nfilter returns the minimum number of bits needed to encode local gray-level\ndistribution.

\n

skimage.rank.entropy returns the local entropy on a given structuring\nelement. The following example shows applies this filter on 8- and 16-bit\nimages.

\n
\n

Note

\n

to better use the available image bit, the function returns 10x entropy for\n8-bit images and 1000x entropy for 16-bit images.

\n
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom skimage import data\nfrom skimage.filter.rank import entropy\nfrom skimage.morphology import disk\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimage = data.camera()\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))\n\nfig.colorbar(ax1.imshow(image, cmap=plt.cm.gray), ax=ax1)\nax1.set_title('Image')\nax1.axis('off')\n\nfig.colorbar(ax2.imshow(entropy(image, disk(5)), cmap=plt.cm.jet), ax=ax2)\nax2.set_title('Entropy')\nax2.axis('off')", "metadata": {} }, { "source": "
\n

Implementation

\n

The central part of the skimage.rank filters is build on a sliding window\nthat updates the local gray-level histogram. This approach limits the algorithm\ncomplexity to O(n) where n is the number of image pixels. The complexity is\nalso limited with respect to the structuring element size.

\n

In the following we compare the performance of different implementations\navailable in skimage.

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfrom time import time\n\nfrom scipy.ndimage.filters import percentile_filter\nfrom skimage.morphology import dilation\nfrom skimage.filter.rank import median, maximum\n\n\ndef exec_and_timeit(func):\n \"\"\"Decorator that returns both function results and execution time.\"\"\"\n def wrapper(*arg):\n t1 = time()\n res = func(*arg)\n t2 = time()\n ms = (t2 - t1) * 1000.0\n return (res, ms)\n return wrapper\n\n\n@exec_and_timeit\ndef cr_med(image, selem):\n return median(image=image, selem=selem)\n\n\n@exec_and_timeit\ndef cr_max(image, selem):\n return maximum(image=image, selem=selem)\n\n\n@exec_and_timeit\ndef cm_dil(image, selem):\n return dilation(image=image, selem=selem)\n\n\n@exec_and_timeit\ndef ndi_med(image, n):\n return percentile_filter(image, 50, size=n * 2 - 1)", "metadata": {} }, { "source": "
\n

Comparison between

\n\n

on increasing structuring element size:

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\na = data.camera()\n\nrec = []\ne_range = range(1, 20, 2)\nfor r in e_range:\n elem = disk(r + 1)\n rc, ms_rc = cr_max(a, elem)\n rcm, ms_rcm = cm_dil(a, elem)\n rec.append((ms_rc, ms_rcm))\n\nrec = np.asarray(rec)\n\nfig, ax = plt.subplots()\nax.set_title('Performance with respect to element size')\nax.set_ylabel('Time (ms)')\nax.set_xlabel('Element radius')\nax.plot(e_range, rec)\nax.legend(['filter.rank.maximum', 'morphology.dilate'])", "metadata": {} }, { "source": "
\n

and increasing image size:

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nr = 9\nelem = disk(r + 1)\n\nrec = []\ns_range = range(100, 1000, 100)\nfor s in s_range:\n a = (np.random.random((s, s)) * 256).astype(np.uint8)\n (rc, ms_rc) = cr_max(a, elem)\n (rcm, ms_rcm) = cm_dil(a, elem)\n rec.append((ms_rc, ms_rcm))\n\nrec = np.asarray(rec)\n\nfig, ax = plt.subplots()\nax.set_title('Performance with respect to image size')\nax.set_ylabel('Time (ms)')\nax.set_xlabel('Image size')\nax.plot(s_range, rec)\nax.legend(['filter.rank.maximum', 'morphology.dilate'])", "metadata": {} }, { "source": "
\n

Comparison between:

\n\n

on increasing structuring element size:

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\na = data.camera()\n\nrec = []\ne_range = range(2, 30, 4)\nfor r in e_range:\n elem = disk(r + 1)\n rc, ms_rc = cr_med(a, elem)\n rndi, ms_ndi = ndi_med(a, r)\n rec.append((ms_rc, ms_ndi))\n\nrec = np.asarray(rec)\n\nfig, ax = plt.subplots()\nax.set_title('Performance with respect to element size')\nax.plot(e_range, rec)\nax.legend(['filter.rank.median', 'scipy.ndimage.percentile'])\nax.set_ylabel('Time (ms)')\nax.set_xlabel('Element radius')", "metadata": {} }, { "source": "
\n

Comparison of outcome of the three methods:

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nfig, ax = plt.subplots()\nax.imshow(np.hstack((rc, rndi)))\nax.set_title('filter.rank.median vs. scipy.ndimage.percentile')\nax.axis('off')", "metadata": {} }, { "source": "
\n

and increasing image size:

\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nr = 9\nelem = disk(r + 1)\n\nrec = []\ns_range = [100, 200, 500, 1000]\nfor s in s_range:\n a = (np.random.random((s, s)) * 256).astype(np.uint8)\n (rc, ms_rc) = cr_med(a, elem)\n rndi, ms_ndi = ndi_med(a, r)\n rec.append((ms_rc, ms_ndi))\n\nrec = np.asarray(rec)\n\nfig, ax = plt.subplots()\nax.set_title('Performance with respect to image size')\nax.plot(s_range, rec)\nax.legend(['filter.rank.median', 'scipy.ndimage.percentile'])\nax.set_ylabel('Time (ms)')\nax.set_xlabel('Image size')", "metadata": {} }, { "source": "
\n
\n", "cell_type": "markdown", "metadata": {} }, { "cell_type": "code", "language": "python", "outputs": [], "collapsed": false, "input": "\nplt.show()", "metadata": {} } ], "metadata": {} } ], "metadata": { "name": "" } }