/*********************************************************************
 * 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 "Hough.h"
#include "../Gaussian.h"

namespace prapi { namespace transforms {
	
	template <> util::Matrix<double> CircularHough<double, graphics::Point<double> >::getTransformedImage(const util::Matrix<graphics::Point<double> >& 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<Point>&): End radius must be larger that start radius.");
		int radii = int(double(radEnd - _iStartRadius) / _dRadiusStep) + 1;

		//Gaussian weighting functions
		util::List<util::List<double> > gaussians(radii);
		//Arc lengths
		util::List<int> arcLengths(radii);

		if (_dGradientError > 0)
			{
				for (int i=0; i<radii; i++)
					{
						int arcLength = int(_dGradientError/M_PI*8*(_iStartRadius+i*_dRadiusStep)+0.5);
						util::List<double> gaussian(Gaussian::create1DGaussian(3,arcLength));
						gaussians += gaussian / util::Math::sum(gaussian);
						arcLengths += arcLength;
					}
			}
		
		util::Matrix<double> result(rows * radii, cols);
		graphics::MemoryGraphics<double> transform(result);
		graphics::Rectangle<int> rect(0, 0, cols, rows);
		transform.setClip(&rect);

		int r,c;
		AllItems(r,c,mat)
			{
				graphics::Point<double> point(mat(r,c));
				if (point.y < _dThreshold)
					continue;
				//If an estimate of the gradient error is given, draw
				//Gaussian-weighted arcs on the transformation domain
				if (_dGradientError > 0)
					{
						for (int i=0; i<radii; i++)
							{
								//Move the origin of the graphics context
								int radius = _iStartRadius + int(i*_dRadiusStep);
								double step = M_PI/(4*radius);
								double angle = point.x - _dGradientError;
								for (int j=0; j<arcLengths[i]; j++, angle+=step)
									{
										double si = sin(angle);
										double co = cos(angle);
										transform.setColor(absolute(point.y) * gaussians[i][j]);
										transform.modifyPixel<std::plus<double> >(c+radius*co, r+radius*si);
										transform.modifyPixel<std::plus<double> >(c-radius*co, r-radius*si);
									}
								transform.translate(0, rows);
							}
					}
				//If an esimate is not given, each edge point in the input
				//adds two points to the transform. (In the gradient direction
				//and its opposite.)
				else
					{
						double si = sin(point.x);
						double co = cos(point.x);
						transform.setColor(absolute(point.y));
						for (int i=0; i<radii; i++)
							{
								double radius = _iStartRadius + double(i)*_dRadiusStep;
								transform.modifyPixel<std::plus<double> >(c+radius*co, r+radius*si);
								transform.modifyPixel<std::plus<double> >(c-radius*co, r-radius*si);
								//transform.modifyPixel<std::plus<double> >(c, r, absolute(point.y));
								//if (!(c%4) && !(r%4))
								//	transform.drawLine(double(c), double(r),radius*co + c,radius*si + r,100);
								transform.translate(0, rows);
							}
					}
				//Shift drawing window back to origin
				transform.translate(0, -radii*rows);
			}
		return result;
	}
}}