1 files changed, 185 insertions, 0 deletions

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;
+}