Squashed commit of the sb-vbs branch.

Includes the SD-VBS benchmarks modified to:
- Use libextra to loop as realtime jobs
- Preallocate memory before starting their main computation
- Accept input via stdin instead of via argc

Does not include the SD-VBS matlab code.

Fixes libextra execution in LITMUS^RT.

11 files changed, 1360 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;
+}
+
+
+
+
+