/*! 
*   \file  CFeatureExtraction.cpp
*   \brief Defines class to extract features.
*    
*   \details  
*   CFeatureExtraction extracts the requested features. 
*
*   \date      October 24, 2011
*   \copyright eBay Research Labs.
*/

#include <iomanip> // For setprecision.

// Include openCV-specific headers.
#include "opencv2/core/core.hpp"
#include <highgui.h>

#include "CFeatureExtraction.h"

#include "CFST.h"
#include "CMD5AndUrl.h"
#include "UtilitiesWeb.h"

#include "Base64.h"
#include <sstream>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/iterators/binary_from_base64.hpp>
#include <boost/archive/iterators/base64_from_binary.hpp>
#include <boost/archive/iterators/transform_width.hpp>
#include "CVMatSerialization.h"

// Namespaces.
using namespace std;
using namespace cv;


// *************************************************************************************************************
// ********************************************   Public Methods   *********************************************
// *************************************************************************************************************


/*!
*   \brief  Destructor.
*/
CFeatureExtraction::~CFeatureExtraction()
{
    delete m_sig;
}

/*!
*   \brief  Extract necessary features from the input image.
*   \param [in]  imgIn     Input image.
*   \param [in]  ftrMapper    reference to CFeatureMapper object in order to extract Matsuda code words.
*   \param [in]  flagCrop     Flag to indicate whether to crop image using ROI rectangle. Default = true.
*   \return                Cropped/rescaled version of imgIn.
*/
CFeatureExtraction::CFeatureExtraction(const Mat &imgIn, CFeatureMapper &ftrMapper, bool flagCrop)
{
    m_sig = new CSignature();
    m_sig->m_version = 3.0; // Current version.

    m_md5 = Mantis::createMD5HashForImage(imgIn);

    //*************************************************************************************
    // Step 1: Crop and resize image.
    //*************************************************************************************
    
    Mat  rgb;
    Rect roi;
    
    if(flagCrop)
    {
        DefaultRectROI(imgIn.size(), roi);
    }

    ResizeImageWithROI(imgIn, rgb, roi, flagCrop); // Resize image if the ROI is big. No data copy when not scaled.

    //*************************************************************************************
    // Step 2: Extract features based on color
    //*************************************************************************************
    Mat hsv;
    cvtColor(rgb, hsv, CV_BGR2HSV);

    m_sig->m_histColor = ClrHSV_BWValHist(hsv, m_sig->m_fracColor);    

    //*************************************************************************************
    // Step 3: Extract features based on pattern/texture.
    //*************************************************************************************    
    Mat gray;
    cvtColor(rgb, gray, CV_RGB2GRAY);
    //OrientHist(gray, m_sig->m_histOrient, m_sig->m_fracEdgePixels, m_sig->m_maxEdgeMag);  
    OrientHist(rgb, m_sig->m_histOrient, m_sig->m_fracEdgePixels, m_sig->m_maxEdgeMag);  

    int dimBlocks[] = {2, 2}; // rows x columns
    int numScales   = 4;
    int numAngles   = 8;
    CFST fst(dimBlocks, numScales, numAngles);
    Mat ftrs;
    fst.CalculateFeatures(gray, m_sig->m_FST);

    //*************************************************************************************
    // Step 4: Extract Matsuda code word for color distribution.
    //************************************************************************************* 
    vector <Rect> boundBoxes;
    vector <int>  codeWords;
    boundBoxes.push_back(Rect(0, 0, rgb.cols, rgb.rows)); 
    ftrMapper.ComputeCodeWords(rgb, boundBoxes, codeWords);
    m_sig->m_matsudaHue = codeWords[0];
}

// This limited constructor only creates color signature (m_sig->m_fracColor, m_histColor)
CFeatureExtraction::CFeatureExtraction(const Mat& imgIn, bool flagCrop)
{
    m_sig = new CSignature();
    m_sig->m_version = 3.0; // Current version.

    m_md5 = Mantis::createMD5HashForImage(imgIn);

    //*************************************************************************************
    // Step 1: Crop and resize image.
    //*************************************************************************************
    
    Mat  rgb;
    Rect roi;
    
    if(flagCrop)
    {
        DefaultRectROI(imgIn.size(), roi);
    }

    ResizeImageWithROI(imgIn, rgb, roi, flagCrop); // Resize image if the ROI is big. No data copy when not scaled.

    //*************************************************************************************
    // Step 2: Extract features based on color
    //*************************************************************************************
    Mat hsv;
    cvtColor(rgb, hsv, CV_BGR2HSV);

    m_sig->m_histColor = ClrHSV_BWValHist(hsv, m_sig->m_fracColor);    
}

CFeatureExtraction::CFeatureExtraction() : m_sig(0)
{

}


/*!
*   \brief  Writes features to strean in JSON  format.
*   \param [inout]   stream.
*/
std::ostream&  CFeatureExtraction::SerializeSignatureJSON(std::ostream& stream)
{
  if (m_sig)
    {
      stream << "{ \"MD5\": \"" 
	     <<  mantis_base64_encode((unsigned char*)m_md5.data(), m_md5.length()) 
	     << "\", \"signature\": \""
	     << SerializeSignature() << "\" }";
    }
  return stream;
}

/*!
*   \brief  Writes features to a stream in JSON format with text only.
*   \param [inout]  stream.
*/
std::ostream& CFeatureExtraction::SerializeSignatureJSONText(std::ostream& stream, bool contentOnly /* = false */)
{
  if (m_sig)
  {
    if (!contentOnly)
    {
      stream << "{ ";
    }
    // The consumer of this data requires that 0.0f be represented as "0.0" in the output. I could achieve this
    // by using stream.setf(ios_base::showpoint), but it would bloat up the output unnecessarily. So instaed I
    // just treat zeros specially.
    stream << "\"fracColor\": ";
    if (m_sig->m_fracColor == 0)
      stream << "0.0";
    else stream << m_sig->m_fracColor;
    stream << ", " << " \"histColor\": [ ";
    const float* M = m_sig->m_histColor.ptr<float>(0);
    for(int i = 0; i < m_sig->m_histColor.rows; i++)
    {
      stream << (i > 0 ? ", " : "");
      if (M[i] == 0)
        stream << "0.0";
      else stream << M[i];
    }
    stream << " ]";
    if (!contentOnly)
    {
      stream << " }";
    }
  }
  return stream;
}

std::string  CFeatureExtraction::SerializeSignature()
{
  std::stringstream s( std::ios_base::out | std::ios_base::binary);
  s << std::noskipws;
  boost::archive::binary_oarchive oa(s);

  oa << m_sig->m_version;
  oa << m_sig->m_matsudaHue;
  oa << m_sig->m_fracColor;
  oa << m_sig->m_maxEdgeMag;
  oa << m_sig->m_histColor;
  oa << m_sig->m_histOrient;
  oa << m_sig->m_fracEdgePixels;
  oa << m_sig->m_FST;
  s.flush();

  std::string sbuf(s.str());
  std::string ret = mantis_base64_encode((unsigned char*)sbuf.data(), sbuf.length());
//  cerr << "serialized string is " << ret << endl;
  return ret;
}

/*!
*   \brief  Writes features to strean in JSON  format.
*   \param [inout]   stream.
*/
std::ostream&  CFeatureExtraction::SerializeSignatureXML(std::ostream& stream)
{
  if (m_sig)
    {
      stream << "<MD5>" 
	     <<  mantis_base64_encode((unsigned char*)m_md5.data(), m_md5.length()) 
	     << "</MD5>/n<signature>"
	     <<  SerializeSignature() << "</signature>";
    }
  return stream;
}



void CFeatureExtraction::DeserializeSignature(const std::string& b64E)
{
  //  cerr << "string is " << b64E << endl; 
  std::stringstream is(std::ios_base::in | std::ios_base::out | std::ios_base::binary ) ;
  is << std::noskipws;
  is << mantis_base64_decode(b64E);

  boost::archive::binary_iarchive ia(is);
  delete m_sig;
  m_sig = new CSignature();

  ia >> m_sig->m_version;
  ia >> m_sig->m_matsudaHue;
  ia >> m_sig->m_fracColor;
  ia >> m_sig->m_maxEdgeMag;
  ia >> m_sig->m_histColor;
  ia >> m_sig->m_histOrient;
  ia >> m_sig->m_fracEdgePixels;
  ia >> m_sig->m_FST;
  
}


// *************************************************************************************************************
// ********************************************   Private Methods   ********************************************
// *************************************************************************************************************


/*!
*   \brief  If crop ROI is too big in terms of #pixels, the image will be resized by this function.
*
*   \param [in]  sz    Size of input image.
*   \param [out] roi   Default rectangular ROI, based on image size.
*/
void CFeatureExtraction::DefaultRectROI(Size sz, Rect &roi)
{
    // TODO: Needs more experimentation.

    // Women's clothing:
    roi.x      = cvRound(0.38*sz.width);
    roi.y      = cvRound(0.30*sz.height);
    roi.width  = cvRound(0.24*sz.width);
    roi.height = cvRound(0.40*sz.height);                  
}


/*!
*   \brief  If crop ROI is too big in terms of #pixels, the image will be resized by this function.
*
*   \param [in]  imgIn     Input image.
*   \param [out] imgOut    Resized image.
*   \param [in]  roi       Rectangle ROI.
*   \param [in]  flagCrop     Flag to indicate whether to crop image using ROI rectangle. Default = true.
*   \param [in]  maxDim    Max dimension (width or height) of rectangular ROI. Default = 128.
*/
void CFeatureExtraction::ResizeImageWithROI(const Mat &imgIn, Mat &imgOut, Rect &roi, bool flagCrop, int maxDim)
{
    float h, w;

    if(flagCrop)
    {
        w = (float) roi.width;
        h = (float) roi.height;
    }
    else
    {
        // Don't crop.
        w = (float) imgIn.cols;
        h = (float) imgIn.rows;
    }

    float fx, fy;
    fx = fy = 1.0;

    if (w>maxDim)
    {
        fx = (float) maxDim / w;
    }
                 
    if (h>maxDim)
    {
        fy = (float) maxDim / h;
    } 
    if (fx<fy) 
    {
        fy = fx;
    }
    else
    {
        fx = fy;
    }

    if (fx<1.0)
    {        
        if(flagCrop)
        {
             resize(imgIn(roi), imgOut, Size(), fx, fy, INTER_LINEAR);   // TODO: Anti-alias based on scale factor.
        }
        else
        {
            resize(imgIn, imgOut, Size(), fx, fy, INTER_LINEAR);  // TODO: Anti-alias based on scale factor.
        }
    }
    else
    {
        // Don't scale. 
        if(flagCrop)
        {
            imgOut = imgIn(roi);  // No copy of data.  
        }
        else
        {
            imgOut = imgIn;  // No copy of data.  
        }
    }
}


/*!
*   \brief  Calculate 1D histogram of values in a Mat variable.
*
*   \param [in] data2D     Input matrix.
*   \param [in] range      Range of values in data2D to be binned.
*   \param [in] numBins    #bins in histogram.
*   \return                1-D histogram.
*/
Mat CFeatureExtraction::Hist1D(const Mat& data2D, const float *range, const int numBins)
{
    int          histSize[]  = {numBins};
    const float* rangesAll[] = {range};
    int          planes[]    = {0};
    Mat          hist;

    calcHist(&data2D, 1, planes, Mat(), hist, 1, histSize, rangesAll);

    return hist;
} // End of function: Hist1D


/*!
*   \brief  Calculates HSV histogram for color and V histogram for BW.
*
*   This partitions pixels into color and b/w. Then finds 1D histograms for 
*   hue, saturation and value. 1D histogram of value is extracted for b/w pixels.
*   feature subsets are then weighted. the histograms are concated into one 
*   single column.
*
*   \param [in]  hsv     Input HSV image.
*   \param [out] fracClr Fraction of color pixels. Value in [0,1].
*   \return              Clr-HSV and BW-Val histogram. 
*/
Mat CFeatureExtraction::ClrHSV_BWValHist(const Mat &hsv, float &fracClr)
{
    CV_Assert( !hsv.data || hsv.type() == CV_8UC3 );
    int numBinsHue    = 24;
    int numBinsSatVal = 8;
    int numBinsVal    = 8;
    float wts[]       = {0.4F, 0.2F, 0.1F, 0.3F}; // Feature weights: clrHue, clrSat, clrVal, bwVal. This must sum to 1.0.
    int numBins       = numBinsHue+2*numBinsSatVal+numBinsVal;
    Mat histHueVal    = Mat::zeros(numBins, 1, CV_32F);
    
    int index;
    double scaleHue    = numBinsHue/180.0;
    double scaleSatVal = numBinsSatVal/256.0;
    double scaleVal    = numBinsVal/256.0;

    int   hue, sat, val;

    Rect clrHueROI;
    clrHueROI.x      = 0;
    clrHueROI.y      = 0;
    clrHueROI.width  = 1;
    clrHueROI.height = numBinsHue;  

    Rect clrSatROI;
    clrSatROI.x      = 0;
    clrSatROI.y      = numBinsHue;
    clrSatROI.width  = 1;
    clrSatROI.height = numBinsSatVal;  

    Rect clrValROI;
    clrValROI.x      = 0;
    clrValROI.y      = numBinsHue+numBinsSatVal;
    clrValROI.width  = 1;
    clrValROI.height = numBinsSatVal;  
                
    Rect bwValROI;
    bwValROI.x      = 0;
    bwValROI.y      = numBinsHue+2*numBinsSatVal;
    bwValROI.width  = 1;
    bwValROI.height = numBinsVal;
    
    Mat histClrHue =  histHueVal(clrHueROI);
    Mat histClrSat =  histHueVal(clrSatROI);
    Mat histClrVal =  histHueVal(clrValROI);
    Mat histBWVal  =  histHueVal(bwValROI);

    MatConstIterator_< Vec<uchar, 3> > hsv_iter;
    int numClr = 0;
    for (hsv_iter = hsv.begin<Vec<uchar, 3> >(); hsv_iter != hsv.end<Vec<uchar, 3> >(); hsv_iter++)
    {
        hue = (*hsv_iter)[0]; 
        sat = (*hsv_iter)[1]; 
        val = (*hsv_iter)[2]; 

        
        if ((float) sat * (float) val < 15.3*255.0) // 6% of 255 is 15.3
        {
            // Hue is not reliable. Update "value" histogram.
            index = (int) (val*scaleVal); 
            histBWVal.at<float>(index)++;
        }
        else
        {
            // Hue is reliable. Update "hue" histogram.
            index = (int) (hue*scaleHue);
            histClrHue.at<float>(index)++;

            index = (int) (sat*scaleSatVal);
            histClrSat.at<float>(index)++;

            index = (int) (val*scaleSatVal);
            histClrVal.at<float>(index)++;

            numClr++;
        }    
    }  

    // Apply weights to subsets of features.
    histClrHue *=  wts[0];
    histClrSat *=  wts[1];
    histClrVal *=  wts[2];
    histBWVal  *=  wts[3];
   
    fracClr = (float) numClr / (float) hsv.total();

    return histHueVal;
}


/*!
*   \brief  Calculates orientation histogram of edges.
*
*   Then finds 1D histogram of orientation of edges detected by Canny detector. 
*   Contributions are weighted based on edge magnitude.
*
*   \param [in]  numAngBins    number of bins for orientation in [0, 180).
*   \param [in]  dx            X-derivative.
*   \param [in]  dy            Y-derivative.
*   \param [in]  imgEdge       Edge map.
*   \param [out] histOrient    Histogram of orientation, weighed by gradient strength.
*                              This needs to be allocated and initialized to zero values.
*   \param [out] numEdgePixels Number of edge pixels.
*   \param [out] maxEdgeMag    Maximum magnitude of edges.
*/
void CFeatureExtraction::GradientHist(int numAngBins, const Mat &dx, const Mat &dy, const Mat &imgEdge, Mat &histOrient, int &numEdgePixels, float &maxEdgeMag)
{
    // Assume histOrient is reset to have zero values.
    uchar* pEdge = imgEdge.data;
    MatConstIterator_<float> dx_iter = dx.begin<float>();
    MatConstIterator_<float> dy_iter = dy.begin<float>();

    float angle, mag;            
    int   index;
    maxEdgeMag = 0;
    // Get weighted histogram of gradient, weighted by edge strength.
    for(; dx_iter != dx.end<float>(); ++pEdge, ++dx_iter, ++dy_iter)
    {
        if ((int)(*pEdge)<1)
            continue;

        angle = fastAtan2(*dy_iter, *dx_iter); // Angle will be in [0, 360].
        mag   = sqrt((*dy_iter)*(*dy_iter) + (*dx_iter)*(*dx_iter))/2;

        if (maxEdgeMag < mag) maxEdgeMag = mag;

        if (angle > 180.0)  angle = 360-angle; // Range: [0, 180];
        if (angle == 180.0) angle = 0.0; // Range: [0, 180);

        index       = cvRound(angle/180.0*(numAngBins-1)); // Range: [0, numAngBins-1]
        histOrient.at<float>(index) += mag;// Weighting by edge strength.     
    }
    numEdgePixels = countNonZero(imgEdge); // TODO: Experiment with other thresholds.
}


/*!
*   \brief  Calculates orientation histogram of edges.
*
*   Then finds 1D histogram of orientation of edges detected by Canny detector. 
*   Contributions are weighted based on edge magnitude.
*
*   \param [in]  rgb            Rgb image.
*   \param [out] fracEdgePixels Fraction of edge pixels in each vertical partition of image. Value in [0,1].
*   \param [out] histOrient     Orientation histogram. 
*   \param [out] maxEdgeMag     Maximum magnitude of edges.
*/
void CFeatureExtraction::OrientHist(const Mat &rgb, Mat &histOrient, Mat &fracEdgePixels, float &maxEdgeMag)
{
  CV_Assert( rgb.data && rgb.type() == CV_8UC3 );

  int height     = rgb.rows;
  int width      = rgb.cols;

  // Split RGB channels.
  int numChannels = rgb.channels();
  vector <Mat> channels;
  split(rgb, channels);

  // Iterate through channels to get gradient. Use max-norm across channels to get the final gradient.
  Mat dx, dy, imgEdge;
  dx      = Mat::zeros(height, width, CV_32F); 
  dy      = Mat::zeros(height, width, CV_32F); 
  imgEdge = Mat::zeros(height, width, CV_8U); 
  for (int ci=0; ci<numChannels; ci++)
    {
      Mat gray;
      GaussianBlur(channels[ci], gray, Size(9, 9), 3, 3);  // Smooth image first.

      // Get edge image.
      Mat imgEdgeCurrent;
      Canny(gray, imgEdgeCurrent, 5, 15, 3); // TODO: Choose best thresholds. 
      float numEdgePixels = (float) countNonZero(imgEdgeCurrent);

      // Get x, y derivatives.
      Mat dxCurrent, dyCurrent;
      Sobel(gray, dxCurrent, CV_32F, 1, 0, CV_SCHARR);
      Sobel(gray, dyCurrent, CV_32F, 0, 1, CV_SCHARR);

      // Use max-norm.
      dx      = max(dx, dxCurrent);
      dy      = max(dy, dyCurrent);
      imgEdge = max(imgEdge, imgEdgeCurrent); // TODO: Non-maximum suppression.
    }

  // Get orientation of gradient.
  int numAngBins = 10; 
  int numRowBins = 2; // # Vertical partitions

  histOrient     = Mat::zeros(numAngBins*numRowBins, 1, CV_32F); 
  fracEdgePixels = Mat::zeros(numRowBins, 1, CV_32F); 

  int avgNumRows = height/numRowBins;
  int numPixels  = (int) rgb.total();

  float maxEdgeMagCurrBlock = 0;
  for (int i=0; i<numRowBins; i++)
    {
      Rect rectImg;
      rectImg.x      = 0;
      rectImg.y      = i*avgNumRows;
      rectImg.width  = width;

      if (i == numRowBins-1) // Last partition.
        {
	  rectImg.height = avgNumRows + (height - numRowBins*avgNumRows);
        }
      else
        {
	  rectImg.height = avgNumRows;
        }

      if (rectImg.height<1)
	continue;

      Rect rectHist(0, i*numAngBins, 1, numAngBins);
      int numEdgePixels = 0;
       
      // FIX: Added to fix compilation issue on ubuntu
      const Mat cdx = dx(rectImg);
      const Mat cdy = dy(rectImg);
      const Mat cimgEdge  = imgEdge(rectImg);
      Mat   chistOrient   = histOrient(rectHist);
      GradientHist(numAngBins, cdx, cdy, cimgEdge, chistOrient, numEdgePixels, maxEdgeMagCurrBlock);
      if (maxEdgeMag < maxEdgeMagCurrBlock) maxEdgeMag = maxEdgeMagCurrBlock;
      fracEdgePixels.at<float>(i) = (float) numEdgePixels / (float) numPixels;
    }

  return;
}