|| || || || || || || ||

---

Sift.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_Sift_h
#define Impala_Core_Array_Sift_h

#include "Core/Array/Octaves.h"
#include "Core/Keypoint/DoG.h"


namespace Impala{
namespace Core{
namespace KeyPoint{

class Sift{

private:
    ILOG_VAR_DEC;
    Array2dScalarUInt8* mSrc;
    Detector* mDetector;
    DetectorType mDetectorType;
    Array2dScalarReal64** mGradMag;
    Array2dScalarReal64** mGradAng;
    Array2dScalarReal64** mGaussians;
    Array2dScalarReal64* mGaussian;
    KeyPointTableType* mKeyPoints;
    FeatureTable* mFeatureTable;
    int mLevelCnt;
    int mOctaveCnt;
    int mDescriptorWindowSize;
    int mDescriptorBinCount;
    int mDescriptorRegDiv;
    int mDescriptorLength;
    
    Octaves* mOctaves;
public:


    Sift(Array2dScalarUInt8* src, int Octaves, int Levels, KeyPoint::DetectorType dT = DOG )
    :mSrc(src)
    {
        mDetectorType = dT;
        Init();
    }
    Octaves* GetOctaves(){
        return mOctaves;
    }
    Detector* GetDetector(){
        return mDetector;    
    }

    void Init(){
        CmdOptions& options = CmdOptions::GetInstance();
        //The window size taken for the desccriptor
        mDescriptorWindowSize=options.GetInt("dws",16);
        //Default is 8 bins from every region
        mDescriptorBinCount=options.GetInt("dbc",8);
        //Default is a 4x4 region
        mDescriptorRegDiv=options.GetInt("drv",4);
        ILOG_INFO("Descriptor Window Size : "<<mDescriptorWindowSize);
        ILOG_INFO("Descriptor Histogram Bin Count : "<<mDescriptorBinCount);

        mDescriptorLength = mDescriptorBinCount *
                            mDescriptorRegDiv   * mDescriptorRegDiv;


        mOctaves = new Impala::Core::Array::Octaves(mSrc,-1,3);
        switch(mDetectorType){
            case DOG:
                mDetector = new DoG(mOctaves);
        }
        mDetector->FindKeyPoints();
        mLevelCnt = mOctaves->GetLevelCount();
        mOctaveCnt = mOctaves->GetOctaveCount();


        mGradMag = new Array2dScalarReal64*[mOctaveCnt*mLevelCnt];
        mGradAng = new Array2dScalarReal64*[mOctaveCnt*mLevelCnt];
        mGaussians = new Array2dScalarReal64*[mOctaveCnt*mLevelCnt];
        for(int i=0;i<mOctaveCnt*mLevelCnt;i++){
            mGaussians[i] =0;
            mGradMag[i] = 0;
            mGradAng[i] = 0;
        }

        BuildGaussians();
        mGaussian = MakeGaussian2d( 1.5 * mDescriptorWindowSize, 
                                    1.5 * mDescriptorWindowSize,0,0,3);
        mKeyPoints=new KeyPointTableType(0);
        
    }
    void ExtractDescriptors()
    {
        ComputeGradientField();
        AssignOrientation();
        Extract();
    }

    Array2dScalarReal64* GetGradMag(int o,int l)
    {
        if(l>mLevelCnt-1){
            ILOG_INFO("Clipping, only "<<mLevelCnt<< " gradmag. at each octave!");
            l=mLevelCnt-1;
        }
        return mGradMag[o*mLevelCnt+l];
    }
    
    Array2dScalarReal64* GetGradAng(int o,int l)
    {
        if(l>mLevelCnt-1){
            ILOG_INFO("Clipping, only "<<mLevelCnt<< " gradang. at each octave!");
            l=mLevelCnt-1;
        }
        return mGradAng[o*mLevelCnt+l];
    }
    void ComputeGradientField(){
        double H[4];
        double G[2];
        double S[2];

        ILOG_INFO("Computing gradient field");

        //For every level, construct the gradient magnitude and angle map
        for(int o=0;o<mOctaveCnt;o++){
            for(int l=0;l<mLevelCnt;l++)
            {
                Array2dScalarReal64* L = mOctaves->GetOctave(o,l);
                Set(mGradMag[o*mLevelCnt+l],L);
                SetVal(mGradMag[o*mLevelCnt+l],0);
                Set(mGradAng[o*mLevelCnt+l],mGradMag[o*mLevelCnt+l]);
                //dd stands for derivative distance..
                int dd;
                //The sigma to be calculated in the gradient comes from 
                //the difference of gaussians, for which peaks are analysed
                //in the sigma neighborhood. because of these two shifts, we
                //have to use l+2
                Real64 sigma=mOctaves->GetSigma(o,l+2);
                // If resampling is enabled, we have already reduced 
                // the scale by two at every octave(doubling of sigma). 
                // Therefore the derivative is the immediate neighbors.
                //
                // Otherwise jump as much as 1 sigma for a derivative,
                // since the images are not scaled down as well.
                if(mOctaves->IsResamplingEnabled())
                    dd=1;
                else
                    dd=sigma;
                
                for(int i=dd;i<L->CW()-dd;i++){
                for(int j=dd;j<L->CH()-dd;j++){
                    //this is the first derivative wrt. x
                    //[ -1 0 1 ]
                    G[0] = L->Value(i+dd,j) - L->Value(i-dd,j);

                    //this is the first derivative wrt. y
                    G[1] = L->Value(i,j+dd) - L->Value(i,j-dd);
                    
                    mGradMag[o*mLevelCnt+l]->SetValue(sqrt(G[0]*G[0]+G[1]*G[1])
                                                   ,i,j);
                    mGradAng[o*mLevelCnt+l]->SetValue(atan(G[1]/G[0]),i,j);


                }
                }
                /*
                ILOG_DEBUG("Octave("<<o<<","<<l<<") : ("
                            <<PixMin( mGradMag[o*mLevels+l]) <<","
                            <<PixMax( mGradMag[o*mLevels+l]) <<") ("
                            <<PixMin( mGradAng[o*mLevels+l]) <<","
                            <<PixMax( mGradAng[o*mLevels+l]) <<")");*/
                Normalize(mGradMag[o*mLevelCnt+l]);

            }
        }
    }
    void Normalize(Array2dScalarReal64* s){
        Real64 min,max;
        PixMinMax(s,&min,&max);
        AddVal(s,s,-min);
        MulVal(s,s,1/(max-min));
    }

    void BuildGaussians(){
        Real64 sigma;
        int index;
        for (int o=0;o<mOctaveCnt;o++){
            for(int l=0;l<mLevelCnt;l++)
            {
                sigma = mOctaves->GetSigma(o,l+2);
                index = o*mLevelCnt+l;

                mGaussians[index] = MakeGaussian2d(1.5 * sigma,
                                                     1.5 * sigma,0,0,3);
                Normalize(mGaussians[index]);
                /*
                ILOG_DEBUG("Gaussian Created");
                ILOG_DEBUG("sigmas = ( "<<sigma<<" , "<<sigma*1.5<<" )");
                ILOG_DEBUG("sizes  = ( "<<mGaussians[index]->CW()<<" x "
                                        <<mGaussians[index]->CH()<<" )");
                                        */
            
            }
        
        }
    }
    void AssignOrientation()
    {
        
        double sigma;
        int range=0;
        int k=0;
        int o,l,i,j,index;
        Array2dScalarReal64 *GA, *GM;
        Array2dScalarReal64 *Gaussian;

        mKeyPoints=mDetector->GetKeyPoints();
        int OrigSize = mKeyPoints->Size();
        while(k<OrigSize)
        {
            sigma= mKeyPoints->Get1(k);
            i = mKeyPoints->Get2(k);
            j = mKeyPoints->Get3(k);
            o = mKeyPoints->Get4(k);
            l = mKeyPoints->Get5(k);
            index=o*mLevelCnt+l;

            GA=mGradAng[index];
            GM=mGradMag[index];
            Gaussian=mGaussians[index];
            //Weight the gradient magnitudes with a gaussian,
            // with a sigma 1.5 times the original,
            // and centered around the keypoint
            range = Gaussian->CW()/2;
            
            double histAng[36];
            for(int h=0;h<36;h++)
                histAng[h]=0;
            double incr=PI/36;

            for(int x=-range;x<range+1;x++)
            {
                for(int y=-range;y<range+1;y++)
                {
                    if((i+x<0)||(i+x>GA->CW())||
                       (j+y<0)||(j+y>GA->CH()))
                        continue;
                    if((x+range>=Gaussian->CW()||y+range>=Gaussian->CH()))
                        continue;
                    int bin=floor(18+GA->Value(i+x,j+y)/incr);
                    if(bin>36)
                        continue;
                    histAng[bin]+=Gaussian->Value(x+range,y+range) *
                                               GM->Value(i+x,j+y);
                }
            
            }
            //Find the maximum in the histogram
            double max_val=0;
            int max_bin=0;
            for(int b=0;b<36;b++)
            {
                if(histAng[b]>max_val)
                {
                    max_val = histAng[b];
                    max_bin=b;
                }
            }
            mKeyPoints->Set6(k,-PI/2+incr*(max_bin+.5));
            //Update the Orientation of the keypoint to max_bin
            double low,high;
            low=-PI/2+incr*max_bin;
            high=-PI/2+incr*(max_bin+1);
            //ILOG_DEBUG("Orientation("<<max_bin<<"):"<<histAng[max_bin]
            //        <<" Values btw "<<low<<"-"<<high);

            //Find other values which are greater than 80% of highest peak
            double max80= .8 * max_val;
            for(int b=1;b<35;b++)
            {
                if(b==max_bin)
                    continue;
                if(histAng[b]>max80)
                {   
                    //Check if it is a peak
                    if((histAng[b]>histAng[b-1])&&(histAng[b]>histAng[b+1])){
                        //Create a new keypoint at the same location 
                        //with this orientation
                        //
                        //ILOG_DEBUG("Create New Keypoint With Orientation:"<<b<<"-"<<histAng[b]);
                        mKeyPoints->Add(sigma,i,j,o,l,-PI/2+incr*(b+.5));
                    }
                }
            }
            //find 3 histogram values closest to each peak, and fit a parabola
            //for better accuracy
            
            k++;
        }
        ILOG_INFO("Total "<<mKeyPoints->Size()-OrigSize
                        <<" new orientations assigned!");
        return;
    }
    void Extract()
    {
        ILOG_DEBUG("Extracting Descriptors")
        int k=0;

        Array2dScalarReal64 *GM,*GA,*Gaussian;
        int o,l,i,j,index;
        Real64 ori, sigma;
        int range = mGaussian->CW()/2;
        int range_ext= range*sqrt(2)+1;
        KeyPointTableType* oldKeyPoints=mKeyPoints;
        int size = oldKeyPoints->Size();
        int current = 0;
        mFeatureTable = new FeatureTable(FeatureDefinition("SIFT"),size,mDescriptorLength);


        mKeyPoints=new KeyPointTableType(0);
        while(k<size)
        {
            if(!((oldKeyPoints->Get2(k)==0)&&(oldKeyPoints->Get3(k)==0))){
                sigma= oldKeyPoints->Get1(k);
                i = oldKeyPoints->Get2(k);
                j = oldKeyPoints->Get3(k);
                o = oldKeyPoints->Get4(k);
                l = oldKeyPoints->Get5(k);
                ori = oldKeyPoints->Get6(k);
                index=o*mLevelCnt+l;
                mKeyPoints->Add(sigma,i,j,o,l,ori);
                GA=mGradAng[index];
                GM=mGradMag[index];

                //Weight the gradient magnitudes with a gaussian,
                // with a sigma 1.5 * descriptor window size
                // and centered around the keypoint

                
                //Rotate
                Array2dScalarReal64 *srcGA=ArrayCreate<Array2dScalarReal64>
                                (mDescriptorWindowSize,mDescriptorWindowSize);
                Array2dScalarReal64 *srcGM=ArrayCreate<Array2dScalarReal64>
                                (mDescriptorWindowSize,mDescriptorWindowSize);
                Array2dScalarReal64 *dstGA=0; 
                Array2dScalarReal64 *dstGM=0; 
                SetVal(srcGA,0);
                SetVal(srcGM,0);


                //mHistogram Bin Count
                double histAng[16*mDescriptorBinCount];
                for(int h=0;h<16*mDescriptorBinCount;h++)
                    histAng[h]=0;
                double incr=PI/mDescriptorBinCount;
                bool clipped=false;

                for(int x = -mDescriptorWindowSize/2;
                        x < mDescriptorWindowSize/2; x++)
                {
                    for(int y =-mDescriptorWindowSize/2;
                            y <mDescriptorWindowSize/2; y++)
                    {
                        if((i+x<0)||(i+x>GA->CW())||
                           (j+y<0)||(j+y>GA->CH())){
                           clipped=true;
                            continue;
                        }
                        srcGA->SetValue(GA->Value(i+x,j+y),
                                      x+mDescriptorWindowSize/2,
                                      y+mDescriptorWindowSize/2);
                        srcGM->SetValue(GM->Value(i+x,j+y),
                                      x+mDescriptorWindowSize/2,
                                      y+mDescriptorWindowSize/2);
                        
                    }
            
                }

                Rotate(dstGA,srcGA,ori*180/PI,Geometry::LINEAR,true,0);
                AddVal(dstGA,dstGA,ori);

                Rotate(dstGM,srcGM,ori*180/PI,Geometry::LINEAR,true,0);

                for(int hx=0;hx<mDescriptorRegDiv;hx++){//hx and hy are divs
                for(int hy=0;hy<mDescriptorRegDiv;hy++){
                for(int ix=0;ix<dstGA->CW()/mDescriptorRegDiv;ix++){
                    for(int iy=0;iy<dstGA->CH()/mDescriptorRegDiv;iy++){
                        int x=hx*mDescriptorRegDiv+ix;
                        int y=hy*mDescriptorRegDiv+iy;
                        if(dstGA->Value(x,y)>PI/2)
                            dstGA->SetValue(dstGA->Value(x,y)-PI/2,x,y);
                        if(dstGA->Value(x,y)<0)
                            dstGA->SetValue(dstGA->Value(x,y)+PI/2,x,y);

                        int bin=floor(mDescriptorBinCount/2+
                              dstGA->Value(x,y)/incr);

                        if(bin>mDescriptorBinCount)
                            continue;
                        if(bin<0)
                            continue;
                        //The descriptor is stored as 4x4 regions, and 8
                        //histogram bins

                        histAng[mDescriptorRegDiv*mDescriptorBinCount*hx+
                                mDescriptorBinCount*hy+
                                bin]
                                += mGaussian->Value(x,y) * dstGM->Value(x,y);

                    }
                    
                }
                }
                }
                mFeatureTable->GetColumn1()->Set(current++,MakeQuid(0,0,0,k));
                //the size of the feature depends on the number of binws in the
                //histograms, and the square of the region( default:8bins-4 x 4)
                mFeatureTable->GetColumn2()->AddVector( histAng,
                                                        mDescriptorLength );

            }
           k++; 
        }
        ILOG_INFO("Total Keypoint Count:"<<mKeyPoints->Size());
        
    
    }
    KeyPointTableType* GetKeypoints(){
        return mKeyPoints;
    } 



};

ILOG_VAR_INIT(Sift, Impala.Core.Array);

}}}//Namespace



#endif

---