/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * This program 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 2 of the License, or (at
 * your option) any later version.
 * 
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/*! \file ImageStats.hpp */
//! ImageStats.hpp

#ifndef _ImageStats_hpp
#define _ImageStats_hpp

//#include<ImLib3D/Image3D.hpp>
#include<ImLib3D/ImageProcessor.hpp>
#include<map>
#include<algorithm>


namespace IP3D
{
/*! #&Average&#     Computes the average and the variance of an images.
        Only for positions where mask is not 0.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Average       type:Im3DValue
     </UL>
~*/
	template <class ImageType>
	void Average(const ImageType& src,typename ImageType::value_type &average)
	{
		typename ImageType::const_iteratorFastMasked p;
		average=src.GetZero();
		int ct=0;
		for(p=src.begin();p!=src.end();++p){ct++;average+=*p;}
		average/=ct;
	}
};
namespace IP3D
{
/*! #&AverageAndVariance&#     Computes the average and the variance of an images.
        Only for positions where mask is not 0.
     <UL>Args:
      <LI>Source image  type:ImageType
      <LI>Average       type:double
      <LI>Variance      type:double
     </UL>
~*/
	template <class ImageType>
	void AverageAndVariance(const ImageType& src,double &average,double &variance)
	{
		typename ImageType::const_iteratorFastMasked p;
		average=0;variance=0;
		int ct=0;
		for(p=src.begin();p!=src.end();++p){ct++;average+=*p;variance+=(*p)*(*p);}
		average/=ct;
		variance=variance/ct-average*average;
	}
};

namespace IP3D
{
/*! #&RobustAverageAndVariance&#     Computes the average and the variance of an images using a robust technique.
        Robust (outlier rejection) technique is: median for average, 1.4826 median 
        absolute deviation for std dev.
     <UL>Args:
      <LI>Source image  type:Image3Df
      <LI>Average       type:double
      <LI>Variance      type:double
     </UL>
~*/
		void RobustAverageAndVariance(const Image3Df &src,double &average,double &variance);
};


//!  Median value of image, or more generaly: portion of histogram
/*!  
 */
namespace IP3D
{
/*! #&Median&#     Median value of an image, or more generaly: portion of histogram.
        This is NOT a median filter (see MedianFilter). For 0<=p<=1, Image I. 
        Computes value v of image for which Porba(I(P)<v)=p. For exmaple, for p=.5, 
        v=median value of I. Uses mask.
     <UL>Args:
      <LI>Source image.      type:ImageType
      <LI>Histogram portion  type:double
      <LI>Variance           type:double
     </UL>
~*/
	template <class ImageType>
	void Median(const ImageType& src,double &res,double part=.5/*! .5 for median*/,double precision=-.1)
	{
		
		vector<typename ImageType::value_type> values;
		values.reserve(src.GetNVoxels());
		for(typename ImageType::const_iteratorFastMasked p=src.begin();p!=src.end();++p)
		{
			values.push_back(*p);
		}
  		int n=values.size();
		sort(values.begin(),values.end());
//  		mprintf("CMedian result:%f %f %f %f %f\n",values[(int)(n*0)],values[(int)(n*.3)],
//  			   values[(int)(n*.5)],values[(int)(n*.66)],values[(int)(n-1)]);
		res=(double)values[(int)(n*part)];
	}
};

//! ProbabilityDistribution
class ProbabilityDistribution 
{
	typedef  float Value;
	vector<double> bins;
	float vmin,vmax;	

	// cumulative prob distribution
	bool cumulIsComputed;
	vector<double> cumul;
	bool cumulInverseIsComputed;
	map<double,int> cumulInverse;

	void SetRange(const vector<Value> &samples);

public:
	static int ScottsRule(vector<float>& vals); // for finding number of bins automatically

	vector<double> GetBins() { return bins;}	
	int GetNBins(){return bins.size();}
	void GetRange(float& rmin, float& rmax){rmin=vmin;rmax=vmax;}
	double GetMax(){return *max_element(bins.begin(), bins.end());}
	double GetMaxPos(){return BinValue(max_element(bins.begin(), bins.end())-bins.begin());}
	double operator[](int bin){return bins[bin];}
	double Density(int bin){return bins[bin]*bins.size()/(vmax-vmin);}
	//! Find value of begining of bin 
	double BinValue(int binNum)
	{
		return vmin+binNum*(vmax-vmin)/bins.size();
	}

	//! Like BinNUm with floating point result
	double BinPos(Value v)
	{
		return (bins.size()*(v-vmin)/(vmax-vmin));
	}

	//! Bin index number corresponding to value v
	int BinNum(Value v)
	{
		int res=(int)(BinPos(v));
		if(res>=(int)bins.size() || res<0){ThrowError("ProbabilityDistribution::BinNum out of range:%f-> %d (min=%f max=%f)\n",v,res,vmin,vmax);}
		return res;
	}
protected:
	void ComputeCumul()
	{
		cumul.clear();
		cumul.assign(bins.size(),0);
		double v=0;
		for(uint i=0;i<bins.size();i++)
		{
			v+=bins[i];
			cumul[i]=v;
		}
		cumulIsComputed=true;
	}
	void ComputeInverseCumul()
	{		
		cumulInverse.clear();
		if(!cumulIsComputed){ComputeCumul();}
		for(uint i=0;i<cumul.size();i++)
		{
			cumulInverse[cumul[i]]=i;
		}
		cumulInverseIsComputed=true;
	}
public:
	//! inverse of cumulative prob distrib
	double CumulInverse(double v)
	{
		if(!cumulInverseIsComputed){ComputeInverseCumul();}
		map<double,int>::iterator res=cumulInverse.upper_bound(v);
		if(res==cumulInverse.end()){return vmax;}
		int cbin=(*res).second-1;
		return cbin*(vmax-vmin)/bins.size()+vmin;
	}

	bool NextIndex(vector<int> &indexes)
	{
		for(uint i=0;i<indexes.size();i++)
		{
			indexes[i]++;
			if(indexes[i]<(int)bins.size()){return true;}
			indexes[i]=0;
		}
		return false;
	}

	class OtsuClass
	{
	public:
		int id;
		double moment0;
		double moment1;
		double moment2 ;
		int j0;
		int j1;
		OtsuClass():
			id(-1),moment0(0),moment1(0),moment2(0),j0(0),j1(0)
		{;}
		
	};

	OtsuClass ComputeClassStats(vector<OtsuClass> &classes,int classNum,vector<int> &indexes);
	void OtsuThreshholds(int nbClasses,vector<double> &res);
	void SetWithSamples(const vector<Value> &samples,int nbBins=1000);
	//! Like SetWithSamples, but uses linear interpolation for smoother distribution
	void SetWithSignal(const vector<Value> &samples,int nbBins=1000);
	//! Uses by default linear interpolation in 3 directions. type=1:faster type=2=fastest
	enum SetWithImage3DType {LINEAR3,LINEAR1,SAMPLES,FAST=LINEAR1,FASTEST=SAMPLES};
	void SetWithImage3D(const Image3Df &src,SetWithImage3DType type=LINEAR3,int nbBins=1000);
	void SetWithImage3Dlinear3(const Image3Df &src,int nbBins);
	void ShowWithGnuPlot();
	
	ProbabilityDistribution():cumulIsComputed(false),cumulInverseIsComputed(false)
	{}
};

#endif // _ImageStats_hpp

