diff --git a/CMakeLists.txt b/CMakeLists.txt index 73b4a48f..41014cfd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,7 +100,7 @@ endif() set(TARGETS meshoptimizer) if(MESHOPT_BUILD_DEMO) - add_executable(demo demo/main.cpp demo/tests.cpp tools/objloader.cpp) + add_executable(demo demo/main.cpp demo/nanite.cpp demo/tests.cpp tools/objloader.cpp) target_link_libraries(demo meshoptimizer) endif() diff --git a/Makefile b/Makefile index 56b00326..04b2b9f0 100644 --- a/Makefile +++ b/Makefile @@ -35,6 +35,11 @@ ifdef BASISU endif endif +ifdef METIS + $(DEMO_OBJECTS): CXXFLAGS+=-DMETIS + $(DEMO): LDFLAGS+=-lmetis +endif + WASMCC?=$(WASI_SDK)/bin/clang++ WASIROOT?=$(WASI_SDK)/share/wasi-sysroot @@ -108,6 +113,9 @@ check: $(DEMO) dev: $(DEMO) $(DEMO) -d $(files) +nanite: $(DEMO) + $(DEMO) -n $(files) + format: clang-format -i $(LIBRARY_SOURCES) $(DEMO_SOURCES) $(GLTFPACK_SOURCES) diff --git a/demo/main.cpp b/demo/main.cpp index 0507ea91..6af34228 100644 --- a/demo/main.cpp +++ b/demo/main.cpp @@ -138,21 +138,21 @@ Mesh parseObj(const char* path, double& reindex) return result; } -void dumpObj(const Mesh& mesh, bool recomputeNormals = false) +void dumpObj(const std::vector& vertices, const std::vector& indices, bool recomputeNormals = false) { std::vector normals; if (recomputeNormals) { - normals.resize(mesh.vertices.size() * 3); + normals.resize(vertices.size() * 3); - for (size_t i = 0; i < mesh.indices.size(); i += 3) + for (size_t i = 0; i < indices.size(); i += 3) { - unsigned int a = mesh.indices[i], b = mesh.indices[i + 1], c = mesh.indices[i + 2]; + unsigned int a = indices[i], b = indices[i + 1], c = indices[i + 2]; - const Vertex& va = mesh.vertices[a]; - const Vertex& vb = mesh.vertices[b]; - const Vertex& vc = mesh.vertices[c]; + const Vertex& va = vertices[a]; + const Vertex& vb = vertices[b]; + const Vertex& vc = vertices[c]; float nx = (vb.py - va.py) * (vc.pz - va.pz) - (vb.pz - va.pz) * (vc.py - va.py); float ny = (vb.pz - va.pz) * (vc.px - va.px) - (vb.px - va.px) * (vc.pz - va.pz); @@ -160,7 +160,7 @@ void dumpObj(const Mesh& mesh, bool recomputeNormals = false) for (int k = 0; k < 3; ++k) { - unsigned int index = mesh.indices[i + k]; + unsigned int index = indices[i + k]; normals[index * 3 + 0] += nx; normals[index * 3 + 1] += ny; @@ -169,9 +169,9 @@ void dumpObj(const Mesh& mesh, bool recomputeNormals = false) } } - for (size_t i = 0; i < mesh.vertices.size(); ++i) + for (size_t i = 0; i < vertices.size(); ++i) { - const Vertex& v = mesh.vertices[i]; + const Vertex& v = vertices[i]; float nx = v.nx, ny = v.ny, nz = v.nz; @@ -193,9 +193,21 @@ void dumpObj(const Mesh& mesh, bool recomputeNormals = false) fprintf(stderr, "vn %f %f %f\n", nx, ny, nz); } - for (size_t i = 0; i < mesh.indices.size(); i += 3) + for (size_t i = 0; i < indices.size(); i += 3) { - unsigned int a = mesh.indices[i], b = mesh.indices[i + 1], c = mesh.indices[i + 2]; + unsigned int a = indices[i], b = indices[i + 1], c = indices[i + 2]; + + fprintf(stderr, "f %d %d %d\n", a + 1, b + 1, c + 1); + } +} + +void dumpObj(const char* section, const std::vector& indices) +{ + fprintf(stderr, "o %s\n", section); + + for (size_t j = 0; j < indices.size(); j += 3) + { + unsigned int a = indices[j], b = indices[j + 1], c = indices[j + 2]; fprintf(stderr, "f %d %d %d\n", a + 1, b + 1, c + 1); } @@ -905,6 +917,18 @@ void shadow(const Mesh& mesh) (end - start) * 1000); } +static int follow(int* parents, int index) +{ + while (index != parents[index]) + { + int parent = parents[index]; + parents[index] = parents[parent]; + index = parent; + } + + return index; +} + void meshlets(const Mesh& mesh, bool scan) { const size_t max_vertices = 64; @@ -940,6 +964,7 @@ void meshlets(const Mesh& mesh, bool scan) double avg_vertices = 0; double avg_triangles = 0; size_t not_full = 0; + size_t not_connected = 0; for (size_t i = 0; i < meshlets.size(); ++i) { @@ -948,14 +973,37 @@ void meshlets(const Mesh& mesh, bool scan) avg_vertices += m.vertex_count; avg_triangles += m.triangle_count; not_full += m.vertex_count < max_vertices; + + // union-find vertices to check if the meshlet is connected + int parents[max_vertices]; + for (unsigned int j = 0; j < m.vertex_count; ++j) + parents[j] = int(j); + + for (unsigned int j = 0; j < m.triangle_count * 3; ++j) + { + int v0 = meshlet_triangles[m.triangle_offset + j]; + int v1 = meshlet_triangles[m.triangle_offset + j + (j % 3 == 2 ? -2 : 1)]; + + v0 = follow(parents, v0); + v1 = follow(parents, v1); + + parents[v0] = v1; + } + + int roots = 0; + for (unsigned int j = 0; j < m.vertex_count; ++j) + roots += follow(parents, j) == int(j); + + assert(roots != 0); + not_connected += roots > 1; } avg_vertices /= double(meshlets.size()); avg_triangles /= double(meshlets.size()); - printf("Meshlets%c: %d meshlets (avg vertices %.1f, avg triangles %.1f, not full %d) in %.2f msec\n", + printf("Meshlets%c: %d meshlets (avg vertices %.1f, avg triangles %.1f, not full %d, not connected %d) in %.2f msec\n", scan ? 'S' : ' ', - int(meshlets.size()), avg_vertices, avg_triangles, int(not_full), (end - start) * 1000); + int(meshlets.size()), avg_vertices, avg_triangles, int(not_full), int(not_connected), (end - start) * 1000); float camera[3] = {100, 100, 100}; @@ -1155,6 +1203,8 @@ void provoking(const Mesh& mesh) int(mesh.indices.size() / 3), int(pcount), double(pcount) / double(bestv) * 100.0 - 100.0, (end - start) * 1000); } +void nanite(const std::vector& vertices, const std::vector& indices); // nanite.cpp + bool loadMesh(Mesh& mesh, const char* path) { double start = timestamp(); @@ -1344,6 +1394,15 @@ void processDev(const char* path) simplifyAttr(mesh); } +void processNanite(const char* path) +{ + Mesh mesh; + if (!loadMesh(mesh, path)) + return; + + nanite(mesh.vertices, mesh.indices); +} + int main(int argc, char** argv) { void runTests(); @@ -1360,16 +1419,17 @@ int main(int argc, char** argv) if (strcmp(argv[1], "-d") == 0) { for (int i = 2; i < argc; ++i) - { processDev(argv[i]); - } + } + else if (strcmp(argv[1], "-n") == 0) + { + for (int i = 2; i < argc; ++i) + processNanite(argv[i]); } else { for (int i = 1; i < argc; ++i) - { process(argv[i]); - } runTests(); } diff --git a/demo/nanite.cpp b/demo/nanite.cpp new file mode 100644 index 00000000..a52a61a3 --- /dev/null +++ b/demo/nanite.cpp @@ -0,0 +1,584 @@ +// This is a playground for experimenting with algorithms necessary for Nanite like hierarchical clustering +// The code is not optimized, not robust, and not intended for production use. +// It optionally supports METIS for clustering and partitioning, with an eventual goal of removing this code +// in favor of meshopt algorithms. + +// For reference, see the original Nanite paper: +// Brian Karis. Nanite: A Deep Dive. 2021 + +#ifndef _CRT_SECURE_NO_WARNINGS +#define _CRT_SECURE_NO_WARNINGS +#endif + +#include "../src/meshoptimizer.h" + +#include +#include +#include + +#include +#include + +#ifdef METIS +#include +#endif + +#ifndef TRACE +#define TRACE 0 +#endif + +struct Vertex +{ + float px, py, pz; + float nx, ny, nz; + float tx, ty; +}; + +struct LODBounds +{ + float center[3]; + float radius; + float error; +}; + +struct Cluster +{ + std::vector indices; + + LODBounds self; + LODBounds parent; +}; + +const size_t kClusterSize = 128; + +static LODBounds bounds(const std::vector& vertices, const std::vector& indices, float error) +{ + meshopt_Bounds bounds = meshopt_computeClusterBounds(&indices[0], indices.size(), &vertices[0].px, vertices.size(), sizeof(Vertex)); + + LODBounds result; + result.center[0] = bounds.center[0]; + result.center[1] = bounds.center[1]; + result.center[2] = bounds.center[2]; + result.radius = bounds.radius; + result.error = error; + return result; +} + +static LODBounds boundsMerge(const std::vector& clusters, const std::vector& group) +{ + LODBounds result = {}; + + // we approximate merged bounds center as weighted average of cluster centers + // (could also use bounds() center, but we can't use bounds() radius so might as well just merge manually) + float weight = 0.f; + for (size_t j = 0; j < group.size(); ++j) + { + result.center[0] += clusters[group[j]].self.center[0] * clusters[group[j]].self.radius; + result.center[1] += clusters[group[j]].self.center[1] * clusters[group[j]].self.radius; + result.center[2] += clusters[group[j]].self.center[2] * clusters[group[j]].self.radius; + weight += clusters[group[j]].self.radius; + } + + if (weight > 0) + { + result.center[0] /= weight; + result.center[1] /= weight; + result.center[2] /= weight; + } + + // merged bounds must strictly contain all cluster bounds + result.radius = 0.f; + for (size_t j = 0; j < group.size(); ++j) + { + float dx = clusters[group[j]].self.center[0] - result.center[0]; + float dy = clusters[group[j]].self.center[1] - result.center[1]; + float dz = clusters[group[j]].self.center[2] - result.center[2]; + result.radius = std::max(result.radius, clusters[group[j]].self.radius + sqrtf(dx * dx + dy * dy + dz * dz)); + } + + // merged bounds error must be conservative wrt cluster errors + result.error = 0.f; + for (size_t j = 0; j < group.size(); ++j) + result.error = std::max(result.error, clusters[group[j]].self.error); + + return result; +} + +static float boundsError(const LODBounds& bounds, float x, float y, float z) +{ + float dx = bounds.center[0] - x, dy = bounds.center[1] - y, dz = bounds.center[2] - z; + float d = sqrtf(dx * dx + dy * dy + dz * dz) - bounds.radius; + return d <= 0 ? FLT_MAX : bounds.error / d; +} + +#ifdef METIS +static void clusterizeMetisRec(std::vector& result, const std::vector& indices, const std::vector& triidx, const std::vector& triadj) +{ + assert(triadj.size() == triidx.size() * 3); + + if (triidx.size() <= kClusterSize) + { + Cluster cluster; + for (size_t i = 0; i < triidx.size(); ++i) + { + cluster.indices.push_back(indices[triidx[i] * 3 + 0]); + cluster.indices.push_back(indices[triidx[i] * 3 + 1]); + cluster.indices.push_back(indices[triidx[i] * 3 + 2]); + } + + cluster.parent.error = FLT_MAX; + result.push_back(cluster); + return; + } + + std::vector xadj(triidx.size() + 1); + std::vector adjncy; + std::vector part(triidx.size()); + + for (size_t i = 0; i < triidx.size(); ++i) + { + for (int j = 0; j < 3; ++j) + if (triadj[i * 3 + j] != -1) + adjncy.push_back(triadj[i * 3 + j]); + + xadj[i + 1] = adjncy.size(); + } + + int options[METIS_NOPTIONS]; + METIS_SetDefaultOptions(options); + options[METIS_OPTION_SEED] = 42; + options[METIS_OPTION_UFACTOR] = triidx.size() > 8 * kClusterSize ? 100 : 1; + + int nvtxs = int(triidx.size()); + int ncon = 1; + int nparts = 2; + int edgecut = 0; + + int nvtxspad = (nvtxs + kClusterSize - 1) / kClusterSize * kClusterSize; + int idealcut = ((nvtxspad / kClusterSize) / 2) * kClusterSize; + float partw[2] = {idealcut / float(nvtxspad), 1.f - idealcut / float(nvtxspad)}; + + int r = METIS_PartGraphRecursive(&nvtxs, &ncon, &xadj[0], &adjncy[0], NULL, NULL, NULL, &nparts, partw, NULL, options, &edgecut, &part[0]); + assert(r == METIS_OK); + (void)r; + + int partsize[2] = {}; + std::vector partoff(part.size()); + for (size_t i = 0; i < part.size(); ++i) + partoff[i] = partsize[part[i]]++; + + for (int p = 0; p < 2; ++p) + { + std::vector partidx, partadj; + partidx.reserve(partsize[p]); + partadj.reserve(partsize[p] * 3); + + for (size_t i = 0; i < triidx.size(); ++i) + { + if (part[i] != p) + continue; + + partidx.push_back(triidx[i]); + + for (int j = 0; j < 3; ++j) + { + if (triadj[i * 3 + j] >= 0 && part[triadj[i * 3 + j]] == p) + partadj.push_back(partoff[triadj[i * 3 + j]]); + else + partadj.push_back(-1); + } + } + + clusterizeMetisRec(result, indices, partidx, partadj); + } +} + +static std::vector clusterizeMetis(const std::vector& vertices, const std::vector& indices) +{ + std::vector shadowib(indices.size()); + meshopt_generateShadowIndexBuffer(&shadowib[0], &indices[0], indices.size(), &vertices[0].px, vertices.size(), sizeof(float) * 3, sizeof(Vertex)); + + std::map, unsigned int> edges; + + for (size_t i = 0; i < indices.size(); ++i) + { + unsigned int v0 = shadowib[i + 0]; + unsigned int v1 = shadowib[i + (i % 3 == 2 ? -2 : 1)]; + + // we don't track adjacency fully on non-manifold edges for now + edges[std::make_pair(v0, v1)] = unsigned(i / 3); + } + + std::vector triadj(indices.size(), -1); + + for (size_t i = 0; i < indices.size(); i += 3) + { + unsigned int v0 = shadowib[i + 0], v1 = shadowib[i + 1], v2 = shadowib[i + 2]; + + std::map, unsigned int>::iterator oab = edges.find(std::make_pair(v1, v0)); + std::map, unsigned int>::iterator obc = edges.find(std::make_pair(v2, v1)); + std::map, unsigned int>::iterator oca = edges.find(std::make_pair(v0, v2)); + + triadj[i + 0] = oab != edges.end() ? int(oab->second) : -1; + triadj[i + 1] = obc != edges.end() ? int(obc->second) : -1; + triadj[i + 2] = oca != edges.end() ? int(oca->second) : -1; + } + + std::vector triidx(indices.size() / 3); + for (size_t i = 0; i < indices.size(); i += 3) + triidx[i / 3] = int(i / 3); + + std::vector result; + clusterizeMetisRec(result, indices, triidx, triadj); + + return result; +} +#endif + +static std::vector clusterize(const std::vector& vertices, const std::vector& indices) +{ +#ifdef METIS + static const char* metis = getenv("METIS"); + if (metis && atoi(metis) >= 2) + return clusterizeMetis(vertices, indices); +#endif + + const size_t max_vertices = 192; // TODO: depends on kClusterSize, also may want to dial down for mesh shaders + const size_t max_triangles = kClusterSize; + + size_t max_meshlets = meshopt_buildMeshletsBound(indices.size(), max_vertices, max_triangles); + + std::vector meshlets(max_meshlets); + std::vector meshlet_vertices(max_meshlets * max_vertices); + std::vector meshlet_triangles(max_meshlets * max_triangles * 3); + + meshlets.resize(meshopt_buildMeshlets(&meshlets[0], &meshlet_vertices[0], &meshlet_triangles[0], &indices[0], indices.size(), &vertices[0].px, vertices.size(), sizeof(Vertex), max_vertices, max_triangles, 0.f)); + + std::vector clusters(meshlets.size()); + + for (size_t i = 0; i < meshlets.size(); ++i) + { + const meshopt_Meshlet& meshlet = meshlets[i]; + + meshopt_optimizeMeshlet(&meshlet_vertices[meshlet.vertex_offset], &meshlet_triangles[meshlet.triangle_offset], meshlet.triangle_count, meshlet.vertex_count); + + // note: for now we discard meshlet-local indices; they are valuable for shader code so in the future we should bring them back + clusters[i].indices.resize(meshlet.triangle_count * 3); + for (size_t j = 0; j < meshlet.triangle_count * 3; ++j) + clusters[i].indices[j] = meshlet_vertices[meshlet.vertex_offset + meshlet_triangles[meshlet.triangle_offset + j]]; + + clusters[i].parent.error = FLT_MAX; + } + + return clusters; +} + +#ifdef METIS +static std::vector > partitionMetis(const std::vector& clusters, const std::vector& pending) +{ + std::vector > result; + + std::map, std::vector > edges; + + for (size_t i = 0; i < pending.size(); ++i) + { + const Cluster& cluster = clusters[pending[i]]; + + for (size_t j = 0; j < cluster.indices.size(); ++j) + { + int v0 = cluster.indices[j + 0]; + int v1 = cluster.indices[j + (j % 3 == 2 ? -2 : 1)]; + + std::vector& list = edges[std::make_pair(std::min(v0, v1), std::max(v0, v1))]; + if (list.empty() || list.back() != int(i)) + list.push_back(int(i)); + } + } + + std::map, int> adjacency; + + for (std::map, std::vector >::iterator it = edges.begin(); it != edges.end(); ++it) + { + const std::vector& list = it->second; + + for (size_t i = 0; i < list.size(); ++i) + for (size_t j = i + 1; j < list.size(); ++j) + adjacency[std::make_pair(std::min(list[i], list[j]), std::max(list[i], list[j]))]++; + } + + std::vector xadj(pending.size() + 1); + std::vector adjncy; + std::vector adjwgt; + std::vector part(pending.size()); + + for (size_t i = 0; i < pending.size(); ++i) + { + for (std::map, int>::iterator it = adjacency.begin(); it != adjacency.end(); ++it) + if (it->first.first == int(i)) + { + adjncy.push_back(it->first.second); + adjwgt.push_back(it->second); + } + else if (it->first.second == int(i)) + { + adjncy.push_back(it->first.first); + adjwgt.push_back(it->second); + } + + xadj[i + 1] = adjncy.size(); + } + + int options[METIS_NOPTIONS]; + METIS_SetDefaultOptions(options); + options[METIS_OPTION_SEED] = 42; + options[METIS_OPTION_UFACTOR] = 100; + + int nvtxs = int(pending.size()); + int ncon = 1; + int nparts = int(pending.size() + 3) / 4; + int edgecut = 0; + + if (nparts <= 1) + { + // not sure why this is a special case that we need to handle but okay metis + result.push_back(pending); + } + else + { + int r = METIS_PartGraphKway(&nvtxs, &ncon, &xadj[0], &adjncy[0], NULL, NULL, &adjwgt[0], &nparts, NULL, NULL, options, &edgecut, &part[0]); + assert(r == METIS_OK); + (void)r; + + result.resize(nparts); + for (size_t i = 0; i < part.size(); ++i) + result[part[i]].push_back(pending[i]); + } + + return result; +} +#endif + +static std::vector > partition(const std::vector& clusters, const std::vector& pending) +{ +#ifdef METIS + static const char* metis = getenv("METIS"); + if (metis && atoi(metis) >= 1) + return partitionMetis(clusters, pending); +#endif + + std::vector > result; + + size_t last_indices = 0; + + // rough merge; while clusters are approximately spatially ordered, this should use a proper partitioning algorithm + for (size_t i = 0; i < pending.size(); ++i) + { + if (result.empty() || last_indices + clusters[pending[i]].indices.size() > kClusterSize * 4 * 3) + { + result.push_back(std::vector()); + last_indices = 0; + } + + result.back().push_back(pending[i]); + last_indices += clusters[pending[i]].indices.size(); + } + + return result; +} + +static std::vector simplify(const std::vector& vertices, const std::vector& indices, size_t target_count, float* error = NULL) +{ + if (target_count > indices.size()) + return indices; + + std::vector lod(indices.size()); + unsigned int options = meshopt_SimplifyLockBorder | meshopt_SimplifySparse | meshopt_SimplifyErrorAbsolute; + lod.resize(meshopt_simplify(&lod[0], &indices[0], indices.size(), &vertices[0].px, vertices.size(), sizeof(Vertex), target_count, FLT_MAX, options, error)); + + return lod; +} + +void dumpObj(const std::vector& vertices, const std::vector& indices, bool recomputeNormals = false); +void dumpObj(const char* section, const std::vector& indices); + +void nanite(const std::vector& vertices, const std::vector& indices) +{ + static const char* metis = getenv("METIS"); + if (metis && atoi(metis)) + { +#ifdef METIS + printf("using metis for %s\n", atoi(metis) >= 2 ? "both clustering and partition" : "partition only"); +#else + printf("ERROR: build does not have metis available\n"); +#endif + } + + static const char* dump = getenv("DUMP"); + + // initial clusterization splits the original mesh + std::vector clusters = clusterize(vertices, indices); + for (size_t i = 0; i < clusters.size(); ++i) + clusters[i].self = bounds(vertices, clusters[i].indices, 0.f); + + printf("lod 0: %d clusters, %d triangles\n", int(clusters.size()), int(indices.size() / 3)); + + std::vector pending(clusters.size()); + for (size_t i = 0; i < clusters.size(); ++i) + pending[i] = int(i); + +#ifndef NDEBUG + std::vector > dag_debug; +#endif + + int depth = 0; + + // merge and simplify clusters until we can't merge anymore + while (pending.size() > 1) + { + std::vector > groups = partition(clusters, pending); + pending.clear(); + + std::vector retry; + + size_t triangles = 0; + size_t stuck_triangles = 0; + int single_clusters = 0; + int stuck_clusters = 0; + int full_clusters = 0; + + if (dump && depth == atoi(dump)) + dumpObj(vertices, std::vector()); + + // every group needs to be simplified now + for (size_t i = 0; i < groups.size(); ++i) + { + if (groups[i].empty()) + continue; // metis shortcut + + if (groups[i].size() == 1) + { +#if TRACE + printf("stuck cluster: singleton with %d triangles\n", int(clusters[groups[i][0]].indices.size() / 3)); +#endif + + if (dump && depth == atoi(dump)) + dumpObj("cluster", clusters[groups[i][0]].indices); + + single_clusters++; + stuck_clusters++; + stuck_triangles += clusters[groups[i][0]].indices.size() / 3; + retry.push_back(groups[i][0]); + continue; + } + + std::vector merged; + for (size_t j = 0; j < groups[i].size(); ++j) + merged.insert(merged.end(), clusters[groups[i][j]].indices.begin(), clusters[groups[i][j]].indices.end()); + + if (dump && depth == atoi(dump)) + dumpObj("group", merged); + + float error = 0.f; + std::vector simplified = simplify(vertices, merged, kClusterSize * 2 * 3, &error); + if (simplified.size() > merged.size() * 0.85f || simplified.size() > kClusterSize * 3 * 3) + { +#if TRACE + printf("stuck cluster: simplified %d => %d over threshold\n", int(merged.size() / 3), int(simplified.size() / 3)); +#endif + stuck_clusters++; + stuck_triangles += merged.size() / 3; + for (size_t j = 0; j < groups[i].size(); ++j) + retry.push_back(groups[i][j]); + continue; // simplification is stuck; abandon the merge + } + + // enforce bounds and error monotonicity + // note: it is incorrect to use the precise bounds of the merged or simplified mesh, because this may violate monotonicity + LODBounds groupb = boundsMerge(clusters, groups[i]); + groupb.error += error; // this may overestimate the error, but we are starting from the simplified mesh so this is a little more correct + + std::vector split = clusterize(vertices, simplified); + + // update parent bounds and error for all clusters in the group + // note that all clusters in the group need to switch simultaneously so they have the same bounds + for (size_t j = 0; j < groups[i].size(); ++j) + { + assert(clusters[groups[i][j]].parent.error == FLT_MAX); + clusters[groups[i][j]].parent = groupb; + } + +#ifndef NDEBUG + // record DAG edges for validation during the cut + for (size_t j = 0; j < groups[i].size(); ++j) + for (size_t k = 0; k < split.size(); ++k) + dag_debug.push_back(std::make_pair(groups[i][j], int(clusters.size()) + int(k))); +#endif + + for (size_t j = 0; j < split.size(); ++j) + { + split[j].self = groupb; + + clusters.push_back(split[j]); // std::move + pending.push_back(int(clusters.size()) - 1); + + triangles += split[j].indices.size() / 3; + full_clusters += split[j].indices.size() == kClusterSize * 3; + } + } + + depth++; + printf("lod %d: simplified %d clusters (%d full, %.1f tri/cl), %d triangles; stuck %d clusters (%d single), %d triangles\n", depth, + int(pending.size()), full_clusters, pending.empty() ? 0 : double(triangles) / double(pending.size()), int(triangles), stuck_clusters, single_clusters, int(stuck_triangles)); + + if (triangles < stuck_triangles / 3) + break; + + pending.insert(pending.end(), retry.begin(), retry.end()); + } + + size_t lowest_triangles = 0; + for (size_t i = 0; i < clusters.size(); ++i) + if (clusters[i].parent.error == FLT_MAX) + lowest_triangles += clusters[i].indices.size() / 3; + + printf("lowest lod: %d triangles\n", int(lowest_triangles)); + + // for testing purposes, we can compute a DAG cut from a given viewpoint and dump it as an OBJ + float maxx = 0.f, maxy = 0.f, maxz = 0.f; + for (size_t i = 0; i < vertices.size(); ++i) + { + maxx = std::max(maxx, vertices[i].px * 2); + maxy = std::max(maxy, vertices[i].py * 2); + maxz = std::max(maxz, vertices[i].pz * 2); + } + float threshold = 3e-3f; + + std::vector cut; + for (size_t i = 0; i < clusters.size(); ++i) + if (boundsError(clusters[i].self, maxx, maxy, maxz) <= threshold && boundsError(clusters[i].parent, maxx, maxy, maxz) > threshold) + cut.insert(cut.end(), clusters[i].indices.begin(), clusters[i].indices.end()); + +#ifndef NDEBUG + for (size_t i = 0; i < dag_debug.size(); ++i) + { + int j = dag_debug[i].first, k = dag_debug[i].second; + float ej = boundsError(clusters[j].self, maxx, maxy, maxz); + float ejp = boundsError(clusters[j].parent, maxx, maxy, maxz); + float ek = boundsError(clusters[k].self, maxx, maxy, maxz); + + assert(ej <= ek); + assert(ejp >= ej); + } +#endif + + printf("cut (%.3f): %d triangles\n", threshold, int(cut.size() / 3)); + + if (dump && -1 == atoi(dump)) + { + dumpObj(vertices, cut); + + for (size_t i = 0; i < clusters.size(); ++i) + if (boundsError(clusters[i].self, maxx, maxy, maxz) <= threshold && boundsError(clusters[i].parent, maxx, maxy, maxz) > threshold) + dumpObj("cluster", clusters[i].indices); + } +}