Source code for
This file is located in the toolbox/demo folder in VLFeat package.
function vl_demo_kdtree_sift
% VL_DEMO_KDTREE_SIFT
% Demonstrates the use of a kd-tree forest to match SIFT
% features. If FLANN is present, this function runs a comparison
% against it.
% AUTORIGHS
rand('state',0) ;
randn('state',0);
do_median = 0 ;
do_mean = 1 ;
% try to setup flann
if ~exist('flann_search', 'file')
if exist(fullfile(vl_root, 'opt', 'flann', 'build', 'matlab'))
addpath(fullfile(vl_root, 'opt', 'flann', 'build', 'matlab')) ;
end
end
do_flann = exist('nearest_neighbors') == 3 ;
if ~do_flann
warning('FLANN not found. Comparison disabled.') ;
end
maxNumComparisonsRange = [1 10 50 100 200 300 400] ;
numTreesRange = [1 2 5 10] ;
% get data (SIFT features)
im1 = imread(fullfile(vl_root, 'data', 'roofs1.jpg')) ;
im2 = imread(fullfile(vl_root, 'data', 'roofs2.jpg')) ;
im1 = single(rgb2gray(im1)) ;
im2 = single(rgb2gray(im2)) ;
[f1,d1] = vl_sift(im1,'firstoctave',-1,'floatdescriptors','verbose') ;
[f2,d2] = vl_sift(im2,'firstoctave',-1,'floatdescriptors','verbose') ;
% add some noise to make matches unique
d1 = single(d1) + rand(size(d1)) ;
d2 = single(d2) + rand(size(d2)) ;
% match exhaustively to get the ground truth
elapsedDirect = tic ;
D = vl_alldist(d1,d2) ;
[drop, best] = min(D, [], 1) ;
elapsedDirect = toc(elapsedDirect) ;
for ti=1:length(numTreesRange)
for vi=1:length(maxNumComparisonsRange)
v = maxNumComparisonsRange(vi) ;
t = numTreesRange(ti) ;
if do_median
tic ;
kdtree = vl_kdtreebuild(d1, ...
'verbose', ...
'thresholdmethod', 'median', ...
'numtrees', t) ;
[i, d] = vl_kdtreequery(kdtree, d1, d2, ...
'verbose', ...
'maxcomparisons',v) ;
elapsedKD_median(vi,ti) = toc ;
errors_median(vi,ti) = sum(double(i) ~= best) / length(best) ;
errorsD_median(vi,ti) = mean(abs(d - drop) ./ drop) ;
end
if do_mean
tic ;
kdtree = vl_kdtreebuild(d1, ...
'verbose', ...
'thresholdmethod', 'mean', ...
'numtrees', t) ;
%kdtree = readflann(kdtree, '/tmp/flann.txt') ;
%checkx(kdtree, d1, 1, 1) ;
[i, d] = vl_kdtreequery(kdtree, d1, d2, ...
'verbose', ...
'maxcomparisons', v) ;
elapsedKD_mean(vi,ti) = toc ;
errors_mean(vi,ti) = sum(double(i) ~= best) / length(best) ;
errorsD_mean(vi,ti) = mean(abs(d - drop) ./ drop) ;
end
if do_flann
tic ;
[i, d] = flann_search(d1, d2, 1, struct('algorithm','kdtree', ...
'trees', t, ...
'checks', v));
ifla = i ;
elapsedKD_flann(vi,ti) = toc;
errors_flann(vi,ti) = sum(i ~= best) / length(best) ;
errorsD_flann(vi,ti) = mean(abs(d - drop) ./ drop) ;
end
end
end
figure(1) ; clf ;
leg = {} ;
hnd = [] ;
sty = {{'color','r'},{'color','g'},...
{'color','b'},{'color','c'},...
{'color','k'}} ;
for ti=1:length(numTreesRange)
s = sty{mod(ti,length(sty))+1} ;
if do_median
h1=loglog(elapsedDirect ./ elapsedKD_median(:,ti),100*errors_median(:,ti),'-*',s{:}) ; hold on ;
leg{end+1} = sprintf('VLFeat median (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h1 ;
end
if do_mean
h2=loglog(elapsedDirect ./ elapsedKD_mean(:,ti), 100*errors_mean(:,ti), '-o',s{:}) ; hold on ;
leg{end+1} = sprintf('VLFeat (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h2 ;
end
if do_flann
h3=loglog(elapsedDirect ./ elapsedKD_flann(:,ti), 100*errors_flann(:,ti), '+--',s{:}) ; hold on ;
leg{end+1} = sprintf('FLANN (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h3 ;
end
end
set([hnd], 'linewidth', 2) ;
xlabel('speedup over linear search (log times)') ;
ylabel('percentage of incorrect matches (%)') ;
h=legend(hnd, leg{:}, 'location', 'southeast') ;
set(h,'fontsize',8) ;
grid on ;
axis square ;
vl_demo_print('kdtree_sift_incorrect',.6) ;
figure(2) ; clf ;
leg = {} ;
hnd = [] ;
for ti=1:length(numTreesRange)
s = sty{mod(ti,length(sty))+1} ;
if do_median
h1=loglog(elapsedDirect ./ elapsedKD_median(:,ti),100*errorsD_median(:,ti),'*-',s{:}) ; hold on ;
leg{end+1} = sprintf('VLFeat median (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h1 ;
end
if do_mean
h2=loglog(elapsedDirect ./ elapsedKD_mean(:,ti), 100*errorsD_mean(:,ti), 'o-',s{:}) ; hold on ;
leg{end+1} = sprintf('VLFeat (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h2 ;
end
if do_flann
h3=loglog(elapsedDirect ./ elapsedKD_flann(:,ti), 100*errorsD_flann(:,ti), '+--',s{:}) ; hold on ;
leg{end+1} = sprintf('FLANN (%d tr.)', numTreesRange(ti)) ;
hnd(end+1) = h3 ;
end
end
set([hnd], 'linewidth', 2) ;
xlabel('speedup over linear search (log times)') ;
ylabel('relative overestimation of minmium distannce (%)') ;
h=legend(hnd, leg{:}, 'location', 'southeast') ;
set(h,'fontsize',8) ;
grid on ;
axis square ;
vl_demo_print('kdtree_sift_distortion',.6) ;
% --------------------------------------------------------------------
function checkx(kdtree, X, t, n, mib, mab)
% --------------------------------------------------------------------
if nargin <= 4
mib = -inf * ones(size(X,1),1) ;
mab = +inf * ones(size(X,1),1) ;
end
lc = kdtree.trees(t).nodes.lowerChild(n) ;
uc = kdtree.trees(t).nodes.upperChild(n) ;
if lc < 0
for i=-lc:-uc-1
di = kdtree.trees(t).dataIndex(i) ;
if any(X(:,di) > mab)
error('a') ;
end
if any(X(:,di) < mib)
error('b') ;
end
end
return
end
i = kdtree.trees(t).nodes.splitDimension(n) ;
v = kdtree.trees(t).nodes.splitThreshold(n) ;
mab_ = mab ;
mab_(i) = min(mab(i), v) ;
checkx(kdtree, X, t, lc, mib, mab_) ;
mib_ = mib ;
mib_(i) = max(mib(i), v) ;
checkx(kdtree, X, t, uc, mib_, mab) ;
% --------------------------------------------------------------------
function kdtree = readflann(kdtree, path)
% --------------------------------------------------------------------
data = textread(path)' ;
for i=1:size(data,2)
nodeIds = data(1,:) ;
ni = find(nodeIds == data(1,i)) ;
if ~isnan(data(2,i))
% internal node
li = find(nodeIds == data(4,i)) ;
ri = find(nodeIds == data(5,i)) ;
kdtree.trees(1).nodes.lowerChild(ni) = int32(li) ;
kdtree.trees(1).nodes.upperChild(ni) = int32(ri) ;
kdtree.trees(1).nodes.splitThreshold(ni) = single(data(2,i)) ;
kdtree.trees(1).nodes.splitDimension(ni) = single(data(3,i)+1) ;
else
di = data(3,i) + 1 ;
kdtree.trees(1).nodes.lowerChild(ni) = int32(- di) ;
kdtree.trees(1).nodes.upperChild(ni) = int32(- di - 1) ;
end
kdtree.trees(1).dataIndex = uint32(1:kdtree.numData) ;
end