/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001-2003 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.19 $
 *********************************************************************/

#ifndef _CIRCULARLOCALSAMPLER_H
#define _CIRCULARLOCALSAMPLER_H

#include "../graphics/Point.h"
#include "../ImageTransform.h"
#include <math.h>

#define interpolateAtPtr(u,i,c) ((_dpMultipliers[(i<<2)] != 1.0) ? \
(*(u)*_dpMultipliers[(i<<2)] + *(u+1)*_dpMultipliers[(i<<2)+1] + *(u+c)*_dpMultipliers[(i<<2)+2] + *(u+c+1)*_dpMultipliers[(i<<2)+3] + 1e-10) : \
*(u))

#define intRound(u,i,c) int(interpolateAtPtr(u,i,c)+0.5)
#define round(x) int(x + ((x>=0)? 0.5 : -0.5))

namespace prapi { namespace texture {

	/**
	 * CircularLocalSampler is a common superclass for all image
	 * transforms that work with a circular neighborhood. "Circular
	 * neighborhood" is used to denote a situation where, instead of the
	 * traditional rectangular one, neighborhood pixels are defined to
	 * be the ones that lie at a certain distance from the center. The
	 * distance is also called "predicate". The number of samples at
	 * this distance and the predicate itself can be dynamically
	 * changed. In digital images, all pixels in a circular neighborhood
	 * do not necessarily match the pixel grid. Pixel values at these
	 * positions are obtained with bilinear interpolation or, if the
	 * interpolation flag is set to false, from the pixel nearest to the
	 * exact position.
	 **/
	template <class T, class U> class CircularLocalSampler : public ImageTransform<T,U>
	{
	public:
		/**
		 * Create a new CircularLocalSampler instance. Make sure the
		 * compiler understands you are giving the predicate as a floating
		 * point value. Otherwise, it may have difficulties in detecting
		 * whether this or the other constructor should be used. Example:
		 * <pre>
		 * new CircularLocalSampler<int,int>(8,1.0);
		 * </pre>
		 *
		 * @param predicate the sampling radius
		 * @param samples the number of samples in the circular neighborhood
		 **/
		CircularLocalSampler(unsigned int samples, double predicate);
		/**
		 * Create a new CircularLocalSampler instance. Make sure the
		 * compiler understands you are giving the predicate as an
		 * unsigned integer Otherwise, it may have difficulties in
		 * detecting whether this or the other constructor should be used. 
		 * Example:
		 * <pre>
		 * new CircularLocalSampler<int,int>(8,1u);
		 * </pre>
		 *
		 * @param samples the number of samples in the circular neighborhood
		 * @param predicate the sampling radius
		 * @param interpolated if true, bilinear interpolation is in use.
		 * Otherwise, nearest neighbor is used.
		 **/
		CircularLocalSampler(unsigned int samples = 8, unsigned int predicate = 1, bool interpolated = true);

		virtual ~CircularLocalSampler() { delete[] _pPoints; delete[] _pOffsets; delete[] _dpMultipliers; }
		/**
		 * Set the predicate value. This will update the precalculated
		 * point table. Example:
		 * <pre>
		 * sampler.setPredicate(3u);
		 * </pre>
		 *
		 * @param p the new predicate. Must be larger than 0.
		 **/
		void setPredicate(double p) { _dPredicate = p; _uiPredicate = (unsigned int)ceil(p); updatePoints(); }
		/**
		 * Set the predicate value, i.e. the radius of the neighborhood. 
		 * This will update the precalculated point table. Example:
		 * <pre>
		 * sampler.setPredicate(3u);
		 * </pre>
		 *
		 * @param p the new predicate. Must be larger than 0.
		 **/
		void setPredicate(unsigned int p) { _dPredicate = p; _uiPredicate = p; updatePoints(); }
		/**
		 * Set the predicate value. This will update the precalculated
		 * point table.
		 *
		 * @param p the new predicate. Must be larger than 0.
		 **/
		void setPredicate(int p) { _dPredicate = p; _uiPredicate = p; updatePoints(); }
		/**
		 * Get the radius of the circular neighborhood.
		 **/
		double getPredicate(void) { return _dPredicate; }

		/**
		 * Set the interpolation flag.
		 * @param intp if true, points that do not match pixels exactly
		 *        are interpolated. If false, nearest pixels are used.
		 **/
		void setInterpolated(bool intp) { _bInterpolated = intp; }
		/**
		 * Get the interpolation flag.
		 **/
		bool isInterpolated(void) { return _bInterpolated; }

		/**
		 * Set the number of samples in the circular neighborhood. The
		 * precalculated point table will be updated.
		 **/
		void setSamples(unsigned int samples);
		/**
		 * Get the number of samples in the circular neighborhood.
		 **/
		unsigned int getSamples() { return _uiSamples; }

		/**
		 * Set the parameters of the circular neighborhood.
		 *
		 * @param samples the number of neighborhood samples (e.g. 8u)
		 * @param the radius of the neighborhood (e.g. 1.5)
		 **/
		void setNeighborhood(unsigned int samples, double predicate)
		{
			_dPredicate = predicate;
			_uiPredicate = (unsigned int)ceil(predicate);
			setSamples(samples);
		}

		/**
		 * Set the parameters of the circular neighborhood.
		 *
		 * @param samples the number of neighborhood samples (e.g. 8u)
		 * @param the radius of the neighborhood (e.g. 2u)
		 **/
		void setNeighborhood(unsigned int samples, unsigned int predicate)
		{
			_dPredicate = predicate;
			_uiPredicate = predicate;
			setSamples(samples);
		}

		/**
		 * Subclasses must implement this method to return a transformed
		 * value for a center pixel given a known circular neighborhood. 
		 * The default implementation returns 0. This method is called for
		 * each pixel in an matrix from getTransformedImage(Matrix<U>&). 
		 * For example, the local variance operator returns the variance
		 * of the values in <i>surrounding</i>. LBP, in turn, compares
		 * each value in <i>surrounding</i> to the value of the center
		 * pixel, and builds up a binary code that it then returns.
		 *
		 * @param center the value of the current center pixel
		 * @param surrounding the values of the surrounding pixels starting
		 *        from horizontal x-axis direction (3 o'clock) and proceeding
		 *        clockwise. The length of this list is equal to the number
		 *        of samples set with the setSamples() method.
		 * @return a transformed value for the center pixel.
		 **/
		virtual T getValue(U center, util::List<U>& surrounding) { return 0; }

		/**
		 * Goes through each matrix pixel, calculates the neighborhood
		 * values using the defined sample count and predicate, either
		 * with or without interpolation. For each neighborhood, calls
		 * getValue(U, List<U>&) and stores the returned value into the
		 * transformed matrix. The size of the resulting matrix is smaller
		 * than the input due to border effects. If you set the radius of
		 * the neighborhood to 3, the width and height of the result
		 * matrix will be 4 (=(3-1)*2) pixels smaller than those of the
		 * input.
		 *
		 * @param mat the input matrix
		 * @return the transformed matrix
		 **/
		virtual util::Matrix<T> getTransformedImage(const util::Matrix<U>& mat) throw (ImageTransformException&);

		/**
		 * Get the angle at which the first neighbor is located. Angle is
		 * represented as radians in counter-clockwise direction. The
		 * direction of the positive x axis is 0. This value defaults to
		 * 0.
		 **/
		double getStartAngle() const { return _dStartAngle; }
		/**
		 * Set the angle of the first neighbor point.
		 **/
		void setStartAngle(double angle) { _dStartAngle = angle; updatePoints(); }

	protected:
		/**
		 * The current predicate (radius), i.e. the distance of the
		 * neighborhood from its center.
		 **/
		double _dPredicate;
		/**
		 * The predicate rounded to the closest integer larger than
		 * _dPredicate.
		 **/
		unsigned int _uiPredicate;
		/**
		 * The interpolation flag.
		 **/
		bool _bInterpolated;
		/**
		 * The number of samples in the local neighborhood.
		 **/
		unsigned int _uiSamples;
		/**
		 * The angle of the first neighbor.
		 **/
		double _dStartAngle;
		/**
		 * A precalculated table of interpolation points.
		 **/
		prapi::graphics::Point<int> *_pPoints;
		/**
		 * A precalculated table of interpolation offsets.
		 **/
		prapi::graphics::Point<double> *_pOffsets;
		/**
		 * Precalculated values for interpolation multiplication.
		 **/
		double *_dpMultipliers;

		/**
		 * Update the precalculated point value tables.
		 **/
		void updatePoints(void);
	};
	

	template <class T, class U> CircularLocalSampler<T,U>::CircularLocalSampler(unsigned int samples, double p) :
		_dPredicate(p), _uiPredicate((unsigned int)ceil(p)), _bInterpolated(true), _uiSamples(samples), _dStartAngle(0.0),
		_pPoints(NULL), _pOffsets(NULL), _dpMultipliers(NULL)
	{
		setSamples(samples);
	}

	template <class T, class U> CircularLocalSampler<T,U>::CircularLocalSampler(unsigned int samples,
																																							unsigned int p,
																																							bool intp) :
		_dPredicate(p), _uiPredicate(p), _bInterpolated(intp), _uiSamples(samples), _dStartAngle(0.0),
		_pPoints(NULL), _pOffsets(NULL), _dpMultipliers(NULL)
	{
		setSamples(samples);
	}

	template <class T, class U> void CircularLocalSampler<T,U>::setSamples(unsigned int samples)
	{
		delete[] _pPoints; delete[] _pOffsets; delete[] _dpMultipliers;
		_pPoints = new prapi::graphics::Point<int>[samples];
		_pOffsets = new prapi::graphics::Point<double>[samples];
		_dpMultipliers = new double[samples*4];
		_uiSamples = samples;
		updatePoints();
	}

	template <class T, class U> void CircularLocalSampler<T,U>::updatePoints(void)
	{
		double step = 2 * M_PI / _uiSamples;
		for (unsigned int i=0;i<_uiSamples;i++)
			{
				double tmpX = _dPredicate * cos(i * step + _dStartAngle);
				double tmpY = _dPredicate * sin(i * step + _dStartAngle);
				_pPoints[i].x = (int)tmpX;
				_pPoints[i].y = (int)tmpY;
				_pOffsets[i].x = tmpX - _pPoints[i].x;
				_pOffsets[i].y = tmpY - _pPoints[i].y;
				if (_pOffsets[i].x < 1.0e-10 && _pOffsets[i].x > -1.0e-10) //rounding error
					_pOffsets[i].x = 0;
				if (_pOffsets[i].y < 1.0e-10 && _pOffsets[i].y > -1.0e-10) //rounding error
					_pOffsets[i].y = 0;
				
				if (tmpX < 0 && _pOffsets[i].x != 0)
					{
						_pPoints[i].x -= 1;
						_pOffsets[i].x += 1;
					}
				if (tmpY < 0 && _pOffsets[i].y != 0)
					{
						_pPoints[i].y -= 1;
						_pOffsets[i].y += 1;
					}
				//cerr << i << ": (" << _pPoints[i].x << ", " << _pOffsets[i].x << ") (" << _pPoints[i].y << ", " << _pOffsets[i].y << ")" << endl;

				double dx = 1-_pOffsets[i].x;
				double dy = 1-_pOffsets[i].y;

				_dpMultipliers[i*4+0] = dx*dy;
				_dpMultipliers[i*4+1] = _pOffsets[i].x*dy;
				_dpMultipliers[i*4+2] = dx*_pOffsets[i].y;
				_dpMultipliers[i*4+3] = _pOffsets[i].x*_pOffsets[i].y;
				//cerr << "   " << multipliers[i][0] << ", " << multipliers[i][1] << ", " << multipliers[i][2] << ", " << multipliers[i][3] << endl;
			}
	}

	//  	template <class T, class U>
	//  	template <class V> inline V interpolateAtPtr(V* upperLeft, int i, int columns)
	//  	{
	//  		return T(*upperLeft*multipliers[i][0] +
	//  						 *(upperLeft+1)*multipliers[i][1] +
	//  						 *(upperLeft+columns)*multipliers[i][2] +
	//  						 *(upperLeft+columns+1)*multipliers[i][3]);
	//  	}


	template <class T, class U> util::Matrix<T> CircularLocalSampler<T,U>::getTransformedImage(const util::Matrix<U>& mat)
		throw (ImageTransformException&)
	{
		util::List<U> values(_uiSamples);
		values.setLength(_uiSamples);
		
		int rows = mat.getRows(), columns = mat.getColumns();
		const U* center = mat.getData() + _uiPredicate * (1 + columns);
		const U* ptr[_uiSamples];
		
		int pred2 = _uiPredicate << 1;
		util::Matrix<T> result(rows-pred2,columns-pred2);
		T *resultPtr = result.getData();

		if (!_bInterpolated)
			{
				for (unsigned int i=0;i<_uiSamples;i++)
					{
						int xc = round(_pPoints[i].x + _pOffsets[i].x);
						int yc = round(_pPoints[i].y + _pOffsets[i].y);
						ptr[i] = center + xc + yc * columns;
					}

				for (int r=0;r<rows-pred2;r++)
					{
						for (int c=0;c<columns-pred2;c++)
							{
								for (unsigned int i=0;i<_uiSamples;i++)
									{
										values[i] = *ptr[i];
										ptr[i]++;
									}

								*resultPtr = getValue(*center, values);

								center++;
								resultPtr++;
							}
						for (unsigned int i=0;i<_uiSamples;i++)
							ptr[i] += pred2;
						center += pred2;
					}
			}
		else
			{
				for (unsigned int i=0;i<_uiSamples;i++)
					ptr[i] = center + _pPoints[i].x + _pPoints[i].y * columns;
				
				for (int r=0;r<rows-pred2;r++)
					{
						for (int c=0;c<columns-pred2;c++)
							{
								for (unsigned int i=0;i<_uiSamples;i++)
									{
										values[i] = U(interpolateAtPtr(ptr[i],i,columns));
										ptr[i]++;
									}
									
								*resultPtr = getValue(*center, values);

								center++;
								resultPtr++;
							}
						for (unsigned int i=0;i<_uiSamples;i++)
							ptr[i] += pred2;
						center += pred2;
					}
			}
		return result;
	}
}}

#endif