#include "dc_simplify.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#include <random>
// -------------------------------------------------------------------------
// Symmetric 4x4 QEM matrix (stored as 10 unique values)
// -------------------------------------------------------------------------
struct QEM {
// Upper triangle of the symmetric 4x4 matrix:
// [ a00 a01 a02 a03 ]
// [ a11 a12 a13 ]
// [ a22 a23 ]
// [ a33 ]
float m[10] = {};
QEM() = default;
// Create QEM from a plane (nx, ny, nz, d) where nx*x + ny*y + nz*z + d = 0
QEM(float nx, float ny, float nz, float d) {
m[0] = nx * nx; m[1] = nx * ny; m[2] = nx * nz; m[3] = nx * d;
m[4] = ny * ny; m[5] = ny * nz; m[6] = ny * d;
m[7] = nz * nz; m[8] = nz * d;
m[9] = d * d;
}
QEM operator+(const QEM& o) const {
QEM r;
for (int i = 0; i < 10; i++) r.m[i] = m[i] + o.m[i];
return r;
}
QEM& operator+=(const QEM& o) {
for (int i = 0; i < 10; i++) m[i] += o.m[i];
return *this;
}
// Evaluate error at position (x, y, z): v^T Q v
float evaluate(float x, float y, float z) const {
// v = [x, y, z, 1]
return m[0]*x*x + 2*m[1]*x*y + 2*m[2]*x*z + 2*m[3]*x
+ m[4]*y*y + 2*m[5]*y*z + 2*m[6]*y
+ m[7]*z*z + 2*m[8]*z
+ m[9];
}
};
// -------------------------------------------------------------------------
// Edge hash for deduplication
// -------------------------------------------------------------------------
struct Edge {
int v0, v1; // v0 < v1
bool operator==(const Edge& o) const { return v0 == o.v0 && v1 == o.v1; }
};
struct EdgeHash {
size_t operator()(const Edge& e) const {
return std::hash<int64_t>()(
(static_cast<int64_t>(e.v0) << 32) | static_cast<int64_t>(e.v1));
}
};
// -------------------------------------------------------------------------
// dc_simplify_mca
// -------------------------------------------------------------------------
int dc_simplify_mca(DCMeshData& mesh, int target_tri_count,
int k_candidates, float max_error) {
int vert_count = static_cast<int>(mesh.positions.size()) / 3;
int tri_count = static_cast<int>(mesh.indices.size()) / 3;
if (tri_count <= target_tri_count || vert_count < 4) return 0;
// Build vertex positions array (working copy)
std::vector<float> pos(mesh.positions); // x,y,z flat
std::vector<float> nrm(mesh.normals); // x,y,z flat
// Build triangle list (working copy — indices into pos)
std::vector<int> tris(mesh.indices);
// Per-vertex: which triangles reference it
std::vector<std::unordered_set<int>> vert_tris(vert_count);
for (int t = 0; t < tri_count; t++) {
vert_tris[tris[t*3+0]].insert(t);
vert_tris[tris[t*3+1]].insert(t);
vert_tris[tris[t*3+2]].insert(t);
}
// Compute per-vertex QEM from incident face planes
std::vector<QEM> vertex_qem(vert_count);
for (int t = 0; t < tri_count; t++) {
int i0 = tris[t*3+0], i1 = tris[t*3+1], i2 = tris[t*3+2];
float ax = pos[i1*3+0] - pos[i0*3+0];
float ay = pos[i1*3+1] - pos[i0*3+1];
float az = pos[i1*3+2] - pos[i0*3+2];
float bx = pos[i2*3+0] - pos[i0*3+0];
float by = pos[i2*3+1] - pos[i0*3+1];
float bz = pos[i2*3+2] - pos[i0*3+2];
// Cross product = face normal (unnormalized)
float nx = ay*bz - az*by;
float ny = az*bx - ax*bz;
float nz = ax*by - ay*bx;
float len = std::sqrt(nx*nx + ny*ny + nz*nz);
if (len < 1e-12f) continue;
nx /= len; ny /= len; nz /= len;
float d = -(nx*pos[i0*3+0] + ny*pos[i0*3+1] + nz*pos[i0*3+2]);
QEM q(nx, ny, nz, d);
vertex_qem[i0] += q;
vertex_qem[i1] += q;
vertex_qem[i2] += q;
}
// Remap table: vertex i -> actual vertex (follows collapse chains)
std::vector<int> remap(vert_count);
for (int i = 0; i < vert_count; i++) remap[i] = i;
auto resolve = [&](int v) -> int {
while (remap[v] != v) v = remap[v];
return v;
};
// Build edge set
std::unordered_set<Edge, EdgeHash> edge_set;
for (int t = 0; t < tri_count; t++) {
int v[3] = { tris[t*3+0], tris[t*3+1], tris[t*3+2] };
for (int e = 0; e < 3; e++) {
int a = v[e], b = v[(e+1)%3];
if (a > b) std::swap(a, b);
edge_set.insert({a, b});
}
}
std::vector<Edge> edges(edge_set.begin(), edge_set.end());
edge_set.clear();
// Alive triangle tracking
std::vector<bool> tri_alive(tri_count, true);
int alive_tris = tri_count;
// RNG for random edge selection
std::mt19937 rng(42);
int collapses = 0;
while (alive_tris > target_tri_count && !edges.empty()) {
// Pick K random candidate edges
int best_idx = -1;
float best_error = max_error;
float best_pos[3] = {};
int k = std::min(k_candidates, static_cast<int>(edges.size()));
for (int c = 0; c < k; c++) {
int idx = static_cast<int>(rng() % edges.size());
Edge& e = edges[idx];
int va = resolve(e.v0);
int vb = resolve(e.v1);
// Skip degenerate edges (already collapsed)
if (va == vb) continue;
// Evaluate collapse: place new vertex at midpoint
float mx = (pos[va*3+0] + pos[vb*3+0]) * 0.5f;
float my = (pos[va*3+1] + pos[vb*3+1]) * 0.5f;
float mz = (pos[va*3+2] + pos[vb*3+2]) * 0.5f;
QEM combined = vertex_qem[va] + vertex_qem[vb];
float err = std::fabs(combined.evaluate(mx, my, mz));
if (err < best_error) {
best_error = err;
best_idx = idx;
best_pos[0] = mx;
best_pos[1] = my;
best_pos[2] = mz;
}
}
if (best_idx < 0) break; // no valid candidate below max_error
// Collapse the best edge: merge vb into va
Edge& best_edge = edges[best_idx];
int va = resolve(best_edge.v0);
int vb = resolve(best_edge.v1);
if (va == vb) {
// Remove degenerate edge and continue
edges[best_idx] = edges.back();
edges.pop_back();
continue;
}
// Move va to optimal position
pos[va*3+0] = best_pos[0];
pos[va*3+1] = best_pos[1];
pos[va*3+2] = best_pos[2];
// Average normals
nrm[va*3+0] = (nrm[va*3+0] + nrm[vb*3+0]) * 0.5f;
nrm[va*3+1] = (nrm[va*3+1] + nrm[vb*3+1]) * 0.5f;
nrm[va*3+2] = (nrm[va*3+2] + nrm[vb*3+2]) * 0.5f;
float nl = std::sqrt(nrm[va*3+0]*nrm[va*3+0] + nrm[va*3+1]*nrm[va*3+1] + nrm[va*3+2]*nrm[va*3+2]);
if (nl > 1e-12f) {
nrm[va*3+0] /= nl; nrm[va*3+1] /= nl; nrm[va*3+2] /= nl;
}
// Merge QEMs
vertex_qem[va] += vertex_qem[vb];
// Remap vb -> va
remap[vb] = va;
// Transfer vb's triangles to va
for (int t : vert_tris[vb]) {
vert_tris[va].insert(t);
}
vert_tris[vb].clear();
// Update triangle indices and remove degenerate triangles
std::vector<int> to_remove;
for (int t : vert_tris[va]) {
if (!tri_alive[t]) continue;
// Remap all vertices in this triangle
for (int j = 0; j < 3; j++) {
tris[t*3+j] = resolve(tris[t*3+j]);
}
// Check for degenerate (two or more vertices the same)
int a = tris[t*3+0], b = tris[t*3+1], c = tris[t*3+2];
if (a == b || b == c || a == c) {
tri_alive[t] = false;
alive_tris--;
to_remove.push_back(t);
}
}
for (int t : to_remove) {
vert_tris[va].erase(t);
}
// Remove the collapsed edge
edges[best_idx] = edges.back();
edges.pop_back();
collapses++;
}
// Rebuild mesh from surviving triangles
// Compact vertices: only include those referenced by alive triangles
std::unordered_map<int, int> vert_map;
DCMeshData result;
for (int t = 0; t < tri_count; t++) {
if (!tri_alive[t]) continue;
int v[3];
for (int j = 0; j < 3; j++) {
int orig = resolve(tris[t*3+j]);
auto it = vert_map.find(orig);
if (it == vert_map.end()) {
int new_idx = static_cast<int>(result.positions.size()) / 3;
vert_map[orig] = new_idx;
result.positions.push_back(pos[orig*3+0]);
result.positions.push_back(pos[orig*3+1]);
result.positions.push_back(pos[orig*3+2]);
result.normals.push_back(nrm[orig*3+0]);
result.normals.push_back(nrm[orig*3+1]);
result.normals.push_back(nrm[orig*3+2]);
v[j] = new_idx;
} else {
v[j] = it->second;
}
}
result.indices.push_back(v[0]);
result.indices.push_back(v[1]);
result.indices.push_back(v[2]);
}
mesh = std::move(result);
return collapses;
}
// ---------------------------------------------------------------------------
// meshoptimizer-based simplification
// ---------------------------------------------------------------------------
#include <meshoptimizer.h>
void dc_simplify_meshopt(DCMeshData& mesh, int target_tri_count, bool lock_border) {
int vert_count = static_cast<int>(mesh.positions.size() / 3);
int idx_count = static_cast<int>(mesh.indices.size());
if (vert_count == 0 || idx_count == 0) return;
// meshoptimizer uses unsigned int indices
std::vector<unsigned int> indices(idx_count);
for (int i = 0; i < idx_count; i++) {
indices[i] = static_cast<unsigned int>(mesh.indices[i]);
}
size_t target_idx_count = static_cast<size_t>(target_tri_count) * 3;
float target_error = 1e-2f;
unsigned int options = lock_border ? meshopt_SimplifyLockBorder : 0;
std::vector<unsigned int> result_indices(idx_count);
size_t result_count = meshopt_simplify(
result_indices.data(),
indices.data(), idx_count,
mesh.positions.data(), vert_count, sizeof(float) * 3,
target_idx_count, target_error, options, nullptr);
result_indices.resize(result_count);
// Remap to compact vertex set
std::vector<unsigned int> remap(vert_count, ~0u);
unsigned int new_vert_count = 0;
for (unsigned int idx : result_indices) {
if (remap[idx] == ~0u) {
remap[idx] = new_vert_count++;
}
}
std::vector<float> new_positions(new_vert_count * 3);
std::vector<float> new_normals(new_vert_count * 3);
bool has_mat = ((int)mesh.materials.size() == vert_count); // route ② texturing
bool has_proj = ((int)mesh.proj.size() == vert_count); // Stage-2 projection candidate 0
bool has_proj2 = ((int)mesh.proj2.size() == vert_count); // Stage-2 projection candidate 1
bool has_curv = ((int)mesh.curv.size() == vert_count); // curvature attribute
std::vector<float> new_materials(has_mat ? new_vert_count : 0);
std::vector<float> new_proj(has_proj ? new_vert_count : 0);
std::vector<float> new_proj2(has_proj2 ? new_vert_count : 0);
std::vector<float> new_curv(has_curv ? new_vert_count : 0);
for (int v = 0; v < vert_count; v++) {
if (remap[v] != ~0u) {
unsigned int nv = remap[v];
new_positions[nv * 3 + 0] = mesh.positions[v * 3 + 0];
new_positions[nv * 3 + 1] = mesh.positions[v * 3 + 1];
new_positions[nv * 3 + 2] = mesh.positions[v * 3 + 2];
new_normals[nv * 3 + 0] = mesh.normals[v * 3 + 0];
new_normals[nv * 3 + 1] = mesh.normals[v * 3 + 1];
new_normals[nv * 3 + 2] = mesh.normals[v * 3 + 2];
if (has_mat) new_materials[nv] = mesh.materials[v];
if (has_proj) new_proj[nv] = mesh.proj[v];
if (has_proj2) new_proj2[nv] = mesh.proj2[v];
if (has_curv) new_curv[nv] = mesh.curv[v];
}
}
mesh.positions = std::move(new_positions);
mesh.normals = std::move(new_normals);
if (has_mat) mesh.materials = std::move(new_materials);
if (has_proj) mesh.proj = std::move(new_proj);
if (has_proj2) mesh.proj2 = std::move(new_proj2);
if (has_curv) mesh.curv = std::move(new_curv);
mesh.indices.resize(result_count);
for (size_t i = 0; i < result_count; i++) {
mesh.indices[i] = static_cast<int>(remap[result_indices[i]]);
}
}
void dc_optimize_vertex_cache(DCMeshData& mesh) {
int idx_count = static_cast<int>(mesh.indices.size());
int vert_count = static_cast<int>(mesh.positions.size() / 3);
if (idx_count == 0 || vert_count == 0) return;
std::vector<unsigned int> indices(idx_count);
for (int i = 0; i < idx_count; i++) {
indices[i] = static_cast<unsigned int>(mesh.indices[i]);
}
meshopt_optimizeVertexCache(indices.data(), indices.data(), idx_count, vert_count);
for (int i = 0; i < idx_count; i++) {
mesh.indices[i] = static_cast<int>(indices[i]);
}
}