/*! 
*   \file  CFeatureMapper.cpp
*   \brief Defines class to transform query signature.
*    
*   \details  
*   One approach to get complementary items is to map the input query signature
*   to the complementary space. This class facilitates this transformation. 
*
*   \date      July 17, 2012 
*   \copyright eBay Research Labs.
*/

// Include boost-specific headers.
#include "boost/regex.hpp"

// Include standard C++ headers.
#include <cmath>
#include <cstdlib> // For rand.
#include <iterator>
#include <algorithm>
#include <functional> // For bind1st
#include <iostream>
#include <sstream>
#include <fstream>


#include <numeric> // For accumulate.
#include <cstdlib> // For rand.
#include <ctime>   // For time.

// Include openCV-specific headers.
#include "opencv2/objdetect/objdetect.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <opencv2/contrib/contrib.hpp>

#include "CFeatureMapper.hpp"
#include "CFeatureExtraction.h"

// Namespaces.
using namespace std;
using namespace cv;
namespace bfs = boost::filesystem;
using namespace boost::lambda;

//#define SAVESWATCHFTRS

#define SQUAREDIM    30
#define SQUARESHIFT  (SQUAREDIM/2)

inline double drand()
{
    return (double) rand() / (double) RAND_MAX;
}

/*!
*   \brief Constructor.
*
*   \param [in] codebookApproach     Approach to derive code book. Default = 0 (Matsuda patterns).
*   \param [in] numParts             #parts to interpret.
*   \param [in] numWords             #words in vocabulary (code book). 
*/
CFeatureMapper::CFeatureMapper(int codebookApproach, int numParts, int numWords)
{
    m_numHueBins = 24; // TODO: Avoid this parameter.
    SelectCodeBookExtractionApproach(codebookApproach, numParts, numWords);
}


/*!
*   \brief Destructor.
*/
CFeatureMapper::~CFeatureMapper()
{

}


//**********************************************************************************************
// Public methods.
//**********************************************************************************************
 
/*!
*   \brief Learn mapping models.
*
*   \param [in]  contextRootDir  Image to be displayed.
*   \return      Boolean status for success (true = successful). 
*/
bool CFeatureMapper::LearnFeatureMappers(bfs::path contextRootDir)
{
    // TODO: Currently implements GMM model. 
    // This will be extended to take advantage of parameters set by SelectCodeBookExtractionApproach.

    // Get list of contexts. Context can be author, season, occasion.
    vector < pair < string, bfs::path > > contextList;
    GetSubdirectories(contextRootDir, contextList);
    int numContexts = (int) contextList.size();
    int numPatterns = (int) m_matsudaPatternLookup.size();

    // TODO: Multi-threaded.
    // Iterate through context directories.
    for (int i=0; i<numContexts; i++)
    {
        TickMeter timer;
        timer.start();

        // Images will be expected under "images" folder for each context. 
        // Annotation files will be expected under "annotations" folder for each context. 
        bfs::path imageRootDir       = contextList[i].second / "csImages"; //csImages
        bfs::path annotationsRootDir = contextList[i].second / "annotations";
        
#ifdef SAVESWATCHFTRS
        // Create directories to save feature files for swatches defined by each bounding box.
        vector <bfs::path> bboxFtrsDir;
        char suffix[8];
        for (int iPart=0; iPart<m_numParts; iPart++)
        {
            sprintf(suffix, "rect%d", iPart);
            bboxFtrsDir.push_back(contextList[i].second / suffix);
            bfs::create_directories(bboxFtrsDir[iPart]);
        }
#endif

        cerr << "Start learning model for: " << contextList[i].first << endl;
        cerr << "using annotations from: "   << annotationsRootDir   << endl;
        cerr << "and images from : "         << imageRootDir         << endl;

        // Get list of annotation files.
        vector < bfs::path > fileList;
        GetFilesInDirectory(annotationsRootDir, fileList, "txt");
        int numFiles = (int) fileList.size();

        // Extract features for each image. Pack each feature vector in a row of features.
        Mat featuresMat = cv::Mat::zeros(numFiles, m_numParts, CV_64FC1);        
        int numValid    = 0; // Maintain a count of valid feature vectors.

        Mat histBoW;
        if (m_numParts==1)
        {
            histBoW = cv::Mat::zeros(1, numPatterns, CV_32F);
        }
        else if (m_numParts==2)
        {
            histBoW = cv::Mat::zeros(numPatterns, numPatterns, CV_32F);
        }
        else
        {
            throw string("#Parts can only be 1 or 2"); 
        }

        // Iterate through annotation files.
        for (int j=0; j<numFiles; j++) 
        {
            if (j%1000 == 0)
                cerr << "File: " << j << ":: " << fileList[j] << endl;

            // Read annotation rectangles.
            vector <Rect> boundBoxes;
            if (!ReadAnnotationsFromFile(fileList[j], boundBoxes))
                continue; // Annotations wer not read successfully. Skip this entry.

            // Get image file name, based on annotation file name.
            // Name of annotation file for image named 0001.abc is 0001.abc.txt;
            string fileNameImage         = fileList[j].stem().string();        
            string fileNameImageFullPath = (imageRootDir / fileNameImage).string();
            if (!bfs::exists(fileNameImageFullPath))
            {
                continue; // Corresponding image file does not exist. Skip this entry.
            }

            // Read image.
            Mat img = imread(fileNameImageFullPath);
            if (img.empty())
            {
                continue; // Empty image. Skip this entry.
            }

            #ifdef _DEBUG
            PlotAnnotations(img, boundBoxes, "Read Image"); waitKey(500);
            #endif

            // Extract codewords for the selected image + annotations.
            vector <int> codeWords;
            if (!ComputeCodeWords(img, boundBoxes, codeWords))
            {
                continue; // Failed to extract codewords. Skip this entry.
            }                  

#ifdef SAVESWATCHFTRS
            // Extract feature for swatches defined by each bounding box. Also, save feature files.
            for (int rectIter = 0; rectIter < m_numParts; rectIter++)
            {
                string fileNameFtrs  = (bboxFtrsDir[rectIter] / fileNameImage).string();
                fileNameFtrs        += ".xml";
                Rect roi(boundBoxes[rectIter]);

                // Use SQUAREDIM, SQUARESHIFT
                if (roi.width<1)
                {
                    roi.x      = max(0, roi.x-SQUARESHIFT);
                    roi.y      = max(0, roi.y-SQUARESHIFT);
                    roi.width  = min(SQUAREDIM, img.cols-roi.x);
                    roi.height = min(SQUAREDIM, img.rows-roi.y);
                }

                CFeatureExtraction fx;
                fx.ExtractFeatures(img(roi), *this, false);
                fx.WriteFeatures(fileNameFtrs);
            }
#endif


            for( int partIter = 0; partIter < m_numParts; partIter++ )
                featuresMat.at<double> (numValid, partIter) = codeWords[partIter];

            // Accumulate histogram.
            if (m_numParts==1)
            {
                histBoW.at<float>(codeWords[0])++;
            }
            else
            {
                histBoW.at<float>(codeWords[1], codeWords[0])++;
            }
            //cerr << histBoW.at<float>(codeWords[1], codeWords[0]) << endl;
            numValid++;
        }
        GaussianBlur(histBoW, histBoW, Size(35, 35), 3);  // Smooth the histogram.

        // Take sqrt.
        MatIterator_<float> iter = histBoW.begin<float>();
        for (; iter != histBoW.end<float>(); iter++)
        {
            *iter = sqrt(*iter);
        }

        // Normalize.
        double maxVal;
        minMaxLoc(histBoW, 0, &maxVal, 0, 0);
        if (maxVal>0)
            histBoW /= maxVal;
         
        bfs::path fileNameFtrs(contextRootDir);
        fileNameFtrs /= "ftrs_" + contextList[i].first + ".xml";

        cerr << "*** " << numValid << " out of " << numFiles << " files were valid. ***\n" << endl;
        cerr << "Writing features to: " << fileNameFtrs.string() << endl;

        FileStorage fsFtrs(fileNameFtrs.string(), FileStorage::WRITE);
        fsFtrs << "f" << featuresMat;
        fsFtrs.release();

        // Fit GMM model using EM.
        cv::EM emTrainer(m_numWords, EM::COV_MAT_DIAGONAL, TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, EM::DEFAULT_MAX_ITERS, FLT_EPSILON) );
        emTrainer.train(featuresMat);
		cerr << "Finished learning model for: " << contextList[i].first << endl;
        

        // Create name for model file. Add prefix "modelComp_". 
        // Models files will be available in contextRootDir. They will be named after context name.
        bfs::path fileNameGMM(contextRootDir);
        fileNameGMM /= "modelComp_" + contextList[i].first + ".xml";

        FileStorage fs(fileNameGMM.string(), FileStorage::WRITE);
        emTrainer.write(fs);
        fs.release();

        
        // Write lookup table for codeWord occurrence.
        fs.open(fileNameGMM.string(), FileStorage::APPEND);
        fs << "LookupCodeWordOccurrence" << histBoW;
        if (countNonZero(histBoW) < histBoW.total())
            fs << "UseModel" << "GMM"; 
        else
            fs << "UseModel" << "LookupTable"; 

        fs.release();
        
        // REMOVE
        //FileStorage fsTest(fileNameGMM.string(), FileStorage::READ);
        //FileNode    gmmNodeTest = fsTest["StatModel.EM"];
        ////if (gmmNodeTest.empty())
        ////{
        ////    cerr << "Failed to Load Context Model from " << fileNameGMM << endl;
        ////    return false; 
        ////}

        //cv::EM   gmm;
        //gmm.read(gmmNodeTest);
        //if (gmm.isTrained())
        //{
        //    cerr << "Loaded Context Model from " << fileNameGMM << endl;
        //}
        //else
        //{
        //    cerr << "Failed to Load Context Model from " << fileNameGMM << endl;
        //    return  false;
        //}


        cerr << "Saved model to " << fileNameGMM << endl;
	timer.stop();
        cerr << "Total time = " << timer.getTimeMilli() << " ms\n" << endl << endl;
    } // for (int i=0; i<numContexts; i++)

    return true;
}


/*!
*   \brief Loads GMM models from the input XML files.
*
*   \param [in]  modelRootDir   Directory containing model files. All XML files in this 
*                               directory are assumed to be model files.
*   \return      Boolean status for success (true = successful). 
*/
bool CFeatureMapper::LoadFeatureMappers(bfs::path modelRootDir)
{
    // TODO: Merge all models into a single xml file. Currently, assume only one set of model files (GMM only or lookup table only).
    // TODO: This will be extended to take advantage of parameters set by SelectCodeBookExtractionApproach.
    
    // Reset.
    m_lookupCollection.clear();
    m_gmmCollection.clear(); 
    m_modelIDCollection.clear();

    vector < bfs::path > modelFiles;
    GetFilesInDirectory(modelRootDir, modelFiles, "xml");

    for (int i=0; i<modelFiles.size(); i++)
    {
        cerr << "Loading Context Model from " << modelFiles[i].string() << endl << flush;	  
        FileStorage fs(modelFiles[i].string(), FileStorage::READ);
        FileNode    gmmNode = fs["StatModel.EM"];

	
        if (gmmNode.empty())
        {
	  //	  cerr << "Warning: Empty GMM node, model :" << modelFiles[i] << endl << flush;
        }

        cv::EM   gmm;
        gmm.read(gmmNode);
        if (gmm.isTrained())
        {
	  cerr << "Loaded Context Model from " << modelFiles[i] << endl << flush;
        }
        else
        {
	  cerr << "Failed to Load Context Model from " << modelFiles[i] << endl << flush;
	  return  false;
        }

        m_gmmCollection.push_back(gmm);

        Mat histBoW;
        fs["LookupCodeWordOccurrence"] >> histBoW; 

        string modelString;
        fs["UseModel"] >> modelString; 

        int modelID = (strcmp(modelString.c_str(), "GMM")==0) ? 0 : 1; // 0 = GMM, 1 = table lookup.

        m_lookupCollection.push_back(histBoW);
        m_modelIDCollection.push_back(modelID);

        fs.release();        
    }   

	return true;
}


/*!
*   \brief Transforms the query signature to the space of given context.
*
*   \param [in, out]  querySig    Query signature. This will be modified to transformed query.
*   \param [in]       queryPartID 0-based PartID of query. TODO: Currently accepts only 0 or 1.
*   \param [in]       context     Context identifier for transformation. 
*                                 Accepted values: <1 (identity map), 1 (color science), >1 (data science)
*   \return           Boolean status for success (true = successful). 
*/
bool CFeatureMapper::TransformQuerySignatureUsingContext(CSignature& querySig, int queryPartID, int context)
{
    if (context<1)
        return true; // Identity map. Do nothing.

	vector < double > query, mappedQuery;
	float hueSum = 0;
    query.reserve(m_numHueBins);

    MatConstIterator_<float> iter = querySig.m_histColor.begin<float>(); //CV_32F
	// Copy from Mat to vector<float>. 
    // TODO: Avoid this step. Just use Mat.
    for(; iter != querySig.m_histColor.end<float>(); iter++)
	{
		query.push_back(*iter);
		hueSum += *iter;
	}
    Vec2i mappedCodeword;

#ifdef USEPATTERN
	TransformQuerySignatureUsingContext(query, queryPartID, querySig.m_fracColor, querySig.m_fracEdgePixels, querySig.m_maxEdgeMag, context, mappedQuery, mappedCodeword);
#endif // #ifdef USEPATTERN
	
    // TODO: Clean up code.
    if (context==3)
    {
        // Update complete signature (not just hue) from vector<float> to Mat. 
	    for( int histIter = 0; histIter < (int) mappedQuery.size(); histIter++ )
	    {
		    querySig.m_histColor.at<float> (histIter) = (float) mappedQuery[histIter];
	    }
        return true;
    }

    // Update signature from vector<float> to Mat. This updates hue.
	for( int histIter = 0; histIter < m_numHueBins; histIter++ )
	{
		querySig.m_histColor.at<float> (histIter) = (float) (mappedQuery[histIter] * hueSum);
	}

    // Update saturation, value (for color pixels) and value (for gray pixels).
    int compPartID = 1-queryPartID;

    int numBinsHue    = 24;
    int numBinsSatVal = 8;
    int numBinsVal    = 8;

    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 roiMat = querySig.m_histColor(clrSatROI);
    ConvertToComplementaryHistogram(roiMat, mappedCodeword[compPartID]);

    roiMat = querySig.m_histColor(clrValROI);
    ConvertToComplementaryHistogram(roiMat, mappedCodeword[compPartID]);

    roiMat = querySig.m_histColor(bwValROI);
    ConvertToComplementaryHistogram(roiMat,  mappedCodeword[compPartID]);

	return true;
}


/*!
*   \brief Compute codewords for the given image in regions defined by input rectangles.
*
*   \param [in]  src         Input image.
*   \param [in]  boundBoxes  Set of rectangular bounsing boxes to define swatches.
*   \param [out] codeWords   Set of codewords estimated for the input swatch regions.
*   \return  Boolean status for success (true = successful). 
*/
bool CFeatureMapper::ComputeCodeWords(Mat& src, vector <Rect>& boundBoxes, vector<int>& codeWords)
{
    int          selHist;
    int          channels[] = {0};
    int          hbins      = m_numHueBins;
    int          histSize[] = {hbins};    
    float        hranges[]  = {0, 180}; // hue varies from 0 to 179, see cvtColor
    const float* ranges[]   = {hranges};

    Mat hsv;
    cvtColor(src, hsv, CV_BGR2HSV); // Convert from BGR space to HSV.

    // Extract BoW for each part (bounding box).
    for (int rectIter = 0; rectIter < boundBoxes.size(); rectIter++)
    {
        Rect roi(boundBoxes[rectIter]);

        // Use SQUAREDIM, SQUARESHIFT
        if (roi.width<1)
        {
            roi.x      = max(0, roi.x-SQUARESHIFT);
            roi.y      = max(0, roi.y-SQUARESHIFT);
            roi.width  = min(SQUAREDIM, src.cols-roi.x);
            roi.height = min(SQUAREDIM, src.rows-roi.y);
        }

        // Use only the ROI as defined by the rectangle.
        Mat currHSV = hsv(roi);

        // Get Hue histogram.
        MatND histHue;
        calcHist(&currHSV, 1, channels, Mat(), // do not use mask
                 histHue, 1, histSize, ranges);

        vector< double > inputPattern;
		double histSum = 0;
        for( int binIter = 0; binIter < m_numHueBins; binIter++)
		{
            inputPattern.push_back( histHue.at<float> (binIter, 0) );
			histSum = histSum + histHue.at<float> (binIter, 0);
		}

		for( int binIter = 0; binIter < m_numHueBins; binIter++)
			inputPattern[binIter] = inputPattern[binIter] / histSum;

        selHist = FindClosestMatsudaPattern(inputPattern);
        codeWords.push_back(selHist);
    }
    return true;
}


//**********************************************************************************************
// Private methods.
//**********************************************************************************************
 
/*!
*   \brief Selects approach to derive code book.
*
*   \param [in] codebookApproach     Approach to derive code book. Default = 0 (Matsuda patterns).
*   \param [in] numParts             #parts to interpret.
*   \param [in] numWords             #words in vocabulary (code book). 
*/
void CFeatureMapper::SelectCodeBookExtractionApproach(int codebookApproach, int numParts, int numWords) 
{
    // TODO: Ignore interpreting codebookApproach for now. Matsuda patterns will be used.
    // ComputeVQCodebook, ComputeFeatureVQ

    m_codebookApproach = codebookApproach;
	m_numParts         = numParts;
    m_numWords         = numWords;
    m_numSectors       = 18; 

    CreateMatsudaCodebook();
}


 /*!
*   \brief Prepares book keeping variables for Matsuda codeword lookup.
*
*   \return Flag for success operation. 
*/
bool CFeatureMapper::CreateMatsudaCodebook()
{
	m_matsudaPatternLookup.clear();
	m_matsudaWeights.clear();

    const  int numPatterns = 7;

    // Pattern types: L, I, Y, X, i, V, T respectively.
    // Define starting locations (angle in degrees) of sector 1 and their widths.
    const int    start1[numPatterns]     = {333, 351, 351, 333, 351, 333, 270};
    const int    width1[numPatterns]     = { 54,  18,  18,  54,  18,  54, 180};

    // Define starting locations (angle in degrees) of sector 2 and their widths. 0 entries indicate that there is only one sector.
    const int    start2[numPatterns]     = { 81, 171, 153, 153,   0,   0,   0};
    const int    width2[numPatterns]     = { 18,  18,  54,  54,   0,   0,   0}; 
	const double inpWeights[numPatterns] = {1.2, 1.2, 1.2, 1.2,   1,   1,   1};
	

	//int  start1[7]  =      {351, 333, 333, 270, 351, 351, 333  };
    //int  start2[7]  =      {0,   0,   81, 0, 171, 153, 153};
    //int  width1[7]  =      {18,  54, 54, 180,  18,  18,  54};
    //int  width2[7]  =      {0,   0,  18,  0,  18,  54,  54 };
	//double inpWeights[7] = {.95,   .85,  .8, .5, .9,  .8,  .7 };
	//int numPatterns = 7;

    // Get bondaries of hue bins.
    vector <double>           hueBinEdges(m_numHueBins+1, 0); // Initalize with 0 values.    
    vector <double>::iterator iterBin = hueBinEdges.begin();
    double binWidth   = 360.0 / (double) m_numHueBins;
    double accumulate = 0;
    for (int bin=1; bin<m_numHueBins; bin++)
    {
        accumulate      += binWidth;
        hueBinEdges[bin] = accumulate;
    } 
    hueBinEdges[m_numHueBins] = 360; // Last bin.

    // Each Matsuda pattern (say "T") can occur with several variations (m_numSectors) by shifting the sectors.
    // Get representative histograms for each such pattern.
    double widthOfSector = 360.0 / (double) m_numSectors;
    for (int row=0; row<numPatterns; row++)
    {
        // Define sectors. There can be a maximum of 2 sectors for complementary recommendation. 
        // These values will be offset by width of sector.
        double s1 = start1[row];
        double w1 = width1[row];
        double s2 = start2[row];
        double w2 = width2[row];

        // Iterate for various shifts.
        for (int col=0; col<m_numSectors; col++)
        {
            vector <double> currentHist(m_numHueBins, 0); // Initalize with 0 values.

            // The first sector.
            IntersectMatsudaSectorWithHueHist(s1, w1, hueBinEdges, currentHist);
            s1 += widthOfSector;
            if (s1>360) s1 = s1-360; // Wrap around.

            if (w2>0)
            {
                // There are 2 sectors. 
                IntersectMatsudaSectorWithHueHist(s2, w2, hueBinEdges, currentHist);
                s2 += widthOfSector;
                if (s2>360) s2 = s2-360; // Wrap around.
            }

            double histSum = std::accumulate(currentHist.begin(), currentHist.end(), 0.0);
            if( histSum == 0 )
            {    
				cerr << "SOMETHING IS WEIRD. HISTSUM IS ZERO."<< endl;
				//cin.get();
			}
            else
            {
                // Normalize histogram.
                for( int iter = 0; iter < currentHist.size(); iter++)
			    {
                    currentHist[iter] /= histSum;
			    }
            }

            // Add representative histogram for future lookup.
			m_matsudaPatternLookup.push_back(currentHist);
			m_matsudaWeights.push_back(inpWeights[row]);
        }
    }

    return true;
}


/*!
*   \brief Maps the input signature to a Matsuda codeword.
*
*   \param [in] inputPattern   Input signature.
*   \return Matsuda codeword for the given signature. 
*/
int CFeatureMapper::FindClosestMatsudaPattern(vector <double>& inputPattern)
{
    double  inpHistSum   = 0;
    double  maxIntersect = 0;

	int numPatterns = (int) m_matsudaPatternLookup.size();
	int selHist     = 0;

	for( int iter = 0; iter < inputPattern.size(); iter++)
		inpHistSum = inpHistSum + inputPattern[iter];
	
	if( inpHistSum == 0 )
    {
		cerr << "Something is Wrong "<< endl;
        return -1;
    }
	
	for( int iter = 0; iter < inputPattern.size(); iter++)
		inputPattern[iter] = inputPattern[iter] / inpHistSum;

    for (int histIter = 0; histIter < numPatterns; histIter++)
    {
        double histInter = 0;
        for (int binIter = 0; binIter < m_numHueBins; binIter++)
			histInter = histInter + min(inputPattern[binIter], m_matsudaPatternLookup[histIter][binIter]);
		
		histInter = histInter * m_matsudaWeights[histIter];
        if( histInter > maxIntersect )
        {
            maxIntersect = histInter;
            selHist      = histIter;
        }
    }
    return selHist;
}


/*!
*   \brief Matsuda patterns contain sectors in hue wheel. This function represents such a sector as a hue hostogram.
*
*   \param [in]   start        Starting angle (degrees) of sector in hue wheel.
*   \param [in]   width        Angular width of sector.
*   \param [in]   hueBinEdges  Edge of hue bins in degrees.
*   \param [out]  overlap      Representative histogram for the given sector. Care is taken to handle franctional overlaps.
*/
void CFeatureMapper::IntersectMatsudaSectorWithHueHist(double start, double width, const vector <double> &hueBinEdges, vector <double> &overlap)
{
    int    numBins = (int) hueBinEdges.size() - 1;
    double scale   = (double) numBins / 360.0;

    //(int) (binWidth*value)
    double end = start+width;
    int i, i1, i2;
    i1 = (int) (scale*start); // Bin index where the starting point lies.
    if (end<360)
    {                
        i2 = (int) (scale*end); // Bin index where the ending point lies.

        // Non-boundary bins have 100% overlap with Matsuda sector.
        for (i=i1+1; i<i2; i++)
            overlap[i] = 1;

        // Boundary bins have partial overlap with Matsuda sector.
        overlap[i1] = (hueBinEdges[i1+1] - start) * scale;
        overlap[i2] = (end - hueBinEdges[i2]) * scale;
    }
    else
    {
        // Fill till 360 first.
        // Non-boundary bins have 100% overlap with Matsuda sector.
        for (i=i1+1; i<m_numHueBins; i++)
            overlap[i] = 1;   

        end  = end - 360;
        i2 = (int) (scale*end); // Bin index where the ending point lies.

        // Non-boundary bins have 100% overlap with Matsuda sector.
        for (i=0; i<i2; i++)
            overlap[i] = 1;   

        // Boundary bins have partial overlap with Matsuda sector.
        overlap[i1] = (hueBinEdges[i1+1] - start) * scale;
        overlap[i2] = (end - hueBinEdges[i2]) * scale;                                 
    }
}


/*!
*   \brief ContextID is passed through the main interface. This is internally mapped to a model ID.
*
*   \param [in]  context   Context identifier for transformation. 
*                          Accepted values: <1 (identity map), 1 (color science), >1 (data science)
*   \return      Mapped model ID. Here is the current interpretation of possible values:
*                0 - Similar, 1 - Color science, 2 - TAR, 3 - Rule-based (hybrid), 4 - Street, 
*                5 - Nicols, 6 - Paris
*/
int CFeatureMapper::ConvertContextIDToModelID(int context)
{
    return (context<4) ? 0 : context-4;
}


/*!
*   \brief Convert given histogram to its complement.
*
*   \param [in, out]  hgram      Histogram that needs to be mapped to its complement.
*   \param [in] mappedCodeWord   Codeword of desired (mapped) histogram.
*   \return  Boolean status for success (true = successful). 
*/
bool CFeatureMapper::ConvertToComplementaryHistogram(Mat &hgram, int mappedCodeWord)
{
    // Note: The current implementation uses Matsuda codewords.
    
    double sumOrig = (double) sum(hgram)[0];
    if (sumOrig<=0)
        return true; // Nothing to do.
     
    int patternID  = (mappedCodeWord<0) ? -1 : mappedCodeWord / m_numSectors; // ID of Matsuda pattern. There are 7 patterns: L, I, Y, X, i, v, T respectively.

    // Estimate number of dominant bins. TODO: Use EstimateNumDominantBins
    int   numDominant    = 0;
    float threshDominant = (float) 0.3 * (float) sumOrig; // Some heuristics to define "large" value. TODO: Needs more experimentation.    
    for (MatConstIterator_<float> iter = hgram.begin<float>(); iter != hgram.end<float>(); iter++)
    {
        if (*iter > threshDominant) numDominant++; // Count bins with "large" values.
    }

    // Get peak value of histogram.
    double maxVal;
    Point  maxLoc;
    minMaxLoc(hgram, 0, &maxVal, 0, &maxLoc);

    // Patterns L, I, Y, X have discontinuous spread of hue.
    // Patterns i, v, T have continuous spread of hue.
    bool  wideHue = (patternID >= 0) && (patternID < 4); // Indices 0, 1, 2, 3 correspond to L, I, Y, X respectively. 

    if (wideHue || (patternID<0 && numDominant>0))
    {
        // wideHue => recommend more variety.
        // (patternID<0 && numDominant>0) => There are few peaks. Recommend something with more variety. Take arithmetic complement.
        hgram = maxVal - hgram;
    }
    else
    {
        // Recommendation should have peaky distribution.
        hgram.at<float>(maxLoc.x) = (float) 1.5 * (float) maxVal; // Emphasize peak.
    }

    // Preserve original sum of elements.
    double sumModified = (double) sum(hgram)[0];
    hgram /= sumOrig/sumModified;

	return true;
}


/*!
*   \brief Estimate number of dominant bins in histogram. This is a coarse estimation.
*
*   \param [in]  queryFirst       Iterator pointing to the first entry of color histogram of query.
*   \param [in]  len              Number of bins in histogram.
*   \param [in]  scaleMaxThresh   SCale factor to scale th epeak value, to be used as threshold to define dominant bins.
*   \return      Estimated number of dominant colors.
*/
int  CFeatureMapper::EstimateNumDominantBins(const vector<double>::const_iterator queryFirst, int len, double scaleMaxThresh)
{
    // Get peak value.
    const  vector<double>::const_iterator maxValIter = max_element(queryFirst, queryFirst+len);
    double maxVal = *maxValIter;

    if (maxVal==0.0)
        return 0;

    // Count #bins above a threshold which is a fixed fraction of maximum value.
    // Distribution with tall peaks will have few such bins.
    double threshDominant = scaleMaxThresh*maxVal;
    return (int) count_if(queryFirst, queryFirst+len, bind1st(greater<double>(), threshDominant)); 
}


/*!
*   \brief Given a query, this function will classify whether it is solid/patterned and bw/colorful.
*
*   \param [in]  queryFirst           Iterator pointing to the first entry of color histogram of query.
*   \param [in]  len                  Length of color histogram.
*   \param [in]  queryFracClr         Fraction of colorful pixels to B/W pixels.
*   \param [in]  queryFracEdgePixels  Fraction of edge pixels.
*   \param [in]  queryMaxEdgeMag      Maximum edge magnitude in query. See CSignature.h
*   \param [in]  codeword             Codeword assigned to query signature. Example: Matsuda codeword.
*   \param [out] solid                Variable indicating whether input query is solid color or not.
*   \param [out] bw                   Variable indicating whether input query is b/w or colorful or mixed.
*/
void CFeatureMapper::IsSolidOrBW(const vector<double>::const_iterator queryFirst, int len, float queryFracClr, Mat &queryFracEdgePixels, float queryMaxEdgeMag, int codeword, int &solid, int &bw)
{
    // -----------------------------------
    //  Var/Val     -1      0       1
    // -----------------------------------
    //  BW          Color   Mixed   B/W
    //  Solid       Comp    Mixed   Solid
    // -----------------------------------

    // Estimate number of dominant bins for each subgroup (hue, sat, val, bw).
    int numDominantBins[4];
    numDominantBins[0] = EstimateNumDominantBins(queryFirst,    24);
    numDominantBins[1] = EstimateNumDominantBins(queryFirst+24,  8);
    numDominantBins[2] = EstimateNumDominantBins(queryFirst+32,  8);
    numDominantBins[3] = EstimateNumDominantBins(queryFirst+40,  8);

    // Classify whether solid color or not.
    float score = 0;
    if (queryFracEdgePixels.empty())
    {
        // Only color features are available.
        score = - 0.443726F*queryFracClr 
                - 0.140190F*numDominantBins[0] - 0.341343F*numDominantBins[1] - 0.063948F*numDominantBins[2] - 0.298221F*numDominantBins[3]  
                + 2.650507F;
    }
    else
    {
        // Both color and pattern features are available.
        score = 1.000314F*queryFracClr - 7.805258F*queryFracEdgePixels.at<float>(0) - 20.290347F*queryFracEdgePixels.at<float>(1) - 0.001610F*queryMaxEdgeMag 
                - 0.062721F*numDominantBins[0] - 0.223601F*numDominantBins[1] - 0.015656F*numDominantBins[2] - 0.107449F*numDominantBins[3]  
                + 3.767735F;
    }

    // Classify solid or not.
    solid = score>0 ? 1 : 0;


    // Classify b/w or not.
    // TODO: Use other factors such as: FracEdgePixels.
    if (queryFracClr>0.9)
    {
        // Colorful.
        bw    = -1;
    }
    else if (queryFracClr<0.1)
    {
        // Mostly B/W.
        //if (solid) // Reclassify.
        //{
        //    solid = (EstimateNumDominantBins(queryFirst+40, 8)<4) ? 1 : 0; // Just use bins corresponding to "Value" of B/W. 
        //}
        bw    = 1;
    }
    else
    {
        // Mixed.
        solid = 0;
        bw    = 0;
    }
}


/*!
*   \brief Given a histogram, this will generate a histogram with desired properties.
*
*   \param [in]   queryFirst         Iterator pointing to the first entry of histogram.
*   \param [in]   solid              Variable indicating whether input query is solid color or not. See IsSolidOrBW for meaning of values.
*   \param [in]   bw                 Variable indicating whether input query is bw or colorful. See IsSolidOrBW for meaning of values.
*   \param [out]  hgramOut           Output histogram. This is assumed to be initialized with 0 values.
*   \param [in]   bwbias             Intensity bias of BW channel. Bright when close to 1. This is ignored when negative. Range: [0, 1] or -1. Default = -1.
*   \param [in]   maxNumDominant     Desired maximum number of dominant colors. Default = 3.

*/
void CFeatureMapper::GenerateColorHistogram(const vector<double>& hgramIn, int solid, int bw, vector<double>& hgramOut, float bwbias, int maxNumDominant)
{
    vector<double>  hgramInCopy(hgramIn.begin(), hgramIn.end());
    vector<double>::const_iterator first = hgramInCopy.begin();

    int len = (int) hgramIn.size();

    // Initialize output with 0.
    vector <double> out;
    out.assign(len, 0); 

    vector <double>::iterator firstOut = out.begin();
    if (solid<0)
    {
        // To use complementary color. Copy only hue bins. TODO: May need to copy all bins.
        copy(hgramOut.begin(), hgramOut.begin()+24, hgramInCopy.begin());
    }

    switch (bw)
    {
        case -1: // Color.
            GenerateHistogram(first,   24, solid,    firstOut, -1, maxNumDominant); // Hue.
            GenerateHistogram(first+24, 8, solid, firstOut+24, -1, maxNumDominant); // Saturation.

            if (find_if(first+32, first+40, bind1st(greater<double>(), 0)) == first+40)
            {                            
                // No color pixels. Use B/W instead of Val. Reverse the order to get better complement (ex: dark against bright).
                vector<double>::const_reverse_iterator r = hgramInCopy.rbegin();
                vector<double> reverseVal(8);
                for (int i=0; i<8; i++)
                {
                    reverseVal[i] = *(r++);
                }
                GenerateHistogram(reverseVal.begin(), 8, solid, firstOut+32, -1, maxNumDominant);
            }
            else
            {
                GenerateHistogram(first+32, 8, solid, firstOut+32, -1, maxNumDominant); // Value.
            }
            break;

        case 0: // Mixed.  
            solid = 0;          
            GenerateHistogram(first,   24, solid, firstOut,        -1, maxNumDominant); // Hue.
            GenerateHistogram(first+24, 8, solid, firstOut+24,     -1, maxNumDominant); // Saturation.
            GenerateHistogram(first+32, 8, solid, firstOut+32,     -1, maxNumDominant); // Value.
            GenerateHistogram(first+40, 8, solid, firstOut+40, bwbias, maxNumDominant); // Value of BW.
            break;

        case 1: // B/W.
            // Update bwbias when recommendation is solid (1 or -1).
            
            if (solid!=0)
            {
                double sum = std::accumulate(first, first+len, 0.0);
                if (sum>0)
                {
                    // Get a measure of how "bright" the query image is.
                    double dark       = std::accumulate(first+32, first+36, 0) + std::accumulate(first+40, first+44, 0.0); 
                    double light      = std::accumulate(first+36, first+40, 0) + std::accumulate(first+44, first+48, 0.0);
                    double brightness = light / (dark+light);
                        
                    if (std::abs(brightness-bwbias)<0.05)
                    {
                        // Brightness of query image and bwbias are similar
                        if (std::abs(bwbias-0.5)<0.2)
                        {
                            // Medium brightness. Recommend pure black or pure white.
                            bwbias = rand() > RAND_MAX/2;
                        }
                        else
                        {
                            // When brightness of query image and bwbias are similar, 
                            // flip bwbias.
                            bwbias = 1-bwbias;
                        }
                    }
                } // if (sum>0)
            } // if (solid!=0)

            if (find_if(first+40, first+48, bind1st(greater<double>(), 0)) == first+48)
            {
                // No B/W pixels. Use Val instead of B/W.
                GenerateHistogram(first+32, 8, solid, firstOut+40, bwbias, maxNumDominant);
            }
            else
            {
                GenerateHistogram(first+40, 8, solid, firstOut+40, bwbias, maxNumDominant);
            }
            break;
    }

    hgramOut.assign(out.begin(), out.end());
}


/*!
*   \brief Given a histogram subset, this will generate a histogram with desired properties.
*
*   \param [in]   queryFirst         Iterator pointing to the first entry of histogram subset (say hue histogram).
*   \param [in]   len                Length of histogram subset.
*   \param [in]   solid              Variable indicating whether input query is solid color or not. See IsSolidOrBW for meaning of values.
*   \param [out]  hgramOut           Output histogram. This is assumed to be initialized with 0 values.
*   \param [in]   bwbias             Intensity bias of BW channel. Bright when close to 1. This is ignored when negative. Range: [0, 1] or -1. Default = -1.
*   \param [in]   maxNumDominant     Desired maximum number of dominant colors. Default = 3.

*/
void CFeatureMapper::GenerateHistogram(const vector<double>::const_iterator queryFirst, int len, int solid, vector<double>::iterator hgramOut, float bwbias, int maxNumDominant)
{
    // Note: hgramOut is assumed to be initialized with 0. sum(hgramOut) will be 1.

    const  vector<double>::const_iterator maxValIter = max_element(queryFirst, queryFirst+len);
    int    indMax = (int) distance(queryFirst, maxValIter); // Location of peak.
    double maxVal = *maxValIter;


    // Estimate number of dominant bins.
    double threshDominant = 0.25*maxVal; // Threshold used to estimate dominant bins.    
    vector <int> indDominant;
    for (int i=0; i<len; i++)
    {
        if (queryFirst[i]>threshDominant)
        {
            indDominant.push_back(i);
        }
    }
    int numDominant = (int) indDominant.size();

    double sum             = std::accumulate(queryFirst, queryFirst+len, 0.0); // Sum of histogram.
    double boostMax        = 2.0; // Bin value to boost the dominant bin.
    bool   needToNormalize = true; // Flag to indicate whether histogram needs to be normalized to have unit sum.
    if (solid==0)
    {
        // Mixed.
        // Start with random histogram, to give more variety. 
        // Limit the number of bins with positive values. Otherwise, recommended clothing is too colorful.
        int n0 = std::min(len, maxNumDominant + rand() % (len/2));
        for (int i=0; i<n0; i++)
        {
            hgramOut[rand() % len] = rand() / (double) RAND_MAX;
        }
        //generate_n(hgramOut, len, drand); // Pick bins at random.

        maxNumDominant = std::min(maxNumDominant, numDominant);
        if (maxNumDominant>0)
        {
            // Add bias to use randomly picked dominant bins.
            int numSelected = 1 + rand() % maxNumDominant;
            int count       = 0;
            // Select up to numSelected dominant dominant at random.
            for (int i=0; i<maxNumDominant; i++)
            {
                if (rand() %2)
                {
                    count++;
                    hgramOut[indDominant[i]] = boostMax * rand() / (double) RAND_MAX;
                    if (count==numSelected) break;
                }
            }
        }
    }
    else
    {
        // Solid. 1 = common color, -1 = complement color.
        if (bwbias>=0)
        {
            // Recommend solid B/W.
            if (rand() <= RAND_MAX /(1+numDominant)) // With equal chance.
            {
                // Pick intensity purely based on bwbias.
                indMax = std::min(len-1, (int) ceil(bwbias*len));
            }
            else
            {
                // Randomly pick one of the dominant bins instead of just max bin.
                indMax = indDominant[rand() % numDominant];
            }
        }
        else
        {
            // Recommend solid color (since bwbias<0).
            if (maxVal==0)
            {
                // Pick a random bin.
                indMax = rand() % len;
            }
            else
            {
                // Randomly pick one of the dominant bins instead of just max bin.
                indMax = indDominant[rand() % numDominant];
            }
        }    
        hgramOut[indMax] = (sum>0) ? sum : 1;
        needToNormalize  = false;  // Preserve sum. Do not normalize to unity.  
    }

    // Scale so that sum is 1.
    if (needToNormalize)
    {        
        double sumOut = std::accumulate(hgramOut, hgramOut+len, 0.0);
        if (sumOut>0 && sumOut!=1)
        {
            double scale = 1 / sumOut;
            for (int i=0; i<len; i++)
                hgramOut[i] *= scale;
        }
    }
}


/*!
*   \brief Maps the query codeword based on input context for the given query partID.
*
*   \param [in]  queryCodeWord    Query signature. This will be modified to transformed query.
*   \param [in]  queryPartID      0-based PartID of query. TODO: Currently accepts only 0 or 1.
*   \param [in]  context          Context identifier for transformation. 
*                                 Accepted values: <1 (identity map), 1 (color science), >1 (data science)
*   \param [in] mappedCodeWord    Codeword mapped to space defined by context and query part.
*   \return     Boolean status for success (true = successful). 
*/
bool CFeatureMapper::MapCodeWordVector(const Vec2i& queryCodeWord, int queryPartID, int context, Vec2i& mappedCodeWord)
{
    mappedCodeWord = queryCodeWord; // Initialize

    if (context<1) // Identity map. Nothing to do. Just return.
        return true;

    int compPartID      = 1-queryPartID;

    int indexCollection = ConvertContextIDToModelID(context);
    int modelID         = m_modelIDCollection[indexCollection]; // 0 = GMM, 1 = table lookup.
    if (modelID>0) 
    {
        // Use lookup table to map codeword.
        Mat    histBoW      = m_lookupCollection[indexCollection]; // Table lookup.
        int    compCodeWord = (int) queryCodeWord[queryPartID];    // Codeword of known part (query).

        // Determine codeword of unknown part (the complementary part) by picking 
        // the codeword with maximum entry in histBoW. 
        // Whether to search along row or column depends on partID of query.
        // Codewords for part 0 vary along columns (x-axis). Codewords for part 1 vary along columns (y-axis).
        Point  maxLoc;
        if (queryPartID==0)
        {
            // Search along column.
            minMaxLoc(histBoW.col(queryCodeWord[queryPartID]), 0, 0, 0, &maxLoc);
            compCodeWord = maxLoc.y;
        }
        else
        {
            // Search along row.
            minMaxLoc(histBoW.row(queryCodeWord[queryPartID]), 0, 0, 0, &maxLoc);
            compCodeWord = maxLoc.x;
        }
        mappedCodeWord[compPartID]  = compCodeWord;
    }
    else
    {
        // Use GMM for pairwise occurrence, to map codeword.
        Mat featuresMat = cv::Mat::zeros( 1, m_numParts, CV_64FC1 );
	    cv::EM emTrainer = m_gmmCollection[indexCollection];
	    Vec2d  samples;
	    int    numPatterns = (int) m_matsudaPatternLookup.size();
        int    maxIndex    = 0;
	    double maxlik      = -10000;
        double lik;

	    featuresMat.at<double> (0,0) = queryCodeWord[0];
        featuresMat.at<double> (0,1) = queryCodeWord[1];

	    for( int iter = 0; iter < numPatterns; iter++ )
	    {
		    featuresMat.at<double> (0,1) = iter;
		    samples = emTrainer.predict(featuresMat);
		    lik     = samples[0];
#ifdef GMM_LEARN
		    cerr << "For sample " << iter << " likelihood corresponding to " << query[0] << " is " << lik << endl;
#endif		
		    if( maxlik < lik )
            {
                maxlik = lik;
                maxIndex = iter;
            }
	    }
	    mappedCodeWord[compPartID]  = maxIndex;
#ifdef GMM_LEARN
	    cerr << "Input Matsuda is " << query[0] << "Mapped Matsuda is " << maxIndex << endl;
#endif
    }

    return true;
}


/*!
*   \brief Read annotation rectangles from the given text file.
*
*   \param [in]   fileName     Name of annotation file.
*   \param [out]  boundBoxes   List of rectangles.
*   \return       Boolean status for success (true = successful). 
*/
bool CFeatureMapper::ReadAnnotationsFromFile(boost::filesystem::path fileName, vector <Rect>& boundBoxes)
{
    // Example content of annotation file (only ineter:
    // 2                                    <= #people
    // 2 135 141 36 40 135 282 36 66 1      <= #rectangles followed by rectagles
    // 2 128 194 26 30 122 236 26 48 0      <= #rectangles followed by rectagles

    boundBoxes.clear();
    
    std::ifstream fs(fileName.string().c_str());
    string singleLine;

    getline(fs, singleLine);
    if (singleLine.empty())
        return false;

    int numPeople = atoi(singleLine.c_str()); 

    for (int i=0; i<numPeople; i++)
    {
        getline(fs, singleLine);
        if (singleLine.empty())
            return false;

        // Read line contents into a vector of integers. Note: This will crash if it encounters non-integers!
        istringstream iss(singleLine);
        vector <int>  values;
        copy(istream_iterator<int>(iss), istream_iterator<int>(), back_inserter(values));

        if (values[0] != m_numParts) // Each part needs a rectangle.
            return false;

        if (values.size() != 1+values[0]*4) // 4 = #parameter for rectangle.
            return false;

        // Valid entry. Add all rectangles.
        for (int j=0; j<values[0]; j++)
        {
            int k = 4*j+1;
            boundBoxes.push_back(Rect(values[k], values[k+1], values[k+2], values[k+3]));
        }
    }

    return true;
}


/*!
*   \brief Transforms the query signature to the space of given context.
*
*   \param [in, out]  querySig             Query signature. This will be modified to transformed query.
*   \param [in]       queryPartID          0-based PartID of query. TODO: Currently accepts only 0 or 1.
*   \param [in]       queryFracClr         Fraction of color pixels in query. See CSignature.h
*   \param [in]       queryFracEdgePixels  Fraction of edge pixels in query. See CSignature.h
*   \param [in]       queryMaxEdgeMag      Maximum edge magnitude in query. See CSignature.h
*   \param [in]       context              Context identifier for transformation. 
*                                          Accepted values: <1 (identity map), 1 (color science), >1 (data science)
*   \param [out]      mappedQuery          Query signature after mapping based on desired context.
*   \param [out]      mappedCodeword       Codeword of mapped signature.
*   \return           Boolean status for success (true = successful). 
*/
bool CFeatureMapper::TransformQuerySignatureUsingContext(vector<double>& query, int queryPartID, float queryFracClr, Mat &queryFracEdgePixels, float queryMaxEdgeMag, int context, vector<double>& mappedQuery, Vec2i &mappedCodeword)
{
    mappedQuery.clear();
    int compPartID    = 1-queryPartID;

    // Some initializations so that empty values are not returned.
    mappedCodeword[0] = -1;
    mappedCodeword[1] = -1;
    if(context==1)
	{
        // Color Science.

		//int firstHalt  = (int) ceil(0.33 * m_numHueBins);
        //int secondHalt = (int) ceil(0.66 * m_numHueBins);
		int halfHalt     = (int) ceil(0.50 * m_numHueBins);
        int iHalf;

        mappedQuery.assign(m_numHueBins, 0); // Initialize with 0.

		for(int binIter = 0; binIter < m_numHueBins; binIter++)
		{
			iHalf              = (binIter + halfHalt) % (int) m_numHueBins; // CircularIndex1Based( binIter + halfHalt, (int) m_numHueBins );
			mappedQuery[iHalf] = mappedQuery[iHalf] + query[binIter];
		}
	}
    else if( context == 2 ) 
    {       
        // TAR (Texture Agnostic Retrieval)

        Mat inpQuery = Mat::zeros((int) query.size(), 1, CV_64FC1);
        
        // This is an option to either sample subarrays (for instance H, S, and V sepatrately) or full arrays
        // (bag of words) ... here we are assuming color histograms and hence splitArray is set to true
        bool splitArray = true;
        
        // Conversion to matrix form for ease of use with OpenCV's random number generators
        for(int rowIter = 0; rowIter < inpQuery.rows; rowIter++)
        {
            inpQuery.at<double>(rowIter,0) = query[rowIter];
            mappedQuery.push_back(0);
        }
        Mat outQuery = Mat::zeros(inpQuery.rows, inpQuery.cols, CV_64FC1);
        cv :: theRNG () = cv :: RNG ( time (0) );
        randu(outQuery, Scalar(0), Scalar(1) );
        
        mappedQuery.assign(inpQuery.rows, 0); // Initialize with 0.
        // srand((unsigned int) time(0));
        srand(1);
        if( splitArray == true )
        {
            // Note: 
            // 24 Hue bins for color pixels.
            // 8 Sat bins for color pixels.
            // 8 Val bins for color pixels.
            // 8 Val bins for gray pixels.

            // Pick a random hue bin and make it 1.
            int ind = (int) ((double) rand() / (double) (RAND_MAX+1) * 24);
            mappedQuery[ind] = 1;

            // Pick a random sat bin and make it 1.
            ind = (int) ((double) rand() / ((double) RAND_MAX+1) * 8);
            mappedQuery[24+ind] = 1;

            // Pick a random val bin and make it 1.
            ind = (int) ((double) rand() / ((double) RAND_MAX+1) * 8);
            mappedQuery[32+ind] = 1;

            // Pick a random b/w val bin and make it 1.
            ind = (int) ((double) rand() / ((double) RAND_MAX+1) * 8);
            mappedQuery[40+ind] = 1;
        }
        else
        {
            // Pick a random bin and make it 1.
            int ind = (int) ((double) rand() / (double) (RAND_MAX+1) * query.size());
            mappedQuery[ind] = 1;
        }
    }
	else
	{
        // Data science.
		Vec2i queryCodeword;
		int selHist;
		// This must give out a 1 Dimensional Vector ... for example the input comes in "queryCodeword" for a shirt and the mapped
		// codeword now is a query for a pant
		selHist = FindClosestMatsudaPattern(query);	

		// This is the final mapped query.
		queryCodeword[queryPartID] = selHist;
		queryCodeword[compPartID]  = selHist;
		
		// Function call to the mapper
		MapCodeWordVector(queryCodeword, queryPartID, context, mappedCodeword);

        if (context==3)
            mappedQuery.assign((int) query.size(), 0); // Initialize with 0.

		// Return the histogram corresponding to mapped codeword
		for( int binIter = 0; binIter < m_matsudaPatternLookup[0].size(); binIter++)
		    mappedQuery.push_back( m_matsudaPatternLookup[ mappedCodeword[compPartID] ][ binIter ] );

        float bwbias         = (compPartID==0) ? 1.0F : 0.0F; // 0  = query image is for torso (ex: tops&blouses); 1 = query image is for legs (ex: pants, skirts)
        int   maxNumDominant = 2;
        if (context==3)
            TransformQuerySignatureUsingContextRule(query, queryPartID, queryFracClr, queryFracEdgePixels, queryMaxEdgeMag, context, mappedCodeword, mappedQuery, bwbias, maxNumDominant);
	}
	return true;
}


/*!
*   \brief Transforms the query signature to the space of given context, using a rule-based approach.
*
*   \param [in, out]  querySig             Query signature. This will be modified to transformed query.
*   \param [in]       queryPartID          0-based PartID of query. TODO: Currently accepts only 0 or 1.
*   \param [in]       queryFracClr         Fraction of color pixels in query. See CSignature.h
*   \param [in]       queryFracEdgePixels  Fraction of edge pixels in query. See CSignature.h
*   \param [in]       queryMaxEdgeMag      Maximum edge magnitude in query. See CSignature.h
*   \param [in]       context              Context identifier for transformation. 
*                                          Accepted values: <1 (identity map), 1 (color science), >1 (data science)
*   \param [out]      mappedQuery          Query signature after mapping based on desired context.
*   \param [out]      mappedCodeword       Codeword of mapped signature.
*   \param [in]       bwbias               Intensity bias of BW channel. Bright when close to 1. This is ignored when negative. Range: [0, 1] or -1. Default = -1.
*   \param [in]       maxNumDominant       Desired maximum number of dominant colors. Default = 3.
*   \return           Boolean status for success (true = successful). 
*/
void  CFeatureMapper::TransformQuerySignatureUsingContextRule(const vector<double>& query, int queryPartID, float queryFracClr, Mat &queryFracEdgePixels, float queryMaxEdgeMag, int context, const Vec2i &mappedCodeword, vector<double> &mappedQuery, float bwbias, int maxNumDominant)
{    
    bool  isColorMapped   = rand() > RAND_MAX/4;
    int   biasSolidThresh = (int) (0.3*RAND_MAX); // Prob(solid) = biasSolidThresh
    int   fracRand        = RAND_MAX/2;
    //int   compPartID      = 1-queryPartID;

    int solid;
    int bw          = 1;
    int solidMapped = 1;
    int bwMapped    = 1;

    IsSolidOrBW(query.begin(), (int) query.size(), queryFracClr, queryFracEdgePixels, queryMaxEdgeMag, mappedCodeword[queryPartID], solid, bw);

    switch (bw)
    {
        case -1: // Color.
            if (isColorMapped)
            {
                // Recommend items with color.
                bwMapped = -1;
                if (solid)
                {
                    // Any comp-solid color or mixed colors will work.     
                    if  (rand()>biasSolidThresh)                
                    {
                        solidMapped = 0;
                    }
                    else
                    {
                        bwMapped    = 0; // Note: There is a risk that complementary color is the same as dominant color. Hence use bwMapped = 0.
                        solidMapped = -1; 
                    }
                }
                else
                {
                    // Use solid color.
                    solidMapped = 1;
                }
            }
            else
            {
                // Recommend b/w items.
                bwMapped = 1;
                if (solid)
                {
                    // Any solid color or mixed colors will work.
                    solidMapped = rand()>biasSolidThresh ? 0 : 1;
                }
                else
                {
                    // Use any solid or comp-solid.
                    solidMapped = rand()>fracRand ? -1 : 1;
                }
            }
            break;

        //--------------------------------------------------------------------------------

        case 0: // Mixed.
            
            if (isColorMapped)
            {
                 // Recommend items with color.
                
                // Any solid.
                bwMapped    = -1;
                solidMapped = 1;
            }
            else
            {
                // Recommend b/w items.

                // Any comp-solid.
                bwMapped    = 1;
                solidMapped = -1;
            }
            break;

        //--------------------------------------------------------------------------------

        case 1: // B/W.
            if (solid==0) isColorMapped = false;
            if (isColorMapped)
            {
                // Recommend items with color.
                bwMapped  = -1;
                if (solid)
                {
                    // Any solid color or mixed colors will work.
                    solidMapped = rand()>biasSolidThresh ? 0 :  1;
                }
                else
                {
                    // Only solid will work.
                    solidMapped = 1;                    
                }
            }
            else
            {
                // Recommend b/w items.
                bwMapped = 1;

                if (solid)
                {
                    // Any comp-solid color or mixed colors will work.                     
                    if  (rand()>biasSolidThresh)                
                    {
                        solidMapped = 0;
                    }
                    else
                    {
                        bwMapped    = 0; // Note: There is a risk that complementary color is the same as dominant color. Hence use bwMapped = 0.
                        solidMapped = -1; 
                    }
                }
                else
                {
                    // Use Solid.
                    solidMapped = 1;
                }
            }
            break;
    }
    cerr << "Solid = " << solid << " BW = " << bw << endl;
    cerr << "SolidMapped = " << solidMapped << " BWMapped = " << bwMapped << endl;
    GenerateColorHistogram(query, solidMapped, bwMapped, mappedQuery, bwbias, maxNumDominant);
}

//**********************************************************************************************
// Utilities.
//**********************************************************************************************

/*!
*   \brief Displays the image along with annotated rectangles.
*
*   \param [in]  img         Image to be displayed.
*   \param [in]  boundBoxes  List of rectangles.
*   \param [in]  figName     Name of figure. Default = "Annotations".
*/
void CFeatureMapper::PlotAnnotations(Mat &img, vector <Rect>& boundBoxes, string figName)
{
    Mat temp = img.clone(); // Make a copy since retangle() will modify the image.
    for (int i=0; i<boundBoxes.size(); i++)
    {
        cerr << boundBoxes[i].x << ", " << boundBoxes[i].y << ", " << boundBoxes[i].width << ", " << boundBoxes[i].height << endl;
        rectangle(temp, Point(boundBoxes[i].x, boundBoxes[i].y), Point(boundBoxes[i].x+boundBoxes[i].width-1, boundBoxes[i].y+boundBoxes[i].height-1), Scalar(0, 0, 200), 2, 8); 
    }
    imshow(figName, temp);
}


/*!
*   \brief Given a 1-based index, this function will return the 1-based index with the given period.
*
*   \param [in]  position    1-based linear index.
*   \param [in]  period      Length of period.
*   \return  1-based periodic index.
*/
 int CFeatureMapper::CircularIndex1Based(int position, int period)
{
	int index = position % period;
	if( index == 0 )
		index = period;

	return index;
}


/*!
*   \brief Given a directory, this function will return the list of files.
*
*   \param [in]  rootDir     Root directory.
*   \param [out] fileList    Children of rootDir that are files.
*   \param [out] extension   Optional filter for file extension (case will be ignored). 
*                           Default filter returns only .jpg files.
*/
void CFeatureMapper::GetFilesInDirectory(bfs::path rootDir, vector < bfs::path > &fileList, string extension)
{
    // Note: Also see bfs::recursive_directory_iterator.

    fileList.clear(); // Initialize output.
    if( !bfs::exists(rootDir) || !bfs::is_directory(rootDir) )
    {
        // rootDir does not exist or is not a directory. 
        return;
    }

    // Use regex to apply pattern. Case-insensitive.
    boost::regex pattern(extension, boost::regex::icase);   

    // Look at children of rootDir.
    bfs::directory_iterator dir_iter(rootDir), dir_end;
    for( ; dir_iter != dir_end; ++dir_iter )
    {
        // Check if child is a file.
        if( bfs::is_regular_file(dir_iter->status()) )
        {
            //cerr << "********" << *dir_iter << ": " << dir_iter->path().filename() << "********" << endl;
            // Apply regex filter to file name.
            boost::smatch matches;
            if( boost::regex_search(dir_iter->path().filename().string(), matches, pattern) )
            {
                // Add this file to the list.
                fileList.push_back(*dir_iter);
            }            
        }
    }
}  // End of function: GetFilesInDirectory

/*!
*   \brief Given a directory, this function will return the list of immediate sub-folders.
*
*   \param [in]  rootDir     Root directory.
*   \param [out] subdirList  Children of rootDir that are directories.
*/
void CFeatureMapper::GetSubdirectories(bfs::path rootDir, vector < pair < string, bfs::path > > &subdirList)
{
    subdirList.clear(); // Initialize output.
    if( !bfs::exists(rootDir) || !bfs::is_directory(rootDir) )
    {
        // rootDir does not exist or is not a directory. 
        return;
    }

    // Look at children of rootDir.
    bfs::directory_iterator dir_iter(rootDir), dir_end;
    for( ; dir_iter != dir_end; ++dir_iter )
    {
        // Check if child is a directory.
        if( bfs::is_directory(*dir_iter) )
        {
            //cerr << "********" << *dir_iter << ": " << dir_iter->path().filename() << "********" << endl;
            string    name(dir_iter->path().filename().string());
            subdirList.push_back(make_pair(name, *dir_iter));
        }
    }
} // End of function: GetSubdirectories

//**********************************************************************************************

/*
int CFeatureMapper::ComputeVQCodebook( string& currDataDir )
{
    int numClusters = 50;
    int numDims     = 24;
    int NUM_FILES = 0;
    Mat hsv;
    int hbins = numDims, hSize[] = {hbins}, hchannels[] = {0};
    int sbins = numDims/3, sSize[] = {sbins}, schannels[] = {1};
    int vbins = numDims/3, vSize[] = {vbins}, vchannels[] = {2};
    
	
	float hrangesIn[] = { 0, 180 };
	float srangesIn[] = { 0, 1 };
	float vrangesIn[] = { 0, 1 };
    
	const float* hranges[] = { hrangesIn };
	const float* sranges[] = { srangesIn };
	const float* vranges[] = { vrangesIn };


    // Calculate the number of files in the directory
    for (bfs::directory_iterator itr( currDataDir.c_str() ); itr!=bfs::directory_iterator(); ++itr)
        if( !strcmp( (itr->path().extension().string()).c_str(), ".jpg") )
            NUM_FILES++;
    
    // Initialize feature matrix and bounding boxes for that directory
    Mat featureMatrix = cv::Mat::zeros( NUM_FILES, numDims, CV_64FC1 );
    Mat bestLabels = cv::Mat::zeros( NUM_FILES, 1, CV_64FC1 );
    Mat centers    = cv::Mat::zeros( numClusters, numDims, CV_64FC1 ); 
    // For each image, parse bounding box and compute bow feature vector
    int fileCtr = 0;
    for (bfs::directory_iterator itr( currDataDir.c_str() ); itr!=bfs::directory_iterator(); ++itr)
    { 
        if( !strcmp( (itr->path().extension().string()).c_str(), ".jpg") )
        {
            Mat inpImage = imread(itr->path().string(), -1);
            cvtColor(inpImage, hsv, CV_BGR2HSV);
            MatND hHist, sHist, vHist;
            calcHist( &hsv, 1, hchannels, Mat(), hHist, 1, hSize, hranges, true, false );
            calcHist( &hsv, 1, schannels, Mat(), sHist, 1, sSize, sranges, true, false );
            calcHist( &hsv, 1, vchannels, Mat(), vHist, 1, vSize, vranges, true, false );
			double hueSum=0, satSum=0, valSum=0;

			for( int binIter = 0; binIter < numDims; binIter++)
				hueSum = hueSum + hHist.at<float> (binIter, 0);
            for( int binIter = 0; binIter < numDims; binIter++)
                featureMatrix.at<double> (fileCtr, binIter) = hHist.at<float> (binIter, 0)/hueSum;

			for( int binIter = 0; binIter < sbins; binIter++)
				satSum = satSum + sHist.at<float> (binIter, 0);
			for( int binIter = 0; binIter < sbins; binIter++)
                featureMatrix.at<double> (fileCtr, binIter+sbins) = sHist.at<float> (binIter, 0)/satSum;

			for( int binIter = 0; binIter < vbins; binIter++)
				valSum = valSum + vHist.at<float> (binIter, 0);
			for( int binIter = 0; binIter < vbins; binIter++)
                featureMatrix.at<double> (fileCtr, binIter+sbins+vbins) = vHist.at<float> (binIter, 0)/valSum;
            fileCtr = fileCtr + 1;
        }
    }
    
    // Perform the Vector Quantizations
    kmeans(featureMatrix, numClusters, bestLabels, TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, EM::DEFAULT_MAX_ITERS, FLT_EPSILON), 10, KMEANS_RANDOM_CENTERS, centers );
    
    // Form a new representation for Matsuda Codebook
    for(int rowIter = 0; rowIter < featureMatrix.rows; rowIter++ )
    {
        vector< double > currentHist;
        for(int binIter = 0; binIter < numDims; binIter++)
            currentHist.push_back(0);
        double histSum = 0;
        for( int iter = 0; iter < featureMatrix.cols; iter++)
            histSum = histSum + featureMatrix.at<double> (rowIter, iter);
        if( histSum == 0 )
        {    
            cerr << "SOMETHING IS WEIRD HISTSUM IS ZERO"<< endl;
            cin.get();
        }
        for( int iter = 0; iter < currentHist.size(); iter++)
            currentHist[iter] = currentHist[iter]/histSum;
        m_matsudaPatternLookup.push_back(currentHist);
    }
    
    return 1;
}


// **********************************************************************************************
int CFeatureMapper::ComputeFeatureVQ(Mat& src, vector<int>& bowFeat, Mat& inpRects)
{
    int selHist;
	Mat hsv;
	int numDims = (int) m_matsudaPatternLookup[0].size();
	int hbins = numDims, hSize[] = {hbins}, hchannels[] = {0};
    int sbins = numDims/3, sSize[] = {sbins}, schannels[] = {1};
    int vbins = numDims/3, vSize[] = {vbins}, vchannels[] = {2};
    
	float hrangesIn[] = { 0, 180 };
	float srangesIn[] = { 0, 1 };
	float vrangesIn[] = { 0, 1 };
    
	const float* hranges[] = { hrangesIn };
	const float* sranges[] = { srangesIn };
	const float* vranges[] = { vrangesIn };

	Mat featureMatrix = Mat::zeros(m_numParts, numDims+numDims/3+numDims/3, CV_64FC1);

	cvtColor(src, hsv, CV_BGR2HSV);

    for(  int rectIter = 0; rectIter < m_numParts; rectIter++)
    {
        Mat currMask = cv::Mat::zeros(src.rows, src.cols, CV_8UC1), 
			visMask = cv::Mat::zeros(src.rows, src.cols, CV_8UC1);
			

        for( int rowIter = 0; rowIter < currMask.rows; rowIter++ )
            for( int colIter = 0; colIter < currMask.cols; colIter++ )
                if( colIter > inpRects.at<double> (rectIter,0) && colIter < (inpRects.at<double> (rectIter,0)+inpRects.at<double> (rectIter,2)) && rowIter > inpRects.at<double> (rectIter,1) && rowIter < (inpRects.at<double> (rectIter,1)+inpRects.at<double> (rectIter,3)) )
					currMask.at<unsigned char> (rowIter, colIter) = 1;
				
			
			MatND hHist, sHist, vHist;
            calcHist( &hsv, 1, hchannels, currMask, hHist, 1, hSize, hranges, true, false );
            calcHist( &hsv, 1, schannels, currMask, sHist, 1, sSize, sranges, true, false );
            calcHist( &hsv, 1, vchannels, currMask, vHist, 1, vSize, vranges, true, false );
			double hueSum=0, satSum=0, valSum=0;

			for( int binIter = 0; binIter < numDims; binIter++)
				hueSum = hueSum + hHist.at<float> (binIter, 0);
            for( int binIter = 0; binIter < numDims; binIter++)
                featureMatrix.at<double> (rectIter, binIter) = hHist.at<float> (binIter, 0)/hueSum;

			for( int binIter = 0; binIter < sbins; binIter++)
				satSum = satSum + sHist.at<float> (binIter, 0);
			for( int binIter = 0; binIter < sbins; binIter++)
                featureMatrix.at<double> (rectIter, binIter+sbins) = sHist.at<float> (binIter, 0)/satSum;

			for( int binIter = 0; binIter < vbins; binIter++)
				valSum = valSum + vHist.at<float> (binIter, 0);
			for( int binIter = 0; binIter < vbins; binIter++)
                featureMatrix.at<double> (rectIter, binIter+sbins+vbins) = vHist.at<float> (binIter, 0)/valSum;


        vector< double > inputPattern;
		double histSum = 0;
        for( int binIter = 0; binIter < m_numHueBins; binIter++)
            inputPattern.push_back( featureMatrix.at<double> (rectIter, binIter) );

        selHist = FindClosestMatsudaPattern(inputPattern);
        bowFeat.push_back(selHist);
    }
    return 1;
}

*/