mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-08 11:26:12 +08:00
Convert template matching Pull Request to skimage
This commit takes the template matching implementation from holtzhau/template (Pull Request #13) and converts the code to use the new package name (scikits.image --> skimage).
This commit is contained in:
Tony S Yu
@@ -0,0 +1,70 @@
|
||||
"""
|
||||
=================
|
||||
Template Matching
|
||||
=================
|
||||
|
||||
In this example, we use template matching to identify the occurrence of an
|
||||
object in an image. The ``match_template`` function uses normalised correlation
|
||||
techniques to find instances of the "target image" in the "test image".
|
||||
|
||||
The output of ``match_template`` is an image where we can easily identify peaks
|
||||
by eye. Nevertheless, this example concludes with a simple peak extraction
|
||||
algorithm to quantify the locations of matches.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from skimage.detection import match_template
|
||||
from numpy.random import randn
|
||||
import matplotlib.pyplot as plt
|
||||
import math
|
||||
|
||||
# We first construct a simple image target:
|
||||
size = 100
|
||||
target = np.tri(size) + np.tri(size)[::-1]
|
||||
target = target.astype(np.float32)
|
||||
|
||||
plt.gray()
|
||||
plt.imshow(target)
|
||||
plt.title("Target image")
|
||||
plt.axis('off')
|
||||
|
||||
# place target in an image at two positions, and add noise.
|
||||
image = np.zeros((400, 400), dtype=np.float32)
|
||||
target_positions = [(50, 50), (200, 200)]
|
||||
for x, y in target_positions:
|
||||
image[x:x+size, y:y+size] = target
|
||||
image += randn(400, 400)*2
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(image)
|
||||
plt.title("Test image")
|
||||
plt.axis('off')
|
||||
|
||||
# Match the template.
|
||||
result = match_template(image, target, method='norm-corr')
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(result)
|
||||
plt.title("Result from ``match_template``")
|
||||
plt.axis('off')
|
||||
|
||||
plt.show()
|
||||
|
||||
# peak extraction algorithm.
|
||||
delta = 5
|
||||
found_positions = []
|
||||
for i in range(50):
|
||||
index = np.argmax(result)
|
||||
y, x = np.unravel_index(index, result.shape)
|
||||
if not found_positions:
|
||||
found_positions.append((x, y))
|
||||
for position in found_positions:
|
||||
distance = math.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2)
|
||||
if distance > delta:
|
||||
found_positions.append((x, y))
|
||||
result[y, x] = 0
|
||||
if len(found_positions) == len(target_positions):
|
||||
break
|
||||
|
||||
assert np.all(found_positions == target_positions)
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
from template import match_template
|
||||
@@ -0,0 +1,219 @@
|
||||
"""template.py - Template matching
|
||||
"""
|
||||
import cython
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
import cv
|
||||
from scipy.signal import fftconvolve
|
||||
|
||||
cdef extern from "math.h":
|
||||
double sqrt(double x)
|
||||
double fabs(double x)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image):
|
||||
"""
|
||||
Calculate the summed integral image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Source image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray, dtype=np.double_t
|
||||
Summed integral image.
|
||||
"""
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1]))
|
||||
cdef double s
|
||||
cdef int x, y
|
||||
cdef int width, height
|
||||
height = image.shape[0]
|
||||
width = image.shape[1]
|
||||
ii[0, 0] = image[0, 0]
|
||||
|
||||
for y in range(1, height):
|
||||
ii[y, 0] = image[y, 0] + ii[y - 1, 0]
|
||||
|
||||
for x in range(1, width):
|
||||
s = 0
|
||||
for y in range(0, height):
|
||||
s += image[y, x]
|
||||
ii[y, x] = s + ii[y, x - 1]
|
||||
|
||||
return ii
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image):
|
||||
"""
|
||||
Calculate the squared integral image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Source image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray, dtype=np.double_t
|
||||
Squared integral image.
|
||||
"""
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1]))
|
||||
cdef double s
|
||||
cdef int x, y
|
||||
cdef int width, height
|
||||
height = image.shape[0]
|
||||
width = image.shape[1]
|
||||
ii2[0, 0] = image[0, 0] * image[0, 0]
|
||||
|
||||
for y in range(1, height):
|
||||
ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0]
|
||||
|
||||
for x in range(1, width):
|
||||
s = 0
|
||||
for y in range(0, height):
|
||||
s += image[y, x] * image[y, x]
|
||||
ii2[y, x] = s + ii2[y, x - 1]
|
||||
|
||||
return ii2
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image):
|
||||
"""
|
||||
Calculate the summed and sqared integral image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Source image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : tuple (ndarray, ndarray) of type np.double_t
|
||||
Summed and squared integral image.
|
||||
"""
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1]))
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1]))
|
||||
cdef double s, s2
|
||||
cdef int x, y
|
||||
cdef int width, height
|
||||
height = image.shape[0]
|
||||
width = image.shape[1]
|
||||
ii[0, 0] = image[0, 0]
|
||||
ii2[0, 0] = image[0, 0] * image[0, 0]
|
||||
|
||||
for y in range(1, height):
|
||||
ii[y, 0] = image[y, 0] + ii[y - 1, 0]
|
||||
ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0]
|
||||
|
||||
for x in range(1, width):
|
||||
s = 0
|
||||
s2 = 0
|
||||
for y in range(0, height):
|
||||
s += image[y, x]
|
||||
s2 += image[y, x] * image[y, x]
|
||||
ii[y, x] = s + ii[y, x - 1]
|
||||
ii2[y, x] = s2 + ii2[y, x - 1]
|
||||
|
||||
return ii, ii2
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef double sum_integral(np.ndarray[np.double_t, 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.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
sat : ndarray of double_t
|
||||
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 double 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, int num_type):
|
||||
# convolve the image with template by frequency domain multiplication
|
||||
cdef np.ndarray[np.double_t, ndim=2] result
|
||||
result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), mode="valid"), dtype=np.double)
|
||||
# calculate squared integral images used for normalization
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sum
|
||||
cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sqr
|
||||
if num_type == 1:
|
||||
integral_sum, integral_sqr = integral_images(image)
|
||||
else:
|
||||
integral_sqr = integral_image_sqr(image)
|
||||
|
||||
# use inversed area for accuracy
|
||||
cdef double inv_area = 1.0 / (template.shape[0] * template.shape[1])
|
||||
# calculate template norm according to the following:
|
||||
# variance ** 2 = 1/K Sigma[(x_k - mean) ** 2] = 1/K Sigma[x_k ** 2] - mean ** 2
|
||||
cdef double template_norm
|
||||
cdef double template_mean = np.mean(template)
|
||||
|
||||
if num_type == 0:
|
||||
template_norm = sqrt((np.std(template) ** 2 + template_mean ** 2)) / sqrt(inv_area)
|
||||
else:
|
||||
template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area)
|
||||
|
||||
# define window of template size in squared integral image
|
||||
cdef int i, j
|
||||
cdef double num, window_sum2, window_mean2, normed, t,
|
||||
# move window through convolution results, normalizing in the process
|
||||
for i in range(result.shape[0] - 1):
|
||||
for j in range(result.shape[1] - 1):
|
||||
num = result[i, j]
|
||||
window_mean2 = 0
|
||||
if num_type == 1:
|
||||
t = sum_integral(integral_sum, i, j, i + template.shape[0], j + template.shape[1])
|
||||
window_mean2 = t * t * inv_area
|
||||
num -= t*template_mean
|
||||
|
||||
# calculate squared template window sum in the image
|
||||
window_sum2 = sum_integral(integral_sqr, i, j, i + template.shape[0], j + template.shape[1])
|
||||
normed = sqrt(window_sum2 - window_mean2) * template_norm
|
||||
# enforce some limits
|
||||
if fabs(num) < normed:
|
||||
num /= normed
|
||||
elif fabs(num) < normed*1.125:
|
||||
if num > 0:
|
||||
num = 1
|
||||
else:
|
||||
num = -1
|
||||
else:
|
||||
num = 0
|
||||
result[i, j] = num
|
||||
# zero boundaries
|
||||
for i in range(result.shape[0]):
|
||||
result[i, -1] = 0
|
||||
for j in range(result.shape[1]):
|
||||
result[-1, j] = 0
|
||||
return result
|
||||
@@ -0,0 +1,33 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
import shutil
|
||||
import hashlib
|
||||
|
||||
from skimage._build import cython
|
||||
|
||||
base_path = os.path.abspath(os.path.dirname(__file__))
|
||||
|
||||
def configuration(parent_package='', top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
|
||||
|
||||
config = Configuration('transform', parent_package, top_path)
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_template.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_template', sources=['_template.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(maintainer = 'Scikits.Image Developers',
|
||||
author = 'Scikits.Image Developers',
|
||||
maintainer_email = 'scikits-image@googlegroups.com',
|
||||
description = 'Transforms',
|
||||
url = 'http://stefanv.github.com/scikits.image/',
|
||||
license = 'SciPy License (BSD Style)',
|
||||
**(configuration(top_path='').todict())
|
||||
)
|
||||
@@ -0,0 +1,75 @@
|
||||
"""template.py - Template matching
|
||||
"""
|
||||
import numpy as np
|
||||
import cv
|
||||
import _template
|
||||
|
||||
#XXX add to opencv backend once backend system in place
|
||||
def match_template_cv(image, template, out=None, method="norm-coeff"):
|
||||
"""Finds a template in an image using normalized correlation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Image to process.
|
||||
template : array_like, dtype=float
|
||||
Template to locate.
|
||||
out: array_like, dtype=float, optional
|
||||
Optional destination.
|
||||
Returns
|
||||
-------
|
||||
output : ndarray, dtype=float
|
||||
Correlation results between 0.0 and 1.0, maximum indicating the most probable match.
|
||||
"""
|
||||
if out == None:
|
||||
out = np.empty((image.shape[0] - template.shape[0] + 1,image.shape[1] - template.shape[1] + 1), dtype=image.dtype)
|
||||
if method == "norm-corr":
|
||||
cv.MatchTemplate(image, template, out, cv.CV_TM_CCORR_NORMED)
|
||||
elif method == "norm-corr":
|
||||
cv.MatchTemplate(image, template, out, cv.CV_TM_CCOEFF_NORMED)
|
||||
else:
|
||||
raise ValueError("Unknown template method: %s" % method)
|
||||
return out
|
||||
|
||||
|
||||
def match_template(image, template, method="norm-coeff"):
|
||||
"""Finds a template in an image using normalized correlation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Image to process.
|
||||
template : array_like, dtype=float
|
||||
Template to locate.
|
||||
method: str (default 'norm-coeff')
|
||||
The correlation method used in scanning.
|
||||
T represents the template, I the image and R the result.
|
||||
The summation is done over x' = 0..w-1 and y' = 0..h-1 of the template.
|
||||
'norm-coeff':
|
||||
R(x, y) = Sigma(x',y')[T(x', y').I(x + x', y + y')] / N
|
||||
N = sqrt(Sigma(x',y')[T(x', y')**2].Sigma(x',y')[I(x + x', y + y')**2])
|
||||
'norm-corr':
|
||||
R(x,y) = Sigma(x',y)[T'(x', y').I'(x + x', y + y')] / N
|
||||
N = sqrt(Sigma(x',y)[T'(x', y')**2].Sigma(x',y')[I'(x + x', y + y')**2])
|
||||
where:
|
||||
T'(x, y) = T(x', y') - 1/(w.h).Sigma(x'',y'')[T(x'', y'')]
|
||||
I'(x + x', y + y') = I(x + x', y + y') -
|
||||
1/(w.h).Sigma(x'',y'')[I(x + x'', y + y'')]
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray, dtype=float
|
||||
Correlation results between 0.0 and 1.0, maximum indicating the most
|
||||
probable match.
|
||||
"""
|
||||
if method == "norm-corr":
|
||||
method_num = 0
|
||||
elif method == "norm-coeff":
|
||||
method_num = 1
|
||||
else:
|
||||
raise ValueError("Unknown template method: %s" % method)
|
||||
return _template.match_template(image, template, method_num)
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,46 @@
|
||||
import os.path
|
||||
import numpy as np
|
||||
from numpy.testing import *
|
||||
from skimage import data_dir
|
||||
from skimage.detection import *
|
||||
from numpy.random import randn
|
||||
|
||||
def test_template():
|
||||
size = 100
|
||||
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
|
||||
image += randn(400, 400)*2
|
||||
|
||||
for method in ["norm-corr", "norm-coeff"]:
|
||||
result = match_template(image, target, method=method)
|
||||
delta = 5
|
||||
found_positions = []
|
||||
# find the targets
|
||||
for i in range(50):
|
||||
index = np.argmax(result)
|
||||
y, x = np.unravel_index(index, result.shape)
|
||||
if not found_positions:
|
||||
found_positions.append((x, y))
|
||||
for position in found_positions:
|
||||
distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2)
|
||||
if distance > delta:
|
||||
found_positions.append((x, y))
|
||||
result[y, x] = 0
|
||||
if len(found_positions) == len(target_positions):
|
||||
break
|
||||
|
||||
for x, y in target_positions:
|
||||
print x, y
|
||||
found = False
|
||||
for position in found_positions:
|
||||
distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2)
|
||||
if distance < delta:
|
||||
found = True
|
||||
assert found
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
@@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_subpackage('draw')
|
||||
config.add_subpackage('feature')
|
||||
config.add_subpackage('measure')
|
||||
config.add_subpackage('detection')
|
||||
|
||||
def add_test_directories(arg, dirname, fnames):
|
||||
if dirname.split(os.path.sep)[-1] == 'tests':
|
||||
|
||||
Reference in New Issue
Block a user