#define _USE_MATH_DEFINES 
#include <iostream>
#include <math.h>
#include <random>

struct vec3
{
	float x, y, z;

	vec3 operator+(const vec3& other)
	{
		return { x + other.x, y + other.y, z + other.z };
	}

	vec3 operator-(const vec3& other)
	{
		return { x - other.x, y - other.y, z - other.z };
	}

	vec3 operator-()
	{
		return { -x, -y, -z };
	}

	vec3 operator*(const float val)
	{
		return { x * val, y * val,z * val };
	}

	vec3 operator/(const vec3& other)
	{
		return { x / other.x, y / other.y, z / other.z };
	}

	vec3 operator/(const float val)
	{
		return { x / val, y / val,z / val };
	}

	vec3 normalize()
	{
		float norm = std::sqrt(x * x + y * y + z * z);
		return { x / norm, y / norm, z / norm };
	}
};

inline vec3 operator*(const vec3& a, const float b)
{
	return { a.x * b, a.y * b, a.z * b };
}

inline vec3 operator-(const vec3& a, const vec3& b)
{
	return { a.x - b.x, a.y - b.y, a.z - b.z };
}

float dot(const vec3& a, const vec3& b)
{
	return { a.x * b.x + a.y * b.y + a.z * b.z };
}

vec3 cross(const vec3& a, const vec3& b)
{
	return { a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x };
}

vec3 min(const vec3& a, const vec3& b)
{
	return { std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z) };
}

float min_component(const vec3& a)
{
	return std::min(a.x, std::min(a.y, a.z));
}

vec3 max(const vec3& a, const vec3& b)
{
	return { std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z) };
}

float max_component(const vec3& a)
{
	return std::max(a.x, std::max(a.y, a.z));
}

vec3 abs(const vec3& a)
{
	return { std::abs(a.x), std::abs(a.y), std::abs(a.z) };
}

float length(const vec3& a)
{
	return std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
}

vec3 cos(const vec3& a)
{
	return { std::cos(a.x), std::cos(a.y), std::cos(a.z) };
}

vec3 sin(const vec3& a)
{
	return { std::sin(a.x), std::sin(a.y), std::sin(a.z) };
}

constexpr size_t renderWidth = 1024;
constexpr size_t renderHeight = 1024;

constexpr vec3 startBounding{ -M_PI, -M_PI, -M_PI };
constexpr vec3 endBounding{ 399.f * M_PI, 399.f * M_PI, 99.f * M_PI };

inline float evalBoundingBox(const vec3& p)
{
	vec3 q = abs(p - startBounding * .5f - endBounding * .5f) - (endBounding * .5f - startBounding * .5f);
	return length(max(q, { .0f, .0f, .0f })) + std::min(.0f, max_component(q));
}

inline float eval(const vec3& p)
{
	float t0 = evalBoundingBox(p);
	float t1 = std::abs(std::cos(p.x) + std::cos(p.y) + std::cos(p.z));	// inverted sign for keeping S(p) > 0 outside the surface
	return std::max(t0, t1);
}

inline float evalTest(const vec3& p)
{
	return -(std::cos(p.x) + std::cos(p.y) + std::cos(p.z));
}

inline bool testIntersection(const vec3& p, float scaling = 1.001f)
{
	float t0 = -(std::cos(p.x) + std::cos(p.y) + std::cos(p.z));
	float t1 = -(std::cos(p.x + t0 * scaling) + std::cos(p.y + t0 * scaling) + std::cos(p.z + t0 * scaling));
	return t0 * t1 < .0f;
}

inline std::pair<float, float> eval(const vec3& p, const vec3& dir)
{
	vec3 cosf = cos(p);
	vec3 sinf = sin(p);
	float t0 = evalBoundingBox(p);
	return { std::max(t0, -(cosf.x + cosf.y + cosf.z)), dot(-sinf, dir) };
}

int main(int argc, char** argv)
{
	vec3 eye = { 199.f * M_PI, 199.f * M_PI, -50.f },
		center = { 199.f * M_PI, 199.f * M_PI, -1.f },
		up = { .0f, 1.f, .0f };
	float fov = 60.f;
	float near = .001f, far = 400.f * M_PI;

	vec3 direction = (center - eye).normalize();
	vec3 right = cross(up, direction);

	std::default_random_engine e(time(nullptr));
	std::uniform_real_distribution<float> u(0, 1);
	std::pair<float, float> uv = { u(e) * 2 - 1, u(e) * 2 - 1 };

	vec3 rayDir = (direction + right * uv.first + up * uv.second).normalize();
	
	float proj = 1.f / dot(rayDir, direction);
	vec3 t0 = (-eye + startBounding) / rayDir;
	vec3 t1 = (-eye + endBounding) / rayDir;
	vec3 tmin = min(t0, t1), tmax = max(t0, t1);
	float tNear = std::max(max_component(tmin), near * proj), tFar = std::min(min_component(tmax), far * proj);

	float t = tNear;
	int i = 0;
	bool hit = false;
	vec3 p, endPoint;
	std::vector<float> ST_t;
	while (t < tFar && tNear < tFar && tFar > .0f && i <= 1024)
	{
		p = eye + rayDir * t;
		float val = eval(p);
		//bool test = testIntersection(p);

		endPoint = eye + rayDir * (t + val);
		auto evalStart = eval(p, rayDir);
		auto evalEnd = eval(endPoint, rayDir);

		std::initializer_list testPoints = { evalStart.first, evalStart.first + evalStart.second * val * .5f, evalEnd.first, evalEnd.first + evalEnd.second * val * .5f };
		float test[] = { std::min(testPoints), std::max(testPoints) };

		if (val < .001f)
		{
			vec3 temp = p, temp2 = endPoint;
			for (int i = 0; i < 0; ++i)
			{
				vec3 midPoint = (p + endPoint) * .5f;
				auto evalMid = eval(midPoint, rayDir);
				std::initializer_list testPoints = { evalStart.first, evalStart.first + evalStart.second * val * .5f, evalMid.first, evalMid.first + evalMid.second * val * .5f };
				float test[] = { std::min(testPoints), std::max(testPoints) };
				bool testRes = test[0] * test[1] < .0f;
				p = testRes ? p : midPoint;
				evalStart = testRes ? evalStart : evalMid;
				endPoint = testRes ? midPoint : endPoint;
				evalEnd = testRes ? evalMid : evalEnd;
				val *= .5f;
			}
			p = (p + endPoint) * .5f;

			std::cout << "Hit after " << i << " steps with p = (" << p.x << ", " << p.y << ", " << p.z << ")!!!" << std::endl;
			std::cout << "\tOriginal p = (" << temp.x << ", " << temp.y << ", " << temp.z << ")" << std::endl;
			std::cout << "\tOriginal endPos = (" << temp2.x << ", " << temp2.y << ", " << temp2.z << ")" << std::endl;
			std::cout << "\tError: " << evalTest(temp2) << std::endl;
			std::cout << "\tError with p: " << evalTest(temp) << std::endl;
			hit = true;
			break;
		}

		t += val;
		i++;

		ST_t.push_back(t);
	}

	if (!hit)
		std::cout << "No hit!" << std::endl;

	t = tNear;
	float s;
	i = 0;
	vec3 startPoint;
	std::vector<float> RM_t;
	while (t < tFar && tNear < tFar && tFar > .0f)
	{
		startPoint = eye + rayDir * t;
		auto evalStart = eval(startPoint, rayDir);
		if (abs(evalStart.first) > .02f)
			s = .01f * 2;
		else if (abs(evalStart.first) < .01f)
		{
			if (abs(evalStart.second) < .005f)
				s = .01f * .25f;
			else
				s = .01f * .5f;
		}
		else
			s = .01f;

		vec3 endPoint = eye + rayDir * (t + s);
		auto evalEnd = eval(endPoint, rayDir);

		std::initializer_list testPoints = { evalStart.first, evalStart.first + evalStart.second * s * .5f, evalEnd.first, evalEnd.first + evalEnd.second * s * .5f };
		float test[] = { std::min(testPoints), std::max(testPoints) };

		if (test[0] * test[1] < .0f)
		{
			std::cout << "Hit after " << i << " steps with p = (" << startPoint.x << ", " << startPoint.y << ", " << startPoint.z << ")!!!" << std::endl;
			std::cout << "\tError: " << evalTest(startPoint) << std::endl;
			break;
		}

		t += s;
		i++;

		RM_t.push_back(t);
	}

	std::cout << "Error: " << length(abs(startPoint - endPoint)) << std::endl;
	if (length(abs(startPoint - endPoint)) > .01f)
	{
		std::cout << "Trace(Sphere Tracing): " << std::endl;
		std::cout << "\tRay steps: ";
		for (auto& i : ST_t)
			std::cout << i << ", ";
		std::cout << std::endl << std::endl;

		std::cout << "Trace(Nvidia's Ray Marching): " << std::endl;
		std::cout << "\tRay steps: ";
		for (auto& i : RM_t)
			std::cout << i << ", ";
		std::cout << std::endl << std::endl;
	}

	return 0;
}