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

---

KfrMotionExtractor.h

Go to the documentation of this file.

#ifndef Impala_Core_VideoSet_KfrMotionExtractor_h
#define Impala_Core_VideoSet_KfrMotionExtractor_h

// created by Jun Wu, on Dec 18, 2009
//
// purpose: extract motion feature before and after keyframes.
// Note: 
// 1) we maintaince a window, the center frame is key-frame
// 2) within it, we also have a sliding buffer, to compute the energy distribution function
// 3) considering all siding buffer as a sequence, we use the average as the final output
//      window size: 7, 9, 11, 13, ...
//      default buffer size: = 5
// 4) arguments: requires --srcWindow "size;9;rgb2ooo"
// 5) to represent dominant motion (horizontal,vertical,radial), we use three motion filters

#include "Core/VideoSet/Reporter.h"
#include "Core/Vector/VectorTem.h"
#include "Persistency/FeatureTableRepository.h"
#include "Core/Stream/SeqConvKernel.h"
#include "Core/Array/RGB2Gray.h"

// comment this line when commit to impala
//#define __BITMAP_FUNC__

#ifdef __BITMAP_FUNC__
#include "Core/Feature/Bitmap.h"
#endif

// motion estimation only supports 16x16 block.
// Note: do NOT change marco MVBLOCKSIZE
#define MVBLOCKSIZE   16 

#define MAX_DISP    10
#define TPL_NUM      3

#ifndef BYTE
#define BYTE unsigned char
#endif

#ifndef BOOL
#define BOOL bool
#endif

#ifndef UINT
#define UINT unsigned int
#endif


namespace Impala
{
namespace Core
{
namespace VideoSet
{

class KfrMotionExtractor : public Listener
{

protected:

    typedef Vector::VectorTem<Real64> VectorReal64;

public:

    KfrMotionExtractor(Reporter* reporter, 
        CmdOptions& options, CString featureName, int binCount)
    {
        mReporter = reporter;
        mFeatureName = featureName;
        mBinCount = binCount;
        Feature::FeatureDefinition def(GetFeatureName());
        mHistogramTable = new Feature::FeatureTable(def, 1000, GetBinCount());

        m_bNormalize = TRUE;
        m_verbose = FALSE;
        m_verbose2 = FALSE;
    }

    KfrMotionExtractor(Reporter* reporter, CmdOptions& options)
    {
        mReporter = reporter;
        mFeatureName = "kfrmotion";
        mBinCount = 2 * TPL_NUM;
        Feature::FeatureDefinition def(GetFeatureName());
        mHistogramTable = new Feature::FeatureTable(def, 1000, GetBinCount());

        m_bNormalize = TRUE;

        m_verbose = FALSE;
        m_verbose2 = FALSE;
        //m_verbose = TRUE;
        //m_verbose2 = TRUE;
        
    }

    virtual 
    ~KfrMotionExtractor()
    {
    }

    String
    GetFeatureName() const
    {
        return mFeatureName;
    }

    virtual void
    HandleNewFrame(VideoSet* vs, int fileId, Stream::RgbDataSrc* src)
    {
        if (!src)
        {
            ILOG_ERROR("Not a valid source for fileId " << fileId);
            return;
        }

        Array::Array2dVec3UInt8* im = Array::ArrayCreate<Array::Array2dVec3UInt8>
            (src->FrameWidth(), src->FrameHeight(), 0, 0, src->DataPtr(), true);
        if (!im)
        {
            ILOG_ERROR("Couldn't load image for fileId " << fileId);
            return;
        }

        if (m_verbose)
        {
            std::cout << "Current Frame No. = " << src->FrameNr() << std::endl;
        }

#ifdef __BITMAP_FUNC__
        Core::Feature::Bitmap bitmap;
        bitmap.SaveRgb2BitmapFile("test_rgb_00.bmp",src->DataPtr(), src->FrameWidth(), src->FrameHeight() );
#endif

        /*VectorReal64 histogram = ComputeHistogram(im);
        Quid quid = vs->GetQuidFrame(fileId, src->FrameNr());
        mHistogramTable->Add(quid, histogram);*/

        double feature[2*TPL_NUM];
        memset(feature,0,sizeof(double)*2*TPL_NUM);

        // For temporal filtering use something like this instead:
        // (requires --srcWindow "size;5;rgb2ooo")
        typedef Stream::RgbDataSrcWindow<Stream::WindowPrepRgb2Ooo> WindowType;
        //WindowType* window = dynamic_cast<WindowType*>(src);
        WindowType* window = (WindowType*)src;
        if (window)
        {
            if (m_verbose)
            {
                std::cout << "window size = " << window->WindowSize() << std::endl;
            }

            //Array::Array2dScalarReal64* tGauss = 0;
            //double temporalSigma = 0.75;
            //tGauss = Array::MakeGaussian1d(temporalSigma, 0, 3.0, 100);
            //Array::Array2dVec3Real64* res = 0;
            //SeqConvKernel(res, window, tGauss);
            //Array::Array2dVec3UInt8* src;
            //Array::InvWiccest(res, src, false, false, 0.75, 3, true);


            // test looping for all frames within the window
            //for (int i=0 ; i<window->WindowSize(); i++)
            //{
            //    Array::Array2dVec3UInt8* cur = Array::ArrayCreate<Array::Array2dVec3UInt8>
            //        (window->FrameWidth(), window->FrameHeight(), 0, 0, window->DataPtrWindow(i), true);
            //    if (!cur)
            //    {
            //        ILOG_ERROR("Couldn't load image for fileId " << fileId);
            //        return;
            //    }

            //    VectorReal64 cur_histogram = ComputeHistogram(cur);
            //    Quid quid = vs->GetQuidFrame(fileId, window->FrameNr() - window->WindowSize()/2 + i);
            //    mHistogramTable->Add(quid, cur_histogram);
            //}
            

            int bWidth = window->FrameWidth()/MVBLOCKSIZE;
            int bHeight = window->FrameHeight()/MVBLOCKSIZE;
            
            // Note: we have tested these pairs: (10,7) (7,10) (10,8) (8,10) (22,18) (20,15)
            //int bWidth = 10;
            //int bHeight = 7;

            for (int t=0;t<TPL_NUM;t++)
            {
                m_Weight[t] = new double[bWidth * bHeight];
                memset(m_Weight[t], 0, sizeof(double)* bWidth * bHeight);

                if (m_verbose2)
                {
                    DumpWeightMatrix(m_Weight[t], bWidth, bHeight, t);
                }
            }

            m_EnergyBuffer = new BYTE[bWidth * bHeight];
            m_EnergyMap = new double[bWidth * bHeight];
            m_Buffer = new BYTE[bWidth * bHeight * 2];

            // Initialization
            LoadEnergyTemplate(bWidth, bHeight);

            // calculate the mean of weight for all templates
            CalcTPLMean(bWidth, bHeight);

            m_MVseq = new BYTE[bWidth * bHeight * 2 * (window->WindowSize()-1)];
            memset(m_MVseq, 0, sizeof(BYTE)*bWidth * bHeight * 2 * (window->WindowSize()-1) );
            //DumpMVFseq(m_MVseq, bWidth, bHeight, window->WindowSize()-1);

            // compute motion features between adjacent frames
            // if window size as 9: 0-1,1-2,2-3,3-4,4-5,5-6,6-7,7-8,8-9
            // based on these motion vectors field, we compute the energy distribution
            for (int i=1 ; i<window->WindowSize() ; i++)
            {
                // previous frame: pre
                Array::Array2dVec3UInt8* pre_rgb = Array::ArrayCreate<Array::Array2dVec3UInt8>
                    (window->FrameWidth(), window->FrameHeight(), 0, 0, window->DataPtrWindow(i-1), true);
                if (!pre_rgb)
                {
                    ILOG_ERROR("Couldn't load image for fileId " << fileId);
                    return;
                }

                // current frame: cur
                Array::Array2dVec3UInt8* cur_rgb = Array::ArrayCreate<Array::Array2dVec3UInt8>
                    (window->FrameWidth(), window->FrameHeight(), 0, 0, window->DataPtrWindow(i), true);
                if (!cur_rgb)
                {
                    ILOG_ERROR("Couldn't load image for fileId " << fileId);
                    return;
                }

#ifdef __BITMAP_FUNC__
                //Core::Feature::Bitmap bitmap;
                //bitmap.SaveRgb2BitmapFile("test_rgb_0.bmp",cur_rgb->CPB(), cur_rgb->CW(), cur_rgb->CH() );
#endif
                //Array2dScalarReal64* pre = 0;
                //Array2dScalarReal64* cur = 0;
                Array::Array2dScalarUInt8* pre = 0;
                Array::Array2dScalarUInt8* cur = 0;
                RGB2Gray(pre,pre_rgb);
                RGB2Gray(cur,cur_rgb);

                // motion vector computing between (pre,cur)
                Int64 body = pre->BW() + pre->BH();

#ifdef __BITMAP_FUNC__
                //Core::Feature::Bitmap bitmap;
                //bitmap.SaveGray2BitmapFile("test_gray_0.bmp",cur->CPB(), cur->CW(), cur->CH() );
#endif
                CalcMotionVectors((unsigned char*)pre->CPB(),(unsigned char*)cur->CPB(),window->FrameWidth(),window->FrameHeight(),m_MVseq, i-1);

                if (m_verbose2) DumpMVF(m_MVseq, bWidth, bHeight, i-1);

                // DO remember: release the memory
                delete pre_rgb;
                delete cur_rgb;
                delete pre;
                delete cur;

            }
            //DumpMVFseq(m_MVseq, bWidth, bHeight, window->WindowSize()-1);

            // -----------------------------------------------------
            // #frame within the window:  window->WindowSize()
            // Motion Vector Field:       window->WindowSize() - 1
            // Energy Distribution Curve: window->WindowSize() - 2
            // -----------------------------------------------------

            // compute energy based on weighted templates
            // using a sliding buffer within a window
            int result = ComputeEnergy(window->WindowSize()-1, bWidth, bHeight, m_MVseq, feature);

            if (0 !=  result) return;

            for (int i=0;i<TPL_NUM;i++)
                delete []m_Weight[i];
            delete []m_EnergyBuffer;
            delete []m_EnergyMap;
            delete []m_Buffer;
            delete []m_MVseq;
     
        }
        delete im;

        VectorReal64 mvfeature(2*TPL_NUM);
        for (int x=0; x<2*TPL_NUM; x++)
        {
            mvfeature[x] = feature[x];
        }
        Quid quid = vs->GetQuidFrame(fileId, src->FrameNr());
        mHistogramTable->Add(quid, mvfeature);

        ILOG_DEBUG("... done keyframe " << src->FrameNr());

    }

    virtual void
    HandleDoneFile(VideoSet* vs, int fileId, Stream::RgbDataSrc* src)
    {
        bool binary = true;
        Feature::FeatureDefinition def = mHistogramTable->GetFeatureDefinition();
        typedef Persistency::FeatureLocator FeatureLocator;
        FeatureLocator loc(vs->GetLocator(), false, false, "Keyframes",
                           def.AsString(), vs->GetContainerFile(fileId));
        Persistency::FeatureTableRepository().Add(loc, mHistogramTable);
        mHistogramTable->SetSize(0);
    }

    // output all the MVF (Motion Vector Field) in the sequence
    void DumpMVFseq(BYTE* mvseq, int bWidth, int bHeight, int N)
    {
        //memset(mvseq, 0, sizeof(BYTE)*bWidth * bHeight * 2 * (window->WindowSize()-1) );

        for (int w=0; w<N; w++)
        {
            printf("MVF (%dx%d) for window %d:\n",bWidth,bHeight,w);
            for (int j=0; j<bHeight; j++)
            {
                for (int i=0; i<bWidth; i++)
                {
                    BYTE sx = mvseq[sizeof(BYTE)*bWidth*bHeight*2* w + (j*bWidth + i)*2 +0];
                    BYTE sy = mvseq[sizeof(BYTE)*bWidth*bHeight*2* w + (j*bWidth + i)*2 +1];
                    printf("(%d,%d) ",sx,sy);
                }
                printf("\n");
            }
            printf("\n");
        }
    }

    // output the index-th MVF (Motion Vector Field) in the sequence
    void DumpMVF(BYTE* mvseq, int bWidth, int bHeight, int index)
    {
        printf("MVF (%dx%d) for window %d:\n",bWidth,bHeight,index);
        for (int j=0; j<bHeight; j++)
        {
            for (int i=0; i<bWidth; i++)
            {
                BYTE sx = mvseq[sizeof(BYTE)*bWidth*bHeight*2* index + (j*bWidth + i)*2 +0];
                BYTE sy = mvseq[sizeof(BYTE)*bWidth*bHeight*2* index + (j*bWidth + i)*2 +1];
                printf("(%d,%d) ",sx,sy);
            }
            printf("\n");
        }
        printf("\n");
    }

    void DumpEnergyMap(int bWidth, int bHeight, char* comments="")
    {
        printf("Energy Map (%dx%d): Window Index -> %s\n", bWidth, bHeight, comments);
        for (int j=0; j<bHeight; j++)
        {
            for (int i=0; i<bWidth; i++)
            {
                printf("%.1f ",m_EnergyMap[j*bWidth + i]);
            }
            printf("\n");
        }
        printf("\n");
    }

    void DumpEnergyBuffer(int bWidth, int bHeight)
    {
        printf("Energy Buffer (%dx%d):\n",bWidth,bHeight);
        for (int j=0; j<bHeight; j++)
        {
            for (int i=0; i<bWidth; i++)
            {
                printf("%d ",m_EnergyBuffer[j*bWidth + i]);
            }
            printf("\n");
        }
        printf("\n");
    }

    //
    // compute the energy distribution map based on sliding buffer,
    // CurBufSize: = 5 (default)
    // WindowLen > 5: such as 7,9,11,13,... 
    int ComputeEnergy(int WindowLen, int bWidth, int bHeight, BYTE* mvseq, double* feature, int CurBufSize = 5)
    {
        IsInShot = FALSE;

        // ZeroMemory for m_EnergyMap[] and m_EnergyBuffer[]
        OnShotStart(bWidth, bHeight);

        for (int t=0; t<2*TPL_NUM; t++)
        {
            feature[t] = 0.0;
        }

        //int CurBufSize = 5;
        if (WindowLen < CurBufSize)
        {
            printf("Error: WindowLen (%d) is smaller than CurBufSize (%d)\n", WindowLen, CurBufSize);
            return 1;
        }
        int num = WindowLen - CurBufSize;
        for (int index = 0; index < 1+num; index++)
        {
            // Just the fill the buffer according to the motion vector information
            FillBuffer(CurBufSize, bWidth, bHeight, mvseq + sizeof(BYTE)*bWidth*bHeight*2*index);

            int i=0;
            float Temp1Energy[TPL_NUM];
            float Temp2Energy[TPL_NUM];
            float diffEnergy[TPL_NUM];
            memset(diffEnergy,0,TPL_NUM*sizeof(float));

            for (i=0;i<TPL_NUM;i++)
            {
                // Use the weight template to calculate the total energy value
                float Energy = 0.0;
                if (m_bNormalize)
                {
                    Energy = CalcEnergy(i, bWidth, bHeight);
                } else {
                    Energy = CalcEnergy(i, bWidth, bHeight, FALSE);
                }

                // update Temp1Energy[] and Temp2Energy[] for late compute their diffs
                if (index==0)
                    Temp1Energy[i]=Energy;
                else if (index==1)
                    Temp2Energy[i]=Energy;
                else //index>1
                {
                    Temp1Energy[i] = Temp2Energy[i];
                    Temp2Energy[i] = Energy;
                }

                if (m_verbose) printf("%.5f\t", Energy);

                //sum up all the values, to compute averages later.
                feature[i] += Energy;
            }
            //printf("\n");

            // Here there are TPL_NUM pairs of data, we get the diffs of them
            // first diff energy has been initialized as zero
            if (index>0)
            {
                for (i=0;i<TPL_NUM;i++)
                {
                    if (m_bNormalize)
                    {
                        // [-1,1] --> [0,1]
                        diffEnergy[i] = fabs(Temp2Energy[i] - Temp1Energy[i]);
                        //diffEnergy[i] = Normalization(diffEnergy[i],-1.0,1.0);
                    } else {
                        diffEnergy[i] = Temp2Energy[i] - Temp1Energy[i];
                    }

                    if (m_verbose) printf("%.5f\t", diffEnergy[i]);

                    //sum up all the values, to compute averages later.
                    feature[TPL_NUM+i] += diffEnergy[i];

                }
                if (m_verbose) printf("\n");

            } else {

                //index==0
                if (m_verbose) printf("diff1\tdiff2\tdiff3\n");

                for (i=0;i<TPL_NUM;i++)
                {
                    //if (m_verbose) printf("%.5f\t", diffEnergy[i]);
                    feature[TPL_NUM+i] = 0.0;
                }
            }
        }
        
        OnShotEnd(bWidth, bHeight);

        // calculate the averages of the feature within current Window buffer
        for (int t=0; t<TPL_NUM; t++)
        {
            feature[t] /= num+1;
            feature[TPL_NUM+t] /= num; // one less than three original values
        }

        if (m_verbose)
        {
            printf("\nThe Average as the final feature:\n");
            for (int t=0; t<2*TPL_NUM; t++)
            {
                printf("%.5f ",feature[t]);
            }
            printf("\n\n");
        }

        return 0;

    }

    // template: m_Weight[0] - m_Weight[2]:
    // 
    //Template Weigth Matrix (10x7)_0 is:
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    // 1  5  9 13 17 21 25 29 33 37
    //
    //Template Weigth Matrix (10x7)_1 is:
    // 1  1  1  1  1  1  1  1  1  1
    // 6  6  6  6  6  6  6  6  6  6
    //11 11 11 11 11 11 11 11 11 11
    //16 16 16 16 16 16 16 16 16 16
    //21 21 21 21 21 21 21 21 21 21
    //26 26 26 26 26 26 26 26 26 26
    //31 31 31 31 31 31 31 31 31 31
    //
    //Template Weigth Matrix (10x7)_2 is:
    // 1  1  1  1  1  1  1  1  1  1
    // 1 17 17 17 17 17 17 17 17  1
    // 1 17 33 33 33 33 33 33 17  1
    // 1 17 33 49 49 49 49 33 17  1
    // 1 17 33 33 33 33 33 33 17  1
    // 1 17 17 17 17 17 17 17 17  1
    // 1  1  1  1  1  1  1  1  1  1

    // Till now, we have verified pairs for (bWidth,bHeight) as follows:
    //
    //    (10,7) (7,10) (10,8) (8,10) (22,18) (20,15)
    //

    // initialize template buffer, and record the (min,max) values
    void LoadEnergyTemplate(int bWidth, int bHeight)
    {
        int i=0,j=0;

        // pan: horizontal template
        int base = 1;
        int delta = 4;
        for (i=0; i<bWidth; i++)
            for (j=0; j<bHeight; j++)
            {
                m_Weight[0][j*bWidth+i] = base + delta*i;
            }
        m_TPLMin[0] = base;
        m_TPLMax[0] = base + delta*(bWidth-1);

        // tilt: vertical tempate
        delta = 5;
        for (j=0; j<bHeight; j++)   
            for (i=0; i<bWidth; i++)
            {
                m_Weight[1][j*bWidth+i] = base + delta*j;
            }
        m_TPLMin[1] = base;
        m_TPLMax[1] = base + delta*(bHeight-1);

        // zoom: radial tempate
        delta=16;
        int k=0;
        int dim = (bWidth<=bHeight) ? bWidth: bHeight;
        for (k=0; k<dim/2; k++)
        {
            // first column: left
            i=0+k;
            for (j=0+k; j<bHeight-k; j++)
            {
                m_Weight[2][j*bWidth+i] = base + delta*k;
            }

            // last colum: right
            i=bWidth-1-k;
            for (j=0+k; j<bHeight-k; j++)
            {
                m_Weight[2][j*bWidth+i] = base + delta*k;
            }

            // first row: top
            j=0+k;
            for (i=0+k; i<bWidth-k; i++)
            {
                m_Weight[2][j*bWidth+i] = base + delta*k;
            }

            // last row: bottom
            j=bHeight-1-k;
            for (i=0+k; i<bWidth-k; i++)
            {
                m_Weight[2][j*bWidth+i] = base + delta*k;
            }

            //DumpWeightMatrix(m_Weight[2], bWidth, bHeight, 2);

        }

        // if dim is odd, one row or one colum will NOT be process.
        // here, we will do this processing for the left one.
        m_TPLMin[2] = base;
        m_TPLMax[2] = base + delta*k;
        if (dim%2==1)
        {
            if (bWidth >= bHeight)
            {
                // middle row
                j=0+k;
                for (i=0+k; i<bWidth-k; i++)
                {
                    m_Weight[2][j*bWidth+i] = base + delta*k;
                }
            }
            else
            {
                // middle column
                i=0+k;
                for (j=0+k; j<bHeight-k; j++)
                {
                    m_Weight[2][j*bWidth+i] = base + delta*k;
                }
            }

            //DumpWeightMatrix(m_Weight[2], bWidth, bHeight, 2);
        }

        // output the weight matrix
        if (m_verbose2)
        {
            for (int t=0;t<TPL_NUM;t++)
            {
                DumpWeightMatrix(m_Weight[t], bWidth, bHeight, t);
            }
        }

    }

    void DumpWeightMatrix(double* weight, int bWidth, int bHeight, int index)
    {
        // output the weight matrix
        printf("Template Weigth Matrix (%dx%d)_%d is:\n",bWidth,bHeight, index);

        int i,j;
        for (j=0; j<bHeight; j++)
        {
            for (i=0; i<bWidth; i++)
            {
                printf("%2.d ",(int)weight[j*bWidth+i]);
                //printf("%.2f ",weight[j*bWidth+i]);
            }
            printf("\n");
        }
        printf("(Min, Max) = (%.2f, %.2f)\n", m_TPLMin[index], m_TPLMax[index]);
        printf("\n");
    }

    void CalcTPLMean(int bWidth, int bHeight)
    {
        for (int indexOfTPL=0;indexOfTPL<TPL_NUM;indexOfTPL++)
        {
            m_TPLMean[indexOfTPL] = 0;
            for (int i = 0; i < bWidth * bHeight; i ++)
            m_TPLMean[indexOfTPL] += m_Weight[indexOfTPL][i];
        }

        if (m_verbose2)
        {
            printf("Template Mean Values TPL[0...2]: ");
            for (int indexOfTPL=0;indexOfTPL<TPL_NUM;indexOfTPL++)
            {
                printf("%.2f ",m_TPLMean[indexOfTPL]);
            }
            printf("\n\n");
        }

    }

    void OnShotStart(int bWidth, int bHeight)
    {
        if (!IsInShot)
        {
            for (int i = 0; i < bWidth * bHeight; i++)
                m_EnergyMap[i] = 1.0;
            memset(m_EnergyBuffer, 255, bWidth * bHeight);

            IsInShot = TRUE;
        }

        //DumpEnergyBuffer(bWidth, bHeight);

    }

    void OnShotEnd(int bWidth, int bHeight)
    {
        if (IsInShot)
        {
            IsInShot = FALSE;
            memset(m_EnergyBuffer, 0, bWidth * bHeight);
        }

    }


    // fill and update the energy map based on MVF sequence
    void FillBuffer(int WindowLen, int bWidth, int bHeight, BYTE* mvseq)
    {
        double * Add = new double[bWidth * bHeight];
        if (IsInShot)
        {
            UINT BufferLen = bWidth * bHeight;
            double * Temp = new double[BufferLen];
            if (Temp == NULL || Add == NULL)
            {
                if (Temp != NULL) delete[] Temp;
                if (Add != NULL) delete[] Add;
                printf("Not Enough Memory! No Map Generated!");
                return;
            }

            // Initialize the map value to 1.0
            for (int i = 0; i < BufferLen; i++)
            {
                m_EnergyMap[i] = 1.0;
                Add[i] = 0.0;
            }

            //DumpEnergyMap(bWidth, bHeight, "Initialization");

            for (int w = 0; w < WindowLen; w ++)
            {
                memset(Temp, 0, BufferLen * sizeof(double));

                // read the motion vectors from the "mvseq" buffer
                memcpy(m_Buffer, mvseq + sizeof(BYTE)*bWidth*bHeight*2* w, BufferLen * 2);

                // loop for the energy map :(20*15, for example)
                for (int i = 0; i < bWidth; i++)
                    for (int j = 0; j < bHeight; j++)
                    {
                        double CurValue = m_EnergyMap[j * bWidth + i];
                        
                        int cx, cy;
                        long Ptr = 2 * (j * bWidth + i);

                        // (cx,cy) is the motion vector
                        cx = (char)m_Buffer[Ptr];
                        cy = (char)m_Buffer[Ptr + 1];    
                        
                        char SignX, SignY;
                        SignX = (cx >= 0)? 1:-1;
                        SignY = (cy >= 0)? 1:-1;
                        
                        int BaseX,BaseY;
                        BaseX = i + (int)(cx / MVBLOCKSIZE);
                        BaseY = j + (int)(cy / MVBLOCKSIZE);
                        
                        double Factor;
                        if (BaseX >= bWidth || BaseY >= bHeight || BaseX < 0 || BaseY < 0) 
                        {
                            if (abs(cx) >= abs(cy))
                            {
                                Factor = (double)(cy) / (double)(cx);
                                cx = ((SignX > 0)?bWidth:0) * MVBLOCKSIZE - MVBLOCKSIZE / 2;
                                cy = Factor * cx + j * MVBLOCKSIZE - Factor * i * MVBLOCKSIZE;
                                if (cy <= -MVBLOCKSIZE || cy >= bHeight * MVBLOCKSIZE)
                                {
                                    cy = ((SignY > 0)?bHeight:0) * MVBLOCKSIZE - MVBLOCKSIZE / 2;
                                    cx = (cy - j * MVBLOCKSIZE + Factor * i * MVBLOCKSIZE) / Factor;
                                }
                            }
                            else
                            {
                                Factor = (double)(cx) / (double)(cy);
                                cy = ((SignY > 0)?bHeight:0) * MVBLOCKSIZE - MVBLOCKSIZE / 2;
                                cx = Factor * cy + i * MVBLOCKSIZE - Factor * j * MVBLOCKSIZE;
                                if (cx <= - MVBLOCKSIZE || cx >= bWidth * MVBLOCKSIZE)
                                {
                                    cx = ((SignX > 0)?bWidth:0) * MVBLOCKSIZE - MVBLOCKSIZE / 2;
                                    cy = (cx - i * MVBLOCKSIZE + Factor * j * MVBLOCKSIZE) / Factor;
                                }
                            }
                            cx -= i * MVBLOCKSIZE;
                            cy -= j * MVBLOCKSIZE;
                            BaseX = i + (int)(cx / MVBLOCKSIZE);
                            BaseY = j + (int)(cy / MVBLOCKSIZE);
                        }
                        
                        double mX, mY;
                        mX = (double)(abs(cx % MVBLOCKSIZE))/MVBLOCKSIZE;
                        mY = (double)(abs(cy % MVBLOCKSIZE))/MVBLOCKSIZE;
                        if (BaseX + SignX >= bWidth || BaseX + SignX < 0) mX = 0.0;
                        if (BaseY + SignY >= bHeight || BaseY + SignY < 0) mY = 0.0;
                        
                        Temp[BaseY * bWidth + BaseX] += (1 - mX) * (1 - mY) * CurValue;
                        
                        if (BaseX + SignX < bWidth && BaseX + SignX >= 0)
                            Temp[BaseY * bWidth + BaseX + SignX] += mX * (1 - mY) * CurValue;
                        if (BaseY + SignY < bHeight && BaseY + SignY >= 0)
                        {
                            Temp[(BaseY + SignY) * bWidth + BaseX] += (1 - mX) * mY * CurValue;
                            if (BaseX + SignX < bWidth && BaseX + SignX >= 0)
                                Temp[(BaseY + SignY) * bWidth + BaseX + SignX] += mX * mY * CurValue;
                        }
                        
                    }

                double * Exch;
                Exch = Temp;
                Temp = m_EnergyMap;
                m_EnergyMap = Exch;

                for (int k = 0; k < bWidth * bHeight; k ++)
                    Add[k] += m_EnergyMap[k];

                //char strWinLen[256];
                //_itoa(w,strWinLen,10);
                //DumpEnergyMap(bWidth, bHeight, strWinLen);

            }
            delete[] Temp;

        }

        // Average the energy map within the window length
        for (int i = 0; i < bWidth * bHeight; i ++)
        {
            m_EnergyMap[i] = Add[i] / WindowLen;
            m_EnergyBuffer[i] = (BYTE)(m_EnergyMap[i] * 255);
        }
        delete[] Add;

        // dump out the final energy map within current sliding window
        if (m_verbose2)
        {
            DumpEnergyMap(bWidth, bHeight, "Average within Window (size=4)");
            //DumpEnergyBuffer(bWidth, bHeight);
        }

    }

    // based on three different templates
    double CalcEnergy(const int indexOfTPL, int bWidth, int bHeight, BOOL bNormalized=TRUE)
    {
        double Energy = 0;
        for (int i = 0; i < bWidth * bHeight; i ++)
            Energy += m_EnergyMap[i] * m_Weight[indexOfTPL][i];

        if (bNormalized){
            //normalize the energy between [0,1]
            double min = m_TPLMin[indexOfTPL]*bWidth*bHeight;
            double max = m_TPLMax[indexOfTPL]*bWidth*bHeight;
            Energy = Normalization(Energy,min,max);
            
        } else {
            // Note: this case is OK only if all the video files are with the same size.
            //       Or, videos with different (bWidth,bHeight)s can NOT be compared each other.
            //
            // only adjust according to the mean of the template
            Energy = Energy - m_TPLMean[indexOfTPL];
        }
        return Energy;
    }

    // [min, max] --> [0,1]
    double Normalization(const double x, const double min, const double max)
    {
        return (x-min)/(max-min);
    }

    // the main function for the motion vector computering
    // pre, cur: are both gray images
    // ------------------------------------------------------
    BOOL CalcMotionVectors(unsigned char* pre, unsigned char* cur, int Width, int Height, BYTE* mvseq, int index)
    {
        // bWidth, bHeight: size of block grid
        int bWidth  = Width / MVBLOCKSIZE;
        int bHeight = Height / MVBLOCKSIZE;

        int i, j;
        int x, y, x_modify, y_modify;
        char sx, sy;

        if (pre == NULL || cur == NULL)
            return TRUE;

        int count=0;
        for (j = 0, y = 0; j < bHeight; j ++, y += 16)
        {
            for (i = 0, x = 0; i < bWidth; i ++, x += 16)
            {
                x_modify = x;
                y_modify = y;
                DiamondSearch16(pre, cur, Width, Height, &x_modify, &y_modify);
                //FullSearch16(pre, cur, Width, Height, &x_modify, &y_modify);

                sx = x_modify - x;
                sy = y_modify - y;

                count += abs(sx)+abs(sy);

                //m_MVseq = new BYTE[bWidth * bHeight * 2 * (window->WindowSize()-1)];
                mvseq[bWidth*bHeight*2*index + 2*(j*bWidth+i) + 0] = sx;
                mvseq[bWidth*bHeight*2*index + 2*(j*bWidth+i) + 1] = sy;

            }
        }

        return TRUE;

    }

    // Use diamond search method to find the best matching block postion for this block
    // return the value of : &x_modify, &y_modify 
    // ------------------------------------------------
    int DiamondSearch16(unsigned char *refframe, unsigned char *desframe, 
                        int framewidth, int frameheight, int *cx, int *cy)
    {
        static char FirstDiamondX[] = { 0,  0, -2,  2, -1,  1, -1,  1};
        static char FirstDiamondY[] = {-2,  2,  0,  0, -1, -1,  1,  1};
        static char FinalDiamondX[] = {-1,  0,  1,  0}; 
        static char FinalDiamondY[] = { 0, -1,  0,  1}; 

        static char CornerDiamondX[] = {0,  0,  1,   1,  2}; 
        static char CornerDiamondY[] = {2, -2,  1,  -1,  0}; 

        static char EdgeDiamondX[] = { 0,  1,  2}; 
        static char EdgeDiamondY[] = { 2,  1,  0}; 

        int minX = 0, minY = 0, minerror;
        int curX, curY;
        int error, signX, signY;
        int sX, sY;

        int i;
        unsigned char * block;

        // the postion of current block for the current frame!
        // this positon will not change during the search process
        block = desframe + * cy * framewidth + * cx; 

        curX = * cx; curY = * cy;
        if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
        {
            * cx = * cy = 0;
            return -1;
        }

        // Initialize the minerror value
        minerror = GetMotionError16(refframe, block, framewidth, curX, curY);

        sX = * cx;
        sY = * cy;
        for (i = 0; i < 8; i ++)
        {
            curX = sX + FirstDiamondX[i]; 
            curY = sY + FirstDiamondY[i];
            if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
                continue;

            error = GetMotionError16(refframe, block, framewidth, curX, curY);
            if (error < minerror)
            {
                minX = FirstDiamondX[i];
                minY = FirstDiamondY[i];
                minerror = error;
            }
        }

        while (1)
        {
            if (!minX && !minY) break;

            sX += minX;
            sY += minY;

            if (abs(sX - * cx) > MAX_DISP) break;
            if (abs(sY - * cy) > MAX_DISP) break;
            
            if (!minY)
            {
                if (minX > 0) signX = 1;
                else signX = -1;

                minX = minY = 0;
                for (i = 0; i < 5; i ++)
                {
                    curX = sX + CornerDiamondX[i] * signX; 
                    curY = sY + CornerDiamondY[i];
                    if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
                        continue;

                    error = GetMotionError16(refframe, block, framewidth, curX, curY);
                    if (error < minerror)
                    {
                        minX = curX - sX;
                        minY = curY - sY;
                        minerror = error;
                    }
                }
                continue;
            }
            
            if (!minX)
            {
                if (minY > 0) signY = 1;
                else signY = -1;

                minX = minY = 0;
                for (i = 0; i < 5; i ++)
                {
                    curX = sX + CornerDiamondY[i]; 
                    curY = sY + CornerDiamondX[i] * signY;
                    if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
                        continue;

                    error = GetMotionError16(refframe, block, framewidth, curX, curY);
                    if (error < minerror)
                    {
                        minX = curX - sX;
                        minY = curY - sY;
                        minerror = error;
                    }
                }
                continue;
            }

            if (minX > 0) signX = 1;
            else signX = -1;
            if (minY > 0) signY = 1;
            else signY = -1;
            
            minX = minY = 0;
            for (i = 0; i < 3; i ++)
            {
                curX = sX + EdgeDiamondX[i] * signX; 
                curY = sY + EdgeDiamondY[i] * signY;
                if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
                    continue;

                error = GetMotionError16(refframe, block, framewidth, curX, curY);
                if (error < minerror)
                {
                    minX = curX - sX;
                    minY = curY - sY;
                    minerror = error;
                }
            }
        }

        for (i = 0; i < 4; i ++)
        {
            curX = sX + FinalDiamondX[i]; 
            curY = sY + FinalDiamondY[i];
            if (curX < 0 || curX > framewidth - 16 || curY < 0 || curY > frameheight - 16) 
                continue;

            error = GetMotionError16(refframe, block, framewidth, curX, curY);
            if (error < minerror)
            {
                minX = FinalDiamondX[i];
                minY = FinalDiamondY[i];
                minerror = error;
            }
        }

        * cx = sX + minX;
        * cy = sY + minY;
        return minerror;

    }

    int FullSearch16(unsigned char *refframe, unsigned char *desframe, int framewidth, int frameheight, int *cx, int *cy)
    {
        int minerror;
        int curX, curY;
        int error, oX, oY;
        int sx, sy;

        int i, j;
        unsigned char * block;
        block = desframe + * cy * framewidth + * cx;

        // first
        oX = curX = * cx; 
        oY = curY = * cy;
        if (curX < 0 || curX > framewidth - 16|| curY < 0 || curY > frameheight - 16) 
        {
            * cx = * cy = 0;
            return -1;
        }
        minerror = GetMotionError16(refframe, block, framewidth, curX, curY);

        for (i = 1; i <= MAX_DISP; i ++)
        {
            sx = sy = -i;
            for (j = 0; j < 8 * i; j ++)
            {
                curX = oX + sx;
                curY = oY + sy;

                if (!(curX < 0 || curX > framewidth - 16|| curY < 0 || curY > frameheight - 16))
                {
                    error = GetMotionError16(refframe, block, framewidth, curX, curY);
                    if (error < minerror)
                    {
                        minerror = error;
                        * cx = curX;
                        * cy = curY;
                    }
                }

                if (j < 2 * i) sx++;
                else if (j < 4 * i) sy++;
                else if (j < 6 * i) sx--;
                else sy--;
            }
        }

        return minerror;

    }


    // ------------------------------------------------------------------------------------
    //    refframe:   the former frame buffer  (with fx,fy )
    //    block   :   the current block position in the current frame buffer (not change!)
    //  framewidth:   the whole frame width
    //    fx ,fy  :   the ref block position       (fx<=Width,fy<=Height)
    // ------------------------------------------------------------------------------------
    int GetMotionError16(unsigned char * refframe, unsigned char * block, 
                         int framewidth, int fx, int fy)
    {
        int i,j;
        int SAD = 0;
        BYTE* refBlock = refframe + fy*framewidth + fx;

        // i: fx direction (width)
        // j: fy direction (height)
        for (j=0;j<MVBLOCKSIZE;j++)
            for (i=0;i<MVBLOCKSIZE;i++)
            {
                int diff = block[j*framewidth+i] - refBlock[j*framewidth+i];
                diff = abs(diff);
                SAD+=diff;
            }
        return SAD;
    }


protected:

    int
    GetBinCount() const
    {
        return mBinCount;
    }

    //virtual VectorReal64
    //ComputeHistogram(Array::Array2dVec3UInt8* image) = 0;

    // compute rgb histogram
    /*VectorReal64
    ComputeHistogram(Array::Array2dVec3UInt8* image)
    {
        UInt8* data = image->CPB(0, 0);
        const int nPix = image->CW() * image->CH();

        VectorReal64 histogram(GetBinCount());
        Real64* histR = histogram.GetData();
        Real64* histG = histR + BINS_PER_CHANNEL;
        Real64* histB = histG + BINS_PER_CHANNEL;

        for (int i = 0; i < GetBinCount(); i++)
            histR[i] = 0;

        for (int j = 0; j < nPix; j++) 
        {
            const int r = ComputeBin((int) data[0], 0, 255, BINS_PER_CHANNEL);
            histR[r] += 1;

            const int g = ComputeBin((int) data[1], 0, 255, BINS_PER_CHANNEL);
            histG[g] += 1;

            const int b = ComputeBin((int) data[2], 0, 255, BINS_PER_CHANNEL);
            histB[b] += 1;

            data += 3;
        }

        const Real64 factor = nPix;
        Vector::DivAssign(histogram, factor);

        return histogram;
    }*/

    template <class T>
    int
    ComputeBin(const T val, const T low, const T high, const int nrBins) const
    {
        const int bin = (nrBins * (val - low)) / (high - low);
        if (bin >= 0 && bin < nrBins)
            return bin;

        if (val == high)
            return nrBins - 1;

        ILOG_WARN("Unexpected value: " << val << 
            " is not in range [" << low << ", " << high << "]");
        return (bin < 0) ? 0 : nrBins-1;
    }

private:

    //static const int BINS_PER_CHANNEL = 32;

    Reporter*              mReporter;
    String                 mFeatureName;
    int                    mBinCount;
    Feature::FeatureTable* mHistogramTable;

    // calc the energy map from motion vectors
    double*                m_Weight[TPL_NUM];
    double                 m_TPLMean[TPL_NUM];
    double                 m_TPLMin[TPL_NUM];
    double                 m_TPLMax[TPL_NUM];
    double*                m_EnergyMap;

    BOOL                   m_bNormalize;

    BOOL                   IsInShot;
    BYTE*                  m_MVseq;   // store all the motion vectors field within the window
    BYTE*                  m_Buffer;
    BYTE*                  m_EnergyBuffer;

    BOOL                   m_verbose;
    BOOL                   m_verbose2; // for dedug only; more detailed info

    ILOG_VAR_DECL;
};

ILOG_VAR_INIT(KfrMotionExtractor, Impala.Core.VideoSet);

} // namespace VideoSet
} // namespace Core
} // namespace Impala

#endif

---