/**
 * @file VCRA-SZCZERBACKI-splatting.cpp
 * @brief Point Based Rendering Tool.
 *                                                                    
 * Point Based Rendering Tool (using OpenGL and GLUT).                 \n
 *                                                                     \n
 * Using partially the splatting code implemented by Jussi Rasanen for \n
 * his Master Thesis "Surface Splatting: Theory, Extensions and        \n
 * Implementation", 2002, at Helsinki University of Technology         \n
 *                                                                     \n
 * Course:      Visao Computacional e Realidade Aumentada (PUC-RIO).   \n
 * Instructor:  Marcelo Gattass.                                       \n
 * Student:     Ricardo Szczerbacki.                                   \n
 * Version:     1.00.                                                  \n
 * 2006.2                                                  \n
 */

#include "splatting.h"

#include "analytic.inl"

extern Mesh    *lwo;
extern Color   *buffer;
extern ABuffer *abuffer;

extern std::vector<Splat>	splats;

void loadLwoFile(char *filename) 
{
  glutSetCursor(GLUT_CURSOR_WAIT);

  lwo = new Mesh( filename );
	lwo->setParametrization();
	lwo->subdivide( subdivision );
	splats.resize(lwo->getVertexCount());
	for(int i=0;i<splats.size();++i)
	{
		float maxl = 0.0f,l;
		int valence, left, right;
		bool boundary;
		const std::vector<Mesh::ventry> &vertexlist = lwo->get1Neighbourhood( i, valence, boundary, left, right );
		for(int j=0;j<vertexlist.size();++j)
		{
			l = (lwo->getVertex(vertexlist[j].m_vertex).m_position - lwo->getVertex(i).m_position).length();
			maxl = max2(l,maxl);
		}

		splats[i].m_position = lwo->getVertex(i).m_position;
		splats[i].m_normal = lwo->getVertex(i).m_normal;
		splats[i].m_size = maxl;
    splats[i].m_color.set(1.0f,0.85f,0.65f, 1.0f);
	}

  glutSetCursor(GLUT_CURSOR_INHERIT);
}

void poissonDiscDistribution( Vector2 *points, int count )
{
	SS_ASSERT(count > 0);

	float threshold = 1.0f;

	//repeat until we have all subpixels distributed
	int j=0;
	while(j < count)
	{
		//try to place a random sample in a unit circle so that it is not closer
		// to earlier samples than the threshold 
		int tries = 100;
		int k;
		Vector2 p;
		do
		{
			//make a random sample inside unit circle
			p = randInCircle();

			//test whether the sample is closer to points already placed than the threshold
			for(k=0;k<j;k++)
			{
				if( (p - points[k]).length() < threshold )
				{	//too close, break
					tries--;
					break;
				}
			}
		} while(k < j && tries > 0);

		if( tries == 0 )
		{
			//if did not succeed, lower the threshold and redistribute all points
			threshold *= 0.98f;	//the closer the number is to 1, the slower the convergence
			j = 0;
		}
		else
		{
			//if succeeded, add the sample
			points[j] = p;
			j++;
		}
	}
}

float cumulativeGaussian[gausssamples];

void initCumulativeGaussian()
{
	//init cumulative gaussian
	float x = 0.0f;
	float xstep = prefiltercutoff / (float)gausssamples;
	x += xstep/2.0f;	//calculate x at the center of each slice
	cumulativeGaussian[0] = exp(-0.5f*x*x)*xstep;
	for(int i=1;i<gausssamples;i++)
	{
		x += xstep;
		cumulativeGaussian[i] = cumulativeGaussian[i-1] + exp(-0.5f*x*x)*xstep;
	}
	//normalize
	for(int i=0;i<gausssamples;i++)
	{
		cumulativeGaussian[i] /= cumulativeGaussian[gausssamples-1];
	}
}

float inverseCumulativeGaussian( float r )
{
	int i;
	for(i=0;i<gausssamples-1;i++)
	{
		if( r > cumulativeGaussian[i] && r < cumulativeGaussian[i+1] )
			break;
	}
	return (float)i/(float)gausssamples;
}

void makeJitter(int subpixels)
{
	SS_ASSERT(subpixels > 0);

	if( prefilterJitter ) delete[] prefilterJitter;
	if( lensJitter) delete[] lensJitter;

	prefilterJitter = new Vector2[subpixels*jitterPatterns];
	lensJitter = new Vector2[subpixels*jitterPatterns];

	svPrint("poisson disk distribution with %d points.",subpixels);
	Timer t;
	t.sample();
	initCumulativeGaussian();
	for(int p=0;p<jitterPatterns;p++)
	{
		poissonDiscDistribution(&prefilterJitter[p*subpixels],subpixels);
		for(int j=0;j<subpixels;j++)
		{
			Vector2 jit = prefilterJitter[p*subpixels+j];
			float r = jit.length();
			if( r > 0.0f ) jit *= 1.0f/r;
			jit *= inverseCumulativeGaussian(r) * prefiltersize;
			prefilterJitter[p*subpixels+j] = jit;
		}

		poissonDiscDistribution(&lensJitter[p*subpixels],subpixels);
		for(int j=0;j<subpixels;j++)
		{
			lensJitter[p*subpixels+j] *= lensRadius;
		}
	}
	t.sample();
	svPrint(" took %f seconds\n",t.getDelta());
}