From 9fa6000344633c1728feae31222e8083aadf2259 Mon Sep 17 00:00:00 2001
From: Andre Anjos <andre.anjos@idiap.ch>
Date: Tue, 8 Nov 2016 15:15:19 +0100
Subject: [PATCH] Improved documentation on MiuraMatch algorithm; Use scipy
convolution
---
bob/bio/vein/algorithm/MiuraMatch.py | 89 +++++++++---------
bob/bio/vein/tests/test.py | 136 +++++++++++++++++++--------
2 files changed, 142 insertions(+), 83 deletions(-)
diff --git a/bob/bio/vein/algorithm/MiuraMatch.py b/bob/bio/vein/algorithm/MiuraMatch.py
index 56b995c..86e0263 100644
--- a/bob/bio/vein/algorithm/MiuraMatch.py
+++ b/bob/bio/vein/algorithm/MiuraMatch.py
@@ -1,32 +1,46 @@
#!/usr/bin/env python
# vim: set fileencoding=utf-8 :
-import bob.sp
-import bob.ip.base
-
import numpy
-import math
import scipy.signal
+import bob.ip.base
from bob.bio.base.algorithm import Algorithm
class MiuraMatch (Algorithm):
- """Finger vein matching: match ratio
+ """Finger vein matching: match ratio via cross-correlation
+
+ The method is based on "cross-correlation" between a model and a probe image.
+ It convolves the binary image(s) representing the model with the binary image
+ representing the probe (rotated by 180 degrees), to evaluate how they
+ cross-correlate. If the model and probe are very similar, the output of the
+ correlation corresponds to a single scalar and approaches a maximum. The
+ value is then normalized by the sum of the pixels lit in both binary images.
+ Therefore, the output of this method is a floating-point number in the range
+ :math:`[0, 0.5]`. The higher, the better match.
+
+ In case model and probe represent images from the same vein structure, but
+ are misaligned, the output is not guaranteed to be accurate. To mitigate this
+ aspect, Miura et al. proposed to add a *small** erosion factor to the model
+ image, assuming not much information is available on the borders (``ch``, for
+ the vertical direction and ``cw``, for the horizontal direction). This allows
+ the convolution to yield searches for different areas in the probe image. The
+ maximum value is then taken from the resulting operation. The convolution
+ result is normalized by the pixels lit in both the eroded model image and the
+ matching pixels on the probe that yield the maximum on the resulting
+ convolution.
Based on N. Miura, A. Nagasaka, and T. Miyatake. Feature extraction of finger
vein patterns based on repeated line tracking and its application to personal
identification. Machine Vision and Applications, Vol. 15, Num. 4, pp.
194--203, 2004
-
Parameters:
- ch (:py:class:`int`, optional): Maximum search displacement in y-direction.
- Different default values based on the different features.
+ ch (:py:class:`int`, optional): Maximum search displacement in y-direction.
- cw (:py:class:`int`, optional): Maximum search displacement in x-direction.
- Different default values based on the different features.
+ cw (:py:class:`int`, optional): Maximum search displacement in x-direction.
"""
@@ -57,39 +71,21 @@ class MiuraMatch (Algorithm):
return numpy.array(enroll_features)
- def convfft(self, t, a):
- # Determine padding size in x and y dimension
- size_t = numpy.array(t.shape)
- size_a = numpy.array(a.shape)
- outsize = size_t + size_a - 1
-
- # Determine 2D cross correlation in Fourier domain
- taux = numpy.zeros(outsize)
- taux[0:size_t[0],0:size_t[1]] = t
- Ft = bob.sp.fft(taux.astype(numpy.complex128))
- aaux = numpy.zeros(outsize)
- aaux[0:size_a[0],0:size_a[1]] = a
- Fa = bob.sp.fft(aaux.astype(numpy.complex128))
+ def score(self, model, probe):
+ """Computes the score between the probe and the model.
- convta = numpy.real(bob.sp.ifft(Ft*Fa))
+ Parameters:
- [w, h] = size_t-size_a+1
- output = convta[size_a[0]-1:size_a[0]-1+w, size_a[1]-1:size_a[1]-1+h]
+ model (numpy.ndarray): The model of the user to test the probe agains
- return output
+ probe (numpy.ndarray): The probe to test
- def score(self, model, probe):
- """
- Computes the score of the probe and the model.
+ Returns:
- **Parameters:**
+ score (float): Value between 0 and 0.5, larger value means a better match
- score : :py:class:`float`
- Value between 0 and 0.5, larger value is better match
"""
- #print model.shape
- #print probe.shape
I=probe.astype(numpy.float64)
@@ -99,22 +95,31 @@ class MiuraMatch (Algorithm):
n_models = model.shape[0]
scores = []
+
for i in range(n_models):
+
+ # erode model by (ch, cw)
R=model[i,:].astype(numpy.float64)
h, w = R.shape
crop_R = R[self.ch:h-self.ch, self.cw:w-self.cw]
+
+ # rotate input image
rotate_R = numpy.zeros((crop_R.shape[0], crop_R.shape[1]))
bob.ip.base.rotate(crop_R, rotate_R, 180)
- #FFT for scoring!
- #Nm=bob.sp.ifft(bob.sp.fft(I)*bob.sp.fft(rotate_R))
- Nm = self.convfft(I, rotate_R)
- #Nm2 = scipy.signal.convolve2d(I, rotate_R, 'valid')
+ # convolve model and probe using FFT/IFFT.
+ #Nm = utils.convfft(I, rotate_R) #drop-in replacement for scipy method
+ Nm = scipy.signal.convolve2d(I, rotate_R, 'valid')
+
+ # figures out where the maximum is on the resulting matrix
t0, s0 = numpy.unravel_index(Nm.argmax(), Nm.shape)
+
+ # this is our output
Nmm = Nm[t0,s0]
- #Nmm = Nm.max()
- #mi = numpy.argwhere(Nmm == Nm)
- #t0, s0 = mi.flatten()[:2]
+
+ # normalizes the output by the number of pixels lit on the input
+ # matrices, taking into consideration the surface that produced the
+ # result (i.e., the eroded model and part of the probe)
scores.append(Nmm/(sum(sum(crop_R)) + sum(sum(I[t0:t0+h-2*self.ch, s0:s0+w-2*self.cw]))))
return numpy.mean(scores)
diff --git a/bob/bio/vein/tests/test.py b/bob/bio/vein/tests/test.py
index da4715c..b2c9530 100644
--- a/bob/bio/vein/tests/test.py
+++ b/bob/bio/vein/tests/test.py
@@ -23,7 +23,8 @@ import bob.io.base
import bob.io.matlab
import bob.io.image
-from ..preprocessor import utils
+from ..preprocessor import utils as preprocessor_utils
+from ..algorithm import utils as algorithm_utils
def F(parts):
@@ -46,7 +47,7 @@ def test_finger_crop():
preprocess = FingerCrop(fingercontour='leemaskMatlab', padding_width=0)
preproc, mask = preprocess(img)
- #utils.show_mask_over_image(preproc, mask)
+ #preprocessor_utils.show_mask_over_image(preproc, mask)
mask_ref = bob.io.base.load(output_fvr_filename).astype('bool')
preproc_ref = bob.core.convert(bob.io.base.load(output_img_filename),
@@ -54,8 +55,8 @@ def test_finger_crop():
assert numpy.mean(numpy.abs(mask - mask_ref)) < 1e-2
- # Very loose comparison!
- #utils.show_image(numpy.abs(preproc.astype('int16') - preproc_ref.astype('int16')).astype('uint8'))
+ # Very loose comparison!
+ #preprocessor_utils.show_image(numpy.abs(preproc.astype('int16') - preproc_ref.astype('int16')).astype('uint8'))
assert numpy.mean(numpy.abs(preproc - preproc_ref)) < 1.3e2
@@ -160,12 +161,12 @@ def test_assert_points():
# Tests that point assertion works as expected
area = (10, 5)
inside = [(0,0), (3,2), (9, 4)]
- utils.assert_points(area, inside) #should not raise
+ preprocessor_utils.assert_points(area, inside) #should not raise
def _check_outside(point):
# should raise, otherwise it is an error
try:
- utils.assert_points(area, [point])
+ preprocessor_utils.assert_points(area, [point])
except AssertionError as e:
assert str(point) in str(e)
else:
@@ -180,22 +181,22 @@ def test_fix_points():
# Tests that point clipping works as expected
area = (10, 5)
inside = [(0,0), (3,2), (9, 4)]
- fixed = utils.fix_points(area, inside)
+ fixed = preprocessor_utils.fix_points(area, inside)
assert numpy.array_equal(inside, fixed), '%r != %r' % (inside, fixed)
- fixed = utils.fix_points(area, [(-1, 0)])
+ fixed = preprocessor_utils.fix_points(area, [(-1, 0)])
assert numpy.array_equal(fixed, [(0, 0)])
- fixed = utils.fix_points(area, [(10, 0)])
+ fixed = preprocessor_utils.fix_points(area, [(10, 0)])
assert numpy.array_equal(fixed, [(9, 0)])
- fixed = utils.fix_points(area, [(0, 5)])
+ fixed = preprocessor_utils.fix_points(area, [(0, 5)])
assert numpy.array_equal(fixed, [(0, 4)])
- fixed = utils.fix_points(area, [(10, 5)])
+ fixed = preprocessor_utils.fix_points(area, [(10, 5)])
assert numpy.array_equal(fixed, [(9, 4)])
- fixed = utils.fix_points(area, [(15, 12)])
+ fixed = preprocessor_utils.fix_points(area, [(15, 12)])
assert numpy.array_equal(fixed, [(9, 4)])
@@ -204,7 +205,7 @@ def test_poly_to_mask():
# Tests we can generate a mask out of a polygon correctly
area = (10, 9) #10 rows, 9 columns
polygon = [(2, 2), (2, 7), (7, 7), (7, 2)] #square shape, (y, x) format
- mask = utils.poly_to_mask(area, polygon)
+ mask = preprocessor_utils.poly_to_mask(area, polygon)
nose.tools.eq_(mask.dtype, numpy.bool)
# This should be the output:
@@ -223,7 +224,7 @@ def test_poly_to_mask():
assert numpy.array_equal(mask, expected)
polygon = [(3, 2), (5, 7), (8, 7), (7, 3)] #trapezoid, (y, x) format
- mask = utils.poly_to_mask(area, polygon)
+ mask = preprocessor_utils.poly_to_mask(area, polygon)
nose.tools.eq_(mask.dtype, numpy.bool)
# This should be the output:
@@ -250,7 +251,7 @@ def test_mask_to_image():
nose.tools.eq_(sample.dtype, numpy.bool)
def _check_uint(n):
- conv = utils.mask_to_image(sample, 'uint%d' % n)
+ conv = preprocessor_utils.mask_to_image(sample, 'uint%d' % n)
nose.tools.eq_(conv.dtype, getattr(numpy, 'uint%d' % n))
target = [0, (2**n)-1]
assert numpy.array_equal(conv, target), '%r != %r' % (conv, target)
@@ -261,7 +262,7 @@ def test_mask_to_image():
_check_uint(64)
def _check_float(n):
- conv = utils.mask_to_image(sample, 'float%d' % n)
+ conv = preprocessor_utils.mask_to_image(sample, 'float%d' % n)
nose.tools.eq_(conv.dtype, getattr(numpy, 'float%d' % n))
assert numpy.array_equal(conv, [0, 1.0]), '%r != %r' % (conv, target)
@@ -272,7 +273,7 @@ def test_mask_to_image():
# This should be unsupported
try:
- conv = utils.mask_to_image(sample, 'int16')
+ conv = preprocessor_utils.mask_to_image(sample, 'int16')
except TypeError as e:
assert 'int16' in str(e)
else:
@@ -292,15 +293,13 @@ def test_jaccard_index():
[True, False],
])
- nose.tools.eq_(utils.jaccard_index(a, b), 1.0/4.0)
- nose.tools.eq_(utils.jaccard_index(a, a), 1.0)
- nose.tools.eq_(utils.jaccard_index(b, b), 1.0)
- nose.tools.eq_(utils.jaccard_index(a, numpy.ones(a.shape, dtype=bool)),
- 2.0/4.0)
- nose.tools.eq_(utils.jaccard_index(a, numpy.zeros(a.shape, dtype=bool)), 0.0)
- nose.tools.eq_(utils.jaccard_index(b, numpy.ones(b.shape, dtype=bool)),
- 3.0/4.0)
- nose.tools.eq_(utils.jaccard_index(b, numpy.zeros(b.shape, dtype=bool)), 0.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(a, b), 1.0/4.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(a, a), 1.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(b, b), 1.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(a, numpy.ones(a.shape, dtype=bool)), 2.0/4.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(a, numpy.zeros(a.shape, dtype=bool)), 0.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(b, numpy.ones(b.shape, dtype=bool)), 3.0/4.0)
+ nose.tools.eq_(preprocessor_utils.jaccard_index(b, numpy.zeros(b.shape, dtype=bool)), 0.0)
def test_intersection_ratio():
@@ -316,18 +315,73 @@ def test_intersection_ratio():
[True, False],
])
- nose.tools.eq_(utils.intersect_ratio(a, b), 1.0/2.0)
- nose.tools.eq_(utils.intersect_ratio(a, a), 1.0)
- nose.tools.eq_(utils.intersect_ratio(b, b), 1.0)
- nose.tools.eq_(utils.intersect_ratio(a, numpy.ones(a.shape, dtype=bool)), 1.0)
- nose.tools.eq_(utils.intersect_ratio(a, numpy.zeros(a.shape, dtype=bool)), 0)
- nose.tools.eq_(utils.intersect_ratio(b, numpy.ones(b.shape, dtype=bool)), 1.0)
- nose.tools.eq_(utils.intersect_ratio(b, numpy.zeros(b.shape, dtype=bool)), 0)
-
- nose.tools.eq_(utils.intersect_ratio_of_complement(a, b), 1.0/2.0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(a, a), 0.0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(b, b), 0.0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(a, numpy.ones(a.shape, dtype=bool)), 1.0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(a, numpy.zeros(a.shape, dtype=bool)), 0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(b, numpy.ones(b.shape, dtype=bool)), 1.0)
- nose.tools.eq_(utils.intersect_ratio_of_complement(b, numpy.zeros(b.shape, dtype=bool)), 0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(a, b), 1.0/2.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(a, a), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(b, b), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(a, numpy.ones(a.shape, dtype=bool)), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(a, numpy.zeros(a.shape, dtype=bool)), 0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(b, numpy.ones(b.shape, dtype=bool)), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio(b, numpy.zeros(b.shape, dtype=bool)), 0)
+
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(a, b), 1.0/2.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(a, a), 0.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(b, b), 0.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(a, numpy.ones(a.shape, dtype=bool)), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(a, numpy.zeros(a.shape, dtype=bool)), 0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(b, numpy.ones(b.shape, dtype=bool)), 1.0)
+ nose.tools.eq_(preprocessor_utils.intersect_ratio_of_complement(b, numpy.zeros(b.shape, dtype=bool)), 0)
+
+
+def test_convolution():
+
+ # A test for convolution performance. Convolutions are used on the Miura
+ # Match algorithm, therefore we want to make sure we can perform them as fast
+ # as possible.
+ import scipy.signal
+
+ Y = 250
+ X = 600
+ CH = 1
+ CW = 1
+
+ def gen_ab():
+ a = numpy.random.randint(256, size=(Y, X)).astype(float)
+ b = numpy.random.randint(256, size=(Y-CH, X-CW)).astype(float)
+ return a, b
+
+ def utils_function(a, b):
+ return algorithm_utils.convfft(a, b)
+
+ def scipy_function(a, b):
+ return scipy.signal.convolve2d(a, b, 'valid')
+
+ def scipy2_function(a, b):
+ return scipy.signal.fftconvolve(a, b, 'valid')
+
+ a, b = gen_ab()
+ assert numpy.allclose(utils_function(a, b), scipy_function(a, b))
+ assert numpy.allclose(scipy_function(a, b), scipy2_function(a, b))
+
+ import time
+
+ start = time.clock()
+ N = 10
+ for i in range(N):
+ a, b = gen_ab()
+ utils_function(a, b)
+ total = time.clock() - start
+ print('utils, %d iterations - %.2e per iteration' % (N, total/N))
+
+ start = time.clock()
+ for i in range(N):
+ a, b = gen_ab()
+ scipy_function(a, b)
+ total = time.clock() - start
+ print('scipy, %d iterations - %.2e per iteration' % (N, total/N))
+
+ start = time.clock()
+ for i in range(N):
+ a, b = gen_ab()
+ scipy2_function(a, b)
+ total = time.clock() - start
+ print('scipy2, %d iterations - %.2e per iteration' % (N, total/N))
--
GitLab