8 files changed, 887 insertions, 0 deletions

diff --git a/SD-VBS/benchmarks/tracking/src/c/calcAreaSum.c b/SD-VBS/benchmarks/tracking/src/c/calcAreaSum.c
new file mode 100644
index 0000000..40f03dc
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/calcAreaSum.c
@@ -0,0 +1,81 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** Compute the sum of pixels over pixel neighborhood.
+    Neighborhood = winSize*winSize 
+    This will be useful when we compute the displacement
+    of the neighborhood across frames instead of tracking each pixel.
+
+    Input:  src Image
+            # rows, cols, size of window
+    Output: windowed-sum image, ret.
+
+    Example:
+    
+    winSize = 4, cols = 8, rows = 16
+    nave_half = 2, nave = 4
+    Say, the first row of the image is:
+                3 8 6 2 4 8 9 5
+    a1 =    0 0 3 8 6 2 4 8 9 5 0 0 (append nave_half zeros to left and right border)
+    asum =  (sum the first nave # pixels - 0 0 3 8 ) = 11
+    ret(0,0) =  11
+    For ret(0,1), we need to move the window to the right by one pixel and subtract
+    from a1sum the leftmost pixel. So, we add value 6 and subtract value at a1(0,0), = 0 here.
+    ret(0,1) = 17 = a1sum
+    For ret(0,2), a1sum - a1(0,1) + a1(2+nave) = 17 - 0 + 2 = 19 = a1sum
+    For ret(0,3), a1sum - a1(0,2) + a1(3+nave) = 19 - 3 + 4 = 20 = a1sum
+
+    We proceed this way for all the rows and then perform summantion across all cols.
+**/
+F2D* calcAreaSum(F2D* src, int cols, int rows, int winSize)
+{
+    int nave, nave_half, i, j, k;
+    F2D *ret, *a1;
+    float a1sum;
+    
+    nave = winSize;
+    nave_half = floor((nave+1)/2);
+
+    ret = fMallocHandle(rows, cols);
+    a1 = fSetArray(1, cols+nave,0);
+
+    for(i=0; i<rows; i++)
+    {
+        for(j=0; j<cols; j++) {
+            asubsref(a1,j+nave_half) = subsref(src,i,j);
+                }
+        a1sum = 0;
+        for(k=0; k<nave; k++) {
+            a1sum += asubsref(a1,k);
+                }
+        for(j=0; j<cols; j++)
+        {
+            subsref(ret,i,j) = a1sum;
+            a1sum += asubsref(a1,j+nave) - asubsref(a1,j);
+        }
+    }
+    fFreeHandle(a1);
+
+    a1 = fSetArray(1, rows+nave,0);
+    for(i=0; i<cols; i++)
+    {
+        for(j=0; j<rows; j++) {
+            asubsref(a1,j+nave_half) = subsref(ret,j,i);
+                }
+        a1sum = 0;
+        for(k=0; k<nave; k++) {
+            a1sum += asubsref(a1,k);
+                }
+        for(j=0; j<rows; j++)
+        {
+            subsref(ret,j,i) = a1sum;
+            a1sum += asubsref(a1,j+nave) - asubsref(a1,j);
+        }
+    }
+    fFreeHandle(a1);
+
+    return ret; 
+}
diff --git a/SD-VBS/benchmarks/tracking/src/c/calcGoodFeature.c b/SD-VBS/benchmarks/tracking/src/c/calcGoodFeature.c
new file mode 100644
index 0000000..0ef1862
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/calcGoodFeature.c
@@ -0,0 +1,78 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** Computes lambda matrix, strength at each pixel
+    
+    det = determinant( [ IverticalEdgeSq IhorzVertEdge; IhorzVertEdge IhorizontalEdgeSq] ) ;
+    tr = IverticalEdgeSq + IhorizontalEdgeSq;
+    lamdba = det/tr;
+
+    Lambda is the measure of the strength of pixel
+    neighborhood. By strength we mean the amount of 
+    edge information it has, which translates to
+    sharp features in the image.
+
+    Input:  Edge images - vertical and horizontal
+            Window size (neighborhood size)
+    Output: Lambda, strength of pixel neighborhood
+
+    Given the edge images, we compute strength based
+    on how strong the edges are within each neighborhood.
+
+**/
+
+F2D* calcGoodFeature(F2D* verticalEdgeImage, F2D* horizontalEdgeImage, int cols, int rows, int winSize)
+{
+    int i, j, k, ind;
+    F2D *verticalEdgeSq, *horizontalEdgeSq, *horzVertEdge;
+    F2D *tr, *det, *lambda;
+    F2D *cummulative_verticalEdgeSq, *cummulative_horzVertEdge, *cummulative_horizontalEdgeSq;
+
+    verticalEdgeSq = fMallocHandle(rows, cols);
+    horzVertEdge = fMallocHandle(rows, cols);
+    horizontalEdgeSq = fMallocHandle(rows, cols);
+    
+    for( i=0; i<rows; i++)
+    {
+        for( j=0; j<cols; j++)
+        {
+            subsref(verticalEdgeSq,i,j) = subsref(verticalEdgeImage,i,j) * subsref(verticalEdgeImage,i,j);
+            subsref(horzVertEdge,i,j) = subsref(verticalEdgeImage,i,j) * subsref(horizontalEdgeImage,i,j);
+            subsref(horizontalEdgeSq,i,j) = subsref(horizontalEdgeImage,i,j) * subsref(horizontalEdgeImage,i,j);
+        }
+    }
+    
+    cummulative_verticalEdgeSq = calcAreaSum(verticalEdgeSq, cols, rows, winSize); 
+    cummulative_horzVertEdge = calcAreaSum(horzVertEdge, cols, rows, winSize);
+    cummulative_horizontalEdgeSq = calcAreaSum(horizontalEdgeSq, cols, rows, winSize);    
+     
+    tr = fMallocHandle(rows, cols);
+    det = fMallocHandle(rows, cols);
+    lambda = fMallocHandle(rows, cols);
+    
+    for( i=0; i<rows; i++)
+    {
+        for( j=0; j<cols; j++)
+        {
+            subsref(tr,i,j) = subsref(cummulative_verticalEdgeSq,i,j) + subsref(cummulative_horizontalEdgeSq,i,j);
+            subsref(det,i,j) = subsref(cummulative_verticalEdgeSq,i,j) * subsref(cummulative_horizontalEdgeSq,i,j) - subsref(cummulative_horzVertEdge,i,j) * subsref(cummulative_horzVertEdge,i,j);
+            subsref(lambda,i,j) = ( subsref(det,i,j) / (subsref(tr,i,j)+0.00001) ) ;
+        }
+    }
+   
+    fFreeHandle(verticalEdgeSq);
+    fFreeHandle(horzVertEdge);
+    fFreeHandle(horizontalEdgeSq);
+
+    fFreeHandle(cummulative_verticalEdgeSq);
+    fFreeHandle(cummulative_horzVertEdge);
+    fFreeHandle(cummulative_horizontalEdgeSq);
+    
+    fFreeHandle(tr);
+    fFreeHandle(det);
+     
+    return lambda;
+}
diff --git a/SD-VBS/benchmarks/tracking/src/c/calcPyrLKTrack.c b/SD-VBS/benchmarks/tracking/src/c/calcPyrLKTrack.c
new file mode 100644
index 0000000..3718673
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/calcPyrLKTrack.c
@@ -0,0 +1,179 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** calcPyrLKTrack tries to find the displacement (translation in x and y) of features computed in the previous frame.
+    To compute the displacement, we interpolate the existing feature position around the pixel neighborhood.
+    For better results, we perform interpolatations across all pyramid levels.
+
+    Input:  Previous frame blurred images, vertical and horizontal
+            Current frame edge images, vertical and horizontal
+            Current frame blurred images, for level 0 and 1
+            Number of features from previous frame
+            Window size
+            Accuracy
+            Number of iterations to compute motion vector
+            Features from the previous frame
+            currentFrameFeatures, Populate the current frame features' displacements
+    Output: Number of valid features
+
+**/
+
+I2D* calcPyrLKTrack(F2D* previousImageBlur_level1, F2D* previousImageBlur_level2, F2D* vertEdge_level1, F2D* vertEdge_level2, F2D* horzEdge_level1, F2D* horzEdge_level2, F2D* currentImageBlur_level1, F2D* currentImageBlur_level2, F2D* previousFrameFeatures, int nFeatures, int winSize, float accuracy, int max_iter, F2D* currentFrameFeatures)
+{
+    int idx, level, pLevel, i, j, k, winSizeSq;
+    I2D *valid;
+    F2D *rate, *iPatch, *jPatch, *iDxPatch;
+    F2D *iDyPatch;
+    float tr, x, y, dX, dY, c_xx, c_yy, c_xy;
+    int imgSize_1, imgSize_2;
+    float mX, mY, dIt, eX, eY, c_det;
+    I2D *imgDims;
+
+    imgDims = iMallocHandle(2, 2);
+    subsref(imgDims,0,0) = previousImageBlur_level1->height;
+    subsref(imgDims,0,1) = previousImageBlur_level1->width;
+    subsref(imgDims,1,0) = previousImageBlur_level2->height;
+    subsref(imgDims,1,1) = previousImageBlur_level2->width;
+
+    pLevel = 2;
+    rate = fMallocHandle(1, 6);
+
+    asubsref(rate,0) = 1;
+    asubsref(rate,1) = 0.5;
+    asubsref(rate,2) = 0.25;
+    asubsref(rate,3) = 0.125;
+    asubsref(rate,4) = 0.0625;
+    asubsref(rate,5) = 0.03125;
+    
+    winSizeSq = 4*winSize*winSize;
+    valid = iSetArray(1,nFeatures, 1);
+
+    /** For each feature passed from previous frame, compute the dx and dy, the displacements **/
+    for(i=0; i<nFeatures; i++)
+    {
+        dX = 0;
+        dY = 0;
+
+        /** Compute the x and y co-ordinate values at "pLevel" **/
+        x = subsref(previousFrameFeatures,0,i) * asubsref(rate,pLevel);
+        y = subsref(previousFrameFeatures,1,i) * asubsref(rate,pLevel);
+        c_det = 0;
+
+        /** For each pyramid level, try to find correspondence.
+            We look for the correspondence in a given window size
+            , (winSize x winSize) neighborhood **/
+        for(level = pLevel-1; level>=0; level--)
+        {
+            x = x+x;
+            y = y+y;
+            dX = dX + dX;
+            dY = dY + dY;
+            imgSize_1 = subsref(imgDims,level,0);
+            imgSize_2 = subsref(imgDims,level,1);
+
+            c_xx = 0;
+            c_xy = 0;
+            c_yy = 0;
+        
+            if((x-winSize)<0 || (y-winSize)<0 || (y+winSize+1)>=imgSize_1 || (x+winSize+1)>=imgSize_2)
+            {
+                asubsref(valid,i) = 0;
+                break;
+            }
+
+            /** Perform interpolation. Use co-ord from previous
+                frame and use the images from current frame **/
+
+            if(level ==0)
+            {
+                iPatch = getInterpolatePatch(previousImageBlur_level1, imgSize_2, x, y, winSize);
+                iDxPatch = getInterpolatePatch(vertEdge_level1, imgSize_2, x, y, winSize);
+                iDyPatch = getInterpolatePatch(horzEdge_level1, imgSize_2, x, y, winSize);
+            }
+            if(level ==1)
+            {
+                iPatch = getInterpolatePatch(previousImageBlur_level2, imgSize_2, x, y, winSize);
+                iDxPatch = getInterpolatePatch(vertEdge_level2, imgSize_2, x, y, winSize);
+                iDyPatch = getInterpolatePatch(horzEdge_level2, imgSize_2, x, y, winSize);
+            }
+
+            /** Compute feature strength in similar way as calcGoodFeature **/
+            for(idx=0; idx<winSizeSq; idx++)
+            {
+                c_xx += asubsref(iDxPatch,idx) * asubsref(iDxPatch,idx);
+                c_xy += asubsref(iDxPatch,idx) * asubsref(iDyPatch,idx);
+                c_yy += asubsref(iDyPatch,idx) * asubsref(iDyPatch,idx);
+            }
+            
+            c_det = (c_xx * c_yy -c_xy * c_xy);
+            tr = c_xx + c_yy;
+
+            if((c_det/(tr+0.00001)) < accuracy)
+            {
+                asubsref(valid,i) = 0;
+                fFreeHandle(iPatch);
+                fFreeHandle(iDxPatch);
+                fFreeHandle(iDyPatch);
+                break;
+            }
+            c_det = 1/c_det;
+
+            /** We compute dX and dY using previous frame and current
+                frame images. For this, the strength is computed at each
+                pixel in the new image.  **/
+            for(k=0; k<max_iter; k++)
+            {
+                if( (x+dX-winSize)<0 || (y+dY-winSize)<0 || (y+dY+winSize+1)>=imgSize_1 || (x+dX+winSize+1)>=imgSize_2)
+                {
+                    asubsref(valid,i) = 0;
+                    break;
+                }
+                
+                if(level == 0)
+                    jPatch = getInterpolatePatch(currentImageBlur_level1, imgSize_2, x+dX, y+dY, winSize);
+                if(level == 1)
+                    jPatch = getInterpolatePatch(currentImageBlur_level2, imgSize_2, x+dX, y+dY, winSize);
+                
+                eX = 0;
+                eY = 0;
+                for(idx=0; idx<winSizeSq; idx++)
+                {
+                    dIt = asubsref(iPatch,idx) - asubsref(jPatch,idx);
+                    eX += dIt * asubsref(iDxPatch,idx);
+                    eY += dIt * asubsref(iDyPatch,idx);
+                }
+                
+                mX = c_det * (eX * c_yy - eY * c_xy);
+                mY = c_det * (-eX * c_xy + eY * c_xx);
+                dX = dX + mX;
+                dY = dY + mY;
+            
+                if( (mX*mX+mY*mY) < accuracy)
+                {
+                    fFreeHandle(jPatch);
+                    break;
+                }
+                
+                fFreeHandle(jPatch);
+            }
+
+            fFreeHandle(iPatch);
+            fFreeHandle(iDxPatch);
+            fFreeHandle(iDyPatch);
+        }
+
+        subsref(currentFrameFeatures,0,i) = subsref(previousFrameFeatures,0,i) + dX;
+        subsref(currentFrameFeatures,1,i) = subsref(previousFrameFeatures,1,i) + dY;
+        
+    }
+
+    fFreeHandle(rate);
+    iFreeHandle(imgDims);
+    return valid;
+}
+
+
+
diff --git a/SD-VBS/benchmarks/tracking/src/c/fillFeatures.c b/SD-VBS/benchmarks/tracking/src/c/fillFeatures.c
new file mode 100644
index 0000000..42de169
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/fillFeatures.c
@@ -0,0 +1,67 @@
+/********************************
+
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** Find the position and values of the top N_FEA features 
+        from the lambda matrix **/
+
+F2D* fillFeatures(F2D* lambda, int N_FEA, int win)
+{
+    int i, j, k, l;
+    int rows = lambda->height;
+    int cols = lambda->width;
+    F2D* features;
+
+        features = fSetArray(3, N_FEA, 0);
+
+        /** init array **/
+    for(i=0; i<N_FEA; i++)
+    {
+        subsref(features, 0, i) = -1.0;
+        subsref(features, 1, i) = -1.0;
+        subsref(features, 2, i) = 0.0;
+    }
+    
+
+        /** 
+                Find top N_FEA values and store them in 
+        features array along with row and col information 
+                It should be possible to make this algorithm better
+                if we use a pointer-based data structure,
+                but have not implemented due to MATLAB compatibility
+        **/ 
+
+        for (i=win; i<rows-win; i++)
+        {
+                for (j=win; j<cols-win; j++)
+                {
+            float currLambdaVal = subsref(lambda,i,j);
+                        if (subsref(features, 2, N_FEA-1) > currLambdaVal)
+                                continue;
+                        
+                        for (k=0; k<N_FEA; k++) 
+                        {
+                                if (subsref(features, 2, k) < currLambdaVal) 
+                                {
+                                        /** shift one slot **/
+                                        for (l=N_FEA-1; l>k; l--)
+                                        {
+                                                subsref(features, 0, l) = subsref(features, 0, l-1);
+                                                subsref(features, 1, l) = subsref(features, 1, l-1);
+                                                subsref(features, 2, l) = subsref(features, 2, l-1);
+                                        }       
+                                        
+                                subsref(features, 0, k) = j * 1.0;
+                                subsref(features, 1, k) = i * 1.0;
+                                subsref(features, 2, k) = currLambdaVal;
+                                        break;
+                                }
+                        }
+                }
+        }
+
+    return features;
+}
diff --git a/SD-VBS/benchmarks/tracking/src/c/getANMS.c b/SD-VBS/benchmarks/tracking/src/c/getANMS.c
new file mode 100644
index 0000000..7ecce8b
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/getANMS.c
@@ -0,0 +1,156 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** ANMS - Adaptive Non-Maximal Suppression
+    This function takes features as input and suppresses those which
+    are close to each other (within the SUPPRESSION_RADIUS) and have
+    similar feature strength
+**/
+F2D *getANMS (F2D *points, float r)
+{
+    float MAX_LIMIT = 100000;
+    F2D *suppressR;
+    float C_ROBUST = 1.0;
+    F2D *srtdPnts;
+    int n;
+    I2D *srtdVIdx, *supId;
+    float t, t1, r_sq;
+    F2D *tempF, *srtdV, *interestPnts;
+    int i, j, validCount, cnt, end, k;
+    int iter, rows, cols;
+    F2D *temp;
+    int supIdPtr = 0;
+    
+    /** Concatenate x,y,v to form points matrix **/
+    r_sq = r * r;
+    n = points->height;
+
+    srtdVIdx = iMallocHandle(points->height, 1);
+    for (i = 0; i < srtdVIdx->height; i++) {
+        asubsref(srtdVIdx,i) = i;
+        }
+    srtdPnts = fMallocHandle (srtdVIdx->height, points->width);
+    for (i = 0; i < srtdVIdx->height; i++) {
+        for(j=0; j<points->width; j++) {
+            subsref(srtdPnts,i,j) = subsref(points, asubsref(srtdVIdx,i), j);
+                }
+        }
+    temp = fSetArray (1, 3, 0);
+    suppressR = fSetArray(n, 1, MAX_LIMIT);
+    validCount = n;
+    iter = 0;
+    
+    /** Allocate supId for #validCount and fill the values of
+        supId with the indices where suppressR>r_sq **/
+    k = 0;
+    supId = iMallocHandle(validCount, 1);
+    for (i = 0; i < (suppressR->height*suppressR->width); i++)
+    {
+        if ( asubsref(suppressR,i) > r_sq)
+        {
+            asubsref(supId,k++) = i;
+        }
+    }
+         
+    /** While number of features not-inspected is >0,  **/
+    while (validCount > 0)
+    {
+        F2D *tempp, *temps;
+        
+        /** Inspect the strongest feature point in srtdPnts
+            The index of that feature is in supId and the
+            index values in supId are arranged in descending order **/
+        asubsref(temp,0) = subsref(srtdPnts, asubsref(supId,0), 0);
+        asubsref(temp,1) = subsref(srtdPnts, asubsref(supId,0), 1);
+        asubsref(temp,2) = subsref(srtdPnts, asubsref(supId,0), 2);
+       
+        /** Stacking up the interestPnts matrix with top features
+            post suppression **/
+        if(iter == 0)
+            interestPnts = fDeepCopy(temp);
+        else
+        {
+            tempp = fDeepCopy(interestPnts);
+            fFreeHandle(interestPnts); 
+            interestPnts = ffVertcat(tempp, temp);
+            fFreeHandle(tempp);
+        }
+        iter++;
+ 
+        tempp = fDeepCopy(srtdPnts);
+        temps = fDeepCopy(suppressR);
+
+        fFreeHandle(srtdPnts);
+        fFreeHandle(suppressR);
+
+        /** Remove the feature that has been added to interestPnts **/
+        srtdPnts = fMallocHandle(supId->height-1, 3);
+        suppressR = fMallocHandle(supId->height-1, 1);
+        k=0;
+
+        for(i=1; i<(supId->height); i++)    /** Filling srtdPnts after removing the feature that was added to interestPnts**/
+        {
+            subsref(srtdPnts,(i-1),0) = subsref(tempp, asubsref(supId,i) ,0);
+            subsref(srtdPnts,(i-1),1) = subsref(tempp, asubsref(supId,i) ,1);
+            subsref(srtdPnts,(i-1),2) = subsref(tempp, asubsref(supId,i) ,2);
+            subsref(suppressR,(i-1),0) = subsref(temps, asubsref(supId,i) ,0);
+        }
+         
+        fFreeHandle(tempp);
+        fFreeHandle(temps);
+        rows = interestPnts->height-1;
+        cols = interestPnts->width;
+
+        /** For each feature, find how robust it is compared to the one in interestPnts **/
+        for (i = 0; i < srtdPnts->height; i++)
+            {
+            t = 0;
+                t1 = 0;
+
+                if ((C_ROBUST * subsref(interestPnts,rows,2)) >= subsref(srtdPnts, i,2))
+                {
+                        t = subsref(srtdPnts, i,0) - subsref(interestPnts,rows,0);
+                        t1 = subsref(srtdPnts, i,1) - subsref(interestPnts,rows,1);
+                        t = t * t + t1 * t1;
+                t1 = 0;
+            }
+
+                if ((C_ROBUST * subsref(interestPnts,rows,2)) < subsref(srtdPnts, i,2))
+                t1 = 1 * MAX_LIMIT;
+
+                if ( asubsref(suppressR, i) > (t + t1))
+            {
+                    asubsref(suppressR, i) = t + t1;
+            }  
+        }
+        
+        /** Inspect the new suppressR to find how many valid features left **/
+        validCount=0;
+        for (i = 0; i < suppressR->height; i++) {
+                if ( asubsref(suppressR,i) > r_sq) {
+                validCount++;
+                        } 
+                }
+        k = 0;
+        iFreeHandle(supId);
+        /** Allocate supId for #validCount and fill the values of
+        supId with the indices where suppressR>r_sq **/
+        supId = iMallocHandle(validCount, 1);
+        for (i = 0; i < suppressR->height*suppressR->width; i++) {
+            if ( asubsref(suppressR,i) > r_sq)
+                asubsref(supId,k++) = i;
+                }
+    }
+    
+    iFreeHandle(supId);
+    iFreeHandle(srtdVIdx);
+    fFreeHandle(srtdPnts);
+    fFreeHandle(temp);
+    fFreeHandle(suppressR);
+
+    return interestPnts;
+}
+
diff --git a/SD-VBS/benchmarks/tracking/src/c/getInterpolatePatch.c b/SD-VBS/benchmarks/tracking/src/c/getInterpolatePatch.c
new file mode 100644
index 0000000..99e95f9
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/getInterpolatePatch.c
@@ -0,0 +1,43 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+
+/** Perform simple interpolation around  2*winSize*2*winSize neighbourhood **/
+F2D* getInterpolatePatch(F2D* src, int cols, float centerX, float centerY, int winSize)
+{
+    F2D *dst;
+    float a, b, a11, a12, a21, a22;
+    int i, j, k, srcIdx, dstIdx;
+    int srcIdxx, dstIdxx;
+
+    a = centerX - floor(centerX);
+    b = centerY - floor(centerY);
+
+    a11 = (1-a)*(1-b);
+    a12 = a*(1-b);
+    a21 = (1-a)*b;
+    a22 = a*b;
+
+    dst = fSetArray(1,2*winSize*2*winSize, 0);
+
+
+    for(i=-winSize; i<winSize; i++)
+    {
+        srcIdxx = floor(centerY) + i;
+        dstIdxx = i+winSize;
+
+        for(j=-winSize; j<(winSize); j++)
+        {
+            srcIdx = srcIdxx * cols + floor(centerX) + j;
+            dstIdx = dstIdxx * 2 * winSize + j + winSize;
+            asubsref(dst,dstIdx) = asubsref(src,srcIdx)*a11 + asubsref(src,srcIdx+1)*a12 + asubsref(src,srcIdx+cols)*a21 + asubsref(src,srcIdx+1+cols)*a22;
+        }
+    }
+
+    return dst;
+}
+
+
+
diff --git a/SD-VBS/benchmarks/tracking/src/c/script_tracking.c b/SD-VBS/benchmarks/tracking/src/c/script_tracking.c
new file mode 100644
index 0000000..bb48ace
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/script_tracking.c
@@ -0,0 +1,263 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#include "tracking.h"
+#include <malloc.h>
+#include "extra.h"
+#define TRACKING_MEM 1<<29
+
+int main(int argc, char* argv[])
+{
+    SET_UP
+    mallopt(M_TOP_PAD, TRACKING_MEM);
+    mallopt(M_MMAP_MAX, 0);
+    int i, j, k, N_FEA, WINSZ, LK_ITER, rows, cols;
+    int endR, endC;
+    F2D *blurredImage, *previousFrameBlurred_level1, *previousFrameBlurred_level2, *blurred_level1, *blurred_level2;
+    F2D *verticalEdgeImage, *horizontalEdgeImage, *verticalEdge_level1, *verticalEdge_level2, *horizontalEdge_level1, *horizontalEdge_level2, *interestPnt;
+    F2D *lambda, *lambdaTemp, *features;
+    I2D *Ic, *status;
+    float SUPPRESION_RADIUS;
+    F2D *newpoints;
+
+    int numFind, m, n;
+    F2D *np_temp;
+
+    char im1[100];
+    int counter=2;
+    float accuracy = 0.03;
+    int count;
+    
+
+    N_FEA = 1600;
+    WINSZ = 4;
+    SUPPRESION_RADIUS = 10.0;
+    LK_ITER = 20;
+
+    #ifdef test
+        WINSZ = 2;
+        N_FEA = 100;
+        LK_ITER = 2;
+        counter = 2;
+        accuracy = 0.1;
+    #endif
+    #ifdef sim_fast
+        WINSZ = 2;
+        N_FEA = 100;
+        LK_ITER = 2;
+        counter = 4;
+    #endif
+    #ifdef sim
+        WINSZ = 2;
+        N_FEA = 200;
+        LK_ITER = 2;
+        counter = 4;
+    #endif
+    #ifdef sqcif
+        WINSZ = 8;
+        N_FEA = 500;
+        LK_ITER = 15;
+        counter = 2;
+    #endif
+    #ifdef qcif
+        WINSZ = 12;
+        N_FEA = 400;
+        LK_ITER = 15;
+        counter = 4;
+    #endif
+    #ifdef cif
+        WINSZ = 20;
+        N_FEA = 500;
+        LK_ITER = 20;
+        counter = 4;
+    #endif
+    #ifdef vga
+        WINSZ = 32;
+        N_FEA = 400;
+        LK_ITER = 20;
+        counter = 4;
+    #endif
+    #ifdef wuxga
+        WINSZ = 64;
+        N_FEA = 500;
+        LK_ITER = 20;
+        counter = 4;
+    #endif
+    #ifdef fullhd
+        WINSZ = 48;
+        N_FEA = 500;
+        LK_ITER = 20;
+        counter = 4;
+    #endif
+
+    I2D* images[counter];
+    /** Read input image **/
+    for(count=1; count<=counter; count++)
+    {
+        /** Read image **/
+        printf("Input image %d: ", count);
+        scanf("%s", im1);
+        images[count - 1] = readImage(im1);
+        if(count == 1) Ic = readImage(im1);
+    }
+
+
+    rows = Ic->height;
+    cols = Ic->width;
+    
+    printf("start\n");   
+    for_each_job{
+        
+        /** IMAGE PRE-PROCESSING **/
+    
+        /** Blur the image to remove noise - weighted avergae filter **/
+        blurredImage = imageBlur(Ic);
+    
+        /** Scale down the image to build Image Pyramid. We find features across all scales of the image **/
+        blurred_level1 = blurredImage;                   /** Scale 0 **/
+        blurred_level2 = imageResize(blurredImage);      /** Scale 1 **/
+    
+    
+        /** Edge Images - From pre-processed images, build gradient images, both horizontal and vertical **/
+        verticalEdgeImage = calcSobel_dX(blurredImage);
+        horizontalEdgeImage = calcSobel_dY(blurredImage);
+    
+        /** Edge images are used for feature detection. So, using the verticalEdgeImage and horizontalEdgeImage images, we compute feature strength
+            across all pixels. Lambda matrix is the feature strength matrix returned by calcGoodFeature **/
+        
+        lambda = calcGoodFeature(verticalEdgeImage, horizontalEdgeImage, verticalEdgeImage->width, verticalEdgeImage->height, WINSZ);
+        endR = lambda->height;
+        endC = lambda->width;
+        lambdaTemp = fReshape(lambda, endR*endC, 1);
+    
+        /** We sort the lambda matrix based on the strengths **/
+        /** Fill features matrix with top N_FEA features **/
+        fFreeHandle(lambdaTemp);
+        lambdaTemp = fillFeatures(lambda, N_FEA, WINSZ);
+        features = fTranspose(lambdaTemp);
+        
+        /** Suppress features that have approximately similar strength and belong to close neighborhood **/
+        interestPnt = getANMS(features, SUPPRESION_RADIUS);
+        
+        /** Refill interestPnt in features matrix **/
+        fFreeHandle(features);
+        features = fSetArray(2, interestPnt->height, 0);
+        for(i=0; i<2; i++) {
+            for(j=0; j<interestPnt->height; j++) {
+                subsref(features,i,j) = subsref(interestPnt,j,i); 
+                }
+        } 
+         
+    
+        fFreeHandle(verticalEdgeImage);
+        fFreeHandle(horizontalEdgeImage);
+        fFreeHandle(interestPnt);
+        fFreeHandle(lambda);
+        fFreeHandle(lambdaTemp);
+        /** Until now, we processed base frame. The following for loop processes other frames **/
+        for(count=1; count<=counter; count++)
+        {
+            /** Read image **/
+            //sprintf(im1, "%s/%d.bmp", argv[1], count);
+            //Ic = readImage(im1);
+            I2D* Icc = images[count-1];
+            rows = Icc->height;
+            cols = Icc->width;
+            
+        
+            /** Blur image to remove noise **/
+            blurredImage = imageBlur(Icc);
+            previousFrameBlurred_level1 = fDeepCopy(blurred_level1);
+            previousFrameBlurred_level2 = fDeepCopy(blurred_level2);
+            
+            fFreeHandle(blurred_level1);
+            fFreeHandle(blurred_level2);
+        
+            /** Image pyramid **/
+            blurred_level1 = blurredImage;
+            blurred_level2 = imageResize(blurredImage);
+        
+            /** Gradient image computation, for all scales **/
+            verticalEdge_level1 = calcSobel_dX(blurred_level1);   
+            horizontalEdge_level1 = calcSobel_dY(blurred_level1); 
+            
+            verticalEdge_level2 = calcSobel_dX(blurred_level2); 
+            horizontalEdge_level2 = calcSobel_dY(blurred_level2);
+            
+            newpoints = fSetArray(2, features->width, 0);
+                
+            /** Based on features computed in the previous frame, find correspondence in the current frame. "status" returns the index of corresponding features **/
+            status = calcPyrLKTrack(previousFrameBlurred_level1, previousFrameBlurred_level2, verticalEdge_level1, verticalEdge_level2, horizontalEdge_level1, horizontalEdge_level2, blurred_level1, blurred_level2, features, features->width, WINSZ, accuracy, LK_ITER, newpoints);
+            fFreeHandle(verticalEdge_level1);
+            fFreeHandle(verticalEdge_level2);
+            fFreeHandle(horizontalEdge_level1);
+            fFreeHandle(horizontalEdge_level2);
+            fFreeHandle(previousFrameBlurred_level1);
+            fFreeHandle(previousFrameBlurred_level2);
+            //printf("height = %d, width = %d, numFind = %d\n", newpoints->height, newpoints->width);
+            
+            /** Populate newpoints with features that had correspondence with previous frame features **/
+            np_temp = fDeepCopy(newpoints);
+            if(status->width > 0 )
+            {
+                k = 0;
+                numFind=0;
+                for(i=0; i<status->width; i++)
+                {
+                    if( asubsref(status,i) == 1)
+                        numFind++;
+                }
+                fFreeHandle(newpoints);
+                newpoints = fSetArray(2, numFind, 0);
+        
+                for(i=0; i<status->width; i++)
+                {
+                    if( asubsref(status,i) == 1)
+                    {
+                        subsref(newpoints,0,k) = subsref(np_temp,0,i);
+                        subsref(newpoints,1,k++) = subsref(np_temp,1,i);
+                    }
+                }    
+            }    
+            
+            iFreeHandle(status);
+            fFreeHandle(np_temp);
+            fFreeHandle(features);
+        
+            /** Populate newpoints into features **/
+            features = fDeepCopy(newpoints);
+            fFreeHandle(newpoints);
+        }
+    }
+    printf("end..\n");
+#ifdef CHECK   
+    /* Self checking */
+    {
+        int ret=0;
+        float tol = 2.0;
+#ifdef GENERATE_OUTPUT
+        fWriteMatrix(features, argv[1]);
+#endif
+        ret = fSelfCheck(features, "expected_C.txt", tol); 
+        if (ret == -1)
+            printf("Error in Tracking Map\n");
+    }
+#endif
+    
+    fFreeHandle(blurred_level1);
+    fFreeHandle(blurred_level2);
+    fFreeHandle(features);
+
+    for(count=1; count<=counter; count++)
+    {
+        free(images[count - 1] );
+    }
+    iFreeHandle(Ic);
+    WRITE_TO_FILE
+    return 0;
+
+}
+
+
+
diff --git a/SD-VBS/benchmarks/tracking/src/c/tracking.h b/SD-VBS/benchmarks/tracking/src/c/tracking.h
new file mode 100644
index 0000000..24b8b2d
--- /dev/null
+++ b/SD-VBS/benchmarks/tracking/src/c/tracking.h
@@ -0,0 +1,20 @@
+/********************************
+Author: Sravanthi Kota Venkata
+********************************/
+
+#ifndef _SCRIPT_TRACK_
+#define _SCRIPT_TRACK_
+
+#include "sdvbs_common.h"
+
+F2D* calcAreaSum(F2D* src, int cols, int rows, int winSize);
+F2D* calcGoodFeature(F2D* dX, F2D* dY, int cols, int rows, int winSize);
+I2D* calcPyrLKTrack(F2D* ip1, F2D* ip2, F2D* idxp1, F2D* idxp2, F2D* idyp1, F2D* idyp2, F2D* jp1, F2D* jp2, F2D* fPnt, int nFeatures, int winSize, float accuracy, int maxIter, F2D* newPnt);
+F2D *getANMS (F2D *points, float r);
+F2D* getInterpolatePatch(F2D* src, int cols, float centerX, float centerY, int winSize);
+float polynomial(int d, F2D* a, F2D* b, int dim);
+I2D* sortInd(F2D* input, int dim);
+F2D* fillFeatures(F2D* lambda, int N_FEA, int win);
+
+#endif
+