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

---

RegionDescriptor.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_RegionDescriptor_h
#define Impala_Core_Feature_RegionDescriptor_h

#include "Core/Array/MakeCircle.h"
#include "Core/Array/ProjectRange.h"
#include "Core/Array/Resample.h"
#include "Core/Histogram/MakeHistogram1d.h"
#include "Core/Histogram/Histogram1dSet.h"
#include "Core/Array/Rgb2Ooo.h"
#include "Core/Array/RGB2rg.h"
#include "Core/Array/MakeFromData.h"
#include "Core/Array/ColorMomentsKoen.h"
#include "Core/Array/Atan2.h"
#include "Core/Array/Sqrt.h"
#include "Core/Array/GaussDerivative.h"
#include "Core/Array/AddVal.h"

#include "Core/Vector/VectorSet.h"
#include "Core/Histogram/MakeHistogram1dSet.h"
#include "Core/Feature/PointDescriptorTable.h"
#include "Persistency/RepositoryInFileSystem.h"

namespace Impala
{
namespace Core
{
namespace Feature
{


Array::Array2dScalarReal64*
CreateCircleMask(int patchSize) 
{
    using namespace Impala::Core::Array;
    Array2dScalarReal64* circleMask = 0;
    circleMask = ArrayCreate<Array2dScalarReal64>(patchSize, patchSize);
    MakeCircle(circleMask, patchSize / 2, patchSize / 2, Real64(patchSize) / 2.0, 1, 0);
    return circleMask;
}


void
ComputeRegionDescriptor(std::vector<Real64>& vec, Array::Array2dVec3UInt8* src,
                        String descriptor, bool useCircularMask, int srcCenterX,
                        int srcCenterY)
{
    using namespace Impala::Core::Array;
    using namespace Impala::Core::Histogram;

    /* note: ensure that vec is empty when you call this function */
    //WriteRaw(src, String("src.raw"), 0);

    /* certain descriptors can only be computed over square images */
    bool imgIsSquare = (src->CW() == src->CH());
    int binCount = 15;
    Real64 uniformPatchSize = 61;           // when making the src uniformly sized; resize to 61x61

    /* support bin count changing */
    CmdOptions& options = CmdOptions::GetInstance();
    binCount = options.GetInt("descriptorBinCount");

    Array2dScalarReal64* circleMask = 0;
    if(useCircularMask) {
        if(!imgIsSquare) {
            std::cerr << "[ERROR] Using circularMask on a non-square source!" << std::endl;
            throw "Using circularMask on a non-square source!";
        }
        circleMask = CreateCircleMask(src->CW());
    }
    
    if(descriptor == "opponenthistogram") {
        /***** OPPONENT COLOR HISTOGRAM *****/
        vec.reserve(binCount * 3);
    
        // color conversion
        Array2dVec3Real64* srcOpponent = 0;
        Rgb2Ooo(srcOpponent, src);
        //WriteRaw(srcOpponent, String("srcopp.raw"), 0);
        
        double lower[3] = {0, -70, -70};
        double upper[3] = {270, 70, 70};

        if(useCircularMask) {
            AddVal(circleMask, circleMask, -1);
            MulVal(circleMask, circleMask, 1000);   // now it is 0 in the circle and -1000 elsewhere
        }

        for (int d=1 ; d<=src->ElemSize() ; d++)
        {
            Array2dScalarReal64* comp = 0;
            ProjectRange(comp, srcOpponent, d);
            
            if(useCircularMask) {
                /* add them together, making values outside the circle very negative (i.e. an outlier) */
                Add(comp, comp, circleMask);
                //WriteRaw(comp, String("x.raw"), 0);
            }
            
            Histogram1dTem<Real64> hist(lower[d-1], upper[d-1], binCount, 0);
            MakeHistogram1d(&hist, comp);
            hist.Normalize();
            if(hist.TotalWeight() == 0) {
                /* debug code */
                hist.Dump(true);
                std::cout << "[WARNING] Opponent color histogram 100% outliers: " << srcCenterX << " " << srcCenterY << std::endl;
                Persistency::File file = Persistency::RepositoryGetFile("outlier.raw", true, false);
                WriteRaw(comp, file, 1);
            }
            delete comp;
            for(int i = 0; i < hist.Size(); i++) {
                vec.push_back(hist.Elem(i));
            }
        }
        delete srcOpponent;
    } else if(descriptor == "nrghistogram") {
        /***** rgb HISTOGRAM (normalized rgb) *****/
        vec.reserve(binCount * 2);

        // color conversion
        Array2dVec3Real64* src_rg = 0;
        RGB2rg(src_rg, src);
    
        double lower[3] = {0, 0, 0};
        double upper[3] = {1, 1, 1};

        if(useCircularMask) {
            AddVal(circleMask, circleMask, -1);
            MulVal(circleMask, circleMask, 1000);   // now it is 0 in the circle and -1000 elsewhere
        }

        for (int d=1 ; d<=2 ; d++)
        {
            Array2dScalarReal64* comp = 0;
            ProjectRange(comp, src_rg, d);
            
            if(useCircularMask) {
                /* add them together, making values outside the circle very negative (i.e. an outlier) */
                Add(comp, comp, circleMask);
                //WriteRaw(comp, String("x.raw"), 0);
            }
            
            Histogram1dTem<Real64> hist(lower[d-1], upper[d-1], binCount, 0);
            MakeHistogram1d(&hist, comp);
            hist.Normalize();
            //hist.Dump(true);
            //std::cout << comp->CW() << std::endl;
            delete comp;
            for(int i = 0; i < hist.Size(); i++) {
                vec.push_back(hist.Elem(i));
            }
        }
        delete src_rg;
    } else if(descriptor == "rgbhistogram") {
        /***** RGB HISTOGRAM *****/
        vec.reserve(binCount * 3);
    
        double lower[3] = {0, 0, 0};
        double upper[3] = {256, 256, 256};

        if(useCircularMask) {
            AddVal(circleMask, circleMask, -1);
            MulVal(circleMask, circleMask, 1000);   // now it is 0 in the circle and -1000 elsewhere
        }

        for (int d=1 ; d<=src->ElemSize() ; d++)
        {
            Array2dScalarReal64* comp = 0;
            ProjectRange(comp, src, d);
            
            if(useCircularMask) {
                /* add them together, making values outside the circle very negative (i.e. an outlier) */
                Add(comp, comp, circleMask);
                //WriteRaw(comp, String("x.raw"), 0);
            }
            
            Histogram1dTem<Real64> hist(lower[d-1], upper[d-1], binCount, 0);
            MakeHistogram1d(&hist, comp);
            hist.Normalize();
            //hist.Dump(true);
            //std::cout << comp->CW() << std::endl;
            delete comp;
            for(int i = 0; i < hist.Size(); i++) {
                vec.push_back(hist.Elem(i));
            }
        }
    } else if(descriptor == "transformedcolorhistogram") {
        /***** TRANSFORMED RGB HISTOGRAM *****/
        vec.reserve(binCount * 3);
    
        double lower[3] = {-2, -2, -2};
        double upper[3] = { 2,  2,  2};

        for (int d=1 ; d<=src->ElemSize() ; d++)
        {
            Array2dScalarReal64* comp = 0;
            ProjectRange(comp, src, d);
            
            Real64 pixelCount = comp->CW() * comp->CH();
            if(useCircularMask) {
                /* make values outside the circle zero */
                Mul(comp, comp, circleMask);
                pixelCount = PixSum(circleMask);
            }

            Real64 mu = PixSum(comp) / pixelCount;
            AddVal(comp, comp, -mu);
            if(useCircularMask) {
                /* make values outside the circle zero */
                Mul(comp, comp, circleMask);
            }
            //std::cout << "0-mean: " << mu << " " << PixSum(comp) << " " << PixSum(comp) / pixelCount << std::endl;

            Array2dScalarReal64* stddev = 0;
            Mul(stddev, comp, comp);    // square it -> (x - mu)^2
            if(useCircularMask) {
                /* make values outside the circle zero */
                Mul(stddev, stddev, circleMask);
            }
            Real64 standardDeviation = sqrt(PixSum(stddev) / pixelCount);
            MulVal(comp, comp, 1.0 / standardDeviation);
            //std::cout << PixSum(comp) / pixelCount << std::endl;
            
            Histogram1dTem<Real64> hist(lower[d-1], upper[d-1], binCount, 0);
            MakeHistogram1d(&hist, comp);
            hist.Normalize();
            //hist.Dump(true);
            delete comp;
            delete stddev;
            for(int i = 0; i < hist.Size(); i++) {
                vec.push_back(hist.Elem(i));
            }
        }
//    } else if(descriptor == "nu-colormoments") {
//        /***** non-uniform COLOR MOMENTS *****/
//        ColorMoments cm;
//        cm.GetMoments(vec, src, 0, 2, useCircularMask, circleMask, false);
//    } else if(descriptor == "nu-spatialcolormoments") {
//        /***** non-uniform SPATIAL COLOR MOMENTS *****/
//        ColorMoments cm;
//        cm.GetMoments(vec, src, 1, 2, useCircularMask, circleMask, false);
//    } else if(descriptor == "nupv-colormoments") {
//        /***** non-uniform position-variant COLOR MOMENTS *****/
//        ColorMoments cm;
//        cm.SetCoordinateShift(srcCenterX, srcCenterY);
//        cm.GetMoments(vec, src, 0, 2, useCircularMask, circleMask, false);
//    } else if(descriptor == "nupv-spatialcolormoments") {
//        /***** non-uniform position-variant SPATIAL COLOR MOMENTS *****/
//        ColorMoments cm;
//        cm.SetCoordinateShift(srcCenterX, srcCenterY);
//        cm.GetMoments(vec, src, 1, 2, useCircularMask, circleMask, false);
//    } else if(descriptor == "colormoments") {
//        /***** COLOR MOMENTS *****/
//        Array2dVec3UInt8* uniformSrc = 0;
//        Resample(uniformSrc, src, uniformPatchSize, uniformPatchSize);
//        if(useCircularMask) delete circleMask;
//        circleMask = CreateCircleMask(uniformSrc->CW());
//        
//        ColorMoments cm;
//        cm.GetMoments(vec, uniformSrc, 0, 2, useCircularMask, circleMask, false);
//        delete uniformSrc;
    } else if(descriptor == "colormoments") {
        /***** SPATIAL COLOR MOMENTS *****/
        Array2dVec3UInt8* uniformSrc = 0;
        Resample(uniformSrc, src, uniformPatchSize, uniformPatchSize);
        if(useCircularMask) delete circleMask;
        circleMask = CreateCircleMask(uniformSrc->CW());

        ColorMoments cm;
        cm.GetMoments(vec, uniformSrc, 1, 2, useCircularMask, circleMask, false);
        delete uniformSrc;
//    } else if(descriptor == "spatialcolormoments3") {
//        /***** SPATIAL COLOR MOMENTS up to degree 3 *****/
//        Array2dVec3UInt8* uniformSrc = 0;
//        Resample(uniformSrc, src, uniformPatchSize, uniformPatchSize);
//        if(useCircularMask) delete circleMask;
//        circleMask = CreateCircleMask(uniformSrc->CW());
//
//        ColorMoments cm;
//        cm.GetMoments(vec, uniformSrc, 1, 3, useCircularMask, circleMask, false);
//        delete uniformSrc;
    } else if(descriptor == "colormomentinvariants") {
        /***** COLOR MOMENT INVARIANTS *****/
        Array2dVec3UInt8* uniformSrc = 0;
        Resample(uniformSrc, src, uniformPatchSize, uniformPatchSize);
        if(useCircularMask) delete circleMask;
        circleMask = CreateCircleMask(uniformSrc->CW());

        ColorMoments cm;
        cm.GetInvariants(vec, uniformSrc, useCircularMask, circleMask, false);
        delete uniformSrc;
    } else if(descriptor == "huehistogram") {
        /***** HUE DESCRIPTOR *****/
        /* by Joost van de Weijer, ECCV2006. Not included: Gaussian weighting mask */

        // default to the setting used by the paper
        if(binCount == 15) binCount = 37;
        
        //H=atan2((patches_R+patches_G-2*patches_B),sqrt(3)*(patches_R-patches_G))+pi;
        Array2dScalarReal64* H = 0;

        Array2dScalarReal64* R = 0;
        Array2dScalarReal64* G = 0;
        Array2dScalarReal64* B = 0;
        ProjectRange(R, src, 1);
        ProjectRange(G, src, 2);
        ProjectRange(B, src, 3);

        Array2dScalarReal64* X = 0;
        Add(X, R, G);
        MulVal(B, B, -1);
        Add(X, X, B);
        Add(X, X, B);
        Array2dScalarReal64* Y = 0;
        MulVal(G, G, -1);
        Add(Y, R, G);
        MulVal(Y, Y, sqrt(3.0));
        Atan2(H, X, Y);
        Real64 myPI = 3.1415927;
        AddVal(H, H, myPI);
        delete X; X = 0;
        delete Y; Y = 0;

        /* reset to correct values */
        ProjectRange(R, src, 1);
        ProjectRange(G, src, 2);
        ProjectRange(B, src, 3);
        
        //saturation=sqrt(2/3*(patches_R.^2+patches_G.^2+patches_B.^2-patches_R.*(patches_G+patches_B)-patches_G.*patches_B)+0.01);
        Array2dScalarReal64* Saturation = 0;
        // patches_R.^2+patches_G.^2+patches_B.^2
        Array2dScalarReal64* R2 = 0;
        Array2dScalarReal64* G2 = 0;
        Array2dScalarReal64* B2 = 0;
        Mul(R2, R, R);
        Mul(G2, G, G);
        Mul(B2, B, B);
        Add(X, R2, G2);
        Add(X, X, B2);
        delete R2; delete G2; delete B2;
        
        // previous-patches_R.*(patches_G+patches_B)
        Add(Y, G, B);
        Mul(Y, R, Y);
        MulVal(Y, Y, -1);
        Add(X, X, Y);
        
        // -patches_G.*patches_B
        Mul(Y, G, B);
        MulVal(Y, Y, -1);
        Add(X, X, Y);
        
        MulVal(X, X, 2.0 / 3.0);
        AddVal(X, X, 0.01);
        Sqrt(Saturation, X);

        delete R; delete G; delete B; delete X; delete Y;

        Histogram1dTem<Real64> hist(0, 2 * myPI, binCount, 0);
        if(useCircularMask) {
            Mul(Saturation, Saturation, circleMask);   // everything not in the mask will not be included
        }
        MakeHistogram1d(&hist, H, Saturation);
        hist.Normalize();
        //hist.Dump(true);
        for(int i = 0; i < hist.Size(); i++) {
            vec.push_back(hist.Elem(i));
        }
        delete H; delete Saturation;
        
    } else {
        std::cerr << "[ERROR] Invalid descriptor name: " << descriptor << std::endl;
        throw "Invalid descriptor name";
    }
    if(circleMask) delete circleMask;
}


void
MakePatch(Array::Array2dVec3UInt8*& patch, Array::Array2dVec3UInt8* img, int x,
          int y, int pixelsAround)
{
    Geometry::Rectangle regionRect(x - pixelsAround, y - pixelsAround,
                                   x + pixelsAround, y + pixelsAround);
    if((regionRect.mLeft   < - img->BW()) || 
       (regionRect.mTop    < - img->BH()) ||
       (regionRect.mRight  >= img->CW() + img->BW()) ||
       (regionRect.mBottom >= img->CH() + img->BH()) ) {
       
        //std::cerr << "[WARNING] Decreased pixelsAround to 120 in MakePatch to avoid going beyond borders of image, value was " << pixelsAround << std::endl;
        
        pixelsAround = 120;    // this should always work
        regionRect.mLeft = x - pixelsAround;
        regionRect.mTop = y - pixelsAround;
        regionRect.mRight = x + pixelsAround;
        regionRect.mBottom = y + pixelsAround;
        /*
        std::cerr << "[WARNING] Segfault possible due to too small border in MakePatch!" << std::endl;
        std::cerr << regionRect.mLeft   <<" "<< - img->BW() <<" "<< 
                     regionRect.mTop    <<" "<< - img->BH() <<" "<<
                     regionRect.mRight  <<" "<< img->CW() + img->BW() <<" "<<
                     regionRect.mBottom <<" "<< img->CH() + img->BH() << std::endl;
        std::cerr << "x=" << x << " y=" << y << " size=" << pixelsAround << std::endl; */
    }
    patch = MakeRoi(img, regionRect);
}


void
CalculateRegionDescriptors(Array::Array2dVec3UInt8* inputNoBorder,
                           PointDescriptorTable* pointData,
                           String descriptor)
{
    using namespace Impala::Core::Array;
    using namespace Impala::Core::Geometry;

    // start of descriptor computation code
    Real64 scaleMultiplier = 10;     // used to be 4
    Real64 extraBorder = 120;
    
    Array2dVec3UInt8* input =
        ArrayCreate<Array2dVec3UInt8>(inputNoBorder->CW(), inputNoBorder->CH(), extraBorder, extraBorder);
    MakeFromData2<Array2dVec3UInt8,Array2dVec3UInt8>(input, inputNoBorder);
    /* now mirror the data into the border */
    Pattern::FuncBorderMirror2d(input, extraBorder, extraBorder);

    /* support bin count changing */
    CmdOptions& options = CmdOptions::GetInstance();
    if(options.GetInt("descriptorPreSmoothing") != 0) {
        /* pre-smooth the input */
        Real64 sigmaSmooth = 1.0;

        std::vector<Array2dScalarReal64*> resList;
        for (int d=1 ; d<=input->ElemSize() ; d++)
        {
            Array2dScalarReal64* comp = 0;
            ProjectRange(comp, input, d);
            GaussDerivative(comp, comp, sigmaSmooth, 0, 0, 3.0);
            resList.push_back(comp);
            delete comp;
        }
        delete input;
        input = 0;
        MakeFrom3Images(input, resList[0], resList[1], resList[2]);
        ArrayListDelete(&resList);
    }

    //ArraySystem::Instance().MarkMemoryUsage(true);          // DEBUG
    std::vector<Real64> tmpDescriptor;
    for(int i = 0; i < pointData->Size(); i++)
    {
        /* assumes a circle */

        /* this code makes the (very strong) assumption that (X,Y) are integers, which is true in the interest point datastructure */
        int pixelsAround = floor(pointData->GetScale(i) * scaleMultiplier + 0.5);
        Array2dVec3UInt8* patch = 0;
        int x = floor(pointData->GetX(i) + 0.5);
        int y = floor(pointData->GetY(i) + 0.5);

        MakePatch(patch, input, x, y, pixelsAround);
        
        /* compute a descriptor over the patch */
        tmpDescriptor.clear();
        ComputeRegionDescriptor(tmpDescriptor, patch, descriptor, true, x, y);
        if(pointData->GetDescriptorLength() != tmpDescriptor.size())
        {
            pointData->SetDescriptorLength(tmpDescriptor.size());
        }
        Real64* ptr = pointData->GetDescriptorData(i);
        for(int j = 0; j < tmpDescriptor.size(); j++)
        {
            ptr[j] = tmpDescriptor[j];
        }
        delete patch;
        
//        } else if(point->geometryType() == "RectangleZ") {
//            /* cast to a rectangle */
//            InterestRectangleZ* rectangle(dynamic_cast<InterestRectangleZ*>(point));
//
//            /* this code makes the (very strong) assumption that (X,Y) are integers, which is true in the interest point datastructure */
//            Array2dVec3UInt8* patch = MakeRoi(input, rectangle->GetRect());
//                   
//            /* compute a descriptor over the patch; do *not* use a circular mask */
//            // do not clear it first, then I cannot concatenate descriptors!
//            //point->mDescriptor.clear();
//
//            ComputeRegionDescriptor(rectangle->mDescriptor, patch, descriptor, false, rectangle->mX, rectangle->mY);
//
//            delete patch;
//        } else {
//            throw "[CalculateDescriptors] Unknown geometry type: " + point->geometryType();
//        }
        //ArraySystem::Instance().MemoryUsageSinceMark(true);     // DEBUG
    }
    delete input;
}

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

#endif

---