initial sd-vbs

initial sd-vbs

add sd-vbs

sd-vbs

41 files changed, 4120 insertions, 0 deletions

diff --git a/SD-VBS/benchmarks/sift/src/matlab/diffss.m b/SD-VBS/benchmarks/sift/src/matlab/diffss.m
new file mode 100644
index 0000000..1e65390
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/diffss.m
@@ -0,0 +1,72 @@
+function dss = diffss(ss)
+% DIFFSS  Difference of scale space
+%   DSS=DIFFSS(SS) returns a scale space DSS obtained by subtracting
+%   consecutive levels of the scale space SS.
+%
+%   In SIFT, this function is used to compute the difference of
+%   Gaussian scale space from the Gaussian scale space of an image.
+%
+%   See also GAUSSIANSS(), PDF:SIFT.USER.SS.
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+dss.smin = ss.smin ;
+dss.smax = ss.smax-1 ;
+dss.omin = ss.omin ;
+dss.O = ss.O ;
+dss.S = ss.S ;
+dss.sigma0 = ss.sigma0 ;
+
+for o=1:dss.O
+  % Can be done at once, but it turns out to be faster
+  % this way
+  [M,N,S] = size(ss.octave{o}) ;
+  dss.octave{o} = zeros(M,N,S-1) ;
+  for s=1:S-1
+    dss.octave{o}(:,:,s) = ...
+        ss.octave{o}(:,:,s+1) - ss.octave{o}(:,:,s) ;
+  end
+end
+
+%    for i=1:32
+%        for j=1:32
+%            fprintf(1, '%f\n', dss.octave{1}(i,j,1));
+%        end
+%    end
diff --git a/SD-VBS/benchmarks/sift/src/matlab/gaussianss.m b/SD-VBS/benchmarks/sift/src/matlab/gaussianss.m
new file mode 100644
index 0000000..c92383c
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/gaussianss.m
@@ -0,0 +1,216 @@
+function SS = gaussianss(I,sigman,O,S,omin,smin,smax,sigma0)
+% GAUSSIANSS
+%   SS = GAUSSIANSS(I,SIGMAN,O,S,OMIN,SMIN,SMAX,SIGMA0) returns the
+%   Gaussian scale space of image I. Image I is assumed to be
+%   pre-smoothed at level SIGMAN. O,S,OMIN,SMIN,SMAX,SIGMA0 are the
+%   parameters of the scale space as explained in PDF:SIFT.USER.SS.
+%
+%   See also DIFFSS(), PDF:SIFT.USER.SS.
+
+% History
+%  4-15-2006  Fixed some comments
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+% --------------------------------------------------------------------
+%                                                  Check the arguments
+% --------------------------------------------------------------------
+% --------------------------------------------------------------------
+%                                                           Do the job
+% --------------------------------------------------------------------
+
+% Scale multiplicative step
+k = 2^(1/S) ;
+
+dsigma0 = sigma0 * sqrt(1 - 1/k^2) ; % Scale step factor
+sigman  = 0.5 ;                      % Nominal smoothing of the image
+
+% Scale space structure
+SS.O          = O ;
+SS.S          = S ;
+SS.sigma0     = sigma0 ;
+SS.omin       = omin ;
+SS.smin       = smin ;
+SS.smax       = smax ;
+
+% If mino < 0, multiply the size of the image.
+% (The rest of the code is consistent with this.)
+
+
+if omin < 0
+        for o=1:-omin
+                I = doubleSize(I) ;
+        end
+elseif omin > 0
+        for o=1:omin
+                I = halveSize(I) ;
+        end
+end
+
+[M,N] = size(I) ;
+
+% Index offset
+so = -smin+1 ;
+
+% --------------------------------------------------------------------
+%                                                         First octave
+% --------------------------------------------------------------------
+%
+% The first level of the first octave has scale index (o,s) =
+% (omin,smin) and scale coordinate
+%
+%    sigma(omin,smin) = sigma0 2^omin k^smin 
+%
+% The input image I is at nominal scale sigman. Thus in order to get
+% the first level of the pyramid we need to apply a smoothing of
+%  
+%   sqrt( (sigma0 2^omin k^smin)^2 - sigman^2 ).
+%
+% As we have pre-scaled the image omin octaves (up or down,
+% depending on the sign of omin), we need to correct this value
+% by dividing by 2^omin, getting
+%e
+%   sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 )
+%
+
+SS.octave{1} = zeros(M,N,smax-smin+1) ;
+SS.octave{1}(:,:,1)  = imsmooth(I, ...
+        sqrt((sigma0*k^smin)^2  - (sigman/2^omin)^2)) ;
+
+temp = sqrt((sigma0*k^smin)^2  - (sigman/2^omin)^2) ;
+
+for s=smin+1:smax
+        % Here we go from (omin,s-1) to (omin,s). The extra smoothing
+        % standard deviation is
+        %
+        %  (sigma0 2^omin 2^(s/S) )^2 - (simga0 2^omin 2^(s/S-1/S) )^2
+        %
+        % Aftred dividing by 2^omin (to take into account the fact
+  % that the image has been pre-scaled omin octaves), the 
+  % standard deviation of the smoothing kernel is
+  %
+        %   dsigma = sigma0 k^s sqrt(1-1/k^2)
+  %
+        dsigma = k^s * dsigma0 ;
+        SS.octave{1}(:,:,s +so) = ...
+                imsmooth(squeeze(...
+                SS.octave{1}(:,:,s-1 +so)...
+                ), dsigma )  ;
+end
+
+% --------------------------------------------------------------------
+%                                                        Other octaves
+% --------------------------------------------------------------------
+
+for o=2:O  
+        % We need to initialize the first level of octave (o,smin) from
+        % the closest possible level of the previous octave. A level (o,s)
+  % in this octave corrsponds to the level (o-1,s+S) in the previous
+  % octave. In particular, the level (o,smin) correspnds to
+  % (o-1,smin+S). However (o-1,smin+S) might not be among the levels
+  % (o-1,smin), ..., (o-1,smax) that we have previously computed.
+  % The closest pick is
+  %
+        %                       /  smin+S    if smin+S <= smax
+        % (o-1,sbest) , sbest = |
+        %                       \  smax      if smin+S > smax
+        %
+        % The amount of extra smoothing we need to apply is then given by
+        %
+        %  ( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(sbest/S - 1) )^2
+        %
+  % As usual, we divide by 2^o to cancel out the effect of the
+  % downsampling and we get
+  %
+        %  ( sigma 0 k^smin )^2 - ( sigma0 2^o k^(sbest - S) )^2
+  %
+        sbest = min(smin + S, smax) ;
+        TMP = halveSize(squeeze(SS.octave{o-1}(:,:,sbest+so))) ;
+
+        target_sigma = sigma0 * k^smin ;
+          prev_sigma = sigma0 * k^(sbest - S) ;
+          
+        if(target_sigma > prev_sigma)
+          temp = sqrt(target_sigma^2 - prev_sigma^2);
+          TMP = imsmooth(TMP, temp ) ;                               
+        end
+        [M,N] = size(TMP) ;
+        
+        SS.octave{o} = zeros(M,N,smax-smin+1) ;
+        SS.octave{o}(:,:,1) = TMP ;
+
+        for s=smin+1:smax
+                % The other levels are determined as above for the first octave.                
+                dsigma = k^s * dsigma0 ;
+                SS.octave{o}(:,:,s +so) = ...
+                        imsmooth(squeeze(...
+                        SS.octave{o}(:,:,s-1 +so)...
+                        ), dsigma)  ;    
+        end
+    
+        
+end
+
+% -------------------------------------------------------------------------
+%                                                       Auxiliary functions
+% -------------------------------------------------------------------------
+function J = doubleSize(I)
+[M,N]=size(I) ;
+J = zeros(2*M,2*N) ;
+J(1:2:end,1:2:end) = I ;
+J(2:2:end-1,2:2:end-1) = ...
+        0.25*I(1:end-1,1:end-1) + ...
+        0.25*I(2:end,1:end-1) + ...
+        0.25*I(1:end-1,2:end) + ...
+        0.25*I(2:end,2:end) ;
+J(2:2:end-1,1:2:end) = ...
+        0.5*I(1:end-1,:) + ...
+    0.5*I(2:end,:) ;
+J(1:2:end,2:2:end-1) = ...
+        0.5*I(:,1:end-1) + ...
+    0.5*I(:,2:end) ;
+
+function J = halveSize(I)
+J=I(1:2:end,1:2:end) ;
+%[M,N] = size(I) ;
+%m=floor((M+1)/2) ;
+%n=floor((N+1)/2) ;
+%J = I(:,1:2:2*n) + I(:,2:2:2*n+1) ;
+%J = 0.25*(J(1:2:2*m,:)+J(2:2:2*m+1,:)) ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m b/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m
new file mode 100644
index 0000000..55cb708
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m
@@ -0,0 +1,52 @@
+function I = imreadbw(file) 
+% IMREADBW  Reads an image as gray-scale
+%   I=IMREADBW(FILE) reads the image FILE and converts the result to a
+%   gray scale image (with DOUBLE storage class anr range normalized
+%   in [0,1]).
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+I=im2double(imread(file)) ;
+
+if(size(I,3) > 1)
+  I = rgb2gray( I ) ;
+end
+
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c
new file mode 100644
index 0000000..5d6a660
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c
@@ -0,0 +1,115 @@
+/** file:        imsmooth.c
+ ** author:      Andrea Vedaldi
+ ** description: Smooth an image.
+ **/
+
+#include"mex.h"
+#include<stdlib.h>
+#include<string.h>
+#include<math.h>
+#include<assert.h>
+
+#define greater(a,b) ((a) > (b))
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+const double win_factor = 1.5 ;
+const int nbins = 36 ;
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int M,N ;
+  double* I_pt ;
+  double* J_pt ;
+  double s ;
+  enum {I=0,S} ;
+  enum {J=0} ;
+
+  /* ------------------------------------------------------------------
+  **                                                Check the arguments
+  ** --------------------------------------------------------------- */ 
+  if (nin != 2) {
+    mexErrMsgTxt("Exactly two input arguments required.");
+  } else if (nout > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  if (!mxIsDouble(in[I]) || 
+      !mxIsDouble(in[S]) ) {
+    mexErrMsgTxt("All arguments must be real.") ;
+  }
+  
+  if(mxGetNumberOfDimensions(in[I]) > 2||
+     mxGetNumberOfDimensions(in[S]) > 2) {
+    mexErrMsgTxt("I must be a two dimensional array and S a scalar.") ;
+  }
+  
+  if(max(mxGetM(in[S]),mxGetN(in[S])) > 1) {
+    mexErrMsgTxt("S must be a scalar.\n") ;
+  }
+
+  M = mxGetM(in[I]) ;
+  N = mxGetN(in[I]) ;
+
+  out[J] = mxCreateDoubleMatrix(M, N, mxREAL) ;
+  
+  I_pt   = mxGetPr(in[I]) ;
+  J_pt   = mxGetPr(out[J]) ;
+  s      = *mxGetPr(in[S]) ;
+
+  /* ------------------------------------------------------------------
+  **                                                         Do the job
+  ** --------------------------------------------------------------- */ 
+  if(s > 0.01) {
+    
+    int W = (int) ceil(4*s) ;
+    int i ; 
+    int j ;
+    double* g0 = (double*) mxMalloc( (2*W+1)*sizeof(double) ) ;
+    double* buffer = (double*) mxMalloc( M*N*sizeof(double) ) ;
+    double acc=0.0 ;
+    
+    for(j = 0 ; j < 2*W+1 ; ++j) {
+      g0[j] = exp(-0.5 * (j - W)*(j - W)/(s*s)) ;
+      acc += g0[j] ;
+    }
+    for(j = 0 ; j < 2*W+1 ; ++j) {
+      g0[j] /= acc ;
+    }
+    
+    /*
+    ** Convolve along the columns
+    **/
+    for(j = 0 ; j < N ; ++j) {
+      for(i = 0 ; i < M ; ++i) {
+        double* start = I_pt + max(i-W,0) + j*M ;        
+        double* stop = I_pt + min(i+W,M-1) + j*M + 1 ;
+        double* g = g0 + max(0, W-i) ;
+        acc = 0.0 ;
+        while(stop != start) acc += (*g++) * (*start++) ;
+        *buffer++ = acc ;
+      }
+    }
+    buffer -= M*N ;
+    
+    /*
+    ** Convolve along the rows
+    **/
+    for(j = 0 ; j < N ; ++j) {
+      for(i = 0 ; i < M ; ++i) {
+        double* start = buffer + i + max(j-W,0)*M ;
+        double* stop  = buffer + i + min(j+W,N-1)*M + M ;
+        double* g = g0 + max(0, W-j) ;
+        acc = 0.0 ;
+        while(stop != start) { acc += (*g++) * (*start) ; start+=M ;}
+        *J_pt++ = acc ;
+      }
+    }    
+    mxFree(buffer) ;
+    mxFree(g0) ;
+  } else {
+    memcpy(J_pt, I_pt, sizeof(double)*M*N) ;   
+  }
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64
new file mode 100755
index 0000000..0997262
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx
new file mode 100755
index 0000000..a4a7a20
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c b/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c
new file mode 100644
index 0000000..e0adeff
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c
@@ -0,0 +1,103 @@
+#include<stdlib.h>
+#include<string.h>
+#include<math.h>
+#include "imsmooth.h"
+
+#define greater(a,b) ((a) > (b))
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+const double win_factor = 1.5 ;
+const int nbins = 36 ;
+
+float* imsmooth(float* I_pt_, float dsigma)
+{
+  int M,N ;
+  float *I_pt ;
+  float* J_pt_, *J_pt ;
+  float* out_, *out;
+  float s ;
+  enum {I=0,S} ;
+  enum {J=0} ;
+
+  /* ------------------------------------------------------------------
+  **                                                Check the arguments
+  ** --------------------------------------------------------------- */ 
+
+  M = (int)in_[0];
+  N = (int)in_[1];
+
+  out_ = fMallocHandle(M, N);
+    J_pt_ = &(out_[0]);
+    J_pt = &(J_pt_[4]);
+    s = dsigma;
+    
+
+  /* ------------------------------------------------------------------
+  **                                                         Do the job
+  ** --------------------------------------------------------------- */ 
+  if(s > 0.01) 
+  {
+    
+    int W = (int) ceil(4*s) ;
+    int i ; 
+    int j ;
+    float* g0_, *g0, *buffer_, *buffer, acc;
+    g0_ = fMallocHandle(1, 2*W+1);
+    g0 = &(g0_[4]);
+    buffer_ = fMallocHandle(M, N);
+    acc=0.0 ;
+    
+    for(j = 0 ; j < 2*W+1 ; ++j) 
+    {
+      g0[j] = exp(-0.5 * (j - W)*(j - W)/(s*s)) ;
+      acc += g0[j] ;
+    }
+
+    for(j = 0 ; j < 2*W+1 ; ++j) 
+    {
+      g0[j] /= acc ;
+    }
+    
+    /*
+    ** Convolve along the columns
+    **/
+
+    for(j = 0 ; j < N ; ++j) 
+    {
+      for(i = 0 ; i < M ; ++i) 
+      {
+        float* start = I_pt + max(i-W,0) + j*M ;        
+        float* stop = I_pt + min(i+W,M-1) + j*M + 1 ;
+        float* g = g0 + max(0, W-i) ;
+        acc = 0.0 ;
+        while(stop != start) acc += (*g++) * (*start++) ;
+        *buffer++ = acc ;
+      }
+    }
+    buffer -= M*N ;
+    
+    /*
+    ** Convolve along the rows
+    **/
+    for(j = 0 ; j < N ; ++j) 
+    {
+      for(i = 0 ; i < M ; ++i) 
+      {
+        float* start = buffer + i + max(j-W,0)*M ;
+        float* stop  = buffer + i + min(j+W,N-1)*M + M ;
+        float* g = g0 + max(0, W-j) ;
+        acc = 0.0 ;
+        while(stop != start) { acc += (*g++) * (*start) ; start+=M ;}
+        *J_pt++ = acc ;
+      }
+    }    
+    free(buffer) ;
+    free(g0) ;
+  } 
+  else 
+  {
+    memcpy(J_pt, I_pt, sizeof(double)*M*N) ;   
+  }
+  return J_pt_;
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/mexutils.c b/SD-VBS/benchmarks/sift/src/matlab/mexutils.c
new file mode 100644
index 0000000..45843cf
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/mexutils.c
@@ -0,0 +1,98 @@
+/* file:        mexutils.c 
+** author:      Andrea Vedaldi
+** description: Utility functions to write MEX files.
+**/
+
+#include"mex.h"
+#include<math.h>
+
+#undef M_PI
+#define M_PI 3.14159265358979
+
+/** @biref Is real scalar?
+ **
+ ** @return @c true if the array @a A is a real scalar.
+ **/
+int
+uIsRealScalar(const mxArray* A)
+{
+  return 
+    mxIsDouble(A) && 
+    !mxIsComplex(A) &&
+    mxGetNumberOfDimensions(A) == 2 &&
+    mxGetM(A) == 1 &&
+    mxGetN(A) == 1 ;
+}
+
+/** @brief Is real matrix?
+ **
+ ** The function checks wether the argument @a A is a real matrix.  In
+ ** addition, if @a M >= 0, it checks wether the number of rows is
+ ** equal to @a M and, if @a N >= 0, if the number of columns is equal
+ ** to @a N.
+ **
+ ** @param M number of rows.
+ ** @param N number of columns.
+ ** @return @c true if the array is a real matrix with the specified format.
+ **/
+int
+uIsRealMatrix(const mxArray* A, int M, int N)
+{
+  return  
+    mxIsDouble(A) &&
+    !mxIsComplex(A) &&
+    mxGetNumberOfDimensions(A) == 2 &&
+    ((M>=0)?(mxGetM(A) == M):1) &&
+    ((N>=0)?(mxGetN(A) == N):1) ;   
+}
+
+/** @brief Is real vector?
+ **
+ ** The function checks wether the argument  @a V is a real vector. By
+ ** definiton, a  matrix is a vector  if one of its  dimension is one.
+ ** In addition, if  @a D >= 0, it checks wether  the dimension of the
+ ** vecotr is equal to @a D.
+ **
+ ** @param D lenght of the vector.
+ ** @return @c true if the array is a real vector of the specified dimension.
+ **/
+int
+uIsRealVector(const mxArray* V, int D) 
+{
+  int M = mxGetM(V) ;
+  int N = mxGetN(V) ;
+  int is_vector = (N == 1) || (M == 1) ;
+  
+  return   
+    mxIsDouble(V) &&
+    !mxIsComplex(V) &&
+    mxGetNumberOfDimensions(V) == 2 &&
+    is_vector &&
+    ( D < 0 || N == D || M == D) ;
+}
+
+
+/** @brief Is a string?
+ **
+ ** The function checks wether the array @a S is a string. If
+ ** @a L is non-negative, it also check wether the strign has
+ ** length @a L.
+ **
+ ** @return @a c true if S is a string of the specified length.
+ **/
+int
+uIsString(const mxArray* S, int L)
+{
+  int M = mxGetM(S) ;
+  int N = mxGetN(S) ;
+
+  return 
+    mxIsChar(S) &&
+    M == 1 &&
+    (L < 0 || N == L) ;
+}
+
+/**
+ **
+ **/
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m b/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m
new file mode 100644
index 0000000..b6324e9
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m
@@ -0,0 +1,124 @@
+function h=plotmatches(I1,I2,P1,P2,matches,varargin)
+% PLOTMATCHES  Plot keypoint matches
+%   PLOTMATCHES(I1,I2,P1,P2,MATCHES) plots the two images I1 and I2
+%   and lines connecting the frames (keypoints) P1 and P2 as specified
+%   by MATCHES.
+%
+%   P1 and P2 specify two sets of frames, one per column. The first
+%   two elements of each column specify the X,Y coordinates of the
+%   corresponding frame. Any other element is ignored.
+%
+%   MATCHES specifies a set of matches, one per column. The two
+%   elementes of each column are two indexes in the sets P1 and P2
+%   respectively.
+%
+%   The images I1 and I2 might be either both grayscale or both color
+%   and must have DOUBLE storage class. If they are color the range
+%   must be normalized in [0,1].
+%
+%   The function accepts the following option-value pairs:
+%
+%   'Stacking' ['h']
+%      Stacking of images: horizontal [h], vertical [v], diagonal
+%      [h], overlap ['o']
+%
+%   See also PLOTSIFTDESCRIPTOR(), PLOTSIFTFRAME(), PLOTSS().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+% --------------------------------------------------------------------
+%                                                  Check the arguments
+% --------------------------------------------------------------------
+
+stack='h' ;
+
+for k=1:2:length(varargin)
+  switch varargin{k}
+    case 'Stacking'
+      stack=varargin{k+1} ;
+  end
+end
+ 
+% --------------------------------------------------------------------
+%                                                           Do the job
+% --------------------------------------------------------------------
+
+[M1,N1,K1]=size(I1) ;
+[M2,N2,K2]=size(I2) ;
+
+switch stack
+  case 'h'
+    N3=N1+N2 ;
+    M3=max(M1,M2) ;
+    oj=N1 ;
+    oi=0 ;
+  case 'v'
+    M3=M1+M2 ;
+    N3=max(N1,N2) ;
+    oj=0 ;
+    oi=M1 ;    
+  case 'd'
+    M3=M1+M2 ;
+    N3=N1+N2 ;
+    oj=N1 ;
+    oi=M1 ;
+  case 'o'
+    M3=max(M1,M2) ;
+    N3=max(N1,N2) ;
+    oj=0;
+    oi=0;
+end
+
+I=zeros(M3,N3,K1) ;
+I(1:M1,1:N1,:) = I1 ;
+I(oi+(1:M2),oj+(1:N2),:) = I2 ;
+         
+axes('Position', [0 0 1 1]) ;
+imagesc(I) ; colormap gray ; hold on ; axis image ; axis off ;
+drawnow ;
+
+K = size(matches, 2) ;
+nans = NaN * ones(1,K) ;
+
+x = [ P1(1,matches(1,:)) ; P2(1,matches(2,:))+oj ; nans ] ;
+y = [ P1(2,matches(1,:)) ; P2(2,matches(2,:))+oi ; nans ] ;
+
+h = line(x(:)', y(:)') ;
+set(h,'Marker','.','Color','g') ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m b/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m
new file mode 100644
index 0000000..2f4af50
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m
@@ -0,0 +1,169 @@
+function h=plotsiftdescriptor(d,f)
+% PLOTSIFTDESCRIPTOR   Plot SIFT descriptor
+%   PLOTSIFTDESCRIPTOR(D) plots the SIFT descriptors D, stored as
+%   columns of the matrix D. D has the same format used by SIFT().
+%
+%   PLOTSIFTDESCRIPTOR(D,F) plots the SIFT descriptors warped to the
+%   SIFT frames F, specified as columns of the matrix F. F has
+%   the same format used by SIFT().
+%
+%   H=PLOTSIFTDESCRIPTOR(...) returns the handle H to the line drawing
+%   representing the descriptors.
+%
+%   REMARK. Currently the function supports only descriptors with 4x4
+%   spatial bins and 8 orientation bins (Lowe's default.)
+%
+%   See also PLOTSIFTFRAME(), PLOTMATCHES(), PLOTSS().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+lowe_compatible = 1 ;
+
+% --------------------------------------------------------------------
+%                                                  Check the arguments
+% --------------------------------------------------------------------
+
+if(size(d,1) ~= 128)
+  error('D should be a 128xK matrix (only standard descriptors accepted)') ;
+end
+
+if nargin > 1
+  if(size(f,1) ~= 4)
+    error('F should be a 4xK matrix');
+  end
+  
+  if(size(f,2) ~= size(f,2))
+    error('D and F must have the same number of columns') ;
+  end
+end
+
+% Descriptors are often non-double numeric arrays
+d = double(d) ;
+K = size(d,2) ;
+if nargin < 2
+  f = repmat([0;0;1;0],1,K) ;
+end
+
+maginf = 3.0 ;
+NBP=4 ;
+NBO=8 ;
+
+% --------------------------------------------------------------------
+%                                                           Do the job
+% --------------------------------------------------------------------
+
+xall=[] ;
+yall=[] ;
+
+for k=1:K
+  SBP = maginf * f(3,k) ;
+  th=f(4,k) ;
+  c=cos(th) ;
+  s=sin(th) ;
+
+  [x,y] = render_descr(d(:,k)) ;
+  xall = [xall SBP*(c*x-s*y)+f(1,k)] ;
+  yall = [yall SBP*(s*x+c*y)+f(2,k)] ;
+end
+
+h=line(xall,yall) ;
+
+% --------------------------------------------------------------------
+%                                                     Helper functions
+% --------------------------------------------------------------------
+
+% Renders a single descriptor
+function [x,y] = render_descr( d )
+
+lowe_compatible=1; 
+NBP=4 ;
+NBO=8 ;
+
+[x,y] = meshgrid(-NBP/2:NBP/2,-NBP/2:NBP/2) ;
+
+% Rescale d so that the biggest peak fits inside the bin diagram
+d = 0.4 * d / max(d(:)) ;
+
+% We have NBP*NBP bins to plot. Here are the centers:
+xc = x(1:end-1,1:end-1) + 0.5 ;
+yc = y(1:end-1,1:end-1) + 0.5 ;
+
+% We swap the order of the bin diagrams because they are stored row
+% major into the descriptor (Lowe's convention that we follow.)
+xc = xc' ;
+yc = yc' ;
+
+% Each bin contains a star with eight tips
+xc = repmat(xc(:)',NBO,1) ;
+yc = repmat(yc(:)',NBO,1) ;
+
+% Do the stars 
+th=linspace(0,2*pi,NBO+1) ;
+th=th(1:end-1) ;
+if lowe_compatible
+  xd = repmat(cos(-th), 1, NBP*NBP ) ;
+  yd = repmat(sin(-th), 1, NBP*NBP ) ;
+else
+  xd = repmat(cos(th), 1, NBP*NBP ) ;
+  yd = repmat(sin(th), 1, NBP*NBP ) ;
+end
+xd = xd .* d(:)' ;
+yd = yd .* d(:)' ;
+
+% Re-arrange in sequential order the lines to draw
+nans = NaN * ones(1,NBP^2*NBO) ;
+x1 = xc(:)' ;
+y1 = yc(:)' ;
+x2 = x1 + xd ;
+y2 = y1 + yd ;
+xstars = [x1;x2;nans] ;
+ystars = [y1;y2;nans] ;
+
+% Horizontal lines of the grid
+nans = NaN * ones(1,NBP+1);
+xh = [x(:,1)' ; x(:,end)' ; nans] ;
+yh = [y(:,1)' ; y(:,end)' ; nans] ;
+
+% Verical lines of the grid
+xv = [x(1,:) ; x(end,:) ; nans] ;
+yv = [y(1,:) ; y(end,:) ; nans] ;
+
+x=[xstars(:)' xh(:)' xv(:)'] ;
+y=[ystars(:)' yh(:)' yv(:)'] ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m b/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m
new file mode 100644
index 0000000..984dfc9
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m
@@ -0,0 +1,119 @@
+function h=plotsiftframe(frames,labels)
+% PLOTSIFTFRAME  Plot SIFT frame
+%   H=PLOTSIFTFRAME(FRAMES) plots the SIFT frames FRAMES and returns
+%   and handle H to the resulting line set. FRAMES has the same format
+%   used by SIFT().
+%
+%   PLOTSIFTFRAME(FRAMES, LABELS) displays nearby the frame centers
+%   the indexes specified by the vector LABELS. This operation is slow
+%   for large sets of frames.
+%
+%   A SIFT frame is denoted by a circle, representing its support, and
+%   one of its radii, representing its orientation. The support is a
+%   disk with radius equal to six times the scale SIGMA of the
+%   frame. If the standard parameters are used for the detector, this
+%   corresponds to four times the standard deviation of the Gaussian
+%   window that has been uses to estimate the orientation, which is in
+%   fact equal to 1.5 times the scale SIGMA.
+%
+%   This function is considerably more efficient if called once on a
+%   whole set of frames as opposed to multiple times, one for each
+%   frame.
+%
+%   See also PLOTMATCHES(), PLOTSIFTDESCRIPTOR(), PLOTSS().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+% --------------------------------------------------------------------
+%                                                  Check the arguments
+% --------------------------------------------------------------------
+
+if size(frames,1) ~= 4
+        error('FRAMES should be a 4xK matrix') ;
+end
+
+K = size(frames,2) ;
+
+if nargin > 1
+        putlabels = 1 ;
+end
+
+% --------------------------------------------------------------------
+%                                                          Do the work
+% --------------------------------------------------------------------
+
+hold on ;
+K=size(frames,2) ;
+thr=linspace(0,2*pi,40) ;
+
+allx = nan*ones(1, 40*K+(K-1)) ;
+ally = nan*ones(1, 40*K+(K-1)) ;
+
+allxf = nan*ones(1, 3*K) ;
+allyf = nan*ones(1, 3*K) ;
+
+putlabel=0 ;
+
+for k=1:K
+        xc=frames(1,k) ;
+        yc=frames(2,k) ;
+        r=1.5*4*frames(3,k) ;
+        th=frames(4,k) ;
+        
+        x = r*cos(thr) + xc ;
+        y = r*sin(thr) + yc ;
+                
+        allx((k-1)*(41) + (1:40)) = x ;
+        ally((k-1)*(41) + (1:40)) = y ;
+
+        allxf((k-1)*3 + (1:2)) = [xc xc+r*cos(th)] ;
+        allyf((k-1)*3 + (1:2)) = [yc yc+r*sin(th)] ;
+        
+        if putlabel
+                x=xc+r ;
+                y=yc ;
+                h=text(x+2,y,sprintf('%d',labels(k))) ;
+                set(h,'Color',[1 0 0]) ;
+                plot(x,y,'r.') ;
+        end
+
+end
+
+h=line([allx nan allxf], [ally nan allyf], 'Color','g','LineWidth',3) ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotss.m b/SD-VBS/benchmarks/sift/src/matlab/plotss.m
new file mode 100644
index 0000000..17868b6
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/plotss.m
@@ -0,0 +1,68 @@
+function plotss(ss,field)
+% PLOTSS  Plot scale space
+%   PLOTSS(SS) plots all octaves of the scale space SS.
+%
+%   See also GAUSSIANSS(), DIFFSS().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+if nargin > 2
+    error('Too many arguments.') ;
+end
+
+omin = ss.omin ;
+smin = ss.smin ;
+nlevels = ss.smax-ss.smin+1 ;
+
+for oi=1:ss.O
+        for si=1:nlevels
+                tightsubplot(nlevels, ss.O, nlevels*(oi-1)+si) ;
+                s = si-1 + smin ;
+                o = oi-1 + omin ;
+                sigma = ss.sigma0 * 2^(s/ss.S + o) ;
+                F=squeeze(ss.octave{oi}(:,:,si)) ;
+                [M,N]=size(F) ;
+                imagesc(squeeze(ss.octave{oi}(:,:,si))) ;       axis image ; axis off ;
+                h=text(M/10,N/20,sprintf('(o,s)=(%d,%d), sigma=%f',o,s,sigma)) ;
+                set(h,'BackgroundColor','w','Color','k') ;
+        end
+end
+
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m b/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m
new file mode 100644
index 0000000..c157507
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m
@@ -0,0 +1,23 @@
+function script_run_profile(dataDir, resultDir, type, common, tooldir)
+
+path(path,common); 
+sift_compile;
+
+I1=readImage([dataDir, '/1.bmp']) ; 
+
+rows = size(I1,1);
+cols = size(I1,2);
+fprintf(1,'Input size\t\t- (%dx%d)\n', rows, cols);
+
+I1=I1-min(I1(:)) ;
+I1=I1/max(I1(:)) ;
+
+%% Timing
+start = photonStartTiming;
+frames1 = sift( I1) ;
+stop = photonEndTiming;
+elapsed = photonReportTiming(start, stop);
+photonPrintTiming(elapsed);
+
+fWriteMatrix(frames1, dataDir);
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift.m b/SD-VBS/benchmarks/sift/src/matlab/sift.m
new file mode 100644
index 0000000..ee79371
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift.m
@@ -0,0 +1,110 @@
+% Input:
+%   I          the input image, with pixel values normalize to lie betwen [0,1].
+%
+% Output:
+%   features    a structure which contains the following fields:
+%       pos     an n*2 matrix containing the (x,y) coordinates of the keypoints
+%               stored in rows.
+%       scale   an n*3 matrix with rows describing the scale of each keypoint
+%               (i.e., first column specifies the octave, second column specifies the interval, and
+%               third column specifies sigma).
+%       orient  a n*1 vector containing the orientations of the keypoints [-pi,pi).
+%       desc    an n*128 matrix with rows containing the feature descriptors 
+%               corresponding to the keypoints.
+
+function [frames]=sift(I)
+
+[M,N,C] = size(I) ;
+
+% Lowe's choices
+S=3 ;
+omin=-1 ;
+O=floor(log2(min(M,N)))-omin-4 ; % Up to 16x16 images
+sigma0=1.6*2^(1/S) ;
+sigman=0.5 ;
+thresh = 0.04 / S / 2 ;
+r = 10 ;
+
+NBP = 4 ;
+NBO = 8 ;
+magnif = 3.0 ;
+
+% Parese input
+compute_descriptor = 0 ;
+discard_boundary_points = 1 ;
+verb = 0 ;
+
+smin = -1;
+smax = S+1;
+intervals = smax - smin + 1;
+
+
+% --------------------------------------------------------------------
+%                                                           Parameters
+% --------------------------------------------------------------------
+
+oframes = [] ;
+frames = [] ;
+descriptors = [] ;
+
+% --------------------------------------------------------------------
+%                                         SIFT Detector and Descriptor
+% --------------------------------------------------------------------
+
+% Compute scale spaces
+%if verb > 0, fprintf('SIFT: computing scale space...') ; end
+
+gss = gaussianss(I,sigman,O,S,omin,-1,S+1,sigma0) ;
+dogss = diffss(gss) ;
+frames = [];
+
+%% To maintain consistency with C code. Once C code is ready, this will be uncommented.
+for o=1:O                               %for o=1:O 
+        
+        % Local maxima of the DOG octave
+        % The 80% tricks discards early very weak points before refinement.
+         
+    idx = siftlocalmax(  dogss.octave{o}, 0.8*thresh  ) ;
+        idx = [idx , siftlocalmax( - dogss.octave{o}, 0.8*thresh)] ;  
+   
+        K=length(idx) ; 
+        [i,j,s] = ind2sub( size( dogss.octave{o} ), idx ) ;
+    
+    y=i-1 ;
+        x=j-1 ;
+        s=s-1+dogss.smin ;
+    oframes = [x(:)';y(:)';s(:)'] ;
+
+%    fWriteMatrix(oframes, '../data/sim');
+        
+        % Remove points too close to the boundary
+        if discard_boundary_points
+                oframes = filter_boundary_points(size(dogss.octave{o}), oframes ) ;
+        end
+                
+        % Refine the location, threshold strength and remove points on edges
+        oframes = siftrefinemx(...
+                oframes, ...
+                dogss.octave{o}, ...
+                dogss.smin, ...
+                thresh, ...
+                r) ;
+    
+   frames = [frames, oframes];     
+ 
+end
+end
+
+%% --------------------------------------------------------------------
+%%                                                              Helpers
+%% --------------------------------------------------------------------
+function oframes=filter_boundary_points(sz, oframes)
+
+sel=find(...
+        oframes(1,:) > 3 & ...
+        oframes(1,:) < sz(2)-3 & ...
+        oframes(2,:) > 3 & ...
+        oframes(2,:) < sz(1)-3 ) ;
+
+oframes=oframes(:,sel) ;
+end
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m b/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m
new file mode 100644
index 0000000..78c0034
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m
@@ -0,0 +1,60 @@
+function sift_compile(type)
+% SIFT_COMPILE  Compile MEX files
+%   Compiling under Windows requires at least Visual C 6 or LCC. You
+%   might try other compilers, but most likely you will need to edit
+%   this file.
+
+siftroot = fileparts(which('siftcompile')) ;
+opts = { '-O', '-I.' } ;
+%opts = { opts{:}, '-v' } ;
+
+if nargin < 1
+    type = 'visualc' ;
+end
+
+switch computer
+  case 'PCWIN'
+    warning('NOTE: compiling has been tested only with Visual C 6-7 and LCC') ;
+    switch type
+      case 'visualc'
+        lib{1}=[matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'] ;
+        lib{2}=[matlabroot '\extern\lib\win32\microsoft\msvc60\libmwlapack.lib'];
+        lib{3}=[matlabroot '\extern\lib\win32\microsoft\msvc71\libmwlapack.lib'];
+      case 'lcc'
+        lib{1}=[matlabroot '\extern\lib\win32\lcc\libmwlapack.lib'] ;
+    end
+    found=0;
+    for k=1:length(lib)
+      fprintf('Trying LAPACK lib ''%s''\n',lib{k}) ;
+      found=exist(lib{k}) ;
+      if found ~= 0
+        break ;
+      end
+    end
+    if found == 0
+      error('Could not find LAPACK library. Please edit this M-file to fix the issue.');
+    end
+    opts = {opts{:}, '-DWINDOWS'} ;
+    opts = {opts{:}, lib{k}} ;
+    
+  case 'MAC'
+    opts = {opts{:}, '-DMACOSX'} ;
+    opts = {opts{:}, 'CFLAGS=\$CFLAGS -faltivec'} ;
+  
+  case 'GLNX86'
+    opts = {opts{:}, '-DLINUX' } ;
+        
+  otherwise
+    error(['Unsupported architecture ', computer, '. Please edit this M-mfile to fix the issue.']) ;    
+end
+
+mex('imsmooth.c',opts{:}) ;
+mex('siftlocalmax.c',opts{:}) ;
+mex('siftrefinemx.c',opts{:}) ;
+mex('siftormx.c',opts{:}) ;
+mex('siftdescriptor.c',opts{:}) ;
+mex('siftmatch.c',opts{:}) ;
+
+    
+
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m b/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m
new file mode 100644
index 0000000..a6f0c66
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m
@@ -0,0 +1,110 @@
+cd % SIFT_DEMO2  Demonstrate SIFT code (2)
+%   This is similar to SIFT_DEMO().
+%
+%   See also SIFT_DEMO().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+I1=imreadbw('data/landscape-a.jpg') ; % I1=I1(1:2:end,:) ;
+I2=imreadbw('data/landscape-b.jpg') ; % I2=I2(1:2:end,:) ;
+I1c=double(imread('data/landscape-a.jpg'))/255.0 ;
+I2c=double(imread('data/landscape-b.jpg'))/255.0 ;
+
+I1=imsmooth(I1,.1) ;
+I2=imsmooth(I2,.1) ;
+
+I1=I1-min(I1(:)) ;
+I1=I1/max(I1(:)) ;
+I2=I2-min(I2(:)) ;
+I2=I2/max(I2(:)) ;
+
+S=3 ;
+
+fprintf('Computing frames and descriptors.\n') ;
+[frames1,descr1,gss1,dogss1] = sift( I1, 'Verbosity', 1, 'Threshold', ...
+                                     0.005, 'NumLevels', S ) ;
+[frames2,descr2,gss2,dogss2] = sift( I2, 'Verbosity', 1, 'Threshold', ...
+                                     0.005, 'NumLevels', S ) ;
+
+figure(11) ; clf ; plotss(dogss1) ; colormap gray ;
+figure(12) ; clf ; plotss(dogss2) ; colormap gray ;
+drawnow ;
+
+figure(2) ; clf ;
+tightsubplot(1,2,1) ; imagesc(I1) ; colormap gray ; axis image ;
+hold on ;
+h=plotsiftframe( frames1 ) ; set(h,'LineWidth',2,'Color','g') ;
+h=plotsiftframe( frames1 ) ; set(h,'LineWidth',1,'Color','k') ;
+
+tightsubplot(1,2,2) ; imagesc(I2) ; colormap gray ; axis image ;
+hold on ;
+h=plotsiftframe( frames2 ) ; set(h,'LineWidth',2,'Color','g') ;
+h=plotsiftframe( frames2 ) ; set(h,'LineWidth',1,'Color','k') ;
+
+fprintf('Computing matches.\n') ;
+% By passing to integers we greatly enhance the matching speed (we use
+% the scale factor 512 as Lowe's, but it could be greater without
+% overflow)
+descr1=uint8(512*descr1) ;
+descr2=uint8(512*descr2) ;
+tic ; 
+matches=siftmatch( descr1, descr2, 3 ) ;
+fprintf('Matched in %.3f s\n', toc) ;
+
+figure(3) ; clf ;
+plotmatches(I1c,I2c,frames1(1:2,:),frames2(1:2,:),matches,...
+  'Stacking','v') ;
+drawnow ;
+
+% Movie
+figure(4) ; set(gcf,'Position',[10 10 1024 512]) ;
+figure(4) ; clf ;
+tightsubplot(1,1);
+imagesc(I1) ; colormap gray ; axis image ; hold on ;
+h=plotsiftframe( frames1 ) ; set(h,'LineWidth',1,'Color','g') ;
+h=plot(frames1(1,:),frames1(2,:),'r.') ;
+MOV(1)=getframe ;
+
+figure(4) ; clf ;
+tightsubplot(1,1);
+imagesc(I2) ; colormap gray ; axis image ; hold on ;
+h=plotsiftframe( frames2 ) ; set(h,'LineWidth',1,'Color','g') ;
+h=plot(frames2(1,:),frames2(2,:),'r.') ;
+MOV(2)=getframe ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css
new file mode 100644
index 0000000..b9bd8e3
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css
@@ -0,0 +1,159 @@
+/** file:        default.css
+ ** author:      Andrea Vedaldi
+ ** description: Default CSS sylesheet for SIFT_GENDOC.PL
+ **/
+
+/* AUTORIGHTS */
+
+h1
+{
+  color: #3a3a3a ;
+}
+
+pre
+{
+  font-family: monaco, courier, monospace ;
+  font-size: 14px ;
+  color: #3a3a3a ;
+}
+
+body
+{
+  margin: 10px ;
+  padding: 1em ;
+  right: 0px ;
+  font-family: arial, sans-serif ;
+
+  color: black ;
+  background-color: #fffffa ; 
+}
+
+/* ------------------------------------------------------------------
+ *                                                             Module
+ * --------------------------------------------------------------- */
+
+div.module
+{
+  margin-bottom: 1em ;
+  //background-color: #ffffff ;
+}
+
+div.module h1
+{
+  margin: 0px ;
+  margin-bottom: 0.5em ;
+  padding: 0px ;
+  font-size: 1.5em ;
+  border-bottom: 2px #3a3a3a solid ;
+}
+
+div.module pre
+{
+  padding-left: 1em ;
+}
+
+div.module div.footer
+{
+  clear: both ;
+}
+
+div.module div.footer :link
+{
+  color: white ;
+  text-decoration: none ;
+}
+
+/* ------------------------------------------------------------------
+ *                                                       Module index
+ * --------------------------------------------------------------- */
+
+div.module div.index
+{
+  font-family: sans-serif ;
+  font-size: 0.8em ;
+
+  width: 15em ;
+  float: right ;
+  border: 1px #ccc0c0 solid ;
+  background-color: #fcf0f0 ;
+}
+
+div.module div.index h1
+{
+  font-size: 1em ;  
+  font-weight: bold ;
+  text-align: center ;
+  border: none ;
+  padding: 0 ;
+  margin: 0 ;
+}
+
+div.module div.index ul
+{
+  list-style-type: square ;
+  list-style-position: inside ;
+  color: #3a3a3a ;
+  padding: 0 ;
+  margin: 0 ;
+  padding: 0.3em ;
+}
+
+div.module div.index ul li
+{
+  margin-bottom: 0.1em ;
+}
+
+
+/* ------------------------------------------------------------------
+ *                                                              Mfile
+ * --------------------------------------------------------------- */
+
+div.mfile
+{
+  background-color: white ;
+  margin-bottom: 1em ;
+  border: 1px #aaaaaa solid ;
+}
+
+div.mfile h1
+{
+  margin: 0 ;
+  padding: 0 ;
+  padding-left: 10px ;
+  font-size: 1.3em ;  
+  background-color: #f0f0ff ;
+}
+
+div.mfile h1 span.name
+{
+  font-family: monaco, courier, fixed ;
+}
+
+div.mfile h1 span.brief
+{
+  padding-left: 10px ; 
+  font-size: 0.9em ;
+  font-weight: normal ;  
+}
+
+div.mfile pre
+{
+  padding-left: 1em ;
+}
+
+div.mfile div.footer
+{
+  font-size: 0.8em ;
+  padding: 0.3em ;
+}
+
+div.mfile div.footer a
+{
+  text-decoration: none ;
+  color: #777777 ;
+}
+
+div.mfile div.footer a:hover
+{
+  text-decoration: underline ;
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m
new file mode 100644
index 0000000..417f9c5
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m
@@ -0,0 +1,50 @@
+function sift_gendoc
+% SIFT_GENDOC  Generate documentation
+%   This function extracts the documentation from the MEX files and
+%   creates an HTML manual.
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+rootdir   = fileparts(which('siftgendoc')) ;
+d=pwd;
+cd([rootdir]) ;
+res = perl([rootdir '/siftgendoc.pl'],[rootdir '/dom/index.html']) ;
+cd(d) ;
+fprintf(res) ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl
new file mode 100755
index 0000000..6bc0b5e
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl
@@ -0,0 +1,296 @@
+#!/usr/bin/perl -w 
+## file:        sift_gendoc.pl
+## author:      Andrea Vedaldi
+## description: Summarize MATLAB M-Files docs
+
+# AUTORIGHTS
+# Copyright (c) 2006 The Regents of the University of California.
+# All Rights Reserved.
+# 
+# Created by Andrea Vedaldi
+# UCLA Vision Lab - Department of Computer Science
+# 
+# Permission to use, copy, modify, and distribute this software and its
+# documentation for educational, research and non-profit purposes,
+# without fee, and without a written agreement is hereby granted,
+# provided that the above copyright notice, this paragraph and the
+# following three paragraphs appear in all copies.
+# 
+# This software program and documentation are copyrighted by The Regents
+# of the University of California. The software program and
+# documentation are supplied "as is", without any accompanying services
+# from The Regents. The Regents does not warrant that the operation of
+# the program will be uninterrupted or error-free. The end-user
+# understands that the program was developed for research purposes and
+# is advised not to rely exclusively on the program for any reason.
+# 
+# This software embodies a method for which the following patent has
+# been issued: "Method and apparatus for identifying scale invariant
+# features in an image and use of same for locating an object in an
+# image," David G. Lowe, US Patent 6,711,293 (March 23,
+# 2004). Provisional application filed March 8, 1999. Asignee: The
+# University of British Columbia.
+# 
+# IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+# FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+# INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+# ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+# CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+# A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+# BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+# MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+# Debugging level
+$verb = 0 ;
+
+# PDF document location
+$pdfdoc = 'sift.pdf' ;
+
+# This is the queue of directories to process.
+# It gets filled recursively.
+@dir_fifo = ('.') ;
+
+# This will hold an entry for each m-file found
+%mfiles = () ;
+
+# This will hold an entry for each module found 
+%modules = () ;
+
+# #ARGV is the index of the last element, which is 0
+if ($#ARGV == 0) {
+    open(FOUT,">$ARGV[0]") ;
+    print STDERR "Writing to file '$ARGV[0]'.\n" if $verb ;
+} else {
+    *FOUT= *STDOUT; 
+    print STDERR "Using standard output.\n" if $verb ; 
+}
+
+# Each module is a subdirectory. The subdirectory is used
+# as the module ID.
+
+while ($module_path = shift(@dir_fifo)) { 
+    print STDERR "=> '$module_path'\n" if $verb ;
+    
+    # get a first version of the module name
+    $module_path =~ m/.*\/([^\\]+)$/ ;
+    $module_name = $1 ;
+
+    # start a new module    
+    $module = {
+        'id'          => $module_path,
+        'path'        => $module_path,
+        'name'        => $module_name,
+        'mfiles'      => [],
+        'description' => ""
+    } ;
+    
+    # .................................................................
+    opendir(DIRHANDLE, $module->{'path'}) ;
+    FILE: foreach (sort readdir(DIRHANDLE)) {
+        $name = $_ ;
+        $path = $module->{'path'} . "/" . $_ ;
+
+        # push if a directory
+        $_ = $path ;
+        if (-d) {
+            next FILE if /(\.$|\.\.$)/ ;
+            push( @dir_fifo, "$_" ) ;
+            next FILE ;
+        }
+        
+        # parse if .m and not test_
+        next FILE unless (/.+\.m$/) ;
+        next FILE if (/test_.x*/) ;
+        $name =~ m/(.*).m/ ;
+        $name = $1 ;
+        print STDERR "  . m-file '$name'\n" if $verb ;
+        my ($brief,$description) = get_comment($path) ;
+
+        # topic description?
+        if (/overview/) {
+            print STDERR "    * module description\n" if $verb ;
+            $module->{'id'}          = $name ;
+            $module->{'name'}        = $brief ;
+            $module->{'description'} = $description ;
+            next FILE ;
+        }
+
+        # use names as IDs
+        $id = $name ;
+        $id =~ tr/A-Z/a-z/ ;
+
+        # create a new mfile object
+        $mfile = {
+            'id'          => $id,
+            'path'        => $path,
+            'name'        => $name,
+            'brief'       => $brief,
+            'module'      => $module,
+            'description' => $description
+        } ;
+        
+        # add a reference to this mfile
+        # object to the global mfile list
+        $mfiles{$id} = $mfile ;
+
+        # add a reference to this mfile
+        # object to the current module
+        push( @{$module->{'mfiles'}}, $mfile) ;
+    }
+    closedir(DIRHANDLE) ;
+    # ................................................................
+
+    # add a reference to the current module to the global
+    # module list
+    $modules{$module->{'id'}} = $module ;
+}
+
+# ....................................................................
+#                                                  write documentation
+# ....................................................................
+
+print FOUT <<EOF;
+<html>
+  <head>
+    <link href="default.css" rel="stylesheet" type="text/css"/>
+  </head>
+  <body>
+EOF
+
+# sort modules by path
+sub criterion { $modules{$a}->{'path'} cmp $modules{$b}->{'path'} ; }
+
+MODULE:
+foreach $id ( sort criterion keys %modules ) {
+    my $module = $modules{$id} ;   
+    my $rich_description = make_rich($module->{'description'}) ;
+
+    next MODULE if $#{$module->{'mfiles'}} < 0 and length($rich_description) == 0;
+
+    print FOUT <<EOF;
+    <div class='module' id='$module->{"id"}'>
+      <h1>$module->{'name'}</h1>
+      <div class='index'>
+        <h1>Module contents</h1>
+        <ul>
+EOF
+    foreach( @{$module->{'mfiles'}} ) {
+        print FOUT "        <li><a href='#$_->{'id'}'>" 
+                 . "$_->{'name'}</a></li>\n" ;
+    }
+    print FOUT <<EOF;
+        </ul>
+      </div>
+      <pre>
+      $rich_description
+      </pre>
+      <div class="footer">
+      </div>
+    </div>
+EOF
+}
+
+foreach $id (sort keys %mfiles) {
+    my $mfile = $mfiles{$id} ;
+    my $rich_description = make_rich($mfile->{'description'}) ;
+
+    print FOUT <<EOF;
+    <div class="mfile" id='$mfile->{"id"}'>
+      <h1>
+       <span class="name">$mfile->{"name"}</span>
+       <span class="brief">$mfile->{"brief"}</span>
+      </h1>
+      <pre>
+$rich_description
+      </pre>
+      <div class="footer">
+        <a href="#$mfile->{'module'}->{'id'}">
+        $mfile->{'module'}->{'name'}
+        </a>
+      </div>
+    </div>
+EOF
+}
+
+print FOUT "</body></html>" ;
+
+# Close file
+close FOUT ;
+
+
+# -------------------------------------------------------------------------
+sub get_comment {
+# -------------------------------------------------------------------------
+    local $_ ;
+    my $full_name = $_[0] ;
+    my @comment = () ;
+    
+    open IN,$full_name ;
+  SCAN: 
+    while( <IN> ) {
+        next if /^function/ ;
+        last SCAN unless ( /^%/ ) ;
+        push(@comment, substr("$_",1)) ;
+    }
+    close IN ;
+
+    my $brief = "" ;
+    if( $#comment >= 0 && $comment[0] =~ m/^\s*\w+\s+(.*)$/ ) {
+        $brief = $1 ;
+        splice (@comment, 0, 1) ;
+    }
+
+    # from the first line
+    return ($brief, join("",@comment)) ;
+}
+
+# -------------------------------------------------------------------------
+sub make_rich {
+# -------------------------------------------------------------------------
+    local $_ = $_[0] ;
+    s/([A-Z]+[A-Z0-9_]*)\(([^\)]*)\)/${\make_link($1,$2)}/g ;
+    s/PDF:([A-Z0-9_\-:.]+[A-Z0-9])/${\make_pdf_link($1)}/g ;
+    return $_ ;
+}
+
+# -------------------------------------------------------------------------
+sub make_link {
+# -------------------------------------------------------------------------
+    local $_ ;
+    my $name = $_[0] ;
+    my $arg  = $_[1] ;
+    my $id   = $name ;
+
+    # convert name to lower case and put into $_
+    $id =~ tr/A-Z/a-z/ ;
+    
+    # get mfile
+    my $mfile = $mfiles{$id} ;
+    my $module = $modules{$id} ;
+    
+    # return as appropriate    
+    if($mfile) {
+        return "<a href='#$id'>" . $name . "</a>" . "(" . $arg . ")" ;
+    } elsif($module) {
+      return "<a class='module' href='#$id'>" . $name . 
+          "</a>" . "(" . $arg . ")" ;
+    } else {
+        return $name . "(" . $arg .")" ;
+    }
+}
+
+
+# -------------------------------------------------------------------------
+sub make_pdf_link {
+# -------------------------------------------------------------------------
+    local $_ ;
+    my $name = $_[0] ;
+    my $id   = $name ;
+
+    # convert name to lower case and put into $_
+    $id =~ tr/A-Z/a-z/ ;
+    
+    return "<a href='${pdfdoc}#$id'>PDF:$1</a>" ;
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m b/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m
new file mode 100644
index 0000000..71557e4
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m
@@ -0,0 +1,33 @@
+% SIFT_OVERVIEW  Scale-Invariant Feature Transfrom
+%
+%  This is a MATLAB/C implementation of SIFT detector and descriptor
+%  [1]. You can:
+%
+%  * Use SIFT() to detect the SIFT frames (keypoints) of a given image
+%    and compute their descriptors. Then you can use SIFTMATCH() to
+%    match the descriptors.
+%
+%  * Use PLOTSS(), PLOTSIFTDESCRIPTOR(), PLOTSIFTFRAME(),
+%    PLOTMATCHES() to visualize the results.
+%
+%  As SIFT is implemented by several reusable M and MEX files, you can
+%  also run portions of the algorithm, or change them. Specifically,
+%  you can:
+%
+%  * Use SIFTDESCRIPTOR() to compute the SIFT descriptor from a list
+%    of frames and a scale space or plain image.
+%
+%  * Use GAUSSIANSS() and DIFFSS() to compute the Gaussian and DOG
+%    scale spaces.
+%
+%  * Use SIFTLOCALMAX(), SIFTREFINEMX(), SIFTORMX() to manually
+%    extract the SIFT frames from the DOG scale space. More in
+%    general, you can use SIFTLOCALMAX() to find maximizers of any
+%    multi-dimensional arrays.
+%
+%  REFERENCES
+%  [1] D. G. Lowe, "Distinctive image features from scale-invariant
+%      keypoints," IJCV, vol. 2, no. 60, pp. 91 110, 2004.
+%
+%  See also PDF:SIFT.INTRODUCTION.
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c
new file mode 100644
index 0000000..63f4830
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c
@@ -0,0 +1,524 @@
+/* file:        siftdescriptor
+** author:      Andrea Vedaldi
+** description: Compute SIFT descriptors
+**/
+
+/* AUTORIGHTS
+Copyright (c) 2006 The Regents of the University of California.
+All Rights Reserved.
+
+Created by Andrea Vedaldi
+UCLA Vision Lab - Department of Computer Science
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for educational, research and non-profit purposes,
+without fee, and without a written agreement is hereby granted,
+provided that the above copyright notice, this paragraph and the
+following three paragraphs appear in all copies.
+
+This software program and documentation are copyrighted by The Regents
+of the University of California. The software program and
+documentation are supplied "as is", without any accompanying services
+from The Regents. The Regents does not warrant that the operation of
+the program will be uninterrupted or error-free. The end-user
+understands that the program was developed for research purposes and
+is advised not to rely exclusively on the program for any reason.
+
+This software embodies a method for which the following patent has
+been issued: "Method and apparatus for identifying scale invariant
+features in an image and use of same for locating an object in an
+image," David G. Lowe, US Patent 6,711,293 (March 23,
+2004). Provisional application filed March 8, 1999. Asignee: The
+University of British Columbia.
+
+IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+*/
+
+/*
+  REMARKS. The use of strcasecmp makes the function POSIX but not ANSI
+  compliant. When compling with Altivec, GCC Altivec extensions are
+  supported.
+*/
+
+#define LOWE_COMPATIBLE
+
+#include"mexutils.c"
+#include<stdlib.h>
+#include<math.h>
+
+#ifdef WINDOWS
+#include<string.h>
+#ifndef __cplusplus
+#define sqrtf(x)    ((float)sqrt((double)(x)))
+#define powf(x,y)   ((float)pow((double)(x),(double)(y)))
+#define fabsf(x)    ((float)fabs((double)(x)))
+#define sinf(x)     ((float)sin((double)(x)))
+#define cosf(x)     ((float)cos((double)(x)))
+#define expf(x)     ((float)exp((double)(x)))
+#define atan2f(x,y) ((float)atan2((double)(x),(double)(y)))
+#endif
+#else
+#include<strings.h>
+#endif
+
+/* Altivec and Accelerate support.
+ * Very crude at this time.
+ */
+#if defined( MACOSX ) && defined( __ALTIVEC__ )
+#include<Accelerate/Accelerate.h>
+typedef union 
+{
+  float x[4] ;
+  vFloat vec ;
+} float4 ;
+#endif
+
+#define greater(a,b) a > b
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+
+enum {SCALESPACE, NOSCALESPACE} ;
+
+enum {PROP_MAGNIF=0,
+      PROP_NBP,
+      PROP_NBO,
+      PROP_UNKNOWN} ;
+
+char const * properties [4] = 
+  { "Magnif",
+    "NumSpatialBins",
+    "NumOrientBins",
+     0L
+  } ;
+
+/** Fast fmodf for 2*PI
+ **/
+/*inline*/
+float fast_mod(float th)
+{
+  while(th < 0) th += 2*M_PI ;
+  while(th > 2*M_PI) th -= 2*M_PI ;
+  return th ;
+}
+
+/** Fast floor. Equivalent to (int) floor(x)
+ **/
+/*inline*/
+int fast_floor(float x)
+{
+  return (int)( x - ((x>=0)?0:1) ) ; 
+}
+
+/** Normalizes in norm L_2 a descriptor.
+ **/
+void
+normalize_histogram(float* L_begin, float* L_end)
+{
+  float* L_iter ;
+  float norm=0.0 ;
+
+  for(L_iter = L_begin; L_iter != L_end ; ++L_iter)
+    norm += (*L_iter) * (*L_iter) ;
+
+  norm = sqrtf(norm) ;
+  /*  mexPrintf("%f\n",norm) ;*/
+
+  for(L_iter = L_begin; L_iter != L_end ; ++L_iter)
+    *L_iter /= norm ;
+}
+
+/** @brief MATLAB Driver.
+ **/
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int M,N,S=0,smin=0,K,num_levels=0 ;
+  const int* dimensions ;
+  const double* P_pt ;
+  const double* G_pt ;
+  float* descr_pt ;
+  float* buffer_pt ;
+  float sigma0 ;
+  float magnif = 3.0f ; /* Spatial bin extension factor. */
+  int NBP = 4 ;         /* Number of bins for one spatial direction (even). */
+  int NBO = 8 ;         /* Number of bins for the ortientation. */
+  int mode = NOSCALESPACE ;
+  int buffer_size=0;
+
+  enum {IN_G=0,IN_P,IN_SIGMA0,IN_S,IN_SMIN} ;
+  enum {OUT_L=0} ;
+
+  /* ------------------------------------------------------------------
+  **                                                Check the arguments
+  ** --------------------------------------------------------------- */ 
+ 
+  if (nin < 3) {
+    mexErrMsgTxt("At least three arguments are required") ;
+  } else if (nout > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+                
+  if( !uIsRealScalar(in[IN_SIGMA0]) ) {
+    mexErrMsgTxt("SIGMA0 should be a real scalar") ;
+  }
+        
+  if(!mxIsDouble(in[IN_G]) ||
+     mxGetNumberOfDimensions(in[IN_G]) > 3) {
+    mexErrMsgTxt("G should be a real matrix or 3-D array") ;
+  }
+  
+  sigma0 = (float) *mxGetPr(in[IN_SIGMA0]) ;
+  
+  dimensions = mxGetDimensions(in[IN_G]) ;
+  M = dimensions[0] ;
+  N = dimensions[1] ;
+  G_pt = mxGetPr(in[IN_G]) ;
+  
+  P_pt = mxGetPr(in[IN_P]) ;    
+  K = mxGetN(in[IN_P]) ;
+  
+  if( !uIsRealMatrix(in[IN_P],-1,-1)) {
+    mexErrMsgTxt("P should be a real matrix") ;
+  }
+
+  if ( mxGetM(in[IN_P])  == 4) {
+    /* Standard (scale space) mode */ 
+    mode = SCALESPACE ;
+    num_levels = dimensions[2] ;
+    
+    if(nin < 5) {
+      mexErrMsgTxt("Five arguments are required in standard mode") ;
+    }
+    
+    if( !uIsRealScalar(in[IN_S]) ) {
+      mexErrMsgTxt("S should be a real scalar") ;
+    }
+    
+    if( !uIsRealScalar(in[IN_SMIN]) ) {
+      mexErrMsgTxt("SMIN should be a real scalar") ;
+    }
+    
+    if( !uIsRealMatrix(in[IN_P],4,-1)) {
+      mexErrMsgTxt("When the e mode P should be a 4xK matrix.") ;
+    }
+    
+    S = (int)(*mxGetPr(in[IN_S])) ;
+    smin = (int)(*mxGetPr(in[IN_SMIN])) ;
+    
+  } else if (  mxGetM(in[IN_P])  == 3 ) {
+    mode = NOSCALESPACE ;
+    num_levels = 1 ;
+    S      = 1 ;
+    smin   = 0 ;
+  } else {
+    mexErrMsgTxt("P should be either a 3xK or a 4xK matrix.") ;
+  }
+
+  /* Parse the property-value pairs */
+  {
+    char str [80] ;
+    int arg = (mode == SCALESPACE) ? IN_SMIN + 1 : IN_SIGMA0 + 1 ;
+
+    while(arg < nin) {
+      int k ;
+
+      if( !uIsString(in[arg],-1) ) {
+        mexErrMsgTxt("The first argument in a property-value pair"
+                     " should be a string") ;
+      }
+      mxGetString(in[arg], str, 80) ;
+
+#ifdef WINDOWS      
+      for(k = 0 ; properties[k] && strcmpi(str, properties[k]) ; ++k) ;
+#else
+      for(k = 0 ; properties[k] && strcasecmp(str, properties[k]) ; ++k) ;
+#endif
+
+      switch (k) {
+      case PROP_NBP:
+        if( !uIsRealScalar(in[arg+1]) ) {
+          mexErrMsgTxt("'NumSpatialBins' should be real scalar") ;
+        }
+        NBP = (int) *mxGetPr(in[arg+1]) ;
+        if( NBP <= 0 || (NBP & 0x1) ) {
+          mexErrMsgTxt("'NumSpatialBins' must be positive and even") ;
+        }
+        break ;
+
+      case PROP_NBO:
+        if( !uIsRealScalar(in[arg+1]) ) {
+          mexErrMsgTxt("'NumOrientBins' should be a real scalar") ;
+        }
+        NBO = (int) *mxGetPr(in[arg+1]) ;
+        if( NBO <= 0 ) {
+          mexErrMsgTxt("'NumOrientlBins' must be positive") ;
+        }
+        break ;
+
+      case PROP_MAGNIF:
+        if( !uIsRealScalar(in[arg+1]) ) {
+          mexErrMsgTxt("'Magnif' should be a real scalar") ;
+        }
+        magnif = (float) *mxGetPr(in[arg+1]) ;
+        if( magnif <= 0 ) {
+          mexErrMsgTxt("'Magnif' must be positive") ;
+        }
+        break ;
+
+      case PROP_UNKNOWN:
+        mexErrMsgTxt("Property unknown.") ;
+        break ;
+      }
+      arg += 2 ;
+    }
+  }
+  
+  /* -----------------------------------------------------------------
+   *                                   Pre-compute gradient and angles
+   * -------------------------------------------------------------- */
+  /* Alloc two buffers and make sure their size is multiple of 128 for
+   * better alignment (used also by the Altivec code below.)
+   */
+  buffer_size = (M*N*num_levels + 0x7f) & (~ 0x7f) ;
+  buffer_pt = (float*) mxMalloc( sizeof(float) * 2 * buffer_size ) ;
+  descr_pt  = (float*) mxCalloc( NBP*NBP*NBO*K,  sizeof(float)  ) ;
+
+  {
+    /* Offsets to move in the scale space. */
+    const int yo = 1 ;
+    const int xo = M ;
+    const int so = M*N ;
+    int x,y,s ;
+
+#define at(x,y) (*(pt + (x)*xo + (y)*yo))
+
+    /* Compute the gradient */
+    for(s = 0 ; s < num_levels ; ++s) {
+      const double* pt = G_pt + s*so ;
+      for(x = 1 ; x < N-1 ; ++x) {
+        for(y = 1 ; y < M-1 ; ++y) {
+          float Dx = 0.5 * ( at(x+1,y) - at(x-1,y) ) ;
+          float Dy = 0.5 * ( at(x,y+1) - at(x,y-1) ) ;
+          buffer_pt[(x*xo+y*yo+s*so) + 0          ] = Dx ;
+          buffer_pt[(x*xo+y*yo+s*so) + buffer_size] = Dy ;
+        }
+      }
+    }
+    
+    /* Compute angles and modules */
+    {
+      float* pt = buffer_pt ;
+      int j = 0 ;
+      while (j < N*M*num_levels) {
+
+#if defined( MACOSX ) && defined( __ALTIVEC__ )
+        if( ((unsigned int)pt & 0x7) == 0 && j+3 < N*M*num_levels ) {
+          /* If aligned to 128 bit and there are at least 4 pixels left */
+          float4 a, b, c, d ;
+          a.vec = vec_ld(0,(vector float*)(pt              )) ;
+          b.vec = vec_ld(0,(vector float*)(pt + buffer_size)) ;
+          c.vec = vatan2f(b.vec,a.vec) ;
+          a.x[0] = a.x[0]*a.x[0]+b.x[0]*b.x[0] ;
+          a.x[1] = a.x[1]*a.x[1]+b.x[1]*b.x[1] ;
+          a.x[2] = a.x[2]*a.x[2]+b.x[2]*b.x[2] ;
+          a.x[3] = a.x[3]*a.x[3]+b.x[3]*b.x[3] ;
+          d.vec = vsqrtf(a.vec) ;
+          vec_st(c.vec,0,(vector float*)(pt + buffer_size)) ;
+          vec_st(d.vec,0,(vector float*)(pt              )) ;
+          j += 4 ;
+          pt += 4 ;
+        } else {
+#endif
+          float Dx = *(pt              ) ;
+          float Dy = *(pt + buffer_size) ;
+          *(pt              ) = sqrtf(Dx*Dx + Dy*Dy) ;
+          *(pt + buffer_size) = atan2f(Dy, Dx) ;
+          j += 1 ;
+          pt += 1 ;
+#if defined( MACOSX ) && defined( __ALTIVEC__ )
+        }
+#endif
+
+      }
+    }
+  }
+
+  /* -----------------------------------------------------------------
+   *                                                        Do the job
+   * -------------------------------------------------------------- */ 
+  if(K > 0) {    
+    int p ;
+
+    /* Offsets to move in the buffer */
+    const int yo = 1 ;
+    const int xo = M ;
+    const int so = M*N ;
+
+    /* Offsets to move in the descriptor. */
+    /* Use Lowe's convention. */
+    const int binto = 1 ;
+    const int binyo = NBO * NBP ;
+    const int binxo = NBO ;
+    const int bino  = NBO * NBP * NBP ;
+
+    for(p = 0 ; p < K ; ++p, descr_pt += bino) {
+      /* The SIFT descriptor is a  three dimensional histogram of the position
+       * and orientation of the gradient.  There are NBP bins for each spatial
+       * dimesions and NBO  bins for the orientation dimesion,  for a total of
+       * NBP x NBP x NBO bins.
+       *
+       * The support  of each  spatial bin  has an extension  of SBP  = 3sigma
+       * pixels, where sigma is the scale  of the keypoint.  Thus all the bins
+       * together have a  support SBP x NBP pixels wide  . Since weighting and
+       * interpolation of  pixel is used, another  half bin is  needed at both
+       * ends of  the extension. Therefore, we  need a square window  of SBP x
+       * (NBP + 1) pixels. Finally, since the patch can be arbitrarly rotated,
+       * we need to consider  a window 2W += sqrt(2) x SBP  x (NBP + 1) pixels
+       * wide.
+       */      
+      const float x = (float) *P_pt++ ;
+      const float y = (float) *P_pt++ ;
+      const float s = (float) (mode == SCALESPACE) ? (*P_pt++) : 0.0 ;
+      const float theta0 = (float) *P_pt++ ;
+
+      const float st0 = sinf(theta0) ;
+      const float ct0 = cosf(theta0) ;
+      const int xi = (int) floor(x+0.5) ; /* Round-off */
+      const int yi = (int) floor(y+0.5) ;
+      const int si = (int) floor(s+0.5) - smin ;
+      const float sigma = sigma0 * powf(2, s / S) ;
+      const float SBP = magnif * sigma ;
+      const int W = (int) floor( sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5) ;      
+      int bin ;
+      int dxi ;
+      int dyi ;
+      const float* pt ;
+      float* dpt ;
+
+      /* Check that keypoints are within bounds . */
+
+      if(xi < 0   || 
+         xi > N-1 || 
+         yi < 0   || 
+         yi > M-1 ||
+         ((mode==SCALESPACE) && 
+          (si < 0   ||
+           si > dimensions[2]-1) ) )
+        continue ;
+
+      /* Center the scale space and the descriptor on the current keypoint. 
+       * Note that dpt is pointing to the bin of center (SBP/2,SBP/2,0).
+       */
+      pt  = buffer_pt + xi*xo + yi*yo + si*so ;
+      dpt = descr_pt + (NBP/2) * binyo + (NBP/2) * binxo ;
+     
+#define atd(dbinx,dbiny,dbint) (*(dpt + (dbint)*binto + (dbiny)*binyo + (dbinx)*binxo))
+      
+      /*
+       * Process each pixel in the window and in the (1,1)-(M-1,N-1)
+       * rectangle.
+       */
+      for(dxi = max(-W, 1-xi) ; dxi <= min(+W, N-2-xi) ; ++dxi) {
+        for(dyi = max(-W, 1-yi) ; dyi <= min(+W, M-2-yi) ; ++dyi) {
+
+          /* Compute the gradient. */
+          float mod   = *(pt + dxi*xo + dyi*yo + 0           ) ;
+          float angle = *(pt + dxi*xo + dyi*yo + buffer_size ) ;
+#ifdef LOWE_COMPATIBLE
+          float theta = fast_mod(-angle + theta0) ;
+#else
+          float theta = fast_mod(angle - theta0) ;
+#endif
+          /* Get the fractional displacement. */
+          float dx = ((float)(xi+dxi)) - x;
+          float dy = ((float)(yi+dyi)) - y;
+
+          /* Get the displacement normalized w.r.t. the keypoint orientation
+           * and extension. */
+          float nx = ( ct0 * dx + st0 * dy) / SBP ;
+          float ny = (-st0 * dx + ct0 * dy) / SBP ; 
+          float nt = NBO * theta / (2*M_PI) ;
+
+          /* Get the gaussian weight of the sample. The gaussian window
+           * has a standard deviation of NBP/2. Note that dx and dy are in
+           * the normalized frame, so that -NBP/2 <= dx <= NBP/2. */
+          const float wsigma = NBP/2 ;
+          float win =  expf(-(nx*nx + ny*ny)/(2.0 * wsigma * wsigma)) ;
+
+          /* The sample will be distributed in 8 adijacient bins. 
+           * Now we get the ``lower-left'' bin. */
+          int binx = fast_floor( nx - 0.5 ) ;
+          int biny = fast_floor( ny - 0.5 ) ;
+          int bint = fast_floor( nt ) ;
+          float rbinx = nx - (binx+0.5) ;
+          float rbiny = ny - (biny+0.5) ;
+          float rbint = nt - bint ;
+          int dbinx ;
+          int dbiny ;
+          int dbint ;
+
+          /* Distribute the current sample into the 8 adijacient bins. */
+          for(dbinx = 0 ; dbinx < 2 ; ++dbinx) {
+            for(dbiny = 0 ; dbiny < 2 ; ++dbiny) {
+              for(dbint = 0 ; dbint < 2 ; ++dbint) {
+                
+                if( binx+dbinx >= -(NBP/2) &&
+                    binx+dbinx <   (NBP/2) &&
+                    biny+dbiny >= -(NBP/2) &&
+                    biny+dbiny <   (NBP/2) ) {
+                  float weight = win 
+                    * mod 
+                    * fabsf(1 - dbinx - rbinx)
+                    * fabsf(1 - dbiny - rbiny)
+                    * fabsf(1 - dbint - rbint) ;
+
+                  atd(binx+dbinx, biny+dbiny, (bint+dbint) % NBO) += weight ;
+                }
+              }            
+            }
+          }
+        }  
+      }
+
+      {
+        /* Normalize the histogram to L2 unit length. */        
+        normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ;
+        
+        /* Truncate at 0.2. */
+        for(bin = 0; bin < NBO*NBP*NBP ; ++bin) {
+          if (descr_pt[bin] > 0.2) descr_pt[bin] = 0.2;
+        }
+        
+        /* Normalize again. */
+        normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ;
+      }
+    }
+  }
+
+  /* Restore pointer to the beginning of the descriptors. */
+  descr_pt -= NBO*NBP*NBP*K ;
+
+  {
+    int k ;
+    double* L_pt ;
+    out[OUT_L] = mxCreateDoubleMatrix(NBP*NBP*NBO, K, mxREAL) ;
+    L_pt = mxGetPr(out[OUT_L]) ;
+    for(k = 0 ; k < NBP*NBP*NBO*K ; ++k) {
+      L_pt[k] = descr_pt[k] ;
+    }
+  }
+
+  mxFree(descr_pt) ;  
+  mxFree(buffer_pt) ;
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m
new file mode 100644
index 0000000..7764897
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m
@@ -0,0 +1,79 @@
+% SIFTDESCRIPTOR  Compute SIFT descriptors
+%   DESCR = SIFTDESCRIPTOR(G, P, SIGMA0, S, MINS) computes the SIFT
+%   descriptors DESCR of the SIFT frames P defined on the octave G of
+%   a Gaussian scale space. SIGMA0, S and MINS are the the parameter
+%   of the scale space as explained in PDF:SIFT.USER.SS. P has one
+%   column per frame, specifiying the center X1,X2, the scal index s
+%   and the orientation THETA of the frame.
+%
+%   As the function operates on a single octave, in order to process
+%   features spanning several octaves, one must group them and call
+%   SIFTDESCRIPTOR() once per each octave.
+%    
+%   DESCR = SIFTDESCRIPTOR(I, P, SIGMA0) operates on a plain image
+%   I. The image I is assumed to be pre-smoothed at scale SIGMA0 and P
+%   is a matrix with a column per frame, specifying the center X1,X2
+%   and the orientation THETA (but NOT the scale).
+%
+%   REMARK. The scale parameter s in P is the scale index, NOT the
+%   scale coordinate (see the PDF doc. for a discussion).
+%
+%   Other parameters can be specfied as option-value paris. These
+%   are
+%
+%   'Magnif' [3.0]
+%      Frame magnification factor. Each spatial bin of the SIFT
+%      histogram has an exentsion equal to magnif * sigma, where
+%      magnif is the frame magnification factor and sigma is the scale
+%      of the frame.
+%
+%   'NumSpatialBins' [4]
+%      This parameter specifies the number of spatial bins in each
+%      spatial direction x1 and x2. It must be a positive and even
+%      number.
+%
+%   'NumOrientBins' [8]
+%      This parameter specifies the number of orietnation bins. It
+%      must be a positive number.
+%
+%   See also SIFT(), GAUSSIANSS(), DIFFSS(), SIFTLOCALMAX(), 
+%            PDF:SIFT.USER.DESCRIPTOR.
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64
new file mode 100755
index 0000000..f094f56
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx
new file mode 100755
index 0000000..4eb94d8
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c
new file mode 100644
index 0000000..3646692
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c
@@ -0,0 +1,291 @@
+/* file:        siftlocalmax.c
+** author:      Andrea Vedaldi
+** description: Find local maximizer of multi-dimensional array.
+**/
+
+/* AUTORIGHTS
+Copyright (c) 2006 The Regents of the University of California.
+All Rights Reserved.
+
+Created by Andrea Vedaldi
+UCLA Vision Lab - Department of Computer Science
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for educational, research and non-profit purposes,
+without fee, and without a written agreement is hereby granted,
+provided that the above copyright notice, this paragraph and the
+following three paragraphs appear in all copies.
+
+This software program and documentation are copyrighted by The Regents
+of the University of California. The software program and
+documentation are supplied "as is", without any accompanying services
+from The Regents. The Regents does not warrant that the operation of
+the program will be uninterrupted or error-free. The end-user
+understands that the program was developed for research purposes and
+is advised not to rely exclusively on the program for any reason.
+
+This software embodies a method for which the following patent has
+been issued: "Method and apparatus for identifying scale invariant
+features in an image and use of same for locating an object in an
+image," David G. Lowe, US Patent 6,711,293 (March 23,
+2004). Provisional application filed March 8, 1999. Asignee: The
+University of British Columbia.
+
+IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+*/
+
+#include"mex.h"
+#include<mexutils.c>
+#include<stdlib.h>
+
+/** Matlab driver.
+ **/
+#define greater(a,b) ((a) > (b)+threshold)
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int M, N ;
+  const double* F_pt ;
+  int ndims ; 
+  int pdims = -1 ;
+  int* offsets ;
+  int* midx ;
+  int* neighbors ;
+  int nneighbors ;
+  int* dims ;
+  enum {F=0,THRESHOLD,P} ;
+  enum {MAXIMA=0} ;
+  double threshold = - mxGetInf() ; 
+
+
+  /* ------------------------------------------------------------------
+   *                                                Check the arguments
+   * --------------------------------------------------------------- */
+  if(nin > 1) {
+    if(!uIsRealScalar(in[THRESHOLD])) {
+      mexErrMsgTxt("THRESHOLD must be a real scalar.") ;
+    }
+    threshold = *mxGetPr(in[THRESHOLD]) ;
+  }
+    
+  ndims = mxGetNumberOfDimensions(in[F]) ;
+
+if(ndims !=3)
+    printf("Error Error: NDIMS not equal to 3. Not handled by C version\n");
+
+  {
+    /* We need to make a copy because in one special case (see below)
+       we need to adjust dims[].
+    */
+    int d ;
+    const int* const_dims = (int*) mxGetDimensions(in[F]) ;
+    dims = mxMalloc(sizeof(int)*ndims) ;
+    for(d=0 ; d < ndims ; ++d)
+    {
+/*        printf("Const Dimensions = %d\n", const_dims[d]);
+*/
+         dims[d] = const_dims[d] ;
+    }
+  }
+  M = dims[0] ;
+  N = dims[1] ;
+  F_pt = mxGetPr(in[F]) ;
+
+  /* 
+     If there are only two dimensions and if one is singleton, then
+     assume that a vector has been provided as input (and treat this
+     as a COLUMN matrix with p=1). We do this because Matlab does not
+     distinguish between vectors and 1xN or Mx1 matrices and because
+     the cases 1xN and Mx1 are trivial (the result is alway empty).
+  */
+  if((ndims == 2) && (pdims < 0) && (M == 1 || N == 1)) {
+
+    printf("ERROR ERROR: Entered a different loop here. Not handled by C version");
+
+    pdims = 1 ;
+    M = (M>N)?M:N ;
+    N = 1 ;
+    dims[0]=M ;
+    dims[1]=N ;
+  }
+
+  /* search the local maxima along the first p dimensions only */
+  if(pdims < 0)
+  {
+    pdims = ndims ;
+    }
+  
+  if(pdims > ndims) {
+    mxFree(dims) ;
+    mexErrMsgTxt("P must not be greater than the number of dimensions") ;
+  }
+
+  /* ------------------------------------------------------------------
+   *                                                         Do the job
+   * --------------------------------------------------------------- */  
+  {
+    int maxima_size = M*N ;
+    int* maxima_start = (int*) mxMalloc(sizeof(int) * maxima_size) ;
+    int* maxima_iterator = maxima_start ;
+    int* maxima_end = maxima_start + maxima_size ;
+    int i,j,h,o ;
+    const double* pt = F_pt ;
+    
+    /* Compute the offsets between dimensions. */
+    offsets = (int*) mxMalloc(sizeof(int) * ndims) ;
+    offsets[0] = 1 ;
+
+    for(h = 1 ; h < ndims ; ++h)
+    {
+      offsets[h] = offsets[h-1]*dims[h-1] ;
+/*      printf("%d:%d\t%d\n", h, offsets[h], dims[h-1]);
+*/ 
+    }
+
+    /* Multi-index. */
+    midx = (int*) mxMalloc(sizeof(int) * ndims) ;
+    for(h = 0 ; h < ndims ; ++h)
+      midx[h] = 1 ;
+
+    /* Neighbors. */
+    nneighbors = 1 ;
+    o=0 ;
+    for(h = 0 ; h < pdims ; ++h) {
+      nneighbors *= 3 ;
+      midx[h] = -1 ;
+      o -= offsets[h] ;
+    }
+    nneighbors -= 1 ;
+    neighbors = (int*) mxMalloc(sizeof(int) * nneighbors) ;
+
+    /* Precompute offsets from offset(-1,...,-1,0,...0) to
+     * offset(+1,...,+1,0,...,0). */
+    i = 0 ;
+
+    while(true) {
+      if(o != 0)
+        neighbors[i++] = o ;
+      h = 0 ;
+      while( o += offsets[h], (++midx[h]) > 1 ) {
+        o -= 3*offsets[h] ;
+        midx[h] = -1 ;
+        if(++h >= pdims)
+          goto stop ;
+      }
+    }
+  stop: ;
+        
+    /* Starts at the corner (1,1,...,1,0,0,...0) */
+    for(h = 0 ; h < pdims ; ++h) {
+      midx[h] = 1 ;
+      pt += offsets[h] ;
+/*      printf("%d:%x\t%d\n", h, pt, offsets[h]);
+*/
+    }
+  
+    for(h = pdims ; h < ndims ; ++h) {
+      midx[h] = 0 ;
+    }
+
+    /* ---------------------------------------------------------------
+     *                                                            Loop
+     * ------------------------------------------------------------ */
+
+    /*
+      If any dimension in the first P is less than 3 elements wide
+      then just return the empty matrix (if we proceed without doing
+      anything we break the carry reporting algorithm below).
+    */
+    for(h=0 ; h < pdims ; ++h)
+      if(dims[h] < 3) goto end ;
+    
+     while(true) {
+ 
+       /* Propagate carry along multi index midx */
+       h = 0 ;
+       while((midx[h]) >= dims[h] - 1) {
+/*         pt += 2*offsets[h] ;  skip first and last el. */
+         midx[h] = 1 ;
+         if(++h >= pdims) 
+           goto next_layer ;
+         ++midx[h] ;
+       }
+      
+      /*  Scan neighbors */
+      {
+        double v = *pt ;
+        bool is_greater = (v >= threshold) ;
+      /* printf("%f\n", v);
+        */
+        i = 0  ;
+        while(is_greater && i < nneighbors)
+        {
+          is_greater &= v > *(pt + neighbors[i++]) ;
+        }
+        
+        /* Add the local maximum */
+        if(is_greater) {
+          /* Need more space? */
+          if(maxima_iterator == maxima_end) {
+            maxima_size += M*N ;
+            maxima_start = (int*) mxRealloc(maxima_start, 
+                                            maxima_size*sizeof(int)) ;
+            maxima_end = maxima_start + maxima_size ;
+            maxima_iterator = maxima_end - M*N ;
+          }
+          
+          *maxima_iterator++ = pt - F_pt + 1 ;
+        }
+        
+        /* Go to next element */
+        pt += 1 ;
+        ++midx[0] ;
+        continue ;
+        
+      next_layer: ;
+        if( h >= ndims )
+          goto end ;
+        
+        while((++midx[h]) >= dims[h]) {
+          midx[h] = 0 ;
+          if(++h >= ndims)
+            goto end ;
+        }
+      }
+    }    
+  end:;
+    /* Return. */
+    {
+      double* M_pt ;
+/*      printf("%x\t%x\t%d\n", maxima_iterator, maxima_start, maxima_iterator-maxima_start);
+*/  
+      out[MAXIMA] = mxCreateDoubleMatrix
+        (1, maxima_iterator-maxima_start, mxREAL) ;
+      maxima_end = maxima_iterator ;
+      maxima_iterator = maxima_start ;
+      M_pt = mxGetPr(out[MAXIMA]) ;
+      while(maxima_iterator != maxima_end) 
+      {
+        *M_pt++ = *maxima_iterator++ ;
+      }
+    }
+    
+    /* Release space. */
+    mxFree(offsets) ;
+    mxFree(neighbors) ;
+    mxFree(midx) ;
+    mxFree(maxima_start) ;     
+  }
+  mxFree(dims) ;
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m
new file mode 100644
index 0000000..cee30ce
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m
@@ -0,0 +1,33 @@
+% SIFTLOCALMAX  Find local maximizers
+%   SEL=SIFTLOCALMAX(F) returns the indexes of the local maximizers of
+%   the Q-dimensional array F.
+%
+%   A local maximizer is an element whose value is greater than the
+%   value of all its neighbors.  The neighbors of an element i1...iQ
+%   are the subscripts j1...jQ such that iq-1 <= jq <= iq (excluding
+%   i1...iQ itself).  For example, if Q=1 the neighbors of an element
+%   are its predecessor and successor in the linear order; if Q=2, its
+%   neighbors are the elements immediately to its north, south, west,
+%   est, north-west, north-est, south-west and south-est
+%   (8-neighborhood).
+%
+%   Points on the boundary of F are ignored (and never selected as
+%   local maximizers).
+%
+%   SEL=SIFTLOCALMAX(F,THRESH) accepts an element as a mazimizer only
+%   if it is at least THRES greater than all its neighbors.
+%
+%   SEL=SIFTLOCALMAX(F,THRESH,P) look for neighbors only in the first
+%   P dimensions of the Q-dimensional array F. This is useful to
+%   process F in ``slices''.
+%
+%   REMARK.  Matrices (2-array) with a singleton dimension are
+%   interpreted as vectors (1-array). So for example SIFTLOCALMAX([0 1
+%   0]) and SIFTLOCALMAX([0 1 0]') both return 2 as an aswer. However,
+%   if [0 1 0] is to be interpreted as a 1x2 matrix, then the correct
+%   answer is the empty set, as all elements are on the boundary.
+%   Unfortunately MATLAB does not distinguish between vectors and
+%   2-matrices with a singleton dimension.  To forece the
+%   interpretation of all matrices as 2-arrays, use
+%   SIFTLOCALMAX(F,TRESH,2) (but note that in this case the result is
+%   always empty!).
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64
new file mode 100755
index 0000000..2b8e264
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx
new file mode 100755
index 0000000..613a74e
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c
new file mode 100644
index 0000000..1150176
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c
@@ -0,0 +1,250 @@
+/* file:        siftmatch.c
+** author:      Andrea Vedaldi
+** description: SIFT descriptor matching.
+**/
+
+/* AUTORIGHTS
+Copyright (c) 2006 The Regents of the University of California.
+All Rights Reserved.
+
+Created by Andrea Vedaldi
+UCLA Vision Lab - Department of Computer Science
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for educational, research and non-profit purposes,
+without fee, and without a written agreement is hereby granted,
+provided that the above copyright notice, this paragraph and the
+following three paragraphs appear in all copies.
+
+This software program and documentation are copyrighted by The Regents
+of the University of California. The software program and
+documentation are supplied "as is", without any accompanying services
+from The Regents. The Regents does not warrant that the operation of
+the program will be uninterrupted or error-free. The end-user
+understands that the program was developed for research purposes and
+is advised not to rely exclusively on the program for any reason.
+
+This software embodies a method for which the following patent has
+been issued: "Method and apparatus for identifying scale invariant
+features in an image and use of same for locating an object in an
+image," David G. Lowe, US Patent 6,711,293 (March 23,
+2004). Provisional application filed March 8, 1999. Asignee: The
+University of British Columbia.
+
+IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+*/
+
+#include"mexutils.c"
+
+#include<stdlib.h>
+#include<string.h>
+#include<math.h>
+
+#define greater(a,b) ((a) > (b))
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+#define TYPEOF_mxDOUBLE_CLASS double
+#define TYPEOF_mxSINGLE_CLASS float
+#define TYPEOF_mxINT8_CLASS   char
+#define TYPEOF_mxUINT8_CLASS  unsigned char
+
+#define PROMOTE_mxDOUBLE_CLASS double
+#define PROMOTE_mxSINGLE_CLASS float
+#define PROMOTE_mxINT8_CLASS   int
+#define PROMOTE_mxUINT8_CLASS  int
+
+#define MAXVAL_mxDOUBLE_CLASS mxGetInf()
+#define MAXVAL_mxSINGLE_CLASS ((float)mxGetInf())
+#define MAXVAL_mxINT8_CLASS   0x7fffffff
+#define MAXVAL_mxUINT8_CLASS  0x7fffffff
+
+typedef struct 
+{ 
+  int k1 ; 
+  int k2 ;
+  double score ;
+} Pair ;
+
+/* 
+ * This macro defines the matching function for abstract type; that
+ * is, it is a sort of C++ template.  This is also a good illustration
+ * of why C++ is preferable for templates :-)
+ */
+#define _COMPARE_TEMPLATE(MXC)                                          \
+  Pair* compare_##MXC (Pair* pairs_iterator,                            \
+                       const TYPEOF_##MXC * L1_pt,                      \
+                       const TYPEOF_##MXC * L2_pt,                      \
+                       int K1, int K2, int ND, float thresh)            \
+  {                                                                     \
+    int k1, k2 ;                                                        \
+    const PROMOTE_##MXC maxval = MAXVAL_##MXC ;                         \
+    for(k1 = 0 ; k1 < K1 ; ++k1, L1_pt += ND ) {                        \
+                                                                        \
+      PROMOTE_##MXC best = maxval ;                                     \
+      PROMOTE_##MXC second_best = maxval ;                              \
+      int bestk = -1 ;                                                  \
+                                                                        \
+      /* For each point P2[k2] in the second image... */                \
+      for(k2 =  0 ; k2 < K2 ; ++k2, L2_pt += ND) {                      \
+                                                                        \
+        int bin ;                                                       \
+        PROMOTE_##MXC acc = 0 ;                                         \
+        for(bin = 0 ; bin < ND ; ++bin) {                               \
+          PROMOTE_##MXC delta =                                         \
+            ((PROMOTE_##MXC) L1_pt[bin]) -                              \
+            ((PROMOTE_##MXC) L2_pt[bin]) ;                              \
+          acc += delta*delta ;                                          \
+        }                                                               \
+                                                                        \
+        /* Filter the best and second best matching point. */           \
+        if(acc < best) {                                                \
+          second_best = best ;                                          \
+          best = acc ;                                                  \
+          bestk = k2 ;                                                  \
+        } else if(acc < second_best) {                                  \
+          second_best = acc ;                                           \
+        }                                                               \
+      }                                                                 \
+                                                                        \
+      L2_pt -= ND*K2 ;                                                  \
+                                                                        \
+      /* Lowe's method: accept the match only if unique. */             \
+      if(thresh * (float) best <= (float) second_best &&                \
+         bestk != -1) {                                                 \
+        pairs_iterator->k1 = k1 ;                                       \
+        pairs_iterator->k2 = bestk ;                                    \
+        pairs_iterator->score = best ;                                  \
+        pairs_iterator++ ;                                              \
+      }                                                                 \
+    }                                                                   \
+                                                                        \
+    return pairs_iterator ;                                             \
+  }                                                                     \
+
+_COMPARE_TEMPLATE( mxDOUBLE_CLASS )
+_COMPARE_TEMPLATE( mxSINGLE_CLASS )
+_COMPARE_TEMPLATE( mxINT8_CLASS   )
+_COMPARE_TEMPLATE( mxUINT8_CLASS  )
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int K1, K2, ND ;
+  void* L1_pt ;
+  void* L2_pt ;
+  double thresh = 1.5 ;
+  mxClassID data_class ;
+  enum {L1=0,L2,THRESH} ;
+  enum {MATCHES=0,D} ;
+
+  /* ------------------------------------------------------------------
+  **                                                Check the arguments
+  ** --------------------------------------------------------------- */ 
+  if (nin < 2) {
+    mexErrMsgTxt("At least two input arguments required");
+  } else if (nout > 2) {
+    mexErrMsgTxt("Too many output arguments");
+  }
+    
+  if(!mxIsNumeric(in[L1]) ||
+     !mxIsNumeric(in[L2]) ||
+     mxGetNumberOfDimensions(in[L1]) > 2 ||
+     mxGetNumberOfDimensions(in[L2]) > 2) {
+    mexErrMsgTxt("L1 and L2 must be two dimensional numeric arrays") ;
+  }
+  
+  K1 = mxGetN(in[L1]) ;
+  K2 = mxGetN(in[L2]) ;
+  ND = mxGetM(in[L1]) ;
+  
+  if(mxGetM(in[L2]) != ND) {
+    mexErrMsgTxt("L1 and L2 must have the same number of rows") ;
+  }
+
+  data_class = mxGetClassID(in[L1]) ;
+  if(mxGetClassID(in[L2]) != data_class) {
+    mexErrMsgTxt("L1 and L2 must be of the same class") ;
+  }
+
+  L1_pt = mxGetData(in[L1]) ;
+  L2_pt = mxGetData(in[L2]) ;
+
+  if(nin == 3) {
+    if(!uIsRealScalar(in[THRESH])) {
+      mexErrMsgTxt("THRESH should be a real scalar") ;
+    }
+    thresh = *mxGetPr(in[THRESH]) ;
+  } else if(nin > 3) {
+    mexErrMsgTxt("At most three arguments are allowed") ;
+  }
+ 
+  /* ------------------------------------------------------------------
+  **                                                         Do the job
+  ** --------------------------------------------------------------- */ 
+  {
+    Pair* pairs_begin = (Pair*) mxMalloc(sizeof(Pair) * (K1+K2)) ;
+    Pair* pairs_iterator = pairs_begin ;
+
+
+#define _DISPATCH_COMPARE( MXC )                                        \
+    case MXC :                                                          \
+      pairs_iterator = compare_##MXC(pairs_iterator,                    \
+                                     (const TYPEOF_##MXC*) L1_pt,       \
+                                     (const TYPEOF_##MXC*) L2_pt,       \
+                                     K1,K2,ND,thresh) ;                 \
+    break ;                                                             \
+    
+    switch (data_class) {
+    _DISPATCH_COMPARE( mxDOUBLE_CLASS ) ;
+    _DISPATCH_COMPARE( mxSINGLE_CLASS ) ; 
+    _DISPATCH_COMPARE( mxINT8_CLASS   ) ;
+    _DISPATCH_COMPARE( mxUINT8_CLASS  ) ;
+    default :
+      mexErrMsgTxt("Unsupported numeric class") ;
+      break ;
+    }
+
+    /* ---------------------------------------------------------------
+     *                                                        Finalize
+     * ------------------------------------------------------------ */
+    {
+      Pair* pairs_end = pairs_iterator ;
+      double* M_pt ;
+      double* D_pt = NULL ;
+
+      out[MATCHES] = mxCreateDoubleMatrix
+        (2, pairs_end-pairs_begin, mxREAL) ;
+
+      M_pt = mxGetPr(out[MATCHES]) ;
+
+      if(nout > 1) {
+        out[D] = mxCreateDoubleMatrix(1,
+                                      pairs_end-pairs_begin,
+                                      mxREAL) ;
+        D_pt = mxGetPr(out[D]) ;
+      }
+
+      for(pairs_iterator = pairs_begin ; 
+          pairs_iterator < pairs_end  ; 
+          ++pairs_iterator) {
+        *M_pt++ = pairs_iterator->k1 + 1 ;
+        *M_pt++ = pairs_iterator->k2 + 1 ;
+        if(nout > 1) {
+          *D_pt++ = pairs_iterator->score ;
+        }
+      }
+    }
+    mxFree(pairs_begin) ;
+  }
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m
new file mode 100644
index 0000000..d4398c1
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m
@@ -0,0 +1,60 @@
+% SIFTMATCH  Match SIFT features
+%   MATCHES=SIFTMATCH(DESCR1, DESCR2) matches the two sets of SIFT
+%   descriptors DESCR1 and DESCR2.
+%
+%   The function uses the same algorithm suggested by D. Lowe [1] to
+%   reject matches that are too ambiguous.
+%
+%   SIFTMATCH(DESCR1, DESCR2, THRESH) uses [1] with the specified
+%   threshold THRESH. A descriptor D1 is matched to a descriptor D2
+%   only if the distance d(D1,D2) multiplied by THRESH is not greather
+%   than the distance of D1 to all other descriptors. The default
+%   value of THRESH is 1.5.
+%   
+%   The storage class of the descriptors can be either DOUBLE, FLOAT,
+%   INT8 or UINT8. Usually interger classes are faster.
+%   
+%   [1] D. G. Lowe, 
+%       `Distinctive image features from scale-invariant keypoints,'
+%       IJCV, vol. 2, no. 60, pp. 91–110, 2004.
+%
+%   See also SIFT(), SIFTDESCRIPTOR().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64
new file mode 100755
index 0000000..e58eaf5
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx
new file mode 100755
index 0000000..d5e32f7
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.c b/SD-VBS/benchmarks/sift/src/matlab/siftormx.c
new file mode 100644
index 0000000..40380a2
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftormx.c
@@ -0,0 +1,251 @@
+/* file:        siftormx.c
+** author:      Andrea Vedaldi
+** description: Computes peaks of orientation histogram.
+**/
+
+/* AUTORIGHTS
+Copyright (c) 2006 The Regents of the University of California.
+All Rights Reserved.
+
+Created by Andrea Vedaldi
+UCLA Vision Lab - Department of Computer Science
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for educational, research and non-profit purposes,
+without fee, and without a written agreement is hereby granted,
+provided that the above copyright notice, this paragraph and the
+following three paragraphs appear in all copies.
+
+This software program and documentation are copyrighted by The Regents
+of the University of California. The software program and
+documentation are supplied "as is", without any accompanying services
+from The Regents. The Regents does not warrant that the operation of
+the program will be uninterrupted or error-free. The end-user
+understands that the program was developed for research purposes and
+is advised not to rely exclusively on the program for any reason.
+
+This software embodies a method for which the following patent has
+been issued: "Method and apparatus for identifying scale invariant
+features in an image and use of same for locating an object in an
+image," David G. Lowe, US Patent 6,711,293 (March 23,
+2004). Provisional application filed March 8, 1999. Asignee: The
+University of British Columbia.
+
+IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+*/
+
+#include"mex.h"
+#include<stdlib.h>
+#include<string.h>
+#include<math.h>
+
+#include<mexutils.c>
+
+#define greater(a,b) a > b
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+const double win_factor = 1.5 ;
+#define NBINS 36
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int M,N,S,smin,K ;
+  const int* dimensions ;
+  const double* P_pt ;
+  const double* G_pt ;
+  double* TH_pt ; 
+  double sigma0 ;
+  double H_pt [ NBINS ] ;
+
+  enum {IN_P=0,IN_G,IN_S,IN_SMIN,IN_SIGMA0} ;
+  enum {OUT_Q=0} ;
+
+  /* -----------------------------------------------------------------
+  **                                               Check the arguments
+  ** -------------------------------------------------------------- */ 
+  if (nin != 5) {
+    mexErrMsgTxt("Exactly five input arguments required.");
+  } else if (nout > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  if( !uIsRealScalar(in[IN_S]) ) {
+    mexErrMsgTxt("S should be a real scalar") ;
+  }
+                
+  if( !uIsRealScalar(in[IN_SMIN]) ) {
+    mexErrMsgTxt("SMIN should be a real scalar") ;
+  }
+        
+  if( !uIsRealScalar(in[IN_SIGMA0]) ) {
+    mexErrMsgTxt("SIGMA0 should be a real scalar") ;
+  }
+        
+  if( !uIsRealMatrix(in[IN_P],3,-1)) {
+    mexErrMsgTxt("P should be a 3xK real matrix") ;
+  }
+        
+  if(mxGetNumberOfDimensions(in[IN_G]) != 3) {
+    mexErrMsgTxt("SSO must be a three dimensional array") ;
+  }
+
+  dimensions = mxGetDimensions(in[IN_G]) ;
+  M = dimensions[0] ;
+  N = dimensions[1] ;
+  S = (int)(*mxGetPr(in[IN_S])) ;
+  smin = (int)(*mxGetPr(in[IN_SMIN])) ;
+  sigma0 = *mxGetPr(in[IN_SIGMA0]) ;
+        
+  K = mxGetN(in[IN_P]) ;
+  P_pt = mxGetPr(in[IN_P]) ;
+  G_pt = mxGetPr(in[IN_G]) ;
+        
+
+  /* If the input array is empty, then output an empty array as well. */
+  if(K == 0) {    
+    out[OUT_Q] = mxCreateDoubleMatrix(4,0,mxREAL) ;
+    return ;
+  }
+
+  /* ------------------------------------------------------------------
+   *                                                         Do the job
+   * --------------------------------------------------------------- */ 
+  {
+    int p ;
+    const int yo = 1 ;
+    const int xo = M ;
+    const int so = M*N ;
+
+    int buffer_size = K*4 ;
+    double* buffer_start = (double*) mxMalloc( buffer_size *sizeof(double)) ;
+    double* buffer_iterator = buffer_start ;
+    double* buffer_end = buffer_start + buffer_size ;
+
+    for(p = 0 ; p < K ; ++p, TH_pt += 2) {
+      const double x = *P_pt++ ;
+      const double y = *P_pt++ ;
+      const double s = *P_pt++ ;
+      int xi = ((int) (x+0.5)) ; /* Round them off. */
+      int yi = ((int) (y+0.5)) ;
+      int si = ((int) (s+0.5)) - smin ;
+      int xs ;
+      int ys ;
+      double sigmaw = win_factor * sigma0 * pow(2, ((double)s) / S) ;
+      int W = (int) ceil(3.0 * sigmaw) ;
+      int bin ;
+      const double* pt ;
+
+      /* Make sure that the rounded off keypoint index is within bound.
+       */
+      if(xi < 0   || 
+         xi > N-1 || 
+         yi < 0   || 
+         yi > M-1 || 
+         si < 0   || 
+         si > dimensions[2]-1 ) {
+        mexPrintf("Dropping %d: W %d x %d y %d si [%d,%d,%d,%d]\n",p,W,xi,yi,si,M,N,dimensions[2]) ;
+        continue ;
+      }
+      
+      /* Clear histogram buffer. */
+      {
+        int i ;
+        for(i = 0 ; i < NBINS ; ++i) 
+          H_pt[i] = 0 ;
+      }
+
+      pt = G_pt + xi*xo + yi*yo + si*so ;      
+
+#define at(dx,dy) (*(pt + (dx)*xo + (dy)*yo))
+
+      for(xs = max(-W, 1-xi) ; xs <= min(+W, N -2 -xi) ; ++xs) {
+        for(ys = max(-W, 1-yi) ; ys <= min(+W, M -2 -yi) ; ++ys) {
+          double Dx = 0.5 * ( at(xs+1,ys) - at(xs-1,ys) ) ;
+          double Dy = 0.5 * ( at(xs,ys+1) - at(xs,ys-1) ) ;
+          double dx = ((double)(xi+xs)) - x;
+          double dy = ((double)(yi+ys)) - y;
+                                        
+          if(dx*dx + dy*dy > W*W+0.5) continue ;
+
+          {
+            double win = exp( - (dx*dx + dy*dy)/(2*sigmaw*sigmaw) ) ;
+            double mod = sqrt(Dx*Dx + Dy*Dy) ;
+            double theta = fmod(atan2(Dy, Dx) + 2*M_PI, 2*M_PI) ; 
+            bin = (int) floor( NBINS * theta / (2*M_PI) ) ;
+            H_pt[bin] += mod*win ;        
+          }
+        }
+      }
+                        
+      /* Smooth histogram */                    
+      {
+        int iter, i ;
+        for (iter = 0; iter < 6; iter++) {
+          double prev;
+          prev = H_pt[NBINS-1];
+          for (i = 0; i < NBINS; i++) {
+            float newh = (prev + H_pt[i] + H_pt[(i+1) % NBINS]) / 3.0;
+            prev = H_pt[i] ;
+            H_pt[i] = newh ;
+          }
+        }
+      }
+                        
+      /* Find most strong peaks. */
+      {
+        int i ;
+        double maxh = H_pt[0] ;
+        for(i = 1 ; i < NBINS ; ++i)
+          maxh = max(maxh, H_pt[i]) ;
+                                
+        for(i = 0 ; i < NBINS ; ++i) {
+          double h0 = H_pt[i] ;
+          double hm = H_pt[(i-1+NBINS) % NBINS] ;
+          double hp = H_pt[(i+1+NBINS) % NBINS] ;
+          
+          if( h0 > 0.8*maxh && h0 > hm && h0 > hp ) {
+
+            double di = -0.5 * (hp-hm) / (hp+hm-2*h0) ; /*di=0;*/
+            double th = 2*M_PI*(i+di+0.5)/NBINS ;
+
+            if( buffer_iterator == buffer_end ) {
+              int offset = buffer_iterator - buffer_start ;
+              buffer_size += 4*max(1, K/16) ;
+              buffer_start = (double*) mxRealloc(buffer_start,
+                                                 buffer_size*sizeof(double)) ;
+              buffer_end = buffer_start + buffer_size ;
+              buffer_iterator = buffer_start + offset ;
+            }
+            
+            *buffer_iterator++ = x ;
+            *buffer_iterator++ = y ;
+            *buffer_iterator++ = s ;
+            *buffer_iterator++ = th ;
+          }
+        } /* Scan histogram */
+      } /* Find peaks */
+    }
+
+    /* Save back the result. */
+    {
+      double* result ;
+      int NL = (buffer_iterator - buffer_start)/4 ;
+      out[OUT_Q] = mxCreateDoubleMatrix(4, NL, mxREAL) ;
+      result  = mxGetPr(out[OUT_Q]);
+      memcpy(result, buffer_start, sizeof(double) * 4 * NL) ;
+    }
+    mxFree(buffer_start) ;
+  }
+}
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64
new file mode 100755
index 0000000..ae1f787
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx
new file mode 100755
index 0000000..de90e4b
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftread.m b/SD-VBS/benchmarks/sift/src/matlab/siftread.m
new file mode 100644
index 0000000..76945c8
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftread.m
@@ -0,0 +1,101 @@
+function [frames,descriptors] = siftread(file) 
+% SIFTREAD Read Lowe's SIFT implementation data files
+%   [FRAMES, DESCRIPTORS] = READSIFT(FILE) reads the frames and the
+%   descriptors from the specified file. The function reads files
+%   produced by Lowe's SIFT implementation.
+%     
+%   FRAMES and DESCRIPTORS have the same format used by SIFT(). 
+%
+%   REMARK. Lowe's and our implementations use a silightly different
+%   convention to store the orientation of the frame. When the file
+%   is read, the orientation is changed to match our convention.
+%
+%   See also SIFT().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+verbosity=0 ;
+
+g = fopen(file, 'r');
+if g == -1
+    error(['Could not open file ''', file, '''.']) ;
+end
+[header, count] = fscanf(g, '%d', [1 2]) ;
+if count ~= 2
+    error('Invalid keypoint file header.');
+end
+K = header(1) ; 
+DL = header(2) ;
+
+if(verbosity > 0)
+        fprintf('%d keypoints, %d descriptor length.\n', K, DL) ;
+end
+
+%creates two output matrices
+P = zeros(4,K) ;
+L = zeros(DL,K) ;
+
+%parse tmp.key
+for k = 1:K
+
+        % Record format: i,j,s,th
+        [record, count] = fscanf(g, '%f', [1 4]) ; 
+        if count ~= 4
+                error(...
+                        sprintf('Invalid keypoint file (parsing keypoint %d)',k) );
+        end
+        P(:,k) = record(:) ;
+        
+        % Record format: descriptor
+        [record, count] = fscanf(g, '%d', [1 DL]) ;
+        if count ~= DL
+                error(...
+                        sprintf('Invalid keypoint file (parsing keypoint %d)',k) );
+        end
+        L(:,k) = record(:) ;
+        
+end
+fclose(g) ;
+
+L=double(L) ;
+P(1:2,:)=flipud(P(1:2,:)) ; % i,j -> x,y
+
+frames=[ P(1:2,:) ; P(3,:) ; -P(4,:) ] ;
+descriptors = L ;
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c
new file mode 100644
index 0000000..a1ef8c7
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c
@@ -0,0 +1,342 @@
+/* file:        siftrefinemx.c
+
+Benchmark               - sift
+Data set                - sqcif
+
+                                                                                < M A T L A B >
+                                                                    Copyright 1984-2006 The MathWorks, Inc.
+                                                                          Version 7.3.0.298 (R2006b)
+                                                                                August 03, 2006
+
+ 
+  To get started, type one of these: helpwin, helpdesk, or demo.
+  For product information, visit www.mathworks.com.
+ 
+Warning: Function /h/g2/kvs/checkParallel/sdvbs-svn/common/matlab/randn.m has the same name as a MATLAB builtin. We suggest you rename the function to avoid a potential name conflict.
+> In path at 113
+  In script_run_profile at 3
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Warning: You are using gcc version "4.1.1".  The earliest gcc version supported
+with mex is "3.4.0".  The latest version tested for use with mex is "3.4.5".
+To download a different version of gcc, visit http://gcc.gnu.org 
+Input size              - (96x96)
+** author:      Andrea Vedaldi
+** description: Subpixel localization, thresholding and off-edge test.
+**/
+
+/* AUTORIGHTS
+Copyright (c) 2006 The Regents of the University of California.
+All Rights Reserved.
+
+Created by Andrea Vedaldi
+UCLA Vision Lab - Department of Computer Science
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for educational, research and non-profit purposes,
+without fee, and without a written agreement is hereby granted,
+provided that the above copyright notice, this paragraph and the
+following three paragraphs appear in all copies.
+
+This software program and documentation are copyrighted by The Regents
+of the University of California. The software program and
+documentation are supplied "as is", without any accompanying services
+from The Regents. The Regents does not warrant that the operation of
+the program will be uninterrupted or error-free. The end-user
+understands that the program was developed for research purposes and
+is advised not to rely exclusively on the program for any reason.
+
+This software embodies a method for which the following patent has
+been issued: "Method and apparatus for identifying scale invariant
+features in an image and use of same for locating an object in an
+image," David G. Lowe, US Patent 6,711,293 (March 23,
+2004). Provisional application filed March 8, 1999. Asignee: The
+University of British Columbia.
+
+IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+*/
+
+#include"mex.h"
+
+#include<mexutils.c>
+
+#include<stdlib.h>
+#include<string.h>
+
+/* Prototype of DGESV LAPACK function for the solution of a linear system. */
+#ifdef WINDOWS
+#define DGESV dgesv
+#undef min
+#undef max
+#else
+#define DGESV dgesv_
+#endif
+
+#ifdef WINDOWS
+#ifdef __cplusplus__
+extern "C" {
+    extern int DGESV(int *n, int *nrhs, double *a, int *lda,
+                     int *ipiv, double *b, int *ldb, int *info) ;
+}
+#else
+  extern int DGESV(int *n, int *nrhs, double *a, int *lda,
+                     int *ipiv, double *b, int *ldb, int *info) ;
+#define sqrtf(x) ((float)sqrt((double)x)
+#define powf(x)  ((float)pow((double)x)
+#define fabsf(x) ((float)fabs((double)x)
+#endif
+#else
+extern int DGESV(int *n, int *nrhs, double *a, int *lda,
+                 int *ipiv, double *b, int *ldb, int *info) ;
+#endif
+
+#define greater(a,b) ((a) > (b))
+#define min(a,b) (((a)<(b))?(a):(b))
+#define max(a,b) (((a)>(b))?(a):(b))
+
+const int max_iter = 5 ;
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+  int M,N,S,smin,K ;
+  const int* dimensions ;
+  const double* P_pt ;
+  const double* D_pt ; 
+  double threshold = 0.01 ; /*0.02 ;*/
+  double r = 10.0 ;
+  double* result ;
+  enum {IN_P=0,IN_D,IN_SMIN,IN_THRESHOLD,IN_R} ;
+  enum {OUT_Q=0} ;
+        
+  /* -----------------------------------------------------------------
+  **                                               Check the arguments
+  ** -------------------------------------------------------------- */ 
+  if (nin < 3) {
+    mexErrMsgTxt("At least three input arguments required.");
+  } else if (nout > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  if( !uIsRealMatrix(in[IN_P],3,-1) ) {
+    mexErrMsgTxt("P must be a 3xK real matrix") ;
+  }
+
+  if( !mxIsDouble(in[IN_D]) || mxGetNumberOfDimensions(in[IN_D]) != 3) {
+    mexErrMsgTxt("G must be a three dimensional real array.") ;
+  }
+
+  if( !uIsRealScalar(in[IN_SMIN]) ) {
+    mexErrMsgTxt("SMIN must be a real scalar.") ;
+  }
+
+  if(nin >= 4) {
+    if(!uIsRealScalar(in[IN_THRESHOLD])) {
+      mexErrMsgTxt("THRESHOLD must be a real scalar.") ;
+    }
+    threshold = *mxGetPr(in[IN_THRESHOLD]) ;
+  }
+  
+  if(nin >= 5) {
+    if(!uIsRealScalar(in[IN_R])) {
+      mexErrMsgTxt("R must be a real scalar.") ;
+    }
+    r = *mxGetPr(in[IN_R]) ;    
+  }
+  
+  dimensions = mxGetDimensions(in[IN_D]) ;
+  M = dimensions[0] ;
+  N = dimensions[1] ;
+  S = dimensions[2] ; 
+  smin = (int)(*mxGetPr(in[IN_SMIN])) ;
+  
+  if(S < 3 || M < 3 || N < 3) {
+    mexErrMsgTxt("All dimensions of DOG must be not less than 3.") ;
+  }
+  
+  K = mxGetN(in[IN_P]) ;  
+  P_pt = mxGetPr(in[IN_P]) ;
+  D_pt = mxGetPr(in[IN_D]) ;
+
+  /* If the input array is empty, then output an empty array as well. */
+  if( K == 0) {
+    out[OUT_Q] = mxDuplicateArray(in[IN_P]) ;
+    return ;
+  }
+
+  /* -----------------------------------------------------------------
+   *                                                        Do the job
+   * -------------------------------------------------------------- */
+  {    
+    double* buffer = (double*) mxMalloc(K*3*sizeof(double)) ;
+    double* buffer_iterator = buffer ;
+    int p ;
+    const int yo = 1 ;
+    const int xo = M ;
+    const int so = M*N ;
+    
+/*    printf("Actual values = %d\n\n", K);
+*/  
+    for(p = 0 ; p < K ; ++p) {
+      int x = ((int)*P_pt++) ;
+      int y = ((int)*P_pt++) ;
+/*      printf("%d\t%d\n", ((int)*P_pt), smin);
+*/  
+      int s = ((int)*P_pt++) - smin ;
+      int iter ;
+      double b[3] ;
+        
+      /* Local maxima extracted from the DOG
+       * have coorrinates 1<=x<=N-2, 1<=y<=M-2
+       * and 1<=s-mins<=S-2. This is also the range of the points
+       * that we can refine.
+       */
+/*      printf("%d\t%d\t%d\t%d\t%d\t%d\n", x, N-2,y,M-2, s, S-2);
+*/  
+      if(x < 1 || x > N-2 ||
+         y < 1 || y > M-2 ||
+         s < 1 || s > S-2) {
+        continue ;
+      }
+
+#define at(dx,dy,ds) (*(pt + (dx)*xo + (dy)*yo + (ds)*so))
+
+      {
+        const double* pt = D_pt + y*yo + x*xo + s*so ;      
+  
+        double Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
+        int dx = 0 ;
+        int dy = 0 ;
+/*        printf("%d\t%d\t%d\t%d\t%d\t%d\t%d\t%f\t%f\n",S, y, yo, x, xo, s, so, *D_pt, *pt); 
+*/        
+        for(iter = 0 ; iter < max_iter ; ++iter) {
+
+          double A[3*3] ;          
+          int ipiv[3] ;
+          int n = 3 ;
+          int one = 1 ;
+          int info = 0 ;
+
+#define Aat(i,j) (A[(i)+(j)*3])    
+
+          x += dx ;
+          y += dy ;
+          pt = D_pt + y*yo + x*xo + s*so ;
+          
+          /* Compute the gradient. */
+          Dx = 0.5 * (at(+1,0,0) - at(-1,0,0)) ;
+          Dy = 0.5 * (at(0,+1,0) - at(0,-1,0));
+          Ds = 0.5 * (at(0,0,+1) - at(0,0,-1)) ;
+          
+          /* Compute the Hessian. */
+          Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0 * at(0,0,0)) ;
+          Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0 * at(0,0,0)) ;
+          Dss = (at(0,0,+1) + at(0,0,-1) - 2.0 * at(0,0,0)) ;
+          
+          Dxy = 0.25 * ( at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0) ) ;
+          Dxs = 0.25 * ( at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1) ) ;
+          Dys = 0.25 * ( at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1) ) ;
+          
+          /* Solve linear system. */
+          Aat(0,0) = Dxx ;
+          Aat(1,1) = Dyy ;
+          Aat(2,2) = Dss ;
+          Aat(0,1) = Aat(1,0) = Dxy ;
+          Aat(0,2) = Aat(2,0) = Dxs ;
+          Aat(1,2) = Aat(2,1) = Dys ;
+          
+          b[0] = - Dx ;
+          b[1] = - Dy ;
+          b[2] = - Ds ;
+          
+   /*       DGESV (&n, &one, A, &n, ipiv, b, &n, &info) ;
+    */      
+          /* If the translation of the keypoint is big, move the keypoint
+           * and re-iterate the computation. Otherwise we are all set.
+           */
+          dx= ((b[0] >  0.6 && x < N-2) ?  1 : 0 )
+            + ((b[0] < -0.6 && x > 1  ) ? -1 : 0 ) ;
+          
+          dy= ((b[1] >  0.6 && y < M-2) ?  1 : 0 )
+            + ((b[1] < -0.6 && y > 1  ) ? -1 : 0 ) ;
+          
+          if( dx == 0 && dy == 0 ) break ;
+          
+        }
+                                
+        {
+          double val = at(0,0,0) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ; 
+          double score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ; 
+          double xn = x + b[0] ;
+          double yn = y + b[1] ;
+          double sn = s + b[2] ;
+         
+/*        printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", fabs(val),threshold,score,(r+1)*(r+1)/r,fabs(b[0]), fabs(b[1]), fabs(b[2]),xn,yn,sn,r);
+*/          
+          if(fabs(val) > threshold &&
+             score < (r+1)*(r+1)/r && 
+             score >= 0 &&
+             fabs(b[0]) < 1.5 &&
+             fabs(b[1]) < 1.5 &&
+             fabs(b[2]) < 1.5 &&
+             xn >= 0 &&
+             xn <= N-1 &&
+             yn >= 0 &&
+             yn <= M-1 &&
+             sn >= 0 &&
+             sn <= S-1) {
+            *buffer_iterator++ = xn ;
+            *buffer_iterator++ = yn ;
+            *buffer_iterator++ = sn+smin  ;
+          }
+        }
+      }
+    }      
+
+    /* Copy the result into an array. */
+    {
+        int i;
+      int NL = (buffer_iterator - buffer)/3 ;
+/*      printf("%NL VALUE = %d\t%d\t%d\n", NL, buffer_iterator, buffer);
+*/
+      out[OUT_Q] = mxCreateDoubleMatrix(3, NL, mxREAL) ;
+      result = mxGetPr(out[OUT_Q]);
+        for(i=0; i<(3*NL); i++)
+        {
+            result[i] = buffer[i];
+/*            printf("%f\t", buffer[i]);
+*/  
+        }
+/*        printf("\n");
+      memcpy(result, buffer, sizeof(double) * 3 * NL) ;
+*/
+    }
+    mxFree(buffer) ;
+  }
+  
+}
+
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m
new file mode 100644
index 0000000..222d610
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m
@@ -0,0 +1,61 @@
+% SIFTREFINEMX  Subpixel localization, thresholding and on-edge test
+%   Q = SIFTREFINEMX(P, OCTAVE, SMIN) refines, thresholds and performs
+%   the on-edge test for the SIFT frames P extracted from the DOG
+%   octave OCTAVE with parameter SMIN (see GAUSSIANSS()).
+%
+%   Q = SIFTREFINEMX(P, OCTAVE, SMIN, THRESH, R) specifies custom
+%   values for the local maximum strength threshold THRESH and the
+%   local maximum peakedeness threshold R.
+%
+%   OCTAVE is an octave of the Difference Of Gaussian scale space. P
+%   is a 3xK matrix specifying the indexes (X,Y,S) of the points of
+%   extremum of the octave OCTAVE. The spatial indexes X,Y are integer
+%   with base zero. The scale index S is integer with base SMIN and
+%   represents a scale sublevel in the specified octave.
+%
+%   The function returns a matrix Q containing the refined keypoints.
+%   The matrix has the same format as P, except that the indexes are
+%   now fractional. The function drops the points that do not satisfy
+%   the strength and peakedness tests.
+%
+%   See also SIFT().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64
new file mode 100755
index 0000000..c91556b
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx
new file mode 100755
index 0000000..74b312c
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx
diff --git a/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m b/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m
new file mode 100644
index 0000000..b790ef0
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m
@@ -0,0 +1,151 @@
+function H = tightsubplot(varargin)
+% TIGHTSUBPLOT  Tiles axes without wasting space
+%   H = TIGHTSUBPLOT(K,P) returns an handle to the P-th axis in a
+%   regular grid of K axes. The K axes are numbered from left to right
+%   and from top to bottom.  The function operates similarly to
+%   SUBPLOT(), but by default it does not put any margin between axes.
+%
+%   H = TIGHTSUBPLOT(M,N,P) retursn an handle to the P-th axes in a
+%   regular subdivision with M rows and N columns.
+%
+%   The function accepts the following option-value pairs:
+%
+%   'Spacing' [0]
+%     Set extra spacing between axes.  The space is added between the
+%     inner or outer boxes, depending on the setting below.
+%
+%   'Box' ['inner'] (** ONLY >R14 **)
+%     If set to 'outer', the function displaces the axes by their
+%     outer box, thus protecting title and labels. Unfortunately
+%     MATLAB typically picks unnecessarily large insets, so that a bit
+%     of space is wasted in this case.  If set to 'inner', the
+%     function uses the inner box. This causes the instets of nearby
+%     axes to overlap, but it is very space conservative.
+%
+%   REMARK. While SUBPLOT kills any pre-existing axes that overalps a
+%   new one, this function does not.
+%
+%   See also SUBPLOT().
+
+% AUTORIGHTS
+% Copyright (c) 2006 The Regents of the University of California.
+% All Rights Reserved.
+% 
+% Created by Andrea Vedaldi
+% UCLA Vision Lab - Department of Computer Science
+% 
+% Permission to use, copy, modify, and distribute this software and its
+% documentation for educational, research and non-profit purposes,
+% without fee, and without a written agreement is hereby granted,
+% provided that the above copyright notice, this paragraph and the
+% following three paragraphs appear in all copies.
+% 
+% This software program and documentation are copyrighted by The Regents
+% of the University of California. The software program and
+% documentation are supplied "as is", without any accompanying services
+% from The Regents. The Regents does not warrant that the operation of
+% the program will be uninterrupted or error-free. The end-user
+% understands that the program was developed for research purposes and
+% is advised not to rely exclusively on the program for any reason.
+% 
+% This software embodies a method for which the following patent has
+% been issued: "Method and apparatus for identifying scale invariant
+% features in an image and use of same for locating an object in an
+% image," David G. Lowe, US Patent 6,711,293 (March 23,
+% 2004). Provisional application filed March 8, 1999. Asignee: The
+% University of British Columbia.
+% 
+% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
+% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
+% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
+% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
+% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
+% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
+% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
+% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+sp=0.0 ;
+use_outer=0 ;
+
+% --------------------------------------------------------------------
+%                                                      Parse arguments
+% --------------------------------------------------------------------
+K=varargin{1} ;
+p=varargin{2} ;
+N = ceil(sqrt(K)) ;
+M = ceil(K/N) ;
+
+a=3 ;
+NA = length(varargin) ;
+if NA > 2
+  if isa(varargin{3},'char')
+    % Called with K and p
+  else
+    % Called with M,N and p
+    a = 4 ;
+    M = K ;
+    N = p ;
+    p = varargin{3} ;
+  end
+end
+
+for a=a:2:NA
+  switch varargin{a}
+    case 'Spacing'
+      sp=varargin{a+1} ;
+    case 'Box'      
+      switch varargin{a+1}
+        case 'inner'
+          use_outer = 0 ;
+        case 'outer'
+        if ~strcmp(version('-release'), '14')
+          %warning(['Box option supported only on MATALB 14']) ;
+          continue;
+        end
+        use_outer = 1 ;
+        otherwise
+          error(['Box is either ''inner'' or ''outer''']) ;
+      end
+    otherwise
+      error(['Uknown parameter ''', varargin{a}, '''.']) ;
+  end      
+end
+
+% --------------------------------------------------------------------
+%                                                  Check the arguments
+% --------------------------------------------------------------------
+
+[j,i]=ind2sub([N M],p) ;
+i=i-1 ;
+j=j-1 ;
+
+dt = sp/2 ;
+db = sp/2 ;
+dl = sp/2 ;
+dr = sp/2 ;
+
+pos = [  j*1/N+dl,...
+       1-i*1/M-1/M-db,...
+       1/N-dl-dr,...
+       1/M-dt-db] ;
+
+switch use_outer
+  case 0
+    H = findobj(gcf, 'Type', 'axes', 'Position', pos) ;
+    if(isempty(H))
+      H = axes('Position', pos) ;
+    else
+      axes(H) ;
+    end
+    
+  case 1
+    H = findobj(gcf, 'Type', 'axes', 'OuterPosition', pos) ;
+    if(isempty(H))
+      H = axes('ActivePositionProperty', 'outerposition',...
+               'OuterPosition', pos) ;
+    else
+      axes(H) ;
+    end
+end