/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001 Topi Mäenpää and Jaakko Viertola
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.10 $
 *********************************************************************/

#ifndef _THRESHOLDING_H
#define _THRESHOLDING_H

#include "../ImageTransform.h"
#include <Matrix.h>
#include <Math.h>
#include <math.h>

namespace prapi { namespace texture {

	/**
	 * The base class for thresholding which allows user give
	 * it's own threshold.
	 **/
  template <class T> class Thresholding: public ImageTransform<T,T>
	{
	public:
		/**
		 * Constructor of thresholding.
		 *
		 * @param lowerThreshold is the lower threshold.
		 * @param upperThreshold is the higher threshold.
		 *        Note. If the upperThreshold is zero it's the same as
		 *        there were no upper threshold.
		 **/
		Thresholding(T lowerThreshold,T upperThreshold=0) :
			_lowerThreshold(lowerThreshold),_upperThreshold(upperThreshold){}

		/**
		 * The destructor of thresholding.
		 **/
		~Thresholding(){}

		/**
		 * Makes the thresholding for the picture by given threshold.
		 *
		 * @param mat which wanted to use thresholding.   
		 **/
		util::Matrix<T> getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&);

		/**
		 * Gives the value of threshold.
		 **/
		T getLowerThreshold(void) const {return _lowerThreshold;}
		T getUpperThreshold(void) const {return _upperThreshold;}

		/**
		 * Set the new value for threshold.
		 **/
		void setLowerThreshold(T threshold){_lowerThreshold=threshold;}
		void setUpperThreshold(T threshold){_upperThreshold=threshold;}

	protected:
		/**
		 * Constructor used if need to calculate threshold before
		 * it will be used.
		 **/
		Thresholding() : _lowerThreshold(0),_upperThreshold(0) {}
		/**
		 * _threshold is the value of threshold.
		 **/
		T _lowerThreshold;
		T _upperThreshold;
	};

	/**
	 * MSEdgedetector calculates the threshold by formula sum(Gx"2+Gy"2)*2/sqrt((row-1)*(col-1))
	 * where row and col is the amount of rows and columns.
	 **/
	template <class T> class MSEdgeThresholding: public Thresholding<T>
	{
	public:
		/**
		 * Constructor which calculates the threshold from the given matrix.
		 *
		 * @param matX Gradient matrix X.
		 * @param matY Gradient matrix Y.
		 **/
		MSEdgeThresholding(const util::Matrix<T>& matX, const util::Matrix<T>& matY);
	};

	template <class T> util::Matrix<T> Thresholding<T>::getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&)
	{
		// go throw every pixel in the matrix and check if its high enough.
		util::Matrix<T> result(mat);
		T* data = result.getData();

		// if its not high enough set the valut T 0 (typecast)
		for(int i=0;i<result.getRows()*result.getColumns();i++,data++)
			if(*data<_lowerThreshold || *data>_upperThreshold)*data = T(0);

		return result;
	}

	template <class T> MSEdgeThresholding<T>::MSEdgeThresholding(const util::Matrix<T>& matX, const util::Matrix<T>& matY) :
		 Thresholding<T>()
	{
		int rowsX = matX.getRows();
		int columnsX = matX.getColumns();
		if(rowsX != matY.getRows() || columnsX != matY.getColumns())
			throw util::MatrixException("MSEdgeThresholding<T>::MSEdgeThresholding(const util::Matrix<T>&,const util::Matrix<T>&): The size of matX and matY differs.");

		const T* matXData = matX.getData();
		const T* matYData = matY.getData();
		T squareSum=0;
		
		// first calculate the sum
		for(int i=0;i<rowsX*columnsX;i++,matXData++,matYData++)
			squareSum += T(pow(double(*matXData),2)+pow(double(*matYData),2)); 

		// then multipilacate by 2 and divide by
		squareSum = T(double(squareSum*2)/sqrt(double(rowsX-1)*double(columnsX-1)));

		// set the new value for threshold in base class Thresholding
		Thresholding<T>::_lowerThreshold = squareSum;
	}

	/**
	 * HysteresisThersholding of edge pixels. Starting at pixels with a
	 * value greater than the upperThreshold, trace a connected sequence
	 * of pixels that have a value greater than the lowerThrehsold.
	 **/
	template <class T> class HysteresisThresholding: public Thresholding<T>
	{
	public:
		/**
		 * Constructor which you have to give the lower and upper limits.
		 *
		 * @param lowerThreshold The threshold used when the edge has been found.
		 * @param upperThreshold The threshold used when locating the edges beginning.
		 **/
		HysteresisThresholding(T lowerThreshold, T upperThreshold): Thresholding<T>(lowerThreshold,upperThreshold){}
		/**
		 * Constructor which calculates the threshold from the histogram given in constructor.
		 * The limits are calculated so that the higher limit is 80 % of all the pixels and the
		 * lower limit is first nonzero index (in histogram) + higher limit divided by 2.
		 *
		 * @param histogram The integer histogram (featureVector) of the Matrix.
		 **/
		HysteresisThresholding(const util::List<int>& histogram):Thresholding<T>()
		{estimateThreshold(histogram);}
		/**
		 * Destructor of HysteresisThresholding.
		 **/
		~HysteresisThresholding(){}
		/**
		 * The fuction makes the hysteresis threshold for the Matrix.
		 * Note that if Constructor which calculates the threshold is used
		 * then the threshold are only specified for that specific Matrix and if
		 * other threshold calculation is wanted the estimateThreshold fuction
		 * must be used.
		 *
		 * @param mat The matrix wanted to threshold.
		 **/
		util::Matrix<T> getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&);
		/**
		 * Estimate Theshold estimates the threshold from the histogram of matrix.
		 * The limits are calculated so that the higher limit is 80 % of all the pixels and the
		 * lower limit is first nonzero index (in histogram) + higher limit divided by 2.
		 * The limit values will be saved on the Threshold object and they will be get
		 * with methods getLowerThreshold and getUpperThreshold. The old threshold values
		 * will be destroyed.
		 *
		 * @param histogram The integer histogram (featureVector) of the Matrix.
		 **/
		void estimateThreshold(const util::List<int>& histogram);

	private:
		/**
		 * The traceEdge traces the edge recursively as long as it find edge which is
		 * Bigger than lowerThreshold.
		 *
		 * @param mat The matrix wanted to process.
		 * @param result The result matrix of threshold.
		 * @param i The index of row processed.
		 * @param j The index of column processed.
		 **/
		bool traceEdge(const util::Matrix<T>& mat, util::Matrix<T>& result, int i,int j);
	};

	template <class T> util::Matrix<T> HysteresisThresholding<T>::getTransformedImage(const util::Matrix<T>& mat)
		throw (ImageTransformException&)
	{
		int rows = mat.getRows();
		int cols = mat.getColumns();
		if(rows < 1 || cols <1)
			throw ImageTransformException("HysteresisThreshold<T>::getTransformedImage(const util::Matrix<T>&) : Cannot make threshold for an empty matrix.");
		// make the matrix
		util::Matrix<T> result(rows,cols);
		
		// For each edge with a magnitude above the high threshold, begin
		// tracing edge pixels that are above the low threshold.                
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				if (mat(i,j) >= _upperThreshold)traceEdge(mat,result,i,j);
	
		return result;
	}
	
	template <class T> bool HysteresisThresholding<T>::traceEdge(const util::Matrix<T>& mat, util::Matrix<T>& result, int i, int j)
	{
		int rows = mat.getRows();
		int cols = mat.getColumns();
		bool flag = false;
		if (result(i,j) == 0)
			{
				result(i,j)=mat(i,j);
				flag=false;
				for (int r= -1; r<=1; r++)
					{
						for(int c= -1; c<=1; c++)
							if ((i+r > 0 && i+r < rows)&&(j+c>0 && j+c<cols) && mat(i+r,j+c) >= _lowerThreshold)
								{
									flag=traceEdge(mat,result,i+r, j+c);
									if (flag)break;
								}//if
						if (flag) break;
					}
				return true;
			}
		return false;
	}
	
	template <class T> void HysteresisThresholding<T>::estimateThreshold(const util::List<int>& histogram)
	{
		int len=histogram.getLength();
		int sum = util::Math::sum(histogram); 
		// The high threshold should be > 80 or 90% of the pixels
		int limit = (int)(0.2*double(sum));
		int maxLimit = len-1;
		int count = histogram[maxLimit];
		// find max limit while 
		while (count < limit && !(--maxLimit))count += histogram[maxLimit];
		_upperThreshold=T(maxLimit);
		//then find the lower threshold
		int i=0;
		while (histogram[i]==0) i++;
		_lowerThreshold = T(double(maxLimit+i)/2.0);
	}
	
}}
#endif