Documentation>vl_demo_kdtree_sift

# Source code for ``vl_demo_kdtree_sift.m

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