/**********************************************************************************
*
* errorStatistics.h
*
* Purpose:  This class will compute various statistics of an image.
*           The statistics include: Truncated Average, Least Trimmed Squared, 
* 	    M-Estimation, Maximum, and 95% Maximum
* Input:    image :  this is the image to be analyzed
*
*	    *Note:  At this time, the images must be of floating type, 
*           templating can be done later to accomadate other types
*
*	    In addition, the user can set the percentage for robust HD statistics,
*	    and whether to trim from top or bottom in robust stats . The user also
* 	    has the ability to specify the region.The default is the
*	    largest possible region in the image. 
*	  
* Authors: Eric Billet, Andriy Fedorov 
* PI: Nikos Chrisochoides
**********************************************************************************/
/*
Copyright (c) 2008 College of William and Mary
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

   * Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
   * Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the
distribution.
   * Neither the name of the College of William and Mary nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef __errorStatistics_h
#define __errorStatistics_h

#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkImageConstIterator.h"
#include "itkNumericTraits.h"

template <typename TInputImageType, typename TEdge1>
class errorStatistics
{

	//Typedef statements//
	//const unsigned int Dimension;
        typedef TEdge1 Edge1Type;
        typedef TInputImageType InputImageType;

        //typedef itk::Image<float,3> InputImageType;

	typedef itk::ImageRegionIterator< InputImageType > IteratorType;
	typedef itk::ImageRegionIterator< Edge1Type > EdgeIteratorType;
	
	typedef itk::ImageRegionConstIterator< InputImageType > ConstIteratorType;
	typedef itk::ImageRegionConstIterator< Edge1Type > EdgeConstIteratorType;
	

	private:
		double m_AVG;
		double m_LTS;
		double m_M;
                double m_max;
		double m_Pmax;
		double m_RMS;
		double m_LMS;
		double m_mean;
		double m_Percentage;
		unsigned m_out;//
		char * outFile_name;//
		typename InputImageType::Pointer m_Image1;
		typename Edge1Type::Pointer m_EdgeImage1;
		typename InputImageType::RegionType m_Image1Region;	
		typename Edge1Type::RegionType m_EdgeImage1Region;
		long m_numMeasures;
		//	unsigned m_HDVectorSize;
		short trim_bottom;

	public:

		//Constructors
		errorStatistics();
		errorStatistics(typename InputImageType::Pointer image1, typename Edge1Type::Pointer edge_image);
		errorStatistics(typename InputImageType::Pointer image1,typename Edge1Type::Pointer edge_image,double percentage);
		

		//Functions to set parameters
		void SetImage1(typename InputImageType::Pointer image1);
		void SetImage1Region(typename InputImageType::RegionType region);

		void SetEdgeImage1(typename Edge1Type::Pointer image1);
		void SetEdgeImage1Region(typename Edge1Type::RegionType region);

		void SetPercentage(double p);

		void SetOutput(unsigned out){
			m_out = out;
		}
		void SetOutputFile(char * errorOut){
			outFile_name = errorOut;
		}
		void SetTrim_Bottom(short trim){
			trim_bottom = trim;
		}

		//Function to compute all of the statistics
		void ComputeValues();

		//Functions to return each statistic
		double GetMax(){
			return this->m_max;
		}
		double GetPmax(){
			return this->m_Pmax;
		}
		double GetAVG(){
			return this->m_AVG;
		}
		double GetLTS(){
			return this->m_LTS;
		}
		double GetM(){
			return this->m_M;
		}
		double GetRMS(){ 
			return this->m_RMS;
		}
		double GetLMS(){ 
			return this->m_LMS;
		}
		double GetMean(){ 
			return this->m_mean;
		}
		long GetnumMeasures(){
			return this->m_numMeasures;
		}
		//unsigned GetHDVectorSize();
};

//#ifdef ITK_MANUAL_INSTANTIATION
#include "errorStatistics.txx"
//#endif

/*
//Constructors:
template <typename TInputImageType, typename TEdge1>
errorStatistics<TInputImageType,TEdge1>::errorStatistics(): 
  m_Percentage(.85), m_out(0), trim_bottom(0){
  //m_HDVectorSize = 0;
}

template <typename TInputImageType, typename TEdge1>
errorStatistics<TInputImageType,TEdge1>::errorStatistics(typename InputImageType::Pointer image1, typename Edge1Type::Pointer edge_image): 
  m_Image1(image1),
  m_EdgeImage1(edge_image), m_out(0),
  m_EdgeImage1Region(m_EdgeImage1->GetRequestedRegion()),
  m_Image1Region(m_Image1->GetRequestedRegion()),trim_bottom(0){	
  // m_HDVectorSize = 0;
}

template <typename TInputImageType, typename TEdge1>
errorStatistics<TInputImageType,TEdge1>::errorStatistics(typename InputImageType::Pointer image1,typename Edge1Type::Pointer edge_image, double percentage): 
  m_Image1(image1), 
  m_EdgeImage1(edge_image), m_out(0),
  m_Percentage(percentage),
  m_EdgeImage1Region(m_EdgeImage1->GetRequestedRegion()),
  m_Image1Region(m_Image1->GetRequestedRegion()),trim_bottom(0){
  //m_HDVectorSize = 0;
}


//Member Functions:

//Set Functions:
template <typename TInputImageType, typename TEdge1>
void errorStatistics<TInputImageType,TEdge1>::SetImage1(typename InputImageType::Pointer image1){
	m_Image1 = image1;
}
template <typename TInputImageType, typename TEdge1>
void errorStatistics<TInputImageType,TEdge1>::SetEdgeImage1(typename Edge1Type::Pointer image1){
	m_EdgeImage1 = image1;
}
template <typename TInputImageType, typename TEdge1>
void errorStatistics<TInputImageType,TEdge1>::SetPercentage(double percentage){
	m_Percentage = percentage;
}
template <typename TInputImageType, typename TEdge1>
void errorStatistics<TInputImageType,TEdge1>::SetImage1Region(typename InputImageType::RegionType region){
	m_Image1Region = region;
}
template <typename TInputImageType, typename TEdge1>
void errorStatistics<TInputImageType,TEdge1>::SetEdgeImage1Region(typename Edge1Type::RegionType region){
	m_EdgeImage1Region = region;
}



//ComputeValues()

//////////////////////////////
//Function to get distance metrics
template <typename TInputImageType, typename TEdge1>
void  errorStatistics<TInputImageType,TEdge1>::ComputeValues(){

  //initialize metrics
  double RMS = 0.0;
	
  std::ofstream outFile(outFile_name);

  //First put all distances in a vector
  std::vector<double> allDistancesAcc;
  ConstIteratorType Iter1(m_Image1, m_Image1Region);
  EdgeConstIteratorType eIter(m_EdgeImage1,m_EdgeImage1Region);

  for(Iter1.GoToBegin(),eIter.GoToBegin();!Iter1.IsAtEnd();++Iter1, ++eIter){
    if(eIter.Get()>0){
      double dist = 0;
      dist = Iter1.Get();
	if(dist >= 0){
      	  RMS+=dist*dist;
	  //dist =sqrt(dist);
	  if(m_out)
	    outFile<<dist<<std::endl;
          allDistancesAcc.push_back(dist);
	}
    }
  }

  outFile.close();
  int i;

  double numMeasurements = allDistancesAcc.size();

  std::sort(allDistancesAcc.begin(),allDistancesAcc.end());

  // max
  double max = allDistancesAcc[(unsigned)(numMeasurements-1)];
	
  // Partial max metric
  double Pmax = allDistancesAcc[(unsigned)(numMeasurements*m_Percentage-1)];  

  double rec_perc = 1.-m_Percentage;
  double AvgDistance = 0.0;
  double LTSDistance = 0.0;
  double MDistance = 0.0;
  double mean = 0.0;

  //Get robust metrics
 for(std::vector<double>::const_iterator mI=allDistancesAcc.begin();
    mI!=allDistancesAcc.end();mI++)
    mean += *mI;
 mean /= numMeasurements;


  // IP-A(AVG) 
  // Trim the bottom and top of the measurement vector
  unsigned firstItem;
  firstItem = unsigned(numMeasurements*(1.-m_Percentage)/2.);
  for(std::vector<double>::const_iterator mI=allDistancesAcc.begin()+firstItem;
    mI!=allDistancesAcc.end()-firstItem;mI++)
    AvgDistance += *mI;
 
  AvgDistance /= numMeasurements*m_Percentage;

  // M-estimated Distance
  //
  double thresh = 0;
  if(trim_bottom){
    thresh = allDistancesAcc[unsigned(numMeasurements*(1.-m_Percentage))];
    for(std::vector<double>::const_iterator mI=allDistancesAcc.begin();
      mI!=allDistancesAcc.end();mI++){
        double distance = *mI;
        if(distance<thresh)
          MDistance += thresh;
        else
          MDistance += distance;
        }
    }
  else{
    thresh = allDistancesAcc[unsigned(numMeasurements*m_Percentage-1)];
    for(std::vector<double>::const_iterator mI=allDistancesAcc.begin();
      mI!=allDistancesAcc.end();mI++){
        double distance = *mI;
        if(distance>thresh)
          MDistance += thresh;
        else
          MDistance += distance;
     }
  }
  MDistance /= numMeasurements;

  //IP-LTS
  // Least Trimmed Squares
  //
  if(trim_bottom){
	unsigned lastItem;
  	lastItem = unsigned(numMeasurements*(1.-m_Percentage));
  	for(std::vector<double>::const_iterator mI=allDistancesAcc.begin()+lastItem;
	  mI!=allDistancesAcc.end();mI++){
    	    LTSDistance += *mI;
  	}
 
  LTSDistance /= numMeasurements*m_Percentage;
  }
  else{
	 unsigned lastItem;
  	lastItem = unsigned(numMeasurements*m_Percentage);
  	for(std::vector<double>::const_iterator mI=allDistancesAcc.begin();
    	  mI!=allDistancesAcc.begin()+lastItem;mI++){
    		LTSDistance += *mI;
  	}	
 
  LTSDistance /= numMeasurements*m_Percentage;
  }

  //LMS, least median squares
	double LMS=0.0;
	double min = 2000;
	
	unsigned endItem;
  	endItem = unsigned(numMeasurements*(.5));
	for(int x = 0;x<numMeasurements/2;x++){
		double residual = fabs(allDistancesAcc[x]-allDistancesAcc[x+endItem]);   	
		if(residual==min){
			//find midpoint
			double mid = (allDistancesAcc[x]+allDistancesAcc[x+endItem])/2.0;
			LMS +=mid;
			LMS/=2;
		}
		else if(residual < min){
			min = residual;
			LMS = (allDistancesAcc[x]+allDistancesAcc[x+endItem])/2.0;
		}

	}


  //Set private datamembers
  m_M = MDistance;
  m_LTS =LTSDistance;
  m_AVG = AvgDistance;
  m_Pmax = Pmax;
  m_max = max;
  m_RMS = RMS/numMeasurements;
  m_numMeasures = numMeasurements;
  m_LMS = LMS;
  m_mean = mean;
  return;
}
*/
  


#endif