diff --git a/matlab/README.md b/matlab/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..b2135fe6927006832f3ce00b91d7c7d69d703067
--- /dev/null
+++ b/matlab/README.md
@@ -0,0 +1,69 @@
+# 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.
diff --git a/matlab/compare.py b/matlab/compare.py
new file mode 100644
index 0000000000000000000000000000000000000000..f78bfad4e513da2d96309c49c13238c235b0f63c
--- /dev/null
+++ b/matlab/compare.py
@@ -0,0 +1,54 @@
+#!/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())
diff --git a/matlab/lib/license.txt b/matlab/lib/license.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0225a2b84323c48b07c435b3fda6abc2b533eba4
--- /dev/null
+++ b/matlab/lib/license.txt
@@ -0,0 +1,24 @@
+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.
diff --git a/matlab/lib/miura_match.m b/matlab/lib/miura_match.m
new file mode 100644
index 0000000000000000000000000000000000000000..ccc040659f30e576eff0a3d7a8e8f19558ce9be1
--- /dev/null
+++ b/matlab/lib/miura_match.m
@@ -0,0 +1,52 @@
+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)));
+
+
diff --git a/matlab/lib/miura_max_curvature.m b/matlab/lib/miura_max_curvature.m
new file mode 100644
index 0000000000000000000000000000000000000000..ee83ee23c59eeb1f54cdae2c0b026ab7a52b6ff3
--- /dev/null
+++ b/matlab/lib/miura_max_curvature.m
@@ -0,0 +1,311 @@
+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('/');
diff --git a/matlab/lib/miura_repeated_line_tracking.m b/matlab/lib/miura_repeated_line_tracking.m
new file mode 100644
index 0000000000000000000000000000000000000000..08b3eb55786c9686ef5e4a339fc2c6b16baf12a7
--- /dev/null
+++ b/matlab/lib/miura_repeated_line_tracking.m
@@ -0,0 +1,186 @@
+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
diff --git a/matlab/lib/miura_usage.m b/matlab/lib/miura_usage.m
new file mode 100644
index 0000000000000000000000000000000000000000..51c815c993dddf069ee587fa24529cceb7fd7f64
--- /dev/null
+++ b/matlab/lib/miura_usage.m
@@ -0,0 +1,66 @@
+% 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
diff --git a/matlab/run.sh b/matlab/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b95ff2dd8e7cb6b8ce0d8d0883157d0e5d0842bd
--- /dev/null
+++ b/matlab/run.sh
@@ -0,0 +1,26 @@
+#!/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;"