Coverage for src/bob/bio/vein/extractor/MaximumCurvature.py: 79%
117 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-07-12 23:27 +0200
1#!/usr/bin/env python
2# vim: set fileencoding=utf-8 :
4import math
6import numpy
7import scipy.ndimage
9from bob.bio.base.extractor import Extractor
12class MaximumCurvature(Extractor):
13 """
14 MiuraMax feature extractor.
16 Based on N. Miura, A. Nagasaka, and T. Miyatake, Extraction of Finger-Vein
17 Pattern Using Maximum Curvature Points in Image Profiles. Proceedings on IAPR
18 conference on machine vision applications, 9 (2005), pp. 347--350.
21 Parameters:
23 sigma (:py:class:`int`, optional): standard deviation for the gaussian
24 smoothing kernel used to denoise the input image. The width of the
25 gaussian kernel will be set automatically to 4x this value (in pixels).
27 """
29 def __init__(self, sigma=5):
30 Extractor.__init__(self, sigma=sigma)
31 self.sigma = sigma
33 def detect_valleys(self, image, mask):
34 r"""Detects valleys on the image respecting the mask
36 This step corresponds to Step 1-1 in the original paper. The objective is,
37 for all 4 cross-sections (z) of the image (horizontal, vertical, 45 and -45
38 diagonals), to compute the following proposed valley detector as defined in
39 Equation 1, page 348:
41 .. math::
43 \kappa(z) = \\frac{d^2P_f(z)/dz^2}{(1 + (dP_f(z)/dz)^2)^\\frac{3}{2}}
46 We start the algorithm by smoothing the image with a 2-dimensional gaussian
47 filter. The equation that defines the kernel for the filter is:
49 .. math::
51 \mathcal{N}(x,y)=\\frac{1}{2\pi\sigma^2}e^\\frac{-(x^2+y^2)}{2\sigma^2}
54 This is done to avoid noise from the raw data (from the sensor). The
55 maximum curvature method then requires we compute the first and second
56 derivative of the image for all cross-sections, as per the equation above.
58 We instead take the following equivalent approach:
60 1. construct a gaussian filter
61 2. take the first (dh/dx) and second (d^2/dh^2) deritivatives of the filter
62 3. calculate the first and second derivatives of the smoothed signal using
63 the results from 3. This is done for all directions we're interested in:
64 horizontal, vertical and 2 diagonals. First and second derivatives of a
65 convolved signal
67 .. note::
69 Item 3 above is only possible thanks to the steerable filter property of
70 the gaussian kernel. See "The Design and Use of Steerable Filters" from
71 Freeman and Adelson, IEEE Transactions on Pattern Analysis and Machine
72 Intelligence, Vol. 13, No. 9, September 1991.
75 Parameters:
77 image (numpy.ndarray): an array of 64-bit floats containing the input
78 image
79 mask (numpy.ndarray): an array, of the same size as ``image``, containing
80 a mask (booleans) indicating where the finger is on ``image``.
83 Returns:
85 numpy.ndarray: a 3-dimensional array of 64-bits containing $\kappa$ for
86 all considered directions. $\kappa$ has the same shape as ``image``,
87 except for the 3rd. dimension, which provides planes for the
88 cross-section valley detections for each of the contemplated directions,
89 in this order: horizontal, vertical, +45 degrees, -45 degrees.
91 """
93 # 1. constructs the 2D gaussian filter "h" given the window size,
94 # extrapolated from the "sigma" parameter (4x)
95 # N.B.: This is a text-book gaussian filter definition
96 winsize = numpy.ceil(4 * self.sigma) # enough space for the filter
97 window = numpy.arange(-winsize, winsize + 1)
98 X, Y = numpy.meshgrid(window, window)
99 G = 1.0 / (2 * math.pi * self.sigma**2)
100 G *= numpy.exp(-(X**2 + Y**2) / (2 * self.sigma**2))
102 # 2. calculates first and second derivatives of "G" with respect to "X"
103 # (0), "Y" (90 degrees) and 45 degrees (?)
104 G1_0 = (-X / (self.sigma**2)) * G
105 G2_0 = ((X**2 - self.sigma**2) / (self.sigma**4)) * G
106 G1_90 = G1_0.T
107 G2_90 = G2_0.T
108 hxy = ((X * Y) / (self.sigma**4)) * G
110 # 3. calculates derivatives w.r.t. to all directions of interest
111 # stores results in the variable "k". The entries (last dimension) in k
112 # correspond to curvature detectors in the following directions:
113 #
114 # [0] horizontal
115 # [1] vertical
116 # [2] diagonal \ (45 degrees rotation)
117 # [3] diagonal / (-45 degrees rotation)
118 image_g1_0 = scipy.ndimage.convolve(image, G1_0, mode="nearest")
119 image_g2_0 = scipy.ndimage.convolve(image, G2_0, mode="nearest")
120 image_g1_90 = scipy.ndimage.convolve(image, G1_90, mode="nearest")
121 image_g2_90 = scipy.ndimage.convolve(image, G2_90, mode="nearest")
122 fxy = scipy.ndimage.convolve(image, hxy, mode="nearest")
124 # support calculation for diagonals, given the gaussian kernel is
125 # steerable. To calculate the derivatives for the "\" diagonal, we first
126 # **would** have to rotate the image 45 degrees counter-clockwise (so the
127 # diagonal lies on the horizontal axis). Using the steerable property, we
128 # can evaluate the first derivative like this:
129 #
130 # image_g1_45 = cos(45)*image_g1_0 + sin(45)*image_g1_90
131 # = sqrt(2)/2*fx + sqrt(2)/2*fx
132 #
133 # to calculate the first derivative for the "/" diagonal, we first
134 # **would** have to rotate the image -45 degrees "counter"-clockwise.
135 # Therefore, we can calculate it like this:
136 #
137 # image_g1_m45 = cos(-45)*image_g1_0 + sin(-45)*image_g1_90
138 # = sqrt(2)/2*image_g1_0 - sqrt(2)/2*image_g1_90
139 #
141 image_g1_45 = 0.5 * numpy.sqrt(2) * (image_g1_0 + image_g1_90)
142 image_g1_m45 = 0.5 * numpy.sqrt(2) * (image_g1_0 - image_g1_90)
144 # NOTE: You can't really get image_g2_45 and image_g2_m45 from the theory
145 # of steerable filters. In contact with B.Ton, he suggested the following
146 # material, where that is explained: Chapter 5.2.3 of van der Heijden, F.
147 # (1994) Image based measurement systems: object recognition and parameter
148 # estimation. John Wiley & Sons Ltd, Chichester. ISBN 978-0-471-95062-2
150 # This also shows the same result:
151 # http://www.mif.vu.lt/atpazinimas/dip/FIP/fip-Derivati.html (look for
152 # SDGD)
154 # He also suggested to look at slide 75 of the following presentation
155 # indicating it is self-explanatory: http://slideplayer.com/slide/5084635/
157 image_g2_45 = 0.5 * image_g2_0 + fxy + 0.5 * image_g2_90
158 image_g2_m45 = 0.5 * image_g2_0 - fxy + 0.5 * image_g2_90
160 # ######################################################################
161 # [Step 1-1] Calculation of curvature profiles
162 # ######################################################################
164 # Peak detection (k or kappa) calculation as per equation (1) page 348 on
165 # Miura's paper
166 finger_mask = mask.astype("float64")
168 return numpy.dstack(
169 [
170 (image_g2_0 / ((1 + image_g1_0**2) ** (1.5))) * finger_mask,
171 (image_g2_90 / ((1 + image_g1_90**2) ** (1.5))) * finger_mask,
172 (image_g2_45 / ((1 + image_g1_45**2) ** (1.5))) * finger_mask,
173 (image_g2_m45 / ((1 + image_g1_m45**2) ** (1.5)))
174 * finger_mask,
175 ]
176 )
178 def eval_vein_probabilities(self, k):
179 r"""Evaluates joint vein centre probabilities from cross-sections
181 This function will take $\kappa$ and will calculate the vein centre
182 probabilities taking into consideration valley widths and depths. It
183 aggregates the following steps from the paper:
185 * [Step 1-2] Detection of the centres of veins
186 * [Step 1-3] Assignment of scores to the centre positions
187 * [Step 1-4] Calculation of all the profiles
189 Once the arrays of curvatures (concavities) are calculated, here is how
190 detection works: The code scans the image in a precise direction (vertical,
191 horizontal, diagonal, etc). It tries to find a concavity on that direction
192 and measure its width (see Wr on Figure 3 on the original paper). It then
193 identifies the centers of the concavity and assign a value to it, which
194 depends on its width (Wr) and maximum depth (where the peak of darkness
195 occurs) in such a concavity. This value is accumulated on a variable (Vt),
196 which is re-used for all directions. Vt represents the vein probabilites
197 from the paper.
200 Parameters:
202 k (numpy.ndarray): a 3-dimensional array of 64-bits containing $\kappa$
203 for all considered directions. $\kappa$ has the same shape as
204 ``image``, except for the 3rd. dimension, which provides planes for the
205 cross-section valley detections for each of the contemplated
206 directions, in this order: horizontal, vertical, +45 degrees, -45
207 degrees.
210 Returns:
212 numpy.ndarray: The un-accumulated vein centre probabilities ``V``. This
213 is a 3D array with 64-bit floats with the same dimensions of the input
214 array ``k``. You must accumulate (sum) over the last dimension to
215 retrieve the variable ``V`` from the paper.
217 """
219 V = numpy.zeros(k.shape[:2], dtype="float64")
221 def _prob_1d(a):
222 """Finds "vein probabilities" in a 1-D signal
224 This function efficiently counts the width and height of concavities in
225 the cross-section (1-D) curvature signal ``s``.
227 It works like this:
229 1. We create a 1-shift difference between the thresholded signal and
230 itself
231 2. We compensate for starting and ending regions
232 3. For each sequence of start/ends, we compute the maximum in the
233 original signal
235 Example (mixed with pseudo-code):
237 a = 0 1 2 3 2 1 0 -1 0 0 1 2 5 2 2 2 1
238 b = a > 0 (as type int)
239 b = 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1
241 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1
242 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 (-)
243 -------------------------------------------
244 X 1 0 0 0 0 -1 0 0 0 1 0 0 0 0 0 0 X (length is smaller than orig.)
246 starts = numpy.where(diff > 0)
247 ends = numpy.where(diff < 0)
249 -> now the number of starts and ends should match, otherwise, we must
250 compensate
252 -> case 1: b starts with 1: add one start in begin of "starts"
253 -> case 2: b ends with 1: add one end in the end of "ends"
255 -> iterate over the sequence of starts/ends and find maximums
258 Parameters:
260 a (numpy.ndarray): 1D signal with curvature to explore
263 Returns:
265 numpy.ndarray: 1D container with the vein centre probabilities
267 """
269 b = (a > 0).astype(int)
270 diff = b[1:] - b[:-1]
271 starts = numpy.argwhere(diff > 0)
272 starts += 1 # compensates for shifted different
273 ends = numpy.argwhere(diff < 0)
274 ends += 1 # compensates for shifted different
275 if b[0]:
276 starts = numpy.insert(starts, 0, 0)
277 if b[-1]:
278 ends = numpy.append(ends, len(a))
280 z = numpy.zeros_like(a)
282 if starts.size == 0 and ends.size == 0:
283 return z
285 for start, end in zip(starts, ends):
286 maximum = numpy.argmax(a[int(start) : int(end)])
287 z[start + maximum] = a[start + maximum] * (end - start)
289 return z
291 # Horizontal direction
292 for index in range(k.shape[0]):
293 V[index, :] += _prob_1d(k[index, :, 0])
295 # Vertical direction
296 for index in range(k.shape[1]):
297 V[:, index] += _prob_1d(k[:, index, 1])
299 # Direction: 45 degrees (\)
300 curv = k[:, :, 2]
301 i, j = numpy.indices(curv.shape)
302 for index in range(-curv.shape[0] + 1, curv.shape[1]):
303 V[i == (j - index)] += _prob_1d(curv.diagonal(index))
305 # Direction: -45 degrees (/)
306 # NOTE: due to the way the access to the diagonals are implemented, in this
307 # loop, we operate bottom-up. To match this behaviour, we also address V
308 # through Vud.
309 curv = numpy.flipud(
310 k[:, :, 3]
311 ) # required so we get "/" diagonals correctly
312 Vud = numpy.flipud(V) # match above inversion
313 for index in reversed(range(curv.shape[1] - 1, -curv.shape[0], -1)):
314 Vud[i == (j - index)] += _prob_1d(curv.diagonal(index))
316 return V
318 def connect_centres(self, V):
319 r"""Connects vein centres by filtering vein probabilities ``V``
321 This function does the equivalent of Step 2 / Equation 4 at Miura's paper.
323 The operation is applied on a row from the ``V`` matrix, which may be
324 acquired horizontally, vertically or on a diagonal direction. The pixel
325 value is then reset in the center of a windowing operation (width = 5) with
326 the following value:
328 .. math::
330 b[w] = min(max(a[w+1], a[w+2]) + max(a[w-1], a[w-2]))
333 Parameters:
335 V (numpy.ndarray): The accumulated vein centre probabilities ``V``. This
336 is a 2D array with 64-bit floats and is defined by Equation (3) on the
337 paper.
340 Returns:
342 numpy.ndarray: A 3-dimensional 64-bit array ``Cd`` containing the result
343 of the filtering operation for each of the directions. ``Cd`` has the
344 dimensions of $\kappa$ and $V_i$. Each of the planes correspond to the
345 horizontal, vertical, +45 and -45 directions.
347 """
349 def _connect_1d(a):
350 """Connects centres in the given vector
352 The strategy we use to vectorize this is to shift a twice to the left and
353 twice to the right and apply a vectorized operation to compute the above.
356 Parameters:
358 a (numpy.ndarray): Input 1D array which will be window scanned
361 Returns:
363 numpy.ndarray: Output 1D array (must be writeable), in which we will
364 set the corrected pixel values after the filtering above. Notice that,
365 given the windowing operation, the returned array size would be 4 short
366 of the input array.
368 """
370 return numpy.amin(
371 [
372 numpy.amax([a[3:-1], a[4:]], axis=0),
373 numpy.amax([a[1:-3], a[:-4]], axis=0),
374 ],
375 axis=0,
376 )
378 Cd = numpy.zeros(V.shape + (4,), dtype="float64")
380 # Horizontal direction
381 for index in range(V.shape[0]):
382 Cd[index, 2:-2, 0] = _connect_1d(V[index, :])
384 # Vertical direction
385 for index in range(V.shape[1]):
386 Cd[2:-2, index, 1] = _connect_1d(V[:, index])
388 # Direction: 45 degrees (\)
389 i, j = numpy.indices(V.shape)
390 border = numpy.zeros((2,), dtype="float64")
391 for index in range(-V.shape[0] + 5, V.shape[1] - 4):
392 # NOTE: hstack **absolutately** necessary here as double indexing after
393 # array indexing is **not** possible with numpy (it returns a copy)
394 Cd[:, :, 2][i == (j - index)] = numpy.hstack(
395 [border, _connect_1d(V.diagonal(index)), border]
396 )
398 # Direction: -45 degrees (/)
399 Vud = numpy.flipud(V)
400 Cdud = numpy.flipud(Cd[:, :, 3])
401 for index in reversed(range(V.shape[1] - 5, -V.shape[0] + 4, -1)):
402 # NOTE: hstack **absolutately** necessary here as double indexing after
403 # array indexing is **not** possible with numpy (it returns a copy)
404 Cdud[:, :][i == (j - index)] = numpy.hstack(
405 [border, _connect_1d(Vud.diagonal(index)), border]
406 )
408 return Cd
410 def binarise(self, G):
411 """Binarise vein images using a threshold assuming distribution is diphasic
413 This function implements Step 3 of the paper. It binarises the 2-D array
414 ``G`` assuming its histogram is mostly diphasic and using a median value.
417 Parameters:
419 G (numpy.ndarray): A 2-dimensional 64-bit array ``G`` containing the
420 result of the filtering operation. ``G`` has the dimensions of the
421 original image.
424 Returns:
426 numpy.ndarray: A 2-dimensional 64-bit float array with the same
427 dimensions of the input image, but containing its vein-binarised version.
428 The output of this function corresponds to the output of the method.
430 """
432 median = numpy.median(G[G > 0])
433 Gbool = G > median
434 return Gbool.astype(numpy.float64)
436 def _view_four(self, k, suptitle):
437 """Display four plots using matplotlib"""
439 import matplotlib.pyplot as plt
441 k[k <= 0] = 0
442 k /= k.max()
444 plt.subplot(2, 2, 1)
445 plt.imshow(k[..., 0], cmap="gray")
446 plt.title("Horizontal")
448 plt.subplot(2, 2, 2)
449 plt.imshow(k[..., 1], cmap="gray")
450 plt.title("Vertical")
452 plt.subplot(2, 2, 3)
453 plt.imshow(k[..., 2], cmap="gray")
454 plt.title("+45 degrees")
456 plt.subplot(2, 2, 4)
457 plt.imshow(k[..., 3], cmap="gray")
458 plt.title("-45 degrees")
460 plt.suptitle(suptitle)
461 plt.tight_layout()
462 plt.show()
464 def _view_single(self, k, title):
465 """Displays a single plot using matplotlib"""
467 import matplotlib.pyplot as plt
469 plt.imshow(k, cmap="gray")
470 plt.title(title)
471 plt.tight_layout()
472 plt.show()
474 def __call__(self, image):
475 finger_image = image[0].astype("float64")
476 finger_mask = image[1]
478 # import time
479 # start = time.time()
481 kappa = self.detect_valleys(finger_image, finger_mask)
483 # self._view_four(kappa, "Valley Detectors - $\kappa$")
485 # print('filtering took %.2f seconds' % (time.time() - start))
486 # start = time.time()
488 V = self.eval_vein_probabilities(kappa)
490 # self._view_single(V, "Accumulated Probabilities - V")
492 # print('probabilities took %.2f seconds' % (time.time() - start))
493 # start = time.time()
495 Cd = self.connect_centres(V)
497 # self._view_four(Cd, "Connected Centers - $C_{di}$")
498 # self._view_single(numpy.amax(Cd, axis=2), "Connected Centers - G")
500 # print('connections took %.2f seconds' % (time.time() - start))
501 # start = time.time()
503 retval = self.binarise(numpy.amax(Cd, axis=2))
505 # self._view_single(retval, "Final Binarised Image")
507 # print('binarization took %.2f seconds' % (time.time() - start))
509 return retval