mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-05 09:42:33 +08:00
Merge pull request #100 from tonysyu/skimage-template
ENH: Add template matching.
This commit is contained in:
Stefan van der Walt
+2
-1
@@ -28,7 +28,7 @@
|
||||
- Tony Yu
|
||||
Reading of paletted images; build, bug and doc fixes.
|
||||
Code to generate skimage logo.
|
||||
Otsu thresholding, histogram equalisation, and more.
|
||||
Otsu thresholding, histogram equalisation, template matching, and more.
|
||||
|
||||
- Zachary Pincus
|
||||
Tracing of low cost paths, FreeImage I/O plugin, iso-contours,
|
||||
@@ -46,6 +46,7 @@
|
||||
|
||||
- Pieter Holtzhausen
|
||||
Incorporating CellProfiler's Sobel edge detector, build and bug fixes.
|
||||
Radon transform, template matching.
|
||||
|
||||
- Emmanuelle Guillart
|
||||
Total variation noise filtering, integration of CellProfiler's
|
||||
|
||||
@@ -43,6 +43,9 @@ Library:
|
||||
Extension: skimage.feature._greycomatrix
|
||||
Sources:
|
||||
skimage/feature/_greycomatrix.pyx
|
||||
Extension: skimage.feature._template
|
||||
Sources:
|
||||
skimage/feature/_template.pyx
|
||||
Extension: skimage.io._plugins._colormixer
|
||||
Sources:
|
||||
skimage/io/_plugins/_colormixer.pyx
|
||||
|
||||
@@ -0,0 +1,56 @@
|
||||
"""
|
||||
=================
|
||||
Template Matching
|
||||
=================
|
||||
|
||||
In this example, we use template matching to identify the occurrence of an
|
||||
image patch (in this case, a sub-image centered on a single coin). Here, we
|
||||
return a single match (the exact same coin), so the maximum value in the
|
||||
``match_template`` result corresponds to the coin location. The other coins
|
||||
look similar, and thus have local maxima; if you expect multiple matches, you
|
||||
should use a proper peak-finding function.
|
||||
|
||||
The ``match_template`` function uses fast, normalized cross-correlation [1]_
|
||||
to find instances of the template in the image. Note that the peaks in the
|
||||
output of ``match_template`` correspond to the origin (i.e. top-left corner) of
|
||||
the template.
|
||||
|
||||
.. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and
|
||||
Magic.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from skimage import data
|
||||
from skimage.feature import match_template
|
||||
|
||||
image = data.coins()
|
||||
coin = image[170:220, 75:130]
|
||||
|
||||
result = match_template(image, coin)
|
||||
ij = np.unravel_index(np.argmax(result), result.shape)
|
||||
x, y = ij[::-1]
|
||||
|
||||
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3))
|
||||
|
||||
ax1.imshow(coin)
|
||||
ax1.set_axis_off()
|
||||
ax1.set_title('template')
|
||||
|
||||
ax2.imshow(image)
|
||||
ax2.set_axis_off()
|
||||
ax2.set_title('image')
|
||||
# highlight matched region
|
||||
hcoin, wcoin = coin.shape
|
||||
rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none')
|
||||
ax2.add_patch(rect)
|
||||
|
||||
ax3.imshow(result)
|
||||
ax3.set_axis_off()
|
||||
ax3.set_title('`match_template`\nresult')
|
||||
# highlight matched region
|
||||
ax3.autoscale(False)
|
||||
ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10)
|
||||
|
||||
plt.show()
|
||||
|
||||
@@ -2,3 +2,4 @@ from .hog import hog
|
||||
from .greycomatrix import greycomatrix, greycoprops
|
||||
from .peak import peak_local_max
|
||||
from .harris import harris
|
||||
from .template import match_template
|
||||
|
||||
@@ -0,0 +1,129 @@
|
||||
"""
|
||||
Template matching using normalized cross-correlation.
|
||||
|
||||
We use fast normalized cross-correlation algorithm (see [1]_ and [2]_) to
|
||||
compute match probability. This algorithm calculates the normalized
|
||||
cross-correlation of an image, `I`, with a template `T` according to the
|
||||
following equation::
|
||||
|
||||
sum{ I(x, y) [T(x, y) - <T>] }
|
||||
-------------------------------------------------------
|
||||
sqrt(sum{ [I(x, y) - <I>]^2 } sum{ [T(x, y) - <T>]^2 })
|
||||
|
||||
where `<T>` is the average of the template, and `<I>` is the average of the
|
||||
image *coincident with the template*, and sums are over the template and the
|
||||
image window coincident with the template. Note that the numerator is simply
|
||||
the cross-correlation of the image and the zero-mean template.
|
||||
|
||||
To speed up calculations, we use summed-area tables (a.k.a. integral images) to
|
||||
quickly calculate sums of image windows inside the loop. This step relies on
|
||||
the following relation (see Eq. 10 of [1])::
|
||||
|
||||
sum{ [I(x, y) - <I>]^2 } =
|
||||
sum{ I^2(x, y) } - [sum{ I(x, y) }]^2 / N_x N_y
|
||||
|
||||
(Without this relation, you would need to subtract each image-window mean from
|
||||
the image window *before* squaring.)
|
||||
|
||||
.. [1] Briechle and Hanebeck, "Template Matching using Fast Normalized
|
||||
Cross Correlation", Proceedings of the SPIE (2001).
|
||||
.. [2] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and
|
||||
Magic.
|
||||
"""
|
||||
import cython
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
from scipy.signal import fftconvolve
|
||||
from skimage.transform import integral
|
||||
|
||||
|
||||
cdef extern from "math.h":
|
||||
float sqrt(float x)
|
||||
float fabs(float x)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef float integrate(np.ndarray[float, ndim=2, mode="c"] sat,
|
||||
int r0, int c0, int r1, int c1):
|
||||
"""
|
||||
Using a summed area table / integral image, calculate the sum
|
||||
over a given window.
|
||||
|
||||
This function is the same as the `integrate` function in
|
||||
`skimage.transform.integrate`, but this Cython version significantly
|
||||
speeds up the code.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
sat : ndarray of float
|
||||
Summed area table / integral image.
|
||||
r0, c0 : int
|
||||
Top-left corner of block to be summed.
|
||||
r1, c1 : int
|
||||
Bottom-right corner of block to be summed.
|
||||
|
||||
Returns
|
||||
-------
|
||||
S : int
|
||||
Sum over the given window.
|
||||
"""
|
||||
cdef float S = 0
|
||||
|
||||
S += sat[r1, c1]
|
||||
|
||||
if (r0 - 1 >= 0) and (c0 - 1 >= 0):
|
||||
S += sat[r0 - 1, c0 - 1]
|
||||
|
||||
if (r0 - 1 >= 0):
|
||||
S -= sat[r0 - 1, c1]
|
||||
|
||||
if (c0 - 1 >= 0):
|
||||
S -= sat[r1, c0 - 1]
|
||||
return S
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def match_template(np.ndarray[float, ndim=2, mode="c"] image,
|
||||
np.ndarray[float, ndim=2, mode="c"] template):
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] corr
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] image_sat
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] image_sqr_sat
|
||||
cdef float template_mean = np.mean(template)
|
||||
cdef float template_ssd
|
||||
cdef float inv_area
|
||||
|
||||
image_sat = integral.integral_image(image)
|
||||
image_sqr_sat = integral.integral_image(image**2)
|
||||
|
||||
template -= template_mean
|
||||
template_ssd = np.sum(template**2)
|
||||
# use inversed area for accuracy
|
||||
inv_area = 1.0 / (template.shape[0] * template.shape[1])
|
||||
|
||||
# when `dtype=float` is used, ascontiguousarray returns ``double``.
|
||||
corr = np.ascontiguousarray(fftconvolve(image,
|
||||
template[::-1, ::-1],
|
||||
mode="valid"),
|
||||
dtype=np.float32)
|
||||
|
||||
cdef int i, j
|
||||
cdef float den, window_sqr_sum, window_mean_sqr, window_sum,
|
||||
# move window through convolution results, normalizing in the process
|
||||
for i in range(corr.shape[0]):
|
||||
for j in range(corr.shape[1]):
|
||||
# subtract 1 because `i_end` and `j_end` are used for indexing into
|
||||
# summed-area table, instead of slicing windows of the image.
|
||||
i_end = i + template.shape[0] - 1
|
||||
j_end = j + template.shape[1] - 1
|
||||
|
||||
window_sum = integrate(image_sat, i, j, i_end, j_end)
|
||||
window_mean_sqr = window_sum * window_sum * inv_area
|
||||
window_sqr_sum = integrate(image_sqr_sat, i, j, i_end, j_end)
|
||||
if window_sqr_sum <= window_mean_sqr:
|
||||
corr[i, j] = 0
|
||||
continue
|
||||
|
||||
den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd)
|
||||
corr[i, j] /= den
|
||||
return corr
|
||||
|
||||
@@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_greycomatrix.pyx'], working_path=base_path)
|
||||
cython(['_template.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_greycomatrix', sources=['_greycomatrix.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_template', sources=['_template.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,84 @@
|
||||
"""template.py - Template matching
|
||||
"""
|
||||
import numpy as np
|
||||
import _template
|
||||
|
||||
from skimage.util.dtype import _convert
|
||||
|
||||
|
||||
def match_template(image, template, pad_input=False):
|
||||
"""Match a template to an image using normalized correlation.
|
||||
|
||||
The output is an array with values between -1.0 and 1.0, which correspond
|
||||
to the probability that the template is found at that position.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like
|
||||
Image to process.
|
||||
template : array_like
|
||||
Template to locate.
|
||||
pad_input : bool
|
||||
If True, pad `image` with image mean so that output is the same size as
|
||||
the image, and output values correspond to the template center.
|
||||
Otherwise, the output is an array with shape `(M - m + 1, N - n + 1)`
|
||||
for an `(M, N)` image and an `(m, n)` template, and matches correspond
|
||||
to origin (top-left corner) of the template.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
Correlation results between -1.0 and 1.0. For an `(M, N)` image and an
|
||||
`(m, n)` template, the `output` is `(M - m + 1, N - n + 1)` when
|
||||
`pad_input = False` and `(M, N)` when `pad_input = True`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> template = np.zeros((3, 3))
|
||||
>>> template[1, 1] = 1
|
||||
>>> print template
|
||||
[[ 0. 0. 0.]
|
||||
[ 0. 1. 0.]
|
||||
[ 0. 0. 0.]]
|
||||
>>> image = np.zeros((6, 6))
|
||||
>>> image[1, 1] = 1
|
||||
>>> image[4, 4] = -1
|
||||
>>> print image
|
||||
[[ 0. 0. 0. 0. 0. 0.]
|
||||
[ 0. 1. 0. 0. 0. 0.]
|
||||
[ 0. 0. 0. 0. 0. 0.]
|
||||
[ 0. 0. 0. 0. 0. 0.]
|
||||
[ 0. 0. 0. 0. -1. 0.]
|
||||
[ 0. 0. 0. 0. 0. 0.]]
|
||||
>>> result = match_template(image, template)
|
||||
>>> print np.round(result, 3)
|
||||
[[ 1. -0.125 0. 0. ]
|
||||
[-0.125 -0.125 0. 0. ]
|
||||
[ 0. 0. 0.125 0.125]
|
||||
[ 0. 0. 0.125 -1. ]]
|
||||
>>> result = match_template(image, template, pad_input=True)
|
||||
>>> print np.round(result, 3)
|
||||
[[-0.125 -0.125 -0.125 0. 0. 0. ]
|
||||
[-0.125 1. -0.125 0. 0. 0. ]
|
||||
[-0.125 -0.125 -0.125 0. 0. 0. ]
|
||||
[ 0. 0. 0. 0.125 0.125 0.125]
|
||||
[ 0. 0. 0. 0.125 -1. 0.125]
|
||||
[ 0. 0. 0. 0.125 0.125 0.125]]
|
||||
"""
|
||||
if np.any(np.less(image.shape, template.shape)):
|
||||
raise ValueError("Image must be larger than template.")
|
||||
image = _convert(image, np.float32)
|
||||
template = _convert(template, np.float32)
|
||||
|
||||
if pad_input:
|
||||
pad_size = tuple(np.array(image.shape) + np.array(template.shape) - 1)
|
||||
pad_image = np.mean(image) * np.ones(pad_size, dtype=np.float32)
|
||||
h, w = image.shape
|
||||
i0, j0 = template.shape
|
||||
i0 /= 2
|
||||
j0 /= 2
|
||||
pad_image[i0:i0+h, j0:j0+w] = image
|
||||
image = pad_image
|
||||
result = _template.match_template(image, template)
|
||||
return result
|
||||
|
||||
@@ -0,0 +1,118 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal as assert_close
|
||||
|
||||
from skimage.morphology import diamond
|
||||
from skimage.feature import match_template, peak_local_max
|
||||
|
||||
|
||||
def test_template():
|
||||
size = 100
|
||||
# Type conversion of image and target not required but prevents warnings.
|
||||
image = np.zeros((400, 400), dtype=np.float32)
|
||||
target = np.tri(size) + np.tri(size)[::-1]
|
||||
target = target.astype(np.float32)
|
||||
target_positions = [(50, 50), (200, 200)]
|
||||
for x, y in target_positions:
|
||||
image[x:x + size, y:y + size] = target
|
||||
np.random.seed(1)
|
||||
image += np.random.randn(400, 400) * 2
|
||||
|
||||
result = match_template(image, target)
|
||||
delta = 5
|
||||
|
||||
positions = peak_local_max(result, min_distance=delta)
|
||||
|
||||
if len(positions) > 2:
|
||||
# Keep the two maximum peaks.
|
||||
intensities = result[tuple(positions.T)]
|
||||
i_maxsort = np.argsort(intensities)[::-1]
|
||||
positions = positions[i_maxsort][:2]
|
||||
|
||||
# Sort so that order matches `target_positions`.
|
||||
positions = positions[np.argsort(positions[:, 0])]
|
||||
|
||||
for xy_target, xy in zip(target_positions, positions):
|
||||
yield assert_close, xy, xy_target
|
||||
|
||||
|
||||
def test_normalization():
|
||||
"""Test that `match_template` gives the correct normalization.
|
||||
|
||||
Normalization gives 1 for a perfect match and -1 for an inverted-match.
|
||||
This test adds positive and negative squares to a zero-array and matches
|
||||
the array with a positive template.
|
||||
"""
|
||||
n = 5
|
||||
N = 20
|
||||
ipos, jpos = (2, 3)
|
||||
ineg, jneg = (12, 11)
|
||||
image = np.zeros((N, N))
|
||||
image[ipos:ipos + n, jpos:jpos + n] = 10
|
||||
image[ineg:ineg + n, jneg:jneg + n] = -10
|
||||
|
||||
# white square with a black border
|
||||
template = np.zeros((n+2, n+2))
|
||||
template[1:1+n, 1:1+n] = 1
|
||||
|
||||
result = match_template(image, template)
|
||||
|
||||
# get the max and min results.
|
||||
sorted_result = np.argsort(result.flat)
|
||||
iflat_min = sorted_result[0]
|
||||
iflat_max = sorted_result[-1]
|
||||
min_result = np.unravel_index(iflat_min, result.shape)
|
||||
max_result = np.unravel_index(iflat_max, result.shape)
|
||||
|
||||
# shift result by 1 because of template border
|
||||
assert np.all((np.array(min_result) + 1) == (ineg, jneg))
|
||||
assert np.all((np.array(max_result) + 1) == (ipos, jpos))
|
||||
|
||||
assert np.allclose(result.flat[iflat_min], -1)
|
||||
assert np.allclose(result.flat[iflat_max], 1)
|
||||
|
||||
|
||||
def test_no_nans():
|
||||
"""Test that `match_template` doesn't return NaN values.
|
||||
|
||||
When image values are only slightly different, floating-point errors can
|
||||
cause a subtraction inside of a square root to go negative (without an
|
||||
explicit check that was added to `match_template`).
|
||||
"""
|
||||
np.random.seed(1)
|
||||
image = 10000 + np.random.normal(size=(20, 20))
|
||||
template = np.ones((6, 6))
|
||||
template[:3, :] = 0
|
||||
result = match_template(image, template)
|
||||
assert not np.any(np.isnan(result))
|
||||
|
||||
|
||||
def test_switched_arguments():
|
||||
image = np.ones((5, 5))
|
||||
template = np.ones((3, 3))
|
||||
np.testing.assert_raises(ValueError, match_template, template, image)
|
||||
|
||||
|
||||
def test_pad_input():
|
||||
template = 10.0 * diamond(2)
|
||||
|
||||
image = np.zeros((9, 19))
|
||||
mid = slice(2, 7)
|
||||
image[mid, :3] = -template[:, -3:] # half min template centered at 0
|
||||
image[mid, 4:9] = template # full max template centered at 6
|
||||
image[mid, -9:-4] = -template # full min template centered at 12
|
||||
image[mid, -3:] = template[:, :3] # half max template centered at 18
|
||||
|
||||
result = match_template(image, template, pad_input=True)
|
||||
|
||||
# get the max and min results.
|
||||
sorted_result = np.argsort(result.flat)
|
||||
i, j = np.unravel_index(sorted_result[:2], result.shape)
|
||||
assert_close(j, (12, 0))
|
||||
i, j = np.unravel_index(sorted_result[-2:], result.shape)
|
||||
assert_close(j, (18, 6))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
|
||||
Reference in New Issue
Block a user