52 files changed, 5480 insertions, 0 deletions
diff --git a/SD-VBS/benchmarks/sift/src/c/diffss.c b/SD-VBS/benchmarks/sift/src/c/diffss.c
new file mode 100644
index 0000000..683cbe7
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/diffss.c
@@ -0,0 +1,53 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include "sift.h"
+/**
+    DIFFSS  Difference of scale space
+    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.
+**/
+F2D** diffss(F2D** ss, int num, int intervals)
+{
+    F2D** dss;
+    int o, sizeM, sizeN, s, i, j; 
+    F2D *current, *in1, *in2;
+    dss = malloc(num*intervals*sizeof(F2D*));
+    for(o=0; o<num; o++)
+    {
+        for(s=0; s<(intervals-1); s++)
+        {
+            sizeM = ss[o*intervals+s]->height;
+            sizeN = ss[o*intervals+s]->width;
+            dss[o*intervals+s] = fMallocHandle(sizeM, sizeN);
+            current = dss[o*intervals+s];
+            in1 = ss[o*intervals+s+1];
+            in2 = ss[o*intervals+s];
+            
+            for(i=0; i<sizeM; i++)
+            {
+                for(j=0; j<sizeN; j++)
+                {
+                    subsref(current,i,j) = subsref(in1,i,j) - subsref(in2,i,j);
+                }
+            }
+        }
+    }
+    return dss;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/doubleSize.c b/SD-VBS/benchmarks/sift/src/c/doubleSize.c
new file mode 100644
index 0000000..dd17f77
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/doubleSize.c
@@ -0,0 +1,55 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include "sift.h"
+F2D* doubleSize(F2D* I)
+{
+    F2D *J;
+    int M, N, i, j;
+    M = I->height;
+    N = I->width;
+    J = fSetArray(2*M, 2*N, 0);
+    for(i=0; i<M; i++)
+    {
+        for(j=0; j<N; j++)
+        {
+            subsref(J,2*i,j*2) = subsref(I,i,j);
+        }
+    }
+    
+    for(i=0; i<M-1; i++)
+    {
+        for(j=0; j<N-1; j++)
+        {
+            subsref(J,i*2+1,j*2+1) = (0.25 * (subsref(I,i,j) + subsref(I,i+1,j) + subsref(I,i,j+1) + subsref(I,(i+1),j+1)));
+        }
+    }
+    for(i=0; i<M-1; i++)
+    {
+        for(j=0; j<N; j++)
+        {
+            subsref(J,i*2+1,j*2) = (0.5 * (subsref(I,i,j) + subsref(I,i+1,j)));
+        }
+    }
+   
+    for(i=0; i<M; i++)
+    {
+        for(j=0; j<N-1; j++)
+        {
+            subsref(J,i*2,j*2+1) = (0.5 * (subsref(I,i,j) + subsref(I,i,j+1)));
+        }
+    }
+    
+    return J;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/filterBoundaryPoints.c b/SD-VBS/benchmarks/sift/src/c/filterBoundaryPoints.c
new file mode 100644
index 0000000..c68e717
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/filterBoundaryPoints.c
@@ -0,0 +1,61 @@
+#include "sift.h"
+/**
+    Filter out points on the boundaries. 
+    The function returns oframes, which
+    contains the row, col and the interval number 
+**/
+F2D* filterBoundaryPoints(int M, int N, F2D* oframes) 
+{
+    int i, k=0, m=0;
+        int cnt;
+    I2D* sel;
+    F2D *ret;
+    
+    cnt = 0;
+    for(i=0; i<oframes->width; i++)
+        {
+        if(asubsref(oframes,i)>3 && asubsref(oframes,i)<(N-3) && 
+                        subsref(oframes,1,i)>3 && subsref(oframes,1,i)<(M-3))
+                {
+            cnt++;
+                }
+        }
+    sel = iSetArray(cnt, 1, 0);
+    for(i=0; i<oframes->width; i++)
+    {
+        if(asubsref(oframes,i)>3 && asubsref(oframes,i)<(N-3) && 
+                        subsref(oframes,1,i)>3 && subsref(oframes,1,i)<(M-3))
+        {
+            asubsref(sel,k) = i;
+            k++;
+        }
+        m++;
+    }
+   
+    if( sel->height > 0)
+    {
+        ret = fSetArray(oframes->height, sel->height, 0);
+        {
+            for(i=0; i<sel->height; i++)
+            {
+                subsref(ret,0,i) = subsref(oframes,0,asubsref(sel,i));
+                subsref(ret,1,i) = subsref(oframes,1,asubsref(sel,i));
+                subsref(ret,2,i) = subsref(oframes,2,asubsref(sel,i));
+            }
+        }
+    }
+    else
+        ret = fSetArray(1,1,0);
+    iFreeHandle(sel);
+    return ret;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/gaussianss.c b/SD-VBS/benchmarks/sift/src/c/gaussianss.c
new file mode 100644
index 0000000..52dd95b
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/gaussianss.c
@@ -0,0 +1,185 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include "sift.h"
+F2D* resizeArray(F2D* array, int omin) 
+{
+        F2D* prev = NULL;
+        F2D* current = array;
+        int o;
+    if(omin<0)
+    {
+        for(o=1; o>=-omin; o--)
+        {
+                        prev = current;
+            current = doubleSize(current);
+                        fFreeHandle(prev);
+        }
+    }
+    if(omin>0)
+    {
+        for(o=1; o<= omin; o++)
+                {
+                        prev = current;
+            current = halveSize(current);
+                        fFreeHandle(prev);
+                }
+    }
+        return current;
+}
+/**
+    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.
+**/
+F2D** gaussianss(F2D* array, float sigman, int O, int S, int omin, int smin, int smax, float sigma0)
+{
+   /* We compute the following items in the function
+    1. Smooth input image per octave
+    2. Smooth each octave for different intervals
+    3. Subtract each "interval-1" smooth image from "interval" image per octave. So, per octave, we have "interval" * DOG images.
+    4. So, octave * intervals * image
+    5. Note: At each octave, the image is scaled down in both x and y directions
+    */
+    float k, dsigma0, dsigma;
+    int s, i, j, o, so, M, N, sbest;
+    int intervals = smax-smin+1;
+    float temp, target_sigma, prev_sigma;
+    F2D *TMP, **gss;
+    F2D* I = array;
+    // Scale multiplicative step
+    k = pow(2, (1.0/S));
+    dsigma0 = sigma0 * sqrt(1-(1.0/pow(k,2)));  // Scale step factor
+    
+    // If omin < 0, multiply the size of the image.
+        I = resizeArray(I, omin);
+    M = I->height;
+    N = I->width;
+    so = -smin+1;       // Index offset
+    gss = malloc(O*intervals*sizeof(F2D*));
+    if(gss == NULL)
+    {
+        printf("Could not allocate memory\n");
+        return NULL;
+    }
+/**
+                                                         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
+   sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 )
+**/
+    temp = sqrt(pow((sigma0*pow(k,smin)),2) - pow((sigman/pow(2,omin)),2));
+    {
+        gss[0] = fSetArray(I->height, I->width, 0);
+        imsmooth(I, temp, gss[0] );
+    }
+    for(s=smin; s<smax; s++)
+    {
+        
+    /**
+        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
+        
+            After 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 = pow(k,s+1) * dsigma0;
+        gss[s+so] = fSetArray(gss[s+so-1]->height, gss[s+so-1]->width, 0);
+        imsmooth( gss[(s+so-1)] , dsigma, gss[(s+so)] );
+    }     
+    
+    /** Other octaves **/
+    
+    for(o=1; o<O; 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-1, smax-1);
+        TMP = halveSize( gss[(o-1)*intervals+sbest+so]);
+        target_sigma = sigma0 * pow(k,smin);
+        prev_sigma = sigma0 * pow(k, (sbest+1)-S);
+        temp = sqrt(pow(target_sigma,2) - pow(prev_sigma, 2));
+        if(target_sigma > prev_sigma)
+        {
+            gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0);
+            imsmooth(TMP, temp, gss[o*intervals] );
+        }
+        else
+        {
+            int i;
+            gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0);
+            for(i=0; i<(TMP->height*TMP->width); i++)
+                asubsref(gss[o*intervals],i) = asubsref(TMP,i);
+        }
+        M = TMP->height;
+        N = TMP->width;
+        fFreeHandle(TMP);
+        for(s=smin; s<smax; s++)
+        {
+                    // The other levels are determined as above for the first octave.           
+            dsigma = pow(k,s+1) * dsigma0;
+            gss[o*intervals+s+so] = fSetArray(gss[o*intervals+s-1+so]->height, gss[o*intervals+s-1+so]->width, 0);
+            imsmooth( gss[o*intervals+s-1+so] , dsigma, gss[o*intervals+s+so]);
+        }    
+    }
+    fFreeHandle(I);
+    return gss;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/halveSize.c b/SD-VBS/benchmarks/sift/src/c/halveSize.c
new file mode 100644
index 0000000..fe1e536
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/halveSize.c
@@ -0,0 +1,35 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include "sift.h"
+F2D* halveSize(F2D* I)
+{
+    F2D *J;
+    int i, j, k;
+    int hM, hN;
+    int M, N;
+    M = I->height;
+    N = I->width;
+    
+    hM = (M+1)/2;
+    hN = (N+1)/2;
+    J = fSetArray(hM, hN, 0.0);
+    k = 0;
+    for(i=0; i<M; i+=2)
+    {
+        for(j=0; j<N; j+=2)
+        {
+            asubsref(J,k++) = subsref(I,i,j);
+        }
+    }
+    return J;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/imsmooth.c b/SD-VBS/benchmarks/sift/src/c/imsmooth.c
new file mode 100644
index 0000000..d901418
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/imsmooth.c
@@ -0,0 +1,104 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include "sift.h"
+#include <math.h>
+#include <assert.h>
+const double win_factor = 1.5 ;
+const int nbins = 36 ;
+const float threshold = 0.01;
+/**
+    This function is similar to imageBlur in common/c folder.
+    Here, we can specify the sigma value for the gaussian filter
+    function.
+**/
+void imsmooth(F2D* array, float dsigma, F2D* out)
+{
+  int M,N ;
+  int i,j,k;
+  float s ;
+  /* ------------------------------------------------------------------
+  **                                                Check the arguments
+  ** --------------------------------------------------------------- */ 
+    M = array->height;
+    N = array->width;
+    s = dsigma;
+    
+  /* ------------------------------------------------------------------
+  **                                                         Do the job
+  ** --------------------------------------------------------------- */ 
+  if(s > threshold) 
+  {
+    int W = (int) ceil(4*s) ;
+    float temp[2*W+1];
+    F2D* buffer;
+    float acc = 0.0;
+   
+    buffer = fSetArray(M,N,0);
+    
+    for(j = 0 ; j < (2*W+1) ; ++j) 
+    {
+      temp[j] = (float)(expf(-0.5 * (j - W)*(j - W)/(s*s))) ;
+      acc += temp[j];
+    }
+    for(j = 0 ; j < (2*W+1) ; ++j) 
+    {
+      temp[j] /= acc ;
+    }
+    
+    /*
+    ** Convolve along the columns
+    **/
+    for(j = 0 ; j < M ; ++j) 
+    {
+      for(i = 0 ; i < N ; ++i) 
+      {
+        int startCol = MAX(i-W,0);
+        int endCol = MIN(i+W, N-1);
+        int filterStart = MAX(0, W-i);
+                assert(j < array->height);
+                assert(j < buffer->height);
+                assert(i < buffer->width);
+        for(k=startCol; k<=endCol; k++) {
+                        assert(k < array->width);
+                        assert(filterStart < 2*W+1);
+            subsref(buffer,j,i) += subsref(array, j, k) * temp[filterStart++];
+                }
+      }
+    }
+ 
+    /*
+    ** Convolve along the rows
+    **/
+    for(j = 0 ; j < M ; ++j) 
+    {
+      for(i = 0 ; i < N ; ++i) 
+      {
+        int startRow = MAX(j-W,0);
+        int endRow = MIN(j+W, M-1);
+        int filterStart = MAX(0, W-j);
+        for(k=startRow; k<=endRow; k++)
+            subsref(out,j,i) += subsref(buffer,k,i) * temp[filterStart++];
+      }
+    }  
+  
+    fFreeHandle(buffer);
+      
+  } 
+  else 
+  {
+      for(i=0;i<M*N;i++)
+          asubsref(out, i) = asubsref(array, i);
+  }
+  return;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/script_sift.c b/SD-VBS/benchmarks/sift/src/c/script_sift.c
new file mode 100644
index 0000000..0b2f106
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/script_sift.c
@@ -0,0 +1,83 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include "sift.h"
+#include <malloc.h>
+#include "extra.h"
+#define SIFT_MEM 1<<29
+void normalizeImage(F2D* image)
+{
+        int i;
+        int rows;
+        int cols;
+ 
+    int tempMin = 10000, tempMax = -1;
+        rows = image->height;
+        cols = image->width;
+    for(i=0; i<(rows*cols); i++)
+        if(tempMin > asubsref(image,i))
+            tempMin = asubsref(image,i);
+    for(i=0; i<(rows*cols); i++)
+        asubsref(image,i) = asubsref(image,i) - tempMin;
+    for(i=0; i<(rows*cols); i++)
+        if(tempMax < asubsref(image,i))
+            tempMax = asubsref(image,i);
+    for(i=0; i<(rows*cols); i++)
+        asubsref(image,i) = ( asubsref(image,i) / (tempMax+0.0) );
+}
+int main(int argc, char* argv[])
+{
+    SET_UP
+    mallopt(M_TOP_PAD, SIFT_MEM);
+    mallopt(M_MMAP_MAX, 0);
+    I2D* im;
+    F2D *image;
+    int rows, cols, i, j;
+    F2D* frames;
+    unsigned int* startTime, *endTime, *elapsed;    
+    
+    char imSrc[100];
+    printf("Input image: ");
+    scanf("%s", imSrc);
+    im = readImage(imSrc);
+    image = fiDeepCopy(im);
+    rows = image->height;
+    cols = image->width;
+    printf("start\n");
+    for_each_job{
+        image = fiDeepCopy(im);
+        normalizeImage(image);
+        /** Extract sift features for the normalized image **/
+        frames = sift(image);
+    }
+    printf("end..\n");
+   
+#ifdef CHECK   
+    {
+        int ret=0;
+        float tol = 0.2;
+#ifdef GENERATE_OUTPUT
+        fWriteMatrix(frames, argv[1]);
+#endif
+        ret = fSelfCheck(frames, "expected_C.txt", tol);
+        if (ret == -1)
+            printf("Error in SIFT\n");
+    }
+#endif    
+    iFreeHandle(im);
+    fFreeHandle(frames);
+    WRITE_TO_FILE
+    return 0;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/sift.c b/SD-VBS/benchmarks/sift/src/c/sift.c
new file mode 100644
index 0000000..ea0d308
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/sift.c
@@ -0,0 +1,314 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#include <math.h>
+#include "sift.h"
+/** SIFT- Scale Invariant Feature Transform. This algorithm is based on
+    David Lowe's implementation of sift. So, we will use the parameter
+    values from Lowe's implementation.
+    See: http://www.cs.ubc.ca/~lowe/keypoints/
+    SIFT extracts from an image a collection of frames or keypoints. These
+    are oriented disks attacked to blob-like structures of the image. As
+    the image translates, rotates and scales, the frames track these blobs 
+    and the deformation.
+    'BoundaryPoint' [bool]
+    If set to 1, frames too close to the image boundaries are discarded.
+    
+    'Sigma0' [pixels]
+    Smoothing of the level 0 of octave 0 of the scale space.
+    (Note that Lowe's 1.6 value refers to the level -1 of octave 0.)
+    'SigmaN' [pixels]
+    Nominal smoothing of the input image. Typically set to 0.5.
+    'Threshold'
+    Threshold used to eliminate weak keypoints. Typical values for
+    intensity images in the range [0,1] are around 0.01. Smaller
+    values mean more keypoints.
+**/
+F2D* sift(F2D* I)
+{
+    int rows, cols, K;
+    int subLevels, omin, Octaves, r, NBP, NBO, smin, smax, intervals, o;
+    float sigma, sigman, sigma0, thresh, magnif;
+    int discardBoundaryPoints;
+    F2D *descriptors;
+    F2D **gss, *tfr;
+    F2D **dogss;
+    I2D* i_, *j_, *s_;
+    F2D *oframes, *frames;
+    int i, j, k, m, startRow, startCol, endRow, endCol, firstIn=1;
+    float minVal;
+    I2D* tx1, *ty1, *ts1;
+    I2D* x1, *y1, *s1, *txys;
+    rows = I->height;
+    cols = I->width;
+    /** 
+    Lowe's choices
+    Octaves - octave
+    subLevels - sub-level for image
+    sigma - sigma value for gaussian kernel, for smoothing the image
+    At each successive octave, the data is spatially downsampled by half
+    **/
+    subLevels = 3;
+    omin = -1;
+    minVal = log2f(MIN(rows,cols));
+    Octaves = (int)(floor(minVal))-omin-4;   /* Upto 16x16 images */
+    sigma0 = 1.6 * pow(2, (1.0/subLevels));
+    sigman = 0.5;
+    thresh = (0.04/subLevels)/2;
+    r = 10;
+#ifdef test
+    subLevels = 1;
+    Octaves = 1;
+    sigma0 = pow(0.6 * 2, 1);
+    sigman = 1.0;
+    thresh = (1/subLevels)/2;
+    r = 1;
+#endif
+    NBP = 4;
+    NBO = 8 ;
+    magnif = 3.0 ;
+    discardBoundaryPoints = 1 ;
+    smin = -1;
+    smax = subLevels+1;
+    intervals = smax - smin + 1;
+    /** 
+        We build gaussian pyramid for the input image. Given image I,
+        we sub-sample the image into octave 'Octaves' number of levels. At 
+        each level, we smooth the image with varying sigman values.
+        
+        Gaussiam pyramid can be assumed as a 2-D matrix, where each
+        element is an image. Number of rows corresponds to the number
+        of scales of the pyramid (octaves, "Octaves"). Row 0 (scale 0) is the 
+        size of the actual image, Row 1 (scale 1) is half the actual 
+        size and so on.
+        At each scale, the image is smoothened with different sigma values.
+        So, each row has "intervals" number of smoothened images, starting
+        with least blurred.
+        
+        gss holds the entire gaussian pyramid. 
+    **/
+    gss = gaussianss(I, sigman, Octaves, subLevels, omin, -1, subLevels+1, sigma0);
+    /** 
+        Once we build the gaussian pyramid, we compute DOG, the 
+        Difference of Gaussians. At every scale, we do:
+        dogss[fixedScale][0] = gss[fixedScale][1] - gss[fixedScale][0]
+        Difference of gaussian gives us edge information, at each scale.
+        In order to detect keypoints on the intervals per octave, we 
+        inspect DOG images at highest and lowest scales of the octave, for
+        extrema detection.
+    **/
+    
+    dogss = diffss(gss, Octaves, intervals);
+    /** The extraction of keypoints is carried one octave per time **/
+    for(o=0; o<Octaves; o++)
+    {
+        F2D *temp;
+        F2D *t;
+        F2D *negate;
+        F2D *idxf, *idxft;
+        I2D *idx;
+        int i1, j1;
+        I2D *s, *i, *j, *x, *y;
+        int sizeRows = dogss[o*intervals]->height;
+        int sizeCols = dogss[o*intervals]->width;
+        
+        {
+            int i,j,k=0;
+            temp = fMallocHandle(intervals-1, sizeRows*sizeCols);
+            negate = fMallocHandle(intervals-1, sizeRows*sizeCols);
+            /** 
+                Keypoints are detected as points of local extrema of the
+                DOG pyramid, at a given octave. In softlocalmax function, the
+                keypoints are extracted by looking at 9x9x9 neighborhood samples.
+                
+                We populate temp and negate arrays with the values of the DOG
+                pyramid for a given octave, o. Since we are interested in both
+                local maxima and minima, we compute negate matrix, which is the
+                negated values of the DOG pyramid.
+            
+            **/
+            for(i1=0; i1<(intervals-1); i1++)
+            {
+                for(j=0; j<sizeCols; j++)
+                {
+                    for(i=0; i<sizeRows; i++)
+                    {
+                        asubsref(temp,k) = subsref(dogss[o*intervals+i1],i,j);
+                        asubsref(negate,k++) = -subsref(dogss[o*intervals+i1],i,j);
+                    }
+                }
+            }
+        }
+        
+        /**
+            siftlocalmax returns indices k, that correspond to local maxima and 
+            minima. 
+                The 80% tricks discards early very weak points before refinement.
+        **/
+        idxf = siftlocalmax( temp, 0.8*thresh, intervals, sizeRows, sizeCols);
+        t = siftlocalmax( negate, 0.8*thresh, intervals, sizeRows, sizeCols);
+        
+        idxft = fHorzcat(idxf, t);
+        /**
+            Since indices is the 1-D index of the temp/negate arrays, we compute 
+            the x,y and intervals(s) co-ordinates corresponding to each index.
+        **/
+        
+        idx = ifDeepCopy(idxft);
+        
+        x = iSetArray(idx->height,idx->width,0);
+        y = iSetArray(idx->height,idx->width,0);
+        s_ = iSetArray(idx->height,idx->width,0);
+        
+        {
+            int i, j;
+            for(i=0; i<idx->height; i++)
+             {
+                 for(j=0; j<idx->width; j++)
+                 {
+                    int v, u, w, z;
+                    w = subsref(idx,i,j);
+                    v = ceil((w/(sizeRows*sizeCols)) + 0.5);
+                    u = floor(w/(sizeRows*sizeCols));
+                    z = w - (sizeRows*sizeCols*u);
+                    
+                    /** v is the interval number, s **/
+                    subsref(s_,i,j) = v;
+                    /** row number of the index **/
+                    subsref(y,i,j) = ceil((z / sizeRows)+0.5);
+                    /** col number of the index **/
+                    subsref(x,i,j) = z - (sizeCols * floor(z / sizeRows));
+                 }
+             }
+        }
+        {
+            
+            tx1 = isMinus(x, 1);
+            x1 = iReshape(tx1, 1, (tx1->height*tx1->width));
+            ty1 = isMinus(y, 1);
+            y1 = iReshape(ty1, 1, (ty1->height*ty1->width));
+            ts1 = isPlus(s_, (smin-1));
+            s1 = iReshape(ts1, 1, (ts1->height*ts1->width));
+        
+            txys = iVertcat(y1, x1);
+            i_ = iVertcat(txys, s1);
+        }
+        /**
+            Stack all x,y,s into oframes.
+            Row 0 of oframes = x
+            Row 1 of oframes = y
+            Row 2 of oframes = s
+        **/
+        oframes = fiDeepCopy(i_);
+         
+        {
+            int i,j;
+            F2D* temp;
+            temp = oframes;
+                
+            /** 
+                Remove points too close to the boundary
+            **/
+            
+            if(discardBoundaryPoints)
+                oframes = filterBoundaryPoints(sizeRows, sizeCols, temp);
+            fFreeHandle(temp);
+        }
+        
+        /**
+            Refine the location, threshold strength and remove points on edges
+        **/
+        if( asubsref(oframes,0) != 0)
+        {
+            F2D* temp_;
+            temp_ = fTranspose(oframes);
+            fFreeHandle(oframes);
+            oframes = siftrefinemx(temp_, temp, smin, thresh, r, sizeRows, sizeCols, intervals-1);
+            fFreeHandle(temp_);
+            if( firstIn == 0)
+            {
+                tfr = fDeepCopy(frames);
+                fFreeHandle(frames);
+                frames = fHorzcat(tfr, oframes);
+                fFreeHandle(tfr);
+            }
+            else
+                frames = fDeepCopy(oframes);
+            firstIn = 0;
+        }
+        else if(Octaves == 1)
+            frames = fDeepCopy(oframes);
+        
+        fFreeHandle(oframes);
+        iFreeHandle(y);
+        iFreeHandle(x);
+        iFreeHandle(s_);
+        iFreeHandle(y1);
+        iFreeHandle(x1);
+        iFreeHandle(s1);
+        iFreeHandle(ty1);
+        iFreeHandle(tx1);
+        iFreeHandle(ts1);
+        iFreeHandle(txys);
+        iFreeHandle(i_);
+        iFreeHandle(idx);
+        fFreeHandle(idxf);
+        fFreeHandle(idxft);
+        fFreeHandle(temp);
+        fFreeHandle(t);
+        fFreeHandle(negate);
+    }
+    { int s;
+    for(o=0; o<Octaves; o++)
+    {
+        for(s=0; s<(intervals-1); s++)
+        {
+            fFreeHandle(dogss[o*intervals+s]);
+        }
+    }
+    for(o=0; o<Octaves; o++)
+    {
+        for(s=0; s<(intervals); s++)
+        {
+            fFreeHandle(gss[o*intervals+s]);
+        }
+    }
+    }
+    free(gss);
+    free(dogss);
+   
+    return frames;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/sift.h b/SD-VBS/benchmarks/sift/src/c/sift.h
new file mode 100644
index 0000000..e8a4ee0
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/sift.h
@@ -0,0 +1,27 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+#ifndef _SIFT_
+#define _SIFT_
+#include "sdvbs_common.h"
+#include <assert.h>
+#define GREATER(a,b) ((a) > (b))
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#define MIN(a,b) (((a)<(b))?(a):(b))
+F2D* sift(F2D* I);
+F2D* halveSize(F2D* I);
+F2D** gaussianss(F2D* I, float sigman, int O, int S, int omin, int smin, int smax, float sigma0);
+F2D** diffss(F2D** ss, int O, int intervals);
+F2D* doubleSize(F2D* I);
+void imsmooth(F2D* I_pt, float dsigma, F2D* out);
+F2D* siftlocalmax(F2D* in, float thresh, int intervals, int M, int N);
+F2D* filterBoundaryPoints(int M, int N, F2D* oframes); 
+F2D* siftrefinemx(F2D* oframes, F2D* dogss, int smin, float thresh, int rin, int M, int N, int intervals);
+#endif
diff --git a/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c b/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c
new file mode 100644
index 0000000..d75bcbd
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c
@@ -0,0 +1,247 @@
+#include "sift.h"
+/**
+    SIFTLOCALMAX  Find local maximizers
+   Returns the indexes of the local maximizers of
+   the 3-dimensional array F.
+    Say, we have a 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!).
+**/
+#define NDIMS 3
+F2D* siftlocalmax(F2D* in, float thresh, int intervals, int M, int N)
+{  
+  int ndims, ptoffset=0, maxIter = 0; 
+  int offsets[NDIMS];
+  int midx[NDIMS];
+  int dims[NDIMS];
+  I2D* neighbors ;
+  int nneighbors ;
+  F2D* out;
+    
+ /**  
+    We pass dogss[o], which has >1 intervals
+    If >1 intervals, then the dimensions of in[F] will be 1 x intervals
+    If =1 interval, then the dimensions of in[F] will be the size of the dogss[o] image
+    We will check for this condition.
+ **/
+    
+    ndims = NDIMS;      /* Since we have intervals number of images of size M x N */
+        dims[0] = M;
+        dims[1] = N;
+        dims[2] = intervals-1;
+  /* 
+     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).
+  */
+  
+  
+ 
+  /* ------------------------------------------------------------------
+   *                                                         Do the job
+   * --------------------------------------------------------------- */  
+    int maxima_size = M*N ;
+    I2D* maxima_start = iMallocHandle(1, maxima_size);
+    int* maxima_iterator = maxima_start->data ;
+    int* maxima_end = maxima_start->data + maxima_size ;
+    int i,h,o,j; 
+    F2D* pt;  
+    
+    pt = in;
+    /* Compute the offsets between dimensions. */
+    offsets[0] = 1 ;
+    for(i = 1 ; i < ndims ; ++i)
+    {
+      offsets[i] = offsets[i-1]*dims[i-1] ;
+    }
+    /* Multi-index. */
+    for(i = 0 ; i < ndims ; ++i)
+      midx[i] = 1 ;
+    /* Neighbors. */
+    nneighbors = 1 ;
+    o=0 ;
+    for(i = 0 ; i < ndims ; ++i) 
+    {
+      nneighbors *= 3 ;
+      midx[i] = -1 ;
+      o -= offsets[i] ;
+    }
+    nneighbors -= 1 ;
+    neighbors = iMallocHandle(1,nneighbors);
+    /* Precompute offsets from offset(-1,...,-1,0,...0) to
+     * offset(+1,...,+1,0,...,0). */
+    i = 0 ;
+    while(1) 
+    {
+      if(o != 0)
+      {
+        asubsref(neighbors, i) = o ;
+        i++;
+      }
+      h = 0 ;
+      while( o += offsets[h], (++midx[h]) > 1 ) 
+      {
+        o -= 3*offsets[h] ;
+        midx[h] = -1 ;
+        if(++h >= ndims)
+          goto stop ;
+      }
+    }
+  stop: ;
+    /* Starts at the corner (1,1,...,1,0,0,...0) */
+    for(h = 0 ; h < ndims ; ++h) 
+    {
+      midx[h] = 1 ;
+      ptoffset += offsets[h] ;
+    }
+    /* ---------------------------------------------------------------
+     *                                                            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 < ndims ; ++h)
+      if(dims[h] < 3) 
+                goto end ;
+   
+    while(1) 
+    {
+      /* Propagate carry along multi index midx */
+      
+      h = 0 ;
+      while(midx[h] >= dims[h] - 1) 
+      {
+        midx[h] = 1 ;
+        if(++h >= ndims) 
+          goto next_layer ;
+        ++midx[h] ;
+      }
+      
+      /*  Scan neighbors */
+      {
+        float v = asubsref(pt, ptoffset);  //*pt ;
+        int is_greater = (v>=thresh) ? 1 : 0;
+        
+                assert(ptoffset < pt->width*pt->height);
+        i = 0  ;
+        while(is_greater && i < nneighbors)
+        {
+            if(v > asubsref(pt, ptoffset + asubsref(neighbors,i)))
+                is_greater *= 1;
+            else
+                is_greater *= 0;
+            i = i+1;
+        }
+        
+        /* Add the local maximum */
+        if(is_greater) 
+        {
+          /* Need more space? */
+          if( &(maxima_iterator[maxIter])  == maxima_end) 
+          {
+              int *temp, i, j;
+              maxima_size += M*N ;
+            
+            free(maxima_start); 
+            maxima_start = iMallocHandle(1, maxima_size);
+            maxima_end = maxima_start->data + maxima_size ;
+            maxima_iterator = maxima_end - M*N ;
+            maxIter = 0;
+          }
+          
+          maxima_iterator[maxIter] = ptoffset + 1;
+          maxIter++;
+        }
+        
+        /* Go to next element */
+        ptoffset += 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. */
+    {
+        int i=0;
+        out = fMallocHandle(1, maxIter);
+      
+        for(i=0; i<maxIter; i++)
+            asubsref(out,i) = maxima_iterator[i] ;
+    }
+    
+    /* Release space. */
+    free(neighbors) ;
+    free(maxima_start);
+  
+  return out;
+}
diff --git a/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c b/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c
new file mode 100644
index 0000000..2b10e9d
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c
@@ -0,0 +1,196 @@
+#include "sift.h"
+const int max_iter = 5 ;
+/**
+    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.
+**/
+F2D* siftrefinemx(F2D* oframes, F2D* dogss, int smin, float thresh, int r, int M, int N, int intervals)
+{
+  int S,K ;
+  F2D* out;
+        
+  /* -----------------------------------------------------------------
+  **                                               Check the arguments
+  ** -------------------------------------------------------------- */ 
+ 
+  S = intervals;
+  K = oframes->height;
+  /* If the input array is empty, then output an empty array as well. */
+  if( K == 0) {
+      out = fDeepCopy(oframes);
+    return out;
+  }
+  /* -----------------------------------------------------------------
+   *                                                        Do the job
+   * -------------------------------------------------------------- */
+  {    
+    F2D* buffer = fMallocHandle(K,3);
+    int p ;
+    const int yo = 1 ;
+    const int xo = M ;
+    const int so = M*N;
+    int buffCounter = 0;
+    int pptcount = 0;
+    
+    for(p = 0 ; p < K ; ++p) 
+    {
+        float tx, ty, ts;
+        int x,y,s;
+        int iter ;
+        float b[3] ;
+        tx = asubsref(oframes, pptcount++);
+        ty = asubsref(oframes, pptcount++);
+        ts = asubsref(oframes, pptcount++);
+        x = ceil(tx);
+        y = ceil(ty);
+        s = ceil(ts) - smin; 
+    
+      /* 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.
+       */
+      
+      if(x < 1 || x > N-2 ||
+         y < 1 || y > M-2 ||
+         s < 1 || s > S-2) {
+        continue ;
+      }
+#define at(dx,dy,ds) asubsref(dogss, (dx)*xo + (dy)*yo + (ds)*so)
+      {
+        float Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
+        int dx = 0 ;
+        int dy = 0 ;
+        
+        for(iter = 0 ; iter < max_iter ; ++iter) 
+        {
+          float A[3*3] ; 
+#define Aat(i,j) (A[(i)+(j)*3])    
+          x += dx ;
+          y += dy ;
+          
+          /* Compute the gradient. */
+          Dx = 0.5 * (at(x+1,y+0,s+0) - at(x-1,y+0,s+0)) ;
+          Dy = 0.5 * (at(x+0,y+1,s+0) - at(x+0,y-1,s+0));
+          Ds = 0.5 * (at(x+0,y+0,s+1) - at(x+0,y+0,s-1)) ;
+          
+          /* Compute the Hessian. */
+          Dxx = (at(x+1,y,s) + at(x-1,y,s) - 2.0 * at(x,y,s)) ;
+          Dyy = (at(x,y+1,s) + at(x,y-1,s) - 2.0 * at(x,y,s)) ;
+          Dss = (at(x,y,s+1) + at(x,y,s-1) - 2.0 * at(x,y,s)) ;
+          
+          Dxy = 0.25 * ( at(x+1,y+1,s) + at(x-1,y-1,s) - at(x-1,y+1,s) - at(x+1,y-1,s) ) ;
+          Dxs = 0.25 * ( at(x+1,y,s+1) + at(x-1,y,s-1) - at(x-1,y,s+1) - at(x+1,y,s-1) ) ;
+          Dys = 0.25 * ( at(x,y+1,s+1) + at(x,y-1,s-1) - at(x,y-1,s+1) - at(x,y+1,s-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 ;
+          
+          /* 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 ;
+          
+        }
+                                
+        {
+          float val = at(x,y,s) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ; 
+          float score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ; 
+          float xn = x + b[0] ;
+          float yn = y + b[1] ;
+          float sn = s + b[2] ;
+        
+          if(fabs(val) > thresh &&
+             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) 
+                  {
+            asubsref(buffer,buffCounter++) = xn ;
+            asubsref(buffer,buffCounter++) = yn ;
+            asubsref(buffer,buffCounter++) = sn+smin ;
+          }
+        }
+      }
+    }      
+    /* Copy the result into an array. */
+    {
+      int i, j, k=0;
+      int NL = buffCounter/3;
+      out = fMallocHandle(3, NL);
+        for(i=0; i<NL; i++)
+        {
+            for(j=0; j<3; j++)
+            {
+                subsref(out,j,i) = asubsref(buffer,k);
+                k++;
+            }
+        }
+    }
+    free(buffer) ;
+  }
+    
+    return out;
+}
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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
Binary files differ
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