mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-20 12:40:31 +08:00
Add morphological reconstruction with test and example.
This commit is contained in:
Tony S Yu
@@ -0,0 +1,218 @@
|
||||
"""
|
||||
==============
|
||||
Peak detection
|
||||
==============
|
||||
|
||||
Peak detection (a.k.a. spot detection or particle detection) is a common
|
||||
image analysis step. For example, it's used to detect tracer particles in a
|
||||
flow for particle image velocimetry (i.e. PIV) and to identify features in
|
||||
the Hough transform.
|
||||
|
||||
To simplify plotting code, let's define a simple function that creates a new
|
||||
figure on each call, and removes tick labels.
|
||||
"""
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
plt.rcParams['axes.titlesize'] = 10
|
||||
plt.rcParams['font.size'] = 10
|
||||
def imshow(image, **kwargs):
|
||||
plt.figure(figsize=(2.5, 2.5))
|
||||
plt.imshow(image, **kwargs)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
To explore different peak detection techniques, we use an image of circles with
|
||||
added noise:
|
||||
"""
|
||||
|
||||
from skimage import data
|
||||
img = data.load('noisy_circles.jpg')
|
||||
imshow(img)
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
This image is noisy and has uneven background illumination. The peaks in the
|
||||
image, while readily identified by eye, can be tricky to find algorithmically.
|
||||
The first thing we need to do is remove the high-frequency noise; this can
|
||||
be done with a simple Gaussian filter.
|
||||
"""
|
||||
|
||||
import scipy.ndimage as ndimg
|
||||
img_smooth = ndimg.gaussian_filter(img, 3)
|
||||
|
||||
imshow(img_smooth)
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Thresholding
|
||||
============
|
||||
|
||||
One way to extract the background is to threshold the image.
|
||||
"""
|
||||
|
||||
thresh_value = 100
|
||||
background = img_smooth.copy()
|
||||
background[img_smooth > thresh_value] = 0
|
||||
peaks = img_smooth - background
|
||||
|
||||
"""
|
||||
Here, all pixels values below the threshold value are subtracted from the
|
||||
image. The resulting background image and the extracted peaks are shown below.
|
||||
"""
|
||||
|
||||
imshow(background, vmin=0, vmax=255)
|
||||
plt.title("background image (thresholding)")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
imshow(peaks, vmin=0, vmax=255)
|
||||
plt.title("peaks (thresholding)")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Because of uneven illumination, peaks on the right bleed into each other.
|
||||
Increasing the threshold will fix this problem, but it will also cause some
|
||||
peaks on the left to go undetected.
|
||||
|
||||
Morphological reconstruction
|
||||
============================
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
img_r = np.int32(img_smooth)
|
||||
|
||||
import skimage.morphology as morph
|
||||
h = 20
|
||||
rec = morph.reconstruction(img_r-h, img_r)
|
||||
|
||||
imshow(img_r, vmin=0, vmax=255)
|
||||
plt.title("original (smoothed) image")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
imshow(rec, vmin=0, vmax=255)
|
||||
plt.title("background image (reconstruction)")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
This reconstructed image looks pretty much like the original, except that the
|
||||
peaks in the image are truncated. The reconstructed image can then be
|
||||
subtracted from the original image to reveal the peaks of the image.
|
||||
"""
|
||||
|
||||
imshow(img_r-rec)
|
||||
plt.title("h-dome of image")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The result is known as the h-dome transformation [2]_, which extracts peaks of
|
||||
height `h` from the original image. To better understand what's going on,
|
||||
let's take a 1D slice along the middle of the image (cutting through peaks in
|
||||
the image).
|
||||
"""
|
||||
|
||||
img_slice = img_r[99:100, :]
|
||||
rec_slice = morph.reconstruction(img_slice-h, img_slice)
|
||||
|
||||
"""
|
||||
Plotting the reconstructed image (slice) next to the original image and the
|
||||
seed image shed light on the reconstruction process
|
||||
"""
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.plot(img_slice[0], 'k', label='original image')
|
||||
plt.plot(img_slice[0]-h, '0.5', label='seed image')
|
||||
plt.plot(rec_slice[0], 'r', label='reconstructed')
|
||||
plt.title("image slice")
|
||||
plt.xlabel('x')
|
||||
plt.ylabel('intensity')
|
||||
plt.legend()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Here, you see that morphological reconstruction dilates the seed image (i.e.
|
||||
the `h`-shifted image) until it intersects the mask (original image). Note that
|
||||
the peaks in the original image have very different intensity values (e.g. the
|
||||
peak at x=200 and x=100 differ by about 80). Subtracting the reconstructed
|
||||
image from the original image gives peaks of roughly equal intensity. Thus, the
|
||||
h-dome transformation is quiet effective at removing uneven, dark backgrounds
|
||||
from bright features. The inverse operation---the h-basin
|
||||
transformation---should be used when removing bright backgrounds from dark
|
||||
features.
|
||||
|
||||
|
||||
White tophat
|
||||
============
|
||||
"""
|
||||
selem = morph.disk(10)
|
||||
img_t = np.uint8(img_smooth)
|
||||
opening = morph.greyscale_open(img_t, selem)
|
||||
top_hat = img_t - opening
|
||||
|
||||
imshow(opening, vmin=0, vmax=255)
|
||||
plt.title("Greyscale opening of image")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
|
||||
imshow(top_hat)
|
||||
plt.title("Tophat with disk of r = 10")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
selem = morph.disk(5)
|
||||
top_hat = morph.greyscale_white_top_hat(img_t, selem)
|
||||
|
||||
imshow(top_hat)
|
||||
plt.title("Tophat with disk of r = 5")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
selem = morph.square(20)
|
||||
opening = morph.greyscale_open(img_t, selem)
|
||||
# scikit's top hat filter uses uint8 and doesn't check for over(under)flow.
|
||||
mask = opening > img_t
|
||||
opening[mask] = img_t[mask]
|
||||
top_hat = img_t - opening
|
||||
|
||||
imshow(opening, vmin=0, vmax=255)
|
||||
plt.title("Greyscale opening of image")
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
imshow(top_hat)
|
||||
plt.title("Tophat with square of w = 10")
|
||||
|
||||
plt.show()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
|
||||
References
|
||||
==========
|
||||
|
||||
.. [1] Crocker and Grier, Journal of Colloid and Interface Science (1996)
|
||||
.. [2] Vincent, L., IEEE Transactions on Image Processing (1993)
|
||||
|
||||
"""
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 26 KiB |
@@ -4,3 +4,4 @@ from .ccomp import label
|
||||
from .watershed import watershed, is_local_maximum
|
||||
from .skeletonize import skeletonize, medial_axis
|
||||
from .convex_hull import convex_hull_image
|
||||
from .greyreconstruct import reconstruction
|
||||
|
||||
@@ -0,0 +1,85 @@
|
||||
"""
|
||||
`reconstruction_loop` originally part of CellProfiler, code licensed under both GPL and BSD licenses.
|
||||
|
||||
Website: http://www.cellprofiler.org
|
||||
Copyright (c) 2003-2009 Massachusetts Institute of Technology
|
||||
Copyright (c) 2009-2011 Broad Institute
|
||||
All rights reserved.
|
||||
Original author: Lee Kamentsky
|
||||
|
||||
"""
|
||||
|
||||
from __future__ import division
|
||||
import numpy as np
|
||||
|
||||
cimport numpy as np
|
||||
cimport cython
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1,
|
||||
negative_indices = False,
|
||||
mode = 'c'] avalues,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices = False,
|
||||
mode = 'c'] aprev,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices = False,
|
||||
mode = 'c'] anext,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices = False,
|
||||
mode = 'c'] astrides,
|
||||
np.int32_t current,
|
||||
int image_stride):
|
||||
"""The inner loop for reconstruction"""
|
||||
cdef:
|
||||
np.int32_t neighbor
|
||||
np.uint32_t neighbor_value
|
||||
np.uint32_t current_value
|
||||
np.uint32_t mask_value
|
||||
np.int32_t link
|
||||
int i
|
||||
np.int32_t nprev
|
||||
np.int32_t nnext
|
||||
int nstrides = astrides.shape[0]
|
||||
np.uint32_t *values = <np.uint32_t *>(avalues.data)
|
||||
np.int32_t *prev = <np.int32_t *>(aprev.data)
|
||||
np.int32_t *next = <np.int32_t *>(anext.data)
|
||||
np.int32_t *strides = <np.int32_t *>(astrides.data)
|
||||
|
||||
while current != -1:
|
||||
if current < image_stride:
|
||||
current_value = values[current]
|
||||
if current_value == 0:
|
||||
break
|
||||
for i in range(nstrides):
|
||||
neighbor = current + strides[i]
|
||||
neighbor_value = values[neighbor]
|
||||
# Only do neighbors less than the current value
|
||||
if neighbor_value < current_value:
|
||||
mask_value = values[neighbor + image_stride]
|
||||
# Only do neighbors less than the mask value
|
||||
if neighbor_value < mask_value:
|
||||
# Raise the neighbor to the mask value if
|
||||
# the mask is less than current
|
||||
if mask_value < current_value:
|
||||
link = neighbor + image_stride
|
||||
values[neighbor] = mask_value
|
||||
else:
|
||||
link = current
|
||||
values[neighbor] = current_value
|
||||
# unlink the neighbor
|
||||
nprev = prev[neighbor]
|
||||
nnext = next[neighbor]
|
||||
next[nprev] = nnext
|
||||
if nnext != -1:
|
||||
prev[nnext] = nprev
|
||||
# link the neighbor after the link
|
||||
nnext = next[link]
|
||||
next[neighbor] = nnext
|
||||
prev[neighbor] = link
|
||||
if nnext >= 0:
|
||||
prev[nnext] = neighbor
|
||||
next[link] = neighbor
|
||||
current = next[current]
|
||||
|
||||
@@ -0,0 +1,150 @@
|
||||
"""
|
||||
`reconstruction` originally part of CellProfiler, code licensed under both GPL and BSD licenses.
|
||||
|
||||
Website: http://www.cellprofiler.org
|
||||
Copyright (c) 2003-2009 Massachusetts Institute of Technology
|
||||
Copyright (c) 2009-2011 Broad Institute
|
||||
All rights reserved.
|
||||
Original author: Lee Kamentsky
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
from skimage.filter.rank_order import rank_order
|
||||
|
||||
|
||||
def reconstruction(image, mask, selem=None, offset=None):
|
||||
"""Perform a morphological reconstruction of the image.
|
||||
|
||||
Reconstruction requires a "seed" image and a "mask" image. The seed image
|
||||
gets dilated until it is constrained by the mask. The "seed" and "mask"
|
||||
images will be the minimum and maximum possible values of the reconstructed
|
||||
image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
The seed image.
|
||||
|
||||
mask : ndarray
|
||||
The maximum allowed value at each point.
|
||||
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
|
||||
Returns
|
||||
-------
|
||||
reconstructed : ndarray
|
||||
The result of morphological reconstruction.
|
||||
|
||||
Notes
|
||||
-----
|
||||
The algorithm is taken from:
|
||||
Robinson, "Efficient morphological reconstruction: a downhill filter",
|
||||
Pattern Recognition Letters 25 (2004) 1759-1767.
|
||||
|
||||
Applications for greyscale reconstruction are discussed in:
|
||||
Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis:
|
||||
Applications and Efficient Algorithms", IEEE Transactions on Image
|
||||
Processing (1993)
|
||||
|
||||
Examples
|
||||
--------
|
||||
Uses for greyscale reconstruction are described in Vincent (1993). For
|
||||
example, let's try to extract the features of an image by subtracting a
|
||||
background image created by reconstruction.
|
||||
|
||||
First, create an image where the "bumps" are the features that
|
||||
we want to extract:
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scikits.image.morphology.grey import grey_reconstruction
|
||||
>>> y, x = np.mgrid[:20:0.5, :20:0.5]
|
||||
>>> bumps = np.sin(x) + np.sin(y)
|
||||
|
||||
To create the background image, set the mask image to the original image,
|
||||
and the seed image to the original image with an intensity offset, `h`.
|
||||
|
||||
>>> h = 0.3
|
||||
>>> seed = bumps - h
|
||||
>>> rec = grey_reconstruction(seed, bumps)
|
||||
|
||||
The resulting reconstructed image looks exactly like the original image,
|
||||
but with the peaks of the bumps cut off. Subtracting this reconstructed
|
||||
image from the original image leaves just the peaks of the bumps
|
||||
|
||||
>>> hdome = bumps - rec
|
||||
|
||||
This operation is known as the h-dome of the image, which leaves features
|
||||
of height `h` in the subtracted image. The h-dome transform, and its
|
||||
inverse h-basin, are analogous to the white top-hat and black top-hat
|
||||
transforms, but don't require a structuring element.
|
||||
|
||||
"""
|
||||
assert tuple(image.shape) == tuple(mask.shape)
|
||||
assert np.all(image <= mask)
|
||||
try:
|
||||
from ._morphrec import reconstruction_loop
|
||||
except ImportError:
|
||||
raise ImportError("_morphrec extension not available.")
|
||||
|
||||
if selem is None:
|
||||
selem = np.ones([3]*image.ndim, bool)
|
||||
else:
|
||||
selem = selem.copy()
|
||||
|
||||
if offset == None:
|
||||
assert all([d % 2 == 1 for d in selem.shape]),\
|
||||
"Footprint dimensions must all be odd"
|
||||
offset = np.array([d/2 for d in selem.shape])
|
||||
# Cross out the center of the selem
|
||||
selem[[slice(d,d+1) for d in offset]] = False
|
||||
#
|
||||
# Construct an array that's padded on the edges so we can ignore boundaries
|
||||
# The array is a dstack of the image and the mask; this lets us interleave
|
||||
# image and mask pixels when sorting which makes list manipulations easier
|
||||
#
|
||||
padding = (np.array(selem.shape)/2).astype(int)
|
||||
dims = np.zeros(image.ndim+1,int)
|
||||
dims[1:] = np.array(image.shape)+2*padding
|
||||
dims[0] = 2
|
||||
inside_slices = [slice(p,-p) for p in padding]
|
||||
values = np.ones(dims)*np.min(image)
|
||||
values[[0]+inside_slices] = image
|
||||
values[[1]+inside_slices] = mask
|
||||
#
|
||||
# Create a list of strides across the array to get the neighbors
|
||||
# within a flattened array
|
||||
#
|
||||
value_stride = np.array(values.strides[1:]) / values.dtype.itemsize
|
||||
image_stride = values.strides[0] / values.dtype.itemsize
|
||||
selem_mgrid = np.mgrid[[slice(-o,d - o)
|
||||
for d,o in zip(selem.shape,offset)]]
|
||||
selem_offsets = selem_mgrid[:,selem].transpose()
|
||||
strides = np.array([np.sum(value_stride * selem_offset)
|
||||
for selem_offset in selem_offsets], np.int32)
|
||||
values = values.flatten()
|
||||
value_sort = np.lexsort([-values]).astype(np.int32)
|
||||
#
|
||||
# Make a linked list of pixels sorted by value. -1 is the list terminator.
|
||||
#
|
||||
prev = -np.ones(len(values), np.int32)
|
||||
next = -np.ones(len(values), np.int32)
|
||||
prev[value_sort[1:]] = value_sort[:-1]
|
||||
next[value_sort[:-1]] = value_sort[1:]
|
||||
#
|
||||
# Create a rank-order value array so that the Cython inner-loop
|
||||
# can operate on a uniform data type
|
||||
#
|
||||
values, value_map = rank_order(values)
|
||||
current = value_sort[0]
|
||||
|
||||
reconstruction_loop(values, prev, next, strides, current, image_stride)
|
||||
#
|
||||
# Reshape the values array to the shape of the padded image
|
||||
# and return the unpadded portion of that result
|
||||
#
|
||||
values = value_map[values[:image_stride]]
|
||||
values.shape = np.array(image.shape)+2*padding
|
||||
return values[inside_slices]
|
||||
|
||||
@@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None):
|
||||
cython(['_skeletonize.pyx'], working_path=base_path)
|
||||
cython(['_pnpoly.pyx'], working_path=base_path)
|
||||
cython(['_convex_hull.pyx'], working_path=base_path)
|
||||
cython(['_greyreconstruct.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('ccomp', sources=['ccomp.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
@@ -31,6 +32,8 @@ def configuration(parent_package='', top_path=None):
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_convex_hull', sources=['_convex_hull.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,73 @@
|
||||
"""
|
||||
These tests are originally part of CellProfiler, code licensed under both GPL and BSD licenses.
|
||||
|
||||
Website: http://www.cellprofiler.org
|
||||
Copyright (c) 2003-2009 Massachusetts Institute of Technology
|
||||
Copyright (c) 2009-2011 Broad Institute
|
||||
All rights reserved.
|
||||
Original author: Lee Kamentsky
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
from skimage.morphology.greyreconstruct import reconstruction
|
||||
|
||||
|
||||
def test_zeros():
|
||||
"""Test reconstruction with image and mask of zeros"""
|
||||
assert np.all(reconstruction(np.zeros((5, 7)), np.zeros((5, 7))) == 0)
|
||||
|
||||
|
||||
def test_image_equals_mask():
|
||||
"""Test reconstruction where the image and mask are the same"""
|
||||
assert np.all(reconstruction(np.ones((7, 5)), np.ones((7, 5))) == 1)
|
||||
|
||||
|
||||
def test_image_less_than_mask():
|
||||
"""Test reconstruction where the image is uniform and less than mask"""
|
||||
image = np.ones((5, 5))
|
||||
mask = np.ones((5, 5)) * 2
|
||||
assert np.all(reconstruction(image,mask) == 1)
|
||||
|
||||
|
||||
def test_one_image_peak():
|
||||
"""Test reconstruction with one peak pixel"""
|
||||
image = np.ones((5, 5))
|
||||
image[2, 2] = 2
|
||||
mask = np.ones((5, 5)) * 3
|
||||
assert np.all(reconstruction(image,mask) == 2)
|
||||
|
||||
|
||||
def test_two_image_peaks():
|
||||
"""Test reconstruction with two peak pixels isolated by the mask"""
|
||||
image = np.array([[1, 1, 1, 1, 1, 1, 1, 1],
|
||||
[1, 2, 1, 1, 1, 1, 1, 1],
|
||||
[1, 1, 1, 1, 1, 1, 1, 1],
|
||||
[1, 1, 1, 1, 1, 1, 1, 1],
|
||||
[1, 1, 1, 1, 1, 1, 3, 1],
|
||||
[1, 1, 1, 1, 1, 1, 1, 1]])
|
||||
|
||||
mask = np.array([[4, 4, 4, 1, 1, 1, 1, 1],
|
||||
[4, 4, 4, 1, 1, 1, 1, 1],
|
||||
[4, 4, 4, 1, 1, 1, 1, 1],
|
||||
[1, 1, 1, 1, 1, 4, 4, 4],
|
||||
[1, 1, 1, 1, 1, 4, 4, 4],
|
||||
[1, 1, 1, 1, 1, 4, 4, 4]])
|
||||
|
||||
expected = np.array([[2, 2, 2, 1, 1, 1, 1, 1],
|
||||
[2, 2, 2, 1, 1, 1, 1, 1],
|
||||
[2, 2, 2, 1, 1, 1, 1, 1],
|
||||
[1, 1, 1, 1, 1, 3, 3, 3],
|
||||
[1, 1, 1, 1, 1, 3, 3, 3],
|
||||
[1, 1, 1, 1, 1, 3, 3, 3]])
|
||||
assert np.all(reconstruction(image,mask) == expected)
|
||||
|
||||
|
||||
def test_zero_image_one_mask():
|
||||
"""Test reconstruction with an image of all zeros and a mask that's not"""
|
||||
result = reconstruction(np.zeros((10, 10)), np.ones((10, 10)))
|
||||
assert np.all(result == 0)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
Reference in New Issue
Block a user