src/geometry.cpp

#include <glm/gtc/matrix_transform.hpp>
#include "geometry.h"
#include "basics.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> linePlaneIntersection(
	const Line<3>& line,
	const 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
 */
Plane planeFromTriangle(const Triangle& triangle)
{
	return Plane{normalVector(triangle), triangle.p1};
}

/**
 * @brief Computes the normal vector of a triangle
 * @param triangle
 * @return normal vector
 */
glm::vec3 normalVector(const Triangle& triangle)
{
	return glm::normalize(
		glm::cross(
			triangle.p2 - triangle.p1,
			triangle.p3 - triangle.p1));
}

/**
 * @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
 */
ScalingExtract extractScaling(const glm::mat4& matrix)
{
	ScalingExtract result;
	result.scaling = 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 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> 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> rayLineSegmentIntersection(const Ray<2>& ray, const LineSegment2D& line)
{
	std::optional<glm::vec2> result = lineLineIntersection(
		rayToLine(ray),
		lineFromPoints(line.p1, line.p2));
	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.p1, *result - line.p2);
			if (d2 > 0)
			{
				result.reset();
			}
		}
	}
	return result;
}

std::optional<PointOnRectagle> 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;
}

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

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

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

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

bool isConvex(const Quadrilateral& quad)
{
	glm::vec3 crosses[4] = {
		glm::cross(quad.p4 - quad.p1, quad.p2 - quad.p1),
		glm::cross(quad.p1 - quad.p2, quad.p3 - quad.p2),
		glm::cross(quad.p2 - quad.p3, quad.p4 - quad.p3),
		glm::cross(quad.p3 - quad.p4, quad.p1 - quad.p4),
	};
	return not std::any_of(
		&crosses[1],
		&crosses[4],
		[&crosses](const glm::vec3& vector) {
			return glm::dot(crosses[0], vector) < 1e-6;
		});
}

bool 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 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 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 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];
}