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

using indextype = std::uint16_t;
using indexpair = std::pair<indextype, indextype>;
struct boundaryinfo
{
	indextype third;
	std::size_t triangleid;
};
using boundarymap_t = std::map<indexpair, boundaryinfo>;

struct MergedTriangles
{
	std::vector<Quadrilateral> quadrilaterals;
	std::set<std::size_t> cutTriangles;
};

static bool findQuadsToMerge(
	const boundarymap_t& boundaries,
	auto&& iscut,
	auto&& addMergedQuad
) {
	bool foundQuads = false;
	// Go through triangle boundaries
	for (boundarymap_t::const_iterator it1 = boundaries.begin(); it1 != boundaries.end(); ++it1) {
		const indexpair& pair1 = it1->first;
		const boundaryinfo& boundary1 = it1->second;
		// .. the ones we haven't already been merged...
		if (not iscut(boundary1.triangleid)) {
			// Look for its inverse boundary to find the touching triangle
			const indexpair inverse_boundary_key = std::make_pair(pair1.second, pair1.first);
			const boundarymap_t::const_iterator it2 = boundaries.find(inverse_boundary_key);
			// 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 MergedTriangles mergeTriangles(
	const std::vector<std::uint16_t>& indices,
	std::vector<glm::vec3>::const_iterator polyBegin)
{
	MergedTriangles result;
	const auto iscut = [&result](const std::size_t i){
		return result.cutTriangles.find(i) != result.cutTriangles.end();
	};
	const auto addMergedQuad = [&result, polyBegin](
		boundarymap_t::const_iterator it1,
		boundarymap_t::const_iterator it2
	) {
		const indexpair& common_indices = it1->first;
		const boundaryinfo& boundary_1 = it1->second;
		const boundaryinfo& boundary_2 = it2->second;
		const Quadrilateral quad{
			*(polyBegin + common_indices.first),
			*(polyBegin + boundary_2.third),
			*(polyBegin + common_indices.second),
			*(polyBegin + boundary_1.third),
		};
		if (isConvex(quad)) {
			result.quadrilaterals.push_back(quad);
			result.cutTriangles.insert(boundary_1.triangleid);
			result.cutTriangles.insert(boundary_2.triangleid);
			return true;
		}
		else {
			return false;
		}
	};
	boundarymap_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] = {indices[i + o3], 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;
}

std::vector<PlainPolygonElement> polygonize(
	std::vector<glm::vec3>::const_iterator begin,
	std::vector<glm::vec3>::const_iterator end,
	const glm::mat4& planeMatrix)
{
	std::vector<PlainPolygonElement> result;
	const glm::mat4 inverseGrid = glm::inverse(planeMatrix);
	// mapbox::earcut will always produce a CCW polygon, so if we're drawing
	// a CW polygon, we should invert the result afterwards
	const float shouldInvert = glm::dot(
		glm::vec3{inverseGrid[2]},
		glm::cross(*begin - *(begin + 1), *(begin + 2) - *(begin + 1)));
	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(inverseGrid * glm::vec4{*it, 1});
	}
	const std::vector<indextype> indices = mapbox::earcut<std::uint16_t>(polygons);
	MergedTriangles mergedTriangles = mergeTriangles(indices, begin);
	for (Quadrilateral& quad : mergedTriangles.quadrilaterals) {
		if (shouldInvert < 0) {
			invert(quad);
		}
		result.push_back(quad);
	}
	for (std::size_t i = 0; i < indices.size(); i += 3) {
		if (mergedTriangles.cutTriangles.find(i) == mergedTriangles.cutTriangles.end()) {
			Triangle tri = triangle(
				*(begin + indices[i]),
				*(begin + indices[i + 1]),
				*(begin + indices[i + 2]));
			if (shouldInvert < 0) {
				invert(tri);
			}
			result.push_back(tri);
		}
	}
	return result;
}