diff --git a/bob/bio/vein/extractor/MaximumCurvature.py b/bob/bio/vein/extractor/MaximumCurvature.py index 0d891cfd1fe43e984d48707b901f114dcaa36141..8453b6bd25b02b13af0ba283c0f5649365a4157f 100644 --- a/bob/bio/vein/extractor/MaximumCurvature.py +++ b/bob/bio/vein/extractor/MaximumCurvature.py @@ -3,14 +3,10 @@ import math import numpy - -import bob.core +import scipy.ndimage import bob.io.base - from bob.bio.base.extractor import Extractor -from .. import utils - class MaximumCurvature (Extractor): """ @@ -18,12 +14,15 @@ class MaximumCurvature (Extractor): Based on N. Miura, A. Nagasaka, and T. Miyatake, Extraction of Finger-Vein Pattern Using Maximum Curvature Points in Image Profiles. Proceedings on IAPR - conference on machine vision applications, 9 (2005), pp. 347--350 + conference on machine vision applications, 9 (2005), pp. 347--350. + - **Parameters:** + Parameters: + + sigma (:py:class:int, optional): standard deviation for the gaussian + smoothing kernel used to denoise the input image. The width of the + gaussian kernel will be set automatically to 4x this value (in pixels). - sigma : :py:class:int - Optional: Sigma used for determining derivatives. """ @@ -32,227 +31,477 @@ class MaximumCurvature (Extractor): self.sigma = sigma - def maximum_curvature(self, image, mask): - """Computes and returns the Maximum Curvature features for the given input - fingervein image""" + def detect_valleys(self, image, mask): + """Detects valleys on the image respecting the mask + + This step corresponds to Step 1-1 in the original paper. The objective is, + for all 4 cross-sections (z) of the image (horizontal, vertical, 45 and -45 + diagonals), to compute the following proposed valley detector as defined in + Equation 1, page 348: + + .. math:: + + \kappa(z) = \\frac{d^2P_f(z)/dz^2}{(1 + (dP_f(z)/dz)^2)^\\frac{3}{2}} + + + We start the algorithm by smoothing the image with a 2-dimensional gaussian + filter. The equation that defines the kernel for the filter is: + + .. math:: + + \mathcal{N}(x,y)=\\frac{1}{2\pi\sigma^2}e^\\frac{-(x^2+y^2)}{2\sigma^2} + + + This is done to avoid noise from the raw data (from the sensor). The + maximum curvature method then requires we compute the first and second + derivative of the image for all cross-sections, as per the equation above. + + We instead take the following equivalent approach: + + 1. construct a gaussian filter + 2. take the first (dh/dx) and second (d^2/dh^2) deritivatives of the filter + 3. calculate the first and second derivatives of the smoothed signal using + the results from 3. This is done for all directions we're interested in: + horizontal, vertical and 2 diagonals. First and second derivatives of a + convolved signal + + .. note:: + + Item 3 above is only possible thanks to the steerable filter property of + the gaussian kernel. See "The Design and Use of Steerable Filters" from + Freeman and Adelson, IEEE Transactions on Pattern Analysis and Machine + Intelligence, Vol. 13, No. 9, September 1991. + + + Parameters: - finger_mask = numpy.zeros(mask.shape) - finger_mask[mask == True] = 1 + image (numpy.ndarray): an array of 64-bit floats containing the input + image + mask (numpy.ndarray): an array, of the same size as image, containing + a mask (booleans) indicating where the finger is on image. - winsize = numpy.ceil(4*self.sigma) - x = numpy.arange(-winsize, winsize+1) - y = numpy.arange(-winsize, winsize+1) - X, Y = numpy.meshgrid(x, y) + Returns: - h = (1/(2*math.pi*self.sigma**2))*numpy.exp(-(X**2 + Y**2)/(2*self.sigma**2)) - hx = (-X/(self.sigma**2))*h - hxx = ((X**2 - self.sigma**2)/(self.sigma**4))*h - hy = hx.T - hyy = hxx.T - hxy = ((X*Y)/(self.sigma**4))*h + numpy.ndarray: a 3-dimensional array of 64-bits containing $\kappa$ for + all considered directions. $\kappa$ has the same shape as image, + except for the 3rd. dimension, which provides planes for the + cross-section valley detections for each of the contemplated directions, + in this order: horizontal, vertical, +45 degrees, -45 degrees. - # Do the actual filtering + """ - fx = utils.imfilter(image, hx) - fxx = utils.imfilter(image, hxx) - fy = utils.imfilter(image, hy) - fyy = utils.imfilter(image, hyy) - fxy = utils.imfilter(image, hxy) + # 1. constructs the 2D gaussian filter "h" given the window size, + # extrapolated from the "sigma" parameter (4x) + # N.B.: This is a text-book gaussian filter definition + winsize = numpy.ceil(4*self.sigma) #enough space for the filter + window = numpy.arange(-winsize, winsize+1) + X, Y = numpy.meshgrid(window, window) + G = 1.0 / (2*math.pi*self.sigma**2) + G *= numpy.exp(-(X**2 + Y**2) / (2*self.sigma**2)) - f1 = 0.5*numpy.sqrt(2)*(fx + fy) # \ # - f2 = 0.5*numpy.sqrt(2)*(fx - fy) # / # - f11 = 0.5*fxx + fxy + 0.5*fyy # \\ # - f22 = 0.5*fxx - fxy + 0.5*fyy # // # + # 2. calculates first and second derivatives of "G" with respect to "X" + # (0), "Y" (90 degrees) and 45 degrees (?) + G1_0 = (-X/(self.sigma**2))*G + G2_0 = ((X**2 - self.sigma**2)/(self.sigma**4))*G + G1_90 = G1_0.T + G2_90 = G2_0.T + hxy = ((X*Y)/(self.sigma**4))*G + + # 3. calculates derivatives w.r.t. to all directions of interest + # stores results in the variable "k". The entries (last dimension) in k + # correspond to curvature detectors in the following directions: + # + # [0] horizontal + # [1] vertical + # [2] diagonal \ (45 degrees rotation) + # [3] diagonal / (-45 degrees rotation) + image_g1_0 = scipy.ndimage.convolve(image, G1_0, mode='nearest') + image_g2_0 = scipy.ndimage.convolve(image, G2_0, mode='nearest') + image_g1_90 = scipy.ndimage.convolve(image, G1_90, mode='nearest') + image_g2_90 = scipy.ndimage.convolve(image, G2_90, mode='nearest') + fxy = scipy.ndimage.convolve(image, hxy, mode='nearest') + + # support calculation for diagonals, given the gaussian kernel is + # steerable. To calculate the derivatives for the "\" diagonal, we first + # **would** have to rotate the image 45 degrees counter-clockwise (so the + # diagonal lies on the horizontal axis). Using the steerable property, we + # can evaluate the first derivative like this: + # + # image_g1_45 = cos(45)*image_g1_0 + sin(45)*image_g1_90 + # = sqrt(2)/2*fx + sqrt(2)/2*fx + # + # to calculate the first derivative for the "/" diagonal, we first + # **would** have to rotate the image -45 degrees "counter"-clockwise. + # Therefore, we can calculate it like this: + # + # image_g1_m45 = cos(-45)*image_g1_0 + sin(-45)*image_g1_90 + # = sqrt(2)/2*image_g1_0 - sqrt(2)/2*image_g1_90 + # + + image_g1_45 = 0.5*numpy.sqrt(2)*(image_g1_0 + image_g1_90) + image_g1_m45 = 0.5*numpy.sqrt(2)*(image_g1_0 - image_g1_90) + + # NOTE: You can't really get image_g2_45 and image_g2_m45 from the theory + # of steerable filters. In contact with B.Ton, he suggested the following + # material, where that is explained: Chapter 5.2.3 of van der Heijden, F. + # (1994) Image based measurement systems: object recognition and parameter + # estimation. John Wiley & Sons Ltd, Chichester. ISBN 978-0-471-95062-2 + + # This also shows the same result: + # http://www.mif.vu.lt/atpazinimas/dip/FIP/fip-Derivati.html (look for + # SDGD) + + # He also suggested to look at slide 75 of the following presentation + # indicating it is self-explanatory: http://slideplayer.com/slide/5084635/ + + image_g2_45 = 0.5*image_g2_0 + fxy + 0.5*image_g2_90 + image_g2_m45 = 0.5*image_g2_0 - fxy + 0.5*image_g2_90 img_h, img_w = image.shape #Image height and width - # Calculate curvatures - k = numpy.zeros((img_h, img_w, 4)) - k[:,:,0] = (fxx/((1 + fx**2)**(3/2)))*finger_mask # hor # - k[:,:,1] = (fyy/((1 + fy**2)**(3/2)))*finger_mask # ver # - k[:,:,2] = (f11/((1 + f1**2)**(3/2)))*finger_mask # \ # - k[:,:,3] = (f22/((1 + f2**2)**(3/2)))*finger_mask # / # + # ###################################################################### + # [Step 1-1] Calculation of curvature profiles + # ###################################################################### + + # Peak detection (k or kappa) calculation as per equation (1) page 348 on + # Miura's paper + finger_mask = mask.astype('float64') + + return numpy.dstack([ + (image_g2_0 / ((1 + image_g1_0**2)**(1.5)) ) * finger_mask, + (image_g2_90 / ((1 + image_g1_90**2)**(1.5)) ) * finger_mask, + (image_g2_45 / ((1 + image_g1_45**2)**(1.5)) ) * finger_mask, + (image_g2_m45 / ((1 + image_g1_m45**2)**(1.5))) * finger_mask, + ]) + + + def eval_vein_probabilities(self, k): + '''Evaluates joint vein centre probabilities from cross-sections + + This function will take $\kappa$ and will calculate the vein centre + probabilities taking into consideration valley widths and depths. It + aggregates the following steps from the paper: + + * [Step 1-2] Detection of the centres of veins + * [Step 1-3] Assignment of scores to the centre positions + * [Step 1-4] Calculation of all the profiles + + Once the arrays of curvatures (concavities) are calculated, here is how + detection works: The code scans the image in a precise direction (vertical, + horizontal, diagonal, etc). It tries to find a concavity on that direction + and measure its width (see Wr on Figure 3 on the original paper). It then + identifies the centers of the concavity and assign a value to it, which + depends on its width (Wr) and maximum depth (where the peak of darkness + occurs) in such a concavity. This value is accumulated on a variable (Vt), + which is re-used for all directions. Vt represents the vein probabilites + from the paper. + + + Parameters: + + k (numpy.ndarray): a 3-dimensional array of 64-bits containing $\kappa$ + for all considered directions. $\kappa$ has the same shape as + image, except for the 3rd. dimension, which provides planes for the + cross-section valley detections for each of the contemplated + directions, in this order: horizontal, vertical, +45 degrees, -45 + degrees. + + + Returns: + + numpy.ndarray: The un-accumulated vein centre probabilities V. This + is a 3D array with 64-bit floats with the same dimensions of the input + array k. You must accumulate (sum) over the last dimension to + retrieve the variable V from the paper. + + ''' + + V = numpy.zeros_like(k) + + def _prob_1d(a): + '''Finds "vein probabilities" in a 1-D signal + + This function efficiently counts the width and height of concavities in + the cross-section (1-D) curvature signal s. + + It works like this: + + 1. We create a 1-shift difference between the thresholded signal and + itself + 2. We compensate for starting and ending regions + 3. For each sequence of start/ends, we compute the maximum in the + original signal + + Example (mixed with pseudo-code): + + a = 0 1 2 3 2 1 0 -1 0 0 1 2 5 2 2 2 1 + b = a > 0 (as type int) + b = 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 + + 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 + 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 (-) + ------------------------------------------- + X 1 0 0 0 0 -1 0 0 0 1 0 0 0 0 0 0 X (length is smaller than orig.) + + starts = numpy.where(diff > 0) + ends = numpy.where(diff < 0) + + -> now the number of starts and ends should match, otherwise, we must + compensate + + -> case 1: b starts with 1: add one start in begin of "starts" + -> case 2: b ends with 1: add one end in the end of "ends" + + -> iterate over the sequence of starts/ends and find maximums + + + Parameters: + + a (numpy.ndarray): 1D signal with curvature to explore + + + Returns: + + numpy.ndarray: 1D container with the vein centre probabilities + + ''' + + b = (a > 0).astype(int) + diff = b[1:] - b[:-1] + starts = numpy.argwhere(diff > 0) + starts += 1 #compensates for shifted different + ends = numpy.argwhere(diff < 0) + ends += 1 #compensates for shifted different + if b[0]: starts = numpy.insert(starts, 0, 0) + if b[-1]: ends = numpy.append(ends, len(a)) + + z = numpy.zeros_like(a) + + if starts.size == 0 and ends.size == 0: return z + + for start, end in zip(starts, ends): + maximum = numpy.argmax(a[int(start):int(end)]) + z[start+maximum] = a[start+maximum] * (end-start) + + return z - # Scores - Vt = numpy.zeros(image.shape) - Wr = 0 # Horizontal direction - bla = k[:,:,0] > 0 - for y in range(0,img_h): - for x in range(0,img_w): - if (bla[y,x]): - Wr = Wr + 1 - if ( Wr > 0 and (x == (img_w-1) or not bla[y,x]) ): - if (x == (img_w-1)): - # Reached edge of image - pos_end = x - else: - pos_end = x - 1 - - pos_start = pos_end - Wr + 1 # Start pos of concave - if (pos_start == pos_end): - I=numpy.argmax(k[y,pos_start,0]) - else: - I=numpy.argmax(k[y,pos_start:pos_end+1,0]) - - pos_max = pos_start + I - Scr = k[y,pos_max,0]*Wr - Vt[y,pos_max] = Vt[y,pos_max] + Scr - Wr = 0 + for index in range(k.shape[0]): + V[index,:,0] += _prob_1d(k[index,:,0]) + + # Vertical direction + for index in range(k.shape[1]): + V[:,index,1] += _prob_1d(k[:,index,1]) + + # Direction: 45 degrees (\) + curv = k[:,:,2] + i,j = numpy.indices(curv.shape) + for index in range(-curv.shape[0]+1, curv.shape[1]): + V[i==(j-index),2] += _prob_1d(curv.diagonal(index)) + + # Direction: -45 degrees (/) + # NOTE: due to the way the access to the diagonals are implemented, in this + # loop, we operate bottom-up. To match this behaviour, we also address V + # through Vud. + curv = numpy.flipud(k[:,:,3]) #required so we get "/" diagonals correctly + Vud = numpy.flipud(V) #match above inversion + for index in reversed(range(curv.shape[1]-1, -curv.shape[0], -1)): + Vud[i==(j-index),3] += _prob_1d(curv.diagonal(index)) + + return V + + + def connect_centres(self, V): + """Connects vein centres by filtering vein probabilities V + + This function does the equivalent of Step 2 / Equation 4 at Miura's paper. + The operation is applied on a row from the V matrix, which may be + acquired horizontally, vertically or on a diagonal direction. The pixel + value is then reset in the center of a windowing operation (width = 5) with + the following value: + + .. math:: + + b[w] = min(max(a[w+1], a[w+2]) + max(a[w-1], a[w-2])) + + + Parameters: + + V (numpy.ndarray): The accumulated vein centre probabilities V. This + is a 2D array with 64-bit floats and is defined by Equation (3) on the + paper. + + + Returns: + + numpy.ndarray: A 3-dimensional 64-bit array Cd containing the result + of the filtering operation for each of the directions. Cd has the + dimensions of $\kappa$ and $V_i$. Each of the planes correspond to the + horizontal, vertical, +45 and -45 directions. + + """ + + def _connect_1d(a): + '''Connects centres in the given vector + + The strategy we use to vectorize this is to shift a twice to the left and + twice to the right and apply a vectorized operation to compute the above. + + + Parameters: + + a (numpy.ndarray): Input 1D array which will be window scanned + + + Returns: + + numpy.ndarray: Output 1D array (must be writeable), in which we will + set the corrected pixel values after the filtering above. Notice that, + given the windowing operation, the returned array size would be 4 short + of the input array. + + ''' + + return numpy.amin([numpy.amax([a[3:-1], a[4:]], axis=0), + numpy.amax([a[1:-3], a[:-4]], axis=0)], axis=0) + + + Cd = numpy.zeros(V.shape + (4,), dtype='float64') + + # Horizontal direction + for index in range(V.shape[0]): + Cd[index, 2:-2, 0] = _connect_1d(V[index,:]) # Vertical direction - bla = k[:,:,1] > 0 - for x in range(0,img_w): - for y in range(0,img_h): - if (bla[y,x]): - Wr = Wr + 1 - if ( Wr > 0 and (y == (img_h-1) or not bla[y,x]) ): - if (y == (img_h-1)): - # Reached edge of image - pos_end = y - else: - pos_end = y - 1 - - pos_start = pos_end - Wr + 1 # Start pos of concave - if (pos_start == pos_end): - I=numpy.argmax(k[pos_start,x,1]) - else: - I=numpy.argmax(k[pos_start:pos_end+1,x,1]) - - pos_max = pos_start + I - Scr = k[pos_max,x,1]*Wr - - Vt[pos_max,x] = Vt[pos_max,x] + Scr - Wr = 0 - - # Direction: \ # - bla = k[:,:,2] > 0 - for start in range(0,img_w+img_h-1): - # Initial values - if (start <= img_w-1): - x = start - y = 0 - else: - x = 0 - y = start - img_w + 1 - done = False - - while (not done): - if(bla[y,x]): - Wr = Wr + 1 - - if ( Wr > 0 and (y == img_h-1 or x == img_w-1 or not bla[y,x]) ): - if (y == img_h-1 or x == img_w-1): - # Reached edge of image - pos_x_end = x - pos_y_end = y - else: - pos_x_end = x - 1 - pos_y_end = y - 1 - - pos_x_start = pos_x_end - Wr + 1 - pos_y_start = pos_y_end - Wr + 1 - - if (pos_y_start == pos_y_end and pos_x_start == pos_x_end): - d = k[pos_y_start, pos_x_start, 2] - elif (pos_y_start == pos_y_end): - d = numpy.diag(k[pos_y_start, pos_x_start:pos_x_end+1, 2]) - elif (pos_x_start == pos_x_end): - d = numpy.diag(k[pos_y_start:pos_y_end+1, pos_x_start, 2]) - else: - d = numpy.diag(k[pos_y_start:pos_y_end+1, pos_x_start:pos_x_end+1, 2]) - - I = numpy.argmax(d) - - pos_x_max = pos_x_start + I - pos_y_max = pos_y_start + I - - Scr = k[pos_y_max,pos_x_max,2]*Wr - - Vt[pos_y_max,pos_x_max] = Vt[pos_y_max,pos_x_max] + Scr - Wr = 0 - - if((x == img_w-1) or (y == img_h-1)): - done = True - else: - x = x + 1 - y = y + 1 - - # Direction: / - bla = k[:,:,3] > 0 - for start in range(0,img_w+img_h-1): - # Initial values - if (start <= (img_w-1)): - x = start - y = img_h-1 - else: - x = 0 - y = img_w+img_h-start-1 - done = False - - while (not done): - if(bla[y,x]): - Wr = Wr + 1 - if ( Wr > 0 and (y == 0 or x == img_w-1 or not bla[y,x]) ): - if (y == 0 or x == img_w-1): - # Reached edge of image - pos_x_end = x - pos_y_end = y - else: - pos_x_end = x - 1 - pos_y_end = y + 1 - - pos_x_start = pos_x_end - Wr + 1 - pos_y_start = pos_y_end + Wr - 1 - - if (pos_y_start == pos_y_end and pos_x_start == pos_x_end): - d = k[pos_y_end, pos_x_start, 3] - elif (pos_y_start == pos_y_end): - d = numpy.diag(numpy.flipud(k[pos_y_end, pos_x_start:pos_x_end+1, 3])) - elif (pos_x_start == pos_x_end): - d = numpy.diag(numpy.flipud(k[pos_y_end:pos_y_start+1, pos_x_start, 3])) - else: - d = numpy.diag(numpy.flipud(k[pos_y_end:pos_y_start+1, pos_x_start:pos_x_end+1, 3])) - - I = numpy.argmax(d) - pos_x_max = pos_x_start + I - pos_y_max = pos_y_start - I - Scr = k[pos_y_max,pos_x_max,3]*Wr - Vt[pos_y_max,pos_x_max] = Vt[pos_y_max,pos_x_max] + Scr - Wr = 0 - - if((x == img_w-1) or (y == 0)): - done = True - else: - x = x + 1 - y = y - 1 - - ## Connection of vein centres - Cd = numpy.zeros((img_h, img_w, 4)) - for x in range(2,img_w-3): - for y in range(2,img_h-3): - Cd[y,x,0] = min(numpy.amax(Vt[y,x+1:x+3]), numpy.amax(Vt[y,x-2:x])) # Hor # - Cd[y,x,1] = min(numpy.amax(Vt[y+1:y+3,x]), numpy.amax(Vt[y-2:y,x])) # Vert # - Cd[y,x,2] = min(numpy.amax(Vt[y-2:y,x-2:x]), numpy.amax(Vt[y+1:y+3,x+1:x+3])) # \ # - Cd[y,x,3] = min(numpy.amax(Vt[y+1:y+3,x-2:x]), numpy.amax(Vt[y-2:y,x+1:x+3])) # / # - - #Veins - img_veins = numpy.amax(Cd,axis=2) - - # Binarise the vein image - md = numpy.median(img_veins[img_veins>0]) - img_veins_bin = img_veins > md - - return img_veins_bin.astype(numpy.float64) + for index in range(V.shape[1]): + Cd[2:-2, index, 1] = _connect_1d(V[:,index]) + + # Direction: 45 degrees (\) + i,j = numpy.indices(V.shape) + border = numpy.zeros((2,), dtype='float64') + for index in range(-V.shape[0]+5, V.shape[1]-4): + # NOTE: hstack **absolutately** necessary here as double indexing after + # array indexing is **not** possible with numpy (it returns a copy) + Cd[:,:,2][i==(j-index)] = numpy.hstack([border, + _connect_1d(V.diagonal(index)), border]) + + # Direction: -45 degrees (/) + Vud = numpy.flipud(V) + Cdud = numpy.flipud(Cd[:,:,3]) + for index in reversed(range(V.shape[1]-5, -V.shape[0]+4, -1)): + # NOTE: hstack **absolutately** necessary here as double indexing after + # array indexing is **not** possible with numpy (it returns a copy) + Cdud[:,:][i==(j-index)] = numpy.hstack([border, + _connect_1d(Vud.diagonal(index)), border]) + + return Cd + + + def binarise(self, G): + """Binarise vein images using a threshold assuming distribution is diphasic + + This function implements Step 3 of the paper. It binarises the 2-D array + G assuming its histogram is mostly diphasic and using a median value. + + + Parameters: + + G (numpy.ndarray): A 2-dimensional 64-bit array G containing the + result of the filtering operation. G has the dimensions of the + original image. + + + Returns: + + numpy.ndarray: A 2-dimensional 64-bit float array with the same + dimensions of the input image, but containing its vein-binarised version. + The output of this function corresponds to the output of the method. + + """ + + median = numpy.median(G[G>0]) + Gbool = G > median + return Gbool.astype(numpy.float64) + + + def _view_four(self, k, suptitle): + '''Display four plots using matplotlib''' + + import matplotlib.pyplot as plt + + k[k<=0] = 0 + k /= k.max() + + plt.subplot(2,2,1) + plt.imshow(k[...,0], cmap='gray') + plt.title('Horizontal') + + plt.subplot(2,2,2) + plt.imshow(k[...,1], cmap='gray') + plt.title('Vertical') + + plt.subplot(2,2,3) + plt.imshow(k[...,2], cmap='gray') + plt.title('+45 degrees') + + plt.subplot(2,2,4) + plt.imshow(k[...,3], cmap='gray') + plt.title('-45 degrees') + + plt.suptitle(suptitle) + plt.tight_layout() + plt.show() + + + def _view_single(self, k, title): + '''Displays a single plot using matplotlib''' + + import matplotlib.pyplot as plt + + plt.imshow(k, cmap='gray') + plt.title(title) + plt.tight_layout() + plt.show() def __call__(self, image): - """Reads the input image, extract the features based on Maximum Curvature of the fingervein image, and writes the resulting template""" - finger_image = image[0] #Normalized image with or without histogram equalization + finger_image = image[0] finger_mask = image[1] - return self.maximum_curvature(finger_image, finger_mask) + import time + start = time.time() + + kappa = self.detect_valleys(finger_image, finger_mask) + + #self._view_four(kappa, "Valley Detectors - $\kappa$") + + print('filtering took %.2f seconds' % (time.time() - start)) + start = time.time() + + V = self.eval_vein_probabilities(kappa) + + #self._view_four(V, "Center Probabilities - $V_i$") + #self._view_single(V.sum(axis=2), "Accumulated Probabilities - V") + + print('probabilities took %.2f seconds' % (time.time() - start)) + start = time.time() + + Cd = self.connect_centres(V.sum(axis=2)) + + #self._view_four(Cd, "Connected Centers - $C_{di}$") + #self._view_single(numpy.amax(Cd, axis=2), "Connected Centers - G") + + print('connections took %.2f seconds' % (time.time() - start)) + start = time.time() + + retval = self.binarise(numpy.amax(Cd, axis=2)) + + #self._view_single(retval, "Final Binarised Image") + + print('binarization took %.2f seconds' % (time.time() - start)) + + return retval diff --git a/bob/bio/vein/tests/extractors/image.hdf5 b/bob/bio/vein/tests/extractors/image.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..bafea75e24b675df88c5879c589fec4f2d31ee7f Binary files /dev/null and b/bob/bio/vein/tests/extractors/image.hdf5 differ diff --git a/bob/bio/vein/tests/extractors/mask.hdf5 b/bob/bio/vein/tests/extractors/mask.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..e0770189463109fa56b0b21d6672bc9fb56862c2 Binary files /dev/null and b/bob/bio/vein/tests/extractors/mask.hdf5 differ diff --git a/bob/bio/vein/tests/extractors/mc_bin_matlab.hdf5 b/bob/bio/vein/tests/extractors/mc_bin_matlab.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..df5dac3b258381a796b27087ba934c21a9b5520f Binary files /dev/null and b/bob/bio/vein/tests/extractors/mc_bin_matlab.hdf5 differ diff --git a/bob/bio/vein/tests/extractors/mc_g_matlab.hdf5 b/bob/bio/vein/tests/extractors/mc_g_matlab.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..64aaa64d820e71d74b9c6991ba586f9d97e9005e Binary files /dev/null and b/bob/bio/vein/tests/extractors/mc_g_matlab.hdf5 differ diff --git a/bob/bio/vein/tests/extractors/mc_vt_matlab.hdf5 b/bob/bio/vein/tests/extractors/mc_vt_matlab.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..e0fcd7d21c735cbce780ac541cbbdbb7dc9ccdaa Binary files /dev/null and b/bob/bio/vein/tests/extractors/mc_vt_matlab.hdf5 differ diff --git a/bob/bio/vein/tests/extractors/miuramax_input_fvr.mat b/bob/bio/vein/tests/extractors/miuramax_input_fvr.mat deleted file mode 100644 index 2ab2dfd863819a018b46b64259a4b70b2ab97461..0000000000000000000000000000000000000000 Binary files a/bob/bio/vein/tests/extractors/miuramax_input_fvr.mat and /dev/null differ diff --git a/bob/bio/vein/tests/extractors/miuramax_input_img.mat b/bob/bio/vein/tests/extractors/miuramax_input_img.mat deleted file mode 100644 index 95e874fb7d0429b14d4543cee964fef7eb64abd5..0000000000000000000000000000000000000000 Binary files a/bob/bio/vein/tests/extractors/miuramax_input_img.mat and /dev/null differ diff --git a/bob/bio/vein/tests/extractors/miuramax_output.mat b/bob/bio/vein/tests/extractors/miuramax_output.mat deleted file mode 100644 index 6a3de98b8c500a73078a21dc44e5acc8ddb99563..0000000000000000000000000000000000000000 Binary files a/bob/bio/vein/tests/extractors/miuramax_output.mat and /dev/null differ diff --git a/bob/bio/vein/tests/test.py b/bob/bio/vein/tests/test.py index ada3a28a0943f680f6b0cbb2896a82ab00d825f5..dbc595313a5c0db088594fc6f7ede6437eb832ab 100644 --- a/bob/bio/vein/tests/test.py +++ b/bob/bio/vein/tests/test.py @@ -42,9 +42,15 @@ def test_finger_crop(): img = bob.io.base.load(input_filename) - from bob.bio.vein.preprocessor.FingerCrop import FingerCrop - preprocess = FingerCrop(fingercontour='leemaskMatlab', padding_width=0) - preproc, mask = preprocess(img) + from bob.bio.vein.preprocessor import Preprocessor, LeeMask, \ + HuangNormalization, NoFilter + + processor = Preprocessor( + LeeMask(filter_height=40, filter_width=4), + HuangNormalization(padding_width=0, padding_constant=0), + NoFilter(), + ) + preproc, mask = processor(img) #preprocessor_utils.show_mask_over_image(preproc, mask) mask_ref = bob.io.base.load(output_fvr_filename).astype('bool') @@ -53,7 +59,7 @@ def test_finger_crop(): assert numpy.mean(numpy.abs(mask - mask_ref)) < 1e-2 - # Very loose comparison! + # 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 @@ -62,25 +68,41 @@ def test_max_curvature(): #Maximum Curvature method against Matlab reference - input_img_filename = F(('extractors', 'miuramax_input_img.mat')) - input_fvr_filename = F(('extractors', 'miuramax_input_fvr.mat')) - output_filename = F(('extractors', 'miuramax_output.mat')) - - # Load inputs - input_img = bob.io.base.load(input_img_filename) - input_fvr = bob.io.base.load(input_fvr_filename) + image = bob.io.base.load(F(('extractors', 'image.hdf5'))) + image = image.T + image = image.astype('float64')/255. + mask = bob.io.base.load(F(('extractors', 'mask.hdf5'))) + mask = mask.T + mask = mask.astype('bool') + vt_ref = bob.io.base.load(F(('extractors', 'mc_vt_matlab.hdf5'))) + vt_ref = vt_ref.T + g_ref = bob.io.base.load(F(('extractors', 'mc_g_matlab.hdf5'))) + g_ref = g_ref.T + bin_ref = bob.io.base.load(F(('extractors', 'mc_bin_matlab.hdf5'))) + bin_ref = bin_ref.T # Apply Python implementation from bob.bio.vein.extractor.MaximumCurvature import MaximumCurvature - MC = MaximumCurvature(5) - output_img = MC((input_img, input_fvr)) - - # Load Matlab reference - output_img_ref = bob.io.base.load(output_filename) - - # Compare output of python's implementation to matlab reference - # (loose comparison!) - assert numpy.mean(numpy.abs(output_img - output_img_ref)) < 8e-3 + MC = MaximumCurvature(3) #value used to create references + + kappa = MC.detect_valleys(image, mask) + V = MC.eval_vein_probabilities(kappa) + Vt = V.sum(axis=2) + Cd = MC.connect_centres(Vt) + G = numpy.amax(Cd, axis=2) + bina = MC.binarise(G) + + assert numpy.allclose(Vt, vt_ref, 1e-3, 1e-4), \ + 'Vt differs from reference by %s' % numpy.abs(Vt-vt_ref).sum() + # Note: due to Matlab implementation bug, can only compare in a limited + # range with a 3-pixel around frame + assert numpy.allclose(G[2:-3,2:-3], g_ref[2:-3,2:-3]), \ + 'G differs from reference by %s' % numpy.abs(G-g_ref).sum() + # We require no more than 30 pixels (from a total of 63'840) are different + # between ours and the matlab implementation + assert numpy.abs(bin_ref-bina).sum() < 30, \ + 'Binarized image differs from reference by %s' % \ + numpy.abs(bin_ref-bina).sum() def test_max_curvature_HE(): @@ -91,9 +113,14 @@ def test_max_curvature_HE(): input_img = bob.io.base.load(input_img_filename) # Preprocess the data and apply Histogram Equalization postprocessing (same parameters as in maximum_curvature.py configuration file + postprocessing) - from bob.bio.vein.preprocessor.FingerCrop import FingerCrop - FC = FingerCrop(postprocessing = 'HE') - preproc_data = FC(input_img) + from bob.bio.vein.preprocessor import Preprocessor, LeeMask, \ + HuangNormalization, HistogramEqualization + processor = Preprocessor( + LeeMask(filter_height=40, filter_width=4), + HuangNormalization(padding_width=0, padding_constant=0), + HistogramEqualization(), + ) + preproc_data = processor(input_img) # Extract features from preprocessed and histogram equalized data using MC extractor (same parameters as in maximum_curvature.py configuration file) from bob.bio.vein.extractor.MaximumCurvature import MaximumCurvature @@ -134,9 +161,14 @@ def test_repeated_line_tracking_HE(): input_img = bob.io.base.load(input_img_filename) # Preprocess the data and apply Histogram Equalization postprocessing (same parameters as in repeated_line_tracking.py configuration file + postprocessing) - from bob.bio.vein.preprocessor.FingerCrop import FingerCrop - FC = FingerCrop(postprocessing = 'HE') - preproc_data = FC(input_img) + from bob.bio.vein.preprocessor import Preprocessor, LeeMask, \ + HuangNormalization, HistogramEqualization + processor = Preprocessor( + LeeMask(filter_height=40, filter_width=4), + HuangNormalization(padding_width=0, padding_constant=0), + HistogramEqualization(), + ) + preproc_data = processor(input_img) # Extract features from preprocessed and histogram equalized data using RLT extractor (same parameters as in repeated_line_tracking.py configuration file) from bob.bio.vein.extractor.RepeatedLineTracking import RepeatedLineTracking @@ -182,9 +214,14 @@ def test_wide_line_detector_HE(): input_img = bob.io.base.load(input_img_filename) # Preprocess the data and apply Histogram Equalization postprocessing (same parameters as in wide_line_detector.py configuration file + postprocessing) - from bob.bio.vein.preprocessor.FingerCrop import FingerCrop - FC = FingerCrop(postprocessing = 'HE') - preproc_data = FC(input_img) + from bob.bio.vein.preprocessor import Preprocessor, LeeMask, \ + HuangNormalization, HistogramEqualization + processor = Preprocessor( + LeeMask(filter_height=40, filter_width=4), + HuangNormalization(padding_width=0, padding_constant=0), + HistogramEqualization(), + ) + preproc_data = processor(input_img) # Extract features from preprocessed and histogram equalized data using WLD extractor (same parameters as in wide_line_detector.py configuration file) from bob.bio.vein.extractor.WideLineDetector import WideLineDetector