/*
 * randomdirs.cpp
 *
 *  Created on: Jul 4, 2013
 *      Author: gregor
 */

#include <iostream>
#include <vector>

#include <gtest/gtest.h>

#include "GLTB/exception.h"
#include "GLTB/geom/randomdir.h"

TEST(DiskSampling, testDiskSampling)
{
	size_t numRadialBins = 20;
	size_t numAngularBins = 20;
	size_t numSamples = 100000000;

	// Idea: compute a 2D histogram of radial and angular distributions.
	// Ratios between histogram values should equal ratios between areas covered by bins

	std::vector<std::vector<unsigned long long> > counts(numRadialBins);
	for(size_t i = 0; i < numRadialBins; i++) {
		counts[i].resize(numAngularBins);
		for(size_t j = 0; j < numAngularBins; j++) {
			counts[i][j] = 0;
		}
	}

	gltb::Random rng;

	for(size_t i = 0; i < numSamples; i++) {
        glm::vec2 sample = gltb::getRandomSampleDisk<glm::vec2>(rng);

		double r = glm::length(sample);
		double phi = atan2(sample.y, sample.x);
		if(phi < 0) {
			phi += 2 * M_PI;
		}

		size_t binR = std::min(numRadialBins - 1.0, r * numRadialBins);
		size_t binPhi = std::min(numAngularBins - 1.0, phi / (2 * M_PI) * numAngularBins);

		counts[binR][binPhi]++;
	}

	// compute histogram and related statistics
	for(size_t i = 0; i < numRadialBins; i++) {
		for(size_t j = 0; j < numAngularBins; j++) {
			double innerRadius = i * (1.0 / numRadialBins);
			double outerRadius = (i + 1) * (1.0 / numRadialBins);
			double binArea = (1.0 / numAngularBins) * M_PI * (outerRadius * outerRadius - innerRadius * innerRadius);

			double areaFraction = binArea / M_PI;
			double countFraction = counts[i][j]/(double)numSamples;

			double ratio = countFraction / areaFraction;
//			std::cout << areaFraction << "\t" << countFraction << "\t" << ratio << std::endl;

			ASSERT_GE(ratio, 0.95);
			ASSERT_LE(ratio, 1.05);
		}
	}
}

TEST(HemisphereSampling, testHemisphereSampling)
{
	size_t numThetaBins = 20;
	size_t numPhiBins = 20;
	size_t numSamples = 100000000;

	// Idea: compute a 2D histogram of radial and angular distributions.
	// Ratios between histogram values should equal ratios between areas covered by bins

	std::vector<std::vector<unsigned long long> > counts(numThetaBins);
	for(size_t i = 0; i < numThetaBins; i++) {
		counts[i].resize(numPhiBins);
		for(size_t j = 0; j < numPhiBins; j++) {
			counts[i][j] = 0;
		}
	}

	gltb::Random rng;

	for(size_t i = 0; i < numSamples; i++) {
        glm::vec3 sample = gltb::getRandomDirHemisphere<glm::vec3>(rng);

		double r = glm::length(glm::vec2(sample.x, sample.y));
		double phi = atan2(sample.y, sample.x);
		if(phi < 0) {
			phi += 2 * M_PI;
		}

		size_t binR = std::min(numThetaBins - 1.0, r * numThetaBins);
		size_t binPhi = std::min(numPhiBins - 1.0, phi / (2 * M_PI) * numPhiBins);

		counts[binR][binPhi]++;
	}

	// compute histogram and related statistics
	for(size_t i = 0; i < numThetaBins; i++) {
		for(size_t j = 0; j < numPhiBins; j++) {
			double innerRadius = i * (1.0 / numThetaBins);
			double outerRadius = (i + 1) * (1.0 / numThetaBins);
			double innerTheta = asin(innerRadius);
			double outerTheta = asin(outerRadius);
			double binArea = (1.0 / numPhiBins) * 2 * M_PI * (cos(innerTheta) - cos(outerTheta));

			double areaFraction = binArea / (2 * M_PI);
			double countFraction = counts[i][j]/(double)numSamples;

			double ratio = countFraction / areaFraction;
//			std::cout << areaFraction << "\t" << countFraction << "\t" << ratio << std::endl;

			ASSERT_GE(ratio, 0.95);
			ASSERT_LE(ratio, 1.05);
		}
	}
}

void testSphereSampling(size_t numThetaBins, size_t numPhiBins, size_t numSamples, bool rotated)
{
	// Idea: compute a 2D histogram of radial and angular distributions.
	// Ratios between histogram values should equal ratios between areas covered by bins
	// Note: this test is sort of stupid and mirrors the -z half-space to +z

	std::vector<std::vector<unsigned long long> > counts(numThetaBins);
	for(size_t i = 0; i < numThetaBins; i++) {
		counts[i].resize(numPhiBins);
		for(size_t j = 0; j < numPhiBins; j++) {
			counts[i][j] = 0;
		}
	}

	gltb::Random rng;

	for(size_t i = 0; i < numSamples; i++) {
        glm::vec3 sample = gltb::getRandomDir<glm::vec3>(rng);

		if(rotated) {
			// rotate by swapping coordinates - should not affect result
			float temp;
			temp = sample.x;
			sample.x = sample.z;
			sample.z = temp;
		}

		double r = glm::length(glm::vec2(sample.x, sample.y));
		double phi = atan2(sample.y, sample.x);
		if(phi < 0) {
			phi += 2 * M_PI;
		}

		size_t binR = std::min(numThetaBins - 1.0, r * numThetaBins);
		size_t binPhi = std::min(numPhiBins - 1.0, phi / (2 * M_PI) * numPhiBins);

		counts[binR][binPhi]++;
	}

	// compute histogram and related statistics
	for(size_t i = 0; i < numThetaBins; i++) {
		for(size_t j = 0; j < numPhiBins; j++) {
			double innerRadius = i * (1.0 / numThetaBins);
			double outerRadius = (i + 1) * (1.0 / numThetaBins);
			double innerTheta = asin(innerRadius);
			double outerTheta = asin(outerRadius);
			double binArea = (1.0 / numPhiBins) * 2 * M_PI * (cos(innerTheta) - cos(outerTheta));

			double areaFraction = binArea / (2 * M_PI);
			double countFraction = counts[i][j]/(double)numSamples;

			double ratio = countFraction / areaFraction;
//			std::cout << areaFraction << "\t" << countFraction << "\t" << ratio << std::endl;

			ASSERT_GE(ratio, 0.95);
			ASSERT_LE(ratio, 1.05);
		}
	}
}

TEST(SphereSampling, testSphereSamplingDirect)
{
	testSphereSampling(20, 20, 100000000, false);
}

TEST(SphereSampling, testSphereSamplingRotated)
{
	testSphereSampling(20, 20, 100000000, true);
}

TEST(HemisphereCosineSampling, testHemisphereCosineSampling)
{
	// unimplemented
	ASSERT_TRUE(false && "unimplemented");
}

TEST(Reorient, testFindOrthogonalVector)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 dir = gltb::getRandomDir<glm::vec3>(rng);

        glm::vec3 result = gltb::findOrthogonalVector(dir);

		ASSERT_LE(fabs(glm::dot(result, dir)), 1e-7);
	}
}

TEST(Reorient, testReorientCentralAxisNoOrientation)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 axis = gltb::getRandomDir<glm::vec3>(rng);

		// choose vector from hemisphere to check scalar product afterwards
        glm::vec3 sourceVector = gltb::getRandomDirHemisphere<glm::vec3>(rng);

        glm::vec3 result = gltb::reorientCentralAxis(sourceVector, axis);

		ASSERT_GE(glm::dot(result, axis), 0);
		ASSERT_LE(fabs(glm::dot(sourceVector, glm::vec3(0, 0, 1)) - glm::dot(result, axis)), 5e-7);
	}
}

TEST(Reorient, testReorientCentralAxisFixedOrientation)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 axis = gltb::getRandomDir<glm::vec3>(rng);
        glm::vec3 axis2 = gltb::getRandomDir<glm::vec3>(rng);

		axis2 = glm::normalize(axis2 - glm::dot(axis, axis2) * axis);

		// choose vector from hemisphere to check scalar product afterwards
        glm::vec3 sourceVector = gltb::getRandomDirHemisphere<glm::vec3>(rng);

        glm::vec3 result = gltb::reorientCentralAxis(sourceVector, axis, axis2);

		ASSERT_GE(glm::dot(result, axis), 0);

		ASSERT_LE(fabs(glm::dot(sourceVector, glm::vec3(0, 0, 1)) - glm::dot(result, axis)), 1e-5);
		ASSERT_LE(fabs(glm::dot(sourceVector, glm::vec3(0, 1, 0)) - glm::dot(result, axis2)), 5e-6);
	}
}

TEST(Reorient, testReorientCentralAxisNoOrientationConsistency)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 axis = gltb::getRandomDir<glm::vec3>(rng);

		// choose vector from hemisphere to check scalar product afterwards
        glm::vec3 sourceVector1 = gltb::getRandomDirHemisphere<glm::vec3>(rng);
        glm::vec3 sourceVector2 = gltb::getRandomDirHemisphere<glm::vec3>(rng);
        glm::vec3 sourceVector3 = gltb::getRandomDirHemisphere<glm::vec3>(rng);

        glm::vec3 result1 = gltb::reorientCentralAxis(sourceVector1, axis);
        glm::vec3 result2 = gltb::reorientCentralAxis(sourceVector2, axis);
        glm::vec3 result3 = gltb::reorientCentralAxis(sourceVector3, axis);

		ASSERT_LE(fabs(glm::dot(sourceVector1, glm::vec3(0,0,1)) - glm::dot(result1, axis)), 1e-6);
		ASSERT_LE(fabs(glm::dot(sourceVector2, glm::vec3(0,0,1)) - glm::dot(result2, axis)), 1e-6);
		ASSERT_LE(fabs(glm::dot(sourceVector3, glm::vec3(0,0,1)) - glm::dot(result3, axis)), 1e-6);

		ASSERT_LE(fabs(glm::dot(sourceVector1, sourceVector2) - glm::dot(result1, result2)), 1e-6);
		ASSERT_LE(fabs(glm::dot(sourceVector1, sourceVector3) - glm::dot(result1, result3)), 1e-6);
		ASSERT_LE(fabs(glm::dot(sourceVector2, sourceVector3) - glm::dot(result2, result3)), 1e-6);
	}
}

TEST(Reorient, testReorientCentralAxisNoOrientationReverse)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 axis = gltb::getRandomDir<glm::vec3>(rng);

		// choose vector from hemisphere to check scalar product afterwards
        glm::vec3 sourceVector = gltb::getRandomDirHemisphere<glm::vec3>(rng);

        glm::vec3 result = gltb::reorientCentralAxis(sourceVector, axis);

        glm::vec3 inverse = gltb::reorientCentralAxisReverse(result, axis);

		ASSERT_LE(fabs(sourceVector.x - inverse.x), 3e-7);
		ASSERT_LE(fabs(sourceVector.y - inverse.y), 3e-7);
		ASSERT_LE(fabs(sourceVector.z - inverse.z), 3e-7);
	}
}

TEST(Reorient, testReorientCentralAxisFixedOrientationReverse)
{
	gltb::Random rng;

	size_t numIterations = 10000;

	for(size_t i = 0; i < numIterations; i++) {
        glm::vec3 axis = gltb::getRandomDir<glm::vec3>(rng);
        glm::vec3 axis2 = gltb::getRandomDir<glm::vec3>(rng);

		axis2 = glm::normalize(axis2 - glm::dot(axis, axis2) * axis);

		// choose vector from hemisphere to check scalar product afterwards
        glm::vec3 sourceVector = gltb::getRandomDirHemisphere<glm::vec3>(rng);

        glm::vec3 result = gltb::reorientCentralAxis(sourceVector, axis, axis2);

        glm::vec3 inverse = gltb::reorientCentralAxisReverse(result, axis, axis2);

		ASSERT_LE(fabs(sourceVector.x - inverse.x), 3e-7);
		ASSERT_LE(fabs(sourceVector.y - inverse.y), 3e-7);
		ASSERT_LE(fabs(sourceVector.z - inverse.z), 3e-7);
	}
}

TEST(Reorient, testReorientCentralAxisFixedOrientationReverseRegression)
{
	glm::vec3 normal(0.523393452, -0.330000997, 0.785594463);
	glm::vec3 tangent(0.832214653, 0, -0.554453611);
	glm::vec3 direction(-0.630678296, -0.515583992, -0.580015659);

    glm::vec3 result = gltb::reorientCentralAxis(direction, normal, tangent);

    glm::vec3 inverse = gltb::reorientCentralAxisReverse(result, normal, tangent);

	ASSERT_LE(fabs(direction.x - inverse.x), 3e-7);
	ASSERT_LE(fabs(direction.y - inverse.y), 3e-7);
	ASSERT_LE(fabs(direction.z - inverse.z), 3e-7);
}


int main(int argc, char *argv[])
{
	::testing::InitGoogleTest(&argc, argv);

	return RUN_ALL_TESTS();
}