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

---

DoG.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_DoG_h
#define Impala_Core_Array_DoG_h

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


namespace Impala{
namespace Core{
namespace KeyPoint{

//using namespace Impala::Core::Array;

class DoG:public Detector{
    
private:
    ILOG_VAR_DEC;
    Array::Octaves* mOctaves;
    Array2dScalarReal64** mDoG;
    Array2dScalarReal64** mLevelExtrema;
    Array2dScalarReal64** mOctaveExtrema;

    int mOctaveCnt;
    int mLevelCnt;
    int mScaleCnt;
    int mCancelled;

public:
    DoG(Array::Octaves* o)
    {
        mOctaves=o;
        mOctaveCnt=o->GetOctaveCount();
        mLevelCnt=o->GetLevelCount();
        mScaleCnt=o->GetScaleCount();
        Init();
        
    }
    ~DoG(){
        for(int i=0; i<mOctaveCnt*(mScaleCnt-1);i++)
        {
            delete mDoG[i];
        }
        delete [] mDoG;
        for(int i=0; i<mOctaveCnt*mLevelCnt;i++)
        {
            delete mLevelExtrema[i];
        }
        delete [] mLevelExtrema;
        for(int i=0; i<mOctaveCnt;i++)
        {
            delete mOctaveExtrema[i];
        }
        delete [] mOctaveExtrema;
    }
    void Init(){
        mCancelled=0;

        //mDoG holds the difference of gaussians between successive levels
        mDoG = new Array2dScalarReal64*[mOctaveCnt*(mScaleCnt-1)];
        for(int i=0;i<mOctaveCnt*(mScaleCnt-1);i++)
            mDoG[i]=0;
        

        //mLevelExtrema holds the max/min of successive levels(mScaleCnt-3)
        //mGradMag holds the gradient magnitude of different levels
        //mGradAng holds the gradient orientation of different levels
        mLevelExtrema= new Array2dScalarReal64*[mOctaveCnt*mLevelCnt];

        for(int i=0;i<mOctaveCnt*mLevelCnt;i++){
            mLevelExtrema[i] = 0;
        }
        
        mOctaveExtrema = new  Array2dScalarReal64*[mOctaveCnt];
        for(int i=0;i<mOctaveCnt;i++)
            mOctaveExtrema[i]=0;
        mKeyPoints=new KeyPointTableType(0);
    
    }
    KeyPointTableType* GetKeyPoints(){
        return mKeyPoints;
    }
    
    Array2dScalarReal64* GetDoG(int o,int l)
    {
        if(l>mScaleCnt-2){
            ILOG_INFO("Clipping, only "<<mScaleCnt-1<< " DoG at each octave!");
            l=mScaleCnt-2;
        }
        return mDoG[o*(mScaleCnt-1)+l];
    }
    
    Array2dScalarReal64* GetLevelExtrema(int o,int l)
    {
        if(l>mLevelCnt-1){
            ILOG_INFO("Clipping, only "<<mLevelCnt<< " extrema at each octave!");
            l=mLevelCnt-1;
        }
        return mLevelExtrema[o*mLevelCnt+l];
    }
    Array2dScalarReal64* GetOctaveExtrema(int o)
    {
        if(o>mOctaveCnt){
            ILOG_INFO("Clipping, only "<<mOctaveCnt<< " octaves!");
            o=mOctaveCnt;
        }
        return mOctaveExtrema[o];
    }

    void FindKeyPoints()
    {
        CalculateDoG();
        //CalculateDoGExtrema();
        CalculateOctaveExtrema();
        RejectLowContrast();
        RejectEdgeResponse();
        AcquireKeypoints();
        //Filter points
        LocalizeKeypoints();
    
    
    }
    void CalculateDoG(){
        ILOG_INFO("Calculating DoG ...");
        for(int o=0;o<mOctaveCnt;o++){
            for(int l=0;l<mScaleCnt-1;l++){
                Sub(mDoG[o*(mScaleCnt-1)+l],
                    mOctaves->GetOctave(o,l),
                    mOctaves->GetOctave(o,l+1));
            }
        }
    }

    void CalculateOctaveExtrema()
    {
        ILOG_INFO("Calculating Octave Extrema ...");
        for(int o=0;o<mOctaveCnt;o++)
        {
            mOctaveExtrema[o] = 0;
            for(int l=0;l<mLevelCnt;l++)
            {
                
                mLevelExtrema[o*mLevelCnt+l]=
                    GetLSSE(mDoG[o*(mScaleCnt-1)+l],
                            mDoG[o*(mScaleCnt-1)+l+1],
                            mDoG[o*(mScaleCnt-1)+l+2]);
                if(l==0)
                    Set(mOctaveExtrema[o],mLevelExtrema[o*mLevelCnt+l]);
                else
                    Add(mOctaveExtrema[o],mOctaveExtrema[o],
                            mLevelExtrema[o*mLevelCnt+l]);
            }
        }
    }
    Array2dScalarReal64* GetLSSE(Array2dScalarReal64* u,
                Array2dScalarReal64* c,
                Array2dScalarReal64* l)
    {
        Array2dScalarReal64* res = 0 ;
        Array2dScalarReal64* arrays[3];

        arrays[0]=u;
        arrays[1]=c;
        arrays[2]=l;

        Set(res,c);
        SetVal(res,0);
        /*
        ILOG_DEBUG("LSSE sizes: "<<l->CH()<<"x"<<l->CW()
                                <<" | "<<c->CH()<<"x"<<c->CW()
                                <<" | "<<u->CH()<<"x"<<u->CW());
                                */
        Real64  max;
        int i,j,m,n,a;
        bool is_max;
        for(i=1;i<c->CW()-1;i++){
            for(j=1;j<c->CH()-1;j++)
            {
                max = arrays[1]->Value(i,j);
                is_max = true;
                for(a=0;a<3;a++)
                {
                    for(m=-1;m<2;m++)
                    {
                        for(n=-1;n<2;n++)
                        {
                            if(max <= arrays[a]->Value(i+m,j+n))
                            {
                                if((a==1)&&(m==0)&&(n==0))
                                    continue;
                                //ILOG_DEBUG(a<<"\t"<<m<<"\t"<<n);
                                is_max=false;
                                break;  
                            }
                        }
                        if(!is_max)
                            break;
                    }
                    if(!is_max)
                        break;
                }
                if(is_max)
                    res->SetValue(1,i,j);
            }
        }
       return res;

    }

    void RejectLowContrast()
    {
        ILOG_WARN("Reject low contrast not implemented yet!");
    }
    void RejectEdgeResponse()
    {
        ILOG_WARN("Reject edge response not implemented!");
        return;    
    }
    void AcquireKeypoints()
    {
        ILOG_INFO("Acquiring Keypoints");
        for(int o=0;o<mOctaveCnt;o++)
        {
            for(int l=0;l<mLevelCnt;l++)
            {
                Array2dScalarReal64* arr = mLevelExtrema[o*mLevelCnt+l];
                for(int i=0;i<arr->CW();i++){
                    for(int j=0;j<arr->CH();j++)
                    {
                        if(arr->Value(i,j)!=0)
                        {
                            mKeyPoints->Add(mOctaves->GetSigma(o,l+2),i,j,o,l,0);
                        }
                    }
                }
            }
        }
        ILOG_INFO("Total "<<mKeyPoints->Size()<<" keypoints!");

    }

    void LocalizeKeypoints()
    {   
        int k = 0 ;
        int OrigSize=mKeyPoints->Size();
        ILOG_INFO("Localizing "<<OrigSize<<" Keypoints");
        int new_i;
        int new_j;
        int last=0;
        int count=0;
        mCancelled=0;
        double H[4];
        double G[2];
        double S[2];
        Real64 Dx=0;
        Array2dScalarReal64* Octave;
        //form the Hessian matrix &Gradient vector around each point, 
        //and find the shift from 
        // -inv((Hessian(A))*Grad(A)) : evaluated at the keypoint
        while(k<OrigSize)
        {
            Real64 sigma=mKeyPoints->Get1(k);
            int i = mKeyPoints->Get2(k);
            int j = mKeyPoints->Get3(k);
            int o = mKeyPoints->Get4(k);
            int l = mKeyPoints->Get5(k);


            Octave=mOctaves->GetOctave(o,l);

            int Width=mOctaves->GetOctave(o,l)->CW();
            int Height=mOctaves->GetOctave(o,l)->CH();


            //dd stands for derivative distance..something I just made up now-)
            int dd;

            // 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 we have to jump as much as one sigma for a derivative,
            // since the images are not scaled down as well.
            if(mOctaves->IsResamplingEnabled())
                dd=1;
            else
                dd=sigma;
                
            //ILOG_DEBUG("Using "<<dd<<" as der.distance @ ("<<o<<","<<l<<")");
            if((i+2*dd>=Width)||(j+2*dd>=Height)||(i-2*dd<0)||(j-2*dd<0))
            {
                //ILOG_WARN("Should cancel the keypoint");
                mCancelled++;
                mKeyPoints->Set2(k,0);
                mKeyPoints->Set3(k,0);
                count=0;
                k++;
                last=k;
                continue;
            }
                
        

            //Second row of Hessian
            //this is the first derivative wrt. x
            //[ -1 0 1 ]
            G[0] = Octave->Value(i+dd,j) - Octave->Value(i-dd,j);


            //second derivative wrt. x
            //[ 1 0 -2 0 1 ]
            H[0]= -2*Octave->Value(i,j)  +
                    Octave->Value(i-2*dd,j) +
                    Octave->Value(i+2*dd,j) ;
            /* Second Derivative
            [ 1  0 -1]
            [ 0  0  0]
            [-1  0  1]*/
            //second derivative wrt. y
            H[1]=    Octave->Value(i-dd,j-dd) +
                     Octave->Value(i+dd,j+dd) -
                     Octave->Value(i-dd,j+dd) -
                     Octave->Value(i+dd,j-dd) ;

            //Third row of Hessian
            //this is the first derivative wrt. y
            G[1] = Octave->Value(i,j+dd) - Octave->Value(i,j-dd);

            //second derivative wrt. y
            //[ 1 0 -2 0 1 ]
            H[2]= -2*Octave->Value(i,j)  +
                    Octave->Value(i,j-2*dd) +
                    Octave->Value(i,j+2*dd) ;

            /* Second Derivative
            [ 1  0 -1]
            [ 0  0  0]
            [-1  0  1]*/
            //second derivative wrt. y
            H[3]=   Octave->Value(i-dd,j-dd) +  Octave->Value(i+dd,j+dd) -
                    Octave->Value(i-dd,j+dd) -  Octave->Value(i+dd,j-dd) ;

            Real64 div = H[0]*H[3]-H[1]*H[2];
            S[0]=(-H[3]*G[0]+H[1]*G[1])/div;
            S[1]=(H[2]*G[0]-H[0]*G[1])/div;
    
            Dx= Octave->Value(i,j)+.5*(S[0]*G[0]+S[1]+G[1]);
            if(Dx<.03){
                //ILOG_WARN("Should cancel the keypoint");
                mCancelled++;
                mKeyPoints->Set2(k,0);
                mKeyPoints->Set3(k,0);
                count=0;
                k++;
                last=k;
                continue;
            }
            //advance to the next keypoint
            if(!((S[0]> 0.5)||(S[1]> 0.5)||(S[0]<-.5)||(S[1]<-.5))){
               k++;
               last=k;
               count=0;
            }else{
                new_i=i+S[0];
                new_j=j+S[1];
                if((new_i>0)&&(new_i<Width)){
                    mKeyPoints->Set2(k,new_i);
                }else{
                    //ILOG_INFO("Keypoint "<<k<<" : Location out of border");
                    count=0;
                    k++;
                    last=k;
                    continue;
                }


                if((new_j>0)&&(new_j<Height)){
                    mKeyPoints->Set3(k,new_j);
                }else{
                    //ILOG_INFO("Keypoint "<<k<<" : Location out of border");
                    count=0;
                    k++;
                    last=k;
                    continue;
                }

                if(last==k)
                    count++;
                /*
                ILOG_DEBUG("Keypoint "<<k<<" : Location shifted!("<<count<<")");
                ILOG_DEBUG("Estimated Shift : "<<S[0]<<"x"<<S[1]);
                */

                if(count>10){
                    k++;
                    last=k;
                    count=0;
                }
            }

        }
        ILOG_INFO("Cancel count:"<<mCancelled<<"/"<<OrigSize);

    }

    


};
ILOG_VAR_INIT(DoG, Impala.Core.Array);

}}}
#endif

---