/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2002 Topi Mäenpää
 * 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.1 $
 *********************************************************************/

#ifndef _PEAKDETECTOR_H
#define _PEAKDETECTOR_H

#include <Math.h>
#include <MatrixUtils.h>
#include <SortedList.h>
#include <Pair.h>
#include "binary/Component.h"
#include "binary/Labeling.h"
#include "binary/BinaryMorphology.h"
#include "graphics/Point.h"
#include "ConvolutionMask.h"
#include "Kernel.h"
#include "Gaussian.h"
#include "Heap.h"

#include <functional>

namespace prapi
{
	/**
	 * GlobalPeakDetector is a class for detecting peaks in images, or
	 * in transform domains, like the Hough transform. The detection
	 * method is as follows. First, the input image is globally
	 * thresholded to produce a binary image (hence the name). Second, a
	 * binary opening operation is performed to reduce noise. Third,
	 * connected components are found. Finally, the found components are
	 * sorted according to their magnitude (the magnitude of the input
	 * image at the center of the mass of the component) and returned. 
	 * The user can decide how many peaks are returned at most.
	 **/
	template <class T> class GlobalPeakDetector : virtual public util::Object
	{
	public:
		/**
		 * Create a PeakDetector that binarizes input images using the
		 * given threshold. The maximum number of peaks returned is given
		 * by the maxPeaks parameter, which defaults to ten.
		 *
		 * @param threshold global binarization threshold
		 * @param maxPeaks the maximum number of peaks to find. 
		 * Non-positive number means all.
		 **/
		GlobalPeakDetector(T threshold=0, int maxPeaks=10) :
			_threshold(threshold), _iMaxPeaks(maxPeaks), _bAutoThreshold(false), _uiOpenRadius(1) {}

		/**
		 * Set the radius of the circular structuring element used to
		 * reduce noise. In noise-free cases, this value can be set to
		 * zero, disabling the morphological opening. The size of the
		 * structuring element will be (radius*2+1)-by-(radius*2+1). The
		 * default value is one.
		 **/
		void setOpenMaskRadius(unsigned int radius) { _uiOpenRadius = radius; }

		/**
		 * Get the radius of the structuring element used to reduce noise.
		 **/
		unsigned int getOpenMaskRadius() const { return _uiOpenRadius; }
		
		/**
		 * Set automatic thresholding on. If automatic thresholding is on,
		 * the input image is binarized with the mean of the image plus
		 * two times standard deviation.
		 **/
		void setAutoThreshold(bool autoThreshold) { _bAutoThreshold = autoThreshold; }
		/**
		 * See if the automatic thresholding is on.
		 **/
		void getAutoThreshold() const { return _bAutoThreshold; }
		
		/**
		 * Set the threshold.
		 **/
		void setThreshold(T threshold) { _threshold = threshold; }

		/**
		 * Get the binarization threshold.
		 **/
		T getThreshold() const { return _threshold; }

		/**
		 * Set the maximum number of peaks to return.
		 **/
		void setMaxPeaks(int peaks) { _iMaxPeaks = peaks; }

		/**
		 * Get the maximum number of peaks.
		 **/
		int getMaxPeaks() const { return _iMaxPeaks; }

		/**
		 * Detect peaks in an image or transform domain. The result is a
		 * list of Pair objects that store the magnitude of the peak, and
		 * the properties of the binary connected component that forms the
		 * peak.
		 **/
		util::List<util::Pair<T,binary::Component> > detectPeaks(const util::Matrix<T>& mat) const
			throw (ImageTransformException&);

	private:
		T _threshold;
		int _iMaxPeaks;
		bool _bAutoThreshold;
		unsigned int _uiOpenRadius;
	};


	/**
	 * As the name implies, LocalPeakDetector detects peaks local in
	 * nature. Instead of the global thresholding in GlobalPeakDetector,
	 * this class uses Gaussian fit to find peaks. The width of the
	 * Gaussian function used in fitting can be determined by the user.
	 **/
	class LocalPeakDetector : virtual public util::Object
	{
	public:
		/**
		 * Create a peak detector that detects peaks with the given
		 * "sharpness". The smaller the peak radius, the sharper peaks are
		 * assumed. The radius parameter directly affects the size of the
		 * Gaussian mask used for fitting. The Gaussian is created so that
		 * at least 95% of its mass is within the given radius.
		 *
		 * @param radius the width of the peaks
		 * @param threshold correlation threshold. All peaks with a local
		 * correlation measure less than this value are discarded. The
		 * maximum correlation is 1.0. Note that correlation measures tend
		 * to decrease with mask size, i.e. peak radius.
		 * @param maxPeaks the maximum number of peaks to find. 
		 * Non-positive number means all.
		 **/
		LocalPeakDetector(double radius, double threshold=0, int peaks=10);

		/**
		 * Local peaks are represented as Pair objects that store the
		 * local magnitude of the peak (how well the peak correlates with
		 * the Gaussian), and the x (column) and y (row) coordinates of
		 * the peak.
		 **/
		typedef util::Pair<double,graphics::Point<int> > Peak;

		/**
		 * Set the "sharpness" of the peaks.
		 **/
		void setRadius(double radius);

		/**
		 * Get the current peak radius.
		 **/
		double getRadius() const { return _dRadius; }
		
		/**
		 * Set the threshold.
		 **/
		void setThreshold(double threshold) { _dThreshold = threshold; }

		/**
		 * Get the binarization threshold.
		 **/
		double getThreshold() const { return _dThreshold; }

		/**
		 * Set the maximum number of peaks to return.
		 **/
		void setMaxPeaks(int peaks) { _iMaxPeaks = peaks; }

		/**
		 * Get the maximum number of peaks.
		 **/
		int getMaxPeaks() const { return _iMaxPeaks; }

		/**
		 * Detect peaks in an image or transform domain. The result is a
		 * list of 
		 **/
		template <class T> util::List<Peak> detectPeaks(const util::Matrix<T>& mat) const
			throw (ImageTransformException&);

	private:
		double _dThreshold;
		int _iMaxPeaks;
		double _dRadius;
		unsigned int _uiRadius;
		Matrix<double> _matMask;

		util::List<Peak> detectPeaks(const util::Matrix<double>&, int, int, int, int) const
			throw (ImageTransformException&);
	};


	template <class T> util::List<util::Pair<T,binary::Component> > GlobalPeakDetector<T>::detectPeaks(const util::Matrix<T>& mat) const
			throw (ImageTransformException&)
	{
		using namespace util;
		using namespace prapi::binary;

		T threshold(_threshold);
		if (_bAutoThreshold)
			{
				double mean,variance;
				Math::meanAndVariance(mat, mean, variance);
				threshold = T(mean + 2*sqrt(variance));
			}

		//Binarize
		Matrix<int> binary(MatrixUtils::compare<std::greater<T> >(mat, threshold));

		if (_uiOpenRadius)
			{
				//Open to reduce noise
				ConvolutionMask<int> se(Kernel::createCircular<int>(_uiOpenRadius));
				BinaryMorphology<int> morph(se,BIN_OPEN);
				binary = morph.getTransformedImage(binary);
			}

		//Label connected components
		//Find component parameters
		Labeling<int> labeling(CONNECT_4);
		binary = labeling.getTransformedImage(binary);

		List<Pair<T,Component> > result(256);
		int* data = binary.getData();
		int r,c;
		AllItems(r,c,binary)
			{
				if (*data)
					{
						if (*data > result.getLength())
							result.setLength(*data, Pair<T,Component>(0,Component()));
						int index = *data-1;
						result[index].second().comR += r;
						result[index].second().comC += c;
						result[index].second().mass++;
						if (mat(r,c) > result[index].first())
							result[index].first() = mat(r,c);
					}
				data++;
			}
		for (int i=result.getLength();i--;)
			{
				if (result[i].second().mass)
					{
						result[i].second().comR /= result[i].second().mass;
						result[i].second().comC /= result[i].second().mass;
					}
				result[i].second().labelIndex = i+1;
			}
		Heap<Pair<T,Component>, std::less<Pair<T,Component> > > sorted(result.getLength());
		sorted.addElements(result);
		sorted.sort();
		return sorted(0,_iMaxPeaks);
	}

	template <class T> util::List<LocalPeakDetector::Peak> LocalPeakDetector::detectPeaks(const util::Matrix<T>& mat) const
		throw (ImageTransformException&)
	{
		using namespace util;
		int maskRows = _matMask.getRows(), maskColumns = _matMask.getColumns();
		int rows = mat.getRows(), columns = mat.getColumns();
		return detectPeaks((Matrix<double>)MatrixUtils::extend(mat, _uiRadius), maskRows, maskColumns, rows, columns);
	}
}

#endif