diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 85043f2c..9907ba6f 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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 diff --git a/bento.info b/bento.info index d181464a..07809644 100644 --- a/bento.info +++ b/bento.info @@ -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 diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py new file mode 100644 index 00000000..65b4571b --- /dev/null +++ b/doc/examples/plot_template.py @@ -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() + diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 4e5d6324..70c154b8 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -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 diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx new file mode 100644 index 00000000..b83761a8 --- /dev/null +++ b/skimage/feature/_template.pyx @@ -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) - ] } + ------------------------------------------------------- + sqrt(sum{ [I(x, y) - ]^2 } sum{ [T(x, y) - ]^2 }) + +where `` is the average of the template, and `` 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) - ]^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 + diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 626b2f5a..13d4fae5 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -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 diff --git a/skimage/feature/template.py b/skimage/feature/template.py new file mode 100644 index 00000000..e5ba8c21 --- /dev/null +++ b/skimage/feature/template.py @@ -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 + diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py new file mode 100644 index 00000000..2eedb922 --- /dev/null +++ b/skimage/feature/tests/test_template.py @@ -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() +