/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2015, Itseez Inc, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Itseez Inc or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "dpm_feature.hpp"

using namespace std;

namespace cv
{
namespace dpm
{

Feature::Feature()
{
}

Feature::Feature (PyramidParameter p):params(p)
{
}

void Feature::computeFeaturePyramid(const Mat &imageM, vector< Mat > &pyramid)
{
#ifdef HAVE_TBB
    ParalComputePyramid paralTask(imageM, pyramid, params);
    paralTask.initialize();
    // perform parallel computing
    parallel_for_(Range(0, params.interval), paralTask);
#else
    CV_Assert(params.interval > 0);
    // scale factor between two levels
    params.sfactor = pow(2.0, 1.0/params.interval);
    const Size_<double> imSize = imageM.size();
    params.maxScale = 1 + (int)floor(log(min(imSize.width, imSize.height)/
                (float)(params.binSize*5.0))/log(params.sfactor));

    if (params.maxScale < params.interval)
    {
        CV_Error(CV_StsBadArg, "The image is too small to create a pyramid");
        return;
    }

    pyramid.resize(params.maxScale + params.interval);
    params.scales.resize(params.maxScale + params.interval);

#ifdef _OPENMP
#pragma omp parallel for
#endif
    for (int i = 0; i < params.interval; i++)
    {
        const double scale = (double)(1.0f/pow(params.sfactor, i));
        Mat imScaled;
        resize(imageM, imScaled, imSize * scale);
        // First octave at twice the image resolution
        computeHOG32D(imScaled, pyramid[i], params.binSize/2,
                params.padx + 1, params.pady + 1);

        params.scales[i] = 2*scale;

        // Second octave at the original resolution
        if (i + params.interval <= params.maxScale)
            computeHOG32D(imScaled, pyramid[i+params.interval],
                    params.binSize, params.padx + 1, params.pady + 1);

        params.scales[i+params.interval] = scale;

        // Remaining octaves
        for ( int j = i + params.interval; j < params.maxScale; j += params.interval)
        {
            Mat imScaled2;
            Size_<double> imScaledSize = imScaled.size();
            resize(imScaled, imScaled2, imScaledSize*0.5);
            imScaled = imScaled2;
            computeHOG32D(imScaled2, pyramid[j+params.interval],
                    params.binSize, params.padx + 1, params.pady + 1);
            params.scales[j+params.interval] = params.scales[j]*0.5;
        }
    }
#endif
}

#ifdef HAVE_TBB
ParalComputePyramid::ParalComputePyramid(const Mat &inputImage, \
        vector< Mat > &outputPyramid,\
        PyramidParameter &p):
    imageM(inputImage), pyramid(outputPyramid), params(p)
{
}

void ParalComputePyramid::initialize()
{
    CV_Assert(params.interval > 0);

    // scale factor between two levels
    params.sfactor = pow(2.0, 1.0/params.interval);
    imSize = imageM.size();
    params.maxScale = 1 + (int)floor(log(min(imSize.width, imSize.height)/(float)(params.binSize*5.0))/log(params.sfactor));

    if (params.maxScale < params.interval)
    {
        CV_Error(CV_StsBadArg, "The image is too small to create a pyramid");
        return;
    }

    pyramid.resize(params.maxScale + params.interval);
    params.scales.resize(params.maxScale + params.interval);
}

void ParalComputePyramid::operator() (const Range &range) const
{
    for (int i = range.start; i != range.end; i++)
    {
        const double scale = (double)(1.0f/pow(params.sfactor, i));
        Mat imScaled;
        resize(imageM, imScaled, imSize * scale);

        params.scales[i] = 2*scale;

        // First octave at twice the image resolution
        Feature::computeHOG32D(imScaled, pyramid[i],
                params.binSize/2, params.padx + 1, params.pady + 1);

        // Second octave at the original resolution
        if (i + params.interval <= params.maxScale)
            Feature::computeHOG32D(imScaled, pyramid[i+params.interval],
                    params.binSize, params.padx + 1, params.pady + 1);

        params.scales[i+params.interval] = scale;

        // Remaining octaves
        for ( int j = i + params.interval; j < params.maxScale; j += params.interval)
        {
            Mat imScaled2;
            Size_<double> imScaledSize = imScaled.size();
            resize(imScaled, imScaled2, imScaledSize*0.5);
            imScaled = imScaled2;
            Feature::computeHOG32D(imScaled2, pyramid[j+params.interval],
                    params.binSize, params.padx + 1, params.pady + 1);
            params.scales[j+params.interval] = params.scales[j]*0.5;
        }
    }
}
#endif

void Feature::computeHOG32D(const Mat &imageM, Mat &featM, const int sbin, const int pad_x, const int pad_y)
{
    CV_Assert(pad_x >= 0);
    CV_Assert(pad_y >= 0);
    CV_Assert(imageM.channels() == 3);
    CV_Assert(imageM.depth() == CV_64F);

    // epsilon to avoid division by zero
    const double eps = 0.0001;
    // number of orientations
    const int numOrient = 18;
    // unit vectors to compute gradient orientation
    const double uu[9] = {1.000, 0.9397, 0.7660, 0.5000, 0.1736, -0.1736, -0.5000, -0.7660, -0.9397};
    const double vv[9] = {0.000, 0.3420, 0.6428, 0.8660, 0.9848,  0.9848,  0.8660,  0.6428,  0.3420};

    // image size
    const Size imageSize = imageM.size();
    // block size
    int bW = cvRound((double)imageSize.width/(double)sbin);
    int bH = cvRound((double)imageSize.height/(double)sbin);
    const Size blockSize(bW, bH);
    // size of HOG features
    int oW = max(blockSize.width-2, 0) + 2*pad_x;
    int oH = max(blockSize.height-2, 0) + 2*pad_y;
    Size outSize = Size(oW, oH);
    // size of visible
    const Size visible = blockSize*sbin;

    // initialize historgram, norm, output feature matrices
    Mat histM = Mat::zeros(Size(blockSize.width*numOrient, blockSize.height), CV_64F);
    Mat normM = Mat::zeros(Size(blockSize.width, blockSize.height), CV_64F);
    featM = Mat::zeros(Size(outSize.width*dimHOG, outSize.height), CV_64F);

    // get the stride of each matrix
    const size_t imStride = imageM.step1();
    const size_t histStride = histM.step1();
    const size_t normStride = normM.step1();
    const size_t featStride = featM.step1();

    // calculate the zero offset
    const double* im = imageM.ptr<double>(0);
    double* const hist = histM.ptr<double>(0);
    double* const norm = normM.ptr<double>(0);
    double* const feat = featM.ptr<double>(0);

    for (int y = 1; y < visible.height - 1; y++)
    {
        for (int x = 1; x < visible.width - 1; x++)
        {
            // OpenCV uses an interleaved format: BGR-BGR-BGR
            const double* s = im + 3*min(x, imageM.cols-2) + min(y, imageM.rows-2)*imStride;

            // blue image channel
            double dyb = *(s+imStride) - *(s-imStride);
            double dxb = *(s+3) - *(s-3);
            double vb = dxb*dxb + dyb*dyb;

            // green image channel
            s += 1;
            double dyg = *(s+imStride) - *(s-imStride);
            double dxg = *(s+3) - *(s-3);
            double vg = dxg*dxg + dyg*dyg;

            // red image channel
            s += 1;
            double dy = *(s+imStride) - *(s-imStride);
            double dx = *(s+3) - *(s-3);
            double v = dx*dx + dy*dy;

            // pick the channel with the strongest gradient
            if (vg > v) { v = vg; dx = dxg; dy = dyg; }
            if (vb > v) { v = vb; dx = dxb; dy = dyb; }

            // snap to one of the 18 orientations
            double best_dot = 0;
            int best_o = 0;
            for (int o = 0; o < (int)numOrient/2; o++)
            {
                double dot =  uu[o]*dx + vv[o]*dy;
                if (dot > best_dot)
                {
                    best_dot = dot;
                    best_o = o;
                }
                else if (-dot > best_dot)
                {
                    best_dot = -dot;
                    best_o = o + (int)(numOrient/2);
                }
            }

            // add to 4 historgrams around pixel using bilinear interpolation
            double yp =  ((double)y+0.5)/(double)sbin - 0.5;
            double xp =  ((double)x+0.5)/(double)sbin - 0.5;
            int iyp = (int)floor(yp);
            int ixp = (int)floor(xp);
            double vy0 = yp - iyp;
            double vx0 = xp - ixp;
            double vy1 = 1.0 - vy0;
            double vx1 = 1.0 - vx0;
            v = sqrt(v);

            // fill the value into the 4 neighborhood cells
            if (iyp >= 0 && ixp >= 0)
                *(hist + iyp*histStride + ixp*numOrient + best_o) += vy1*vx1*v;

            if (iyp >= 0 && ixp+1 < blockSize.width)
                *(hist + iyp*histStride + (ixp+1)*numOrient + best_o) += vx0*vy1*v;

            if (iyp+1 < blockSize.height && ixp >= 0)
                *(hist + (iyp+1)*histStride + ixp*numOrient + best_o) += vy0*vx1*v;

            if (iyp+1 < blockSize.height && ixp+1 < blockSize.width)
                *(hist + (iyp+1)*histStride + (ixp+1)*numOrient + best_o) += vy0*vx0*v;

        } // for y
    } // for x

    // compute the energy in each block by summing over orientation
    for (int y = 0; y < blockSize.height; y++)
    {
        const double* src = hist + y*histStride;
        double* dst = norm + y*normStride;
        double const* const dst_end = dst + blockSize.width;
        // for each cell
        while (dst < dst_end)
        {
            *dst = 0;
            for (int o = 0; o < (int)(numOrient/2); o++)
            {
                *dst += (*src + *(src + numOrient/2))*
                    (*src + *(src + numOrient/2));
                src++;
            }
            dst++;
            src += numOrient/2;
        }
    }

    // compute the features
    for (int y = pad_y; y < outSize.height - pad_y; y++)
    {
        for (int x = pad_x; x < outSize.width - pad_x; x++)
        {
            double* dst = feat + y*featStride + x*dimHOG;
            double* p, n1, n2, n3, n4;
            const double* src;

            p = norm + (y - pad_y + 1)*normStride + (x - pad_x + 1);
            n1 = 1.0f / sqrt(*p + *(p + 1) + *(p + normStride) + *(p + normStride + 1) + eps);
            p = norm + (y - pad_y)*normStride + (x - pad_x + 1);
            n2 = 1.0f / sqrt(*p + *(p + 1) + *(p + normStride) + *(p + normStride + 1) + eps);
            p = norm + (y- pad_y + 1)*normStride + x - pad_x;
            n3 = 1.0f / sqrt(*p + *(p + 1) + *(p + normStride) + *(p + normStride + 1) + eps);
            p = norm + (y - pad_y)*normStride + x - pad_x;
            n4 = 1.0f / sqrt(*p + *(p + 1) + *(p + normStride) + *(p + normStride + 1) + eps);

            double t1 = 0.0, t2 = 0.0, t3 = 0.0, t4 = 0.0;

            // contrast-sesitive features
            src = hist + (y - pad_y + 1)*histStride + (x - pad_x + 1)*numOrient;
            for (int o = 0; o < numOrient; o++)
            {
                double val = *src;
                double h1 = min(val*n1, 0.2);
                double h2 = min(val*n2, 0.2);
                double h3 = min(val*n3, 0.2);
                double h4 = min(val*n4, 0.2);
                *(dst++) = 0.5 * (h1 + h2 + h3 + h4);
                src++;
                t1 += h1;
                t2 += h2;
                t3 += h3;
                t4 += h4;
            }

            // contrast-insensitive features
            src =  hist + (y - pad_y + 1)*histStride + (x - pad_x + 1)*numOrient;
            for (int o = 0; o < numOrient/2; o++)
            {
                double sum = *src + *(src + numOrient/2);
                double h1 = min(sum * n1, 0.2);
                double h2 = min(sum * n2, 0.2);
                double h3 = min(sum * n3, 0.2);
                double h4 = min(sum * n4, 0.2);
                *(dst++) = 0.5 * (h1 + h2 + h3 + h4);
                src++;
            }

            // texture features
            *(dst++) = 0.2357 * t1;
            *(dst++) = 0.2357 * t2;
            *(dst++) = 0.2357 * t3;
            *(dst++) = 0.2357 * t4;

            // truncation feature
            *dst = 0;
        }// for x
    }// for y

    // Truncation features
    for (int m = 0; m < featM.rows; m++)
    {
        for (int n = 0; n < featM.cols; n += dimHOG)
        {
            if (m > pad_y - 1 && m < featM.rows - pad_y && n > pad_x*dimHOG - 1 && n < featM.cols - pad_x*dimHOG)
                continue;

            featM.at<double>(m, n + dimHOG - 1) = 1;
        } // for x
    }// for y
}

void Feature::projectFeaturePyramid(const Mat &pcaCoeff, const std::vector< Mat > &pyramid, std::vector< Mat > &projPyramid)
{
    CV_Assert(dimHOG == pcaCoeff.rows);
    dimPCA = pcaCoeff.cols;

    projPyramid.resize(pyramid.size());

    // loop for each level of the pyramid
    for (unsigned int i = 0; i < pyramid.size(); i++)
    {
        Mat orgM = pyramid[i];
        // note that the features are stored in 32-32-32
        int width = orgM.cols/dimHOG;
        int height = orgM.rows;
        // initialize the project feature matrix
        Mat projM = Mat::zeros(height, width*dimPCA, CV_64F);
        //get the pointer of the matrix
        double* const featOrg = orgM.ptr<double>(0);
        double* const featProj = projM.ptr<double>(0);

        // get the stride of each matrix
        const size_t orgStride = orgM.step1();
        const size_t projStride = projM.step1();

        for (int y = 0; y < height; y++)
        {
            for (int x = 0; x < width; x++)
            {
                double* proj = featProj + y*projStride + x*dimPCA;
                // for each pca dimension
                for (int c = 0; c < dimPCA; c++)
                {
                    double* org = featOrg + y*orgStride + x*dimHOG;
                    // dot product 32d HOG feature with the coefficient vector
                    for (int r = 0; r < dimHOG; r++)
                    {
                        *proj += *org * pcaCoeff.at<double>(r, c);
                        org++;
                    }
                    proj++;
                }

            } // for x
        } // for y
        projPyramid[i] = projM;
    } // for each level of the pyramid
}

void Feature::computeLocationFeatures(const int numLevels, Mat &locFeature)
{
    locFeature = Mat::zeros(Size(numLevels, 3), CV_64F);

    int b = 0;
    int e = min(numLevels, params.interval);

    for (int x = b; x < e; x++)
        locFeature.at<double>(0, x) = 1;

    b = e;
    e = min(numLevels, 2*e);

    for (int x = b; x < e; x++)
        locFeature.at<double>(1, x) = 1;

    b = e;
    e = min(numLevels, 3*e);

    for (int x = b; x < e; x++)
        locFeature.at<double>(2, x) = 1;
}
} // namespace dpm
} // namespace cv
