Add matlab reference generation and readme

matlab/README.md

+# Maximum Curvature and Miura Matching (Cross-correlation)
+
+This directory contains scripts to generate references from the Matlab library
+for Repeated Line Tracking (RLT) and Maximum Curvature extraction from B.Ton.
+The original source code download link is here: http://ch.mathworks.com/matlabcentral/fileexchange/35716-miura-et-al-vein-extraction-methods
+
+
+## Usage Instructions
+
+Make sure you have the `matlab` binary available on your path. At Idiap, this
+can be done by executing:
+
+```sh
+$ SETSHELL matlab
+$ which matlab
+/idiap/resource/software/matlab/stable/bin/matlab
+```
+
+Call `run.sh`, that will perform the maximum curvature on the provided image,
+and output various control variables in various HDF5 files:
+
+```sh
+$ run.sh ../bob/bio/vein/tests/extractors/image.hdf5 ../bob/bio/vein/tests/extractors/mask.hdf5 mc
+...
+```
+
+Or, a quicker way, without contaminating your environment:
+
+```sh
+$ setshell.py matlab ./run.sh ../bob/bio/vein/tests/extractors/image.hdf5 ../bob/bio/vein/tests/extractors/mask.hdf5 mc
+...
+```
+
+The program `run.sh` calls the function `lib/maximum_curvature.m`, which
+processes and dumps results to output files.
+
+After running, this program produces 5 files:
+
+* `*_kappa_matlab.hdf5`: contains the kappa variable
+* `*_v_matlab.hdf5`: contains the V variable
+* `*_vt_matlab.hdf5`: contains the accumulated Vt variable
+* `*_cd_matlab.hdf5`: contains the Cd variable
+* `*_g_matlab.hdf5`: contains the accumulated Cd variable called G (paper)
+* `*_bin_matlab.hdf5:`: the final binarised results (G thresholded)
+
+Once you have generated the references, you can compare the Maximum Curvature
+algorithm implemented in Python against the one in Matlab with the following
+command:
+
+```sh
+$ ../bin/python compare.py
+Comparing kappa[0]: 2.51624726501e-14
+Comparing kappa[1]: 2.10662186414e-13
+Comparing kappa[2]: 1.09901515494e-13
+Comparing kappa[3]: 1.0902340027e-13
+Comparing V[0]: 1.04325752096e-14
+Comparing V[1]: 8.5519523202e-14
+Comparing V[2]: 9.20009110336e-05
+Comparing V[3]: 4.02339072347e-14
+Comparing Vt: 9.20009111675e-05
+Comparing Cd[0]: 1.08636347658e-13
+Comparing Cd[1]: 2.8698038577e-14
+Comparing Cd[2]: 3.40774680626e-14
+Comparing Cd[3]: 3.2892922132e-14
+Comparing G: 1.57966982699e-13
+```
+
+The program displays the sum of absolute differences between the Matlab
+reference and Python's.

matlab/compare.py

+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+
+# Run this script to output a debugging comparison of the Python implementation
+# against matlab references you just extracted
+
+import numpy
+import bob.io.base
+import bob.io.matlab
+from bob.bio.vein.extractor import MaximumCurvature
+
+# Load inputs
+image  = bob.io.base.load('../bob/bio/vein/tests/extractors/image.hdf5')
+image  = image.T.astype('float64')/255.
+region = bob.io.base.load('../bob/bio/vein/tests/extractors/mask.hdf5')
+region = region.T.astype('bool')
+
+# Loads matlab references
+kappa_matlab = bob.io.base.load('mc_kappa_matlab.hdf5')
+kappa_matlab = kappa_matlab.transpose(2,1,0)
+V_matlab = bob.io.base.load('mc_v_matlab.hdf5')
+V_matlab = V_matlab.transpose(2,1,0)
+Vt_matlab = bob.io.base.load('mc_vt_matlab.hdf5')
+Vt_matlab = Vt_matlab.T
+Cd_matlab = bob.io.base.load('mc_cd_matlab.hdf5')
+Cd_matlab = Cd_matlab.transpose(2,1,0)
+G_matlab = bob.io.base.load('mc_g_matlab.hdf5')
+G_matlab = G_matlab.T
+
+# Apply Python implementation
+from bob.bio.vein.extractor.MaximumCurvature import MaximumCurvature
+MC = MaximumCurvature(3)
+
+kappa = MC.detect_valleys(image, region) #OK
+V = MC.eval_vein_probabilities(kappa) #OK
+Vt = V.sum(axis=2) #OK
+Cd = MC.connect_centres(Vt) #OK
+G = numpy.amax(Cd, axis=2) #OK
+
+# Compare outputs
+for k in range(4):
+  print('Comparing kappa[%d]: %s' % (k,
+    numpy.abs(kappa[...,k]-kappa_matlab[...,k]).sum()))
+
+for k in range(4):
+  print('Comparing V[%d]: %s' % (k, numpy.abs(V[...,k]-V_matlab[...,k]).sum()))
+
+print('Comparing Vt: %s' % numpy.abs(Vt-Vt_matlab).sum())
+
+for k in range(4):
+  print('Comparing Cd[%d]: %s' % (k,
+    numpy.abs(Cd[2:-3,2:-3,k]-Cd_matlab[2:-3,2:-3,k]).sum()))
+
+print('Comparing G: %s' % numpy.abs(G[2:-3,2:-3]-G_matlab[2:-3,2:-3]).sum())

matlab/lib/license.txt

+Copyright (c) 2012, Bram Ton
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in
+      the documentation and/or other materials provided with the distribution
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.

matlab/lib/miura_match.m

+function score = miura_match(I, R, cw, ch)
+% This is the matching procedure described by Miura et al. in their paper.
+% A small difference is that this matching function calculates the match 
+% ratio instead of the mismatch ratio.
+
+% Parameters:
+%  I  - Input image
+%  R  - Registered template image
+%  cw - Maximum search displacement in x-direction
+%  ch - Maximum search displacement in y-direction
+
+% Returns:
+%  score - Value between 0 and 0.5, larger value is better match
+
+% Reference:
+% Feature extraction of finger vein patterns based on iterative line 
+%    tracking and its application to personal identification
+% N. Miura, A. Nagasaka, and T. Miyatake
+% Syst. Comput. Japan 35 (7 June 2004), pp. 61--71
+% doi: 10.1002/scj.v35:7
+
+% Author:  Bram Ton <b.t.ton@alumnus.utwente.nl>
+% Date:    20th December 2011
+% License: Simplified BSD License
+
+[h w] = size(R); % Get height and width of registered data
+
+% Determine value of match, just cross-correlation, see also xcorr2
+Nm = conv2(I, rot90(R(ch+1:h-ch, cw+1:w-cw),2), 'valid');
+
+% Maximum value of match
+[Nmm,mi] = max(Nm(:)); % (what about multiple maximum values ?)
+[t0,s0] = ind2sub(size(Nm),mi);
+
+% Normalize
+score = Nmm/(sum(sum(R(ch+1:h-ch, cw+1:w-cw))) + sum(sum(I(t0:t0+h-2*ch-1, s0:s0+w-2*cw-1))));
+
+% %% Bram Test
+% Ipad = zeros(h+2*ch,w+2*cw);
+% Ipad(ch+1:ch+h,cw+1:cw+w) = I;
+% 
+% Nm = conv2(Ipad, rot90(R,2), 'valid');
+% 
+% % Maximum value of match
+% [Nmm,mi] = max(Nm(:)); % (what about multiple maximum values ?)
+% [t0,s0] = ind2sub(size(Nm),mi);
+% 
+% % Normalize
+% score = Nmm/(sum(sum(R(ch+1:h-ch, cw+1:w-cw))) + sum(sum(Ipad(t0:t0+h-2*ch-1, s0:s0+w-2*cw-1))));
+% %score = max(max(normxcorr2(R(ch+1:h-ch, cw+1:w-cw),I)));
+
+

matlab/lib/miura_max_curvature.m

+function veins = miura_max_curvature(img, fvr, sigma, stem)
+% Maximum curvature method
+
+% Parameters:
+%  img    - Input vascular image
+%  fvr    - Finger vein region
+%  sigma  - Sigma used for determining derivatives
+%  stem   - (optional) name of the stem where to save partial results
+
+% Returns:
+%  veins - Vein image
+
+% Reference:
+% Extraction of finger-vein patterns using maximum curvature points in
+%   image profiles
+% N. Miura, A. Nagasaka, and T. Miyatake
+% IAPR conference on machine vision applications 9 (2005), pp. 347--350
+
+% Author:  Bram Ton <b.t.ton@alumnus.utwente.nl>
+% Date:    20th December 2011
+% License: Simplified BSD License
+
+% Changelog:
+% 2012/01/10 - Speed enhancement by rewriting diag() functions
+%              with linear indices.
+
+% Construct filter kernels
+winsize = ceil(4*sigma);
+[X,Y] = meshgrid(-winsize:winsize, -winsize:winsize);
+
+h = (1/(2*pi*sigma^2)).*exp(-(X.^2 + Y.^2)/(2*sigma^2));
+hx = (-X/(sigma^2)).*h;
+hxx = ((X.^2 - sigma^2)/(sigma^4)).*h;
+hy = hx';
+hyy = hxx';
+hxy = ((X.*Y)/(sigma^4)).*h;
+
+% Do the actual filtering
+fx  = imfilter(img, hx,  'replicate', 'conv');
+fxx = imfilter(img, hxx, 'replicate', 'conv');
+fy  = imfilter(img, hy,  'replicate', 'conv');
+fyy = imfilter(img, hyy, 'replicate', 'conv');
+fxy = imfilter(img, hxy, 'replicate', 'conv');
+f1  = 0.5*sqrt(2)*(fx + fy); % \
+f2  = 0.5*sqrt(2)*(fx - fy); % /
+f11 = 0.5*fxx + fxy + 0.5*fyy; % \\
+f22 = 0.5*fxx - fxy + 0.5*fyy; % //
+
+[img_h, img_w] = size(img); % Image height and width
+
+%% Calculate curvatures
+k = zeros(img_h, img_w, 4);
+k(:,:,1) = (fxx./((1 + fx.^2).^(3/2))).*fvr; % hor
+k(:,:,2) = (fyy./((1 + fy.^2).^(3/2))).*fvr; % ver
+k(:,:,3) = (f11./((1 + f1.^2).^(3/2))).*fvr; % \
+k(:,:,4) = (f22./((1 + f2.^2).^(3/2))).*fvr; % /
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_kappa_matlab.hdf5'), '/kappa', k);
+end
+
+%% Scores
+V = zeros(img_h, img_w, 4);
+Vt = zeros(img_h, img_w);
+Wr = 0;
+
+% Horizontal direction
+bla = k(:,:,1) > 0;
+for y=1:img_h
+    for x=1:img_w
+        if(bla(y,x))
+            Wr = Wr + 1;
+        end
+
+        if ( Wr > 0 && (x == img_w || ~bla(y,x)) )
+            if (x == img_w)
+                % Reached edge of image
+                pos_end = x;
+            else
+                pos_end = x - 1;
+            end
+
+            pos_start = pos_end - Wr + 1; % Start pos of concave
+            [~, I] = max(k(y, pos_start:pos_end,1));
+            pos_max = pos_start + I - 1;
+            Scr = k(y,pos_max,1)*Wr;
+            V(y,pos_max,1) = V(y,pos_max,1) + Scr;
+            Vt(y,pos_max) = Vt(y,pos_max) + Scr;
+            Wr = 0;
+        end
+    end
+end
+
+% Vertical direction
+bla = k(:,:,2) > 0;
+for x=1:img_w
+    for y=1:img_h
+        if(bla(y,x))
+            Wr = Wr + 1;
+        end
+
+        if ( Wr > 0 && (y == img_h || ~bla(y,x)) )
+            if (x == img_h)
+                % Reached edge of image
+                pos_end = y;
+            else
+                pos_end = y - 1;
+            end
+
+            pos_start = pos_end - Wr + 1; % Start pos of concave
+            [~, I] = max(k(pos_start:pos_end,x,2));
+            pos_max = pos_start + I - 1;
+            Scr = k(pos_max,x,2)*Wr;
+            V(pos_max,x,2) = V(pos_max,x,2) + Scr;
+            Vt(pos_max,x) = Vt(pos_max,x) + Scr;
+            Wr = 0;
+        end
+    end
+end
+
+% Direction: \
+bla = k(:,:,3) > 0;
+for start=1:(img_w+img_h-1)
+    % Initial values
+    if (start <= img_w)
+        x = start;
+        y = 1;
+    else
+        x = 1;
+        y = start - img_w + 1;
+    end
+    done = false;
+
+    while ~done
+        if(bla(y,x))
+            Wr = Wr + 1;
+        end
+
+        if ( Wr > 0 && (y == img_h || x == img_w || ~bla(y,x)) )
+            if (y == img_h || x == img_w)
+                % Reached edge of image
+                pos_x_end = x;
+                pos_y_end = y;
+            else
+                pos_x_end = x - 1;
+                pos_y_end = y - 1;
+            end
+            pos_x_start = pos_x_end - Wr + 1;
+            pos_y_start = pos_y_end - Wr + 1;
+
+            %d = diag(k(pos_y_start:pos_y_end, pos_x_start:pos_x_end, 3));
+            % More efficient implementation than diag(..)
+            d = k(((pos_x_start-1)*img_h + pos_y_start + 2*img_w*img_h):(img_h + 1):((pos_x_end-1)*img_h + pos_y_end + 2*img_w*img_h));
+            [~, I] = max(d);
+            pos_x_max = pos_x_start + I - 1;
+            pos_y_max = pos_y_start + I - 1;
+            Scr = k(pos_y_max,pos_x_max,3)*Wr;
+            V(pos_y_max,pos_x_max,3) = V(pos_y_max,pos_x_max,3) + Scr;
+            Vt(pos_y_max,pos_x_max) = Vt(pos_y_max,pos_x_max) + Scr;
+            Wr = 0;
+        end
+
+        if((x == img_w) || (y == img_h))
+            done = true;
+        else
+            x = x + 1;
+            y = y + 1;
+        end
+    end
+end
+
+% Direction: /
+bla = k(:,:,4) > 0;
+for start=1:(img_w+img_h-1)
+    % Initial values
+    if (start <= img_w)
+        x = start;
+        y = img_h;
+    else
+        x = 1;
+        y = img_w+img_h-start;
+    end
+    done = false;
+
+    while ~done
+        if(bla(y,x))
+            Wr = Wr + 1;
+        end
+
+        if ( Wr > 0 && (y == 1 || x == img_w || ~bla(y,x)) )
+            if (y == 1 || x == img_w)
+                % Reached edge of image
+                pos_x_end = x;
+                pos_y_end = y;
+            else
+                pos_x_end = x - 1;
+                pos_y_end = y + 1;
+            end
+            pos_x_start = pos_x_end - Wr + 1;
+            pos_y_start = pos_y_end + Wr - 1;
+
+            %d = diag(flipud(k(pos_y_end:pos_y_start, pos_x_start:pos_x_end, 4)));
+            % More efficient implementation than diag(flipud(..))
+            d = k(((pos_x_start-1)*img_h + pos_y_start + 3*img_w*img_h):(img_h - 1):((pos_x_end-1)*img_h + pos_y_end + 3*img_w*img_h));
+            [~, I] = max(d);
+            pos_x_max = pos_x_start + I - 1;
+            pos_y_max = pos_y_start - I + 1;
+            Scr = k(pos_y_max,pos_x_max,4)*Wr;
+            V(pos_y_max,pos_x_max,4) = V(pos_y_max,pos_x_max,4) + Scr;
+            Vt(pos_y_max,pos_x_max) = Vt(pos_y_max,pos_x_max) + Scr;
+            Wr = 0;
+        end
+
+        if((x == img_w) || (y == 1))
+            done = true;
+        else
+            x = x + 1;
+            y = y - 1;
+        end
+    end
+end
+
+
+%Vt = V(:,:,1) + V(:,:,2) + V(:,:,3) + V(:,:,4);
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_v_matlab.hdf5'), '/V', V);
+end
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_vt_matlab.hdf5'), '/Vt', Vt);
+end
+
+%% Connection of vein centres
+Cd = zeros(img_h, img_w, 4);
+for x=3:img_w-3
+    for y=3:img_h-3
+        Cd(y,x,1) = min(max(Vt(y,x+1),  Vt(y,x+2))  ,...
+            max(Vt(y,x-1),  Vt(y,x-2)));   % Hor
+        Cd(y,x,2) = min(max(Vt(y+1,x),  Vt(y+2,x))  ,...
+            max(Vt(y-1,x),  Vt(y-2,x)));   % Vert
+        Cd(y,x,3) = min(max(Vt(y-1,x-1),Vt(y-2,x-2)),...
+            max(Vt(y+1,x+1),Vt(y+2,x+2))); % \
+        Cd(y,x,4) = min(max(Vt(y+1,x-1),Vt(y+2,x-2)),...
+            max(Vt(y-1,x+1),Vt(y-2,x+2))); % /
+    end
+end
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_cd_matlab.hdf5'), '/Cd', Cd);
+end
+
+veins = max(Cd,[],3);
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_g_matlab.hdf5'), '/G', veins);
+end
+
+md = median(veins(veins>0));
+veins_bin = veins > md;
+
+if nargin >= 4
+  hdf5write(strcat(stem, '_bin.hdf5'), '/binarised', uint8(veins_bin));
+end
+
+% %% Plot results
+% figure('Name', 'Second order derivatives');
+% subplot(2,2,1);
+%   imshow(fxx, []);
+%   title('Horizontal');
+% subplot(2,2,2);
+%   imshow(fyy, []);
+%   title('Vertical');
+% subplot(2,2,3);
+%   imshow(f11, []);
+%   title('\');
+% subplot(2,2,4);
+%   imshow(f22, []);
+%   title('/');
+%
+% figure('Name', 'Curvatures');
+% subplot(2,2,1);
+%   %imshow(log(k(:,:,1) + 1), []);
+%   imshow(k(:,:,1) > 0, []);
+%   title('Horizontal');
+% subplot(2,2,2);
+%   %imshow(log(k(:,:,2) + 1), []);
+%   imshow(k(:,:,2) > 0, []);
+%   title('Vertical');
+% subplot(2,2,3);
+%   %imshow(log(k(:,:,3) + 1), []);
+%   imshow(k(:,:,3) > 0, []);
+%   title('\');
+% subplot(2,2,4);
+%   %imshow(log(k(:,:,4) + 1), []);
+%   imshow(k(:,:,4) > 0, []);
+%   title('/');
+%
+% figure('Name', 'Scores');
+% subplot(2,2,1);
+%   imshow(V(:,:,1));
+%   title('Horizontal');
+% subplot(2,2,2);
+%   imshow(V(:,:,2));
+%   title('Vertical');
+% subplot(2,2,3);
+%   imshow(V(:,:,3));
+%   title('\');
+% subplot(2,2,4);
+%   imshow(V(:,:,3));
+%   title('/');

matlab/lib/miura_repeated_line_tracking.m

+function veins = miura_repeated_line_tracking(img, fvr, iterations, r, W)
+% Repeated line tracking
+
+% Parameters:
+%  img        - Input vascular image
+%  fvr        - Binary image of the finger region
+%  iterations - Maximum number of iterations
+%  r - Distance between tracking point and cross section of the profile
+%  W - Width of profile
+
+% Returns:
+%  veins - Vein image
+
+% Reference:
+% Feature extraction of finger vein patterns based on repeated line 
+%    tracking and its application to personal identification
+% N. Miura, A. Nagasaka, and T. Miyatake
+% Machine Vision and Applications, Volume 15, Number 4 (2004), pp. 194--203
+% doi: 10.1007/s00138-004-0149-2
+
+% Author:  Bram Ton <b.t.ton@alumnus.utwente.nl>
+% Date:    20th December 2011
+% License: Simplified BSD License
+
+p_lr = 0.5;  % Probability of goin left or right
+p_ud = 0.25; % Probability of going up or down
+% writerObj = VideoWriter('peaks.avi');
+% open(writerObj);
+Tr = zeros(size(img));     % Locus space
+bla = [-1,-1; -1,0; -1,1; 0,-1; 0,0; 0,1; 1,-1; 1,0; 1,1];
+
+% Check if W is even
+if (mod(W,2) == 0)
+    disp('Error: W must be odd')
+end
+ro = round(r*sqrt(2)/2);   % r for oblique directions
+hW = (W-1)/2;              % half width for horz. and vert. directions
+hWo = round(hW*sqrt(2)/2); % half width for oblique directions
+
+% Omit unreachable borders
+fvr(1:r+hW,:) = 0;
+fvr(end-(r+hW-1):end,:) = 0;
+fvr(:,1:r+hW) = 0;
+fvr(:,end-(r+hW-1):end) = 0;
+
+%% Uniformly distributed starting points
+indices = find(fvr > 0);
+a = randperm(length(indices));
+a = a(1:iterations); % Limit to number of iterations
+[ys, xs] = ind2sub(size(img), indices(a));
+%ys = [200;20];% LET OP
+%xs=  [200;200];% LET OP
+
+%% Iterate through all starting points
+for it = 1:size(ys,1)
+    xc = xs(it); % Current tracking point, x
+    yc = ys(it); % Current tracking point, y
+    
+    % Determine the moving-direction attributes
+    % Going left or right ?
+    if rand() >= 0.5
+        Dlr = -1;  % Going left
+    else
+        Dlr = 1; % Going right 
+    end
+
+    % Going up or down ?
+    if rand() >= 0.5
+        Dud = -1;  % Going up
+    else
+        Dud = 1; % Going down 
+    end
+    
+    % Initialize locus-positition table Tc
+    Tc = false(size(img));
+
+    %Dlr = -1; Dud=-1;% LET OP
+    Vl = 1;
+    while Vl > 0;
+        %% Determine the moving candidate point set Nc
+        Nr = false(3);
+        Rnd = rand();
+        %Rnd = 0.8;% LET OP
+        if Rnd < p_lr
+            % Going left or right
+            Nr(:,2+Dlr) = true;
+        elseif (Rnd >= p_lr) && (Rnd < p_lr + p_ud)  
+            % Going up or down
+            Nr(2+Dud,:) = true;
+        else
+            % Going any direction
+            Nr = true(3);
+            Nr(2,2) = false;
+        end
+
+        tmp = find( ~Tc(yc-1:yc+1,xc-1:xc+1) & Nr & fvr(yc-1:yc+1,xc-1:xc+1) );
+        Nc =[xc + bla(tmp,1), yc + bla(tmp,2)];
+        
+        if size(Nc,1)==0
+            Vl=-1;
+            continue
+        end
+
+        %% Detect dark line direction near current tracking point
+        Vdepths = zeros(size(Nc,1),1); % Valley depths
+        for i = 1:size(Nc,1)
+            % Horizontal or vertical 
+            if Nc(i,2) == yc
+                % Horizontal plane
+                yp = Nc(i,2);
+                if Nc(i,1) > xc
+                    % Right direction
+                    xp = Nc(i,1) + r;
+                else
+                    % Left direction
+                    xp = Nc(i,1) - r;
+                end
+
+                Vdepths(i) = img(yp + hW, xp) - ...
+                    2*img(yp,xp) + ...
+                    img(yp - hW, xp);
+            elseif Nc(i,1) == xc
+                % Vertical plane
+                xp = Nc(i,1);
+                if Nc(i,2) > yc
+                    % Down direction
+                    yp = Nc(i,2) + r;
+                else
+                    % Up direction
+                    yp = Nc(i,2) - r;
+                end
+
+                Vdepths(i) = img(yp, xp + hW) - ...
+                    2*img(yp,xp) + ...
+                    img(yp, xp - hW);
+            end
+            
+            % Oblique directions
+            if ((Nc(i,1) > xc) && (Nc(i,2) < yc)) || ((Nc(i,1) < xc) && (Nc(i,2) > yc))
+                % Diagonal, up /
+                if Nc(i,1) > xc && Nc(i,2) < yc
+                    % Top right
+                    xp = Nc(i,1) + ro;
+                    yp = Nc(i,2) - ro;
+                else
+                    % Bottom left
+                    xp = Nc(i,1) - ro;
+                    yp = Nc(i,2) + ro;
+                end
+
+                Vdepths(i) = img(yp - hWo, xp - hWo) - ...
+                    2*img(yp,xp) + ...
+                    img(yp + hWo, xp + hWo);
+            else
+                % Diagonal, down \
+                if Nc(i,1) < xc && Nc(i,2) < yc
+                    % Top left
+                    xp = Nc(i,1) - ro;
+                    yp = Nc(i,2) - ro;
+                else
+                    % Bottom right
+                    xp = Nc(i,1) + ro;
+                    yp = Nc(i,2) + ro;
+                end
+                
+                Vdepths(i) = img(yp + hWo, xp - hWo) - ...
+                    2*img(yp,xp) + ...
+                    img(yp - hWo, xp + hWo);
+            end
+        end %  End search of candidates
+
+        [~, index] = max(Vdepths); %  Determine best candidate
+
+        % Register tracking information
+        Tc(yc, xc) = true;
+
+        % Increase value of tracking space
+        Tr(yc, xc) = Tr(yc, xc) + 1;
+        %writeVideo(writerObj,Tr);
+        
+        % Move tracking point
+        xc = Nc(index, 1);
+        yc = Nc(index, 2); 
+    end
+end
+veins = Tr;
\ No newline at end of file

matlab/lib/miura_usage.m

+% Howto use the miura_* scripts.
+
+img = im2double(imread('finger.png')); % Read the image
+img = imresize(img,0.5);               % Downscale image
+
+% Get the valid region, this is a binary mask which indicates the region of 
+% the finger. For quick testing it is possible to use something like:
+% fvr = ones(size(img));
+% The lee_region() function can be found here:
+% http://www.mathworks.com/matlabcentral/fileexchange/35752-finger-region-localisation
+fvr = lee_region(img,4,40);    % Get finger region
+
+%% Extract veins using maximum curvature method
+sigma = 3; % Parameter
+v_max_curvature = miura_max_curvature(img,fvr,sigma);
+
+% Binarise the vein image
+md = median(v_max_curvature(v_max_curvature>0));
+v_max_curvature_bin = v_max_curvature > md; 
+
+%% Extract veins using repeated line tracking method
+max_iterations = 3000; r=1; W=17; % Parameters
+v_repeated_line = miura_repeated_line_tracking(img,fvr,max_iterations,r,W);
+
+% Binarise the vein image
+md = median(v_repeated_line(v_repeated_line>0));
+v_repeated_line_bin = v_repeated_line > md; 
+
+%% Match
+cw = 80; ch=30;
+% Note that the match score is between 0 and 0.5
+score = miura_match(double(v_repeated_line_bin), double(v_max_curvature_bin), cw, ch);
+fprintf('Match score: %6.4f %%\n', score);
+
+%% Visualise
+% Overlay the extracted veins on the original image
+overlay_max_curvature = zeros([size(img) 3]);
+overlay_max_curvature(:,:,1) = img;
+overlay_max_curvature(:,:,2) = img + 0.4*v_max_curvature_bin;
+overlay_max_curvature(:,:,3) = img;
+
+% Overlay the extracted veins on the original image
+overlay_repeated_line = zeros([size(img) 3]);
+overlay_repeated_line(:,:,1) = img;
+overlay_repeated_line(:,:,2) = img + 0.4*v_repeated_line_bin;
+overlay_repeated_line(:,:,3) = img;
+
+figure;
+subplot(3,2,1)
+  imshow(img,[])
+  title('Original captured image')
+subplot(3,2,2)
+  imshow(fvr)
+  title('Detected finger region')
+subplot(3,2,3)
+  imshow(v_max_curvature_bin)
+  title('Binarised veins extracted by maximum curvature method')
+subplot(3,2,4)
+  imshow(overlay_max_curvature)
+  title('Maximum curvature method')  
+subplot(3,2,5)
+  imshow(v_repeated_line_bin)
+  title('Binarised veins extracted by repeated line tracking method')
+subplot(3,2,6)
+  imshow(overlay_repeated_line)
+  title('Repeated line tracking method')
\ No newline at end of file

matlab/run.sh

+#!/usr/bin/env bash
+
+if [ $# == 0 ]; then
+  echo "usage: $0 input_image.mat input_region.mat output_stem"
+  exit 1
+fi
+
+_matlab=`which matlab`;
+
+if [ -z "${_matlab}" ]; then
+  echo "I cannot find a matlab binary to execute, please setup first"
+  exit 1
+fi
+
+# Does some environment manipulation to get paths right before we start
+_basedir="$(cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd)"
+_matlabpath=${_basedir}/lib
+unset _basedir;
+if [ -n "${MATLABPATH}" ]; then
+  _matlabpath=${_matlabpath}:${MATLABPATH}
+fi
+export MATLABPATH=${_matlabpath};
+unset _matlabpath;
+
+# Calls matlab with our inputs
+${_matlab} -nodisplay -nosplash -nodesktop -r "image = im2double(hdf5read('${1}','/array')); region = double(hdf5read('${2}','/array')); miura_max_curvature(image, region, 3, '${3}'); quit;"