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

---

KoenFIST2.h

Go to the documentation of this file.

#ifndef Impala_Core_Feature_KoenFIST_h
#define Impala_Core_Feature_KoenFIST_h

#include "Core/Array/RecGauss.h"
#include "Core/Array/Sqrt.h"
#include "Core/Matrix/Mat.h"
#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatConcatenate.h"
#include "Core/Feature/PointDescriptorTable.h"
#include "Core/Feature/GetColorChannels.h"
#include "Core/Array/Pow.h"
#include "Core/Array/sRGB2RGB.h"
#ifdef CUDA
#include "Link/Cuda/KoenSIFT.h"
#endif

namespace Impala
{
namespace Core
{
namespace Feature
{


#ifdef PI
#undef PI
#endif
#define PI 3.141592654


/******************* original FIST *********************/

template <int INDEX_SIZE, int ORI_SIZE>
class FISTDescriptor
{
public:
    static const int VecLength = (INDEX_SIZE * INDEX_SIZE * ORI_SIZE);
    bool noGaussianWeighting;
    Real64 mOutputMultiplier;
    int mCustomNorm;
    bool mNoMaxBins;
    
    FISTDescriptor(bool noGaussianWeighting=false)
        : noGaussianWeighting(noGaussianWeighting), mOutputMultiplier(512), mCustomNorm(-1), mNoMaxBins(false)
    {
    }
    
    
    inline Real64
    CanonicalizeAngle(Real64 angle)
    {
        // make sure an angle lies between 0 and 2PI
        while (angle > 2*(Real64)PI)
        {
            angle -= 2*(Real64)PI;
        }
        while (angle < 0.0)
        {
            angle += 2*(Real64)PI;
        }
        return angle;
    }
    
    
    void
    PutInVector(Real64* descriptor,
                      Real64 mag, Real64 ori, Real64 rx, Real64 cx)
    {
        int ri, ci, oi;
        Real64 oval, rfrac, cfrac, ofrac;
       
        ori = CanonicalizeAngle(ori);
        oval = ORI_SIZE * ori / (2*(Real64)PI);
       
        ri = (int)floor(rx);
        ci = (int)floor(cx);
        oi = (int)floor(oval);
    
        // fractional parts of bin placement
        rfrac = rx - ri;
        cfrac = cx - ci;
        ofrac = oval - oi;
        if(!(rfrac >= 0.0 && rfrac <= 1.0 && cfrac >= 0.0 && cfrac <= 1.0 && ofrac >= 0.0 && ofrac <= 1.0))
        {
            std::cerr << "[Failure] Fractions outside valid range " << rfrac << " " << cfrac << " " << ofrac << std::endl;
            return;
        }
        if(oi >= ORI_SIZE) oi = 0;
    
        int rindex = ri;
        Real64 v;
        if((rindex >= 0) && (rindex < INDEX_SIZE))
        {
            int cindex = ci;
            if((cindex >= 0) && (cindex < INDEX_SIZE))
            {
                // case 0
                int oindex = oi;
                v = mag * (1.0 - rfrac) * (1.0 - cfrac) * (1.0 - ofrac);
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
             
                // case 1
                oindex = oi + 1;
                if(oindex >= ORI_SIZE) oindex = 0;
                v = mag * (1.0 - rfrac) * (1.0 - cfrac) * ofrac;
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
            }
             
            cindex = ci + 1;
            if((cindex >= 0) && (cindex < INDEX_SIZE))
            {
                // case 2
                int oindex = oi;
                v = mag * (1.0 - rfrac) * cfrac * (1.0 - ofrac);
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
             
                // case 3
                oindex = oi + 1;
                if(oindex >= ORI_SIZE) oindex = 0;
                v = mag * (1.0 - rfrac) * cfrac * ofrac;
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
            }
        }
        
        rindex = ri + 1;
        if((rindex >= 0) && (rindex < INDEX_SIZE))
        {
            int cindex = ci;
            if((cindex >= 0) && (cindex < INDEX_SIZE))
            {
                // case 4
                int oindex = oi;
                v = mag * rfrac * (1.0 - cfrac) * (1.0 - ofrac);
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
             
                // case 5
                oindex = oi + 1;
                if(oindex >= ORI_SIZE) oindex = 0;
                v = mag * rfrac * (1.0 - cfrac) * ofrac;
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
            }
             
            cindex = ci + 1;
            if((cindex >= 0) && (cindex < INDEX_SIZE))
            {
                // case 6
                int oindex = oi;
                v = mag * rfrac * cfrac * (1.0 - ofrac);
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
             
                // case 7
                oindex = oi + 1;
                if(oindex >= ORI_SIZE) oindex = 0;
                v = mag * rfrac * cfrac * ofrac;
                descriptor[(rindex * INDEX_SIZE + cindex) * ORI_SIZE + oindex] += v;
            }
        }
    }
    
    
    void
    InternalCount(Real64* descriptor, Array::Array2dScalarReal64* magnitude,
                  Array::Array2dScalarReal64* direction, 
                  int r, int c, Real64 rpos, Real64 cpos,
                  Real64 rx, Real64 cx, Real64 pointOrientation)
    {
        // bounds check
        if (r < 0  ||  r >= magnitude->CH()  ||  c < 0  ||  c >= magnitude->CW())
           return;
        
        Real64 mag = 0.0;
        if(noGaussianWeighting)
        {
            mag = magnitude->Val(c, r);
        }
        else
        {
            // find Gaussian weight, depending on distance from center
            const Real64 IndexSigma = 1.0;
            Real64 sigma = IndexSigma * 0.5 * INDEX_SIZE;
            Real64 weight = exp(- (rpos * rpos + cpos * cpos) / (2.0 * sigma * sigma));
            mag = weight * magnitude->Val(c, r);
        }
    
        /* Subtract keypoint orientation to give ori relative to keypoint. */
        Real64 ori = direction->Val(c, r) - pointOrientation;
        ori = CanonicalizeAngle(ori);
        PutInVector(descriptor, mag, ori, rx, cx);
    }
    
    
    /* Add features to vec obtained from sampling the grad and ori images
       for a particular scale.  Location of key is (scale,row,col) with respect
       to images at this scale.  We examine each pixel within a circular
       region containing the keypoint, and distribute the gradient for that
       pixel into the appropriate bins of the index array.
    */
    void
    InternalSample(Real64* descriptor, Array::Array2dScalarReal64* magnitude,
                   Array::Array2dScalarReal64* direction, 
                   Real64 scale, Real64 row, Real64 col, Real64 pointOrientation)
    {
       Real64 rpos, cpos, rx, cx;
    
       // integer version of the patch center position (in image coordinates)
       int irow = (int) (row + 0.5);
       int icol = (int) (col + 0.5);
       Real64 sinus = sin(pointOrientation);
       Real64 cosinus = cos(pointOrientation);
              
       /* The spacing of index samples in terms of pixels at this scale. */
       const int magnificationFactor = 3;
       Real64 spacing = scale * magnificationFactor;
    
       /* Radius of index sample region must extend to diagonal corner of
          index patch plus half sample for interpolation.
       */
       Real64 radius = 1.414 * spacing * INDEX_SIZE / 2.0;
       int iradius = (int) (radius + 0.5);
       
       /* Examine all points from the gradient image that could lie within the
          index square.
       */
       for (int i = -iradius; i <= iradius; i++)
         for (int j = -iradius; j <= iradius; j++) {
             
           /* Rotate sample offset to make it relative to key orientation.
              Uses (row,col) instead of (x,y) coords.  Also, make subpixel
              correction as later image offset must be an integer.  Divide
              by spacing to put in index units.
           */
           rpos = ((cosinus * i +   sinus * j) - (row - irow)) / spacing;
           cpos = ((- sinus * i + cosinus * j) - (col - icol)) / spacing;

           // these above are relative to the center of the patch, and also
           // normalized to a "unit distance" (e.g. we can just put in a 
           // gaussian now)
           // simplified (compensates for deviation from integer positions):
           // rpos = (i - (row - irow)) / spacing;
           // cpos = (j - (col - icol)) / spacing;
             
           /* Compute location of sample in terms of real-valued index array
              coordinates.  Subtract 0.5 so that rx of 1.0 means to put full
              weight on index[1] (e.g., when rpos is 0 and INDEX_SIZE is 3.
           */
           rx = rpos + INDEX_SIZE / 2.0 - 0.5;
           cx = cpos + INDEX_SIZE / 2.0 - 0.5;
    
           /* Test whether this sample falls within boundary of index patch. */
           if((rx > -1.0) && (rx < (Real64) INDEX_SIZE) &&
              (cx > -1.0) && (cx < (Real64) INDEX_SIZE))
             InternalCount(descriptor, magnitude, direction, irow + i, icol + j, rpos, cpos,
                           rx, cx, pointOrientation);
         }
    }


    void
    ComputeGradient(Array::Array2dScalarReal64*& magnitude, 
                    Array::Array2dScalarReal64*& direction, 
                    Array::Array2dScalarReal64* im,
                    Real64 sigma)
    {
        using namespace Impala::Core::Array;

        Real64 precision = 3.0;
    
        Array2dScalarReal64* Lx = 0;
        Array2dScalarReal64* Ly = 0;
        RecGauss(Lx, im, sigma, sigma, 1, 0, precision);
        RecGauss(Ly, im, sigma, sigma, 0, 1, precision);
    
        // compute gradient magnitude
        Array2dScalarReal64* Lx2 = 0;
        Array2dScalarReal64* Ly2 = 0;
        Mul(Lx2, Lx, Lx);
        Mul(Ly2, Ly, Ly);
        Add(magnitude, Lx2, Ly2);
        Sqrt(magnitude, magnitude);
        delete Lx2;
        delete Ly2;
        
        // compute gradient direction
        Atan2(direction, Ly, Lx);
        delete Lx;
        delete Ly;
    }


    void 
    SamplePoint(Real64* descriptor, Array::Array2dScalarReal64* magnitude,
                Array::Array2dScalarReal64* direction, 
                Real64 scale, Real64 y, Real64 x, Real64 pointOrientation)
    {
        for(int j = 0; j < VecLength; j++)
        {
            descriptor[j] = 0.0;
        }

        InternalSample(descriptor, magnitude, direction, scale, y, x, pointOrientation);
        
        if(mCustomNorm < 0)
        {
            // Normalize length of vec
            Real64 sqlen = 0.0;
            for (int j = 0; j < VecLength; j++)
                sqlen += (descriptor[j] * descriptor[j]);
    
            Real64 factor = 1.0 / sqrt(sqlen);
            if(sqlen == 0.0) factor = 1.0;
            //std::cout << x << " " << sqlen << " " << factor << std::endl;
    
            const Real64 maxBinValue = 0.2f;
            for (int j = 0; j < VecLength; j++)
            {
                descriptor[j] = descriptor[j] * factor;
                if(descriptor[j] > maxBinValue)
                {
                    descriptor[j] = maxBinValue;
                }
            }
    
            // Normalize length of vec
            sqlen = 0.0;
            for (int j = 0; j < VecLength; j++)
                sqlen += (descriptor[j] * descriptor[j]);
    
            factor = 1.0 / sqrt(sqlen);
            if(sqlen == 0.0) factor = 1.0;
            //std::cout << x << " " << factor << std::endl;
            
            for (int j = 0; j < VecLength; j++)
                descriptor[j] = static_cast<int>(descriptor[j] * factor * mOutputMultiplier);
        }
        else if(mCustomNorm == 0)
        {
            for (int j = 0; j < VecLength; j++)
                descriptor[j] = descriptor[j] * mOutputMultiplier;
        }
        else
        {
            // Normalize length of vec
            Real64 sqlen = 0.0;
            for (int j = 0; j < VecLength; j++)
                sqlen += pow(descriptor[j], mCustomNorm);
    
            Real64 factor = 1.0 / pow(sqlen, 1.0 / mCustomNorm);
            if(sqlen == 0.0) factor = 1.0;

            if(!mNoMaxBins)
            {
                const Real64 maxBinValue = 0.2f;
                for (int j = 0; j < VecLength; j++)
                {
                    descriptor[j] = descriptor[j] * factor;
                    if(descriptor[j] > maxBinValue)
                    {
                        descriptor[j] = maxBinValue;
                    }
                }
        
                // Normalize length of vec
                sqlen = 0.0;
                for (int j = 0; j < VecLength; j++)
                    sqlen += pow(descriptor[j], mCustomNorm);
        
                factor = 1.0 / pow(sqlen, 1.0 / mCustomNorm);
                if(sqlen == 0.0) factor = 1.0;
            }
    
            for (int j = 0; j < VecLength; j++)
                descriptor[j] = descriptor[j] * factor * mOutputMultiplier;
        }
    }


    Matrix::Mat*
    DoCalculateFISTDescriptors(Array::Array2dScalarReal64* im, 
                               PointDescriptorTable* pointData)
    {
        using namespace Impala::Core::Array;

        /* Scales used by Lowe */
        const int clampScaleCount = 18;
        Real64 clampScales[clampScaleCount+1] = { 1.00794,
                                                  1.26992,
                                                  1.6,
                                                  2.01587,
                                                  2.53984,
                                                  3.2,
                                                  4.03175,
                                                  5.07968,
                                                  6.4,
                                                  8.06349,
                                                  10.1594,
                                                  12.8,
                                                  16.127,
                                                  20.3187,
                                                  25.6,
                                                  32.254,
                                                  40.6375,
                                                  51.2,
                                                  -1};
    
        bool warned = false;
        Matrix::Mat* allDescriptors = new Array2dScalarReal64(
                                           VecLength, pointData->Size(), 0, 0);
        for(int i = 0; i < clampScaleCount; i++) {
            Array2dScalarReal64* magnitude = 0; 
            Array2dScalarReal64* direction = 0;
            bool haveSome = false;
            for(int iter = 0; iter < pointData->Size(); iter++)
            {
                Real64 scale = pointData->GetScale(iter);
                if((clampScales[i] > scale) && (i != 0)) continue;
                if(clampScales[i+1] <= scale) continue;
                haveSome = true;
                break;
            }
            if(haveSome)
            {
                ComputeGradient(magnitude, direction, im, clampScales[i]);
                MulVal(direction, direction, -1);

                #pragma omp parallel for private(iter) shared(pointData, allDescriptors, magnitude, direction)
                for(int iter = 0; iter < pointData->Size(); iter++)
                {
                    Real64 scale = pointData->GetScale(iter);
                
                    // check if this is the scale to compute the descriptor at
                    if((clampScales[i] > scale) && (i != 0)) continue;
                    if(clampScales[i+1] <= scale) continue;

                    Real64 x = pointData->GetX(iter);
                    Real64 y = pointData->GetY(iter);
                    Real64 pointOrientation = pointData->GetOrientation(iter);
    
                    SamplePoint(allDescriptors->CPB(0, iter), magnitude, direction, scale, y, x, pointOrientation);
                }
            }
            /* clean up */
            if(magnitude) delete magnitude;
            if(direction) delete direction;
        }
        return allDescriptors;
    }

    
    Matrix::Mat*
    DoCalculateFISTDescriptorsAndEat(Array::Array2dScalarReal64* image, 
                                     PointDescriptorTable* pointData)
    {
        ILOG_VAR(Core.Feature.DoCalculateFISTDescriptorsAndEat);
        Matrix::Mat* result = 0;
        Timer timer;
#ifdef CUDA
        if(Link::Cuda::CudaUsed())
        {
            result = Link::Cuda::DoCalculateSIFTDescriptors(image, pointData);
        }
        else
#endif
        {
            result = DoCalculateFISTDescriptors(image, pointData);
        }
        ILOG_DEBUG("FIST channel in " << timer.SplitTimeStr());
        delete image;
        return result;
    }
};


bool
StringIsFullyNumeric(String s)
{
    for(int i = 0; i < s.size(); i++)
    {
        if((s[i] >= '0') && (s[i] <= '9'))
            continue;
        return false;
    }
    return true;
}

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

    ILOG_VAR(Core.Feature.CalculateFISTDescriptors);
    FISTDescriptor<4, 8> defaultFIST(false);
    Array2dVec3Real64* input2 = 0;
    if (descriptor.find("fist") != String::npos)
    {
        ILOG_INFO("No max bins in norm");
        defaultFIST.mNoMaxBins = true;
        defaultFIST.mCustomNorm = 2;
    }
    if (descriptor.find("g0sift") != String::npos)
    {
        sRGB2RGB(input2, inputNoBorder);
    }
    else if (descriptor.find("g1sift") != String::npos)
    {
        sRGB2RGB(input2, inputNoBorder);
        Pow(input2, input2, 1.0 / 1.8);
    }
    else if (descriptor.find("g2sift") != String::npos)
    {
        sRGB2RGB(input2, inputNoBorder);
        Pow(input2, input2, 1.0 / 2.0);
    }
    else if (descriptor.find("g3sift") != String::npos)
    {
        sRGB2RGB(input2, inputNoBorder);
        MulVal(input2, input2, 1.0 / 255.0);
        Pow(input2, input2, 1.0 / 1.8);
        MulVal(input2, input2, 255.0);
    }
    else if (descriptor.find("g4sift") != String::npos)
    {
        sRGB2RGB(input2, inputNoBorder);
        MulVal(input2, input2, 1.0 / 255.0);
        Pow(input2, input2, 1.0 / 2.0);
        MulVal(input2, input2, 255.0);
    }
    else
    {
        MulVal(input2, inputNoBorder, 1.0);
    }
    descriptor = StringReplace(descriptor, "g0sift", "");
    descriptor = StringReplace(descriptor, "g1sift", "");
    descriptor = StringReplace(descriptor, "g2sift", "");
    descriptor = StringReplace(descriptor, "g3sift", "");
    descriptor = StringReplace(descriptor, "g4sift", "");
    descriptor = StringReplace(descriptor, "sift", "");
    descriptor = StringReplace(descriptor, "fist", "");
    if(descriptor.size() >= 2)
    {
        if((descriptor[descriptor.size() - 2] == 'n') && (StringIsFullyNumeric(descriptor.substr(descriptor.size() - 1, 1))))
        {
            defaultFIST.mCustomNorm = atoi(descriptor.substr(descriptor.size() - 1, 1));
            descriptor = descriptor.substr(0, descriptor.size() - 2);
        }
        if(descriptor.size() >= 3)
        {
            if((descriptor[descriptor.size() - 3] == 'n') && (StringIsFullyNumeric(descriptor.substr(descriptor.size() - 2, 2))))
            {
                defaultFIST.mCustomNorm = atoi(descriptor.substr(descriptor.size() - 2, 2));
                descriptor = descriptor.substr(0, descriptor.size() - 3);
            }
        }
        ILOG_DEBUG("Set custom norm: " << defaultFIST.mCustomNorm);
    }

    /* assumption: inputNoBorder has RGB input image */
    std::vector<Array2dScalarReal64*> channels =
        Core::Feature::GetColorChannels(input2, descriptor);
    delete input2;
    std::vector<Matrix::Mat*> channelDescriptors;
    for(int i=0 ; i<channels.size() ; i++)
    {
        channelDescriptors.push_back(defaultFIST.DoCalculateFISTDescriptorsAndEat(channels[i], pointData));
    }
    
    if(channelDescriptors.size() == 1)
    {
        // takes ownership of memory
        pointData->ReplaceAllDescriptors(channelDescriptors[0]);
    }
    else
    {
        // now turn them all into a big descriptor
        Matrix::Mat* descriptors =
            Matrix::MatConcatenateHorizontal(channelDescriptors);
        
        // takes ownership of memory
        pointData->ReplaceAllDescriptors(descriptors);

        for(int i=0 ; i<channelDescriptors.size() ; i++)
            delete channelDescriptors[i];
    }

}

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

#endif

---