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

---

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 "Core/Matrix/MatTrace.h"
#include "Core/Matrix/MatEye.h"
#include "Util/Random.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/FeatureTableResult.h"
#include "Core/Feature/PointDescriptorTable.h"
#include "Core/Feature/ReadCodebook.h"
#include "Core/Table/Select.h"
#include "Core/Matrix/MatFunc.h"
#include "Persistency/RepositoryInFileSystem.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 std::vector<Vector::VectorTem<Real64>*> CodebookVectors;
    typedef Persistency::FileLocator FileLocator;
    typedef Persistency::File File;

    bool              mClusterInput;
    FeatureTable*     mClusterInputData;
    String            mClusterInputImage;

    InterestPointFeature(CmdOptions& options) :
        mVerbose(0),
        mStripBorder(0),
        mUseRecGauss(false),
        mKeepLimited(0),
        mCodebook(0),
        mDescriptorReduce(0),
        mExtraBoxes(0),
        mResultOutput(0),
        mClusterInput(false),
        mClusterInputData(0),
        mDynamicPointSelector(false)
    {
        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"))
            mRNG.RandomizeSeed();
        mTimeStats.AddGroup("detector");
        mTimeStats.AddGroup("descriptor");
        mTimeStats.AddGroup("projection");

        String detector    = options.GetString("detector");
        mName              = detector + "-" + options.GetString("descriptor");
        mFeatureDefinition = FeatureDefinition(options.GetString("featurePrefix") + 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";
            }
            else 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";
            }
            else 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";
            }
            else if(config == "pyramid-1x1-2x2-4x4")
            {
                // this is three levels classic
                config = "P1x1#0+P2x2#0+P2x2#1+P2x2#2+P2x2#3+P4x4#0+P4x4#1+P4x4#2+P4x4#3+P4x4#4+P4x4#5+P4x4#6+P4x4#7+P4x4#8+P4x4#9+P4x4#10+P4x4#11+P4x4#12+P4x4#13+P4x4#14+P4x4#15";
            }
            else if(config == "pyramid-1x1-2x2-4x4-1x3")
            {
                // this is three levels classic plus horizontal bars
                config = "P1x1#0+P2x2#0+P2x2#1+P2x2#2+P2x2#3+P4x4#0+P4x4#1+P4x4#2+P4x4#3+P4x4#4+P4x4#5+P4x4#6+P4x4#7+P4x4#8+P4x4#9+P4x4#10+P4x4#11+P4x4#12+P4x4#13+P4x4#14+P4x4#15+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",       "");
        options.AddOption(0, "codebookName",     "name",        "");

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

        /*** 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+...",       "");

        options.AddOption(0, "featurePrefix",   "prefix",       "");

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


    FeatureDefinition
    GetFeatureDefinition()
    {
        return mFeatureDefinition;
    }

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


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


    void
    SetCodebook(String filename, const Persistency::Locator& locator)
    {
        mLocator = locator;
        FileLocator codebookLoc(mLocator, filename);
        File file = RepositoryGetFile(codebookLoc, false, false);
        Util::IOBuffer* buf = file.GetReadBuffer();
        if (!mDescriptorReduceFilename.empty())
        {
            FileLocator fLoc(mLocator, mDescriptorReduceFilename);
            File f = RepositoryGetFile(fLoc, false, false);
            ReadRaw(mDescriptorReduce, f);
            mFeatureDefinition.AddParameter("red", MakeString(Matrix::MatNrCol(mDescriptorReduce)));
        }

        ILOG_DEBUG("SetCodebook after IOBuffer " << buf);
        if (!buf)
        {
            ILOG_ERROR("SetCodebook: IOBuffer is 0");
            return;
        }
        mCodebook = ReadCodebookFromBuffer(buf);
        mCodebook->SetFeatureDefinition(FeatureDefinition(filename));
        delete buf;
        
        ILOG_INFO("Reading codebook: " << codebookLoc <<
                  "; dimensionCount = " << mCodebook->GetFeatureVectorLength()
                  << "; size = " << mCodebook->Size());

        CmdOptions& options = CmdOptions::GetInstance();
        String CM = options.GetString("codebookMode");
        if(CM != "hard")
        {
            CM = StringReplace(CM, "llcdump", "llc");
            if(CM == "hardindex")
            {
                mFeatureDefinition.AddParameter
                    ("cm", CM + "-" + FileNameBase(filename));
            }
            else
            {
                mFeatureDefinition.AddParameter("cm", CM);
            }
        }
        if(options.GetString("codebookSigma") != "")
        {
            mFeatureDefinition.AddParameter
                ("cs", options.GetString("codebookSigma"));
        }
        if(options.GetString("codebookName") != "")
        {
            mFeatureDefinition.AddParameter
                ("cn", options.GetString("codebookName"));
        }

        if(!options.GetString("extraBoxes").empty())
        {
            FileLocator fLoc(mLocator, options.GetString("extraBoxes"));
            File f = RepositoryGetFile(fLoc, false, false);
            ReadRaw(mExtraBoxes, f);
        }
    }
    
    void
    SetResultOutput(FeatureTableResult* result)
    {
        mResultOutput = result;
    }

    bool
    FindInterestPoints(CmdOptions& options, Array::Array2dVec3UInt8* input,
                       Quid quid,     // identifies input image data
                       String outputFilename, bool borderAlreadyRemoved=false,
                       Util::IOBuffer* loadRegionsFile=0)
    {
        ILOG_DEBUG("start FindInterestPoints");
        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");
        mTimeStats.MeasureFirst();
        if(loadRegionsFile)
        {
            // load the points from file
            pointData = PointDescriptorTable::ImportFromBuffer(detectorMode, 
                                          loadRegionsFile, mStripBorder);
        }
        else
        {
            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");
        mTimeStats.MeasureNext();
        if(descriptorMode.size() != 0)
        {
            ComputeDescriptors(pointData, inputNoBorder, descriptorMode);

            ILOG_DEBUG("[Timing] Descriptors finished at " << timer.SplitTimeStr());
            if(mCodebook)
            {
                ILOG_DEBUG("preparing for codebook");
                /* project it onto a codebook */
                Timer timerProjection;
                if(mExtraBoxes)
                {
                    for(int j = 0; j < mExtraBoxes->CH(); j++)
                    {
                        if(*mExtraBoxes->CPB(0, j) == quid)
                        {
                            // it is a box for this image
                            Quid boxQuid = *mExtraBoxes->CPB(1, j);
                            UInt64 xlow = *mExtraBoxes->CPB(2, j);
                            UInt64 ylow = *mExtraBoxes->CPB(3, j);
                            UInt64 xhi = *mExtraBoxes->CPB(4, j);
                            UInt64 yhi = *mExtraBoxes->CPB(5, j);
                            UInt64 width = xhi - xlow;
                            UInt64 height = yhi - ylow;
                            // need to make a copy
                            PointDescriptorTable* table = new PointDescriptorTable(pointData->GetFeatureDefinition());
                            table->SetDescriptorLength(pointData->GetDescriptorLength());
                            for(int i = 0; i < pointData->Size(); i++)
                            {
                                Real64 x = pointData->GetX(i) - xlow;
                                Real64 y = pointData->GetY(i) - ylow;
                                if((x < 0) || (y < 0) || (x > width) || (y > height))
                                {
                                    continue;
                                }
                                table->Add(x, y, 0, 0, 0, pointData->Get6(i));
                            }
                            ILOG_DEBUG("extrabox:"<<xlow << " "<< ylow << " "<<xhi << " " << yhi << " " << table->Size() << " " << pointData->Size());
                            table = ProjectOntoCodebook(table,
                                            options.GetString("codebookMode"),
                                            width, height, boxQuid);
                            delete table;
                        }
                    }
                }
                else
                {
                    mTimeStats.MeasureNext();
                    pointData = ProjectOntoCodebook(pointData,
                                       options.GetString("codebookMode"),
                                       inputNoBorder->CW(), inputNoBorder->CH(),
                                       quid);
                    mTimeStats.MeasureLast();
                    ILOG_INFO("Overall " << mTimeStats.AsString());
                    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");
                    int n = pointData->GetDescriptorLength();
                    for(int i = 0; i < pointData->Size(); i++)
                    {
                        Vector::VectorTem<Real64> descriptor(n, 
                            pointData->GetDescriptorData(i), true);

                        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");
                }
                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(CodebookVectors& codebookVectors, int codeword, double x, 
                  double y, int imageWidth, int imageHeight, 
                  std::vector<Real64>& regionCounts)
    {
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            InterestPointSelector* selector = mPointSelector[ps];
            if(selector->Accept(x, y, imageWidth, imageHeight))
            {
                codebookVectors[ps]->Elem(codeword) =
                    codebookVectors[ps]->Elem(codeword) + 1;
                regionCounts[ps] = regionCounts[ps] + 1;
            }
        }
    }
    
    void
    ProjectUNC(CodebookVectors& codebookVectors, Mat* matrixDescriptors,
               int imageWidth, int imageHeight, std::vector<int>& X, 
               std::vector<int>& Y, std::vector<Real64>& 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 = Min(Int64(10000), Int64(20000000 / mCodebook->Size()));
        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++)
                        {
                            codebookVectors[ps]->Elem(i) =
                                codebookVectors[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
    ProjectLLC(CodebookVectors& codebookVectors, String codebookMode,
               Mat* matrixDescriptors, int imageWidth,
               int imageHeight, std::vector<int>& X, std::vector<int>& Y,
               std::vector<Real64>& regionCounts, Util::IOBuffer* dumpIndices=0)
    {
        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 knn = CmdOptions::GetInstance().GetInt("codebookSigma");
        Real64 beta = 1e-4;
        int step = Min(Int64(10000), Int64(20000000 / mCodebook->Size()));
        Mat* II = MatEye<Mat>(knn);
        MulVal(II, II, beta);
        Mat* onesVector = MatCreate<Mat>(knn, 1);
        SetVal(onesVector, 1.0);
        
        Mat* dumpedIndices = 0;
        if(dumpIndices)
        {
            dumpedIndices = MatCreate<Mat>(MatNrRow(matrixDescriptors), 2*knn+2);
        }
        String config = codebookMode;
        String absolute = StringSplit(config, '-');
        String aggregationMode = StringSplit(config, '-');
        String normMode = StringSplit(config, '-');
        ILOG_INFO("absolute=" << absolute << "; aggregationMode=" <<
                  aggregationMode << "; normMode=" << normMode);
        for(int kkk = 0; kkk < MatNrRow(matrixDescriptors); kkk += step)
        {
            int s = MatNrRow(matrixDescriptors) - kkk;
            if(step < s) s = step;
    
            Mat* partial = MatCreate<Mat>(s, MatNrCol(matrixDescriptors),
                                          MatE(matrixDescriptors, kkk, 0), true);

            // result will be a (#partial-descriptors, codebook size) matrix
            Mat* distances = MatNorm2DistTransposed(partial,
                                                    matrixCodebook);

            ILOG_DEBUG("Euclidean distance matrix took: " <<
                       timerProjection.SplitTimeStr());
                       

            std::vector<std::pair<Real64, int> > distancesAndIndices;
            // find k closest codebook elements
            for(int vj = 0; vj < s; vj++)
            {
                distancesAndIndices.clear();
                distancesAndIndices.reserve(mCodebook->Size());
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    distancesAndIndices.push_back(
                                std::make_pair(*MatE(distances, vj, i), i));
                }
                std::nth_element(distancesAndIndices.begin(), 
                                 distancesAndIndices.begin() + knn, 
                                 distancesAndIndices.end());
                // the first k elements are now the k smallest
                
                // idx = indices[:, descriptorIdx]
                // z = codebook[idx, :] - 
                //          numpy.tile(descriptors[descriptorIdx, :], (knn, 1))
                Mat* z = MatCreate<Mat>(knn, MatNrCol(matrixDescriptors));
                for(int kk = 0; kk < knn; kk++)
                {
                    int closestCodebookIndex = distancesAndIndices[kk].second;
                    for(int i = 0; i < MatNrCol(matrixDescriptors); i++)
                    {
                        *MatE(z, kk, i) = *MatE(matrixCodebook, 
                                                closestCodebookIndex, i) - 
                                          *MatE(partial, vj, i);
                    }
                }
                
                // C = numpy.dot(z, z.transpose());
                Mat* zT = MatTranspose(z);
                Mat* C = MatMul(z, zT);
                delete z; delete zT;
                
                // C = C + II * beta * numpy.trace(C);
                Mat* tmp = 0;
                MulVal(tmp, II, MatTrace(C));
                Add(C, C, tmp);
                delete tmp;
                
                Mat* invertedC = MatInverse(C);
                Mat* w = MatMul(invertedC, onesVector);
                delete invertedC;
                delete C;
                
                // should we take the absolute value of w ??? TODO
                
                Real64 total = 0.0;
                if(absolute == "llca")
                {
                    for(int kk = 0; kk < knn; kk++)
                        total += fabs(*MatE(w, kk, 0));
                    for(int kk = 0; kk < knn; kk++)
                    {
                        *MatE(w, kk, 0) = fabs(*MatE(w, kk, 0)) / total;
                    }
                }
                else
                {
                    for(int kk = 0; kk < knn; kk++)
                        total += *MatE(w, kk, 0);
                    for(int kk = 0; kk < knn; kk++)
                    {
                        *MatE(w, kk, 0) = *MatE(w, kk, 0) / total;
                    }
                }
                
                if(aggregationMode == "avg")
                {
                    for(int ps = 0; ps < mPointSelector.size(); ps++)
                    {
                        if(mPointSelector[ps]->Accept(X[kkk+vj], Y[kkk+vj],
                                                      imageWidth, imageHeight))
                        {
                            for(int kk = 0; kk < knn; kk++)
                            {
                                int closestCodebookIndex = 
                                                distancesAndIndices[kk].second;
                                codebookVectors[ps]->Elem(closestCodebookIndex) =
                                    codebookVectors[ps]->Elem(closestCodebookIndex) + *MatE(w, kk, 0);
                            }
                            regionCounts[ps] = regionCounts[ps] + 1;
                        }
                    }
                }
                else if(aggregationMode == "max" || aggregationMode.empty())
                {
                    for(int ps = 0; ps < mPointSelector.size(); ps++)
                    {
                        if(mPointSelector[ps]->Accept(X[kkk+vj], Y[kkk+vj],
                                                      imageWidth, imageHeight))
                        {
                            for(int kk = 0; kk < knn; kk++)
                            {
                                int closestCodebookIndex = 
                                                distancesAndIndices[kk].second;
                                codebookVectors[ps]->Elem(closestCodebookIndex) =
                                    Max(codebookVectors[ps]->Elem(closestCodebookIndex), *MatE(w, kk, 0));
                            }
                            regionCounts[ps] = 1;
                        }
                    }
                }
                
                // dump
                if(dumpIndices)
                {
                    *MatE(dumpedIndices, kkk+vj, 0) = X[kkk+vj];
                    *MatE(dumpedIndices, kkk+vj, 1) = Y[kkk+vj];
                    for(int kk = 0; kk < knn; kk++)
                    {
                        int closestCodebookIndex = 
                                        distancesAndIndices[kk].second;
                        *MatE(dumpedIndices, kkk+vj, kk*2 + 2) = closestCodebookIndex;
                        *MatE(dumpedIndices, kkk+vj, kk*2 + 3) = *MatE(w, kk, 0);
                    }
                }
                delete w;
            }
            delete partial;
            delete distances;
            ILOG_DEBUG("Codebook projection (piece) took: " <<
                       timerProjection.SplitTimeStr() << " " << mPointSelector.size());
        }
        delete II;
        delete onesVector;
        ILOG_INFO("Codebook projection (LLC/" << MatNrRow(matrixDescriptors) <<
                  ") took: " << timerProjection.SplitTimeStr());
        if(mDescriptorReduce)
        {
            delete matrixDescriptors;
        }
        if(normMode.empty())
        {
            // for backward compatibility
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                regionCounts[ps] = 1;
            }
        }
        else if(normMode == "L2")
        {
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                Real64 tmp = 0.0;
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    tmp += codebookVectors[ps]->Elem(i) * codebookVectors[ps]->Elem(i);
                }
                regionCounts[ps] = sqrt(tmp);
            }
        }
        else if(normMode == "L2a")
        {
            Real64 tmp = 0.0;
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    tmp += codebookVectors[ps]->Elem(i) * codebookVectors[ps]->Elem(i);
                }
            }
            for(int ps = 0; ps < mPointSelector.size(); ps++)
                regionCounts[ps] = sqrt(tmp);
        }
        else if(normMode == "L1")
        {
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                Real64 tmp = 0.0;
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    tmp += codebookVectors[ps]->Elem(i);
                }
                regionCounts[ps] = tmp;
            }
        }
        else if(normMode == "L1a")
        {
            Real64 total = 0.0;
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                Real64 tmp = 0.0;
                for(int i = 0; i < mCodebook->Size(); i++)
                {
                    tmp += codebookVectors[ps]->Elem(i);
                }
                regionCounts[ps] = tmp;
                total += tmp;
            }
            for(int ps = 0; ps < mPointSelector.size(); ps++)
            {
                regionCounts[ps] = total;
            }
        }

        if(dumpIndices)
        {
            WriteRaw(dumpedIndices, dumpIndices, 3);
            delete dumpedIndices;
        }
    }

    void
    ProjectHardMatrix(CodebookVectors& codebookVectors, Mat* matrixDescriptors,
                      int imageWidth, int imageHeight, std::vector<int>& X,
                      std::vector<int>& Y, std::vector<Real64>& 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(codebookVectors, 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;
        if(mDescriptorReduce)
        {
            delete matrixDescriptors;
        }
    }
    
    
    void
    ProjectForest(CodebookVectors& codebookVectors, Mat* descriptors, int imageWidth,
                  int imageHeight, std::vector<int>& X, std::vector<int>& Y,
                  std::vector<Real64>& 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(codebookVectors, 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 quid)
    {
        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(quid) + ".raw";
            FileLocator loc(mLocator, filename);
            File file = RepositoryGetFile(loc, false, false);
            Util::IOBuffer* buf = file.GetReadBuffer();

            // 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<Real64> regionCounts;
        int codebookLength = mCodebook->Size();
        if(codebookMode == "forest")
        {
            codebookLength = GetCodebookLength(mCodebook);
            ILOG_DEBUG("forest codebook length = "<< codebookLength);
        }

        CodebookVectors codebookVectors;
        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;
            }
            codebookVectors.push_back(codebookVector);
            regionCounts.push_back(0);
        }
        ILOG_DEBUG("done making codebook vector");
        if(pointData->Size() == 0)
        {
            if(mResultOutput)
            {
                for(int i = 0; i < codebookVectors.size(); i++)
                {
                    mResultOutput->AddVector(quid, i, *codebookVectors[i]);
                }
            }
            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(codebookVectors, descriptors, imageWidth, imageHeight, 
                       X, Y, regionCounts);
        // hard assignment, matrix style
        if((codebookMode == "hardmatrix") || (codebookMode == "hard"))
            ProjectHardMatrix(codebookVectors, 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;
            for(int i = 0; i < codebookVectors.size(); i++)
                delete codebookVectors[i];
            codebookVectors.clear();
            return pointData;
        }
        if(codebookMode == "forest")
            ProjectForest(codebookVectors, descriptors, imageWidth, imageHeight,
                          X, Y, regionCounts);
        // LLC assignment, dump the indices to disk
        if(StringStartsWith(codebookMode, "llcdump"))
        {
            String folder = "ClusterIndices/" + GetFeatureName();
            String filename = PathJoin(folder, MakeString(quid) + ".raw");
            FileLocator loc(mLocator, filename);
            File file = RepositoryGetFile(loc, true, false);
            Util::IOBuffer* buf = file.GetWriteBuffer();
            codebookMode = StringReplace(codebookMode, "llcdump", "llc");
            ProjectLLC(codebookVectors, codebookMode, descriptors, imageWidth, 
                       imageHeight, X, Y, regionCounts, buf);
            delete buf;
        }
        else if(StringStartsWith(codebookMode, "llc"))
        {
            ProjectLLC(codebookVectors, codebookMode, descriptors, imageWidth,
                       imageHeight, X, Y, regionCounts);
        }
        delete descriptors;

        // normalize so it sums to one (although this is not guaranteed,
        // it depends on what the functions called above put into the
        // regionCount variable)
        for(int ps = 0; ps < mPointSelector.size(); ps++)
        {
            // debugStream << "\n\npyramid part #"<<ps<<"\n";
            Real64 regionCount = regionCounts[ps];
            ILOG_DEBUG("regCount " << regionCount);
            if(regionCount == 0)
                regionCount = 1;
            for(int i = 0; i < codebookLength; i++)
            {
                codebookVectors[ps]->Elem(i) =
                    codebookVectors[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, *(codebookVectors[ps]));
        }

        if(mResultOutput)
        {
            for(int i = 0; i < codebookVectors.size(); i++)
            {
                mResultOutput->AddVector(quid, i, 
                                         *codebookVectors[i]);
            }
        }

        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;
        }
        std::random_shuffle(indices.begin(), indices.end(), mRNG);
        
        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"))
        {
            using namespace Core::Matrix;
            CalculateRegionDescriptors(inputNoBorder, pointData,
                                       "huehistogram");

            // now fix the range of hue...
            Mat* part1 = MatCreate<Mat>(pointData->Size(), pointData->GetDescriptorLength());
            for(int i = 0; i < pointData->Size(); i++)
            {
                Real64* ptr = pointData->GetDescriptorData(i);
                for(int j = 0; j < MatNrCol(part1); j++)
                {
                    *MatE(part1, i, j) = static_cast<int>(ptr[j] * 512);
                }
            }
            
            CalculateFISTDescriptors(inputNoBorder, pointData, "sift");
            Mat* part0 = pointData->GetDescriptorMatrix();
            std::vector<Mat*> matrices;
            matrices.push_back(part0);
            matrices.push_back(part1);
            Matrix::Mat* descriptors =
                Matrix::MatConcatenateHorizontal(matrices);
        
            // takes ownership of memory
            pointData->ReplaceAllDescriptors(descriptors);
            
            delete part1;
       }
        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"))
        {
            CalculateSurfDescriptors
                (inputNoBorder, pointData, descriptorMode, mSurfParams[0],
                 mSurfParams[1], mSurfParams[2]);
        }
        else if(descriptorMode == "none")
        {
            // do nothing, assume the user knows what he is doing
        }
        else
        {
            CalculateRegionDescriptors
                (inputNoBorder, pointData, descriptorMode);
        }
    }

    // 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;

    FeatureTableResult* mResultOutput;
    Array::Array2dScalarUInt64* mExtraBoxes;

    std::vector<InterestPointSelector*> mPointSelector;
    bool                         mDynamicPointSelector;
    String                       mDynamicPointSelectorConfig;
    Persistency::Locator    mLocator;
    std::vector<Real64> mDenseSamplingScales;
    int                 mSurfParams[3];

    Util::Random mRNG;
    Util::TimeStats mTimeStats;

    ILOG_VAR_DEC;
};

ILOG_VAR_INIT(InterestPointFeature, Core.Feature);

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

#endif

---