scikit-image

1
import math
2

3
import numpy as np
4
from scipy.signal import fftconvolve
5

6
from .._shared.utils import check_nD, _supported_float_type
7

8

9
def _window_sum_2d(image, window_shape):
10
    window_sum = np.cumsum(image, axis=0)
11
    window_sum = window_sum[window_shape[0] : -1] - window_sum[: -window_shape[0] - 1]
12

13
    window_sum = np.cumsum(window_sum, axis=1)
14
    window_sum = (
15
        window_sum[:, window_shape[1] : -1] - window_sum[:, : -window_shape[1] - 1]
16
    )
17

18
    return window_sum
19

20

21
def _window_sum_3d(image, window_shape):
22
    window_sum = _window_sum_2d(image, window_shape)
23

24
    window_sum = np.cumsum(window_sum, axis=2)
25
    window_sum = (
26
        window_sum[:, :, window_shape[2] : -1]
27
        - window_sum[:, :, : -window_shape[2] - 1]
28
    )
29

30
    return window_sum
31

32

33
def match_template(
34
    image, template, pad_input=False, mode='constant', constant_values=0
35
):
36
    """Match a template to a 2-D or 3-D image using normalized correlation.
37

38
    The output is an array with values between -1.0 and 1.0. The value at a
39
    given position corresponds to the correlation coefficient between the image
40
    and the template.
41

42
    For `pad_input=True` matches correspond to the center and otherwise to the
43
    top-left corner of the template. To find the best match you must search for
44
    peaks in the response (output) image.
45

46
    Parameters
47
    ----------
48
    image : (M, N[, P]) array
49
        2-D or 3-D input image.
50
    template : (m, n[, p]) array
51
        Template to locate. It must be `(m <= M, n <= N[, p <= P])`.
52
    pad_input : bool
53
        If True, pad `image` so that output is the same size as the image, and
54
        output values correspond to the template center. Otherwise, the output
55
        is an array with shape `(M - m + 1, N - n + 1)` for an `(M, N)` image
56
        and an `(m, n)` template, and matches correspond to origin
57
        (top-left corner) of the template.
58
    mode : see `numpy.pad`, optional
59
        Padding mode.
60
    constant_values : see `numpy.pad`, optional
61
        Constant values used in conjunction with ``mode='constant'``.
62

63
    Returns
64
    -------
65
    output : array
66
        Response image with correlation coefficients.
67

68
    Notes
69
    -----
70
    Details on the cross-correlation are presented in [1]_. This implementation
71
    uses FFT convolutions of the image and the template. Reference [2]_
72
    presents similar derivations but the approximation presented in this
73
    reference is not used in our implementation.
74

75
    References
76
    ----------
77
    .. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light
78
           and Magic.
79
    .. [2] Briechle and Hanebeck, "Template Matching using Fast Normalized
80
           Cross Correlation", Proceedings of the SPIE (2001).
81
           :DOI:`10.1117/12.421129`
82

83
    Examples
84
    --------
85
    >>> template = np.zeros((3, 3))
86
    >>> template[1, 1] = 1
87
    >>> template
88
    array([[0., 0., 0.],
89
           [0., 1., 0.],
90
           [0., 0., 0.]])
91
    >>> image = np.zeros((6, 6))
92
    >>> image[1, 1] = 1
93
    >>> image[4, 4] = -1
94
    >>> image
95
    array([[ 0.,  0.,  0.,  0.,  0.,  0.],
96
           [ 0.,  1.,  0.,  0.,  0.,  0.],
97
           [ 0.,  0.,  0.,  0.,  0.,  0.],
98
           [ 0.,  0.,  0.,  0.,  0.,  0.],
99
           [ 0.,  0.,  0.,  0., -1.,  0.],
100
           [ 0.,  0.,  0.,  0.,  0.,  0.]])
101
    >>> result = match_template(image, template)
102
    >>> np.round(result, 3)
103
    array([[ 1.   , -0.125,  0.   ,  0.   ],
104
           [-0.125, -0.125,  0.   ,  0.   ],
105
           [ 0.   ,  0.   ,  0.125,  0.125],
106
           [ 0.   ,  0.   ,  0.125, -1.   ]])
107
    >>> result = match_template(image, template, pad_input=True)
108
    >>> np.round(result, 3)
109
    array([[-0.125, -0.125, -0.125,  0.   ,  0.   ,  0.   ],
110
           [-0.125,  1.   , -0.125,  0.   ,  0.   ,  0.   ],
111
           [-0.125, -0.125, -0.125,  0.   ,  0.   ,  0.   ],
112
           [ 0.   ,  0.   ,  0.   ,  0.125,  0.125,  0.125],
113
           [ 0.   ,  0.   ,  0.   ,  0.125, -1.   ,  0.125],
114
           [ 0.   ,  0.   ,  0.   ,  0.125,  0.125,  0.125]])
115
    """
116
    check_nD(image, (2, 3))
117

118
    if image.ndim < template.ndim:
119
        raise ValueError(
120
            "Dimensionality of template must be less than or "
121
            "equal to the dimensionality of image."
122
        )
123
    if np.any(np.less(image.shape, template.shape)):
124
        raise ValueError("Image must be larger than template.")
125

126
    image_shape = image.shape
127

128
    float_dtype = _supported_float_type(image.dtype)
129
    image = image.astype(float_dtype, copy=False)
130

131
    pad_width = tuple((width, width) for width in template.shape)
132
    if mode == 'constant':
133
        image = np.pad(
134
            image, pad_width=pad_width, mode=mode, constant_values=constant_values
135
        )
136
    else:
137
        image = np.pad(image, pad_width=pad_width, mode=mode)
138

139
    # Use special case for 2-D images for much better performance in
140
    # computation of integral images
141
    if image.ndim == 2:
142
        image_window_sum = _window_sum_2d(image, template.shape)
143
        image_window_sum2 = _window_sum_2d(image**2, template.shape)
144
    elif image.ndim == 3:
145
        image_window_sum = _window_sum_3d(image, template.shape)
146
        image_window_sum2 = _window_sum_3d(image**2, template.shape)
147

148
    template_mean = template.mean()
149
    template_volume = math.prod(template.shape)
150
    template_ssd = np.sum((template - template_mean) ** 2)
151

152
    if image.ndim == 2:
153
        xcorr = fftconvolve(image, template[::-1, ::-1], mode="valid")[1:-1, 1:-1]
154
    elif image.ndim == 3:
155
        xcorr = fftconvolve(image, template[::-1, ::-1, ::-1], mode="valid")[
156
            1:-1, 1:-1, 1:-1
157
        ]
158

159
    numerator = xcorr - image_window_sum * template_mean
160

161
    denominator = image_window_sum2
162
    np.multiply(image_window_sum, image_window_sum, out=image_window_sum)
163
    np.divide(image_window_sum, template_volume, out=image_window_sum)
164
    denominator -= image_window_sum
165
    denominator *= template_ssd
166
    np.maximum(denominator, 0, out=denominator)  # sqrt of negative number not allowed
167
    np.sqrt(denominator, out=denominator)
168

169
    response = np.zeros_like(xcorr, dtype=float_dtype)
170

171
    # avoid zero-division
172
    mask = denominator > np.finfo(float_dtype).eps
173

174
    response[mask] = numerator[mask] / denominator[mask]
175

176
    slices = []
177
    for i in range(template.ndim):
178
        if pad_input:
179
            d0 = (template.shape[i] - 1) // 2
180
            d1 = d0 + image_shape[i]
181
        else:
182
            d0 = template.shape[i] - 1
183
            d1 = d0 + image_shape[i] - template.shape[i] + 1
184
        slices.append(slice(d0, d1))
185

186
    return response[tuple(slices)]
187