src/geometry.cpp

#include <glm/gtc/matrix_transform.hpp>
#include "geometry.h"
#include "main.h"
#include "ring.h"

/**
 * @brief Computes line-plane intersection
 * @param line
 * @param plane
 * @return point of intersection. Does not return a value if the line is in parallel to the plane.
 */
std::optional<glm::vec3> geom::linePlaneIntersection(
	const geom::Line<3>& line,
	const geom::Plane& plane,
	const float epsilon)
{
	const float denominator = glm::dot(line.direction, plane.normal);
	if (std::abs(denominator) < epsilon)
	{
		return {};
	}
	else
	{
		const float d = glm::dot(plane.anchor - line.anchor, plane.normal) / denominator;
		return line.anchor + d * line.direction;
	}
}

/**
 * @brief Computes the plane of a triangle
 * @param triangle
 * @return plane
 */
geom::Plane geom::planeFromTriangle(const geom::Triangle& triangle)
{
	return geom::Plane{normalVector(triangle), triangle.points[0]};
}

/**
 * @brief Computes the normal vector of a triangle
 * @param triangle
 * @return normal vector
 */
glm::vec3 geom::normalVector(const geom::Triangle& triangle)
{
	return glm::normalize(
		glm::cross(
			triangle.points[1] - triangle.points[0],
			triangle.points[2] - triangle.points[0]));
}

/**
 * @brief Extracts the scaling component of the specified matrix into a vector and returns both the scaling
 * components as well as the unscaled matrix.
 * @param matrix Matrix to compute
 * @return scaling vector and unscaled matrix
 */
geom::ScalingExtract geom::extractScaling(const glm::mat4& matrix)
{
	geom::ScalingExtract result;
	result.scaling = geom::scalingVector(matrix);
	result.unscaled = glm::scale(matrix, 1.0f / result.scaling);
	return result;
}

/**
 * @brief Computes the scaling vector, which contains the scaling of the specified matrix
 * @param matrix
 * @return scaling vector
 */
glm::vec3 geom::scalingVector(const glm::mat4 matrix)
{
	auto component = [](const glm::mat4& matrix, const int i) -> float
	{
		return std::hypot(std::hypot(matrix[i][0], matrix[i][1]), matrix[i][2]);
	};
	return glm::vec3{component(matrix, 0), component(matrix, 1), component(matrix, 2)};
}

std::optional<glm::vec2> geom::lineLineIntersection(const Line<2>& line_1, const Line<2>& line_2)
{
	const float denominator = (line_1.direction.x * line_2.direction.y) - (line_1.direction.y * line_2.direction.x);
	constexpr float epsilon = 1e-6f;
	if (std::abs(denominator) < epsilon)
	{
		return {};
	}
	else
	{
		const glm::vec2 p1 = line_1.anchor + line_1.direction;
		const glm::vec2& p2 = line_1.anchor;
		const glm::vec2 p3 = line_2.anchor + line_2.direction;
		const glm::vec2& p4 = line_2.anchor;
		const float a = glm::determinant(glm::mat2{{p1.x, p2.x}, {p1.y, p2.y}});
		const float b = glm::determinant(glm::mat2{{p3.x, p4.x}, {p3.y, p4.y}});
		const float c_x = glm::determinant(glm::mat2{{p1.x, p2.x}, {1, 1}});
		const float c_y = glm::determinant(glm::mat2{{p1.y, p2.y}, {1, 1}});
		const float d_x = glm::determinant(glm::mat2{{p3.x, p4.x}, {1, 1}});
		const float d_y = glm::determinant(glm::mat2{{p3.y, p4.y}, {1, 1}});
		const float x = glm::determinant(glm::mat2{{a, b}, {c_x, d_x}});
		const float y = glm::determinant(glm::mat2{{a, b}, {c_y, d_y}});
		return glm::vec2{x / denominator, y / denominator};
	}
}

std::optional<glm::vec2> geom::rayLineSegmentIntersection(const Ray<2>& ray, const LineSegment2D& line)
{
	std::optional<glm::vec2> result = lineLineIntersection(
		rayToLine(ray),
		lineFromPoints(line.points[0], line.points[1]));
	if (result.has_value())
	{
		const float d1 = glm::dot(*result - ray.anchor, ray.direction);
		if (d1 < 0)
		{
			result.reset();
		}
		else
		{
			const float d2 = glm::dot(*result - line.points[0], *result - line.points[1]);
			if (d2 > 0)
			{
				result.reset();
			}
		}
	}
	return result;
}

std::optional<geom::PointOnRectagle> geom::rayRectangleIntersection(const Ray<2>& ray, const QRectF& rectangle)
{
	std::optional<glm::vec2> position;
	std::optional<PointOnRectagle> result;
	// Try top
	position = rayLineSegmentIntersection(ray, top(rectangle));
	if (position.has_value())
	{
		result = {*position, RectangleSide::Top};
	}
	else
	{
		// Try bottom
		position = rayLineSegmentIntersection(ray, bottom(rectangle));
		if (position.has_value())
		{
			result = {*position, RectangleSide::Bottom};
		}
		else
		{
			// Try left
			position = rayLineSegmentIntersection(ray, left(rectangle));
			if (position.has_value())
			{
				result = {*position, RectangleSide::Left};
			}
			else
			{
				// Try right
				position = rayLineSegmentIntersection(ray, right(rectangle));
				if (position.has_value())
				{
					result = {*position, RectangleSide::Right};
				}
			}
		}
	}
	return result;
}

geom::LineSegment2D geom::top(const QRectF& rectangle)
{
	return {
		glm::vec2{rectangle.left(), rectangle.top()},
		glm::vec2{rectangle.right(), rectangle.top()}
	};
}

geom::LineSegment2D geom::bottom(const QRectF& rectangle)
{
	return {
		glm::vec2{rectangle.left(), rectangle.bottom()},
		glm::vec2{rectangle.right(), rectangle.bottom()}
	};
}

geom::LineSegment2D geom::left(const QRectF& rectangle)
{
	return {
		glm::vec2{rectangle.left(), rectangle.top()},
		glm::vec2{rectangle.left(), rectangle.bottom()}
	};
}

geom::LineSegment2D geom::right(const QRectF& rectangle)
{
	return {
		glm::vec2{rectangle.right(), rectangle.top()},
		glm::vec2{rectangle.right(), rectangle.bottom()}
	};
}

bool geom::isConvex(const std::vector<glm::vec3>& polygon)
{
	const int n = polygon.size();
	auto polygonRing = iter::ring(polygon, n);
	std::vector<glm::vec3> crosses;
	crosses.resize(n);
	for (int i = 0; i < n; i += 1)
	{
		crosses[i] = glm::cross(polygonRing[i - 1] - polygonRing[i], polygonRing[i + 1] - polygonRing[i]);
	}
	return not std::any_of(
		crosses.begin() + 1,
		crosses.end(),
		[&crosses](const glm::vec3& vector)
		{
			return glm::dot(crosses[0], vector) < 1e-6;
		});
}

/**
 * @brief Determines the winding of a 2d polygon
 * @param polygon
 * @return winding
 */
Winding geom::winding(const QPolygonF &polygon)
{
	// based on https://stackoverflow.com/a/1165943
	double sum = 0.0;
	for (int i = 0; i < polygon.size(); i += 1)
	{
		const QPointF& p1 = polygon[i];
		const QPointF& p2 = polygon[(i + 1) % polygon.size()];
		sum += (p2.x() - p1.x()) * (p2.y() + p1.y());
	}
	return (sum < 0) ? Winding::Anticlockwise : Winding::Clockwise;
}

/**
 * @brief computes the point on a Bezier curve
 * @param curve
 * @param t scalar between 0 and 1, with t=0 being P0 and t=1 being P3
 * @return point on curve
 */
glm::vec3 geom::pointOnCurve(const BezierCurve &curve, float t)
{
	// clamp t as rounding errors might make it slightly out of bounds
	t = std::clamp(t, 0.0f, 1.0f);
	const float t_2 = t * t;
	const float t_3 = t * t * t;
	const float coeffs[3] = {
		-1*t_3  +3*t_2  -3*t  +1,
		+3*t_3  -6*t_2  +3*t,
		-3*t_3  +3*t_2,
	};
	return coeffs[0] * curve[0] + coeffs[1] * curve[1] + coeffs[2] * curve[2] + t_3 * curve[3];
}

/**
 * @brief computes the derivative of a point on a Bezier curve
 * @param curve
 * @param t scalar between 0 and 1, with t=0 being P0 and t=1 being P3
 * @return point on curve
 */
glm::vec3 geom::derivativeOnCurve(const BezierCurve &curve, float t)
{
	// clamp t as rounding errors might make it slightly out of bounds
	t = std::clamp(t, 0.0f, 1.0f);
	const float t_2 = t * t;
	const float coeffs[4] = {
		-3*t_2  + 6*t  -3,
		+9*t_2  -12*t  +3,
		-9*t_2  + 6*t,
		+3*t_2
	};
	return coeffs[0] * curve[0] + coeffs[1] * curve[1] + coeffs[2] * curve[2] + coeffs[3] * curve[3];
}