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

---

Octaves.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_Octaves_h
#define Impala_Core_Array_Octaves_h

#include "Core/Array/Set.h"
#include "Core/Array/Sub.h"
#include "Core/Array/Add.h"
#include "Core/Array/LocalMinMax.h"
#include "Core/Array/Scale.h"
#include "Core/Array/Rotate.h"
#include "Core/Array/PixMinMax.h"
#include "Core/Array/Mul.h"
#include "Core/Array/MulVal.h"
#include "Core/Array/AddVal.h"
#include "Core/Array/WritePng.h"
#include "Core/Array/GaussDerivative.h"
#include "Core/Array/RecGauss.h"
#include "Core/Array/MakeRoi.h"
#include "Core/Array/MakeGaussian2d.h"
#include "Core/Array/Array2dTem.h"
#include "Core/Array/PrintData.h"

#include "Core/Feature/FeatureTable.h"

#include "Core/Histogram/Histogram1dTem.h"
#include "Core/Histogram/MakeHistogram1d.h"

#include "Core/Array/RecGauss.h"
#include "Core/Matrix/MatMulVec.h"

//#include "Core/Feature/SiftFeatureTable.h"
//
#include "Core/Table/TableTem.h"
#include "Core/Column/ColumnTem.h"
#include <math.h>

#ifndef PI
#define PI 3.14159265358979324
#endif

//using namespace Impala::Core::Matrix;
//using namespace Impala::Core::Feature;

namespace Impala{
namespace Core{
namespace Array{
//This class represents the multiresolution octaves of the given array.
//
//Given a number of octaves, and number of scales at each octave, with an
//initial sigma value, It calculates the necessary scalings of the sigma 
//parameter to be used for gaussian blurring
//
//
class Octaves{
public:

private:
    ILOG_VAR_DEC;
    Impala::Real64 mSigma0;
    Impala::Real64 mSigmaN;
    int mOctaves;
    int mMinOctave;
    int mLevels;
    int mScales;
    double * mSigmas;
    double mSigmaK;
    Array2dVec3Real64* mSrcVec3;
    Array2dScalarReal64* mSrcScalar;
    Array2dScalarReal64** mOctaveArray;
    Array2dScalarReal64** mSrcArray;


    bool mRecGauss;
    bool mResample;
    bool mRecursive;

public:
    Octaves(Array2dScalarUInt8* src, int Octaves, int Levels, int MinOctave=0)
    : mOctaves(Octaves),
    mMinOctave(MinOctave),
    mLevels(Levels),
    mRecGauss(true),
    mRecursive(true),
    mSrcScalar(0)
    {
        CmdOptions& options = CmdOptions::GetInstance();
        mRecGauss=options.GetBool("recGauss",false);
        mResample=options.GetBool("resample",false);
        mRecursive=options.GetBool("recursive",false);


        ILOG_INFO("Using RecGauss  : "<<(mRecGauss?"yes":"no"));
        ILOG_INFO("Using Resampling: "<<(mResample?"yes":"no"));
        ILOG_INFO("All Recursive   : "<<(mRecursive?"yes":"no"));
        //Copy the input array to our Real64 array
        Real64 min,max;
        Set(mSrcScalar,src);
        PixMinMax(src,&min,&max);
        ILOG_INFO("Input image Min-Max:"<<min<<" "<<max);
        Normalize(mSrcScalar);
        PixMinMax(mSrcScalar,&min,&max);
        ILOG_INFO("Input image Min-Max:"<<min<<" "<<max);
        Init();
    }
    ~Octaves(){
        delete [] mSigmas;
        
        for(int i=0;i<mScales*mOctaves;i++)
            delete mOctaveArray[i];
        delete [] mOctaveArray;

        for(int i=0;i<mOctaves;i++)
            delete mSrcArray[i];
        delete [] mSrcArray;


    }
    void Normalize(Array2dScalarReal64* s){
        Real64 min,max;
        PixMinMax(s,&min,&max);
        AddVal(s,s,-min);
        MulVal(s,s,1/(max-min));
    }

    int GetOctaveCount(){
        return mOctaves;
    }
    
    int GetLevelCount(){
        return mLevels;
    }

    int GetScaleCount(){
        return mScales;
    }

    void Init(){
        //If a negative count is given for the octaves, we calculate the max
        //possible, deducing it from the image size. A minimum size of 8x8 
        //array is produced from the original array
        if(mOctaves<1)
        {
            mOctaves = Max(1,
            int(std::floor(log2(Min(mSrcScalar->CW(),mSrcScalar->CH()))) 
                           - mMinOctave - 3));
        }
        //We need levels+3 in each octave, since they will first be subtracted,
        //and then the subtration will be used for extrema localization
        //within +/- 1 scales
        mScales=mLevels+3;

        ILOG_INFO("Initializing "<<mOctaves
                <<" octaves, each with "<<mLevels<<" levels and "
                <<mScales<<" scales.");
        ILOG_INFO((mResample?"Resampling Enabled":"Resampling Disabled"));

        
        mSigmas = new double[mOctaves*mScales];
        mOctaveArray = new  Array2dScalarReal64*[mOctaves*mScales];
        for(int i=0;i<mOctaves*mScales;i++)
            mOctaveArray[i]=0;

        //mSrcArray holds the images that will be gaussian smoothed,
        //for first octave it's the original image
        //for octaves below, it is calculated from the gaussian smoothing with
        //2*sigma
        mSrcArray = new  Array2dScalarReal64*[mOctaves];
        for(int i=0;i<mOctaves;i++)
            mSrcArray[i]=0;

        
        ILOG_INFO("Octave structures constructed!");
        BuildOctave();
        ILOG_INFO("Octave structures initialized!");

        
    }

    void BuildOctave(){
        //Calculate sigmas for the first levels of each octave
        //Each sigma is double its counterpart in the octave below
        //
        //This was 0.5 before, for prior smoothing
        mSigmas[0] = 1;

        for(int o=1;o<mOctaves;o++)
        {
            mSigmas[o*mScales] = 2 * mSigmas[(o-1)*mScales];
        }

        //mSigmaK is the ratio of sigmas between successive levels
        mSigmaK = powf(2.0,1.0/mLevels);

        //Calculate the sigma values of other levels.
        //Depending on the number of levels in each octave, we divide the sigmas
        //logarithmically, each scaled by 2^(1/mLevels)
        for(int o=0;o<mOctaves;o++){
            std::ostringstream oss;
            oss<<"Octave "<<o<<" :";
            for(int l=0;l<mScales;l++)
            {
                mSigmas[o*mScales+l] = mSigmas[o*mScales] * powf(mSigmaK,l);
                oss<<" "<< mSigmas[o*mScales+l];
            }
            ILOG_INFO(oss.str());
        }

        //Calculate the gaussian smoothed versions of the image, 
        //given the sigma values for each level
        for(int o=0;o<mOctaves;o++)
        {
            //Calculate the initial image, 

            //If it's not the first octave, the image is resampled from 2*sigma
            //of above octaves first level.
            if(mResample)
            {
                if((o==0))
                {
                    //The first array is the linear upsampling of the original
                    mSrcArray[o] = 0;
                    Scale(mSrcArray[o], mSrcScalar,2,2,Geometry::LINEAR);

                }else{
                    //resample the scale that has double the sigma of 
                    //previous octave's first scale
                    Array2dScalarReal64* src=mOctaveArray[(o-1)*mScales+mLevels];
                    mSrcArray[o] = 0;
                    Scale(mSrcArray[o],src,.5,.5,Geometry::NEAREST);
                }
            }
            else
            {
                //Don't resample the image from source, use directly
                //This causes all the analysis to be done in the same resolution
                if(o==0){
                    mSrcArray[o] = 0;
                    Set(mSrcArray[o], mSrcScalar);
                }else{
                    mSrcArray[o] = 0;
                    Set(mSrcArray[o],mOctaveArray[(o-1)*mScales+mLevels]);
                }

            }
            //The source image for this octave is ready
            //Do any prior smoothing
            

            for(int l=0;l<mScales;l++)
            {
                //ILOG_DEBUG("Octave "<<o<<" Level "<<l<<" Src:"<<mSrcArray[o]);
                if(mRecGauss)
                {
                    if(mRecursive)
                    {
                        //Calculate each gaussian smoothing recursively
                        Array2dScalarReal64* src;
                        //If it's the first level in the octave, 
                        //we should use the array we calculated in the previous
                        //step, else we use the previous level in the octave
                        if(l==0)
                            src = mSrcArray[o];
                        else
                            src = mOctaveArray[o*mScales+l-1];

                        RecGauss(mOctaveArray[o*mScales+l],
                                 src,
                                 mSigmaK,
                                 mSigmaK,0,0,3);
                    }else{
                        //Calculate each gaussian smoothing directly using
                        //recursive gauss
                        RecGauss(mOctaveArray[o*mScales+l],
                             mSrcArray[o],
                             mSigmas[l],
                             mSigmas[l],
                             0,0,3);

                    }
                }
                else
                {
                    if(mRecursive)
                    {
                        Array2dScalarReal64* src; 
                        if(l==0)
                            src = mSrcArray[o];
                        else
                            src = mOctaveArray[o*mScales+l-1];
                        
                        GaussDerivative(mOctaveArray[o*mScales+l],
                             src,
                             mSigmaK,0,0,3);

                    }else{
                        GaussDerivative(mOctaveArray[o*mScales+l],
                             mSrcArray[o],
                             mSigmas[l],
                             0,0,3);
                    }
                }
                Real64 min=0;
                Real64 max=0;
                PixMinMax(mOctaveArray[o*mScales+l],&min,&max);
                //ILOG_DEBUG("Gaussian @ "<<o<<"x"<<l<<" : ("<<min<<","<<max<<")");
            }
        }
    }

    Array2dScalarReal64* GetOctave(int o,int l)
    {
        return mOctaveArray[o*mScales+l];
    }
    

    double GetSigma(int o, int l){
        return mSigmas[o*mScales+l];
    }
    bool IsResamplingEnabled(){
        return mResample;    
    }

};
ILOG_VAR_INIT(Octaves, Impala.Core.Array);

}}}//Namespace



#endif

---