/*********************************************************************
 * 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 $
 *********************************************************************/

#include "PeakDetector.h"

using namespace util;

namespace prapi
{
	LocalPeakDetector::LocalPeakDetector(double radius, double threshold, int peaks) :
		_dThreshold(threshold), _iMaxPeaks(peaks)
	{
		setRadius(radius);
	}

	void LocalPeakDetector::setRadius(double radius)
	{
		_uiRadius = (unsigned int)ceil(radius);
		_dRadius = radius;

		//Create a Gaussian mask in which at least 95% of the
		//total mass lies within 'radius' pixels. The
		//constant x = 1.38591 is selected so that erf(x) = 0.95. 
		//Since uiRadius is rounded up, it is guaranteed that the
		//95% limit is always met.
		_matMask = Gaussian::create2DNormal(_uiRadius, _uiRadius*2+1, _dRadius/1.38591);
	}

	util::List<LocalPeakDetector::Peak> LocalPeakDetector::detectPeaks(const util::Matrix<double>& extended,
																																		 int maskRows, int maskColumns, int rows, int columns) const
		throw (ImageTransformException&)
	{
		Matrix<double> neighborhood(maskRows, maskColumns, false);

		int size = _iMaxPeaks > 0 ? _iMaxPeaks : 256;
		SortedList<Peak> result(size,256);
		if (_iMaxPeaks > 0)
			result.setMaximumSize(_iMaxPeaks);

		Matrix<double> tmp(rows,columns);

		for (int r=0; r<rows; r++)
			for (int c=0; c<columns; c++)
				{
					neighborhood = extended(r,c,maskRows,maskColumns);
					double sum = neighborhood.sum(), correlation = 0.0;
					if (sum != 0)
						{
							neighborhood /= neighborhood.sum();
							neighborhood.dotProduct(_matMask);
							correlation = neighborhood.sum();
						}
					tmp(r,c) = correlation;
					if (correlation >= _dThreshold)
						result += Peak(correlation, graphics::Point<int>(c,r));
				}

		return result;
	}
}