/*
	File:			QMCSampleGen.cc

	Function:		Program to calculate stratified sample positions
					for... well, all sorts of things these days.

	Author(s):		Andrew Willmott

	Copyright:		Copyright (c) 1999, Andrew Willmott

	Notes:			Integrated all sample rules into this file + unified
					Samples.h, Samples.cc.
					Rewrote to use Quasi-Monte Carlo Sampling, 1/99

*/

#include "Geometry.h"
#include <stdio.h>

// --- Quasi-Monte Carlo sampling ---------------------------------------------

// Mapping routines

Point SquareToSphere(Coord c)
// convert c to a unit length screw
// this is not a particularly good mapping -- distortion at the poles.
// need better!
{
    GCLReal 	phirad, w;
	Point		result;
	
	c = 2.0 * c - vl_1;
	phirad =  2.0 * c[1] * vl_pi;
	w = sqrt(1.0 - sqr(c[0]));
	result = Point(w * cos(phirad), w * sin(phirad), c[0]);

	return(result);
}

Coord SquareToTriangle(Coord c)
{
	GCLReal t = sqrt(c[1]);
	Coord	result;	

	result[0] = t * c[0];
	result[1] = 1 - t;

	return(result);
}

Coord SquareToCircle(Coord c)
// From Chiu & Shirley
{
	GCLReal		phi, r;
	Coord		a = 2.0 * c - vl_1;	// a is now on [-1,1] ^ 2 

	if (a[0] > -a[1])    	// region 1 or 2
	{
		if (a[0] > a[1])   	// region 1, also |a| > |b| 

		{
			r = a[0];
			phi = (vl_pi / 4.0) * (a[1] / a[0]);
		}
		else				// region 2, also |b| > |a| 
		{
			r = a[1];
			phi = (vl_pi / 4.0) * (2.0 - (a[0] / a[1]));
		}
	}
	else         			// region 3 or 4
	{
		if (a[0] < a[1])	// region 3, also |a| >= |b|, a != 0 
        { 
            r = -a[0];
            phi = (vl_pi / 4.0) * (4.0 + (a[1] / a[0]));
		}
		else          		// region 4, |b| >= |a|, but a==0 and b==0 could occur.
		{
            r = -a[1];
            if (a[1] != 0)
                phi = (vl_pi / 4.0) * (6.0 - (a[0] / a[1]));
            else
                phi = 0.0;
		}
	}

	return(r * Coord(cos(phi), sin(phi)));
}

// Sample generation routines 

GCLReal HaltonSequence(Int n, Int b)
// return term i of the base b Halton sequence
{
	GCLReal		p, ip, result;
	Int			k, a;
	
	result = 0;
	ip = 1.0 / b;         

	for (p = ip, k = n ; k; p *= ip, k = (Int) (k / b))
		if (a = k % b)
			result += a * p;

	return(result);
}

Void PrintCoords(Coord *coords, Int n, Char *name, Char *brackets)
{
	Int		i;

	cout << "GCLReal " << name << brackets << "[2] =" << endl << "{" << endl;

	for (i = 0; i < n; i++)
	{
		if ((i & 0x3) == 0)
			cout << "    ";
			
		cout << coords[i][0] << ", " << coords[i][1];

		if (i == n - 1)
			cout << endl;
		else
		{
			cout << ", ";
			if ((i & 0x3) == 3)
				cout << endl;
		}
	}
	
	cout << "};" << endl << endl;
}

Void PrintPoints(Point *points, Int n, Char *name, Char *brackets)
{
	Int		i;

	cout << "GCLReal " << name << brackets << "[3] =" << endl << "{" << endl;

	for (i = 0; i < n; i++)
	{
		if ((i & 0x1) == 0)
			cout << "    ";
			
		cout << points[i][0] << ", " << points[i][1] << ", " << points[i][2];

		if (i == n - 1)
			cout << endl;
		else
		{
			cout << ", ";
			if ((i & 0x1) == 1)
				cout << endl;
		}
	}
	
	cout << "};" << endl << endl;
}

Void PrintCoordsArray(Coord *coords, Int n, Char *name)
{
	static Char str[16];

	sprintf(str, "[%d]", n);
	PrintCoords(coords, n, name, str);
}

Void PrintPointsArray(Point *points, Int n, Char *name)
{
	static Char str[16];

	sprintf(str, "[%d]", n);
	PrintPoints(points, n, name, str);
}

Void GenSquareCoords(Int m, Char *name)
{
	Coord		*coords;
	Int			i;

	coords = new Coord[m];

	for (i = 0; i < m; i++)
		coords[i] = Coord(HaltonSequence(i + 1, 2), HaltonSequence(i + 1, 3));

	PrintCoordsArray(coords, m, name);

	delete[] coords;
}

Void GenCubePoints(Int m, Char *name)
{
	Point		*points;
	Int			i;
	
	points = new Point[m];

	for (i = 0; i < m; i++)
		points[i] = Point(HaltonSequence(i + 1, 2), HaltonSequence(i + 1, 3),
			HaltonSequence(i + 1, 5));

	PrintPointsArray(points, m, name);	

	delete[] points;
}

Void GenTriCoords(Int m, Char *name)
{
	Coord		*coords;
	Int			i;

	coords = new Coord[m];

	for (i = 0; i < m; i++)
		coords[i] = SquareToTriangle(Coord(HaltonSequence(i + 1, 2), HaltonSequence(i + 1, 3)));
	
	PrintCoordsArray(coords, m, name);	

	delete[] coords;
}

Void GenCircleCoords(Int m, Char *name)
{
	Coord		*coords;
	Int			i;

	coords = new Coord[m];

	for (i = 0; i < m; i++)
	{
		coords[i] = SquareToCircle(Coord(HaltonSequence(i + 1, 2), HaltonSequence(i + 1, 3)));
	}
	
	PrintCoordsArray(coords, m, name);	

	delete[] coords;
}

Void GenSpherePoints(Int m, Char *name)
{
	Point		*points;
	Int			i;
	
	points = new Point[m];

	for (i = 0; i < m; i++)
		points[i] = SquareToSphere(Coord(HaltonSequence(i + 1, 2), HaltonSequence(i + 1, 3)));

	PrintPointsArray(points, m, name);	

	delete[] points;
}


// --- main -------------------------------------------------------------------


static Char *gHeaderMsg[] = 
{
	"/*",
	"	File:           Samples.h",
	"",
	"	Function:       Gives samples over various unit entities.",
	"*/",
	"",
	"#include \"Geometry.h\"",
	"",
	0
};

Int main(Int argc, Char *argv[])
{
	GCLReal	jitterAmt = 0.8;
	Int		numArrays = 16;
	Int		numGenSamples = 256;
	
	if (argc != 2 || argv[1][1] != 'h')
	{
		cout << "// This file was automatically generated by " << argv[0] << endl
			 << endl;
			 
		cout << "#include \"Samples.h\"" << endl << endl;
		
		cerr << "gen square..." << endl;
		GenSquareCoords(numGenSamples, "tSquareSamples");
		cerr << "gen triangle..." << endl;
		GenTriCoords(numGenSamples, "tTriangleSamples");
		cerr << "gen circle..." << endl;
		GenCircleCoords(numGenSamples, "tCircleSamples");
		cerr << "gen cube..." << endl;
		GenCubePoints(numGenSamples, "tCubeSamples");
		cerr << "done." << endl;
		GenSpherePoints(numGenSamples, "tSphereSamples");
		cerr << "done." << endl;
	}
	else
	{
		Char	**p = gHeaderMsg;

		while (*p)
			printf("%s\n", *p++);

		printf("const Int kMaxStratSamples = %d;\n", numGenSamples);
		printf("\n");

		printf("extern GCLReal tSquareSamples[%d][2];\n", numGenSamples);
		printf("extern GCLReal tTriangleSamples[%d][2];\n", numGenSamples);
		printf("extern GCLReal tCircleSamples[%d][2];\n", numGenSamples);
		printf("extern GCLReal tCubeSamples[%d][3];\n", numGenSamples);
		printf("extern GCLReal tSphereSamples[%d][3];\n", numGenSamples);
	}		

	return(0);
}