/**
 * @file  RegionFeature.h
 * @brief A data colloection of some region to provide statistics for this region.
 *
 * Copyright (c) 2008, 2009, 2012, 2013 University of Pennsylvania.
 *
 * This file is part of DTI-DROID.
 *
 * DTI-DROID is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * DTI-DROID is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with DTI-DROID.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Contact: SBIA Group <sbia-software at uphs.upenn.edu>
 */

#ifndef _DTIDROID_REGIONFEATURE_H
#define _DTIDROID_REGIONFEATURE_H

#include <iostream>
#include <vector>
#include <string>

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>

#include "Histogram.h"


/***
  * class: ImageFeature
  * typename T - data type of the image feature
  * n - dimension of the image feature, or the size of member _data.
***/
template < typename T, int n>
class ImageFeature
{
public:
  ///@{
  /// default constructor
  ImageFeature() { _data = new T[n]; }
  /// copy constructor
  ImageFeature(const ImageFeature<T,n> & rhs)
  {
    _data = new T[n];
    memcpy(_data, rhs.getPtr(), n*sizeof(T));
  }
  
  ///destructor
  virtual ~ImageFeature() { if (_data) delete [] _data; _data=NULL; }
  ///@}
  
  /// set value for a specified position. default the first position
  bool set(T val, int i=0)
  {
    if (i>=n || i<0) { std::cout <<"WARNING: Out of range!" << std::endl; return false;}
    _data[i] = val;
    return true;
  }
  
  /// retrieve the data at specified position. default the first position
  T get(int i = 0) 
  {
    if (i>=n || i<0) { std::cout << "WARNING: Out of range!" << std::endl; return (T)0; }
    return _data[i];
  }
  
  T get(int i = 0) const 
  {
    if (i>=n || i<0) { std::cout << "WARNING: Out of range!" << std::endl; return (T)0; }
    return _data[i];
  }
  
  /// print out the data
  void print()
  {
    for (int i=0; i<n-1; i++)
      std::cout <<(int)(_data[i]) << ", ";
    std::cout << (int)(_data[n-1]) << std::endl;
  }
  
  /// retrieve the data pointer
  T* getPtr() const { return _data;}
  
private:
  T * _data;
};

typedef ImageFeature<float,1>           ScalarFeature;
typedef ImageFeature<float,3>           VectorFeature;
typedef ImageFeature<float,6>           TensorFeature;

/***
  * class: RegionFeature - abstract class for region feature
   * typename T - data type of the image feature
  * n - dimension of the image feature, or the size of member _data.
***/ 
template <typename T, int n>
class RegionFeature: public std::vector< ImageFeature<T,n> >
{
public:
  /// constructor
  RegionFeature() : bHistAvail(false), bStatAvail(false), _histValueLow(0), _histValueHigh(255) 
  {
    _nHistBins.clear();
    _nHistBins.resize(n, 1);
  }
  
  /***
    * function: compute the histogram for this region.
    * param: none
    * return: reference to histogram object.
  ***/
  void computeHistogram();
  
  /*** 
    * function: compute entropy for this region
    * param: none.
    * return: entropy value
  ***/
  double entropy();
  
  /***
    * function: compute geometric moment based on histogram
    * param: desired order for geometric moment
    * return: moment value
  ***/
  double moments(const std::vector<int>& orders); 

  /***
    * function: compute mutual information of 2 region features
    * param: rhs - right hand side region feature. Input region feature 
    *            is assumed one-to-one corresponding to this region feature
    *            in vector index.
    * return: mutual information, normalized to between [0,1]
  ***/
  /// regular MI = H(x)+H(y)-H(x,Y)
  double MI(RegionFeature<T,n>& rhs);
  
  /// normalized MI = 2 - 2(H(x,y)/(H(x)+H(y)))
  double NMI(RegionFeature<T,n>& rhs);

  /// K-L divergence
  double KL_DIVERGENCE(RegionFeature<T,n>&rhs);
  
  /***
    * function: compute moment distance given the order
    * param:  rhs - right hand side region feature. 
    * param: desired order for geometric moment
    * return: distance, normalized between [0,1]
  ***/
  double GM_DIST(RegionFeature<T,n>& rhs, const std::vector<int>& orders);
  
  /// virtual method: compute the statistical information for this region
  virtual void statistics();

  /// print histogram
  void printHistogram()
  {
    if (!bHistAvail) return;
    _regionHistogram.print();
  }
  
  void printHistogram() const
  {
    if (!bHistAvail) return;
    _regionHistogram.print();
  }
  
  /// access the region statistics data
  const Histogram<T,n>& getHistogram() 
  { 
    if (!bHistAvail) computeHistogram();
    return _regionHistogram;
  }
  
  const ImageFeature<T,n>& getMeanFeature() 
  {
    if (!bStatAvail) statistics();
    return _regionMeanFeature;
  }
  
  const ImageFeature<T,n>& getFeatureVariance()
  {
    if (!bStatAvail) statistics();
    return _regionFeatureVariance;
  }
  
  /// set&get histogram value range
  void setHistogramRange(T valLow, T valHigh)
  { _histValueLow = valLow; _histValueHigh = valHigh; }
  
  void getHistogramRange(T *valLow, T *valHigh)
  { *valLow = _histValueLow; *valHigh = _histValueHigh; }
  
  /// set&get histogram bins
  void setHistBins(int nBins) { for (int i=0; i<n; i++) _nHistBins[i]=nBins; }
  void setHistBins(const std::vector<int>& nBins) { _nHistBins = nBins; }
  const std::vector<int>& getHistBins() { return _nHistBins; }
  const std::vector<int>& getHistBins() const { return _nHistBins; }
  
private:
  
  bool bHistAvail;
  bool bStatAvail;
  
  std::vector<int> _nHistBins;
  T _histValueLow, _histValueHigh;
  
  Histogram<T,n> _regionHistogram;

protected: 
 
  ImageFeature<T,n> _regionMeanFeature;
  
  ImageFeature<T,n> _regionFeatureVariance;

};


/***
  * class - VectorRegionFeature
***/
template<typename T>
class VectorRegionFeature: public virtual RegionFeature<T,3>
{
public: 
  /// compute a principal direction of this vector region using PCA
  void computePD();
  
  /// get the principal direction of the region
  const VectorFeature & getPD() const {return _regionPD;}
  
private:
  /// principal direction of the region
  VectorFeature _regionPD;
};

//==================I am a separator:)=========================//

/***
  * Implementation: class - RegionFeature
***/
template <typename T, int n>
void RegionFeature<T,n>::statistics()
{
  for (int d=0; d<n; d++)
  {
    T mean = (T)0.0;
    for (int i=0; i<this->size(); i++)
      mean += this->at(i).get(d);
  
    if (this->size() > 0)  
      mean = (T)((float)mean / (float)this->size());

    _regionMeanFeature.set(mean, d);
  
    T var = (T)0.0;
    for (int i=0; i<this->size(); i++)
      var += (this->at(i).get(d)-mean)*(this->at(i).get(d)-mean);
    
    if (this->size()>1)
      var = (T)((float)var / (float)(this->size()-1));
    
    _regionFeatureVariance.set(var,d);
  }
}

template <typename T, int n>
void RegionFeature<T,n>::computeHistogram()
{
  /// compose a data vector for building histogram
  std::vector< std::vector<T> > histData(n);
  for (int d=0; d<n; d++)
  {
    //cout << "region size: " << size() << endl;
    histData[d].resize(this->size());
    for (int i=0; i<this->size(); i++)
    {
      histData[d][i] = this->at(i).get(d);
    }
  }
  
 // cout << "histData size: " << histData[0].size() << endl;
  
  _regionHistogram.resize(_nHistBins);
  _regionHistogram.setValueRange(_histValueLow, _histValueHigh);
  _regionHistogram.build(histData);
  
  bHistAvail = true;
}

template <typename T, int n>
double RegionFeature<T,n>::entropy()
{
  if (!bHistAvail) computeHistogram();
  
  return _regionHistogram.entropy();
}

template <typename T, int n>
double RegionFeature<T,n>::moments(const std::vector<int>& orders)
{
  if (!bHistAvail) computeHistogram();
  
  _regionHistogram.setRemoveMomentDC(true);
  return _regionHistogram.moments(orders);
}

template <typename T, int n>
double RegionFeature<T,n>::MI(RegionFeature<T,n>& rhs)
{
  if (rhs.size() != this->size())
  {
    std::cout << "ERROR: right-hand-side feature is not corresponding to current feature!" << std::endl;
    return 0.0;
  }
  
  double thisEntropy = entropy();
  double rhsEntropy = rhs.entropy();
  
  /// calculate a joint histogram for these two regions
  Histogram<T,2*n> jointHist;
  
  T valLow, valHigh;
  _regionHistogram.getValueRange(&valLow, &valHigh);
  jointHist.setValueRange(valLow,valHigh);
  
  std::vector<int> nJointHistBins(_nHistBins);
  std::vector<int> rhsHistBins = rhs.getHistBins();
  
  for (int i=0; i<rhsHistBins.size(); i++)
    nJointHistBins.push_back(rhsHistBins[i]);

  //for (int i=0; i<nJointHistBins.size(); i++)
  //  cout << nJointHistBins[i] << "\t" << endl;
    
  std::vector< std::vector<T> > histData(2*n);
  for (int d=0; d<n; d++)
  {
    histData[d].resize(this->size());
    for (int i=0; i<this->size(); i++)
    {
      histData[d][i]=(this->at(i).get(d));
    }
  }
   
  /// the input feature is one-to-one corresponding to current feature, in vetor index
  for (int d=0; d<n; d++)
  {
    histData[n+d].resize(rhs.size());

    for (int i=0; i<rhs.size(); i++)
    { 
      histData[n+d][i]=rhs.at(i).get(d);
    }
  }
  
  jointHist.resize(nJointHistBins);
  jointHist.build(histData);
  //jointHist.print();
  
  double jointEntropy = jointHist.entropy();
  
  //*** DEBUG ***
  /*
  cout << "jointHistogram:" << endl;
  jointHist.print();
  cout << "jointEn: " << jointEntropy << " En A: " << thisEntropy << " En B: " << rhsEntropy << endl;
  */
  //*** DEBUG***
  
  double mi = thisEntropy+rhsEntropy-jointEntropy;
  return mi;
}

template <typename T, int n>
double RegionFeature<T,n>::NMI(RegionFeature<T,n>& rhs)
{
  if (rhs.size() != this->size())
  {
    std::cout << "ERROR: right-hand-side feature is not corresponding to current feature!" << std::endl;
    return 0.0;
  }
  
  double thisEntropy = this->entropy();
  
  double rhsEntropy = rhs.entropy();
  
  /// calculate a joint histogram for these two regions
  Histogram<T,2*n> jointHist;
  
  T valLow, valHigh;
  _regionHistogram.getValueRange(&valLow, &valHigh);
  jointHist.setValueRange(valLow,valHigh);
  
  std::vector<int> nJointHistBins(_nHistBins);
  std::vector<int> rhsHistBins = rhs.getHistBins();
  
  for (int i=0; i<rhsHistBins.size(); i++)
    nJointHistBins.push_back(rhsHistBins[i]);

  //for (int i=0; i<nJointHistBins.size(); i++)
  //  cout << nJointHistBins[i] << "\t" << endl;
    
  std::vector< std::vector<T> > histData(2*n);
  for (int d=0; d<n; d++)
  {
    histData[d].resize(this->size());
    for (int i=0; i<this->size(); i++)
    {
      histData[d][i]=(this->at(i).get(d));
    }
  }
   
  /// the input feature is one-to-one corresponding to current feature, in vetor index
  for (int d=0; d<n; d++)
  {
    histData[n+d].resize(rhs.size());

    for (int i=0; i<rhs.size(); i++)
    { 
      histData[n+d][i]=rhs.at(i).get(d);
    }
  }
  
  jointHist.resize(nJointHistBins);
  jointHist.build(histData);
  //jointHist.print();
  
  double jointEntropy = jointHist.entropy();
  
  //*** DEBUG ***
  //std::cout << "this circular region data:" << std::endl;
  //for (int i=0; i<size(); i++)
  //  at(i).print();
    
  //std::cout << "thisHistogram:" << std::endl;
  //printHistogram();
  //std::cout << "rhsHistogram:" << std::endl;
  //rhs.printHistogram();
  //std::cout << "jointHistogram:" << std::endl;
  //jointHist.print();
  //std::cout << "jointEn: " << jointEntropy << " En A: " << thisEntropy << " En B: " << rhsEntropy << std::endl;
  
  //*** DEBUG***
  
  double mi = 2.0-2.0*jointEntropy/(thisEntropy+rhsEntropy+1e-24);
  if (mi > 1.0) mi = 1.0;
  return mi;
}

template <typename T, int n>
double RegionFeature<T,n>::KL_DIVERGENCE(RegionFeature<T,n>&rhs)
{
  if (n>1) return -1.0;
  
  if (rhs.size() != this->size())
  {
    std::cout << "ERROR: right-hand-side feature is not corresponding to current feature!" << std::endl;
    return -1.0;
  }
  
  // compute cross entropy
  Histogram<T,n> rhsHist = rhs.getHistogram();
  std::vector<int> rhsHistBins = rhs.getHistBins();
  
  if (_nHistBins[0] != rhsHistBins[0])
  {
    std::cout << "ERROR: the histogram bins are different!" << std::endl;
    return -1.0;
  }
  double croEn = 0.0;
  for (int i=0; i<_nHistBins[0]; i++)
    croEn += (-_regionHistogram[i]*log(rhsHist[i]+1e-20)/log(2.0));
    
  double en = _regionHistogram.entropy();
  
  return croEn -en;
}

template <typename T, int n>
double RegionFeature<T,n>::GM_DIST(RegionFeature<T,n>& rhs, const std::vector<int>& orders)
{
  if (rhs.size() != this->size())
  {
    std::cout << "ERROR: right-hand-side feature is not corresponding to current feature!" << std::endl;
    return 0.0;
  }
  
  if (orders.size() != n)
  {
    std::cout << "ERROR: dimension of specified orders wrong!" << std::endl;
    return 0.0;
  }
  
  double m1 = moments(orders);
  double m2 = rhs.moments(orders);
  
  double norm = 1.0;
  for (int i=0; i<n; i++)
    norm *= pow((double)_nHistBins[i], (double)orders[i]);

  //cout << "m1 = " << m1 << "; m2 = " << m2 << "; norm = " << norm << endl;
  
  return fabs(m1-m2)/norm;
}


/***
  * Implementation: class - VectorRegionFeature
***/
template<typename T>
void VectorRegionFeature<T>::computePD()
{
  // compute a covariance matrix
  this->statistics();
  
  double mean[3];
  mean[0] = this->_regionMeanFeature.get(0);
  mean[1] = this->_regionMeanFeature.get(1);
  mean[2] = this->_regionMeanFeature.get(2);
  
  double M[9];
  M[0] = this->_regionFeatureVariance.get(0);
  M[4] = this->_regionFeatureVariance.get(1);
  M[8] = this->_regionFeatureVariance.get(2);
  
  M[1] = 0.0f; M[2] = 0.0f; M[5] = 0.0f;
  for (int i=0; i<this->size(); i++)
  {
    M[1] += (this->at(i).get(0)-mean[0])*(this->at(i).get(1)-mean[1]);
    M[2] += (this->at(i).get(0)-mean[0])*(this->at(i).get(2)-mean[2]);
    M[5] += (this->at(i).get(1)-mean[1])*(this->at(i).get(2)-mean[2]);
  }
   
  if (this->size()>1)
  {
    M[1] /= (double)(this->size()-1);
    M[2] /= (double)(this->size()-1);
    M[5] /= (double)(this->size()-1);
  }
  
  M[3]=M[1]; M[6]=M[2]; M[7]=M[5];
 
  // eigen decomposition
  gsl_vector *eval ;
  gsl_matrix *evec ;
  gsl_eigen_symmv_workspace  *w ;
  
  eval = gsl_vector_alloc(3);
  evec = gsl_matrix_alloc(3, 3);
  w = gsl_eigen_symmv_alloc(3);
  
  gsl_matrix_view m = gsl_matrix_view_array(M, 3, 3);
  gsl_eigen_symmv(&m.matrix, eval, evec, w);
  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
  
  // retrieve PD
  double lb = gsl_vector_get(eval, 0);
  gsl_vector_view evec0 = gsl_matrix_column(evec, 0);
  for (int i=0; i<3; i++)
  {
    double val = gsl_vector_get(&evec0.vector, i);
    _regionPD.set((float)(lb*val), i);
  }
}


#endif // _DTIDROID_REGIONFEATURE_H
