src/triangulate.cpp

#include "src/triangulate.h"
#include "thirdparty/earcut.h"
#include "src/invert.h"

// Make mapbox::earcut work with glm::vec3
template<> struct mapbox::util::nth<0, glm::vec3>
{
	static constexpr float get(const glm::vec3& t) { return t.x; }
};

template<> struct mapbox::util::nth<1, glm::vec3>
{
	static constexpr float get(const glm::vec3& t) { return t.y; }
};

namespace triangulate
{
	struct triangleid
	{
		std::size_t id;
		constexpr auto operator<=>(const triangleid&) const = default;
	};
	//! \brief Index type of vertex
	using index_t = std::uint16_t;
	//! \brief Pair of vertex indices, forming a triangle boundary
	using boundary_t = std::pair<index_t, index_t>;
	//! \brief Additional information about a triangle boundary
	struct triangleinfo
	{
		//! \brief What's the third vertex of this triangle
		index_t third_vertex;
		//! \brief Id of triangle (starting position of the triangle in the result
		//! value of mapbox::earcut)
		triangleid triangleid;
	};
	using triangleinfomap_t = std::map<boundary_t, triangleinfo>;
	struct MergedTriangles
	{
		//! \brief List of merged quadrilaterals
		std::vector<Quadrilateral> quadrilaterals;
		//! \brief Ids of cut triangles
		std::set<triangleid> cutTriangles;
	};
	static MergedTriangles mergeTriangles(
		const std::vector<std::uint16_t>& indices,
		std::vector<glm::vec3>::const_iterator polyBegin);
}

static bool findQuadsToMerge(
	const triangulate::triangleinfomap_t& boundaries,
	auto&& iscut,
	auto&& addMergedQuad
) {
	bool foundQuads = false;
	// Go through triangle boundaries
	for (
		triangulate::triangleinfomap_t::const_iterator it1 = boundaries.begin();
		it1 != boundaries.end();
		++it1
	) {
		const triangulate::boundary_t& boundary = it1->first;
		const triangulate::triangleinfo& earcut_triangle_1 = it1->second;
		// .. the ones we haven't already been merged...
		if (not iscut(earcut_triangle_1.triangleid)) {
			// Look for its inverse boundary to find the touching triangle
			const triangulate::boundary_t inverse_boundary = std::make_pair(boundary.second, boundary.first);
			const triangulate::triangleinfomap_t::const_iterator it2 = boundaries.find(inverse_boundary);
			// Also if that hasn't been cut
			if (it2 != boundaries.end() and not iscut(it2->second.triangleid)) {
				// Continue until no more new quads are found
				foundQuads |= addMergedQuad(it1, it2);
			}
		}
	}
	return foundQuads;
}

static triangulate::MergedTriangles triangulate::mergeTriangles(
	const std::vector<std::uint16_t>& indices,
	std::vector<glm::vec3>::const_iterator polyBegin)
{
	MergedTriangles result;
	const auto iscut = [&result](const triangleid i){
		return result.cutTriangles.find(i) != result.cutTriangles.end();
	};
	const auto addMergedQuad = [&result, polyBegin](
		triangleinfomap_t::const_iterator it1,
		triangleinfomap_t::const_iterator it2
	) {
		const boundary_t& shared_boundary = it1->first;
		const triangleinfo& earcut_triangle_1 = it1->second;
		const triangleinfo& earcut_triangle_2 = it2->second;
		const Quadrilateral quad{
			*(polyBegin + shared_boundary.first),
			*(polyBegin + earcut_triangle_2.third_vertex),
			*(polyBegin + shared_boundary.second),
			*(polyBegin + earcut_triangle_1.third_vertex),
		};
		if (quadrilateral_convexity(quad) == convexity_e::convex) {
			result.quadrilaterals.push_back(quad);
			result.cutTriangles.insert(earcut_triangle_1.triangleid);
			result.cutTriangles.insert(earcut_triangle_2.triangleid);
			return true;
		}
		else {
			return false;
		}
	};
	triangleinfomap_t boundaries;
	for (std::size_t i = 0; i < indices.size(); i += 3) {
		const auto add = [&](const std::size_t o1, const std::size_t o2, const std::size_t o3){
			const auto key = std::make_pair(indices[i + o1], indices[i + o2]);
			boundaries[key] = triangleinfo{
				.third_vertex = indices[i + o3],
				.triangleid = {i},
			};
		};
		add(0, 1, 2);
		add(1, 2, 0);
		add(2, 0, 1);
	}
	bool repeat = true;
	while (repeat) {
		repeat = findQuadsToMerge(boundaries, iscut, addMergedQuad);
	}
	return result;
}

//! \brief Polygonize a polygon into triangles and try merge as many of them
//! as possible into quadrilaterals. The polygons must all reside in the same
//! plane.
std::vector<PlainPolygonElement> polygonize(
	std::vector<glm::vec3>::const_iterator begin,
	std::vector<glm::vec3>::const_iterator end)
{
	std::vector<PlainPolygonElement> result;
	// Transform the polygon into XY plane
	opt<glm::mat3> normal = calculateNormal(begin, end);
	if (normal.has_value()) {
		const glm::mat3 normalInverse = glm::inverse(*normal);
		std::vector<std::vector<glm::vec3>> polygons{1};
		std::vector<glm::vec3>& polygon2d = polygons.back();
		polygon2d.reserve(static_cast<std::size_t>(end - begin));
		for (std::vector<glm::vec3>::const_iterator it = begin; it != end; ++it) {
			polygon2d.push_back(normalInverse * glm::vec4{*it, 1});
		}
		const std::vector<triangulate::index_t> indices = mapbox::earcut<triangulate::index_t>(polygons);
		triangulate::MergedTriangles mergedTriangles = triangulate::mergeTriangles(indices, begin);
		for (Quadrilateral& quad : mergedTriangles.quadrilaterals) {
			result.push_back(quad);
		}
		for (std::size_t i = 0; i < indices.size(); i += 3) {
			if (not mergedTriangles.cutTriangles.contains({i})) {
				Triangle tri = Triangle{
					.p1 = *(begin + indices[i]),
					.p2 = *(begin + indices[i + 1]),
					.p3 = *(begin + indices[i + 2]),
				};
				result.push_back(tri);
			}
		}
	}
	return result;
}