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

#ifndef _HOUGH_H
#define _HOUGH_H

#include "../ImageTransform.h"
#include "../graphics/Graphics.h"
#include "../graphics/Point.h"
#include <Matrix.h>
#include <Math.h>
#include <math.h>
#include <functional>


namespace prapi { namespace transforms {

	/**
	 * Linear Hough transform. Detects lines in images. Lines are
	 * paramerized by their slope (degrees, zero at positive x axis,
	 * growing couter-clockwise) and distance to the origin (pixels),
	 * resulting in a two-dimensional transformation domain. In the
	 * transformation domain, columns represent different slopes and
	 * rows different distances to the origin. The origin of the domain
	 * is at the center of the image. The quantization of the
	 * transformation domain is user-specifiable. The linear Hough
	 * transform can be though of as an image-to-image transform as it
	 * effectively produces an image as output.
	 **/
	template <class T, class U> class LinearHough : public ImageTransform<T,U>
	{
	public:
		/**
		 * Create an instance of the linear Hough transform. There are two
		 * parameters that need to be initialized. The threshold will be
		 * intialized to zero and the size of the transformation result
		 * will be intialized to 180x180.
		 *
		 * @param threshold the value with which the matrix will be thresholded.
		 * @param maxRadius the number of rows in the transformed matrix
		 *                  (quantization can be made for this parameter).
		 * @param maxAngle the number of columns in the transformed matrix
		 *                  (quantization can be made for this parameter).
		 **/
		LinearHough(U threshold = 0, int maxRadius = 180, int maxAngle = 180) :
			_threshold(threshold),_iRadius(maxRadius),_iAngle(maxAngle){}

		virtual ~LinearHough(){}

		/**
		 * Transform an image.
		 **/
		util::Matrix<T> getTransformedImage(const util::Matrix<U>& mat) throw (ImageTransformException&);

		/**
		 * Set the transformed matrix size. (Maximum radius in the
		 * parametric representation.)
		 **/
		void setMaxRadius(int radius){_iRadius=radius;}
		/**
		 * Get the maximum radius.
		 **/
		int getMaxRadius() const {return _iRadius;}
		/**
		 * Set the transformed matrix size. (Maximum angle in the
		 * parametric representation.)
		 **/
		void setMaxAngle(int angle){_iAngle=angle;}
		/**
		 * Get the maximum angle.
		 **/
		int getMaxAngle() const {return _iAngle;}

		/**
		 * Set a new threshold for Hough transform.
		 **/
		void setThreshold(U threshold){_threshold =  threshold;}
		/**
		 * Get the current threshold.
		 **/
		U getThreshold(){return _threshold;}
		
	private:
		/**
		 * The value of threshold.
		 **/
		U _threshold;
		/**
		 * The sizes of Transformed Matrix.
		 **/
		int _iRadius;
		int _iAngle;
		
	};

	
	template <class T, class U> util::Matrix<T> LinearHough<T,U>::getTransformedImage(const util::Matrix<U>& mat)
		throw (ImageTransformException&)
	{
		int rows = mat.getRows();
		int cols = mat.getColumns();
		
		// check the nececary values
		if(_iRadius <= 0 || _iAngle <= 0)
			throw ImageTransformException("LinearHough::getTransformedImage(const Matrix&): The size of the transformation must be greater than zero.");
		else if(rows == 0 || cols == 0)
			throw ImageTransformException("LinearHough::getTransformedImage(const Matrix&): The size of the input matrix must be greater than zero.");
			
		double centerX = cols/2.0; // the center point of X coordinates
		double centerY = rows/2.0; // the center point of Y coordinates
		// the conversion values degs to rads and scale the omega
		// the calculation is = (M_PI/180.0)*(180.0/_iAngle)=M_PI/_iAngle.
		double conversion = M_PI/double(_iAngle); // the conversion value degs to rads
			
		double max = sqrt(rows*rows+cols*cols);
		double scaleRadius = double((_iRadius)/max);

		util::Matrix<T> result(_iRadius,_iAngle);
		//graphics::MemoryGraphics<T> g(result);
		util::List<double> cosAngleValues(_iAngle);
		util::List<double> sinAngleValues(_iAngle); //for optimization

		// calculate first the values of angles
		for(int omega=0; omega<_iAngle; omega++)
			{
				sinAngleValues += sin(double(omega)*conversion-M_PI/2);
				cosAngleValues += cos(double(omega)*conversion-M_PI/2);
			}
		
		// go through the matrix
		const U* data = mat.getData();
		for(int y=0;y<rows;y++)
			for(int x=0;x<cols;x++,data++)
				{
					if(*data > _threshold)
						{
							double coordX = x - centerX;
							double coordY = centerY - y;
							for(int omega=0; omega<_iAngle; omega++)
								{
									int r = int(scaleRadius*(coordX*cosAngleValues[omega] + coordY*sinAngleValues[omega])+0.5);
									if (r >= 0 && r < _iRadius)
										result(r,omega) += T(*data);
								}
						}
				}
		return result;
	}

	/**
	 * Circular Hough transform. The circular Hough transform detects
	 * circles in images. Circles are parametrized with their center
	 * pixel (R,C) and radius (r), resulting in a three-dimensional
	 * transform domain. The transformation method however returns a
	 * two-dimensional matrix. The third dimension is "virtualized" by
	 * copying two-dimensional domains on top of each other as shown in
	 * the following example:
	 *
	 * <pre>
	 * Two-dimensional transform domain:
	 *
	 *     +------> C
	 *     |
	 *     |
	 *     |
	 *     V
	 *     R
	 *
	 * Three-dimensional domain:
	 *
	 *     +------> C
	 *     |
	 * r=1 |
	 *     |
	 *    ...
	 *     |
	 * r=N |
	 *     |
	 *     V
	 *     R
	 * </pre>
	 *
	 * The circular Hough transform takes two forms. In the default
	 * form, input images are represented as gradient magnitude (edge)
	 * images. The transformation cannot be guided by gradient
	 * direction, which results in inaccurate and somewhat slow
	 * operation. If the U template parameter is of type
	 * Point&lt;double&gt;, each pixel in the image is treated as a
	 * gradient angle-magnitude pair. This kind of images can be
	 * obtained for example from DifferentialEdgeDetector. Given the
	 * gradient information, the transform can be optimized quite a lot.
	 **/
	template <class T, class U> class CircularHough : public ImageTransform<T,U>
	{
	public:
		/**
		 * Create an instance of the circular Hough transform. The
		 * transform detects circles with a number of different radii,
		 * specified with the start/endRadius parameters and radiusStep. 
		 * The total number of 2-D transformation domains stacked on each
		 * other is (endRadius - startRadius)/radiusStep.
		 *
		 * @param magnitudeThreshold in transformation, use only pixels
		 * with gradient magnitudes not less than this value.
		 * @param startRadius find circles with radii not smaller than this value
		 * @param endRadius find circles with radii not larger than this
		 * value. If endRadius is set to a number less than 1, all circles
		 * up to the largest possible one are found.
		 * @param radiusStep create a new 2-D transformation domain for
		 * each this many pixels.
		 **/
		CircularHough(T magnitudeThreshold,
									int startRadius=1, int endRadius=-1, double radiusStep=1) :
			_threshold(magnitudeThreshold),
			_iStartRadius(startRadius), _iEndRadius(endRadius), _dRadiusStep(radiusStep) {}

		/**
		 * Set the threshold.
		 **/
		void setThreshold(T threshold) { _threshold = threshold; }
		/**
		 * Get the current threshold.
		 **/
		T getThreshold() const { return _threshold; }

		/**
		 * Set the radii that the transform looks for, and the sampling
		 * step.
		 **/
		void setRadii(int start, int end, double step)
		{
			_iStartRadius = start; _iEndRadius = end;
			_dRadiusStep = step;
		}

		/**
		 * Get the start radius.
		 **/
		int getStartRadius() { return _iStartRadius; }
		/**
		 * Get the end radius.
		 **/
		int getEndRadius() { return _iEndRadius; }
		/**
		 * Get the sampling step.
		 **/
		double getRadiusStep() { return _dRadiusStep; }
		
		/**
		 * Transform a gradient image.
		 **/
		util::Matrix<T> getTransformedImage(const util::Matrix<U>& mat) throw (ImageTransformException&);

	private:
		T _threshold;
		int _iStartRadius, _iEndRadius;
		double _dRadiusStep;
	};

	
	/**
	 * The circular Hough transform for gradient images. The algorithm
	 * used here is a modified version of the one typically found in
	 * computer vision textbooks. Instead of blindly trusting the given
	 * gradient directions, an error term can be given to indicate the
	 * trustworthness of the gradient estimates. Since gradients are
	 * typically measured with 3x3 masks, they cannot be very accurate. 
	 * When a measure of the angular error is involved, the
	 * transformation method assumes the real value of the gradient is
	 * normally distributed around the given angle estimate. The width
	 * of the distribution is dependent on the error estimate.
	 **/
	template <class T> class CircularHough<T, graphics::Point<double> > : public ImageTransform<T, graphics::Point<double> >
	{
	public:
		/**
		 * Create an instance of the circular Hough transform. The
		 * transform detects circles with a number of different radii,
		 * specified with the start/endRadius parameters and radiusStep. 
		 * The total number of 2-D transformation domains stacked on each
		 * other is (endRadius - startRadius)/radiusStep.
		 *
		 * @param magnitudeThreshold in transformation, use only pixels
		 * with gradient magnitudes not less than this value.
		 * @param startRadius find circles with radii not smaller than this value
		 * @param endRadius find circles with radii not larger than this
		 * value. If endRadius is set to a number less than 1, all circles
		 * up to the largest possible one are found.
		 * @param radiusStep create a new 2-D transformation domain for
		 * each this many pixels.
		 * @param gradientError an estimate of the error in your gradient
		 * measurements (in radians). This value defines the range the
		 * real gradients will "most likely" fall into, given a gradient
		 * estimate. The transform assumes the real gradient is somewhere
		 * in the direction of the estimate +- gradientError.
		 **/
		CircularHough(double magnitudeThreshold,
									int startRadius=1, int endRadius=-1, double radiusStep=1,
									double gradientError=0) :
			_dThreshold(magnitudeThreshold),
			_iStartRadius(startRadius), _iEndRadius(endRadius), _dRadiusStep(radiusStep),
			_dGradientError(gradientError) {}

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

		/**
		 * Set the radii that the transform looks for, and the sampling
		 * step.
		 **/
		void setRadii(int start, int end, double step)
		{
			_iStartRadius = start; _iEndRadius = end;
			_dRadiusStep = step;
		}

		/**
		 * Get the start radius.
		 **/
		int getStartRadius() { return _iStartRadius; }
		/**
		 * Get the end radius.
		 **/
		int getEndRadius() { return _iEndRadius; }
		/**
		 * Get the sampling step.
		 **/
		double getRadiusStep() { return _dRadiusStep; }

		/**
		 * Set the gradient error estimate.
		 **/
		void setGradientError(double error) { _dGradientError = error; }

		/**
		 * Get the gradient error estimate.
		 **/
		double getGradientError() { return _dGradientError; }
		
		/**
		 * Transform a gradient image.
		 **/
		util::Matrix<T> getTransformedImage(const util::Matrix<graphics::Point<double> >& mat) throw (ImageTransformException&);

	private:
		double _dThreshold;
		int _iStartRadius, _iEndRadius;
		double _dRadiusStep;
		double _dGradientError;
	};

	template <class T, class U> util::Matrix<T> CircularHough<T,U>::getTransformedImage(const util::Matrix<U>& mat)
		throw (ImageTransformException&)
	{
		int rows =  mat.getRows();
		int cols =  mat.getColumns();
		int radEnd = _iEndRadius > 0 ? _iEndRadius : maximum(rows, cols)/2;
		if (_iStartRadius > _iEndRadius)
			throw ImageTransformException("CircularHough::getTransformedImage(const Matrix&): End radius must be larger that start radius.");
		int radii = int(double(radEnd - _iStartRadius) / _dRadiusStep) + 1;

		util::Matrix<double> result(rows * radii, cols);
		graphics::MemoryGraphics<double> transform(result);
		graphics::Rectangle<int> rect(0, 0, cols, rows);
		transform.setClip(&rect);
		transform.setDrawMode(DRAWMODE_PLUS);

		int r,c;
		AllItems(r,c,mat)
			{
				T magnitude = mat(r,c);
				if (magnitude < _threshold)
					continue;

				transform.setColor(magnitude);
				for (int i=0; i<radii; i++)
					{
						transform.drawCircle((double)c, (double)r, _iStartRadius + _dRadiusStep*i);
						transform.translate(0, rows);
					}
				//Shift drawing window back to origin
				transform.translate(0, -radii*rows);
			}
		return result;
	}
	
}}
#endif