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

---

InterestPointFeature.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_InterestPointFeature_h
#define Impala_Core_Feature_InterestPointFeature_h

#include "Core/Array/MakeRoi.h"
#include "Core/Geometry/Rectangle.h"
#include "Core/Feature/FeatureDefinition.h"
#include "Core/Feature/FeatureTable.h"
#include "Core/Feature/Surf.h"
#include "Core/Feature/RandomForest.h"
#include "Basis/Timer.h"

#include "Core/Matrix/MatNorm2Dist.h"
#include "Core/Matrix/VectorQuantize.h"
#include "Util/RandomGenerator.h"

#include "Core/Geometry/InterestPointList.h"
#include "Core/Geometry/InterestCircle.h"
#include "Core/Geometry/InterestRectangleZ.h"
#include "Core/Feature/RegionDescriptor.h"
#include "Core/Feature/KoenFIST2.h"
#include "Core/Feature/HarrisLaplaceDetector.h"
#include "Core/Feature/LaplacianDetector.h"
#include "Core/Geometry/InterestPointSelector.h"
#include "Core/Geometry/MaskedInterestPointSelector.h"

#include "Core/Feature/PointDescriptorTable.h"

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

namespace Impala
{
namespace Core
{
namespace Feature
{


class InterestPointFeature
{
public:
    typedef Array::Array2dScalarReal64 Array2dScalarReal64;
    typedef Matrix::Mat Mat;
    typedef Geometry::InterestPointSelector InterestPointSelector;
    typedef Geometry::InterestPointList InterestPointList;
    typedef Geometry::InterestPoint InterestPoint;

    // command line parameters
    int         mVerbose;
    int         mStripBorder;
    bool        mUseRecGauss;
    String      mName;         // for display/filename purposes
    FeatureDefinition mFeatureDefinition;
    int         mKeepLimited;
    String      mOutputFormat;
    FeatureTable*     mCodebook;
    Mat*              mDescriptorReduce;
    String            mDescriptorReduceFilename;
    bool              mClusterInput;
    FeatureTable*     mClusterInputData;
    String            mClusterInputImage;
    std::vector<Vector::VectorTem<Real64>*> mLastCodebookVectors;//used in InterestPointProc
    std::vector<InterestPointSelector*> mPointSelector;
    bool                         mDynamicPointSelector;
    String                       mDynamicPointSelectorConfig;
    Impala::Util::Database* mDatabase;
    std::vector<Real64> mDenseSamplingScales;
    int                 mSurfParams[3];

    InterestPointFeature(CmdOptions& options) :
        mVerbose(0),
        mStripBorder(0),
        mUseRecGauss(false),
        mKeepLimited(0),
        mCodebook(0),
        mDescriptorReduce(0),
        mClusterInput(false),
        mClusterInputData(0),
        mDynamicPointSelector(false),
        mDatabase(0)
    {
        mVerbose     = options.GetInt("verb");
        mUseRecGauss = options.GetBool("recGauss");
        mStripBorder = options.GetInt("borderWidth");
        mOutputFormat = options.GetString("outputFormat");

        mClusterInput = options.GetBool("clusterInput");
        mKeepLimited = options.GetInt("keepLimited");
        if (!options.GetBool("noRandomizeSeed"))
            Util::RandomGenerator::GetInstance().RandomizeSeed();

        String detector    = options.GetString("detector");
        mName              = detector + "-" + options.GetString("descriptor");
        mFeatureDefinition = FeatureDefinition(mName);

        String descriptor = options.GetString("descriptor");
        if((descriptor.size() >= 4) &&
           (descriptor.substr(descriptor.size()-4,4) == "surf"))
        {
            Core::Feature::GetDSurfOptions(options, mSurfParams[0],
                                           mSurfParams[1], mSurfParams[2]);
            mFeatureDefinition.AddParameter
                ("sc", MakeString(mSurfParams[0]) + "-" +
                 MakeString(mSurfParams[1]) + "-" + MakeString(mSurfParams[2]));
        }
        if((detector == "harrislaplace") || (detector == "harrislaplace2"))
        {
            if(options.GetString("harrisThreshold") != "1e-9")
            {
                mFeatureDefinition.AddParameter
                    ("ht", options.GetString("harrisThreshold"));
            }
        }
        if((detector == "densesampling") || (detector == "denseall"))
        {
            // detector config
            int spacing = options.GetInt("ds_spacing");
            String scales = options.GetString("ds_scales");

            Util::StringParser sp(scales);
            while(!sp.TheEnd())
            {
                Real64 scale = sp.GetDouble('+', false, true);
                mDenseSamplingScales.push_back(scale);
            }
            if(mDenseSamplingScales.size() == 0)
            {
                // fallback to estimate
                Real64 scale = Real64(spacing) / 5;
                // Most computation code uses 4*scale for the circular region size.
                // Therefore, it makes sense to not provide the 'actual' region of
                // interest here, but something smaller.
                mDenseSamplingScales.push_back(scale);
            }
            String name = MakeString(spacing);
            for(int i = 0; i < mDenseSamplingScales.size(); i++)
            {
                name += "-" + MakeString(mDenseSamplingScales[i]);
            }
            mFeatureDefinition.AddParameter("ds", name);
        }

        mDescriptorReduceFilename = options.GetString("descriptorReduce");
        if(!mDescriptorReduceFilename.empty()) mFeatureDefinition.AddParameter("red", "1");

        if(!options.GetString("descriptorOpponentWeight").empty())
        {
            mFeatureDefinition.AddParameter
                ("ow", options.GetString("descriptorOpponentWeight"));
        }
        if(StringStartsWith(options.GetString("pointSelector"), "dynamic-"))
        {
            // keep mPointSelector empty: we load it on-demand
            mDynamicPointSelector = true;
            String config = options.GetString("pointSelector");
            mFeatureDefinition.AddParameter("ps", config);
            config = StringReplace(config, "dynamic-", "");
            mDynamicPointSelectorConfig = config;
        }
        else if(options.GetString("pointSelector").size() != 0)
        {
            // filter on location (spatial pyramid)
            String config = options.GetString("pointSelector");
            String configName = config;
            if(config == "pyramid-1x1-2x2")
            {
                config = "P1x1#0+P2x2#0+P2x2#1+P2x2#2+P2x2#3";
            }
            if(config == "pyramid-1x1-2x2-3x1")
            {   // this is vertical (!)
                config = "P1x1#0+P2x2#0+P2x2#1+P2x2#2+P2x2#3+P3x1#0+P3x1#1+P3x1#2";
            }
            if(config == "pyramid-1x1-2x2-1x3")
            {
                // this is division into three horizontal bars
                config = "P1x1#0+P2x2#0+P2x2#1+P2x2#2+P2x2#3+P1x3#0+P1x3#1+P1x3#2";
            }
            Util::StringParser sp(config);
            while(!sp.TheEnd())
            {
                String selectorConfig = sp.GetString('+', false);
                if(selectorConfig.size() == 0) break;
                InterestPointSelector* selector =
                    new InterestPointSelector(selectorConfig, mStripBorder);
                mPointSelector.push_back(selector);
            }
            configName = StringReplaceAll(configName, "+", "");
            configName = StringReplaceAll(configName, "#", "x");
            mFeatureDefinition.AddParameter("ps", configName);
        }
        else
        {
            // default: a filter spanning the complete image
            InterestPointSelector* selector =
                new InterestPointSelector("P1x1#0", mStripBorder);
            mPointSelector.push_back(selector);
        }
    }

    ~InterestPointFeature()
    {
        if(mCodebook)
        {
            delete mCodebook;
            mCodebook = 0;
        }
    }

    static void AddCmdOptions()
    {
        CmdOptions& options = CmdOptions::GetInstance();
        options.AddOption(0, "borderWidth",      "size",        "0");
        options.AddOption(0, "recGauss",         "",            "0");
        options.AddOption(0, "loadRegions",      "filename",    "");
        options.AddOption(0, "saveRegions",      "filename",    "");
        options.AddOption(0, "descriptor",       "mode",        "");
        options.AddOption(0, "output",           "filename",    "");
        options.AddOption(0, "outputFormat",     "mode",    "");
        options.AddOption(0, "descriptorBinCount","bins",        "15");
        options.AddOption(0, "descriptorPreSmoothing","",        "0");
        options.AddOption(0, "descriptorOpponentWeight","weight","");

        /* clustering support (i.e. when not in codebook mode) */
        options.AddOption(0, "clusterInput",     "",            "0");
        options.AddOption(0, "keepLimited",      "count",       "-1");
        options.AddOption(0, "noRandomizeSeed",  "",            "0");

        /* codebook support */
        options.AddOption(0, "descriptorReduce", "filename",    "");
        options.AddOption(0, "codebook",         "filename",    "");
        options.AddOption(0, "codebookMode",     "name",        "hard");
        options.AddOption(0, "codebookSigma",    "sigma",       "");

        /* Lazebnik spatial pyramid support */
        options.AddOption(0, "pointSelector",    "config",    "");

        /*** detectors ***/
        options.AddOption(0, "detector",         "mode",        "harrislaplace");

        /* for Harris-laplace detectors */
        options.AddOption(0, "harrisThreshold",  "threshold",   "1e-9");
        options.AddOption(0, "harrisK",          "k",           "0.06");
        options.AddOption(0, "laplaceThreshold", "threshold",   "0.03");

        /* for the dense sampling detector */
        options.AddOption(0, "ds_spacing",  "pixels",                  "6");
        options.AddOption(0, "ds_scales",   "scale1+scale2+...",       "");

        Core::Feature::AddDSurfOptions(options);
    }


    String
    GetFeatureName()
    {
        return mFeatureDefinition.AsString();
    }


    String
    GetFeatureNameShort()
    {
        return mFeatureDefinition.GetName();
    }


    void
    SetCodebook(String filename, Impala::Util::Database* db)
    {
        String codebook;
        Util::IOBuffer* buf = 0;
        if(db)
        {
            /* translate this name into a complete path */
            codebook = db->GetFilePath(filename, false, false);
            ILOG_DEBUG("SetCodebook using DB; " << filename << "; serverName="
                       << codebook);
            /* read it */
            if(codebook.empty())
            {
                ILOG_ERROR("Cannot open codebook: " << filename);
                exit(1);
            }
            buf = db->GetIOBuffer(codebook, true, false, "", 0, true);
            mDatabase = db;
            if(!mDescriptorReduceFilename.empty())
            {
                ReadRaw(mDescriptorReduce, db->GetFilePath(mDescriptorReduceFilename, false, false), db);
                mFeatureDefinition.AddParameter("red", MakeString(Matrix::MatNrCol(mDescriptorReduce)));
            }
        }
        else
        {
            ILOG_DEBUG("Command-line (colorDescriptor binary)");
            codebook = filename;
            buf = new Util::IOBufferFile(filename, true, true);
            if(!mDescriptorReduceFilename.empty())
            {
                Util::IOBuffer* buf2 = new Util::IOBufferFile(mDescriptorReduceFilename, true, true);
                ReadRaw(mDescriptorReduce, buf2);
                mFeatureDefinition.AddParameter("red", MakeString(Matrix::MatNrCol(mDescriptorReduce)));
                delete buf2;
            }
        }
        
        ILOG_DEBUG("SetCodebook after IOBuffer " << buf);
        if (buf)
        {
            buf->Seek(0, SEEK_SET);
            char buf2[2049];
            buf2[2048] = '\0';
            buf->Read(buf2, 2048);
            buf->Seek(0, SEEK_SET);
            String recognizer(buf2);
            if(recognizer.find("nr strings") != String::npos)
            {   // nr strings : then FeatureTableOld
                ILOG_ERROR("Your codebook is in FeatureTableOld format!");
            }
            else
            {
                mCodebook = new FeatureTable(filename);
                Read(mCodebook, buf);
                ILOG_DEBUG("FeatureTable style: " << filename);
            }
            delete buf;
        }
        else
        {
            ILOG_ERROR("SetCodebook: IOBuffer is 0");
            return;
        }
        
        ILOG_INFO("Reading codebook: " << codebook <<
                  "; dimensionCount = " << mCodebook->GetFeatureVectorLength()
                  << "; size = " << mCodebook->Size());

        CmdOptions& options = CmdOptions::GetInstance();
        if(options.GetString("codebookMode") != "hard")
        {
            if(options.GetString("codebookMode") == "hardindex")
            {
                mFeatureDefinition.AddParameter
                    ("cm", options.GetString("codebookMode") + "-" +
                     FileNameBase(codebook));
            }
            else
            {
                mFeatureDefinition.AddParameter
                    ("cm", options.GetString("codebookMode"));
            }
        }
        if(options.GetString("codebookSigma") != "")
        {
            mFeatureDefinition.AddParameter
                ("cs", options.GetString("codebookSigma"));
        }
    }

    bool
    FindInterestPoints(CmdOptions& options, Array::Array2dVec3UInt8* input,
                       Quid quid,     // identifies input image data
                       String outputFilename, bool borderAlreadyRemoved=false,
                       String loadRegionsFilename="")
    {
        ILOG_DEBUG("start FindInterestPoints");
        if(loadRegionsFilename == "")
        {
            // try the command-line argument if it is not already specified
            loadRegionsFilename = options.GetString("loadRegions");
        }
        String saveRegionsFilename = options.GetString("saveRegions");
        String detectorMode = options.GetString("detector");
        String descriptorMode = options.GetString("descriptor");
        Timer timer;
        Array::Array2dVec3UInt8* inputNoBorder = input;
        if(!borderAlreadyRemoved)
        {
            // remove the border, if requested
            Geometry::Rectangle rect(mStripBorder, mStripBorder,
                                     input->CW() - mStripBorder - 1,
                                     input->CH() - mStripBorder - 1, true);
            inputNoBorder = MakeRoi(input, rect);
        }

        PointDescriptorTable* pointData = 0;

        /* Apply detector */
        ILOG_DEBUG("applying detector");
        if(!loadRegionsFilename.empty())
        {
            // load the points from file
            pointData = PointDescriptorTable::ImportFromFile(detectorMode, 
                                          loadRegionsFilename, mStripBorder);
        }

        if(loadRegionsFilename.empty())
        {
            pointData = ApplyDetector(inputNoBorder, detectorMode);
        }
        if(descriptorMode.empty())
        {
            // dump the interest regions to the output stream if we are not
            // going to compute descriptors later
            pointData->ExportToFile(outputFilename, mStripBorder, 
                                    mOutputFormat);
        }
        if(!saveRegionsFilename.empty())
        {
            ILOG_INFO("Writing to: " << saveRegionsFilename);
            pointData->ExportToFile(saveRegionsFilename, mStripBorder,
                                    mOutputFormat);
        }

        ILOG_DEBUG("[Timing] Region loading / applying detector finished in "
                   << timer.SplitTimeStr());

        /* descriptor computation */
        ILOG_DEBUG("computing descriptor");
        if(descriptorMode.size() != 0)
        {
            ComputeDescriptors(pointData, inputNoBorder, descriptorMode);

            ILOG_DEBUG("[Timing] Descriptors finished at " << timer.SplitTimeStr());
            if(mCodebook)
            {
                ILOG_DEBUG("preparing for codebook");
                // clean up previous codebook vectors (if any)
                for(int i = 0; i < mLastCodebookVectors.size(); i++)
                {
                    delete mLastCodebookVectors[i];
                }
                mLastCodebookVectors.clear();
                /* project it onto a codebook */
                Timer timerProjection;
                pointData = ProjectOntoCodebook(pointData,
                                    options.GetString("codebookMode"),
                                    inputNoBorder->CW(), inputNoBorder->CH(),
                                    quid);
                ILOG_DEBUG("Just codebook projection took: " <<
                           timerProjection.SplitTimeStr());
                if(outputFilename.size() != 0)
                {
                    pointData->ExportToFile(outputFilename, mStripBorder, 
                                            mOutputFormat);
                }
            }
            else
            {
                // do reduction
                if((mKeepLimited > 0) && (pointData->Size() > mKeepLimited))
                {
                    ILOG_DEBUG("reducing point count");
                    // we might have to limit the number of points
                    pointData = ReducePointCount(pointData, mKeepLimited);
                }

                if(mClusterInput)
                {
                    ILOG_DEBUG("doing cluster input");
                    InterestPointList outputPointList;
                    pointData->ToPointList(outputPointList);

                    for(InterestPointList::iterator iter = outputPointList.begin();
                        iter != outputPointList.end(); iter++)
                    {
                        InterestPoint* point = *iter;
                        int n = point->mDescriptor.size();
                        Vector::VectorTem<Real64> descriptor(n);
                        for(int i = 0; i < n; i++)
                        {
                            descriptor.Elem(i) = point->mDescriptor[i];
                        }

                        if(!mClusterInputData)
                        {
                            ILOG_DEBUG("initializing mClusterInputData");
                            // initialize the table
                            FeatureDefinition def = GetFeatureName();
                            def.AddParameter("clusterinput", "2");
                            mClusterInputData = new FeatureTable(def, 500, n);
                        }

                        // string representation of point:
                        //point->Serialize(borderOffset+1, borderOffset+1)
                        //mClusterInputData->Add(mClusterInputImage + " " +
                        //point->Serialize(), descriptor); TODO: create proper
                        //CLUSTERINPUT quid class
                        mClusterInputData->Add(mClusterInputData->Size() + 1,
                                               descriptor);
                    }
                    ILOG_DEBUG("done with cluster input");
                    outputPointList.EraseAndDelete();
                }
                else
                {
                    ILOG_DEBUG("writing point list");
                    if(outputFilename.empty())
                    {
                        ILOG_WARN("Warning: no output file to write to. Did you forget to specify --output?");
                    }
                    else
                    {
                        pointData->ExportToFile(outputFilename, mStripBorder, 
                                                mOutputFormat);
                    }
                }
            }
            if(!borderAlreadyRemoved)
            {
                delete inputNoBorder;
            }
        }

        ILOG_DEBUG("Completed " << GetFeatureName() << " in " <<
                   timer.SplitTimeStr());
        delete pointData;

        return true;
    }
private:
    void
    AddToPyramids(int codeword, double x, double y,
                  int imageWidth, int imageHeight, std::vector<int>& regionCounts)
    {
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            InterestPointSelector* selector = mPointSelector[ps];
            if(selector->Accept(x, y, imageWidth, imageHeight))
            {
                mLastCodebookVectors[ps]->Elem(codeword) =
                    mLastCodebookVectors[ps]->Elem(codeword) + 1;
                regionCounts[ps] = regionCounts[ps] + 1;
            }
        }
    }
    
    void
    ProjectUNC(Mat* matrixDescriptors, int imageWidth,
               int imageHeight, std::vector<int>& X, std::vector<int>& Y,
               std::vector<int>& regionCounts)
    {
        using namespace Matrix;
        Timer timerProjection;
        if(mDescriptorReduce)
        {
            matrixDescriptors = MatMul(matrixDescriptors, mDescriptorReduce);
        }
        Mat* matrixCodebook = mCodebook->GetColumn2()->GetStorage();
        ILOG_DEBUG(MatNrRow(matrixDescriptors) << " " <<
                   MatNrCol(matrixDescriptors));
        ILOG_DEBUG(MatNrRow(matrixCodebook) << " " <<
                   MatNrCol(matrixCodebook));
        ILOG_DEBUG("Codebook projection (init) took: " <<
                   timerProjection.SplitTimeStr());

        int step = 10000;
        for(int k = 0; k < MatNrRow(matrixDescriptors); k += step)
        {
            int s = MatNrRow(matrixDescriptors) - k;
            if(step < s) s = step;
    
            Mat* partial = MatCreate<Mat>(s, MatNrCol(matrixDescriptors),
                                          MatE(matrixDescriptors, k, 0), true);

            // result will be a (#partial-descriptors, codebook size) matrix
            Mat* distances = MatNorm2DistTransposed(partial,
                                                    matrixCodebook);
            delete partial;
            
            ILOG_DEBUG("Codebook projection (matrix) took: " <<
                       timerProjection.SplitTimeStr());
            CmdOptions& options = CmdOptions::GetInstance();
            Real64 sigma = options.GetDouble("codebookSigma");

            // pull a Gaussian kernel over the Euclidian distances
            for(int vj = 0; vj < s; vj++)
            {
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    Real64 value = *MatE(distances, vj, i);
                    *MatE(distances, vj, i) = exp((-0.5 * value * value) /
                                                  (sigma * sigma));
                }
            }

            // sum over axis 0 for the normalization factor
            // expected size is (#descriptors, 1)
            Mat* normalizationFactor = MatSumAxis1(distances);
            ILOG_DEBUG("normalizationFactor: " << MatNrRow(normalizationFactor)
                       << " " << MatNrCol(normalizationFactor));

            // perform normalization
            for(int vj = 0; vj < s; vj++)
            {
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    *MatE(distances, vj, i) = *MatE(distances, vj, i) /
                        *MatE(normalizationFactor, vj, 0);
                }
            }
            delete normalizationFactor;

            // remains: sum over the *right* descriptors
            for(int vj = 0; vj < s; vj++)
            {
                for(int ps = 0; ps < mPointSelector.size(); ps++)
                {
                    if(mPointSelector[ps]->Accept(X[k+vj], Y[k+vj],
                                                  imageWidth, imageHeight))
                    {
                        for(int i = 0; i < mCodebook->Size(); i++)
                        {
                            mLastCodebookVectors[ps]->Elem(i) =
                                mLastCodebookVectors[ps]->Elem(i) +
                                *MatE(distances, vj, i);
                        }
                        regionCounts[ps] = regionCounts[ps] + 1;
                    }
                }
            }
            delete distances;
            ILOG_DEBUG("Codebook projection (piece) took: " <<
                       timerProjection.SplitTimeStr() << " " << mPointSelector.size());
        }
        ILOG_INFO("Codebook projection (UNC/" << MatNrRow(matrixDescriptors) <<
                  ") took: " << timerProjection.SplitTimeStr());
        if(mDescriptorReduce)
        {
            delete matrixDescriptors;
        }
    }

    void
    ProjectHardMatrix(Mat* matrixDescriptors,
                      int imageWidth, int imageHeight, std::vector<int>& X,
                      std::vector<int>& Y, std::vector<int>& regionCounts)
    {
        using namespace Matrix;
        Timer timerProjection;
        if(mDescriptorReduce)
        {
            matrixDescriptors = MatMul(matrixDescriptors, mDescriptorReduce);
        }
        Mat* matrixCodebook = mCodebook->GetColumn2()->GetStorage();
        ILOG_DEBUG(MatNrRow(matrixDescriptors) << " " <<
                   MatNrCol(matrixDescriptors));
        ILOG_DEBUG(MatNrRow(matrixCodebook) << " " <<
                   MatNrCol(matrixCodebook));
        ILOG_DEBUG("Codebook projection (init) took: " <<
                   timerProjection.SplitTimeStr());
        // result will be a (#descriptors, 2) matrix
        Mat* assignment =
            Core::Matrix::VectorQuantize(matrixDescriptors, matrixCodebook);
            
        ILOG_DEBUG("Codebook projection (matrix) took: " <<
                   timerProjection.SplitTimeStr());
        for(int vj = 0; vj < MatNrRow(matrixDescriptors); vj++)
        {
            Real64 bestDistance = *MatE(assignment, vj, 1);
            int    bestIndex = static_cast<int>(*MatE(assignment, vj, 0));
            AddToPyramids(bestIndex, X[vj], Y[vj], imageWidth, imageHeight,
                          regionCounts);
            //std::cout << "Codebook projection (inner loop) took: " <<
            //timerProjection.SplitTimeStr() << std::endl;
        }
        ILOG_DEBUG("Codebook projection (matrix-form) took: " <<
                   timerProjection.SplitTimeStr());
        delete assignment;
        if(mDescriptorReduce)
        {
            delete matrixDescriptors;
        }
    }


    void
    ProjectHardIndex(Mat* matrixDescriptors, 
                     PointDescriptorTable* pointData)
    {
        using namespace Matrix;
        Timer timerProjection;
        if(mDescriptorReduce)
        {
            matrixDescriptors = MatMul(matrixDescriptors, mDescriptorReduce);
        }
        Mat* matrixCodebook = mCodebook->GetColumn2()->GetStorage();
        ILOG_DEBUG("Codebook projection (init) took: " <<
                   timerProjection.SplitTimeStr());
        // result will be a (#descriptors, 2) matrix
        Mat* assignment =
            Core::Matrix::VectorQuantize(matrixDescriptors, matrixCodebook);
        ILOG_INFO("Codebook projection (matrix) took: " <<
                  timerProjection.SplitTimeStr());

        Array::Array2dScalarReal64* newDesc = new Array::Array2dScalarReal64(1, 
                                                   pointData->Size(), 0, 0);
        for(int vj = 0; vj < MatNrRow(matrixDescriptors); vj++)
        {
            Real64 bestDistance = *MatE(assignment, vj, 1);
            int    bestIndex = static_cast<int>(*MatE(assignment, vj, 0));
            *MatE(newDesc, vj, 0) = bestIndex;
        }
        // this call also takes ownership
        pointData->ReplaceAllDescriptors(newDesc);
        delete assignment;
        for(int i = 0; i < mLastCodebookVectors.size(); i++)
            delete mLastCodebookVectors[i];
        mLastCodebookVectors.clear();
        if(mDescriptorReduce)
        {
            delete matrixDescriptors;
        }
    }
    
    
    void
    ProjectForest(Mat* descriptors, int imageWidth,
                  int imageHeight, std::vector<int>& X, std::vector<int>& Y,
                  std::vector<int>& regionCounts)
    {
        using namespace Matrix;
        // make trees from the table
        ILOG_DEBUG("creating random forest from feature table");
        RandomForest forest = ReadRandomForest(mCodebook);
        ILOG_DEBUG("codewords: "<<GetCodebookLength(mCodebook)<<", feattable size: "<<
                  mLastCodebookVectors[0]->Size() <<", desc count"<<
                  MatNrRow(descriptors) <<", desc length: "<< MatNrCol(descriptors));
        Timer timerProjection;
        for(int vj = 0; vj < MatNrRow(descriptors); vj++)
        {
            for(int i=0; i<forest.size(); ++i)
            {
                RandomTree* tree = forest[i];
                ILOG_DEBUG("getting codeword");
                int codeword = tree->GetCodeWord(Vector::VectorTem<Real64>(MatNrCol(descriptors), descriptors->CPB(0, vj), true));
                ILOG_DEBUG("got codeword");
                AddToPyramids(codeword, X[vj], Y[vj], imageWidth, imageHeight,
                              regionCounts);
            }
        }
        DeleteForest(forest);
        ILOG_DEBUG("Codebook projection (forest) took: " <<
                   timerProjection.SplitTimeStr());
    }
    
    PointDescriptorTable*
    ProjectOntoCodebook(PointDescriptorTable* pointData, String codebookMode,
                        int imageWidth, int imageHeight, Quid q)
    {
        using namespace Matrix;
        ILOG_DEBUG("ProjectOntoCodebook called");
        PointDescriptorTable* output = new PointDescriptorTable(
                                       pointData->GetFeatureDefinition());
        Array::Array2dScalarInt32* tempDynamic = 0;
        if(mDynamicPointSelector)
        {
            // set the point selector configuration based on the image
            String filename = "PointSelector/" + mDynamicPointSelectorConfig +
                              "/" + MakeString(q) + ".raw";
            Util::IOBuffer* buf = 0;
            if(mDatabase)
            {
                /* translate this name into a complete path */
                filename = mDatabase->GetFilePath(filename, false, false);
                /* read it */
                if(filename.empty())
                {
                    ILOG_ERROR("Cannot open dynamic point selector: " << filename);
                }
                buf = mDatabase->GetIOBuffer(filename, true, false, "", 0, true);
            }
            else
            {
                buf = new Util::IOBufferFile(filename, true, true);
            }

            // load the selectors from file
            ReadRaw(tempDynamic, buf);
            delete buf;
            for(int i = 0; i < tempDynamic->CH(); i++)
            {
                *tempDynamic->CPB(0, i) -= (mStripBorder + 1);
                *tempDynamic->CPB(1, i) -= (mStripBorder + 1);
            }
            
            for(int ps = 0; ps < tempDynamic->CW() - 2; ps++)
            {
                mPointSelector.push_back
                    (new Geometry::MaskedInterestPointSelector(tempDynamic, ps));
            }
        }

        std::vector<int> regionCounts;
        int codebookLength = mCodebook->Size();
        if(codebookMode == "forest")
        {
            codebookLength = GetCodebookLength(mCodebook);
            ILOG_DEBUG("forest codebook length = "<< codebookLength);
        }
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            Vector::VectorTem<Real64>* codebookVector = new Vector::VectorTem<Real64>(codebookLength);
            for(int i = 0; i < codebookLength; i++)
            {
                codebookVector->Elem(i) = 0.0;
            }
            mLastCodebookVectors.push_back(codebookVector);
            regionCounts.push_back(0);
        }
        ILOG_DEBUG("done making codebook vector");
        if(pointData->Size() == 0)
        {
            delete pointData;
            return output;
        }
        int descriptorLength = pointData->GetDescriptorLength();
        std::vector<int> X;
        std::vector<int> Y;
        X.reserve(pointData->Size());
        Y.reserve(pointData->Size());

        ILOG_DEBUG("extracting descriptor data");
        for(int it = 0; it < pointData->Size(); it++)
        {
            X.push_back(static_cast<int>(pointData->GetX(it)));
            Y.push_back(static_cast<int>(pointData->GetY(it)));
        }

        // dimensionality checks
        if((codebookMode != "forest") &&
           (mCodebook->Get2(0).Size() != descriptorLength))
        {
            if(mDescriptorReduce)
            {
                if((mCodebook->Get2(0).Size() != MatNrCol(mDescriptorReduce))
                   || (descriptorLength != MatNrRow(mDescriptorReduce)))
                {
                    ILOG_ERROR("Reduced descriptor dimensionality mismatch");
                    throw "bye";
                }
            }
            else
            {
                ILOG_ERROR("ProjectOntoCodebook found visual word"<<
                       " and descriptor of different sizes!\n\t\t"<<
                       "codebook: "<< mCodebook->Get2(0).Size() <<
                       ", descriptor:"<< descriptorLength);
                throw "Dimensionality mismatch!";
            }
        }

        // create a wrapper with the correct dimensions
        Array2dScalarReal64* descriptors = 
            new Array2dScalarReal64(descriptorLength, pointData->Size(), 0, 0,
            pointData->GetColumn6()->GetStorage()->mData, true);
            
        ILOG_DEBUG("extracted descriptor data");

        // Visual word uncertainty (UNC)
        if(codebookMode == "unc")
            ProjectUNC(descriptors, imageWidth, imageHeight, 
                       X, Y, regionCounts);
        // hard assignment, matrix style
        if((codebookMode == "hardmatrix") || (codebookMode == "hard"))
            ProjectHardMatrix(descriptors, imageWidth, imageHeight,
                              X, Y, regionCounts);
        // hard assignment, matrix style, and only output the indices (special)
        if(codebookMode == "hardindex")
        {
            ProjectHardIndex(descriptors, pointData);
            delete descriptors;
            delete output;
            return pointData;
        }
        if(codebookMode == "forest")
            ProjectForest(descriptors, imageWidth, imageHeight,
                          X, Y, regionCounts);
        delete descriptors;

        // normalize so it sums to one
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            // debugStream << "\n\npyramid part #"<<ps<<"\n";
            int regionCount = regionCounts[ps];
            ILOG_DEBUG("regCount " << regionCount);
            if(regionCount == 0)
                regionCount = 1;
            for(int i = 0; i < codebookLength; i++)
            {
                mLastCodebookVectors[ps]->Elem(i) =
                    mLastCodebookVectors[ps]->Elem(i) / regionCount;
            }
        }
        if(mPointSelector.size() > 0)
            output->SetDescriptorLength(codebookLength);
        output->Reserve(mPointSelector.size(), false);
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            output->Add(0, 0, 0, 0, 0);
            //for(int i = 0; i < codebookLength; i++)
            //    std::cout << mLastCodebookVectors[ps]->Elem(i) << " ";
            output->SetDescriptor(ps, *(mLastCodebookVectors[ps]));
        }

        if(mDynamicPointSelector)
        {
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                delete mPointSelector[ps];
            }
            mPointSelector.clear();
            delete tempDynamic;
        }
        delete pointData;
        return output;
    }


    PointDescriptorTable*
    ReducePointCount(PointDescriptorTable* pointData, int keepCount)
    {
        std::vector<int> indices;
        bool* keepers = new bool[pointData->Size()];
        for(int j = 0; j < pointData->Size(); j++)
        {
            indices.push_back(j);
            keepers[j] = false;
        }
        // do not know how to initialize this random number generator
        //std::random_shuffle(indices.begin(), indices.end());
        Util::RandomGenerator::GetInstance().PermutateVector(indices);
        
        for(int j = 0; j < keepCount; j++)
        {
            keepers[indices[j]] = true;
        }
        
        PointDescriptorTable* output = new PointDescriptorTable(pointData->GetFeatureDefinition());
        output->SetDescriptorLength(pointData->GetDescriptorLength());
        Select(output, pointData, keepers, false);
        delete pointData;
        return output;
    }


    template<class SrcArrayT>
    PointDescriptorTable*
    DenseAllDetector(SrcArrayT* input)
    {
        PointDescriptorTable* output = new PointDescriptorTable();
        int w = input->CW();
        int h = input->CH();
        output->SetEmpty();
        output->Reserve(w * h * mDenseSamplingScales.size(), false);
        for(int y = 0; y < h; y += 1)
        {
            for(int x = 0; x < w; x++)
            {
                for(int q = 0; q < mDenseSamplingScales.size(); q++)
                {
                    output->Add(x, y, mDenseSamplingScales[q], 0, 0);
                }
            }
        }
        ILOG_INFO("Points detected: " << output->Size());
        return output;
    }


    template<class SrcArrayT>
    PointDescriptorTable*
    DenseSamplingDetector(SrcArrayT* input, int spacing)
    {
        PointDescriptorTable* output = new PointDescriptorTable();
        int halfspacing = spacing / 2;
        int w = input->CW();
        int h = input->CH();
        int i = 0;
        for(int y = spacing; y < h - spacing; y += spacing)
        {
            int x = spacing;
            if(i % 2 == 0)
            {
                x += halfspacing;
            }
            for(; x < w - spacing; x += spacing)
            {
                for(int q = 0; q < mDenseSamplingScales.size(); q++)
                {
                    //Real64 scale = sqrt(Real64(spacing * spacing + spacing *
                    //spacing));
                    //scale = scale / 4;
                    output->Add(x, y, mDenseSamplingScales[q], 0, 0);
                }
            }
            i += 1;
        }
        return output;
    }
    
    
    PointDescriptorTable*
    ApplyDetector(Array::Array2dVec3UInt8*  inputNoBorder, String detectorMode)
    {
        CmdOptions& options = CmdOptions::GetInstance();
        Real64 harrisThreshold  = options.GetDouble("harrisThreshold");
        Real64 harrisK          = options.GetDouble("harrisK");
        Real64 laplaceThreshold = options.GetDouble("laplaceThreshold");
        
        PointDescriptorTable* output = 0;

        if((detectorMode == "harrislaplace") || (detectorMode == "harrislaplace2"))
        {
            HarrisLaplaceDetector detector(mUseRecGauss);
            detector.mBackwardsCompatible = false;
            output = detector.HLDetector(inputNoBorder, harrisThreshold,
                                harrisK, laplaceThreshold);
        }
        else if(detectorMode == "colorharrislaplace")
        {
            // recommendation: set harrisThreshold to at least 0.001 in color
            // mode (due to very different ranges!)
            HarrisLaplaceDetector detector(mUseRecGauss);
            detector.mBackwardsCompatible = false;
            output = detector.CHLDetector(inputNoBorder, harrisThreshold,
                                 harrisK, laplaceThreshold);
        }
        else if(detectorMode == "laplacian")
        {
            output = LaplacianDetector(inputNoBorder, true);
        }
        else if(detectorMode == "densesampling")
        {
            int spacing = options.GetInt("ds_spacing");
            output = DenseSamplingDetector(inputNoBorder, spacing);
        }
        else if(detectorMode == "denseall")
        {
            output = DenseAllDetector(inputNoBorder);
        }
        else if(detectorMode == "none")
        {
            // we don't add point in the detector phase, they are added in the
            // descsriptor fase. see Core/Feature/Surf.h for an example
            output = new PointDescriptorTable();
        }
        else
        {
            ILOG_ERROR("[ERROR] Unknown detector mode: " << detectorMode);
            return 0;
        }
        output->SetFeatureDefinition(FeatureDefinition(detectorMode));
        return output;
    }
    
    
    void
    ComputeDescriptors(PointDescriptorTable* pointData,
                       Array::Array2dVec3UInt8*  inputNoBorder,
                       String descriptorMode)
    {
        if((descriptorMode == "huesift") || (descriptorMode == "huefist"))
        {
            CalculateFISTDescriptors(inputNoBorder, pointData, "sift");
            InterestPointList outputPointList;
            pointData->ToPointList(outputPointList);
            CalculateRegionDescriptors(inputNoBorder, outputPointList,
                                       "huehistogram");

            // now fix the range of hue...
            for(InterestPointList::iterator iter = outputPointList.begin();
                iter != outputPointList.end(); iter++)
            {
                InterestPoint* point = *iter;
                for(int i = 128; i < point->mDescriptor.size(); i++)
                {
                    point->mDescriptor[i] =
                        static_cast<int>(point->mDescriptor[i] * 512);
                }
            }
            pointData->SetEmpty();
            pointData->FromPointList(outputPointList);
            outputPointList.EraseAndDelete();
       }
        else if((descriptorMode.size() >= 4) &&
                (descriptorMode.substr(descriptorMode.size()-4,4) == "fist"))
        {
            CalculateFISTDescriptors
                (inputNoBorder, pointData, descriptorMode);
        }
        else if((descriptorMode.size() >= 4) &&
                (descriptorMode.substr(descriptorMode.size()-4,4) == "sift"))
        {
            CalculateFISTDescriptors
                (inputNoBorder, pointData, descriptorMode);
        }
        else if((descriptorMode.size() >= 4) &&
                (descriptorMode.substr(descriptorMode.size()-4,4) == "surf"))
        {
            InterestPointList outputPointList;
            pointData->ToPointList(outputPointList);
            CalculateSurfDescriptors
                (inputNoBorder, outputPointList, descriptorMode, mSurfParams[0],
                 mSurfParams[1], mSurfParams[2]);
            pointData->SetEmpty();
            pointData->FromPointList(outputPointList);
            outputPointList.EraseAndDelete();
        }
        else
        {
            InterestPointList outputPointList;
            pointData->ToPointList(outputPointList);
            CalculateRegionDescriptors
                (inputNoBorder, outputPointList, descriptorMode);
            pointData->SetEmpty();
            pointData->FromPointList(outputPointList);
            outputPointList.EraseAndDelete();
        }
    }


    ILOG_VAR_DEC;
};

ILOG_VAR_INIT(InterestPointFeature, Core.Feature);

} // namespace Koen
} // namespace Sandbox
} // namespace Impala

#endif

---