OpenVulkano/openVulkanoCpp/Scene/Ray.cpp

/*
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/.
 */

#include "Ray.hpp"
#include <vector>

namespace
{
	int SolveQuadraticEquation(float a, float b, float c, float& x0, float& x1)
	{
		const float discr = b * b - 4 * a * c;
		if (discr < 0) return 0;
		if (discr == 0)
		{
			x0 = x1 = (-b) / (2 * a);
			return 1;
		}
		float q = (b > 0) ? -0.5 * (b + std::sqrt(discr)) : -0.5 * (b - std::sqrt(discr));
		x0 = q / a;
		x1 = c / q;
		OpenVulkano::Math::Utils::SortPair(x0, x1);
		return 2;
	}
};

namespace OpenVulkano::Scene
{
	using namespace Math::Utils;

	std::optional<RayHit> Ray::IntersectSphere(const Math::Vector3f& center, float radius) const
	{
		RayHit hitRes;
		if (intersectRaySphere(m_origin, m_dir, center, radius, hitRes.point, hitRes.normal))
		{
			hitRes.distance2 = distance2(m_origin, hitRes.point);
			return hitRes;
		}
		return {};
	}

	std::optional<RayHit> Ray::IntersectTriangle(const Math::Vector3f& v0, const Math::Vector3f& v1,
	                                             const Math::Vector3f& v2) const
	{
		float d;
		if (intersectRayTriangle(m_origin, m_dir, v0, v1, v2, m_baryPos, d) && d >= 0)
		{
			RayHit hitRes;
			hitRes.point = (1.f - m_baryPos.x - m_baryPos.y) * v0 + m_baryPos.x * v1 + m_baryPos.y * v2;
			Math::Vector3f e = v1 - v0;
			Math::Vector3f e2 = v2 - v0;
			// triangle normal. won't work if triangle is smoothly shaded. use other overloaded method instead.
			hitRes.SetDistance(d);
			hitRes.normal = normalize(cross(e, e2));
			return hitRes;
		}
		return {};
	}

	std::optional<RayHit> Ray::IntersectTriangle(const Math::Vector3f& v0, const Math::Vector3f& v1,
	                                             const Math::Vector3f& v2, const Math::Vector3f& n0,
	                                             const Math::Vector3f& n1, const Math::Vector3f& n2) const
	{
		auto hit = IntersectTriangle(v0, v1, v2);
		if (hit && (n0 != n1 || n0 != n2 || n1 != n2))
		{
			// TODO: NEED ACTUAL TESTING
			hit->normal = normalize((1.f - m_baryPos.x - m_baryPos.y) * n0 + m_baryPos.x * n1 + m_baryPos.y * n2);
		}
		return hit;
	}

	std::optional<RayHit> Ray::IntersectQuad(const Math::Vector3f& v0, const Math::Vector3f& v1,
	                                         const Math::Vector3f& v2, const Math::Vector3f& v3) const
	{
		if (const auto hitRes = IntersectTriangle(v0, v1, v2))
		{
			return hitRes;
		}
		if (const auto hitRes = IntersectTriangle(v0, v2, v3))
		{
			return hitRes;
		}
		return {};
	}

	std::optional<RayHit> Ray::IntersectAABB(const Math::AABB2f& bbox) const
	{
		//TODO impl that skips z checks
		return IntersectAABB(Math::AABB({bbox.GetMin(), 0}, {bbox.GetMax(), 0}));
	}

	std::optional<RayHit> Ray::IntersectAABB(const Math::AABB& bbox) const
	{
		RayHit h1, h2;
		const int intersections = this->IntersectAABB(bbox, h1, h2);
		if (intersections == 1) return h1;
		if (intersections == 2) return (h1.distance2 < h2.distance2) ? h1 : h2;
		return {};
	}

	int Ray::IntersectAABB(const Math::AABB& bbox, RayHit& p1, RayHit& p2) const
	{
		auto tmin = (bbox.min - m_origin) / m_dir;
		auto tmax = (bbox.max - m_origin) / m_dir;

		SortPair(tmin.x, tmax.x);
		SortPair(tmin.y, tmax.y);
		SortPair(tmin.z, tmax.z);

		if ((tmin.x > tmax.y) || (tmin.y > tmax.x)) return 0;

		if (tmin.y > tmin.x)
		{
			tmin.x = tmin.y;
		}
		if (tmax.y < tmax.x)
		{
			tmax.x = tmax.y;
		}

		if ((tmin.x > tmax.z) || (tmin.z > tmax.x)) return 0;

		if (tmin.z > tmin.x)
		{
			tmin.x = tmin.z;
		}
		if (tmax.z < tmax.x)
		{
			tmax.x = tmax.z;
		}

		int intersections = 2;
		if (tmin.x < 0)
		{
			if (tmax.x < 0) return 0;
			intersections--;
			tmin.x = tmax.x;
		}

		p1.point = m_origin + tmin.x * m_dir;
		p2.point = m_origin + tmax.x * m_dir;
		p1.distance2 = distance2(m_origin, p1.point);
		p2.distance2 = distance2(m_origin, p2.point);
		p1.normal = p2.normal = Math::Vector3f(0);
		return intersections;
	}

	std::optional<RayHit> Ray::IntersectPlane(const Math::Vector3f& pOrigin, const Math::Vector3f& pNorm) const
	{
		Math::Vector3f norm = normalize(pNorm);
		float d;
		if (intersectRayPlane(m_origin, m_dir, pOrigin, pNorm, d))
		{
			return {{ m_origin + m_dir * d, norm, d, d * d }};
		}
		return {};
	}

	int Ray::IntersectSphere(const Math::Vector3f& center, float radius, RayHit& p1, RayHit& p2) const
	{
		const Math::Vector3f L = m_origin - center;
		const float a = length2(m_dir);
		const float b = 2 * dot(m_dir, L);
		const float c = length2(L) - radius * radius;
		float x1, x2;
		int roots = ::SolveQuadraticEquation(a, b, c, x1, x2);

		if (roots == 0) return 0;
		SortPair(x1, x2);
		bool calcP2 = false;
		if (roots == 1)
		{
			if (x1 < 0) return 0; // ray intersects sphere behind the origin
		}
		else
		{
			if (x1 < 0 && x2 < 0) return 0; // ray intersects sphere behind the origin
			if (x1 >= 0 && x2 >= 0)
			{
				calcP2 = true;
			}
			else
			{
				--roots;
				if (x1 < 0) x1 = x2;
			}
		}
		p1.point = m_origin + x1 * m_dir;
		p1.distance2 = distance2(m_origin, p1.point);
		p1.normal = normalize(p1.point - center);
		if (calcP2)
		{
			p2.point = m_origin + x2 * m_dir;
			p2.distance2 = distance2(m_origin, p2.point);
			p2.normal = normalize(p2.point - center);
		}
		else p2 = p1;
		return roots;
	}

	bool RayHit::operator==(const RayHit& other) const 
	{
		return distance2 == other.distance2 && point == other.point && normal == other.normal;
	}

	bool RayHit::operator!=(const RayHit& other) const { return !(*this == other); }
}