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

---

HarrisLaplaceDetector.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_HarrisLaplaceDetector_h
#define Impala_Core_Feature_HarrisLaplaceDetector_h

#include "Core/Array/Pattern/PatInOutOp.h"
#include "Core/Array/Abs.h"
#include "Core/Array/PixMax.h"
#include "Core/Array/Dilation.h"
#include "Core/Array/Resample.h"
#include "Core/Array/GaussDerivative.h"
#include "Core/Array/RecGauss.h"
#include "Core/Array/Trait/BpoEqual.h"
#include "Core/Array/Harris.h"
#include "Core/Array/RGB2Gray.h"
#include "Core/Array/Rgb2Ooo.h"
#include "Core/Geometry/Rectangle.h"
#include "Core/Array/Laplacian.h"
#include "Basis/Timer.h"

#include "Core/Feature/PointDescriptorTable.h"
#include "Core/Table/Sort.h"

#include <iostream> 
#include <fstream>
#include <algorithm>
#include <iterator>

namespace Impala
{
namespace Core
{
namespace Feature
{


class HarrisLaplaceDetector
{
public:
    typedef Array::Array2dScalarReal64 Array2dScalarReal64;
    typedef Array::Array2dVec3UInt8 Array2dVec3UInt8;
    typedef Array::Array2dVec3Real64 Array2dVec3Real64;

    // command line parameters
    Real64      mPrecision;
    Real64      mSigmaA;
    bool        mUseRecGauss;
    bool        mBackwardsCompatible;
    
    HarrisLaplaceDetector(bool useRecGauss=false) : 
        mPrecision(5.0),     // needs to be this high for better Matlab correspondance
        mSigmaA(3.0), 
        mUseRecGauss(useRecGauss),
        mBackwardsCompatible(true)
    { 
    }

    class MyPoint 
    {
    public:
        Real64 mDetectionScale;
        Real64 mScale;
        Real64 mCornerness;
        int mX;
        int mY;
    
        MyPoint(int x, int y, Real64 cornerness)
            : mX(x), mY(y), mCornerness(cornerness)
        {
        }
    };
    
    typedef std::vector<MyPoint*> MyPointList;
    typedef std::back_insert_iterator<MyPointList> MyPointListBI;

    class ExportPoints
    {
    public:
        typedef Array::Pattern::TagPixOpOut       DirectionCategory;
        typedef Array::Pattern::TagTransVar       TransVarianceCategory;
        typedef Array::Pattern::Tag1Phase         PhaseCategory;
    
        ExportPoints(MyPointListBI& pListBI,
                     Real64 harrisThreshold) :
            mHarrisThreshold(harrisThreshold), mPListBI(pListBI)
        {
        }
    
        void
        DoIt(const Real64& s, int x, int y)
        {
            if(s > mHarrisThreshold) {
                *mPListBI++ = new MyPoint(x, y, s);
            }
        }
    
    private:
        Real64 mHarrisThreshold;
        MyPointListBI& mPListBI;
    };
    
    // finds max points and stores result in dst as well as pList
    void
    DoMaxPoints(Array2dScalarReal64* src,
                MyPointListBI pListBI,
                Real64 harrisThreshold)
    {
        using namespace Impala::Core::Array;
        int cw = 2*mSigmaA + 1;
        int ch = 2*mSigmaA + 1;
        Array2dScalarReal64* ker = ArrayCreate<Array2dScalarReal64>(cw, ch);
        Core::Matrix::MatSet(ker, 0);
    
        Array2dScalarReal64* tmp1 = 0;
        Dilation(tmp1, src, ker);
        Trait::BpoEqual<Array2dScalarReal64, Array2dScalarReal64,
                        Array2dScalarReal64> bpo(0.0000000000000001);
        Pattern::PatBinaryPixOp(tmp1, tmp1, src, bpo); // src == dilated implies local max
        Mul(tmp1, tmp1, src); // obtain orignal values for local max from src
        Real64 theMax = PixMax(tmp1);
        //std::cout << "theMax : " << theMax << std::endl;
    
        ExportPoints expPoints(pListBI, harrisThreshold); // obtain coordinates
        Pattern::PatInOutOp(tmp1, expPoints);
    
        delete tmp1;
        delete ker;
    }
    
    void
    harrisCornerDetector(MyPointList& pointList,
                         Array2dScalarReal64* imIntensity,
                         Real64 sigmaN, Real64 harrisK, Real64 harrisThreshold) 
    {
        Real64 sigmaIntegr = sigmaN;
        Real64 sigmaDiff = 0.7125 * sigmaN;
    
        // compute Lx and Ly Gaussian derivatives
        Array2dScalarReal64* Lx = 0;
        Array2dScalarReal64* Ly = 0;
        if(mUseRecGauss)
        {
            RecGauss(Lx, imIntensity, sigmaDiff, sigmaDiff, 1, 0, mPrecision);
            RecGauss(Ly, imIntensity, sigmaDiff, sigmaDiff, 0, 1, mPrecision);
        } else {
            GaussDerivative(Lx, imIntensity, sigmaDiff, 1, 0, mPrecision);
            GaussDerivative(Ly, imIntensity, sigmaDiff, 0, 1, mPrecision);
        }
        
        // VERY IMPORTANT: SCALE MULTIPLIER!
        Array2dScalarReal64* Lx2 = 0;    // aka mu11
        Mul(Lx2, Lx, Lx);
        MulVal(Lx2, Lx2, sigmaDiff * sigmaDiff);
    
        Array2dScalarReal64* Ly2 = 0;    // aka mu22
        Mul(Ly2, Ly, Ly);
        MulVal(Ly2, Ly2, sigmaDiff * sigmaDiff);
    
        Array2dScalarReal64* LxLy = 0;   // aka mu12/mu21
        Mul(LxLy, Lx, Ly);
        MulVal(LxLy, LxLy, sigmaDiff * sigmaDiff);
    
        delete Lx;
        delete Ly;
    
        // compute Harris  
        Array2dScalarReal64* harrisOutput = 0;
        Array::Harris<Array2dScalarReal64>(harrisOutput, Lx2, Ly2, LxLy,
                                           sigmaIntegr, harrisK, mUseRecGauss);
    
        delete Lx2;
        delete Ly2;
        delete LxLy;
    
        DoMaxPoints(harrisOutput, std::back_inserter(pointList), harrisThreshold);
        
        // write the Harris cornerness to an image on disk
        /*std::string output = "debug.png";
        Array2dVec3UInt8* buf = ArrayCreate<Array2dVec3UInt8>
            (harrisOutput->CW(), harrisOutput->CH());
        GetRgbPixels(harrisOutput, buf->CPB(), std::string("Stretch"));
    
        WritePng(buf, output);
        delete buf;*/
        
        delete harrisOutput;
    }

    // multi-channel version
    template<class SrcArrayT>
    void
    harrisCornerDetector(MyPointList& pointList,
                         SrcArrayT* imColor, Real64 sigmaN, Real64 harrisK,
                         Real64 harrisThreshold) 
    {
        Real64 sigmaIntegr = sigmaN;
        Real64 sigmaDiff = 0.7125 * sigmaN;
        
        Array2dScalarReal64* TLx2 = 0;
        Array2dScalarReal64* TLy2 = 0;
        Array2dScalarReal64* TLxLy = 0;
    
        for (int d=1 ; d<=imColor->ElemSize() ; d++) {
            Array2dScalarReal64* imIntensity = 0;
            ProjectRange(imIntensity, imColor, d);
        
            // compute Lx and Ly Gaussian derivatives
            Array2dScalarReal64* Lx = 0;
            Array2dScalarReal64* Ly = 0;
            if(mUseRecGauss)
            {
                RecGauss(Lx, imIntensity, sigmaDiff, sigmaDiff, 1, 0, mPrecision);
                RecGauss(Ly, imIntensity, sigmaDiff, sigmaDiff, 0, 1, mPrecision);
            } else {
                GaussDerivative(Lx, imIntensity, sigmaDiff, 1, 0, mPrecision);
                GaussDerivative(Ly, imIntensity, sigmaDiff, 0, 1, mPrecision);
            }
            
            delete imIntensity;
            
            // VERY IMPORTANT: SCALE MULTIPLIER!
            Array2dScalarReal64* Lx2 = 0;    // aka mu11
            Mul(Lx2, Lx, Lx);
            MulVal(Lx2, Lx2, sigmaDiff * sigmaDiff);
        
            Array2dScalarReal64* Ly2 = 0;    // aka mu22
            Mul(Ly2, Ly, Ly);
            MulVal(Ly2, Ly2, sigmaDiff * sigmaDiff);
        
            Array2dScalarReal64* LxLy = 0;   // aka mu12/mu21
            Mul(LxLy, Lx, Ly);
            MulVal(LxLy, LxLy, sigmaDiff * sigmaDiff);
        
            delete Lx;
            delete Ly;
            
            if(d == 1) {
                TLx2 = Lx2;
                TLy2 = Ly2;
                TLxLy = LxLy;
            } else {
                Add(TLx2,TLx2,Lx2);
                Add(TLy2,TLy2,Ly2);
                Add(TLxLy,TLxLy,LxLy);
                delete Lx2;
                delete Ly2;
                delete LxLy;
            }
        }
        
        // compute Harris  
        Array2dScalarReal64* harrisOutput = 0;
        Array::Harris<Array2dScalarReal64>(harrisOutput, TLx2, TLy2, TLxLy,
                                           sigmaIntegr, harrisK, mUseRecGauss);
        delete TLx2;
        delete TLy2;
        delete TLxLy;
        DoMaxPoints(harrisOutput, std::back_inserter(pointList), harrisThreshold);
            
        // write the Harris cornerness to an image on disk
        /*std::string output = "debug.png";
        Array2dVec3UInt8* buf = ArrayCreate<Array2dVec3UInt8>
            (harrisOutput->CW(), harrisOutput->CH());
        GetRgbPixels(harrisOutput, buf->CPB(), std::string("Stretch"));
    
        WritePng(buf, output);
        delete buf;*/

        delete harrisOutput;
    }
        
    
    template<class SrcArrayT>
    void
    laplacianScaleSelection(PointDescriptorTable* output, SrcArrayT* imIntensity, MyPointList& pointList, Real64 laplaceThreshold) 
    {
        std::vector<Real64> precomputeScales;
        std::vector<Array2dScalarReal64*> cache;
        for(int i = 0; i < 21; i++) {
            Real64 sigma = pow(sqrt(sqrt(Real64(2))),i) * 0.75;
            precomputeScales.push_back(sigma);
            Array2dScalarReal64* tmp = 0;
            Laplacian(tmp, imIntensity, sigma, 2, mUseRecGauss, mPrecision);
            Abs(tmp, tmp);
            cache.push_back(tmp);

            ILOG_DEBUG("Computing Laplacian scale: " << sigma);
        }
        
        typedef std::pair<int, std::pair<int, int> > MyKey;
        std::map<MyKey, int> alreadyHavePoint;

        for(MyPointList::iterator iter = pointList.begin(); iter != pointList.end(); iter++)
        {
            MyPoint* point(*iter);
            Real64 sigmaN = point->mDetectionScale;

            // we must find the index of the nearest sigma in precomputeScales
            int nearestIndex = 0;
            Real64 smallestDifference = 999;
            for(int i = 0; i < precomputeScales.size(); i++) {
                Real64 difference;
                if(mBackwardsCompatible)
                {
                    // make the current (wrong) behavior explicit
                    difference = static_cast<int>(abs(precomputeScales[i] - sigmaN));
                }
                else
                {
                    difference = fabs(precomputeScales[i] - sigmaN);
                }
                if(difference < smallestDifference) 
                {
                    smallestDifference = difference;
                    nearestIndex = i;
                }
            }
            if(smallestDifference > 0.0001) {
                ILOG_WARN("WARNING: scale mismatch");
            }
    
            std::vector<int> interestingScales;
            std::vector<Real64> laplacianValues;
            //std::cout << "LAP: ";
            for(int i = 0; i < precomputeScales.size(); i++) {
                if((i >= nearestIndex - 4) && (i <= nearestIndex + 4)) {
                    interestingScales.push_back(i);
                    // lookup the Laplacian value for this cached scale
                    laplacianValues.push_back(cache[i]->Value(point->mX, point->mY));
                    //std::cout << cache[i]->Value(point.mX, point.mY) << " ";
                }
            }
            //std::cout << std::endl;
    
            /*std::cout << "plot([";
            for(int i = 0; i < precomputeScales.size(); i++) {
                std::cout << cache[i]->Value(point.mX, point.mY) << " ";
            }
            std::cout << "])" << std::endl;*/
            
            // laplacianValues now contains the values of the Laplacian at interestingScales
    
            // select the local maximum
            int nrFound = 0;
            for(int i = 1; i < interestingScales.size() - 1; i++) {
                if((laplacianValues[i] > laplacianValues[i-1]) &&
                   (laplacianValues[i] > laplacianValues[i+1]))
                {
                    if(laplacianValues[i] > laplaceThreshold) {
                        MyKey key(std::make_pair(point->mX, std::make_pair(point->mY, interestingScales[i])));
                        if(alreadyHavePoint.find(key) == alreadyHavePoint.end())
                        {
                            alreadyHavePoint[key] = 0;
                            output->Add(point->mX, point->mY, precomputeScales[interestingScales[i]], 0, point->mCornerness);
                            nrFound++;
                        }
                    }
                }
            }
        }
        for(int i = 0; i < cache.size(); i++) {
            delete cache[i];
        }
    }

    void
    BoostImage(Array2dVec3Real64*& output, Array2dVec3UInt8* input, Real64 w1, Real64 w2, Real64 w3)
    {
        Array2dVec3Real64* opponentImage = 0;
        Rgb2Ooo(opponentImage, input);

        Array2dScalarReal64* channel1 = 0;
        ProjectRange(channel1, opponentImage, 1);
        MulVal(channel1, channel1, w1);

        Array2dScalarReal64* channel2 = 0;
        ProjectRange(channel2, opponentImage, 2);
        MulVal(channel2, channel2, w2);

        Array2dScalarReal64* channel3 = 0;
        ProjectRange(channel3, opponentImage, 3);
        MulVal(channel3, channel3, w3);
        
        MakeFrom3Images(output, channel1, channel2, channel3);
        delete channel1;
        delete channel2;
        delete channel3;
        delete opponentImage;
    }
    

    PointDescriptorTable*
    CHLDetector(Array2dVec3UInt8* input, Real64 harrisThreshold, Real64 harrisK, Real64 laplaceThreshold)
    {
        Timer timer;
        PointDescriptorTable* output = new PointDescriptorTable();
    
        // boost image
        Array2dVec3Real64* imBoosted = 0;
        BoostImage(imBoosted, input, 0.065, 0.850, 0.524);
      
        //MulVal(imIntensity, imIntensity, 1 / 255.0);
        
        // Harris corner detection
        MyPointList centroids;
        Real64 sigmaZero = 0.75;
        Real64 borderIgnoreSize = 4.0;
        for(int sigmaCounter = 1; sigmaCounter <= 8; sigmaCounter++) 
        {
            Real64 sigmaN = pow(sqrt(Real64(2)), sigmaCounter) * sigmaZero;   
            MyPointList pointList;
            harrisCornerDetector(pointList, imBoosted, sigmaN, harrisK, harrisThreshold);
            for(MyPointList::iterator iter = pointList.begin(); iter != pointList.end(); iter++)
            {
                MyPoint* point(*iter);
                
                point->mDetectionScale = sigmaN;
                //ILOG_DEBUG("   " << point->mX << " " << point->mY << " " << point->mCornerness << " " << point->mDetectionScale);
                if((point->mX < borderIgnoreSize) ||
                   (point->mY < borderIgnoreSize) ||
                   (point->mX >= imBoosted->CW() - borderIgnoreSize) ||
                   (point->mY >= imBoosted->CH() - borderIgnoreSize)) 
                {
                    ILOG_DEBUG("   " << point->mX << " " << point->mY << " BORDER CLOSENESS REJECT");
                    continue;
                }
                centroids.push_back(point);
            }
            ILOG_DEBUG(sigmaN << " " << pointList.size());
        }
        ILOG_INFO("Harris corner detector found " << centroids.size() << " points");
        ILOG_DEBUG("[Timing2] HarrisCornerDetector " << timer.SplitTimeStr());
        
        // Laplacian scale selection
        laplacianScaleSelection(output, imBoosted, centroids, laplaceThreshold);
    
        delete imBoosted;
    
        ILOG_INFO("Laplacian scale selection left " << output->Size() << " points");
        ILOG_DEBUG("[Timing2] LaplacianScaleSelection" << timer.SplitTimeStr());
        
        // now sort the points by cornerness
        SortOnColumn5(output, false);
        return output;
    }
    
    
    template<class SrcArrayT>
    PointDescriptorTable*
    HLDetector(SrcArrayT* input, Real64 harrisThreshold, Real64 harrisK, Real64 laplaceThreshold)
    {
        Timer timer;
        PointDescriptorTable* output = new PointDescriptorTable();
    
        // convert to grey-scale
        Array2dScalarReal64* imIntensity = 0;
        RGB2Gray(imIntensity, input);
      
        // scale intensity to the range 0..1 instead of 0..255 so it is comparable to Matlab
        MulVal(imIntensity, imIntensity, 1 / 255.0);
        
        // Harris corner detection
        MyPointList centroids;
        Real64 sigmaZero = 0.75;
        Real64 borderIgnoreSize = 4.0;
        for(int sigmaCounter = 1; sigmaCounter <= 8; sigmaCounter++) 
        {
            Real64 sigmaN = pow(sqrt(Real64(2)), sigmaCounter) * sigmaZero;   
            MyPointList pointList;
            harrisCornerDetector(pointList, imIntensity, sigmaN, harrisK, harrisThreshold);
            for(MyPointList::iterator iter = pointList.begin(); iter != pointList.end(); iter++)
            {
                MyPoint* point(*iter);
                
                point->mDetectionScale = sigmaN;
                //ILOG_DEBUG("   " << point->mX << " " << point->mY << " " << point->mCornerness << " " << point->mDetectionScale);
                if((point->mX < borderIgnoreSize) ||
                   (point->mY < borderIgnoreSize) ||
                   (point->mX >= imIntensity->CW() - borderIgnoreSize) ||
                   (point->mY >= imIntensity->CH() - borderIgnoreSize)) 
                {
                    //ILOG_DEBUG("   " << point->mX << " " << point->mY << " BORDER CLOSENESS REJECT");
                    continue;
                }
                centroids.push_back(point);
            }
            ILOG_DEBUG(sigmaN << " " << pointList.size());
        }
        int info_harrisCount = centroids.size();
        String info_harrisTiming = timer.SplitTimeStr();
        
        // Laplacian scale selection
        laplacianScaleSelection(output, imIntensity, centroids, laplaceThreshold);
    
        delete imIntensity;
    
        ILOG_INFO("Harris found " << info_harrisCount << " (" << info_harrisTiming << "); Laplacian left " << output->Size() << " (" << timer.SplitTimeStr() << ")");
        
        // now sort the points by cornerness
        SortOnColumn5(output, false);
        return output;
    }
    
    ILOG_VAR_DEC;
};

ILOG_VAR_INIT(HarrisLaplaceDetector, Core.Feature);

} // namespace Feature
} // namespace Core
} // namespace Impala

#endif

---