diff -r 216f02b50b0a -r 9727e545b0bc src/triangulate.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/triangulate.cpp	Sun Jul 03 13:44:11 2022 +0300
@@ -0,0 +1,143 @@
+#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;
+}