Add geom::subdivide for subdividing an indexed mesh

libs/geom/include/psemek/geom/mesh.hpp

@@ -156,4 +156,56 @@ namespace psemek::geom
 		return result;
 	}
 
+	template <typename Vertex, typename Index>
+	void subdivide(std::vector<Vertex> & vertices, std::vector<triangle<Index>> & triangles)
+	{
+		using edge = std::pair<Index, Index>;
+
+		using scalar = std::decay_t<decltype(vertices[0][0])>;
+
+		auto normalise = [](edge e) -> edge
+		{
+			if (e.first > e.second)
+				std::swap(e.first, e.second);
+			return e;
+		};
+
+		std::unordered_map<edge, Index> midpoints;
+
+		auto add_midpoint = [&](Index i0, Index i1)
+		{
+			auto ne = normalise({i0, i1});
+
+			if (!midpoints.contains(ne))
+			{
+				midpoints[ne] = vertices.size();
+				vertices.push_back(lerp(vertices[i0], vertices[i1], scalar{0.5}));
+			}
+		};
+
+		for (auto const & t : triangles)
+		{
+			add_midpoint(t[0], t[1]);
+			add_midpoint(t[1], t[2]);
+			add_midpoint(t[2], t[0]);
+		}
+
+		auto get_midpoint = [&](Index i0, Index i1)
+		{
+			return midpoints.at(normalise({i0, i1}));
+		};
+
+		std::vector<triangle<Index>> new_triangles;
+
+		for (auto const & t : triangles)
+		{
+			new_triangles.push_back({t[0], get_midpoint(t[0], t[1]), get_midpoint(t[2], t[0])});
+			new_triangles.push_back({t[1], get_midpoint(t[1], t[2]), get_midpoint(t[0], t[1])});
+			new_triangles.push_back({t[2], get_midpoint(t[2], t[0]), get_midpoint(t[1], t[2])});
+			new_triangles.push_back({get_midpoint(t[0], t[1]), get_midpoint(t[1], t[2]), get_midpoint(t[2], t[0])});
+		}
+
+		triangles.swap(new_triangles);
+	}
+
 }