template.py 6.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. import math
  2. import numpy as np
  3. from scipy.signal import fftconvolve
  4. from .._shared.utils import check_nD, _supported_float_type
  5. def _window_sum_2d(image, window_shape):
  6. window_sum = np.cumsum(image, axis=0)
  7. window_sum = window_sum[window_shape[0] : -1] - window_sum[: -window_shape[0] - 1]
  8. window_sum = np.cumsum(window_sum, axis=1)
  9. window_sum = (
  10. window_sum[:, window_shape[1] : -1] - window_sum[:, : -window_shape[1] - 1]
  11. )
  12. return window_sum
  13. def _window_sum_3d(image, window_shape):
  14. window_sum = _window_sum_2d(image, window_shape)
  15. window_sum = np.cumsum(window_sum, axis=2)
  16. window_sum = (
  17. window_sum[:, :, window_shape[2] : -1]
  18. - window_sum[:, :, : -window_shape[2] - 1]
  19. )
  20. return window_sum
  21. def match_template(
  22. image, template, pad_input=False, mode='constant', constant_values=0
  23. ):
  24. """Match a template to a 2-D or 3-D image using normalized correlation.
  25. The output is an array with values between -1.0 and 1.0. The value at a
  26. given position corresponds to the correlation coefficient between the image
  27. and the template.
  28. For `pad_input=True` matches correspond to the center and otherwise to the
  29. top-left corner of the template. To find the best match you must search for
  30. peaks in the response (output) image.
  31. Parameters
  32. ----------
  33. image : (M, N[, P]) array
  34. 2-D or 3-D input image.
  35. template : (m, n[, p]) array
  36. Template to locate. It must be `(m <= M, n <= N[, p <= P])`.
  37. pad_input : bool
  38. If True, pad `image` so that output is the same size as the image, and
  39. output values correspond to the template center. Otherwise, the output
  40. is an array with shape `(M - m + 1, N - n + 1)` for an `(M, N)` image
  41. and an `(m, n)` template, and matches correspond to origin
  42. (top-left corner) of the template.
  43. mode : see `numpy.pad`, optional
  44. Padding mode.
  45. constant_values : see `numpy.pad`, optional
  46. Constant values used in conjunction with ``mode='constant'``.
  47. Returns
  48. -------
  49. output : array
  50. Response image with correlation coefficients.
  51. Notes
  52. -----
  53. Details on the cross-correlation are presented in [1]_. This implementation
  54. uses FFT convolutions of the image and the template. Reference [2]_
  55. presents similar derivations but the approximation presented in this
  56. reference is not used in our implementation.
  57. References
  58. ----------
  59. .. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light
  60. and Magic.
  61. .. [2] Briechle and Hanebeck, "Template Matching using Fast Normalized
  62. Cross Correlation", Proceedings of the SPIE (2001).
  63. :DOI:`10.1117/12.421129`
  64. Examples
  65. --------
  66. >>> template = np.zeros((3, 3))
  67. >>> template[1, 1] = 1
  68. >>> template
  69. array([[0., 0., 0.],
  70. [0., 1., 0.],
  71. [0., 0., 0.]])
  72. >>> image = np.zeros((6, 6))
  73. >>> image[1, 1] = 1
  74. >>> image[4, 4] = -1
  75. >>> image
  76. array([[ 0., 0., 0., 0., 0., 0.],
  77. [ 0., 1., 0., 0., 0., 0.],
  78. [ 0., 0., 0., 0., 0., 0.],
  79. [ 0., 0., 0., 0., 0., 0.],
  80. [ 0., 0., 0., 0., -1., 0.],
  81. [ 0., 0., 0., 0., 0., 0.]])
  82. >>> result = match_template(image, template)
  83. >>> np.round(result, 3)
  84. array([[ 1. , -0.125, 0. , 0. ],
  85. [-0.125, -0.125, 0. , 0. ],
  86. [ 0. , 0. , 0.125, 0.125],
  87. [ 0. , 0. , 0.125, -1. ]])
  88. >>> result = match_template(image, template, pad_input=True)
  89. >>> np.round(result, 3)
  90. array([[-0.125, -0.125, -0.125, 0. , 0. , 0. ],
  91. [-0.125, 1. , -0.125, 0. , 0. , 0. ],
  92. [-0.125, -0.125, -0.125, 0. , 0. , 0. ],
  93. [ 0. , 0. , 0. , 0.125, 0.125, 0.125],
  94. [ 0. , 0. , 0. , 0.125, -1. , 0.125],
  95. [ 0. , 0. , 0. , 0.125, 0.125, 0.125]])
  96. """
  97. check_nD(image, (2, 3))
  98. if image.ndim < template.ndim:
  99. raise ValueError(
  100. "Dimensionality of template must be less than or "
  101. "equal to the dimensionality of image."
  102. )
  103. if np.any(np.less(image.shape, template.shape)):
  104. raise ValueError("Image must be larger than template.")
  105. image_shape = image.shape
  106. float_dtype = _supported_float_type(image.dtype)
  107. image = image.astype(float_dtype, copy=False)
  108. pad_width = tuple((width, width) for width in template.shape)
  109. if mode == 'constant':
  110. image = np.pad(
  111. image, pad_width=pad_width, mode=mode, constant_values=constant_values
  112. )
  113. else:
  114. image = np.pad(image, pad_width=pad_width, mode=mode)
  115. # Use special case for 2-D images for much better performance in
  116. # computation of integral images
  117. if image.ndim == 2:
  118. image_window_sum = _window_sum_2d(image, template.shape)
  119. image_window_sum2 = _window_sum_2d(image**2, template.shape)
  120. elif image.ndim == 3:
  121. image_window_sum = _window_sum_3d(image, template.shape)
  122. image_window_sum2 = _window_sum_3d(image**2, template.shape)
  123. template_mean = template.mean()
  124. template_volume = math.prod(template.shape)
  125. template_ssd = np.sum((template - template_mean) ** 2)
  126. if image.ndim == 2:
  127. xcorr = fftconvolve(image, template[::-1, ::-1], mode="valid")[1:-1, 1:-1]
  128. elif image.ndim == 3:
  129. xcorr = fftconvolve(image, template[::-1, ::-1, ::-1], mode="valid")[
  130. 1:-1, 1:-1, 1:-1
  131. ]
  132. numerator = xcorr - image_window_sum * template_mean
  133. denominator = image_window_sum2
  134. np.multiply(image_window_sum, image_window_sum, out=image_window_sum)
  135. np.divide(image_window_sum, template_volume, out=image_window_sum)
  136. denominator -= image_window_sum
  137. denominator *= template_ssd
  138. np.maximum(denominator, 0, out=denominator) # sqrt of negative number not allowed
  139. np.sqrt(denominator, out=denominator)
  140. response = np.zeros_like(xcorr, dtype=float_dtype)
  141. # avoid zero-division
  142. mask = denominator > np.finfo(float_dtype).eps
  143. response[mask] = numerator[mask] / denominator[mask]
  144. slices = []
  145. for i in range(template.ndim):
  146. if pad_input:
  147. d0 = (template.shape[i] - 1) // 2
  148. d1 = d0 + image_shape[i]
  149. else:
  150. d0 = template.shape[i] - 1
  151. d1 = d0 + image_shape[i] - template.shape[i] + 1
  152. slices.append(slice(d0, d1))
  153. return response[tuple(slices)]