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

#ifndef _EDGEDETECTOR_H
#define _EDGEDETECTOR_H

#include "../ImageTransform.h"
#include "../ConvolutionMask.h"
#include "../Image.h"
#include "../dsp/Convolution.h"
#include "../graphics/Point.h"
#include "Thresholding.h"

#include <MatrixUtils.h>
#include <Matrix.h>
#include <math.h>

namespace prapi { namespace texture {

	/**
	 * The EdgeDetectorException.
	 **/
	class EdgeDetectorException : public ImageTransformException 
	{
	public:
		EdgeDetectorException(std::string message) : ImageTransformException(message) {}
	};

	/**
	 * An enumeration for the getGradient function.
	 * <ul>
	 * <li>GRADIENT_X - The gradient in x direction.</li>
	 * <li>GRADIENT_Y - The gradient in y direction.</li>
	 * </ul>
	 **/
	enum Gradient { GRADIENT_X, GRADIENT_Y };

	/**
	 * An enumeration for already implemented, well-known edge
	 * detection masks. These include:
	 * 
	 * <ul>
	 * <li>SOBEL - Sobel in a 3x3 neighborhood</li>
	 * <li>PREWITT - Prewitt in a 3x3 neighborhood</li>
	 * <li>ROBERTS - Roberts in a 2x2 neighborhood</li>
	 * </ul>
	 **/
	enum EdgeMask { SOBEL, PREWITT, ROBERTS };

	/**
	 * The EdgeThinner class is an image transform that makes thin
	 * edges. A non-maximum suppression algorithm is used that looks for
	 * edge pixels. Once it finds one, it tries to find a stronger one
	 * in the direction of the edge gradient at the pixel. If a larger
	 * gradient magnitude value is found within the given radius, the
	 * edge pixel is set to zero.
	 **/
	class EdgeThinner : public ImageTransform<double, graphics::Point<double> >
	{
	public:
		/**
		 * Create a new EdgeThinner that uses the given radius for
		 * non-maximum suppression.
		 **/
		EdgeThinner(unsigned int radius) : _uiRadius(radius) {}
		
		/**
		 * Suppress non-maximum gradien values. Gradient angles must be in
		 * the range [0,2*M_PI). Otherwise the behavior is undefined.
		 *
		 * @param gradients gradient angles and magnitudes
		 * @param radius the search radius
		 **/
		static util::Matrix<double> suppressNonMaxima(const util::Matrix<graphics::Point<double> >& gradients, unsigned int radius);

		util::Matrix<double> getTransformedImage(const util::Matrix<graphics::Point<double> >& gradients)
			throw (ImageTransformException&)
		{
			return suppressNonMaxima(gradients, _uiRadius);
		}

	private:
		unsigned int _uiRadius;
	};

	/**
	 * The EdgeUtils class contains static utility methods for
	 * converting between different edge representations.
	 **/
	class EdgeUtils
	{
	public:
		/**
		 * Given a matrix containing gradient magnitudes in x and y
		 * directions, return the gradient angle for each pixel. Angles
		 * are represented as positive floating point numbers in the range
		 * [0,2*M_PI).
		 *
		 * @param mat gradient magnitudes in x and y directions
		 **/
		static util::Matrix<double> getAngle(const util::Matrix<prapi::graphics::Point<double> >& mat);

		/**
		 * Given a matrix containing gradient magnitudes in x and y
		 * directions, return the gradient magnitude for each pixel. The
		 * magnitude of a gradient is equal to the length of the gradient
		 * vector.
		 *
		 * @param mat gradient magnitudes in x and y directions
		 **/
		static util::Matrix<double> getMagnitude(const util::Matrix<prapi::graphics::Point<double> >& mat);

		/**
		 * Get either of the gradient coordinate layers from a two-layer
		 * gradient image.
		 *
		 * @param mat The gradient image.
		 * @param gradient The gradient which you want to separate (X or Y).
		 **/
		static util::Matrix<double> getGradient(const util::Matrix<prapi::graphics::Point<double> >& mat, Gradient gradient);

		/**
		 * Set a gradient coordinate layer to a value. With this method it
		 * is possible to alter gradient coordinates separately.
		 *
		 * @param matGradient Set gradient in this gradient image
		 * @param mat The gradient values
		 * @param gradient The gradient direction you want to alter
		 **/
		static void setGradient(util::Matrix<prapi::graphics::Point<double> >& matGradient,
														const util::Matrix<double>& mat, Gradient gradient)
			throw (EdgeDetectorException&);
	};
	
	/**
	 * EdgeDectector includes various utility functions for edge detection.
	 **/
  template <class T, class U=bool> class EdgeDetector: public ImageTransform<U,T>
	{
	public:
		/**
		 * Constructor of EdgeDetector.
		 *
		 * @param borderAction for the border handling. 
		 **/
		EdgeDetector(BorderAction borderAction = BORDER_REFLECT): _BorderAction(borderAction) {}
		/**
		 * Destructor of EdgeDetector.
		 **/
		~EdgeDetector(){}
		
	protected:
		/**
		 * BorderAction is for the border handling.
		 **/
		BorderAction _BorderAction;
	};
	

	/**
	 * This is the base class for all EdgeDetectors which are made by using two, X and Y masks.
	 * The base operator are defined in the class and the base operators like Sobel, Prewitt and
	 * Roberts.
	 **/
	template <class T, class U=bool> class DifferentialEdgeDetector: public EdgeDetector<T,U>
	{
	public:
		/**
		 * Constructor which will be used when implemented mask are used.
		 *
		 * @param mask The implemented mask used in EdgeDetection.
		 * @param radius The radius used in suppresNonMaxima. Note that when
		 *               radius is zero suppression isn't made.
		 * @param borderAction How to handel the borders when needeed.
		 **/
		DifferentialEdgeDetector(EdgeMask mask, BorderAction borderAction = BORDER_REFLECT):
			EdgeDetector<T,U>(borderAction), _Mask(mask)
		{ setMask(mask);}

		/**
		 * Constructor which will be used when own mask are used.
		 *
		 * @param maskX The convolution mask which is used calculatin the X gradient.
		 * @param maskY The convolution mask which is used calculatin the Y gradient.
		 * @param radius The radius used in suppresNonMaxima. Note that when
		 *               radius is zero suppression isn't made.
		 * @param borderAction How to handel the borders when needeed.
		 **/
		DifferentialEdgeDetector(ConvolutionMask<double>& maskX, ConvolutionMask<double>& maskY, BorderAction borderAction = BORDER_REFLECT):
			EdgeDetector<T,U>(borderAction), _ConvolutionMaskX(maskX), _ConvolutionMaskY(maskY), _Mask(SOBEL) {}

		/**
		 * The Destructor of Differentian edge detector
		 **/
		~DifferentialEdgeDetector(){}

		/**
		 * Function set's the mask (enumeration Mask) for EdgeDetector.
		 *
		 * @param mask The Mask wanted to use.
		 **/ 
		void setMask(EdgeMask mask);

		/**
		 * Get Convolution mask gives the convolution mask.
		 *
		 * @param gradient The X or Y convolution mask.
		 **/
		ConvolutionMask<double> getConvolutionMask(Gradient gradient);
		
		/**
		 * Calculate the edge gradients in an input image, and return the
		 * gradients in x and y directions.
		 *
		 * @param mat the input image.
		 **/
		util::Matrix<prapi::graphics::Point<double> > getGradientMatrices(const util::Matrix<T>& mat);

		/**
		 * Calculate edge gradients and return them as a gradient
		 * magnitude-angle image. Each pixel in the returned image stores
		 * the gradient angle (Point::x) and the gradient magnitude
		 * (Point::y) of the corresponding pixel of the input image.
		 **/
		util::Matrix<prapi::graphics::Point<double> > getGradientAngleAndMagnitude(const util::Matrix<T>& mat);
		
		/**
		 * Detect edges in an image. The returned image contains the
		 * detected edges.
		 *
		 * @param mat input image
		 * @return edges, typecasted to the type U
		 **/
		util::Matrix<U> getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&);

	protected:
	/**
		 * The constructor for the child classes.
		 **/
		DifferentialEdgeDetector(BorderAction borderAction):
			EdgeDetector<T,U>(borderAction), _ConvolutionMaskX(), _ConvolutionMaskY(),_Mask(SOBEL) {}
		/**
		 * The Convolution mask for the X and Y gradient.
		 **/
		ConvolutionMask<double> _ConvolutionMaskX;
		ConvolutionMask<double> _ConvolutionMaskY;
		
	private:
		EdgeMask _Mask;
	};
	
	/**
	 * This Class makes the Canny edge detection. It uses the gaussian mean and
	 * derivative mask to make the result. First the convolution is made for
	 * the matrix with gaussian mask then the suppress non maxima is made if
	 * the supressRadius is not zero and then usally the Hysteresis threshold is
	 * used when using Canny operator therefore you should use Hysteresis threshold
	 * for the matrix which you get from the getTransformedImage.
	 **/
	template <class T, class U=bool> class Canny: public DifferentialEdgeDetector<T,U>
	{
	public:
		/**
		 * This constructor makes the Canny masks. 
		 *
		 * @param sigma The standard deviation at the gausian masks.
		 * @param supressRadius The radius used in suppresNonMaxima. Note that when
		 *               radius is zero suppression isn't made.
		 * @param borderAction How to handel the borders when needeed.
		 **/
		Canny(double sigma=2.0, int suppressRadius=0, BorderAction borderAction = BORDER_REFLECT);
		/**
		 * The destructor of Canny.
		 **/
		~Canny(){}
		
		/**
		 * Get the standard deviation.
		 **/
		double getSigma(){return _dSigma;}
		/**
		 * Set the standard deviation. When the sigma is set the gaussian masks are
		 * calculated again.
		 **/
		void setSigma(double sigma);
			
	private:
		/**
		 * This fuction makes the gaussian mask used in the Canny operator.
		 *
		 * @param sigma The standard deviation in the gaussian mask.
		 **/
		ConvolutionMask<double> makeCannyMask(double sigma);
		/**
		 * The standard deviation in the gaussian mask.
		 **/
		double _dSigma;
	};
			
	/**
	 * Laplace of gaussian detects edges in a local 3-by-3 neighborhood
	 * using a simple convolution mask:
	 * <pre>
	 *   0 -1  0
	 *  -1  4 -1
	 *   0 -1  0
	 * </pre>
	 **/
	template <class T, class U=bool> class LaplaceOfGaussian: public EdgeDetector<T,U>
	{	
	public:
		/**
		 * Constructor for the LaplaceOfGaussian operator.
		 *
		 * @param borderAction for the border handling. 
		 **/
		 LaplaceOfGaussian(BorderAction borderAction = BORDER_CROP):
			EdgeDetector<T,U>(borderAction){}

		/**
		 * The destructor of LaplaceOfGaussian.
		 */
		~LaplaceOfGaussian(){}
		
		/**
		 * Makes the LaplaceOfGaussian operation for the matrix given in parameter.
		 *
		 * @param mat an input image
		 * @return the new LOG matrix which indicates the edges.
		 *         The type of the matrix is specified by the template parameter U.
		 **/
		util::Matrix<U> getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&);

	};

	
	template<class T,class U> void DifferentialEdgeDetector<T,U>::setMask(EdgeMask mask)
	{
		if(mask == SOBEL)
			{
				double maskFactorsY[] = {-1,-2,-1,0,0,0,1,2,1}; // the sobel mask factors
				double maskFactorsX[] = {-1,0,1,-2,0,2,-1,0,1}; // the sobel mask factors
				ConvolutionMask<double> maskX(3,3,maskFactorsX); // the X mask
				ConvolutionMask<double> maskY(3,3,maskFactorsY); // the Y mask
				_ConvolutionMaskX = maskX;
				_ConvolutionMaskY = maskY;
			}
		else if(mask == PREWITT)
			{
				double maskFactorsY[] = {-1,-1,-1,0,0,0,1,1,1}; // the prewitt mask factors
				double maskFactorsX[] = {-1,0,1,-1,0,1,-1,0,1}; // the prewitt mask factors
				ConvolutionMask<double> maskX(3,3,maskFactorsX); // the X mask
				ConvolutionMask<double> maskY(3,3,maskFactorsY); // the Y mask
				_ConvolutionMaskX = maskX;
				_ConvolutionMaskY = maskY;
			}
		else if(mask == ROBERTS)
			{
				double maskFactorsY[] = {1,0,0,-1}; // the roberts mask factors
				double maskFactorsX[] = {0,1,-1,0}; // the roberts mask factors
				ConvolutionMask<double> maskX(2,2,maskFactorsX); // the X mask
				ConvolutionMask<double> maskY(2,2,maskFactorsY); // the Y mask
				_ConvolutionMaskX = maskX;
				_ConvolutionMaskY = maskY;
			}
	}

	template<class T,class U> ConvolutionMask<double> DifferentialEdgeDetector<T,U>::getConvolutionMask(Gradient gradient)
	{
		if(gradient == GRADIENT_X) return _ConvolutionMaskX;
		return _ConvolutionMaskY;
	}
	
	
	template <class T, class U> util::Matrix<prapi::graphics::Point<double> > DifferentialEdgeDetector<T,U>::getGradientMatrices(const util::Matrix<T>& mat)
	{
		int rows=mat.getRows();
		int cols=mat.getColumns();
		
		util::Matrix<T> tempX(Image::convolve(mat,_ConvolutionMaskX, _BorderAction));
		util::Matrix<T> tempY(Image::convolve(mat,_ConvolutionMaskY, _BorderAction));
		
		util::Matrix<prapi::graphics::Point<double> > result(rows,cols);
		T* matXData = tempX.getData();
		T* matYData = tempY.getData();
		prapi::graphics::Point<double>* resultData = result.getData();

		for(int i=rows*cols; i--; resultData++,matXData++,matYData++)
			{
				resultData->x = double(*matXData);
				resultData->y = double(*matYData);
			}
		
		return result;
	}

	template <class T, class U> util::Matrix<prapi::graphics::Point<double> > DifferentialEdgeDetector<T,U>::getGradientAngleAndMagnitude(const util::Matrix<T>& mat)
	{
		int rows=mat.getRows();
		int cols=mat.getColumns();
		
		util::Matrix<T> tempX(Image::convolve(mat,_ConvolutionMaskX, _BorderAction));
		util::Matrix<T> tempY(Image::convolve(mat,_ConvolutionMaskY, _BorderAction));
		
		util::Matrix<prapi::graphics::Point<double> > result(rows,cols);
		T* matXData = tempX.getData();
		T* matYData = tempY.getData();
		prapi::graphics::Point<double>* resultData = result.getData();

		for(int i=rows*cols; i--; resultData++,matXData++,matYData++)
			{
				double x = double(*matXData), y = double(*matYData);
				resultData->x = atan2(y,x);
				resultData->y = sqrt(x*x+y*y);
			}
		
		return result;
	}
	
	template <class T, class U> util::Matrix<U> DifferentialEdgeDetector<T,U>::getTransformedImage(const util::Matrix<T>& mat)
		throw (ImageTransformException&)
	{
		// first calculate the gradient matrices
		util::Matrix<prapi::graphics::Point<double> > gradients(getGradientMatrices(mat));

		return util::Matrix<U>(EdgeUtils::getMagnitude(gradients));
	}

	template <class T, class U> util::Matrix<U> LaplaceOfGaussian<T,U>::getTransformedImage(const util::Matrix<T>& mat) throw (ImageTransformException&)
	{
		double maskFactors[] = {0,-1,0,-1,4,-1,0,-1,0}; // the LOG mask factors
		ConvolutionMask<double> LOGMask(3,3,maskFactors); // the LOG mask

		// then before returning typycast all items to U
		return util::Matrix<U>(Image::convolve(mat,LOGMask, _BorderAction));
	}

	template <class T,class U> Canny<T,U>::Canny(double sigma, int suppressRadius, BorderAction borderAction):
		DifferentialEdgeDetector<T,U>(suppressRadius,borderAction),_dSigma(sigma)
	{
		_ConvolutionMaskX = makeCannyMask(_dSigma);
		ConvolutionMask<double> tmp(_ConvolutionMaskX.getColumns(),_ConvolutionMaskX.getRows(),_ConvolutionMaskX.getData());
		_ConvolutionMaskY = tmp;
	}

	template <class T,class U> void Canny<T,U>::setSigma(double sigma)
	{
		_dSigma = sigma;
		_ConvolutionMaskX = makeCannyMask(_dSigma);
		ConvolutionMask<double> tmp(_ConvolutionMaskX.getColumns(),_ConvolutionMaskX.getRows(),_ConvolutionMaskX.getData());
		_ConvolutionMaskY = tmp;
	}
	
	template <class T,class U> ConvolutionMask<double> Canny<T,U>::makeCannyMask(double sigma)
	{
		// first make the meanGaussian mask
		int maskRadius=30;
		double squareSigma = sigma*sigma;
		double square2Sigma = 2.0*squareSigma;
		double divider = 6.0*M_PI*squareSigma;
		double i=-double(maskRadius);
		int j=0;

		// find the right size.
		do
			{
				if(exp(-i*i/square2Sigma)>0.0001)break;
				else j++;
				i++;
			}while(i<=maskRadius);

		// check all values which are bigger than 0.0001
		util::List<double> mFactors;
		for(double d=-maskRadius+j;d<=maskRadius-j;d++)
			{		
				double j=d+0.5,k=d-0.5;
				double tmp = (exp(-d*d/square2Sigma)+
											exp(-j*j/square2Sigma)+
											exp(-k*k/square2Sigma)
											)/divider;
				mFactors += tmp;
			}

		int mLen = mFactors.getLength();
		// then make the derivate Gaussian mask
		util::List<double> dFactors;
		for(double j=-double(mLen>>1);j<=mLen>>1;j++)
			dFactors += -j/squareSigma * exp(-j*j/square2Sigma);
		
		// then make the convolution for the factors
		util::List<double> factors(dsp::Convolution::conv(mFactors, dFactors, true));
		// then the convolution mask
		ConvolutionMask<double> filter(1,factors.getLength(),factors.getData());
				
		return filter;
	}
}}
#endif