From f3ffcc648a33418d75c5788752e9f47f7cfd0a75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=91=D1=80=D0=B0=D0=BD=D0=B8=D0=BC=D0=B8=D1=80=20=D0=9A?= =?UTF-8?q?=D0=B0=D1=80=D0=B0=D1=9F=D0=B8=D1=9B?= Date: Sun, 30 Aug 2026 20:49:06 -0700 Subject: [PATCH] Updated meshoptimizer. --- 3rdparty/meshoptimizer/LICENSE.md | 2 +- 3rdparty/meshoptimizer/src/allocator.cpp | 4 + 3rdparty/meshoptimizer/src/clusterizer.cpp | 449 +------ 3rdparty/meshoptimizer/src/indexanalyzer.cpp | 1 - 3rdparty/meshoptimizer/src/indexcodec.cpp | 73 +- 3rdparty/meshoptimizer/src/indexgenerator.cpp | 127 +- 3rdparty/meshoptimizer/src/meshletcodec.cpp | 1051 +++++++++++++++++ 3rdparty/meshoptimizer/src/meshletutils.cpp | 619 ++++++++++ 3rdparty/meshoptimizer/src/meshoptimizer.h | 349 +++++- 3rdparty/meshoptimizer/src/opacitymap.cpp | 557 +++++++++ 3rdparty/meshoptimizer/src/partition.cpp | 54 +- 3rdparty/meshoptimizer/src/quantization.cpp | 51 + 3rdparty/meshoptimizer/src/rasterizer.cpp | 2 +- 3rdparty/meshoptimizer/src/remesher.cpp | 751 ++++++++++++ 3rdparty/meshoptimizer/src/simplifier.cpp | 120 +- 3rdparty/meshoptimizer/src/tangentspace.cpp | 730 ++++++++++++ 3rdparty/meshoptimizer/src/vertexcodec.cpp | 376 +++--- 3rdparty/meshoptimizer/src/vertexfilter.cpp | 62 +- .../meshoptimizer/src/vfetchoptimizer.cpp | 3 - 19 files changed, 4605 insertions(+), 776 deletions(-) create mode 100644 3rdparty/meshoptimizer/src/meshletcodec.cpp create mode 100644 3rdparty/meshoptimizer/src/meshletutils.cpp create mode 100644 3rdparty/meshoptimizer/src/opacitymap.cpp create mode 100644 3rdparty/meshoptimizer/src/remesher.cpp create mode 100644 3rdparty/meshoptimizer/src/tangentspace.cpp diff --git a/3rdparty/meshoptimizer/LICENSE.md b/3rdparty/meshoptimizer/LICENSE.md index a5c3b1ccc..468827f4b 100644 --- a/3rdparty/meshoptimizer/LICENSE.md +++ b/3rdparty/meshoptimizer/LICENSE.md @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2016-2025 Arseny Kapoulkine +Copyright (c) 2016-2026 Arseny Kapoulkine Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/3rdparty/meshoptimizer/src/allocator.cpp b/3rdparty/meshoptimizer/src/allocator.cpp index 6b6083da2..058fbf312 100644 --- a/3rdparty/meshoptimizer/src/allocator.cpp +++ b/3rdparty/meshoptimizer/src/allocator.cpp @@ -1,6 +1,8 @@ // This file is part of meshoptimizer library; see meshoptimizer.h for version/license details #include "meshoptimizer.h" +#include + #ifdef MESHOPTIMIZER_ALLOC_EXPORT meshopt_Allocator::Storage& meshopt_Allocator::storage() { @@ -11,6 +13,8 @@ meshopt_Allocator::Storage& meshopt_Allocator::storage() void meshopt_setAllocator(void* (MESHOPTIMIZER_ALLOC_CALLCONV* allocate)(size_t), void (MESHOPTIMIZER_ALLOC_CALLCONV* deallocate)(void*)) { + assert(allocate && deallocate); + meshopt_Allocator::Storage& s = meshopt_Allocator::storage(); s.allocate = allocate; s.deallocate = deallocate; diff --git a/3rdparty/meshoptimizer/src/clusterizer.cpp b/3rdparty/meshoptimizer/src/clusterizer.cpp index 73cc0ab53..0b630b2cb 100644 --- a/3rdparty/meshoptimizer/src/clusterizer.cpp +++ b/3rdparty/meshoptimizer/src/clusterizer.cpp @@ -8,10 +8,10 @@ // The block below auto-detects SIMD ISA that can be used on the target platform #ifndef MESHOPTIMIZER_NO_SIMD -#if defined(__SSE2__) || (defined(_MSC_VER) && defined(_M_X64)) +#if defined(__SSE2__) || (defined(_MSC_VER) && defined(_M_X64) && !defined(_M_ARM64EC)) #define SIMD_SSE #include -#elif defined(__aarch64__) || (defined(_MSC_VER) && defined(_M_ARM64) && _MSC_VER >= 1922) +#elif defined(__aarch64__) || (defined(_MSC_VER) && (defined(_M_ARM64) || defined(_M_ARM64EC)) && _MSC_VER >= 1922) #define SIMD_NEON #include #endif @@ -19,19 +19,10 @@ // This work is based on: // Graham Wihlidal. Optimizing the Graphics Pipeline with Compute. 2016 -// Matthaeus Chajdas. GeometryFX 1.2 - Cluster Culling. 2016 -// Jack Ritter. An Efficient Bounding Sphere. 1990 -// Thomas Larsson. Fast and Tight Fitting Bounding Spheres. 2008 // Ingo Wald, Vlastimil Havran. On building fast kd-Trees for Ray Tracing, and on doing that in O(N log N). 2006 namespace meshopt { -// This must be <= 256 since meshlet indices are stored as bytes -const size_t kMeshletMaxVertices = 256; - -// A reasonable limit is around 2*max_vertices or less -const size_t kMeshletMaxTriangles = 512; - // We keep a limited number of seed triangles and add a few triangles per finished meshlet const size_t kMeshletMaxSeeds = 256; const size_t kMeshletAddSeeds = 4; @@ -173,116 +164,6 @@ static void clearUsed(short* used, size_t vertex_count, const unsigned int* indi } } -static void computeBoundingSphere(float result[4], const float* points, size_t count, size_t points_stride, const float* radii, size_t radii_stride, size_t axis_count) -{ - static const float kAxes[7][3] = { - // X, Y, Z - {1, 0, 0}, - {0, 1, 0}, - {0, 0, 1}, - - // XYZ, -XYZ, X-YZ, XY-Z; normalized to unit length - {0.57735026f, 0.57735026f, 0.57735026f}, - {-0.57735026f, 0.57735026f, 0.57735026f}, - {0.57735026f, -0.57735026f, 0.57735026f}, - {0.57735026f, 0.57735026f, -0.57735026f}, - }; - - assert(count > 0); - assert(axis_count <= sizeof(kAxes) / sizeof(kAxes[0])); - - size_t points_stride_float = points_stride / sizeof(float); - size_t radii_stride_float = radii_stride / sizeof(float); - - // find extremum points along all axes; for each axis we get a pair of points with min/max coordinates - size_t pmin[7], pmax[7]; - float tmin[7], tmax[7]; - - for (size_t axis = 0; axis < axis_count; ++axis) - { - pmin[axis] = pmax[axis] = 0; - tmin[axis] = FLT_MAX; - tmax[axis] = -FLT_MAX; - } - - for (size_t i = 0; i < count; ++i) - { - const float* p = points + i * points_stride_float; - float r = radii[i * radii_stride_float]; - - for (size_t axis = 0; axis < axis_count; ++axis) - { - const float* ax = kAxes[axis]; - - float tp = ax[0] * p[0] + ax[1] * p[1] + ax[2] * p[2]; - float tpmin = tp - r, tpmax = tp + r; - - pmin[axis] = (tpmin < tmin[axis]) ? i : pmin[axis]; - pmax[axis] = (tpmax > tmax[axis]) ? i : pmax[axis]; - tmin[axis] = (tpmin < tmin[axis]) ? tpmin : tmin[axis]; - tmax[axis] = (tpmax > tmax[axis]) ? tpmax : tmax[axis]; - } - } - - // find the pair of points with largest distance - size_t paxis = 0; - float paxisdr = 0; - - for (size_t axis = 0; axis < axis_count; ++axis) - { - const float* p1 = points + pmin[axis] * points_stride_float; - const float* p2 = points + pmax[axis] * points_stride_float; - float r1 = radii[pmin[axis] * radii_stride_float]; - float r2 = radii[pmax[axis] * radii_stride_float]; - - float d2 = (p2[0] - p1[0]) * (p2[0] - p1[0]) + (p2[1] - p1[1]) * (p2[1] - p1[1]) + (p2[2] - p1[2]) * (p2[2] - p1[2]); - float dr = sqrtf(d2) + r1 + r2; - - if (dr > paxisdr) - { - paxisdr = dr; - paxis = axis; - } - } - - // use the longest segment as the initial sphere diameter - const float* p1 = points + pmin[paxis] * points_stride_float; - const float* p2 = points + pmax[paxis] * points_stride_float; - float r1 = radii[pmin[paxis] * radii_stride_float]; - float r2 = radii[pmax[paxis] * radii_stride_float]; - - float paxisd = sqrtf((p2[0] - p1[0]) * (p2[0] - p1[0]) + (p2[1] - p1[1]) * (p2[1] - p1[1]) + (p2[2] - p1[2]) * (p2[2] - p1[2])); - float paxisk = paxisd > 0 ? (paxisd + r2 - r1) / (2 * paxisd) : 0.f; - - float center[3] = {p1[0] + (p2[0] - p1[0]) * paxisk, p1[1] + (p2[1] - p1[1]) * paxisk, p1[2] + (p2[2] - p1[2]) * paxisk}; - float radius = paxisdr / 2; - - // iteratively adjust the sphere up until all points fit - for (size_t i = 0; i < count; ++i) - { - const float* p = points + i * points_stride_float; - float r = radii[i * radii_stride_float]; - - float d2 = (p[0] - center[0]) * (p[0] - center[0]) + (p[1] - center[1]) * (p[1] - center[1]) + (p[2] - center[2]) * (p[2] - center[2]); - float d = sqrtf(d2); - - if (d + r > radius) - { - float k = d > 0 ? (d + r - radius) / (2 * d) : 0.f; - - center[0] += k * (p[0] - center[0]); - center[1] += k * (p[1] - center[1]); - center[2] += k * (p[2] - center[2]); - radius = (radius + d + r) / 2; - } - } - - result[0] = center[0]; - result[1] = center[1]; - result[2] = center[2]; - result[3] = radius; -} - struct Cone { float px, py, pz; @@ -365,6 +246,14 @@ static bool appendMeshlet(meshopt_Meshlet& meshlet, unsigned int a, unsigned int short& bv = used[b]; short& cv = used[c]; +#ifdef MESHOPTIMIZER_CLUSTERIZER_INDEXLIMIT + // to ensure that after adding up to 2 new vertices, the existing vertex is within span limit (32), we need to reset vertex indices that are out of range + const int span_limit = 32 - 2; + av = (av + span_limit < int(meshlet.vertex_count)) ? -1 : av; + bv = (bv + span_limit < int(meshlet.vertex_count)) ? -1 : bv; + cv = (cv + span_limit < int(meshlet.vertex_count)) ? -1 : cv; +#endif + bool result = false; int used_extra = (av < 0) + (bv < 0) + (cv < 0); @@ -813,7 +702,7 @@ static float boxMerge(BVHBoxT& box, const BVHBox& other) } float sx = box.max[0] - box.min[0], sy = box.max[1] - box.min[1], sz = box.max[2] - box.min[2]; - return sx * sy + sx * sz + sy * sz; + return sx * sy + sy * sz + sz * sx; } #endif @@ -961,7 +850,7 @@ static void bvhPackTail(unsigned char* boundary, const unsigned int* order, size static bool bvhDivisible(size_t count, size_t min, size_t max) { - // count is representable as a sum of values in [min..max] if if it in range of [k*min..k*min+k*(max-min)] + // count is representable as a sum of values in [min..max] if it is in range of [k*min..k*min+k*(max-min)] // equivalent to ceil(count / max) <= floor(count / min), but the form below allows using idiv (see nv_cluster_builder) // we avoid expensive integer divisions in the common case where min is <= max/2 return min * 2 <= max ? count >= min : count % min <= (count / min) * (max - min); @@ -1130,11 +1019,8 @@ size_t meshopt_buildMeshletsBound(size_t index_count, size_t max_vertices, size_ using namespace meshopt; assert(index_count % 3 == 0); - assert(max_vertices >= 3 && max_vertices <= kMeshletMaxVertices); - assert(max_triangles >= 1 && max_triangles <= kMeshletMaxTriangles); - - (void)kMeshletMaxVertices; - (void)kMeshletMaxTriangles; + assert(max_vertices >= 3 && max_vertices <= 256); + assert(max_triangles >= 1 && max_triangles <= 512); // meshlet construction is limited by max vertices and max triangles per meshlet // the worst case is that the input is an unindexed stream since this equally stresses both limits @@ -1154,8 +1040,8 @@ size_t meshopt_buildMeshletsFlex(meshopt_Meshlet* meshlets, unsigned int* meshle assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); assert(vertex_positions_stride % sizeof(float) == 0); - assert(max_vertices >= 3 && max_vertices <= kMeshletMaxVertices); - assert(min_triangles >= 1 && min_triangles <= max_triangles && max_triangles <= kMeshletMaxTriangles); + assert(max_vertices >= 3 && max_vertices <= 256); + assert(min_triangles >= 1 && min_triangles <= max_triangles && max_triangles <= 512); assert(cone_weight >= 0 && cone_weight <= 1); assert(split_factor >= 0); @@ -1348,8 +1234,8 @@ size_t meshopt_buildMeshletsScan(meshopt_Meshlet* meshlets, unsigned int* meshle assert(index_count % 3 == 0); - assert(max_vertices >= 3 && max_vertices <= kMeshletMaxVertices); - assert(max_triangles >= 1 && max_triangles <= kMeshletMaxTriangles); + assert(max_vertices >= 3 && max_vertices <= 256); + assert(max_triangles >= 1 && max_triangles <= 512); meshopt_Allocator allocator; @@ -1385,8 +1271,10 @@ size_t meshopt_buildMeshletsSpatial(struct meshopt_Meshlet* meshlets, unsigned i assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); assert(vertex_positions_stride % sizeof(float) == 0); - assert(max_vertices >= 3 && max_vertices <= kMeshletMaxVertices); - assert(min_triangles >= 1 && min_triangles <= max_triangles && max_triangles <= kMeshletMaxTriangles); + assert(max_vertices >= 3 && max_vertices <= 256); + assert(min_triangles >= 1 && min_triangles <= max_triangles && max_triangles <= 512); + + assert(fill_weight >= 0); if (index_count == 0) return 0; @@ -1476,298 +1364,5 @@ size_t meshopt_buildMeshletsSpatial(struct meshopt_Meshlet* meshlets, unsigned i return meshlet_offset; } -meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride) -{ - using namespace meshopt; - - assert(index_count % 3 == 0); - assert(index_count / 3 <= kMeshletMaxTriangles); - assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); - assert(vertex_positions_stride % sizeof(float) == 0); - - (void)vertex_count; - - size_t vertex_stride_float = vertex_positions_stride / sizeof(float); - - // compute triangle normals and gather triangle corners - float normals[kMeshletMaxTriangles][3]; - float corners[kMeshletMaxTriangles][3][3]; - size_t triangles = 0; - - for (size_t i = 0; i < index_count; i += 3) - { - unsigned int a = indices[i + 0], b = indices[i + 1], c = indices[i + 2]; - assert(a < vertex_count && b < vertex_count && c < vertex_count); - - const float* p0 = vertex_positions + vertex_stride_float * a; - const float* p1 = vertex_positions + vertex_stride_float * b; - const float* p2 = vertex_positions + vertex_stride_float * c; - - float p10[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; - float p20[3] = {p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; - - float normalx = p10[1] * p20[2] - p10[2] * p20[1]; - float normaly = p10[2] * p20[0] - p10[0] * p20[2]; - float normalz = p10[0] * p20[1] - p10[1] * p20[0]; - - float area = sqrtf(normalx * normalx + normaly * normaly + normalz * normalz); - - // no need to include degenerate triangles - they will be invisible anyway - if (area == 0.f) - continue; - - // record triangle normals & corners for future use; normal and corner 0 define a plane equation - normals[triangles][0] = normalx / area; - normals[triangles][1] = normaly / area; - normals[triangles][2] = normalz / area; - memcpy(corners[triangles][0], p0, 3 * sizeof(float)); - memcpy(corners[triangles][1], p1, 3 * sizeof(float)); - memcpy(corners[triangles][2], p2, 3 * sizeof(float)); - triangles++; - } - - meshopt_Bounds bounds = {}; - - // degenerate cluster, no valid triangles => trivial reject (cone data is 0) - if (triangles == 0) - return bounds; - - const float rzero = 0.f; - - // compute cluster bounding sphere; we'll use the center to determine normal cone apex as well - float psphere[4] = {}; - computeBoundingSphere(psphere, corners[0][0], triangles * 3, sizeof(float) * 3, &rzero, 0, 7); - - float center[3] = {psphere[0], psphere[1], psphere[2]}; - - // treating triangle normals as points, find the bounding sphere - the sphere center determines the optimal cone axis - float nsphere[4] = {}; - computeBoundingSphere(nsphere, normals[0], triangles, sizeof(float) * 3, &rzero, 0, 3); - - float axis[3] = {nsphere[0], nsphere[1], nsphere[2]}; - float axislength = sqrtf(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]); - float invaxislength = axislength == 0.f ? 0.f : 1.f / axislength; - - axis[0] *= invaxislength; - axis[1] *= invaxislength; - axis[2] *= invaxislength; - - // compute a tight cone around all normals, mindp = cos(angle/2) - float mindp = 1.f; - - for (size_t i = 0; i < triangles; ++i) - { - float dp = normals[i][0] * axis[0] + normals[i][1] * axis[1] + normals[i][2] * axis[2]; - - mindp = (dp < mindp) ? dp : mindp; - } - - // fill bounding sphere info; note that below we can return bounds without cone information for degenerate cones - bounds.center[0] = center[0]; - bounds.center[1] = center[1]; - bounds.center[2] = center[2]; - bounds.radius = psphere[3]; - - // degenerate cluster, normal cone is larger than a hemisphere => trivial accept - // note that if mindp is positive but close to 0, the triangle intersection code below gets less stable - // we arbitrarily decide that if a normal cone is ~168 degrees wide or more, the cone isn't useful - if (mindp <= 0.1f) - { - bounds.cone_cutoff = 1; - bounds.cone_cutoff_s8 = 127; - return bounds; - } - - float maxt = 0; - - // we need to find the point on center-t*axis ray that lies in negative half-space of all triangles - for (size_t i = 0; i < triangles; ++i) - { - // dot(center-t*axis-corner, trinormal) = 0 - // dot(center-corner, trinormal) - t * dot(axis, trinormal) = 0 - float cx = center[0] - corners[i][0][0]; - float cy = center[1] - corners[i][0][1]; - float cz = center[2] - corners[i][0][2]; - - float dc = cx * normals[i][0] + cy * normals[i][1] + cz * normals[i][2]; - float dn = axis[0] * normals[i][0] + axis[1] * normals[i][1] + axis[2] * normals[i][2]; - - // dn should be larger than mindp cutoff above - assert(dn > 0.f); - float t = dc / dn; - - maxt = (t > maxt) ? t : maxt; - } - - // cone apex should be in the negative half-space of all cluster triangles by construction - bounds.cone_apex[0] = center[0] - axis[0] * maxt; - bounds.cone_apex[1] = center[1] - axis[1] * maxt; - bounds.cone_apex[2] = center[2] - axis[2] * maxt; - - // note: this axis is the axis of the normal cone, but our test for perspective camera effectively negates the axis - bounds.cone_axis[0] = axis[0]; - bounds.cone_axis[1] = axis[1]; - bounds.cone_axis[2] = axis[2]; - - // cos(a) for normal cone is mindp; we need to add 90 degrees on both sides and invert the cone - // which gives us -cos(a+90) = -(-sin(a)) = sin(a) = sqrt(1 - cos^2(a)) - bounds.cone_cutoff = sqrtf(1 - mindp * mindp); - - // quantize axis & cutoff to 8-bit SNORM format - bounds.cone_axis_s8[0] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[0], 8)); - bounds.cone_axis_s8[1] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[1], 8)); - bounds.cone_axis_s8[2] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[2], 8)); - - // for the 8-bit test to be conservative, we need to adjust the cutoff by measuring the max. error - float cone_axis_s8_e0 = fabsf(bounds.cone_axis_s8[0] / 127.f - bounds.cone_axis[0]); - float cone_axis_s8_e1 = fabsf(bounds.cone_axis_s8[1] / 127.f - bounds.cone_axis[1]); - float cone_axis_s8_e2 = fabsf(bounds.cone_axis_s8[2] / 127.f - bounds.cone_axis[2]); - - // note that we need to round this up instead of rounding to nearest, hence +1 - int cone_cutoff_s8 = int(127 * (bounds.cone_cutoff + cone_axis_s8_e0 + cone_axis_s8_e1 + cone_axis_s8_e2) + 1); - - bounds.cone_cutoff_s8 = (cone_cutoff_s8 > 127) ? 127 : (signed char)(cone_cutoff_s8); - - return bounds; -} - -meshopt_Bounds meshopt_computeMeshletBounds(const unsigned int* meshlet_vertices, const unsigned char* meshlet_triangles, size_t triangle_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride) -{ - using namespace meshopt; - - assert(triangle_count <= kMeshletMaxTriangles); - assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); - assert(vertex_positions_stride % sizeof(float) == 0); - - unsigned int indices[kMeshletMaxTriangles * 3]; - - for (size_t i = 0; i < triangle_count * 3; ++i) - { - unsigned int index = meshlet_vertices[meshlet_triangles[i]]; - assert(index < vertex_count); - - indices[i] = index; - } - - return meshopt_computeClusterBounds(indices, triangle_count * 3, vertex_positions, vertex_count, vertex_positions_stride); -} - -meshopt_Bounds meshopt_computeSphereBounds(const float* positions, size_t count, size_t positions_stride, const float* radii, size_t radii_stride) -{ - using namespace meshopt; - - assert(positions_stride >= 12 && positions_stride <= 256); - assert(positions_stride % sizeof(float) == 0); - assert((radii_stride >= 4 && radii_stride <= 256) || radii == NULL); - assert(radii_stride % sizeof(float) == 0); - - meshopt_Bounds bounds = {}; - - if (count == 0) - return bounds; - - const float rzero = 0.f; - - float psphere[4] = {}; - computeBoundingSphere(psphere, positions, count, positions_stride, radii ? radii : &rzero, radii ? radii_stride : 0, 7); - - bounds.center[0] = psphere[0]; - bounds.center[1] = psphere[1]; - bounds.center[2] = psphere[2]; - bounds.radius = psphere[3]; - - return bounds; -} - -void meshopt_optimizeMeshlet(unsigned int* meshlet_vertices, unsigned char* meshlet_triangles, size_t triangle_count, size_t vertex_count) -{ - using namespace meshopt; - - assert(triangle_count <= kMeshletMaxTriangles); - assert(vertex_count <= kMeshletMaxVertices); - - unsigned char* indices = meshlet_triangles; - unsigned int* vertices = meshlet_vertices; - - // cache tracks vertex timestamps (corresponding to triangle index! all 3 vertices are added at the same time and never removed) - unsigned char cache[kMeshletMaxVertices]; - memset(cache, 0, vertex_count); - - // note that we start from a value that means all vertices aren't in cache - unsigned char cache_last = 128; - const unsigned char cache_cutoff = 3; // 3 triangles = ~5..9 vertices depending on reuse - - for (size_t i = 0; i < triangle_count; ++i) - { - int next = -1; - int next_match = -1; - - for (size_t j = i; j < triangle_count; ++j) - { - unsigned char a = indices[j * 3 + 0], b = indices[j * 3 + 1], c = indices[j * 3 + 2]; - assert(a < vertex_count && b < vertex_count && c < vertex_count); - - // score each triangle by how many vertices are in cache - // note: the distance is computed using unsigned 8-bit values, so cache timestamp overflow is handled gracefully - int aok = (unsigned char)(cache_last - cache[a]) < cache_cutoff; - int bok = (unsigned char)(cache_last - cache[b]) < cache_cutoff; - int cok = (unsigned char)(cache_last - cache[c]) < cache_cutoff; - - if (aok + bok + cok > next_match) - { - next = (int)j; - next_match = aok + bok + cok; - - // note that we could end up with all 3 vertices in the cache, but 2 is enough for ~strip traversal - if (next_match >= 2) - break; - } - } - - assert(next >= 0); - - unsigned char a = indices[next * 3 + 0], b = indices[next * 3 + 1], c = indices[next * 3 + 2]; - - // shift triangles before the next one forward so that we always keep an ordered partition - // note: this could have swapped triangles [i] and [next] but that distorts the order and may skew the output sequence - memmove(indices + (i + 1) * 3, indices + i * 3, (next - i) * 3 * sizeof(unsigned char)); - - indices[i * 3 + 0] = a; - indices[i * 3 + 1] = b; - indices[i * 3 + 2] = c; - - // cache timestamp is the same between all vertices of each triangle to reduce overflow - cache_last++; - cache[a] = cache_last; - cache[b] = cache_last; - cache[c] = cache_last; - } - - // reorder meshlet vertices for access locality assuming index buffer is scanned sequentially - unsigned int order[kMeshletMaxVertices]; - - short remap[kMeshletMaxVertices]; - memset(remap, -1, vertex_count * sizeof(short)); - - size_t vertex_offset = 0; - - for (size_t i = 0; i < triangle_count * 3; ++i) - { - short& r = remap[indices[i]]; - - if (r < 0) - { - r = short(vertex_offset); - order[vertex_offset] = vertices[indices[i]]; - vertex_offset++; - } - - indices[i] = (unsigned char)r; - } - - assert(vertex_offset <= vertex_count); - memcpy(vertices, order, vertex_offset * sizeof(unsigned int)); -} - #undef SIMD_SSE #undef SIMD_NEON diff --git a/3rdparty/meshoptimizer/src/indexanalyzer.cpp b/3rdparty/meshoptimizer/src/indexanalyzer.cpp index 87ceeae66..e8e975c27 100644 --- a/3rdparty/meshoptimizer/src/indexanalyzer.cpp +++ b/3rdparty/meshoptimizer/src/indexanalyzer.cpp @@ -74,7 +74,6 @@ meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const unsigned int* ind meshopt_VertexFetchStatistics meshopt_analyzeVertexFetch(const unsigned int* indices, size_t index_count, size_t vertex_count, size_t vertex_size) { - assert(index_count % 3 == 0); assert(vertex_size > 0 && vertex_size <= 256); meshopt_Allocator allocator; diff --git a/3rdparty/meshoptimizer/src/indexcodec.cpp b/3rdparty/meshoptimizer/src/indexcodec.cpp index 7a8fd6867..8dead6343 100644 --- a/3rdparty/meshoptimizer/src/indexcodec.cpp +++ b/3rdparty/meshoptimizer/src/indexcodec.cpp @@ -19,15 +19,9 @@ const int kDecodeIndexVersion = 1; typedef unsigned int VertexFifo[16]; typedef unsigned int EdgeFifo[16][2]; -static const unsigned int kTriangleIndexOrder[3][3] = { - {0, 1, 2}, - {1, 2, 0}, - {2, 0, 1}, -}; - static const unsigned char kCodeAuxEncodingTable[16] = { 0x00, 0x76, 0x87, 0x56, 0x67, 0x78, 0xa9, 0x86, 0x65, 0x89, 0x68, 0x98, 0x01, 0x69, - 0, 0, // last two entries aren't used for encoding + 0, 0 // last two entries aren't used for encoding }; static int rotateTriangle(unsigned int a, unsigned int b, unsigned int c, unsigned int next) @@ -144,19 +138,25 @@ static int getCodeAuxIndex(unsigned char v, const unsigned char* table) return -1; } -static void writeTriangle(void* destination, size_t offset, size_t index_size, unsigned int a, unsigned int b, unsigned int c) +static void* writeTriangle(void* destination, size_t index_size, unsigned int a, unsigned int b, unsigned int c) { if (index_size == 2) { - static_cast(destination)[offset + 0] = (unsigned short)(a); - static_cast(destination)[offset + 1] = (unsigned short)(b); - static_cast(destination)[offset + 2] = (unsigned short)(c); + unsigned short* tri = static_cast(destination); + tri[0] = (unsigned short)(a); + tri[1] = (unsigned short)(b); + tri[2] = (unsigned short)(c); + + return tri + 3; } else { - static_cast(destination)[offset + 0] = a; - static_cast(destination)[offset + 1] = b; - static_cast(destination)[offset + 2] = c; + unsigned int* tri = static_cast(destination); + tri[0] = a; + tri[1] = b; + tri[2] = c; + + return tri + 3; } } @@ -194,6 +194,8 @@ size_t meshopt_encodeIndexBuffer(unsigned char* buffer, size_t buffer_size, cons int fecmax = version >= 1 ? 13 : 15; + static const int rotations[] = {0, 1, 2, 0, 1}; + // use static encoding table; it's possible to pack the result and then build an optimal table and repack // for now we keep it simple and use the table that has been generated based on symbol frequency on a training mesh set const unsigned char* codeaux_table = kCodeAuxEncodingTable; @@ -211,7 +213,7 @@ size_t meshopt_encodeIndexBuffer(unsigned char* buffer, size_t buffer_size, cons if (fer >= 0 && (fer >> 2) < 15) { // note: getEdgeFifo implicitly rotates triangles by matching a/b to existing edge - const unsigned int* order = kTriangleIndexOrder[fer & 3]; + const int* order = rotations + (fer & 3); unsigned int a = indices[i + order[0]], b = indices[i + order[1]], c = indices[i + order[2]]; @@ -247,7 +249,7 @@ size_t meshopt_encodeIndexBuffer(unsigned char* buffer, size_t buffer_size, cons else { int rotation = rotateTriangle(indices[i + 0], indices[i + 1], indices[i + 2], next); - const unsigned int* order = kTriangleIndexOrder[rotation]; + const int* order = rotations + rotation; unsigned int a = indices[i + order[0]], b = indices[i + order[1]], c = indices[i + order[2]]; @@ -413,19 +415,16 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde // since we store 16-byte codeaux table at the end, triangle data has to begin before data_safe_end const unsigned char* code = buffer + 1; - const unsigned char* data = code + index_count / 3; + const unsigned char* code_end = code + index_count / 3; + const unsigned char* data = code_end; + + // each triangle reads at most 16 bytes of data: 1b for codeaux and 5b for each free index const unsigned char* data_safe_end = buffer + buffer_size - 16; const unsigned char* codeaux_table = data_safe_end; - for (size_t i = 0; i < index_count; i += 3) + while (code < code_end) { - // make sure we have enough data to read for a triangle - // each triangle reads at most 16 bytes of data: 1b for codeaux and 5b for each free index - // after this we can be sure we can read without extra bounds checks - if (data > data_safe_end) - return -2; - unsigned char codetri = *code++; if (codetri < 0xf0) @@ -445,6 +444,14 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde { // fifo reads are wrapped around 16 entry buffer unsigned int cf = vertexfifo[(vertexfifooffset - 1 - fec) & 15]; + +#if (defined(__GNUC__) || defined(__clang__)) && (defined(__x86_64__) || defined(__aarch64__)) + // on clang, x86-cmov-conversion pass emits a branch since cf is a load from memory; asm barrier defeats this + // this can be fixed with __builtin_unpredictable, but gcc doesn't support it and has a similar problem due to if-conversion + // additionally, gcc can fuse a/b into one 64-bit load but use 32-bit edgefifo[] stores, which breaks load store forwarding + __asm__("" : "+r"(a) : "r"(cf)); +#endif + c = (fec == 0) ? next : cf; int fec0 = fec == 0; @@ -455,9 +462,13 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde } else { - // fec - (fec ^ 3) decodes 13, 14 into -1, 1 + // make sure we have enough data to read for a triangle; this check covers worst case advance + if (data > data_safe_end) + return -2; + + // fec * 2 - 27 decodes 13, 14 into -1, 1 // note that we need to update the last index since free indices are delta-encoded - last = c = (fec != 15) ? last + (fec - (fec ^ 3)) : decodeIndex(data, last); + last = c = (fec != 15) ? last + (fec * 2 - 27) : decodeIndex(data, last); // push vertex/edge fifo must match the encoding step *exactly* otherwise the data will not be decoded correctly pushVertexFifo(vertexfifo, c, vertexfifooffset); @@ -468,7 +479,7 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde pushEdgeFifo(edgefifo, a, c, edgefifooffset); // output triangle - writeTriangle(destination, i, index_size, a, b, c); + destination = writeTriangle(destination, index_size, a, b, c); } else { @@ -498,7 +509,7 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde next += fec0; // output triangle - writeTriangle(destination, i, index_size, a, b, c); + destination = writeTriangle(destination, index_size, a, b, c); // push vertex/edge fifo must match the encoding step *exactly* otherwise the data will not be decoded correctly pushVertexFifo(vertexfifo, a, vertexfifooffset); @@ -511,6 +522,10 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde } else { + // make sure we have enough data to read for a triangle; this check covers worst case advance + if (data > data_safe_end) + return -2; + // slow path: read a full byte for codeaux instead of using a table lookup unsigned char codeaux = *data++; @@ -539,7 +554,7 @@ int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t inde last = c = decodeIndex(data, last); // output triangle - writeTriangle(destination, i, index_size, a, b, c); + destination = writeTriangle(destination, index_size, a, b, c); // push vertex/edge fifo must match the encoding step *exactly* otherwise the data will not be decoded correctly pushVertexFifo(vertexfifo, a, vertexfifooffset); diff --git a/3rdparty/meshoptimizer/src/indexgenerator.cpp b/3rdparty/meshoptimizer/src/indexgenerator.cpp index 4bf9fccad..1b707bb1b 100644 --- a/3rdparty/meshoptimizer/src/indexgenerator.cpp +++ b/3rdparty/meshoptimizer/src/indexgenerator.cpp @@ -127,6 +127,30 @@ struct VertexCustomHasher } }; +struct TriangleKey +{ + unsigned int a, b, c; + + bool operator==(const TriangleKey& other) const + { + return a == other.a && b == other.b && c == other.c; + } +}; + +struct TriangleKeyHasher +{ + size_t hash(const TriangleKey& k) const + { + // Optimized Spatial Hashing for Collision Detection of Deformable Objects + return (k.a * 73856093) ^ (k.b * 19349663) ^ (k.c * 83492791); + } + + bool equal(const TriangleKey& lhs, const TriangleKey& rhs) const + { + return lhs == rhs; + } +}; + struct EdgeHasher { const unsigned int* remap; @@ -304,6 +328,63 @@ static void generateShadowBuffer(unsigned int* destination, const unsigned int* } } +static size_t filterIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const unsigned int* remap, size_t vertex_count, meshopt_Allocator& allocator) +{ + (void)vertex_count; + + TriangleKeyHasher hasher = {}; + TriangleKey empty = {~0u, ~0u, ~0u}; + + size_t face_count = index_count / 3; + size_t table_size = hashBuckets(face_count); + TriangleKey* table = allocator.allocate(table_size); + memset(table, -1, table_size * sizeof(TriangleKey)); + + size_t write = 0; + + for (size_t i = 0; i < face_count; ++i) + { + unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2]; + assert(a < vertex_count && b < vertex_count && c < vertex_count); + + unsigned int ra = remap[a], rb = remap[b], rc = remap[c]; + + // degenerate triangle + if (ra == rb || ra == rc || rb == rc) + continue; + + // canonicalize triangle wrt winding preserving rotation + if (rb < ra && rb < rc) + { + // abc -> bca + unsigned int t = ra; + ra = rb, rb = rc, rc = t; + } + else if (rc < ra && rc < rb) + { + // abc -> cab + unsigned int t = rc; + rc = rb, rb = ra, ra = t; + } + + TriangleKey key = {ra, rb, rc}; + TriangleKey* entry = hashLookup(table, table_size, hasher, key, empty); + + // duplicate triangle + if (entry->a != ~0u) + continue; + + *entry = key; + + destination[write + 0] = a; + destination[write + 1] = b; + destination[write + 2] = c; + write += 3; + } + + return write; +} + } // namespace meshopt size_t meshopt_generateVertexRemap(unsigned int* destination, const unsigned int* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size) @@ -311,7 +392,6 @@ size_t meshopt_generateVertexRemap(unsigned int* destination, const unsigned int using namespace meshopt; assert(indices || index_count == vertex_count); - assert(!indices || index_count % 3 == 0); assert(vertex_size > 0 && vertex_size <= 256); meshopt_Allocator allocator; @@ -325,7 +405,6 @@ size_t meshopt_generateVertexRemapMulti(unsigned int* destination, const unsigne using namespace meshopt; assert(indices || index_count == vertex_count); - assert(index_count % 3 == 0); assert(stream_count > 0 && stream_count <= 16); for (size_t i = 0; i < stream_count; ++i) @@ -345,7 +424,6 @@ size_t meshopt_generateVertexRemapCustom(unsigned int* destination, const unsign using namespace meshopt; assert(indices || index_count == vertex_count); - assert(!indices || index_count % 3 == 0); assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); assert(vertex_positions_stride % sizeof(float) == 0); @@ -393,8 +471,6 @@ void meshopt_remapVertexBuffer(void* destination, const void* vertices, size_t v void meshopt_remapIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const unsigned int* remap) { - assert(index_count % 3 == 0); - for (size_t i = 0; i < index_count; ++i) { unsigned int index = indices ? indices[i] : unsigned(i); @@ -404,11 +480,49 @@ void meshopt_remapIndexBuffer(unsigned int* destination, const unsigned int* ind } } +size_t meshopt_filterIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride) +{ + using namespace meshopt; + + assert(index_count % 3 == 0); + assert(vertex_size > 0 && vertex_size <= 256); + assert(vertex_size <= vertex_stride); + + meshopt_Allocator allocator; + VertexHasher hasher = {static_cast(vertices), vertex_size, vertex_stride}; + + unsigned int* remap = allocator.allocate(vertex_count); + generateVertexRemap(remap, indices, index_count, vertex_count, hasher, allocator); + + return filterIndexBuffer(destination, indices, index_count, remap, vertex_count, allocator); +} + +size_t meshopt_filterIndexBufferMulti(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, const struct meshopt_Stream* streams, size_t stream_count) +{ + using namespace meshopt; + + assert(index_count % 3 == 0); + assert(stream_count > 0 && stream_count <= 16); + + for (size_t i = 0; i < stream_count; ++i) + { + assert(streams[i].size > 0 && streams[i].size <= 256); + assert(streams[i].size <= streams[i].stride); + } + + meshopt_Allocator allocator; + VertexStreamHasher hasher = {streams, stream_count}; + + unsigned int* remap = allocator.allocate(vertex_count); + generateVertexRemap(remap, indices, index_count, vertex_count, hasher, allocator); + + return filterIndexBuffer(destination, indices, index_count, remap, vertex_count, allocator); +} + void meshopt_generateShadowIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride) { using namespace meshopt; - assert(indices); assert(index_count % 3 == 0); assert(vertex_size > 0 && vertex_size <= 256); assert(vertex_size <= vertex_stride); @@ -423,7 +537,6 @@ void meshopt_generateShadowIndexBufferMulti(unsigned int* destination, const uns { using namespace meshopt; - assert(indices); assert(index_count % 3 == 0); assert(stream_count > 0 && stream_count <= 16); diff --git a/3rdparty/meshoptimizer/src/meshletcodec.cpp b/3rdparty/meshoptimizer/src/meshletcodec.cpp new file mode 100644 index 000000000..f6984a77b --- /dev/null +++ b/3rdparty/meshoptimizer/src/meshletcodec.cpp @@ -0,0 +1,1051 @@ +// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details +#include "meshoptimizer.h" + +#include +#include + +// The block below auto-detects SIMD ISA that can be used on the target platform +#ifndef MESHOPTIMIZER_NO_SIMD + +// The SIMD implementation requires SSE4.1, which can be enabled unconditionally through compiler settings +#if defined(__AVX__) || defined(__SSE4_1__) +#define SIMD_SSE +#endif + +// MSVC supports compiling SSE4.1 code regardless of compile options; we use a cpuid-based scalar fallback +#if !defined(SIMD_SSE) && defined(_MSC_VER) && !defined(__clang__) && (defined(_M_IX86) || (defined(_M_X64) && !defined(_M_ARM64EC))) +#define SIMD_SSE +#define SIMD_FALLBACK +#endif + +// GCC 4.9+ and clang 3.8+ support targeting SIMD ISA from individual functions; we use a cpuid-based scalar fallback +#if !defined(SIMD_SSE) && ((defined(__clang__) && __clang_major__ * 100 + __clang_minor__ >= 308) || (defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ >= 409)) && (defined(__i386__) || defined(__x86_64__)) +#define SIMD_SSE +#define SIMD_FALLBACK +#define SIMD_TARGET __attribute__((target("sse4.1"))) +#endif + +// When targeting AArch64, enable NEON SIMD unconditionally; we do not support SIMD decoding for 32-bit ARM +#if defined(__aarch64__) || (defined(_MSC_VER) && (defined(_M_ARM64) || defined(_M_ARM64EC)) && _MSC_VER >= 1922) +#define SIMD_NEON +#endif + +#if defined(_MSC_VER) && !defined(__clang__) && _MSC_VER > 1930 +#define SIMD_FLATTEN [[msvc::flatten]] +#elif defined(__GNUC__) || defined(__clang__) +#define SIMD_FLATTEN __attribute__((flatten)) +#else +#define SIMD_FLATTEN +#endif + +#ifndef SIMD_TARGET +#define SIMD_TARGET +#endif + +#endif // !MESHOPTIMIZER_NO_SIMD + +#ifdef SIMD_SSE +#include +#endif + +#ifdef SIMD_NEON +#include +#endif + +#if defined(SIMD_SSE) && defined(SIMD_FALLBACK) +#ifdef _MSC_VER +#include // __cpuid +#else +#include // __cpuid +#endif +#endif + +#ifndef TRACE +#define TRACE 0 +#endif + +#if TRACE +#include +#endif + +namespace meshopt +{ + +typedef unsigned int EdgeFifo8[8][2]; + +static int rotateTriangle(unsigned int a, unsigned int b, unsigned int c) +{ + return (a > b && a > c) ? 1 : (b > c ? 2 : 0); +} + +static int getEdgeFifo8(EdgeFifo8 fifo, unsigned int a, unsigned int b, unsigned int c, size_t offset) +{ + for (int i = 0; i < 8; ++i) + { + size_t index = (offset - 1 - i) & 7; + + unsigned int e0 = fifo[index][0]; + unsigned int e1 = fifo[index][1]; + + if (e0 == a && e1 == b) + return (i << 2) | 0; + if (e0 == b && e1 == c) + return (i << 2) | 1; + if (e0 == c && e1 == a) + return (i << 2) | 2; + } + + return -1; +} + +static void pushEdgeFifo8(EdgeFifo8 fifo, unsigned int a, unsigned int b, size_t& offset) +{ + fifo[offset][0] = a; + fifo[offset][1] = b; + offset = (offset + 1) & 7; +} + +static size_t encodeTriangles(unsigned char* codes, unsigned char* extra, const unsigned char* triangles, size_t triangle_count) +{ + EdgeFifo8 edgefifo; + memset(edgefifo, -1, sizeof(edgefifo)); + + size_t edgefifooffset = 0; + + unsigned int next = 0; + + // 4-bit triangle codes give us 16 options that we use as follows: + // 3*2 edge reuse (2 edges * 3 last triangles) * 2 next/explicit = 12 options + // 4 remaining options = next bits; 000, 001, 011, 111. + // triangles are rotated to make next bits line up. + memset(codes, 0, (triangle_count + 1) / 2); + + static const int rotations[] = {0, 1, 2, 0, 1}; + + unsigned char* start = extra; + + for (size_t i = 0; i < triangle_count; ++i) + { +#if TRACE > 1 + unsigned int last = next; +#endif + + int fer = getEdgeFifo8(edgefifo, triangles[i * 3 + 0], triangles[i * 3 + 1], triangles[i * 3 + 2], edgefifooffset); + + if (fer >= 0 && (fer >> 2) < 6) + { + // note: getEdgeFifo8 implicitly rotates triangles by matching a/b to existing edge + const int* order = rotations + (fer & 3); + + unsigned int a = triangles[i * 3 + order[0]], b = triangles[i * 3 + order[1]], c = triangles[i * 3 + order[2]]; + + int fec = (c == next) ? (next++, 0) : 1; + +#if TRACE > 1 + printf("%3d+ | %3d %3d %3d | edge: e%d c%d\n", last, a, b, c, fer >> 2, fec); +#endif + + unsigned int code = (fer >> 2) * 2 + fec; + + codes[i / 2] |= (unsigned char)(code << ((i & 1) * 4)); + + if (fec) + *extra++ = (unsigned char)c; + + pushEdgeFifo8(edgefifo, c, b, edgefifooffset); + pushEdgeFifo8(edgefifo, a, c, edgefifooffset); + } + else + { + // rotate triangles to minimize the need for extra vertices + int rotation = rotateTriangle(triangles[i * 3 + 0], triangles[i * 3 + 1], triangles[i * 3 + 2]); + const int* order = rotations + rotation; + + unsigned int a = triangles[i * 3 + order[0]], b = triangles[i * 3 + order[1]], c = triangles[i * 3 + order[2]]; + + // fe must be continuous: once a vertex is encoded with next, further vertices must also be encoded with next + int fea = (a == next && b == next + 1 && c == next + 2) ? (next++, 0) : 1; + int feb = (b == next && c == next + 1) ? (next++, 0) : 1; + int fec = (c == next) ? (next++, 0) : 1; + + assert(fea == 1 || feb == 0); + assert(feb == 1 || fec == 0); + +#if TRACE > 1 + printf("%3d+ | %3d %3d %3d | restart: %d%d%d\n", last, a, b, c, fea, feb, fec); +#endif + + unsigned int code = 12 + (fea + feb + fec); + + codes[i / 2] |= (unsigned char)(code << ((i & 1) * 4)); + + if (fea) + *extra++ = (unsigned char)a; + if (feb) + *extra++ = (unsigned char)b; + if (fec) + *extra++ = (unsigned char)c; + + pushEdgeFifo8(edgefifo, c, b, edgefifooffset); + pushEdgeFifo8(edgefifo, a, c, edgefifooffset); + } + } + + return extra - start; +} + +static size_t encodeVertices(unsigned char* ctrl, unsigned char* data, const unsigned int* vertices, size_t vertex_count) +{ + // grouped varint, 2 bit per value to indicate 0/1/2/3 byte deltas, with per-group 4-byte fallback + memset(ctrl, 0, (vertex_count + 3) / 4); + + unsigned char* start = data; + + unsigned int last = ~0u; + + for (size_t i = 0; i < vertex_count; i += 4) + { + unsigned int gv[4] = {}; + + for (int k = 0; k < 4 && i + k < vertex_count; ++k) + { + unsigned int d = vertices[i + k] - last - 1; + unsigned int v = (d << 1) ^ (int(d) >> 31); + + gv[k] = v; + last = vertices[i + k]; + } + + // if any value needs 4 bytes, or if *all* values need 3 bytes, we use 4 bytes for all values + // this allows us to encode most 3-byte deltas with 3 bytes which saves space overall + bool use4 = (gv[0] | gv[1] | gv[2] | gv[3]) > 0xffffff || (gv[0] > 0xffff && gv[1] > 0xffff && gv[2] > 0xffff && gv[3] > 0xffff); + + for (int k = 0; k < 4; ++k) + { + unsigned int v = gv[k]; + + // 0/1/2/3 bytes per value, or all 4 values use 4 bytes + int code = use4 ? 3 : (v == 0 ? 0 : (v < 256 ? 1 : (v < 65536 ? 2 : 3))); + + if (code > 0) + *data++ = (unsigned char)(v & 0xff); + if (code > 1) + *data++ = (unsigned char)((v >> 8) & 0xff); + if (code > 2) + *data++ = (unsigned char)((v >> 16) & 0xff); + if (use4) + *data++ = (unsigned char)((v >> 24) & 0xff); + + // split low and high bits into two nibbles for better packing + ctrl[i / 4] |= ((code & 1) << k) | ((code >> 1) << (k + 4)); + } + } + + return data - start; +} + +#if defined(SIMD_FALLBACK) || (!defined(SIMD_SSE) && !defined(SIMD_NEON)) +inline void writeTriangle(unsigned int* triangles, size_t i, unsigned int fifo) +{ + // output triangle is stored without extra edge vertex (0xcbac => 0xcba) + triangles[i] = fifo >> 8; +} + +inline void writeTriangle(unsigned char* triangles, size_t i, unsigned int fifo) +{ + triangles[i * 3 + 0] = (unsigned char)(fifo >> 8); + triangles[i * 3 + 1] = (unsigned char)(fifo >> 16); + triangles[i * 3 + 2] = (unsigned char)(fifo >> 24); +} + +template +static const unsigned char* decodeTriangles(T* triangles, const unsigned char* codes, const unsigned char* extra, const unsigned char* bound, size_t triangle_count) +{ + // branchlessly read next or extra vertex and advance pointers +#define NEXT(var, ec) \ + e = *extra; \ + unsigned int var = (ec) ? e : next; \ + extra += (ec), next += 1 - (ec) + + unsigned int next = 0; + unsigned int fifo[3] = {}; // two edge fifo entries in one uint: 0xcbac + + for (size_t i = 0; i < triangle_count; ++i) + { + if (extra > bound) + return NULL; + + unsigned int code = (codes[i / 2] >> ((i & 1) * 4)) & 0xF; + unsigned int tri; + + if (code < 12) + { + // reuse + unsigned int edge = fifo[code / 4]; + edge >>= (code << 3) & 16; // shift by 16 if bit 1 is set (odd edge for each triangle) + + // 0-1 extra vertices + unsigned int e; + NEXT(c, code & 1); + + // repack triangle into edge format (0xcbac) + tri = ((edge & 0xff) << 16) | (edge & 0xff00) | c | (c << 24); + } + else + { + // restart + int fea = code > 12; + int feb = code > 13; + int fec = code > 14; + + // 0-3 extra vertices + unsigned int e; + NEXT(a, fea); + NEXT(b, feb); + NEXT(c, fec); + + // repack triangle into edge format (0xcbac) + tri = c | (a << 8) | (b << 16) | (c << 24); + } + + writeTriangle(triangles, i, tri); + + fifo[2] = fifo[1]; + fifo[1] = fifo[0]; + fifo[0] = tri; + } + + return extra; + +#undef NEXT +} + +template +static const unsigned char* decodeVertices(V* vertices, const unsigned char* ctrl, const unsigned char* data, const unsigned char* bound, size_t vertex_count) +{ + unsigned int last = ~0u; + + for (size_t i = 0; i < vertex_count; i += 4) + { + if (data > bound) + return NULL; + + unsigned char code4 = ctrl[i / 4]; + + for (int k = 0; k < 4; ++k) + { + int code = ((code4 >> k) & 1) | ((code4 >> (k + 3)) & 2); + int length = code4 == 0xff ? 4 : code; + + // branchlessly read up to 4 bytes + unsigned int mask = (length == 4) ? ~0u : (1 << (8 * length)) - 1; + unsigned int v = (data[0] | (data[1] << 8) | (data[2] << 16) | (data[3] << 24)) & mask; + + // unzigzag + 1 + unsigned int d = (v >> 1) ^ -int(v & 1); + unsigned int r = last + d + 1; + + if (i + k < vertex_count) + vertices[i + k] = V(r); + + data += length; + last = r; + } + } + + return data; +} + +static int decodeMeshlet(void* vertices, void* triangles, const unsigned char* codes, const unsigned char* ctrl, const unsigned char* data, const unsigned char* bound, size_t vertex_count, size_t triangle_count, size_t vertex_size, size_t triangle_size) +{ + if (vertex_size == 4) + data = decodeVertices(static_cast(vertices), ctrl, data, bound, vertex_count); + else + data = decodeVertices(static_cast(vertices), ctrl, data, bound, vertex_count); + if (!data) + return -2; + + if (triangle_size == 4) + data = decodeTriangles(static_cast(triangles), codes, data, bound, triangle_count); + else + data = decodeTriangles(static_cast(triangles), codes, data, bound, triangle_count); + if (!data) + return -2; + + return (data == bound) ? 0 : -3; +} +#endif + +#if defined(SIMD_SSE) || defined(SIMD_NEON) +// SIMD state is stored in a single 16b register as follows: +// 0..5: 6 next extra bytes +// 6..14: 9 bytes = 3 triangles worth of index data +// 15: 'next' byte + +// upon reading each triangle pair we need to transform this state such that the 9 bytes with triangle data contain the newly decoded triangles, +// which is a permutation of original state modulo per-element additions +// this transform can be chained to decode second triangle from original state; we create tables for 256 combinations of two 4-bit triangle codes +// the actual decoding becomes shuffle+add per triangle pair, plus management of extra bytes +static unsigned char kDecodeTableMasks[256][16]; +static unsigned char kDecodeTableExtra[256]; + +// for SIMD vertex decoding we need to unpack 4 values with 0-4 bytes in each +// this can be done with a single control-dependent shuffle per group +static unsigned char kDecodeTableVerts[256][16]; +static unsigned char kDecodeTableLength[256]; + +static bool decodeBuildTables() +{ +#define NEXT(var, ec) \ + shuf[var] = (ec) ? (unsigned char)extra : 15; \ + next[var] = (ec) ? 0 : (unsigned char)nextoff; \ + extra += (ec), nextoff += 1 - (ec) + + // check for SSE4.1 support if we have a fallback path +#if defined(SIMD_SSE) && defined(SIMD_FALLBACK) + int cpuinfo[4] = {}; +#ifdef _MSC_VER + __cpuid(cpuinfo, 1); +#else + __cpuid(1, cpuinfo[0], cpuinfo[1], cpuinfo[2], cpuinfo[3]); +#endif + // bit 19 = SSE4.1 + if ((cpuinfo[2] & (1 << 19)) == 0) + return false; +#endif + + // fill triangle decoding tables for each combination of two triangle codes + for (int code = 0; code < 256; ++code) + { + unsigned char shuf[16] = {}; + unsigned char next[16] = {}; + int extra = 0; + int nextoff = 0; + + // state 0..5 will be refilled every iteration, so we ignore that + // state 6..8 will always contain the last decoded triangle because every triangle shifts fifo equally, so we can decode it independently + shuf[6] = 12; + shuf[7] = 13; + shuf[8] = 14; + + // state 15 will contain next (potentially incremented a few times) + shuf[15] = 15; + + // state 9..11 will contain the first decoded triangle (tri0), which can refer to extra/next and the original triangle history + // state 12..14 will contain the second decoded triangle (tri1); when decoding edge reuse, we need to handle edge 0/1 specially as it was just decoded earlier + for (int k = 0; k < 2; ++k) + { + int tri = (code >> (k * 4)) & 0xf; + + if (tri < 12) + { + if (k == 1 && tri / 4 == 0) + { + // we need to decode one of two edges from the triangle we just decoded earlier + // for that we simply need to copy shuf/next values for the two decoded indices + shuf[9 + k * 3] = shuf[9 + ((tri & 2) ? 2 : 0)]; + next[9 + k * 3] = next[9 + ((tri & 2) ? 2 : 0)]; + + shuf[10 + k * 3] = shuf[9 + ((tri & 2) ? 1 : 2)]; + next[10 + k * 3] = next[9 + ((tri & 2) ? 1 : 2)]; + } + else + { + // reuse: edge comes from the history based on edge index + // note: we reuse with an offset because last triangle in the original history was consumed by tri0 + int trioff = 6 + k * 3 + (2 - tri / 4) * 3; + + // edge cb or ac + shuf[9 + k * 3] = (unsigned char)(trioff + ((tri & 2) ? 2 : 0)); + shuf[10 + k * 3] = (unsigned char)(trioff + ((tri & 2) ? 1 : 2)); + } + + // third vertex is either next or comes from extra + NEXT(11 + k * 3, tri & 1); + } + else + { + // restart: three vertices, each comes from next or extra + int fea = tri > 12; + int feb = tri > 13; + int fec = tri > 14; + + NEXT(9 + k * 3, fea); + NEXT(10 + k * 3, feb); + NEXT(11 + k * 3, fec); + } + } + + // next needs to advance + next[15] = (unsigned char)nextoff; + + // next[0..8] = 0 trivially (never written to); next[9] must also be 0 because nextoff is 0 initially + // shuf[0..5] is not used, which allows us to pack next[10..15] + shuf[6..15] into a single 16-byte entry + assert(next[9] == 0); + memcpy(&kDecodeTableMasks[code][0], &next[10], 6); + memcpy(&kDecodeTableMasks[code][6], &shuf[6], 10); + kDecodeTableExtra[code] = (unsigned char)extra; + } + + // fill vertex decoding tables for each combination of four vertex references + for (unsigned int i = 0; i < 256; ++i) + { + unsigned char shuf[16] = {}; + int offset = 0; + + for (int k = 0; k < 4; ++k) + { + int code = ((i >> k) & 1) | ((i >> (k + 3)) & 2); + int length = i == 0xff ? 4 : code; // 0/1/2/3 bytes, or all 4 bytes if code==0xff + + shuf[k * 4 + 0] = (length > 0) ? (unsigned char)(offset + 0) : 0x80; + shuf[k * 4 + 1] = (length > 1) ? (unsigned char)(offset + 1) : 0x80; + shuf[k * 4 + 2] = (length > 2) ? (unsigned char)(offset + 2) : 0x80; + shuf[k * 4 + 3] = (length > 3) ? (unsigned char)(offset + 3) : 0x80; + + offset += length; + } + + memcpy(kDecodeTableVerts[i], shuf, sizeof(shuf)); + kDecodeTableLength[i] = (unsigned char)offset; + } + + return true; + +#undef NEXT +} + +static bool gDecodeTablesInitialized = decodeBuildTables(); +#endif + +#if defined(SIMD_SSE) +SIMD_TARGET +inline __m128i decodeTriangleGroup(__m128i state, unsigned char code, const unsigned char*& extra) +{ + __m128i shuf = _mm_loadu_si128(reinterpret_cast(kDecodeTableMasks[code])); + __m128i next = _mm_slli_si128(shuf, 10); + + // patch first 6 bytes with current extra and roll state forward + __m128i ext = _mm_loadl_epi64(reinterpret_cast(extra)); + state = _mm_blend_epi16(state, ext, 7); + state = _mm_add_epi8(_mm_shuffle_epi8(state, shuf), next); + + extra += kDecodeTableExtra[code]; + + return state; +} + +SIMD_TARGET +inline __m128i decodeVertexGroup(__m128i last, unsigned char code, const unsigned char*& data) +{ + __m128i word = _mm_loadu_si128(reinterpret_cast(data)); + __m128i shuf = _mm_loadu_si128(reinterpret_cast(kDecodeTableVerts[code])); + + __m128i v = _mm_shuffle_epi8(word, shuf); + + // unzigzag+1 + __m128i xl = _mm_sub_epi32(_mm_setzero_si128(), _mm_and_si128(v, _mm_set1_epi32(1))); + __m128i xr = _mm_srli_epi32(v, 1); + __m128i x = _mm_add_epi32(_mm_xor_si128(xl, xr), _mm_set1_epi32(1)); + + // prefix sum + x = _mm_add_epi32(x, _mm_slli_si128(x, 8)); + x = _mm_add_epi32(x, _mm_slli_si128(x, 4)); + x = _mm_add_epi32(x, _mm_shuffle_epi32(last, 0xff)); + + data += kDecodeTableLength[code]; + + return x; +} +#endif + +#if defined(SIMD_NEON) +SIMD_TARGET +inline uint8x16_t decodeTriangleGroup(uint8x16_t state, unsigned char code, const unsigned char*& extra) +{ + uint8x16_t shuf = vld1q_u8(kDecodeTableMasks[code]); + uint8x16_t next = vextq_u8(vdupq_n_u8(0), shuf, 6); + + // patch first 6 bytes with current extra and roll state forward + uint8x8_t extl = vld1_u8(extra); + uint8x16_t ext = vcombine_u8(extl, vdup_n_u8(0)); + state = vbslq_u8(vcombine_u8(vcreate_u8(0xffffffffffffull), vdup_n_u8(0)), ext, state); + state = vaddq_u8(vqtbl1q_u8(state, shuf), next); + + extra += kDecodeTableExtra[code]; + + return state; +} + +SIMD_TARGET +inline uint32x4_t decodeVertexGroup(uint32x4_t last, unsigned char code, const unsigned char*& data) +{ + uint8x16_t word = vld1q_u8(data); + uint8x16_t shuf = vld1q_u8(kDecodeTableVerts[code]); + + uint32x4_t v = vreinterpretq_u32_u8(vqtbl1q_u8(word, shuf)); + + // unzigzag+1 + uint32x4_t xl = vsubq_u32(vdupq_n_u32(0), vandq_u32(v, vdupq_n_u32(1))); + uint32x4_t xr = vshrq_n_u32(v, 1); + uint32x4_t x = vaddq_u32(veorq_u32(xl, xr), vdupq_n_u32(1)); + + // prefix sum + x = vaddq_u32(x, vextq_u32(vdupq_n_u32(0), x, 2)); + x = vaddq_u32(x, vextq_u32(vdupq_n_u32(0), x, 3)); + x = vaddq_u32(x, vdupq_n_u32(vgetq_lane_u32(last, 3))); + + data += kDecodeTableLength[code]; + + return x; +} +#endif + +#if defined(SIMD_SSE) +#ifdef __GNUC__ +typedef int __attribute__((aligned(1))) unaligned_int; +#else +typedef int unaligned_int; +#endif +#endif + +#if defined(SIMD_SSE) || defined(SIMD_NEON) +SIMD_TARGET +static const unsigned char* decodeTrianglesSimd(unsigned int* triangles, const unsigned char* codes, const unsigned char* extra, const unsigned char* bound, size_t triangle_count) +{ +#if defined(SIMD_SSE) + __m128i repack = _mm_setr_epi8(9, 10, 11, -1, 12, 13, 14, -1, 0, 0, 0, 0, 0, 0, 0, 0); + __m128i state = _mm_setzero_si128(); +#elif defined(SIMD_NEON) + uint8x8_t repack = vcreate_u8(0xff0e0d0cff0b0a09ull); + uint8x16_t state = vdupq_n_u8(0); +#endif + + size_t groups = triangle_count / 2; + + // process all complete groups + for (size_t i = 0; i < groups; ++i) + { + unsigned char code = *codes++; + + if (extra > bound) + return NULL; + + state = decodeTriangleGroup(state, code, extra); + + // write 6 bytes of new triangle data into output, formatted as 8 bytes with 0 padding +#if defined(SIMD_SSE) + __m128i r = _mm_shuffle_epi8(state, repack); + _mm_storel_epi64(reinterpret_cast<__m128i*>(&triangles[i * 2]), r); +#elif defined(SIMD_NEON) + uint32x2_t r = vreinterpret_u32_u8(vqtbl1_u8(state, repack)); + vst1_u32(&triangles[i * 2], r); +#endif + } + + // process a 1 triangle tail; to maintain the memory safety guarantee we have to write a 32-bit element + if (triangle_count & 1) + { + unsigned char code = *codes++; + + if (extra > bound) + return NULL; + + state = decodeTriangleGroup(state, code, extra); + + unsigned int* tail = &triangles[triangle_count & ~1u]; + +#if defined(SIMD_SSE) + __m128i r = _mm_shuffle_epi8(state, repack); + *tail = unsigned(_mm_cvtsi128_si32(r)); +#elif defined(SIMD_NEON) + uint32x2_t r = vreinterpret_u32_u8(vqtbl1_u8(state, repack)); + vst1_lane_u32(tail, r, 0); +#endif + } + + return extra; +} + +SIMD_TARGET +static const unsigned char* decodeTrianglesSimd(unsigned char* triangles, const unsigned char* codes, const unsigned char* extra, const unsigned char* bound, size_t triangle_count) +{ +#if defined(SIMD_SSE) + __m128i state = _mm_setzero_si128(); +#elif defined(SIMD_NEON) + uint8x16_t state = vdupq_n_u8(0); +#endif + + // because the output buffer is guaranteed to have 32-bit aligned size available, we can optimize writes and tail processing + // instead of processing triangles 2 at a time, we process 2 *pairs* at a time (12-byte write) followed by a tail pair, if present + // if the number of triangles mod 4 is 3, we'd normally need to write 12k+9 bytes, but we can instead overwrite up to 3 bytes in the main loop + size_t groups = (triangle_count + 1) / 4; + + // process all complete groups + for (size_t i = 0; i < groups; ++i) + { + unsigned char code0 = *codes++; + unsigned char code1 = *codes++; + + // each triangle pair reads <=6 bytes from extra, so two pairs need <=12 bytes and gap guarantees 16 byte of overread + if (extra > bound) + return NULL; + + state = decodeTriangleGroup(state, code0, extra); + + // write first decoded triangle and first index of second decoded triangle +#if defined(SIMD_SSE) + __m128i r0 = _mm_srli_si128(state, 9); + *reinterpret_cast(&triangles[i * 12]) = _mm_cvtsi128_si32(r0); +#elif defined(SIMD_NEON) + uint8x16_t r0 = vextq_u8(state, vdupq_n_u8(0), 9); + vst1q_lane_u32(reinterpret_cast(&triangles[i * 12]), vreinterpretq_u32_u8(r0), 0); +#endif + + state = decodeTriangleGroup(state, code1, extra); + + // write last two indices of second decoded triangle that we didn't write above plus two new ones + // note that the second decoded triangle has shifted down to 6-8 bytes, hence shift by 7 +#if defined(SIMD_SSE) + __m128i r1 = _mm_srli_si128(state, 7); + _mm_storel_epi64(reinterpret_cast<__m128i*>(&triangles[i * 12 + 4]), r1); +#elif defined(SIMD_NEON) + uint8x16_t r1 = vextq_u8(state, vdupq_n_u8(0), 7); + vst1_u8(&triangles[i * 12 + 4], vget_low_u8(r1)); +#endif + } + + // process a 1-2 triangle tail; to maintain the memory safety guarantee we have to write 1-2 32-bit elements + if (groups * 4 < triangle_count) + { + unsigned char code = *codes++; + + if (extra > bound) + return NULL; + + state = decodeTriangleGroup(state, code, extra); + + unsigned char* tail = &triangles[(triangle_count & ~3u) * 3]; + +#if defined(SIMD_SSE) + __m128i r = _mm_srli_si128(state, 9); + + *reinterpret_cast(tail) = _mm_cvtsi128_si32(r); + if ((triangle_count & 3) > 1) + *reinterpret_cast(tail + 4) = _mm_extract_epi32(r, 1); +#elif defined(SIMD_NEON) + uint8x16_t r = vextq_u8(state, vdupq_n_u8(0), 9); + + vst1q_lane_u32(reinterpret_cast(tail), vreinterpretq_u32_u8(r), 0); + if ((triangle_count & 3) > 1) + vst1q_lane_u32(reinterpret_cast(tail + 4), vreinterpretq_u32_u8(r), 1); +#endif + } + + return extra; +} + +SIMD_TARGET +static const unsigned char* decodeVerticesSimd(unsigned int* vertices, const unsigned char* ctrl, const unsigned char* data, const unsigned char* bound, size_t vertex_count) +{ +#if defined(SIMD_SSE) + __m128i last = _mm_set1_epi32(-1); +#elif defined(SIMD_NEON) + uint32x4_t last = vdupq_n_u32(~0u); +#endif + + size_t groups = vertex_count / 4; + + // process all complete groups + for (size_t i = 0; i < groups; ++i) + { + unsigned char code = *ctrl++; + if (data > bound) + return NULL; + + last = decodeVertexGroup(last, code, data); + +#if defined(SIMD_SSE) + _mm_storeu_si128(reinterpret_cast<__m128i*>(&vertices[i * 4]), last); +#elif defined(SIMD_NEON) + vst1q_u32(&vertices[i * 4], last); +#endif + } + + // process a 1-3 vertex tail; to maintain the memory safety guarantee we have to write individual elements + if (vertex_count & 3) + { + unsigned char code = *ctrl++; + + if (data > bound) + return NULL; + + last = decodeVertexGroup(last, code, data); + + unsigned int* tail = &vertices[vertex_count & ~3u]; + +#if defined(SIMD_SSE) + tail[0] = _mm_cvtsi128_si32(last); + if ((vertex_count & 3) > 1) + tail[1] = _mm_extract_epi32(last, 1); + if ((vertex_count & 3) > 2) + tail[2] = _mm_extract_epi32(last, 2); +#elif defined(SIMD_NEON) + vst1q_lane_u32(&tail[0], last, 0); + if ((vertex_count & 3) > 1) + vst1q_lane_u32(&tail[1], last, 1); + if ((vertex_count & 3) > 2) + vst1q_lane_u32(&tail[2], last, 2); +#endif + } + + return data; +} + +SIMD_TARGET +static const unsigned char* decodeVerticesSimd(unsigned short* vertices, const unsigned char* ctrl, const unsigned char* data, const unsigned char* bound, size_t vertex_count) +{ +#if defined(SIMD_SSE) + __m128i repack = _mm_setr_epi8(0, 1, 4, 5, 8, 9, 12, 13, 0, 0, 0, 0, 0, 0, 0, 0); + __m128i last = _mm_set1_epi32(-1); +#elif defined(SIMD_NEON) + uint32x4_t last = vdupq_n_u32(~0u); +#endif + + // because the output buffer is guaranteed to have 32-bit aligned size available, we can simplify tail processing + // if the number of vertices mod 4 is 3, we'd normally need to write 8+6 bytes, but we can instead overwrite up to 2 bytes in the main loop + size_t groups = (vertex_count + 1) / 4; + + // process all complete groups + for (size_t i = 0; i < groups; ++i) + { + unsigned char code = *ctrl++; + + if (data > bound) + return NULL; + + last = decodeVertexGroup(last, code, data); + +#if defined(SIMD_SSE) + __m128i r = _mm_shuffle_epi8(last, repack); + _mm_storel_epi64(reinterpret_cast<__m128i*>(&vertices[i * 4]), r); +#elif defined(SIMD_NEON) + uint16x4_t r = vmovn_u32(last); + vst1_u16(&vertices[i * 4], r); +#endif + } + + // process a 1-2 vertex tail; to maintain the memory safety guarantee we have to write a 32-bit element + if (groups * 4 < vertex_count) + { + unsigned char code = *ctrl++; + + if (data > bound) + return NULL; + + last = decodeVertexGroup(last, code, data); + + unsigned short* tail = &vertices[vertex_count & ~3u]; + +#if defined(SIMD_SSE) + __m128i r = _mm_shufflelo_epi16(last, 8); + *reinterpret_cast(tail) = _mm_cvtsi128_si32(r); +#elif defined(SIMD_NEON) + uint16x4_t r = vmovn_u32(last); + vst1_lane_u32(reinterpret_cast(tail), vreinterpret_u32_u16(r), 0); +#endif + } + + return data; +} + +template +SIMD_TARGET SIMD_FLATTEN static int +decodeMeshletSimd(void* vertices, void* triangles, const unsigned char* codes, const unsigned char* ctrl, const unsigned char* data, const unsigned char* bound, size_t vertex_count, size_t triangle_count, size_t vertex_size, size_t triangle_size) +{ + assert(gDecodeTablesInitialized); + (void)gDecodeTablesInitialized; + +#ifdef __clang__ + // data is guaranteed to be non-null initially; if decode loops never hit bounds errors, it remains non-null + __builtin_assume(data); +#endif + + // decodes 4 vertices at a time with tail processing; writes up to align(vertex_size * vertex_count, 4) + // raw decoding skips tail processing by rounding up vertex count; it's safe because output buffer is guaranteed to have extra space, and tail control data is 0 + if (vertex_size == 4 || Raw) + data = decodeVerticesSimd(static_cast(vertices), ctrl, data, bound, Raw ? (vertex_count + 3) & ~3 : vertex_count); + else + data = decodeVerticesSimd(static_cast(vertices), ctrl, data, bound, vertex_count); + if (!data) + return -2; + + // decodes 2/4 triangles at a time with tail processing; writes up to align(triangle_size * triangle_count, 4) + // raw decoding skips tail processing by rounding up triangle count; it's safe because output buffer is guaranteed to have extra space, and tail code data is 0 + if (triangle_size == 4 || Raw) + data = decodeTrianglesSimd(static_cast(triangles), codes, data, bound, Raw ? (triangle_count + 1) & ~1 : triangle_count); + else + data = decodeTrianglesSimd(static_cast(triangles), codes, data, bound, triangle_count); + if (!data) + return -2; + + return (data == bound) ? 0 : -3; +} +#endif + +} // namespace meshopt + +size_t meshopt_encodeMeshletBound(size_t max_vertices, size_t max_triangles) +{ + size_t codes_size = (max_triangles + 1) / 2; + size_t extra_size = max_triangles * 3; + + size_t ctrl_size = (max_vertices + 3) / 4; + size_t data_size = (max_vertices + 3) / 4 * 16; // worst case: 16 bytes per vertex group + + size_t gap_size = (codes_size + ctrl_size < 16) ? 16 - (codes_size + ctrl_size) : 0; + + return codes_size + extra_size + ctrl_size + data_size + gap_size; +} + +size_t meshopt_encodeMeshlet(unsigned char* buffer, size_t buffer_size, const unsigned int* vertices, size_t vertex_count, const unsigned char* triangles, size_t triangle_count) +{ + using namespace meshopt; + + assert(triangle_count <= 256 && vertex_count <= 256); + + // 4 bits per triangle + up to three bytes of extra data + unsigned char codes[256 / 2]; + unsigned char extra[256 * 3]; + size_t codes_size = (triangle_count + 1) / 2; + size_t extra_size = encodeTriangles(codes, extra, triangles, triangle_count); + assert(extra_size <= sizeof(extra)); + + // 2 bits per vertex + up to 4 bytes of actual data + unsigned char ctrl[256 / 4]; + unsigned char data[256 * 4]; + size_t ctrl_size = (vertex_count + 3) / 4; + size_t data_size = encodeVertices(ctrl, data, vertices, vertex_count); + assert(data_size <= sizeof(data)); + + // we need to ensure that up to 16 bytes after extra+data are available for SIMD decoding + // to minimize overhead, we place fixed-size codes+control at the end of the buffer + size_t gap_size = (codes_size + ctrl_size < 16) ? 16 - (codes_size + ctrl_size) : 0; + + size_t result = codes_size + extra_size + ctrl_size + data_size + gap_size; + + if (result > buffer_size) + return 0; + + // variable-size data first + memcpy(buffer, data, data_size); + buffer += data_size; + memcpy(buffer, extra, extra_size); + buffer += extra_size; + + // gap (for accelerated decoding) separates variable-size and fixed-size data + memset(buffer, 0, gap_size); + buffer += gap_size; + + // fixed-size data last; it can be located from buffer end during decoding + memcpy(buffer, ctrl, ctrl_size); + buffer += ctrl_size; + memcpy(buffer, codes, codes_size); + buffer += codes_size; + +#if TRACE > 1 + printf("extra:"); + for (size_t i = 0; i < extra_size; ++i) + printf(" %d", extra[i]); + printf("\n"); + + unsigned int minv = ~0u; + for (size_t i = 0; i < vertex_count; ++i) + minv = minv < vertices[i] ? minv : vertices[i]; + + printf("vertices: [%d+]", minv); + for (size_t i = 0; i < vertex_count; ++i) + printf(" %d", vertices[i] - minv); + printf("\n"); +#endif + +#if TRACE + printf("stats: %d vertices, %d triangles => %d bytes (triangles: %d codes, %d extra; vertices: %d control, %d data; %d gap)\n", + int(vertex_count), int(triangle_count), int(result), + int(codes_size), int(extra_size), int(ctrl_size), int(data_size), int(gap_size)); +#endif + + return result; +} + +int meshopt_decodeMeshlet(void* vertices, size_t vertex_count, size_t vertex_size, void* triangles, size_t triangle_count, size_t triangle_size, const unsigned char* buffer, size_t buffer_size) +{ + using namespace meshopt; + + assert(triangle_count <= 256 && vertex_count <= 256); + assert(vertex_size == 4 || vertex_size == 2); + assert(triangle_size == 4 || triangle_size == 3); + + // layout must match encoding + size_t codes_size = (triangle_count + 1) / 2; + size_t ctrl_size = (vertex_count + 3) / 4; + size_t gap_size = (codes_size + ctrl_size < 16) ? 16 - (codes_size + ctrl_size) : 0; + + if (buffer_size < codes_size + ctrl_size + gap_size) + return -2; + + const unsigned char* end = buffer + buffer_size; + const unsigned char* codes = end - codes_size; + const unsigned char* ctrl = codes - ctrl_size; + const unsigned char* data = buffer; + + // gap ensures we have at least 16 bytes available after bound; this allows SIMD decoders to over-read safely + const unsigned char* bound = ctrl - gap_size; + assert(bound >= buffer && bound + 16 <= buffer + buffer_size); + +#if defined(SIMD_FALLBACK) + return (gDecodeTablesInitialized ? decodeMeshletSimd<0> : decodeMeshlet)(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, vertex_size, triangle_size); +#elif defined(SIMD_SSE) || defined(SIMD_NEON) + return decodeMeshletSimd<0>(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, vertex_size, triangle_size); +#else + return decodeMeshlet(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, vertex_size, triangle_size); +#endif +} + +int meshopt_decodeMeshletRaw(unsigned int* vertices, size_t vertex_count, unsigned int* triangles, size_t triangle_count, const unsigned char* buffer, size_t buffer_size) +{ + using namespace meshopt; + + assert(triangle_count <= 256 && vertex_count <= 256); + + // layout must match encoding + size_t codes_size = (triangle_count + 1) / 2; + size_t ctrl_size = (vertex_count + 3) / 4; + size_t gap_size = (codes_size + ctrl_size < 16) ? 16 - (codes_size + ctrl_size) : 0; + + if (buffer_size < codes_size + ctrl_size + gap_size) + return -2; + + const unsigned char* end = buffer + buffer_size; + const unsigned char* codes = end - codes_size; + const unsigned char* ctrl = codes - ctrl_size; + const unsigned char* data = buffer; + + // gap ensures we have at least 16 bytes available after bound; this allows SIMD decoders to over-read safely + const unsigned char* bound = ctrl - gap_size; + assert(bound >= buffer && bound + 16 <= buffer + buffer_size); + +#if defined(SIMD_FALLBACK) + return (gDecodeTablesInitialized ? decodeMeshletSimd<1> : decodeMeshlet)(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, 4, 4); +#elif defined(SIMD_SSE) || defined(SIMD_NEON) + return decodeMeshletSimd<1>(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, 4, 4); +#else + return decodeMeshlet(vertices, triangles, codes, ctrl, data, bound, vertex_count, triangle_count, 4, 4); +#endif +} + +#undef SIMD_SSE +#undef SIMD_NEON +#undef SIMD_FALLBACK +#undef SIMD_FLATTEN +#undef SIMD_TARGET diff --git a/3rdparty/meshoptimizer/src/meshletutils.cpp b/3rdparty/meshoptimizer/src/meshletutils.cpp new file mode 100644 index 000000000..b92332add --- /dev/null +++ b/3rdparty/meshoptimizer/src/meshletutils.cpp @@ -0,0 +1,619 @@ +// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details +#include "meshoptimizer.h" + +#include +#include +#include +#include + +// This work is based on: +// Matthaeus Chajdas. GeometryFX 1.2 - Cluster Culling. 2016 +// Jack Ritter. An Efficient Bounding Sphere. 1990 +// Thomas Larsson. Fast and Tight Fitting Bounding Spheres. 2008 +namespace meshopt +{ + +// This must be <= 256 since meshlet indices are stored as bytes +const size_t kMeshletMaxVertices = 256; + +// A reasonable limit is around 2*max_vertices or less +const size_t kMeshletMaxTriangles = 512; + +static void computeBoundingSphere(float result[4], const float* points, size_t count, size_t points_stride, const float* radii, size_t radii_stride, size_t axis_count, const unsigned int* indices = NULL) +{ + static const float axes[7][3] = { + // X, Y, Z + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1}, + + // XYZ, -XYZ, X-YZ, XY-Z; normalized to unit length + {0.57735026f, 0.57735026f, 0.57735026f}, + {-0.57735026f, 0.57735026f, 0.57735026f}, + {0.57735026f, -0.57735026f, 0.57735026f}, + {0.57735026f, 0.57735026f, -0.57735026f}, + }; + + assert(count > 0); + assert(axis_count <= sizeof(axes) / sizeof(axes[0])); + + size_t points_stride_float = points_stride / sizeof(float); + size_t radii_stride_float = radii_stride / sizeof(float); + + // find extremum points along all axes; for each axis we get a pair of points with min/max coordinates + unsigned int pmin[7], pmax[7]; + float tmin[7], tmax[7]; + + for (size_t axis = 0; axis < axis_count; ++axis) + { + pmin[axis] = pmax[axis] = 0; + tmin[axis] = FLT_MAX; + tmax[axis] = -FLT_MAX; + } + + for (size_t i = 0; i < count; ++i) + { + unsigned int v = indices ? indices[i] : unsigned(i); + const float* p = points + v * points_stride_float; + float r = radii[v * radii_stride_float]; + + for (size_t axis = 0; axis < axis_count; ++axis) + { + const float* ax = axes[axis]; + + float tp = ax[0] * p[0] + ax[1] * p[1] + ax[2] * p[2]; + float tpmin = tp - r, tpmax = tp + r; + + pmin[axis] = (tpmin < tmin[axis]) ? v : pmin[axis]; + pmax[axis] = (tpmax > tmax[axis]) ? v : pmax[axis]; + tmin[axis] = (tpmin < tmin[axis]) ? tpmin : tmin[axis]; + tmax[axis] = (tpmax > tmax[axis]) ? tpmax : tmax[axis]; + } + } + + // find the pair of points with largest distance + size_t paxis = 0; + float paxisdr = 0; + + for (size_t axis = 0; axis < axis_count; ++axis) + { + const float* p1 = points + pmin[axis] * points_stride_float; + const float* p2 = points + pmax[axis] * points_stride_float; + float r1 = radii[pmin[axis] * radii_stride_float]; + float r2 = radii[pmax[axis] * radii_stride_float]; + + float d2 = (p2[0] - p1[0]) * (p2[0] - p1[0]) + (p2[1] - p1[1]) * (p2[1] - p1[1]) + (p2[2] - p1[2]) * (p2[2] - p1[2]); + float dr = sqrtf(d2) + r1 + r2; + + if (dr > paxisdr) + { + paxisdr = dr; + paxis = axis; + } + } + + // use the longest segment as the initial sphere diameter + const float* p1 = points + pmin[paxis] * points_stride_float; + const float* p2 = points + pmax[paxis] * points_stride_float; + float r1 = radii[pmin[paxis] * radii_stride_float]; + float r2 = radii[pmax[paxis] * radii_stride_float]; + + float paxisd = sqrtf((p2[0] - p1[0]) * (p2[0] - p1[0]) + (p2[1] - p1[1]) * (p2[1] - p1[1]) + (p2[2] - p1[2]) * (p2[2] - p1[2])); + float paxisk = paxisd > 0 ? (paxisd + r2 - r1) / (2 * paxisd) : 0.f; + + float center[3] = {p1[0] + (p2[0] - p1[0]) * paxisk, p1[1] + (p2[1] - p1[1]) * paxisk, p1[2] + (p2[2] - p1[2]) * paxisk}; + float radius = paxisdr / 2; + + // iteratively adjust the sphere up until all points fit + for (size_t i = 0; i < count; ++i) + { + unsigned int v = indices ? indices[i] : unsigned(i); + const float* p = points + v * points_stride_float; + float r = radii[v * radii_stride_float]; + + float d2 = (p[0] - center[0]) * (p[0] - center[0]) + (p[1] - center[1]) * (p[1] - center[1]) + (p[2] - center[2]) * (p[2] - center[2]); + float d = sqrtf(d2); + + if (d + r > radius) + { + float k = d > 0 ? (d + r - radius) / (2 * d) : 0.f; + + center[0] += k * (p[0] - center[0]); + center[1] += k * (p[1] - center[1]); + center[2] += k * (p[2] - center[2]); + radius = (radius + d + r) / 2; + } + } + + result[0] = center[0]; + result[1] = center[1]; + result[2] = center[2]; + result[3] = radius; +} + +static meshopt_Bounds computeClusterBounds(const unsigned int* indices, size_t index_count, const unsigned int* corners, size_t corner_count, const float* vertex_positions, size_t vertex_positions_stride) +{ + size_t vertex_stride_float = vertex_positions_stride / sizeof(float); + + // compute triangle normals (.w completes plane equation) + float normals[kMeshletMaxTriangles][4]; + size_t triangles = 0; + + for (size_t i = 0; i < index_count; i += 3) + { + unsigned int a = indices[i + 0], b = indices[i + 1], c = indices[i + 2]; + + const float* p0 = vertex_positions + vertex_stride_float * a; + const float* p1 = vertex_positions + vertex_stride_float * b; + const float* p2 = vertex_positions + vertex_stride_float * c; + + float p10[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]}; + float p20[3] = {p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]}; + + float normalx = p10[1] * p20[2] - p10[2] * p20[1]; + float normaly = p10[2] * p20[0] - p10[0] * p20[2]; + float normalz = p10[0] * p20[1] - p10[1] * p20[0]; + + float area = sqrtf(normalx * normalx + normaly * normaly + normalz * normalz); + + // no need to include degenerate triangles - they will be invisible anyway + if (area == 0.f) + continue; + + normalx /= area; + normaly /= area; + normalz /= area; + + // record triangle normals; normal and corner 0 define a plane equation + normals[triangles][0] = normalx; + normals[triangles][1] = normaly; + normals[triangles][2] = normalz; + normals[triangles][3] = -(normalx * p0[0] + normaly * p0[1] + normalz * p0[2]); + triangles++; + } + + meshopt_Bounds bounds = {}; + + // degenerate cluster, no valid triangles => trivial reject (cone data is 0) + if (triangles == 0) + return bounds; + + const float rzero = 0.f; + + // compute cluster bounding sphere; we'll use the center to determine normal cone apex as well + float psphere[4] = {}; + computeBoundingSphere(psphere, vertex_positions, corner_count, vertex_positions_stride, &rzero, 0, 7, corners); + + float center[3] = {psphere[0], psphere[1], psphere[2]}; + + // treating triangle normals as points, find the bounding sphere - the sphere center determines the optimal cone axis + float nsphere[4] = {}; + computeBoundingSphere(nsphere, normals[0], triangles, sizeof(float) * 4, &rzero, 0, 3); + + float axis[3] = {nsphere[0], nsphere[1], nsphere[2]}; + float axislength = sqrtf(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]); + float invaxislength = axislength == 0.f ? 0.f : 1.f / axislength; + + axis[0] *= invaxislength; + axis[1] *= invaxislength; + axis[2] *= invaxislength; + + // compute a tight cone around all normals, mindp = cos(angle/2) + float mindp = 1.f; + + for (size_t i = 0; i < triangles; ++i) + { + float dp = normals[i][0] * axis[0] + normals[i][1] * axis[1] + normals[i][2] * axis[2]; + + mindp = (dp < mindp) ? dp : mindp; + } + + // fill bounding sphere info; note that below we can return bounds without cone information for degenerate cones + bounds.center[0] = center[0]; + bounds.center[1] = center[1]; + bounds.center[2] = center[2]; + bounds.radius = psphere[3]; + + // degenerate cluster, normal cone is larger than a hemisphere => trivial accept + // note that if mindp is positive but close to 0, the triangle intersection code below gets less stable + // we arbitrarily decide that if a normal cone is ~168 degrees wide or more, the cone isn't useful + if (mindp <= 0.1f) + { + bounds.cone_cutoff = 1; + bounds.cone_cutoff_s8 = 127; + return bounds; + } + + float maxt = 0; + + // we need to find the point on center-t*axis ray that lies in negative half-space of all triangles + for (size_t i = 0; i < triangles; ++i) + { + // dot(center-t*axis-corner, trinormal) = 0 + // dot(center-corner, trinormal) - t * dot(axis, trinormal) = 0 + float dc = center[0] * normals[i][0] + center[1] * normals[i][1] + center[2] * normals[i][2] + normals[i][3]; + float dn = axis[0] * normals[i][0] + axis[1] * normals[i][1] + axis[2] * normals[i][2]; + + // dn should be larger than mindp cutoff above + assert(dn > 0.f); + float t = dc / dn; + + maxt = (t > maxt) ? t : maxt; + } + + // cone apex should be in the negative half-space of all cluster triangles by construction + bounds.cone_apex[0] = center[0] - axis[0] * maxt; + bounds.cone_apex[1] = center[1] - axis[1] * maxt; + bounds.cone_apex[2] = center[2] - axis[2] * maxt; + + // note: this axis is the axis of the normal cone, but our test for perspective camera effectively negates the axis + bounds.cone_axis[0] = axis[0]; + bounds.cone_axis[1] = axis[1]; + bounds.cone_axis[2] = axis[2]; + + // cos(a) for normal cone is mindp; we need to add 90 degrees on both sides and invert the cone + // which gives us -cos(a+90) = -(-sin(a)) = sin(a) = sqrt(1 - cos^2(a)) + bounds.cone_cutoff = sqrtf(1 - mindp * mindp); + + // quantize axis & cutoff to 8-bit SNORM format + bounds.cone_axis_s8[0] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[0], 8)); + bounds.cone_axis_s8[1] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[1], 8)); + bounds.cone_axis_s8[2] = (signed char)(meshopt_quantizeSnorm(bounds.cone_axis[2], 8)); + + // for the 8-bit test to be conservative, we need to adjust the cutoff by measuring the max. error + float cone_axis_s8_e0 = fabsf(bounds.cone_axis_s8[0] / 127.f - bounds.cone_axis[0]); + float cone_axis_s8_e1 = fabsf(bounds.cone_axis_s8[1] / 127.f - bounds.cone_axis[1]); + float cone_axis_s8_e2 = fabsf(bounds.cone_axis_s8[2] / 127.f - bounds.cone_axis[2]); + + // note that we need to round this up instead of rounding to nearest, hence +1 + int cone_cutoff_s8 = int(127 * (bounds.cone_cutoff + cone_axis_s8_e0 + cone_axis_s8_e1 + cone_axis_s8_e2) + 1); + + bounds.cone_cutoff_s8 = (cone_cutoff_s8 > 127) ? 127 : (signed char)(cone_cutoff_s8); + + return bounds; +} + +} // namespace meshopt + +meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride) +{ + using namespace meshopt; + + assert(index_count % 3 == 0); + assert(index_count / 3 <= kMeshletMaxTriangles); + assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); + assert(vertex_positions_stride % sizeof(float) == 0); + + (void)vertex_count; + + unsigned int cache[512]; + memset(cache, -1, sizeof(cache)); + + unsigned int corners[kMeshletMaxTriangles * 3 + 1]; // +1 for branchless slot + size_t corner_count = 0; + + for (size_t i = 0; i < index_count; ++i) + { + unsigned int v = indices[i]; + assert(v < vertex_count); + + unsigned int& c = cache[v & (sizeof(cache) / sizeof(cache[0]) - 1)]; + + // branchless append if vertex isn't in cache + corners[corner_count] = v; + corner_count += (c != v); + c = v; + } + + return computeClusterBounds(indices, index_count, corners, corner_count, vertex_positions, vertex_positions_stride); +} + +meshopt_Bounds meshopt_computeMeshletBounds(const unsigned int* meshlet_vertices, const unsigned char* meshlet_triangles, size_t triangle_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride) +{ + using namespace meshopt; + + assert(triangle_count <= kMeshletMaxTriangles); + assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); + assert(vertex_positions_stride % sizeof(float) == 0); + + (void)vertex_count; + + unsigned int indices[kMeshletMaxTriangles * 3]; + size_t corner_count = 0; + + for (size_t i = 0; i < triangle_count * 3; ++i) + { + unsigned char t = meshlet_triangles[i]; + unsigned int index = meshlet_vertices[t]; + assert(index < vertex_count); + + indices[i] = index; + + // meshlet_vertices[] slice should only contain vertices used by triangle indices, which is the case for any well formed meshlet + corner_count = t >= corner_count ? t + 1 : corner_count; + } + + return computeClusterBounds(indices, triangle_count * 3, meshlet_vertices, corner_count, vertex_positions, vertex_positions_stride); +} + +meshopt_Bounds meshopt_computeSphereBounds(const float* positions, size_t count, size_t positions_stride, const float* radii, size_t radii_stride) +{ + using namespace meshopt; + + assert(positions_stride >= 12 && positions_stride <= 256); + assert(positions_stride % sizeof(float) == 0); + assert((radii_stride >= 4 && radii_stride <= 256) || radii == NULL); + assert(radii_stride % sizeof(float) == 0); + + meshopt_Bounds bounds = {}; + + if (count == 0) + return bounds; + + const float rzero = 0.f; + + float psphere[4] = {}; + computeBoundingSphere(psphere, positions, count, positions_stride, radii ? radii : &rzero, radii ? radii_stride : 0, 7); + + bounds.center[0] = psphere[0]; + bounds.center[1] = psphere[1]; + bounds.center[2] = psphere[2]; + bounds.radius = psphere[3]; + + return bounds; +} + +void meshopt_optimizeMeshletLevel(unsigned int* meshlet_vertices, size_t vertex_count, unsigned char* meshlet_triangles, size_t triangle_count, int level) +{ + using namespace meshopt; + + assert(triangle_count <= kMeshletMaxTriangles); + assert(vertex_count <= kMeshletMaxVertices); + assert(level >= 0 && level <= 9); + + unsigned char* indices = meshlet_triangles; + unsigned int* vertices = meshlet_vertices; + + // cache tracks vertex timestamps (corresponding to triangle index! all 3 vertices are added at the same time and never removed) + unsigned char cache[kMeshletMaxVertices]; + memset(cache, 0, vertex_count); + + // note that we start from a value that means all vertices aren't in cache + unsigned char cache_last = 128; + const unsigned char cache_cutoff = 3; // 3 triangles = ~5..9 vertices depending on reuse + + // vertex valence is used to prioritize triangles for level>0 + // note: we use 8-bit counters for performance; for outlier vertices the valence is incorrect but that just affects the heuristic + unsigned char valence[kMeshletMaxVertices]; + memset(valence, 0, vertex_count); + + for (size_t i = 0; i < triangle_count; ++i) + { + unsigned char a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2]; + assert(a < vertex_count && b < vertex_count && c < vertex_count); + + valence[a]++; + valence[b]++; + valence[c]++; + } + + for (size_t i = 0; i < triangle_count; ++i) + { + int next = -1; + int next_score = -1; + int edges = 0; + + for (size_t j = i; j < triangle_count; ++j) + { + unsigned char a = indices[j * 3 + 0], b = indices[j * 3 + 1], c = indices[j * 3 + 2]; + assert(a < vertex_count && b < vertex_count && c < vertex_count); + + // compute cache distance using unsigned 8-bit subtraction, so cache timestamp overflow is handled gracefully + unsigned char ad = (unsigned char)(cache_last - cache[a]); + unsigned char bd = (unsigned char)(cache_last - cache[b]); + unsigned char cd = (unsigned char)(cache_last - cache[c]); + + int match = (ad < cache_cutoff) + (bd < cache_cutoff) + (cd < cache_cutoff); + + if (level) + { + // prefer low minimum valence + int vmin = valence[a] < valence[b] ? valence[a] : valence[b]; + vmin = valence[c] < vmin ? valence[c] : vmin; + + // prefer vertices with smaller cache distance and valence to improve traversal locality + int score = match * 1024 + (1023 - ad - bd - cd); + score = score * 256 + (255 - vmin); + + next = (score > next_score) ? int(j) : next; + next_score = (score > next_score) ? score : next_score; + + // terminate after finding enough edge matches + if (match >= 2 && ++edges >= level) + break; + } + else + { + int score = match; + + next = (score > next_score) ? int(j) : next; + next_score = (score > next_score) ? score : next_score; + + // settle for a first edge match, which makes the function ~linear in practice + if (match >= 2) + break; + } + } + + assert(next >= 0); + + unsigned char a = indices[next * 3 + 0], b = indices[next * 3 + 1], c = indices[next * 3 + 2]; + + // shift triangles before the next one forward so that we always keep an ordered partition + // note: this could have swapped triangles [i] and [next] but that distorts the order and may skew the output sequence + memmove(indices + (i + 1) * 3, indices + i * 3, (next - i) * 3 * sizeof(unsigned char)); + + indices[i * 3 + 0] = a; + indices[i * 3 + 1] = b; + indices[i * 3 + 2] = c; + + // cache timestamp is the same between all vertices of each triangle to reduce overflow + cache_last++; + cache[a] = cache_last; + cache[b] = cache_last; + cache[c] = cache_last; + + // update vertex valences for scoring heuristic + valence[a]--; + valence[b]--; + valence[c]--; + } + + // rotate triangles to maximize compressibility; only done at level >= 1 for compatibility + if (level >= 1) + { + memset(cache, 0, vertex_count); + + for (size_t i = 0; i < triangle_count; ++i) + { + unsigned char a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2]; + + // if only the middle vertex has been used, rotate triangle to ensure new vertices are always sequential + if (!cache[a] && cache[b] && !cache[c]) + { + // abc -> bca + unsigned char t = a; + a = b, b = c, c = t; + } + else if (!cache[a] && !cache[b] && !cache[c]) + { + // out of three edges, the edge ab can not be reused by subsequent triangles in some encodings + // if subsequent triangles don't share edges ca or bc, we can rotate the triangle to fix this + bool needab = false, needbc = false, needca = false; + + for (size_t j = i + 1; j < triangle_count && j <= i + cache_cutoff; ++j) + { + unsigned char oa = indices[j * 3 + 0], ob = indices[j * 3 + 1], oc = indices[j * 3 + 2]; + + // note: edge comparisons are reversed as reused edges are flipped + needab |= (oa == b && ob == a) || (ob == b && oc == a) || (oc == b && oa == a); + needbc |= (oa == c && ob == b) || (ob == c && oc == b) || (oc == c && oa == b); + needca |= (oa == a && ob == c) || (ob == a && oc == c) || (oc == a && oa == c); + } + + if (needab && !needbc) + { + // abc -> bca + unsigned char t = a; + a = b, b = c, c = t; + } + else if (needab && !needca) + { + // abc -> cab + unsigned char t = c; + c = b, b = a, a = t; + } + else if (!needab && needbc && !needca) + { + // abc -> cab: put bc at most recent FIFO position + unsigned char t = c; + c = b, b = a, a = t; + } + } + + indices[i * 3 + 0] = a, indices[i * 3 + 1] = b, indices[i * 3 + 2] = c; + + cache[a] = cache[b] = cache[c] = 1; + } + } + + // reorder meshlet vertices for access locality assuming index buffer is scanned sequentially + unsigned int order[kMeshletMaxVertices]; + + short remap[kMeshletMaxVertices]; + memset(remap, -1, vertex_count * sizeof(short)); + + size_t vertex_offset = 0; + + for (size_t i = 0; i < triangle_count * 3; ++i) + { + short& r = remap[indices[i]]; + + if (r < 0) + { + r = short(vertex_offset); + order[vertex_offset] = vertices[indices[i]]; + vertex_offset++; + } + + indices[i] = (unsigned char)r; + } + + assert(vertex_offset <= vertex_count); + memcpy(vertices, order, vertex_offset * sizeof(unsigned int)); +} + +void meshopt_optimizeMeshlet(unsigned int* meshlet_vertices, unsigned char* meshlet_triangles, size_t triangle_count, size_t vertex_count) +{ + meshopt_optimizeMeshletLevel(meshlet_vertices, vertex_count, meshlet_triangles, triangle_count, 0); +} + +size_t meshopt_extractMeshletIndices(unsigned int* vertices, unsigned char* triangles, const unsigned int* indices, size_t index_count) +{ + using namespace meshopt; + + assert(index_count % 3 == 0); + assert(index_count / 3 <= kMeshletMaxTriangles); + + size_t unique = 0; + + // direct mapped cache for fast lookups based on low index bits; inspired by vk_lod_clusters from NVIDIA + short cache[1024]; + memset(cache, -1, sizeof(cache)); + + for (size_t i = 0; i < index_count; ++i) + { + unsigned int v = indices[i]; + unsigned int key = v & (sizeof(cache) / sizeof(cache[0]) - 1); + short c = cache[key]; + + // fast path: vertex has been seen before + if (c >= 0 && vertices[c] == v) + { + triangles[i] = (unsigned char)c; + continue; + } + + // fast path: vertex has never been seen before + if (c < 0) + { + assert(unique < kMeshletMaxVertices); + cache[key] = short(unique); + triangles[i] = (unsigned char)unique; + vertices[unique++] = v; + continue; + } + + // slow path: collision with a different vertex, so we need to look through all vertices + int pos = -1; + for (size_t j = 0; j < unique; ++j) + if (vertices[j] == v) + { + pos = int(j); + break; + } + + if (pos < 0) + { + assert(unique < kMeshletMaxVertices); + pos = int(unique); + vertices[unique++] = v; + } + + cache[key] = short(pos); + triangles[i] = (unsigned char)pos; + } + + assert(unique <= kMeshletMaxVertices); + return unique; +} diff --git a/3rdparty/meshoptimizer/src/meshoptimizer.h b/3rdparty/meshoptimizer/src/meshoptimizer.h index d3faf869d..47e782ae4 100644 --- a/3rdparty/meshoptimizer/src/meshoptimizer.h +++ b/3rdparty/meshoptimizer/src/meshoptimizer.h @@ -1,7 +1,7 @@ /** - * meshoptimizer - version 0.25 + * meshoptimizer - version 1.2 * - * Copyright (C) 2016-2025, by Arseny Kapoulkine (arseny.kapoulkine@gmail.com) + * Copyright (C) 2016-2026, by Arseny Kapoulkine (arseny.kapoulkine@gmail.com) * Report bugs and download new versions at https://github.com/zeux/meshoptimizer * * This library is distributed under the MIT License. See notice at the end of this file. @@ -12,7 +12,7 @@ #include /* Version macro; major * 1000 + minor * 10 + patch */ -#define MESHOPTIMIZER_VERSION 250 /* 0.25 */ +#define MESHOPTIMIZER_VERSION 1020 /* 1.2 */ /* If no API is defined, assume default */ #ifndef MESHOPTIMIZER_API @@ -103,6 +103,27 @@ MESHOPTIMIZER_API void meshopt_remapVertexBuffer(void* destination, const void* */ MESHOPTIMIZER_API void meshopt_remapIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const unsigned int* remap); +/** + * Experimental: Filter out redundant triangles from the index buffer and return the number of remaining indices + * Triangles are considered redundant if they are degenerate (two vertices have the same vertex key) or duplicate (matching triangle was present earlier). + * First vertex_size bytes of every vertex are compared for equality; typically vertex_size should be set to the size of the position attribute. + * Note that duplicate triangles with opposite windings are preserved, as they may be needed for double-sided rendering. + * + * destination must contain enough space for the resulting index buffer (index_count elements) + */ +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_filterIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride); + +/** + * Experimental: Filter out redundant triangles from the index buffer and return the number of remaining indices + * Triangles are considered redundant if they are degenerate (two vertices have the same vertex key) or duplicate (matching triangle was present earlier). + * All bytes in specified streams are compared for equality; streams should include attributes relevant for position transform (e.g. bone influences). + * Note that duplicate triangles with opposite windings are preserved, as they may be needed for double-sided rendering. + * + * destination must contain enough space for the resulting index buffer (index_count elements) + * stream_count must be <= 16 + */ +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_filterIndexBufferMulti(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, const struct meshopt_Stream* streams, size_t stream_count); + /** * Generate index buffer that can be used for more efficient rendering when only a subset of the vertex attributes is necessary * All vertices that are binary equivalent (wrt first vertex_size bytes) map to the first vertex in the original vertex buffer. @@ -179,7 +200,7 @@ MESHOPTIMIZER_API size_t meshopt_generateProvokingIndexBuffer(unsigned int* dest /** * Vertex transform cache optimizer * Reorders indices to reduce the number of GPU vertex shader invocations - * If index buffer contains multiple ranges for multiple draw calls, this functions needs to be called on each range individually. + * If index buffer contains multiple ranges for multiple draw calls, this function needs to be called on each range individually. * * destination must contain enough space for the resulting index buffer (index_count elements) */ @@ -198,7 +219,7 @@ MESHOPTIMIZER_API void meshopt_optimizeVertexCacheStrip(unsigned int* destinatio * Vertex transform cache optimizer for FIFO caches * Reorders indices to reduce the number of GPU vertex shader invocations * Generally takes ~3x less time to optimize meshes but produces inferior results compared to meshopt_optimizeVertexCache - * If index buffer contains multiple ranges for multiple draw calls, this functions needs to be called on each range individually. + * If index buffer contains multiple ranges for multiple draw calls, this function needs to be called on each range individually. * * destination must contain enough space for the resulting index buffer (index_count elements) * cache_size should be less than the actual GPU cache size to avoid cache thrashing @@ -208,7 +229,7 @@ MESHOPTIMIZER_API void meshopt_optimizeVertexCacheFifo(unsigned int* destination /** * Overdraw optimizer * Reorders indices to reduce the number of GPU vertex shader invocations and the pixel overdraw - * If index buffer contains multiple ranges for multiple draw calls, this functions needs to be called on each range individually. + * If index buffer contains multiple ranges for multiple draw calls, this function needs to be called on each range individually. * * destination must contain enough space for the resulting index buffer (index_count elements) * indices must contain index data that is the result of meshopt_optimizeVertexCache (*not* the original mesh indices!) @@ -221,7 +242,7 @@ MESHOPTIMIZER_API void meshopt_optimizeOverdraw(unsigned int* destination, const * Vertex fetch cache optimizer * Reorders vertices and changes indices to reduce the amount of GPU memory fetches during vertex processing * Returns the number of unique vertices, which is the same as input vertex count unless some vertices are unused - * This functions works for a single vertex stream; for multiple vertex streams, use meshopt_optimizeVertexFetchRemap + meshopt_remapVertexBuffer for each stream. + * This function works for a single vertex stream; for multiple vertex streams, use meshopt_optimizeVertexFetchRemap + meshopt_remapVertexBuffer for each stream. * * destination must contain enough space for the resulting vertex buffer (vertex_count elements) * indices is used both as an input and as an output index buffer @@ -252,6 +273,7 @@ MESHOPTIMIZER_API size_t meshopt_encodeIndexBufferBound(size_t index_count, size /** * Set index encoder format version (defaults to 1) + * This function is not thread safe and must not be called concurrently with meshopt_encodeIndexBuffer/meshopt_encodeIndexSequence. * * version must specify the data format version to encode; valid values are 0 (decodable by all library versions) and 1 (decodable by 0.14+) */ @@ -295,6 +317,38 @@ MESHOPTIMIZER_API size_t meshopt_encodeIndexSequenceBound(size_t index_count, si */ MESHOPTIMIZER_API int meshopt_decodeIndexSequence(void* destination, size_t index_count, size_t index_size, const unsigned char* buffer, size_t buffer_size); +/** + * Meshlet encoder + * Encodes meshlet data into an array of bytes that is generally smaller and compresses better compared to original. + * Returns encoded data size on success, 0 on error; the only error condition is if buffer doesn't have enough space + * This function encodes a single meshlet; when encoding multiple meshlets, additional headers may be necessary to store vertex/triangle count and encoded size. + * For maximum efficiency the meshlet being encoded should be optimized using meshopt_optimizeMeshletLevel with level 1+ (3 recommended); additionally, vertex reference data should be optimized for locality (fetch). + * + * buffer must contain enough space for the encoded meshlet (use meshopt_encodeMeshletBound to compute worst case size) + * vertices may be NULL, in which case vertex_count must be 0 and only triangle data is encoded + * vertex_count and triangle_count must be <= 256. + */ +MESHOPTIMIZER_API size_t meshopt_encodeMeshlet(unsigned char* buffer, size_t buffer_size, const unsigned int* vertices, size_t vertex_count, const unsigned char* triangles, size_t triangle_count); +MESHOPTIMIZER_API size_t meshopt_encodeMeshletBound(size_t max_vertices, size_t max_triangles); + +/** + * Meshlet decoder + * Decodes meshlet data from an array of bytes generated by meshopt_encodeMeshlet + * Returns 0 if decoding was successful, and an error code otherwise + * The decoder is safe to use for untrusted input, but it may produce garbage data. + * + * vertices must contain enough space for the resulting vertex data, aligned to 4 bytes (align(vertex_count * vertex_size, 4) bytes) + * vertex_size must be 2 (16-bit vertex references) or 4 (32-bit vertex references) + * triangles must contain enough space for the resulting triangle data, aligned to 4 bytes (align(triangle_count * triangle_size, 4) bytes) + * triangle_size must be 3 (8-bit triangle indices) or 4 (32-bit packed triangles, stored as (a) | (b << 8) | (c << 16)) + * vertex_count, triangle_count match those used during encoding exactly; buffer_size must be equal to the encoded size returned by meshopt_encodeMeshlet. + * vertices may be NULL, in which case vertex_count must be 0 and the meshlet must contain just triangle data + * + * When using "raw" decoding (meshopt_decodeMeshletRaw), both vertices and triangles should have available space further aligned to 16 bytes for efficient SIMD decoding. + */ +MESHOPTIMIZER_API int meshopt_decodeMeshlet(void* vertices, size_t vertex_count, size_t vertex_size, void* triangles, size_t triangle_count, size_t triangle_size, const unsigned char* buffer, size_t buffer_size); +MESHOPTIMIZER_API int meshopt_decodeMeshletRaw(unsigned int* vertices, size_t vertex_count, unsigned int* triangles, size_t triangle_count, const unsigned char* buffer, size_t buffer_size); + /** * Vertex buffer encoder * Encodes vertex data into an array of bytes that is generally smaller and compresses better compared to original. @@ -324,6 +378,7 @@ MESHOPTIMIZER_API size_t meshopt_encodeVertexBufferLevel(unsigned char* buffer, /** * Set vertex encoder format version (defaults to 1) + * This function is not thread safe and must not be called concurrently with meshopt_encodeVertexBuffer/meshopt_encodeVertexBufferLevel. * * version must specify the data format version to encode; valid values are 0 (decodable by all library versions) and 1 (decodable by 0.23+) */ @@ -351,16 +406,16 @@ MESHOPTIMIZER_API int meshopt_decodeVertexVersion(const unsigned char* buffer, s * Vertex buffer filters * These functions can be used to filter output of meshopt_decodeVertexBuffer in-place. * - * meshopt_decodeFilterOct decodes octahedral encoding of a unit vector with K-bit (K <= 16) signed X/Y as an input; Z must store 1.0f. + * meshopt_decodeFilterOct decodes octahedral encoding of a unit vector with K-bit signed X/Y as an input; Z must store 1.0f. * Each component is stored as an 8-bit or 16-bit normalized integer; stride must be equal to 4 or 8. W is preserved as is. * - * meshopt_decodeFilterQuat decodes 3-component quaternion encoding with K-bit (4 <= K <= 16) component encoding and a 2-bit component index indicating which component to reconstruct. - * Each component is stored as an 16-bit integer; stride must be equal to 8. + * meshopt_decodeFilterQuat decodes 3-component quaternion encoding with K-bit component encoding and a 2-bit component index indicating which component to reconstruct. + * Each component is stored as a 16-bit integer; stride must be equal to 8. * * meshopt_decodeFilterExp decodes exponential encoding of floating-point data with 8-bit exponent and 24-bit integer mantissa as 2^E*M. * Each 32-bit component is decoded in isolation; stride must be divisible by 4. * - * meshopt_decodeFilterColor decodes YCoCg (+A) color encoding where RGB is converted to YCoCg space with variable bit quantization. + * meshopt_decodeFilterColor decodes RGBA colors from YCoCg (+A) color encoding where RGB is converted to YCoCg space with K-bit component encoding, and A is stored using K-1 bits. * Each component is stored as an 8-bit or 16-bit normalized integer; stride must be equal to 4 or 8. */ MESHOPTIMIZER_API void meshopt_decodeFilterOct(void* buffer, size_t count, size_t stride); @@ -372,20 +427,20 @@ MESHOPTIMIZER_API void meshopt_decodeFilterColor(void* buffer, size_t count, siz * Vertex buffer filter encoders * These functions can be used to encode data in a format that meshopt_decodeFilter can decode * - * meshopt_encodeFilterOct encodes unit vectors with K-bit (K <= 16) signed X/Y as an output. - * Each component is stored as an 8-bit or 16-bit normalized integer; stride must be equal to 4 or 8. W is preserved as is. + * meshopt_encodeFilterOct encodes unit vectors with K-bit (2 <= K <= 16) signed X/Y as an output. + * Each component is stored as an 8-bit or 16-bit normalized integer; stride must be equal to 4 (K <= 8) or 8 (K <= 16). Z will store 1.0f, W is preserved as is. * Input data must contain 4 floats for every vector (count*4 total). * * meshopt_encodeFilterQuat encodes unit quaternions with K-bit (4 <= K <= 16) component encoding. - * Each component is stored as an 16-bit integer; stride must be equal to 8. + * Each component is stored as a 16-bit integer; stride must be equal to 8. * Input data must contain 4 floats for every quaternion (count*4 total). * * meshopt_encodeFilterExp encodes arbitrary (finite) floating-point data with 8-bit exponent and K-bit integer mantissa (1 <= K <= 24). - * Exponent can be shared between all components of a given vector as defined by stride or all values of a given component; stride must be divisible by 4. + * Exponent can be shared between all components of a given vector as defined by stride or all values of a given component; stride must be divisible by 4 (and <= 256). * Input data must contain stride/4 floats for every vector (count*stride/4 total). * - * meshopt_encodeFilterColor encodes RGBA color data by converting RGB to YCoCg color space with variable bit quantization. - * Each component is stored as an 8-bit or 16-bit integer; stride must be equal to 4 or 8. + * meshopt_encodeFilterColor encodes RGBA color data by converting RGB to YCoCg color space with K-bit (2 <= K <= 16) component encoding; A is stored using K-1 bits. + * Each component is stored as an 8-bit or 16-bit integer; stride must be equal to 4 (K <= 8) or 8 (K <= 16). * Input data must contain 4 floats for every color (count*4 total). */ enum meshopt_EncodeExpMode @@ -418,14 +473,18 @@ enum meshopt_SimplifyErrorAbsolute = 1 << 2, /* Remove disconnected parts of the mesh during simplification incrementally, regardless of the topological restrictions inside components. */ meshopt_SimplifyPrune = 1 << 3, - /* Produce more regular triangle sizes and shapes during simplification, at some cost to geometric quality. */ + /* Produce more regular triangle sizes and shapes during simplification, at some cost to geometric and attribute quality. */ meshopt_SimplifyRegularize = 1 << 4, /* Experimental: Allow collapses across attribute discontinuities, except for vertices that are tagged with meshopt_SimplifyVertex_Protect in vertex_lock. */ meshopt_SimplifyPermissive = 1 << 5, + /* Produce more regular triangle sizes and shapes during simplification, at a small cost to geometric and attribute quality. */ + meshopt_SimplifyRegularizeLight = 1 << 6, + /* Experimental: Try to preserve fold lines between opposite-facing triangles, at a small performance cost. */ + meshopt_SimplifyPreserveFolds = 1 << 7, }; /** - * Experimental: Simplification vertex flags/locks, for use in `vertex_lock` arrays in simplification APIs + * Simplification vertex flags/locks, for use in `vertex_lock` arrays in simplification APIs */ enum { @@ -433,6 +492,8 @@ enum meshopt_SimplifyVertex_Lock = 1 << 0, /* Protect attribute discontinuity at this vertex; must be used together with meshopt_SimplifyPermissive option. */ meshopt_SimplifyVertex_Protect = 1 << 1, + /* Increase priority for this vertex, making it more likely that it's preserved during simplification. */ + meshopt_SimplifyVertex_Priority = 1 << 2, }; /** @@ -447,9 +508,9 @@ enum * * destination must contain enough space for the target index buffer, worst case is index_count elements (*not* target_index_count)! * vertex_positions should have float3 position in the first 12 bytes of each vertex - * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation; value range [0..1] + * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation (unless absolute error option is used) * options must be a bitmask composed of meshopt_SimplifyX options; 0 is a safe default - * result_error can be NULL; when it's not NULL, it will contain the resulting (relative) error after simplification + * result_error can be NULL; when it's not NULL, it will contain the resulting (relative/absolute) error after simplification */ MESHOPTIMIZER_API size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error, unsigned int options, float* result_error); @@ -470,10 +531,10 @@ MESHOPTIMIZER_API size_t meshopt_simplify(unsigned int* destination, const unsig * vertex_attributes should have attribute_count floats for each vertex * attribute_weights should have attribute_count floats in total; the weights determine relative priority of attributes between each other and wrt position * attribute_count must be <= 32 - * vertex_lock can be NULL; when it's not NULL, it should have a value for each vertex; 1 denotes vertices that can't be moved - * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation; value range [0..1] + * vertex_lock can be NULL; when it's not NULL, it should have a value for each vertex composed of meshopt_SimplifyVertex_* flags + * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation (unless absolute error option is used) * options must be a bitmask composed of meshopt_SimplifyX options; 0 is a safe default - * result_error can be NULL; when it's not NULL, it will contain the resulting (relative) error after simplification + * result_error can be NULL; when it's not NULL, it will contain the resulting (relative/absolute) error after simplification */ MESHOPTIMIZER_API size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_attributes, size_t vertex_attributes_stride, const float* attribute_weights, size_t attribute_count, const unsigned char* vertex_lock, size_t target_index_count, float target_error, unsigned int options, float* result_error); @@ -493,10 +554,10 @@ MESHOPTIMIZER_API size_t meshopt_simplifyWithAttributes(unsigned int* destinatio * vertex_attributes should have attribute_count floats for each vertex * attribute_weights should have attribute_count floats in total; the weights determine relative priority of attributes between each other and wrt position * attribute_count must be <= 32 - * vertex_lock can be NULL; when it's not NULL, it should have a value for each vertex; 1 denotes vertices that can't be moved - * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation; value range [0..1] + * vertex_lock can be NULL; when it's not NULL, it should have a value for each vertex composed of meshopt_SimplifyVertex_* flags + * target_error represents the error relative to mesh extents that can be tolerated, e.g. 0.01 = 1% deformation (unless absolute error option is used) * options must be a bitmask composed of meshopt_SimplifyX options; 0 is a safe default - * result_error can be NULL; when it's not NULL, it will contain the resulting (relative) error after simplification + * result_error can be NULL; when it's not NULL, it will contain the resulting (relative/absolute) error after simplification */ MESHOPTIMIZER_API size_t meshopt_simplifyWithUpdate(unsigned int* indices, size_t index_count, float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float* vertex_attributes, size_t vertex_attributes_stride, const float* attribute_weights, size_t attribute_count, const unsigned char* vertex_lock, size_t target_index_count, float target_error, unsigned int options, float* result_error); @@ -645,7 +706,7 @@ struct meshopt_Meshlet unsigned int vertex_offset; unsigned int triangle_offset; - /* number of vertices and triangles used in the meshlet; data is stored in consecutive range defined by offset and count */ + /* number of vertices and triangles used in the meshlet; data is stored in consecutive range [offset..offset+count) for vertices and [offset..offset+count*3) for triangles */ unsigned int vertex_count; unsigned int triangle_count; }; @@ -706,6 +767,15 @@ MESHOPTIMIZER_API size_t meshopt_buildMeshletsSpatial(struct meshopt_Meshlet* me */ MESHOPTIMIZER_API void meshopt_optimizeMeshlet(unsigned int* meshlet_vertices, unsigned char* meshlet_triangles, size_t triangle_count, size_t vertex_count); +/** + * Meshlet optimizer + * Reorders meshlet vertices and triangles to maximize locality, with higher levels resulting in smaller compressed size at the cost of optimization time. + * At level 0 the result is equivalent to meshopt_optimizeMeshlet; levels >= 1 may rotate triangle corners to improve compression (which can change provoking vertex and affect OMM data). + * + * level should be in the range [0, 9] with 0 equivalent to meshopt_optimizeMeshlet and 9 being the slowest; the sweet spot for compression ratio is around 3 + */ +MESHOPTIMIZER_API void meshopt_optimizeMeshletLevel(unsigned int* meshlet_vertices, size_t vertex_count, unsigned char* meshlet_triangles, size_t triangle_count, int level); + struct meshopt_Bounds { /* bounding sphere, useful for frustum and occlusion culling */ @@ -743,6 +813,7 @@ struct meshopt_Bounds * * vertex_positions should have float3 position in the first 12 bytes of each vertex * vertex_count should specify the number of vertices in the entire mesh, not cluster or meshlet + * indices should have at most 256 unique vertex indices * index_count/3 and triangle_count must not exceed implementation limits (<= 512) */ MESHOPTIMIZER_API struct meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride); @@ -757,6 +828,17 @@ MESHOPTIMIZER_API struct meshopt_Bounds meshopt_computeMeshletBounds(const unsig */ MESHOPTIMIZER_API struct meshopt_Bounds meshopt_computeSphereBounds(const float* positions, size_t count, size_t positions_stride, const float* radii, size_t radii_stride); +/** + * Extract meshlet-local vertex and triangle indices from absolute cluster indices. + * Fills triangles[] and vertices[] such that vertices[triangles[i]] == indices[i], and returns the number of unique vertices. + * + * vertices must contain enough space for the resulting references (up to 256 elements) + * triangles must contain enough space for local indices (index_count elements) + * indices should have at most 256 unique vertex indices + * index_count/3 must not exceed implementation limits (<= 512) + */ +MESHOPTIMIZER_API size_t meshopt_extractMeshletIndices(unsigned int* vertices, unsigned char* triangles, const unsigned int* indices, size_t index_count); + /** * Cluster partitioner * Partitions clusters into groups of similar size, prioritizing grouping clusters that share vertices or are close to each other. @@ -800,6 +882,118 @@ MESHOPTIMIZER_API void meshopt_spatialSortTriangles(unsigned int* destination, c */ MESHOPTIMIZER_API void meshopt_spatialClusterPoints(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t cluster_size); +/** + * Experimental: Opacity micromap generator (measure) + * Computes a subdivision level for each input triangle, as well as deduplicating the triangles that reference the same UVs to reduce rasterization requests. + * Returns the number of OMM entries. + * + * levels and sources must contain enough space for the worst case output (index_count/3 elements, one per resulting OMM entry) + * levels[i] will contain the subdivision level for entry i, with the total number of entries returned by the function; each entry should be rasterized from triangle index sources[i] + * omm_indices must contain enough space for the resulting OMM indices (index_count/3 elements, one per triangle) + * vertex_uvs should have float2 texture coordinate in the first 8 bytes of each vertex + * max_level specifies the maximum subdivision level (0..12) + * target_edge can be 0; when >0, triangle subdivision is adaptive and targets target_edge^2 texel area + */ +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_opacityMapMeasure(unsigned char* levels, unsigned int* sources, int* omm_indices, const unsigned int* indices, size_t index_count, const float* vertex_uvs, size_t vertex_count, size_t vertex_uvs_stride, unsigned int texture_width, unsigned int texture_height, int max_level, float target_edge); + +/** + * Experimental: Opacity micromap generator (rasterize) + * Rasterizes opacity state for a single triangle entry by sampling the alpha texture, using bilinear filtering and 0.5 alpha cutoff. + * + * result should contain enough space for the output opacity data (which can be computed using meshopt_opacityMapEntrySize) + * level specifies the subdivision level (0..12) + * states should be 2 for 2-state format (opaque/transparent) and 4 for 4-state format (opaque/transparent/unknown) + * uv0/uv1/uv2 should refer to a float2 texture coordinate for each triangle corner; note that micromap data is sensitive to the corner order + * texture_data should point to the alpha channel of the first pixel, encoded as UNORM8 + * texture_stride specifies the distance in bytes between consecutive pixels, e.g. 4 for RGBA input + * texture_pitch specifies the distance in bytes between consecutive rows, e.g. 4*texture_width for tightly packed RGBA input + */ +MESHOPTIMIZER_EXPERIMENTAL void meshopt_opacityMapRasterize(unsigned char* result, int level, int states, const float* uv0, const float* uv1, const float* uv2, const unsigned char* texture_data, size_t texture_stride, size_t texture_pitch, unsigned int texture_width, unsigned int texture_height); +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_opacityMapEntrySize(int level, int states); + +/** + * Experimental: Opacity micromap generator (compact) + * Compacts and deduplicates opacity data, merging identical micromap entries and replacing micromap states with special indices (-4..-1) when possible. + * Returns the number of OMM entries after compaction; the data array should be trimmed using the last offset/size. + * + * data should contain opacity data for all input/output entries + * levels should contain subdivision levels for all input/output entries + * offsets should contain offset into data[] for each entry + * levels[i] and offsets[i] will be updated with post-compaction level/offset for entry i, with the total number of entries returned by the function + * omm_indices should contain indices into the original OMM data, and will be updated with a new index or a special index (-4..-1) when possible + * states should be 2 for 2-state format (opaque/transparent) and 4 for 4-state format (opaque/transparent/unknown) + */ +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_opacityMapCompact(unsigned char* data, size_t data_size, unsigned char* levels, unsigned int* offsets, size_t omm_count, int* omm_indices, size_t triangle_count, int states); + +/** + * Tangent generation options + */ +enum +{ + /* Produce tangents compatible with MikkTSpace (same weighting and fallbacks) at the cost of reduced quality. Not recommended unless normal maps are baked. */ + meshopt_TangentCompatible = 1 << 0, + /* Experimental: For vertices only connected to degenerate triangles, output zero tangents instead of an arbitrary fallback. */ + meshopt_TangentZeroFallback = 1 << 1, +}; + +/** + * Experimental: Tangent space generator + * Computes per-corner tangent vectors; for each corner, computes normalized tangent vector (xyz) and orientation (w, +/-1). + * Bitangent can be reconstructed via cross(normal, tangent.xyz) * tangent.w. + * To apply tangents to the mesh, either deindex and reindex it with the tangent stream, or copy tangents to existing vertex data while duplicating + * vertices with different tangent vectors (e.g. on UV mirror seams). + * Input can be indexed or unindexed (indices=NULL); this does not affect the resulting tangents, but indexed inputs are faster to process. + * + * result must contain enough space for the output tangent data (index_count*4 elements) + * indices can be NULL if the input is unindexed + * vertex_positions should have float3 position in the first 12 bytes of each vertex + * vertex_normals should have unit float3 normal in the first 12 bytes of each vertex + * vertex_uvs should have float2 texture coordinate in the first 8 bytes of each vertex + */ +MESHOPTIMIZER_EXPERIMENTAL void meshopt_generateTangents(float* result, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, const float* vertex_uvs, size_t vertex_uvs_stride, unsigned int options); + +/** + * Experimental: Normal generator + * Computes per-corner normal vectors; normals are averaged across soft edges, using crease angle to determine which edges are hard. + * To apply normals to the mesh, either deindex and reindex it with the normal stream, or copy normals to existing vertex data while duplicating + * vertices with different normal vectors (on normal seams). + * Input can be indexed or unindexed (indices=NULL); this does not affect the resulting normals, but indexed inputs are faster to process. + * + * result must contain enough space for the output normal data (index_count*3 elements) + * indices can be NULL if the input is unindexed + * vertex_positions should have float3 position in the first 12 bytes of each vertex + * crease_angle specifies the angle in radians below which normals are averaged across edges + * smoothing specifies the amount of smoothing to apply to resulting normals as a post-process; larger smoothing values take longer, with the recommended range being [0..5]. + */ +MESHOPTIMIZER_EXPERIMENTAL void meshopt_generateNormals(float* result, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float crease_angle, float smoothing); + +/** + * Experimental: Remeshing flags + */ +enum +{ + /* Thicken thin geometry sheets, producing unique vertex positions for each side of thin geometry. Currently non-functional. */ + meshopt_RemeshThicken = 1 << 0, + /* Produce a two-sided shell that wraps around surfaces of the original mesh, instead of a solid mesh. */ + meshopt_RemeshShell = 1 << 1, + /* Compute optimal output positions that approximate the original surface as closely as possible. */ + meshopt_RemeshSolve = 1 << 2, +}; + +/** + * Experimental: Voxel remesher + * Generates a new mesh that approximates or refines the original mesh, attempting to preserve mesh appearance at a given voxel resolution. + * The original topology is not preserved; features that are closer than the voxel size may be merged, and small gaps in the original mesh may be closed. + * Returns the number of triangles in the new mesh, with destination containing a vertex position for each triangle corner. + * + * destination can be NULL; when it's not NULL, it must contain enough space for the resulting triangle buffer (max_triangle_count * 3 vertices, 3 floats per vertex) + * max_triangle_count is the number of triangles that can be written; when destination is NULL or insufficiently large, the returned value is an upper bound on the number of triangles + * vertex_positions should have float3 position in the first 12 bytes of each vertex + * resolution is the dimension of the internal voxel grid and should be in the range [4, 256] + * options must be a bitmask composed of meshopt_RemeshX options; 0 is a safe default + */ +MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_remesh(float* destination, size_t max_triangle_count, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, unsigned int options); + /** * Quantize a float into half-precision (as defined by IEEE-754 fp16) floating point value * Generates +-inf for overflow, preserves NaN, flushes denormals to zero, rounds to nearest @@ -810,7 +1004,7 @@ MESHOPTIMIZER_API unsigned short meshopt_quantizeHalf(float v); /** * Quantize a float into a floating point value with a limited number of significant mantissa bits, preserving the IEEE-754 fp32 binary representation - * Generates +-inf for overflow, preserves NaN, flushes denormals to zero, rounds to nearest + * Preserves infinities/NaN, flushes denormals to zero, rounds to nearest * Assumes N is in a valid mantissa precision range, which is 1..23 */ MESHOPTIMIZER_API float meshopt_quantizeFloat(float v, int N); @@ -821,11 +1015,31 @@ MESHOPTIMIZER_API float meshopt_quantizeFloat(float v, int N); */ MESHOPTIMIZER_API float meshopt_dequantizeHalf(unsigned short h); +/** + * Experimental: Compute shared exponent suitable for mesh/cluster position quantization + * Given mesh or cluster bounds, compute a shared exponent that can be used to quantize any position inside the bounds to a 24-bit integer grid. + * The resulting output can be stored as a compact bit-stream to be decoded directly in shaders, or to be used as an input to RT BVH builders, + * for example via D3D12_VERTEX_FORMAT_COMPRESSED1 in DXR2 (max_bits=16). + * + * To quantize positions, compute: + * scale = pow(2, exponent) + * iv = int(round(v / scale)) + * The resulting integer can be stored as signed 24-bit, or as an unsigned offset from a signed 24-bit anchor value, shared between all positions. + * + * minv/maxv specify the axis-aligned bounding box of the mesh or cluster; each should refer to a float3 value + * min_exp specifies the minimum value for the returned exponent, limiting precision to reduce size; e.g. min_exp = -10 will produce minimum error of 1mm given metric units + * max_bits specifies the maximum allowed number of bits for the quantized integer range (offset from anchor is an unsigned integer up to 2^max_bits-1) + */ +MESHOPTIMIZER_EXPERIMENTAL int meshopt_computePositionExponent(const float* minv, const float* maxv, int min_exp, int max_bits); + /** * Set allocation callbacks * These callbacks will be used instead of the default operator new/operator delete for all temporary allocations in the library. * Note that all algorithms only allocate memory for temporary use. * allocate/deallocate are always called in a stack-like order - last pointer to be allocated is deallocated first. + * This function is not thread safe and must not be called concurrently with any other meshopt_ function. + * + * In shared library builds, allocations from templated index wrappers in this header will only be redirected if MESHOPTIMIZER_ALLOC_EXPORT is defined. */ MESHOPTIMIZER_API void meshopt_setAllocator(void* (MESHOPTIMIZER_ALLOC_CALLCONV* allocate)(size_t), void (MESHOPTIMIZER_ALLOC_CALLCONV* deallocate)(void*)); @@ -870,6 +1084,10 @@ inline size_t meshopt_generateVertexRemapCustom(unsigned int* destination, const template inline void meshopt_remapIndexBuffer(T* destination, const T* indices, size_t index_count, const unsigned int* remap); template +inline size_t meshopt_filterIndexBuffer(T* destination, const T* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride); +template +inline size_t meshopt_filterIndexBufferMulti(T* destination, const T* indices, size_t index_count, size_t vertex_count, const struct meshopt_Stream* streams, size_t stream_count); +template inline void meshopt_generateShadowIndexBuffer(T* destination, const T* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride); template inline void meshopt_generateShadowIndexBufferMulti(T* destination, const T* indices, size_t index_count, size_t vertex_count, const meshopt_Stream* streams, size_t stream_count); @@ -899,6 +1117,8 @@ template inline size_t meshopt_encodeIndexSequence(unsigned char* buffer, size_t buffer_size, const T* indices, size_t index_count); template inline int meshopt_decodeIndexSequence(T* destination, size_t index_count, const unsigned char* buffer, size_t buffer_size); +template +inline int meshopt_decodeMeshlet(V* vertices, size_t vertex_count, T* triangles, size_t triangle_count, const unsigned char* buffer, size_t buffer_size); inline size_t meshopt_encodeVertexBufferLevel(unsigned char* buffer, size_t buffer_size, const void* vertices, size_t vertex_count, size_t vertex_size, int level); template inline size_t meshopt_simplify(T* destination, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error, unsigned int options = 0, float* result_error = NULL); @@ -917,7 +1137,7 @@ inline size_t meshopt_stripify(T* destination, const T* indices, size_t index_co template inline size_t meshopt_unstripify(T* destination, const T* indices, size_t index_count, T restart_index); template -inline meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const T* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int warp_size, unsigned int buffer_size); +inline meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const T* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int warp_size, unsigned int primgroup_size); template inline meshopt_VertexFetchStatistics meshopt_analyzeVertexFetch(const T* indices, size_t index_count, size_t vertex_count, size_t vertex_size); template @@ -938,6 +1158,12 @@ template inline size_t meshopt_partitionClusters(unsigned int* destination, const T* cluster_indices, size_t total_index_count, const unsigned int* cluster_index_counts, size_t cluster_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t target_partition_size); template inline void meshopt_spatialSortTriangles(T* destination, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride); +template +inline void meshopt_generateTangents(float* result, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, const float* vertex_uvs, size_t vertex_uvs_stride, unsigned int options = 0); +template +inline void meshopt_generateNormals(float* result, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float crease_angle, float smoothing); +template +inline size_t meshopt_remesh(float* destination, size_t max_triangle_count, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, unsigned int options = 0); #endif /* Inline implementation */ @@ -1115,11 +1341,29 @@ template inline void meshopt_remapIndexBuffer(T* destination, const T* indices, size_t index_count, const unsigned int* remap) { meshopt_IndexAdapter in(NULL, indices, indices ? index_count : 0); - meshopt_IndexAdapter out(destination, 0, index_count); + meshopt_IndexAdapter out(destination, NULL, index_count); meshopt_remapIndexBuffer(out.data, indices ? in.data : NULL, index_count, remap); } +template +inline size_t meshopt_filterIndexBuffer(T* destination, const T* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride) +{ + meshopt_IndexAdapter in(NULL, indices, index_count); + meshopt_IndexAdapter out(destination, NULL, index_count); + + return meshopt_filterIndexBuffer(out.data, in.data, index_count, vertices, vertex_count, vertex_size, vertex_stride); +} + +template +inline size_t meshopt_filterIndexBufferMulti(T* destination, const T* indices, size_t index_count, size_t vertex_count, const struct meshopt_Stream* streams, size_t stream_count) +{ + meshopt_IndexAdapter in(NULL, indices, index_count); + meshopt_IndexAdapter out(destination, NULL, index_count); + + return meshopt_filterIndexBufferMulti(out.data, in.data, index_count, vertex_count, streams, stream_count); +} + template inline void meshopt_generateShadowIndexBuffer(T* destination, const T* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size, size_t vertex_stride) { @@ -1163,7 +1407,7 @@ inline size_t meshopt_generateProvokingIndexBuffer(T* destination, unsigned int* meshopt_IndexAdapter out(destination, NULL, index_count); size_t bound = vertex_count + (index_count / 3); - assert(size_t(T(bound - 1)) == bound - 1); // bound - 1 must fit in T + assert(bound == 0 || size_t(T(bound - 1)) == bound - 1); // bound - 1 must fit in T (void)bound; return meshopt_generateProvokingIndexBuffer(out.data, reorder, in.data, index_count, vertex_count); @@ -1255,6 +1499,15 @@ inline int meshopt_decodeIndexSequence(T* destination, size_t index_count, const return meshopt_decodeIndexSequence(destination, index_count, sizeof(T), buffer, buffer_size); } +template +inline int meshopt_decodeMeshlet(V* vertices, size_t vertex_count, T* triangles, size_t triangle_count, const unsigned char* buffer, size_t buffer_size) +{ + char types_valid[(sizeof(V) == 2 || sizeof(V) == 4) && (sizeof(T) == 1 || sizeof(T) == 4) ? 1 : -1]; + (void)types_valid; + + return meshopt_decodeMeshlet(vertices, vertex_count, sizeof(V), triangles, triangle_count, sizeof(T) == 1 ? 3 : 4, buffer, buffer_size); +} + inline size_t meshopt_encodeVertexBufferLevel(unsigned char* buffer, size_t buffer_size, const void* vertices, size_t vertex_count, size_t vertex_size, int level) { return meshopt_encodeVertexBufferLevel(buffer, buffer_size, vertices, vertex_count, vertex_size, level, -1); @@ -1326,17 +1579,17 @@ template inline size_t meshopt_unstripify(T* destination, const T* indices, size_t index_count, T restart_index) { meshopt_IndexAdapter in(NULL, indices, index_count); - meshopt_IndexAdapter out(destination, NULL, (index_count - 2) * 3); + meshopt_IndexAdapter out(destination, NULL, index_count == 0 ? 0 : (index_count - 2) * 3); return meshopt_unstripify(out.data, in.data, index_count, unsigned(restart_index)); } template -inline meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const T* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int warp_size, unsigned int buffer_size) +inline meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const T* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int warp_size, unsigned int primgroup_size) { meshopt_IndexAdapter in(NULL, indices, index_count); - return meshopt_analyzeVertexCache(in.data, index_count, vertex_count, cache_size, warp_size, buffer_size); + return meshopt_analyzeVertexCache(in.data, index_count, vertex_count, cache_size, warp_size, primgroup_size); } template @@ -1419,10 +1672,34 @@ inline void meshopt_spatialSortTriangles(T* destination, const T* indices, size_ meshopt_spatialSortTriangles(out.data, in.data, index_count, vertex_positions, vertex_count, vertex_positions_stride); } + +template +inline void meshopt_generateTangents(float* result, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, const float* vertex_uvs, size_t vertex_uvs_stride, unsigned int options) +{ + meshopt_IndexAdapter in(NULL, indices, indices ? index_count : 0); + + meshopt_generateTangents(result, indices ? in.data : NULL, index_count, vertex_positions, vertex_count, vertex_positions_stride, vertex_normals, vertex_normals_stride, vertex_uvs, vertex_uvs_stride, options); +} + +template +inline void meshopt_generateNormals(float* result, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float crease_angle, float smoothing) +{ + meshopt_IndexAdapter in(NULL, indices, indices ? index_count : 0); + + meshopt_generateNormals(result, indices ? in.data : NULL, index_count, vertex_positions, vertex_count, vertex_positions_stride, crease_angle, smoothing); +} + +template +inline size_t meshopt_remesh(float* destination, size_t max_triangle_count, const T* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, unsigned int options) +{ + meshopt_IndexAdapter in(NULL, indices, index_count); + + return meshopt_remesh(destination, max_triangle_count, in.data, index_count, vertex_positions, vertex_count, vertex_positions_stride, resolution, options); +} #endif /** - * Copyright (c) 2016-2025 Arseny Kapoulkine + * Copyright (c) 2016-2026 Arseny Kapoulkine * * Permission is hereby granted, free of charge, to any person * obtaining a copy of this software and associated documentation diff --git a/3rdparty/meshoptimizer/src/opacitymap.cpp b/3rdparty/meshoptimizer/src/opacitymap.cpp new file mode 100644 index 000000000..0ec5eaab4 --- /dev/null +++ b/3rdparty/meshoptimizer/src/opacitymap.cpp @@ -0,0 +1,557 @@ +// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details +#include "meshoptimizer.h" + +#include +#include +#include + +namespace meshopt +{ + +// opacity micromaps use a "bird" traversal order which recursively subdivides the triangles: +// https://docs.vulkan.org/spec/latest/_images/micromap-subd.svg +// note that triangles 0 and 2 have the same winding as the source triangle, however triangles 1 (flipped) +// and 3 (upright) have flipped winding; this is obvious from the level 2 subdivision in the diagram above +inline size_t getLevelSize(int level, int states) +{ + // 1-bit 2-state or 2-bit 4-state per micro triangle, rounded up to whole bytes + return ((1 << (level * 2)) * (states >> 1) + 7) >> 3; +} + +struct Texture +{ + const unsigned char* data; + size_t stride, pitch; + unsigned int width, height; + float widthf, heightf; // width * 256.f, height * 256.f +}; + +static float sampleTexture(const Texture& texture, float u, float v) +{ + // wrap texture coordinates; floor is expensive so only call it if we're outside of [0, 1] range (+eps) + u = fabsf(u - 0.5f) > 0.5f ? u - floorf(u) : u; + v = fabsf(v - 0.5f) > 0.5f ? v - floorf(v) : v; + + // convert from [0, 1] to 16.8 fixed point coordinates (rounded to nearest subpixel) with texel centers on integer grid + int uf = int(u * texture.widthf - 127.5f); + int vf = int(v * texture.heightf - 127.5f); + + // clamp to avoid extrapolation past left/top edge since we don't wrap across the edge + uf = uf < 0 ? 0 : uf; + vf = vf < 0 ? 0 : vf; + + // x/y are texel coordinates, rx/ry are subpixel offsets + int x = uf >> 8; + int y = vf >> 8; + int rx = uf & 255; + int ry = vf & 255; + + // safeguard: this should not happen but if it ever does, ensure the accesses are inbounds + if (unsigned(x) >= texture.width || unsigned(y) >= texture.height) + return 0.f; + + // clamp the offsets instead of wrapping for simplicity and performance + size_t offset = size_t(y) * texture.pitch + x * texture.stride; + size_t offsetx = (x + 1 < int(texture.width)) ? texture.stride : 0; + size_t offsety = (y + 1 < int(texture.height)) ? texture.pitch : 0; + + unsigned char a00 = texture.data[offset]; + unsigned char a10 = texture.data[offset + offsetx]; + unsigned char a01 = texture.data[offset + offsety]; + unsigned char a11 = texture.data[offset + offsetx + offsety]; + + // bilinear interpolation in integer space: result is 8.16 fixed point + int ax0 = a00 * 256 + (a10 - a00) * rx; + int ax1 = a01 * 256 + (a11 - a01) * rx; + int axy = ax0 * 256 + (ax1 - ax0) * ry; + + return float(axy) * (1.f / (255.f * 65536.f)); +} + +static unsigned int hashUpdate4u(unsigned int h, const unsigned char* key, size_t len) +{ + // MurmurHash2 + const unsigned int m = 0x5bd1e995; + const int r = 24; + + while (len >= 4) + { + unsigned int k; + memcpy(&k, key, sizeof(k)); + + k *= m; + k ^= k >> r; + k *= m; + + h *= m; + h ^= k; + + key += 4; + len -= 4; + } + + return h; +} + +struct TriangleOMM +{ + int uvs[6]; + int level; +}; + +struct TriangleOMMHasher +{ + const TriangleOMM* data; + + size_t hash(unsigned int index) const + { + const TriangleOMM& tri = data[index]; + + return hashUpdate4u(tri.level, reinterpret_cast(tri.uvs), sizeof(tri.uvs)); + } + + bool equal(unsigned int lhs, unsigned int rhs) const + { + const TriangleOMM& lt = data[lhs]; + const TriangleOMM& rt = data[rhs]; + + return lt.level == rt.level && memcmp(lt.uvs, rt.uvs, sizeof(lt.uvs)) == 0; + } +}; + +struct OMMHasher +{ + const unsigned char* data; + const unsigned int* offsets; + const unsigned char* levels; + int states; + + size_t hash(unsigned int index) const + { + const unsigned char* key = data + offsets[index]; + size_t size = getLevelSize(levels[index], states); + + unsigned int h = levels[index]; + + // MurmurHash2 for large keys, simple fold for small; note that size is a power of two + if (size < 4) + h ^= key[0] | (key[size - 1] << 8); + else + h = hashUpdate4u(h, key, size); + + // MurmurHash2 finalizer + h ^= h >> 13; + h *= 0x5bd1e995; + h ^= h >> 15; + return h; + } + + bool equal(unsigned int lhs, unsigned int rhs) const + { + size_t size = getLevelSize(levels[lhs], states); + + return levels[lhs] == levels[rhs] && memcmp(data + offsets[lhs], data + offsets[rhs], size) == 0; + } +}; + +static size_t hashBuckets3(size_t count) +{ + size_t buckets = 1; + while (buckets < count + count / 4) + buckets *= 2; + + return buckets; +} + +template +static T* hashLookup3(T* table, size_t buckets, const Hash& hash, const T& key, const T& empty) +{ + assert(buckets > 0); + assert((buckets & (buckets - 1)) == 0); + + size_t hashmod = buckets - 1; + size_t bucket = hash.hash(key) & hashmod; + + for (size_t probe = 0; probe <= hashmod; ++probe) + { + T& item = table[bucket]; + + if (item == empty) + return &item; + + if (hash.equal(item, key)) + return &item; + + // hash collision, quadratic probing + bucket = (bucket + probe + 1) & hashmod; + } + + assert(false && "Hash table is full"); // unreachable + return NULL; +} + +inline int quantizeSubpixel(float v, unsigned int size) +{ + return int(v * float(int(size) * 4) + (v >= 0 ? 0.5f : -0.5f)); +} + +static int rasterizeEdge(float u0, float v0, float u1, float v1, int edgeres, const Texture& texture) +{ + float edgestep = 1.f / float(edgeres + 1); + + float ud = (u1 - u0) * edgestep, vd = (v1 - v0) * edgestep; + float u = u0, v = v0; + + int mask = 0; + int count = 0; + + for (int i = 0; i < edgeres; ++i) + { + u += ud; + v += vd; + + float a = sampleTexture(texture, u, v); + mask |= (a >= 0.5f) << i; + count += a >= 0.5f; + } + + return mask | (count << 16); +} + +template +static void rasterizeOpacity0(unsigned char* result, size_t index, float a0, float a1, float a2, float ac, int e0, int e1, int e2, int edgeres) +{ + int states = States; + + // basic coverage estimator from center and corner values; trained to minimize error + float coverage = (a0 + a1 + a2) * 0.12f + ac * 0.64f; + + if (edgeres) + { + float edgescale = 1.f / edgeres; + + // if we have edge samples, we can get a better coverage estimate by including them; trained to minimize error + coverage = ac * 0.22f + float((e0 >> 16) + (e1 >> 16) + (e2 >> 16)) * edgescale * 0.23f + (a0 + a1 + a2) * 0.03f; + } + + if (states == 2) + { + result[index / 8] |= (coverage >= 0.5f) << (index % 8); + return; + } + + int transp = (a0 < 0.5f) & (a1 < 0.5f) & (a2 < 0.5f) & (ac < 0.5f); + int opaque = (a0 > 0.5f) & (a1 > 0.5f) & (a2 > 0.5f) & (ac > 0.5f); + + // treat state as known if thresholding of corners & centers against wider bounds is consistent + // for unknown states, we currently use the same formula as the 2-state opacity for better consistency with forced 2-state + int unknown = 2 + (coverage >= 0.5f); + int state = (transp | opaque) ? opaque : unknown; + + if (edgeres && (transp | opaque)) + { + // if we have edge samples, ensure they are consistent too, falling back to unknown if not + int exp = opaque ? (1 << edgeres) - 1 : 0; + int eok = ((e0 & 0xffff) == exp) & ((e1 & 0xffff) == exp) & ((e2 & 0xffff) == exp); + + state = eok ? state : unknown; + } + + result[index / 4] |= state << ((index % 4) * 2); +} + +template +static void rasterizeOpacity1(unsigned char* result, size_t index, int edgeres, const float* c0, const float* c1, const float* c2, const Texture& texture) +{ + // compute each edge midpoint & sample + float c01[3] = {(c0[0] + c1[0]) / 2, (c0[1] + c1[1]) / 2, 0.f}; + float c12[3] = {(c1[0] + c2[0]) / 2, (c1[1] + c2[1]) / 2, 0.f}; + float c20[3] = {(c2[0] + c0[0]) / 2, (c2[1] + c0[1]) / 2, 0.f}; + + c01[2] = sampleTexture(texture, c01[0], c01[1]); + c12[2] = sampleTexture(texture, c12[0], c12[1]); + c20[2] = sampleTexture(texture, c20[0], c20[1]); + + // corner tables for each edge, and corner + edge tables for each triangle + // edges are numbered counter clockwise, 6 outer first, 3 inner last; triangle vertex and edge references are in triangle winding order + static const unsigned char edges[9][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 0}, {5, 1}, {1, 3}, {3, 5}}; + static const unsigned char triangles[4][6] = {{0, 1, 5, 0, 6, 5}, {5, 3, 1, 8, 7, 6}, {1, 2, 3, 1, 2, 7}, {3, 5, 4, 8, 4, 3}}; + + const float* points[] = {c0, c01, c1, c12, c2, c20}; + + int em[9] = {}; + + // sample additional points on the edges to improve state estimation + if (edgeres > 0) + for (size_t i = 0; i < 9; ++i) + em[i] = rasterizeEdge(points[edges[i][0]][0], points[edges[i][0]][1], points[edges[i][1]][0], points[edges[i][1]][1], edgeres, texture); + + for (size_t i = 0; i < 4; ++i) + { + const unsigned char* tri = triangles[i]; + const float* p0 = points[tri[0]]; + const float* p1 = points[tri[1]]; + const float* p2 = points[tri[2]]; + + // compute triangle center & sample + float uc = (p0[0] + p1[0] + p2[0]) * (1.f / 3.f); + float vc = (p0[1] + p1[1] + p2[1]) * (1.f / 3.f); + float ac = sampleTexture(texture, uc, vc); + + // rasterize opacity state based on alpha values in corners and center (and optionally edges) + rasterizeOpacity0(result, index * 4 + i, p0[2], p1[2], p2[2], ac, em[tri[3]], em[tri[4]], em[tri[5]], edgeres); + } +} + +template +static void rasterizeOpacityRec(unsigned char* result, size_t index, int level, int edgeres, const float* c0, const float* c1, const float* c2, const Texture& texture) +{ + if (level == 0) + { + // compute triangle center & sample + float uc = (c0[0] + c1[0] + c2[0]) * (1.f / 3.f); + float vc = (c0[1] + c1[1] + c2[1]) * (1.f / 3.f); + float ac = sampleTexture(texture, uc, vc); + + int e0 = 0, e1 = 0, e2 = 0; + + if (edgeres > 0) + { + // sample additional points on the edges to improve state estimation + e0 = rasterizeEdge(c0[0], c0[1], c1[0], c1[1], edgeres, texture); + e1 = rasterizeEdge(c1[0], c1[1], c2[0], c2[1], edgeres, texture); + e2 = rasterizeEdge(c2[0], c2[1], c0[0], c0[1], edgeres, texture); + } + + // rasterize opacity state based on alpha values in corners and center (and optionally edges) + return rasterizeOpacity0(result, index, c0[2], c1[2], c2[2], ac, e0, e1, e2, edgeres); + } + + // fast path: equivalent to recursive rasterization, but reuses edge data to reduce sample count + if (level == 1 && edgeres > 0) + return rasterizeOpacity1(result, index, edgeres, c0, c1, c2, texture); + + // compute each edge midpoint & sample + float c01[3] = {(c0[0] + c1[0]) / 2, (c0[1] + c1[1]) / 2, 0.f}; + float c12[3] = {(c1[0] + c2[0]) / 2, (c1[1] + c2[1]) / 2, 0.f}; + float c20[3] = {(c2[0] + c0[0]) / 2, (c2[1] + c0[1]) / 2, 0.f}; + + c01[2] = sampleTexture(texture, c01[0], c01[1]); + c12[2] = sampleTexture(texture, c12[0], c12[1]); + c20[2] = sampleTexture(texture, c20[0], c20[1]); + + // recursively rasterize each triangle + // note: triangles 1 and 3 have flipped winding, and 1 is flipped upside down + rasterizeOpacityRec(result, index * 4 + 0, level - 1, edgeres, c0, c01, c20, texture); + rasterizeOpacityRec(result, index * 4 + 1, level - 1, edgeres, c20, c12, c01, texture); + rasterizeOpacityRec(result, index * 4 + 2, level - 1, edgeres, c01, c1, c12, texture); + rasterizeOpacityRec(result, index * 4 + 3, level - 1, edgeres, c12, c20, c2, texture); +} + +static int getSpecialIndex(const unsigned char* data, int level, int states) +{ + int first = data[0] & (states == 2 ? 1 : 3); + int special = -(1 + first); + + // at level 0, every micromap can be converted to a special index + if (level == 0) + return special; + + // at level 1 with 2 states, the byte is partially filled so we need a separate check + if (level == 1 && states == 2) + return (data[0] & 15) == ((-first) & 15) ? special : 0; + + // otherwise we need to check that all bytes are consistent with the first value and we can do this byte-wise + int expected = first * (states == 2 ? 0xff : 0x55); + size_t size = getLevelSize(level, states); + + for (size_t i = 0; i < size; ++i) + if (data[i] != expected) + return 0; + + return special; +} + +} // namespace meshopt + +size_t meshopt_opacityMapMeasure(unsigned char* levels, unsigned int* sources, int* omm_indices, const unsigned int* indices, size_t index_count, const float* vertex_uvs, size_t vertex_count, size_t vertex_uvs_stride, unsigned int texture_width, unsigned int texture_height, int max_level, float target_edge) +{ + using namespace meshopt; + + assert(index_count % 3 == 0); + assert(vertex_uvs_stride >= 8 && vertex_uvs_stride <= 256); + assert(vertex_uvs_stride % sizeof(float) == 0); + assert(unsigned(texture_width - 1) < 16384 && unsigned(texture_height - 1) < 16384); + assert(max_level >= 0 && max_level <= 12); + assert(target_edge >= 0); + + (void)vertex_count; + + meshopt_Allocator allocator; + + size_t vertex_stride_float = vertex_uvs_stride / sizeof(float); + float texture_area = float(texture_width) * float(texture_height); + + // hash map used to deduplicate triangle rasterization requests based on UV + size_t table_size = hashBuckets3(index_count / 3); + unsigned int* table = allocator.allocate(table_size); + memset(table, -1, table_size * sizeof(unsigned int)); + + TriangleOMM* triangles = allocator.allocate(index_count / 3); + TriangleOMMHasher hasher = {triangles}; + + size_t result = 0; + + for (size_t i = 0; i < index_count; i += 3) + { + unsigned int a = indices[i + 0], b = indices[i + 1], c = indices[i + 2]; + assert(a < vertex_count && b < vertex_count && c < vertex_count); + + float u0 = vertex_uvs[a * vertex_stride_float + 0], v0 = vertex_uvs[a * vertex_stride_float + 1]; + float u1 = vertex_uvs[b * vertex_stride_float + 0], v1 = vertex_uvs[b * vertex_stride_float + 1]; + float u2 = vertex_uvs[c * vertex_stride_float + 0], v2 = vertex_uvs[c * vertex_stride_float + 1]; + + int level = max_level; + + if (target_edge > 0) + { + // compute ratio of edge length (in texels) to target and determine subdivision level + float uvarea = fabsf((u1 - u0) * (v2 - v0) - (u2 - u0) * (v1 - v0)) * 0.5f * texture_area; + float ratio = sqrtf(uvarea) / target_edge; + float levelf = log2f(ratio > 1 ? ratio : 1); + + // round to nearest and clamp + level = int(levelf + 0.5f); + level = level < 0 ? 0 : level; + level = level < max_level ? level : max_level; + } + + // deduplicate rasterization requests based on UV + int su0 = quantizeSubpixel(u0, texture_width), sv0 = quantizeSubpixel(v0, texture_height); + int su1 = quantizeSubpixel(u1, texture_width), sv1 = quantizeSubpixel(v1, texture_height); + int su2 = quantizeSubpixel(u2, texture_width), sv2 = quantizeSubpixel(v2, texture_height); + + TriangleOMM tri = {{su0, sv0, su1, sv1, su2, sv2}, level}; + triangles[result] = tri; // speculatively write triangle data to give hasher a way to compare it + + unsigned int* entry = hashLookup3(table, table_size, hasher, unsigned(result), ~0u); + + if (*entry == ~0u) + { + *entry = unsigned(result); + levels[result] = (unsigned char)level; + sources[result] = unsigned(i / 3); + result++; + } + + omm_indices[i / 3] = int(*entry); + } + + return result; +} + +size_t meshopt_opacityMapEntrySize(int level, int states) +{ + assert(level >= 0 && level <= 12); + assert(states == 2 || states == 4); + + return meshopt::getLevelSize(level, states); +} + +void meshopt_opacityMapRasterize(unsigned char* result, int level, int states, const float* uv0, const float* uv1, const float* uv2, const unsigned char* texture_data, size_t texture_stride, size_t texture_pitch, unsigned int texture_width, unsigned int texture_height) +{ + using namespace meshopt; + + assert(level >= 0 && level <= 12); + assert(states == 2 || states == 4); + assert(unsigned(texture_width - 1) < 16384 && unsigned(texture_height - 1) < 16384); + assert(texture_stride >= 1 && texture_stride <= 4); + assert(texture_pitch >= texture_stride * texture_width); + + memset(result, 0, getLevelSize(level, states)); + + Texture texture = {texture_data, texture_stride, texture_pitch, texture_width, texture_height, float(int(texture_width)) * 256.f, float(int(texture_height)) * 256.f}; + + // determine number of edge samples for conservative state estimation + float texture_area = float(int(texture_width)) * float(int(texture_height)); + float uvarea = fabsf((uv1[0] - uv0[0]) * (uv2[1] - uv0[1]) - (uv2[0] - uv0[0]) * (uv1[1] - uv0[1])) * 0.5f * texture_area; + float uvedge = sqrtf(uvarea) / float(1 << level); + + // target ~2px distance between edge samples (assuming equilateral microtriangles) + int edgeres = int(uvedge * 0.75f); + edgeres = edgeres < 0 ? 0 : edgeres; + edgeres = edgeres > 7 ? 7 : edgeres; + + // rasterize all micro triangles recursively, passing corner data down to reduce redundant sampling + float c0[3] = {uv0[0], uv0[1], sampleTexture(texture, uv0[0], uv0[1])}; + float c1[3] = {uv1[0], uv1[1], sampleTexture(texture, uv1[0], uv1[1])}; + float c2[3] = {uv2[0], uv2[1], sampleTexture(texture, uv2[0], uv2[1])}; + + (states == 2 ? rasterizeOpacityRec<2> : rasterizeOpacityRec<4>)(result, 0, level, edgeres, c0, c1, c2, texture); +} + +size_t meshopt_opacityMapCompact(unsigned char* data, size_t data_size, unsigned char* levels, unsigned int* offsets, size_t omm_count, int* omm_indices, size_t triangle_count, int states) +{ + using namespace meshopt; + + assert(states == 2 || states == 4); + + meshopt_Allocator allocator; + + unsigned char* data_old = allocator.allocate(data_size); + memcpy(data_old, data, data_size); + + size_t table_size = hashBuckets3(omm_count); + unsigned int* table = allocator.allocate(table_size); + memset(table, -1, table_size * sizeof(unsigned int)); + + OMMHasher hasher = {data, offsets, levels, states}; + + int* remap = allocator.allocate(omm_count); + + size_t next = 0; + size_t offset = 0; + + for (size_t i = 0; i < omm_count; ++i) + { + int level = levels[i]; + assert(level >= 0 && level <= 12); + + const unsigned char* old = data_old + offsets[i]; + size_t size = getLevelSize(level, states); + assert(offsets[i] + size <= data_size); + + // try to convert to a special index if all micro-triangle states are the same + int special = getSpecialIndex(old, level, states); + if (special < 0) + { + remap[i] = special; + continue; + } + + // speculatively write data to give hasher a way to compare it + memcpy(data + offset, old, size); + offsets[next] = unsigned(offset); + levels[next] = (unsigned char)level; + + unsigned int* entry = hashLookup3(table, table_size, hasher, unsigned(next), ~0u); + + if (*entry == ~0u) + { + *entry = unsigned(next); + next++; + offset += size; + } + + remap[i] = int(*entry); + } + + // remap triangle indices to new indices or special indices + for (size_t i = 0; i < triangle_count; ++i) + { + assert(omm_indices[i] < 0 || unsigned(omm_indices[i]) < omm_count); + omm_indices[i] = omm_indices[i] < 0 ? omm_indices[i] : remap[omm_indices[i]]; + } + + return next; +} diff --git a/3rdparty/meshoptimizer/src/partition.cpp b/3rdparty/meshoptimizer/src/partition.cpp index 4119a53ed..c19cf19f7 100644 --- a/3rdparty/meshoptimizer/src/partition.cpp +++ b/3rdparty/meshoptimizer/src/partition.cpp @@ -268,24 +268,6 @@ static GroupOrder heapPop(GroupOrder* heap, size_t size) return top; } -static unsigned int countShared(const ClusterGroup* groups, int group1, int group2, const ClusterAdjacency& adjacency) -{ - unsigned int total = 0; - - for (int i1 = group1; i1 >= 0; i1 = groups[i1].next) - for (int i2 = group2; i2 >= 0; i2 = groups[i2].next) - { - for (unsigned int adj = adjacency.offsets[i1]; adj < adjacency.offsets[i1 + 1]; ++adj) - if (adjacency.clusters[adj] == unsigned(i2)) - { - total += adjacency.shared[adj]; - break; - } - } - - return total; -} - static void mergeBounds(ClusterGroup& target, const ClusterGroup& source) { float r1 = target.radius, r2 = source.radius; @@ -323,7 +305,7 @@ static float boundsScore(const ClusterGroup& target, const ClusterGroup& source) return mr > 0 ? r1 / mr : 0.f; } -static int pickGroupToMerge(const ClusterGroup* groups, int id, const ClusterAdjacency& adjacency, size_t max_partition_size, bool use_bounds) +static int pickGroupToMerge(const ClusterGroup* groups, int id, const ClusterAdjacency& adjacency, size_t max_partition_size, bool use_bounds, unsigned int* shared_acc, unsigned int& out_shared) { assert(groups[id].size > 0); @@ -331,20 +313,34 @@ static int pickGroupToMerge(const ClusterGroup* groups, int id, const ClusterAdj int best_group = -1; float best_score = 0; + unsigned int best_shared = 0; + // compute shared count based on underlying clusters for every adjacent group + // note: we assume shared_acc only contains zeros on input for (int ci = id; ci >= 0; ci = groups[ci].next) - { for (unsigned int adj = adjacency.offsets[ci]; adj != adjacency.offsets[ci + 1]; ++adj) { int other = groups[adjacency.clusters[adj]].group; - if (other < 0) + if (other >= 0) + shared_acc[other] += adjacency.shared[adj]; + } + + // score each adjacent group and pick the best one; resets shared_acc in the process + for (int ci = id; ci >= 0; ci = groups[ci].next) + for (unsigned int adj = adjacency.offsets[ci]; adj != adjacency.offsets[ci + 1]; ++adj) + { + int other = groups[adjacency.clusters[adj]].group; + if (other < 0 || shared_acc[other] == 0) continue; + // note: we reset shared_acc here to make sure the next invocation sees zero inputs; it's important that this is done before any early outs + unsigned int shared = shared_acc[other]; + shared_acc[other] = 0; + assert(groups[other].size > 0); if (groups[id].size + groups[other].size > max_partition_size) continue; - unsigned int shared = countShared(groups, id, other, adjacency); float other_rsqrt = 1.f / sqrtf(float(int(groups[other].vertices))); // normalize shared count by the expected boundary of each group (+ keeps scoring symmetric) @@ -358,10 +354,11 @@ static int pickGroupToMerge(const ClusterGroup* groups, int id, const ClusterAdj { best_group = other; best_score = score; + best_shared = shared; } } - } + out_shared = best_shared; return best_group; } @@ -513,6 +510,9 @@ size_t meshopt_partitionClusters(unsigned int* destination, const unsigned int* GroupOrder* order = allocator.allocate(cluster_count); size_t pending = 0; + unsigned int* shared_acc = allocator.allocate(cluster_count); + memset(shared_acc, 0, cluster_count * sizeof(unsigned int)); + // create a singleton group for each cluster and order them by priority for (size_t i = 0; i < cluster_count; ++i) { @@ -553,15 +553,13 @@ size_t meshopt_partitionClusters(unsigned int* destination, const unsigned int* if (groups[top.id].size >= target_partition_size) continue; - int best_group = pickGroupToMerge(groups, top.id, adjacency, max_partition_size, /* use_bounds= */ vertex_positions); + unsigned int best_shared = 0; + int best_group = pickGroupToMerge(groups, top.id, adjacency, max_partition_size, /* use_bounds= */ vertex_positions, shared_acc, best_shared); // we can't grow the group any more, emit as is if (best_group == -1) continue; - // compute shared vertices to adjust the total vertices estimate after merging - unsigned int shared = countShared(groups, top.id, best_group, adjacency); - // combine groups by linking them together unsigned int tail = top.id; while (groups[tail].next >= 0) @@ -572,7 +570,7 @@ size_t meshopt_partitionClusters(unsigned int* destination, const unsigned int* // update group sizes; note, the vertex update is a O(1) approximation which avoids recomputing the true size groups[top.id].size += groups[best_group].size; groups[top.id].vertices += groups[best_group].vertices; - groups[top.id].vertices = (groups[top.id].vertices > shared) ? groups[top.id].vertices - shared : 1; + groups[top.id].vertices = (groups[top.id].vertices > best_shared) ? groups[top.id].vertices - best_shared : 1; groups[best_group].size = 0; groups[best_group].vertices = 0; diff --git a/3rdparty/meshoptimizer/src/quantization.cpp b/3rdparty/meshoptimizer/src/quantization.cpp index 149835f50..857197302 100644 --- a/3rdparty/meshoptimizer/src/quantization.cpp +++ b/3rdparty/meshoptimizer/src/quantization.cpp @@ -2,6 +2,8 @@ #include "meshoptimizer.h" #include +#include +#include union FloatBits { @@ -74,3 +76,52 @@ float meshopt_dequantizeHalf(unsigned short h) u.ui = s | r; return u.f; } + +int meshopt_computePositionExponent(const float* minv, const float* maxv, int min_exp, int max_bits) +{ + assert(min_exp >= -126); + assert(max_bits >= 2 && max_bits <= 24); + + int exp = min_exp; + + // compute max absolute component to ensure that individual endpoints fit on a 24-bit signed grid + float maxc = 0.f; + + for (int k = 0; k < 3; ++k) + { + maxc = maxc < fabsf(minv[k]) ? fabsf(minv[k]) : maxc; + maxc = maxc < fabsf(maxv[k]) ? fabsf(maxv[k]) : maxc; + } + + int maxc_exp = 0; + float maxc_fr = frexpf(maxc, &maxc_exp); + + // maxc is representable as 2^(maxc_exp-24) * 24-bit *unsigned* integer + // we have to use a 24-bit *signed* grid, so we have to chop off the last bit + // however, rounding in the corner case (mantissa is 1.111...) may increase the exponent by 1 so we need to offset it + int maxc_off = 23 - (maxc_fr >= 1.f - FLT_EPSILON / 2); + + exp = exp < maxc_exp - maxc_off ? maxc_exp - maxc_off : exp; + + // compute effective range with conservative rounding, to allow for some ambiguity in the caller's rounding direction + float scale = ldexpf(1.f, -exp); + float range = 0.f; + + for (int k = 0; k < 3; ++k) + { + float a = floorf(minv[k] * scale); + float v = ceilf(maxv[k] * scale); + range = range < v - a ? v - a : range; + } + + // range must be representable as max_bits unsigned integer + int range_exp = 0; + float range_fr = frexpf(range, &range_exp); + + exp += (range_exp > max_bits) ? range_exp - max_bits : 0; + + // correct range if we are at the rounding boundary that would push us to overflow max_bits + exp += range_exp > max_bits && range_fr >= 1.f - 1.f / float((1 << max_bits) - 1); + + return exp; +} diff --git a/3rdparty/meshoptimizer/src/rasterizer.cpp b/3rdparty/meshoptimizer/src/rasterizer.cpp index bd788ffdb..0e640bc02 100644 --- a/3rdparty/meshoptimizer/src/rasterizer.cpp +++ b/3rdparty/meshoptimizer/src/rasterizer.cpp @@ -172,7 +172,7 @@ static float transformTriangles(float* triangles, const unsigned int* indices, s extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]); extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]); - float scale = kViewport / extent; + float scale = extent == 0.f ? 0.f : kViewport / extent; for (size_t i = 0; i < index_count; ++i) { diff --git a/3rdparty/meshoptimizer/src/remesher.cpp b/3rdparty/meshoptimizer/src/remesher.cpp new file mode 100644 index 000000000..45e72e49f --- /dev/null +++ b/3rdparty/meshoptimizer/src/remesher.cpp @@ -0,0 +1,751 @@ +// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details +#include "meshoptimizer.h" + +#include +#include +#include +#include + +#ifndef TRACE +#define TRACE 0 +#endif + +#if TRACE +#include +#endif + +// Note: this is only exposed for development purposes; do *not* use +enum +{ + meshopt_RemeshInternalDebug = 1 << 30 +}; + +// This work is based on: +// William Lorensen, Harvey Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. 1987 +// Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997 +namespace meshopt +{ + +// triangles are encoded as three vertex codes in three nibbles +// vertex codes match the bit order XYZ (e.g. 3 maps to X=1 Y=1 Z=0) +struct Case +{ + unsigned char code; + unsigned char altmode; + unsigned short triangles[5]; // up to 4 triangles and a null terminator + unsigned short alternate[5]; // alternate triangulation, restricted to have <= triangles vs primary +}; + +// the base cases establish a canonical triangulation; unlike classical Marching Cubes, we connect cube *corners* instead of edges +// each case may have an alternate triangulation, used based on deciders; crucially, open boundary edges must stay consistent as decider varies per cell +static const Case kBaseCases[] = { + {0x00, 0, {}, {}}, + {0x01, 0, {}, {}}, + {0x03, 0, {}, {}}, + {0x06, 0, {}, {}}, + {0x07, 0, {0x120}, {}}, + {0x0f, 0, {0x130, 0x320}, {}}, + {0x16, 1, {0x124, 0x142}, {}}, + {0x17, 0, {0x124}, {}}, + {0x18, 0, {}, {}}, + {0x19, 1, {0x043, 0x403}, {}}, + {0x1b, 1, {0x134, 0x430}, {0x410, 0x013}}, + {0x1d, 1, {0x243, 0x340}, {0x320, 0x024}}, + {0x1e, 1, {0x142, 0x432, 0x341}, {0x132}}, + {0x1f, 2, {0x412, 0x213}, {0x134, 0x432}}, + {0x3c, 1, {0x345, 0x324, 0x543, 0x423}, {}}, + {0x3d, 1, {0x253, 0x524, 0x350}, {0x240, 0x504, 0x320}}, + {0x3f, 0, {0x253, 0x452}, {}}, + {0x69, 1, {0x035, 0x560, 0x306, 0x536}, {}}, + {0x6b, 1, {0x063, 0x605, 0x365}, {0x305}}, + {0x6f, 1, {0x365, 0x560}, {0x530, 0x036}}, + {0x7e, 0, {0x142, 0x536}, {}}, + {0x7f, 0, {0x365}, {}}, + {0xf0, 0, {0x547, 0x746}, {}}, + {0xff, 0, {}, {}}, +}; + +// for each cube case, a triangulation is derived from the base cases via 90 degree rotations +static unsigned short kTriangleTable[256][2][5]; +static unsigned char kTriangleCount[256]; +static unsigned char kTriangleAlt[256]; + +static bool buildRemeshTables() +{ + // corner remap for 90 degree rotations around X and Y + static const unsigned char rotations[2][8] = {{2, 3, 6, 7, 0, 1, 4, 5}, {4, 0, 6, 2, 5, 1, 7, 3}}; + + unsigned char filled[256] = {}; + + // copy base cases as is; we will generate the rest via rotations + for (size_t i = 0; i < sizeof(kBaseCases) / sizeof(kBaseCases[0]); ++i) + { + const Case& c = kBaseCases[i]; + memcpy(kTriangleTable[c.code][0], c.triangles, sizeof(c.triangles)); + memcpy(kTriangleTable[c.code][1], c.altmode ? c.alternate : c.triangles, sizeof(c.triangles)); + kTriangleAlt[c.code] = c.altmode; + filled[kBaseCases[i].code] = 1; + } + + // propagate base cases iteratively to fill all cube cases via rotations + // note: the worst case is 5 passes but in practice we only need 3 to fill all cases + for (int pass = 0; pass < 3; ++pass) + for (int code = 0; code < 256; ++code) + { + if (filled[code] != 1) + continue; + + for (int r = 0; r < 2; ++r) + { + unsigned char rotated = 0; + for (int i = 0; i < 8; ++i) + rotated |= ((code >> i) & 1) << rotations[r][i]; + + if (filled[rotated]) + continue; + + for (int k = 0; k < 2; ++k) + for (int i = 0; kTriangleTable[code][k][i]; ++i) + { + unsigned short tri = kTriangleTable[code][k][i]; + kTriangleTable[rotated][k][i] = (rotations[r][(tri >> 8) & 0xf] << 8) | (rotations[r][(tri >> 4) & 0xf] << 4) | rotations[r][tri & 0xf]; + } + + kTriangleAlt[rotated] = kTriangleAlt[code]; + filled[rotated] = 1; // mark as pending + } + + filled[code] = 2; // mark as processed + } + + // finalize auxiliary tables + for (int code = 0; code < 256; ++code) + { + assert(filled[code]); + + int count = 0; + while (kTriangleTable[code][0][count]) + count++; + + kTriangleCount[code] = (unsigned char)count; + + // counting pass does not use alternate triangulations, so they need to stay at or below primary triangle count + assert(kTriangleTable[code][1][count] == 0); + + // quadric decider expects a canonical quad encoding (0xabc 0xcbd) + if (kTriangleAlt[code] == 2) + assert(count == 2 && (kTriangleTable[code][0][0] & 0xf0) == (kTriangleTable[code][0][1] & 0xf0) && (kTriangleTable[code][0][0] & 0xf) == (kTriangleTable[code][0][1] >> 8)); + } + + return true; +} + +static bool gRemeshTablesInitialized = buildRemeshTables(); + +struct Voxel +{ + unsigned int coord; + unsigned char octants; + + float px, py, pz; + float w; + + // a00*x^2 + a11*y^2 + a22*z^2 + 2*a10*xy + 2*a20*xz + 2*a21*yz + 2*b0*x + 2*b1*y + 2*b2*z + c + float a00, a11, a22; + float a10, a20, a21; + float b0, b1, b2, c; +}; + +static void voxelAccumulate(Voxel& vox, float px, float py, float pz, float w) +{ + vox.px += px * w; + vox.py += py * w; + vox.pz += pz * w; + vox.w += w; +} + +static void voxelAccumulateQ(Voxel& vox, float a, float b, float c, float d, float w) +{ + float aw = a * w; + float bw = b * w; + float cw = c * w; + float dw = d * w; + + vox.a00 += a * aw; + vox.a11 += b * bw; + vox.a22 += c * cw; + vox.a10 += a * bw; + vox.a20 += a * cw; + vox.a21 += b * cw; + vox.b0 += a * dw; + vox.b1 += b * dw; + vox.b2 += c * dw; + vox.c += d * dw; +} + +static float voxelError(const Voxel& vox, float x, float y, float z) +{ + float rx = (vox.b0 + vox.a10 * y) * 2.f + vox.a00 * x; + float ry = (vox.b1 + vox.a21 * z) * 2.f + vox.a11 * y; + float rz = (vox.b2 + vox.a20 * x) * 2.f + vox.a22 * z; + + return fabsf(vox.c + rx * x + ry * y + rz * z); +} + +static bool voxelSolve(float& rx, float& ry, float& rz, const Voxel& vox, float lambda) +{ + // solve A*p = -b where A is the quadric matrix and b is the linear term + float rw = lambda * vox.w; + float a00 = rw + vox.a00, a11 = rw + vox.a11, a22 = rw + vox.a22; + float a10 = vox.a10, a20 = vox.a20, a21 = vox.a21; + float x0 = lambda * vox.px - vox.b0, x1 = lambda * vox.py - vox.b1, x2 = lambda * vox.pz - vox.b2; + + float eps = 1e-6f * vox.w; + + // LDL decomposition: A = LDL^T + float d0 = a00; + float l10 = a10 / d0; + float l20 = a20 / d0; + + float d1 = a11 - a10 * l10; + float dl21 = a21 - a20 * l10; + float l21 = dl21 / d1; + + float d2 = a22 - a20 * l20 - dl21 * l21; + + // solve L*y = x + float y0 = x0; + float y1 = x1 - l10 * y0; + float y2 = x2 - l20 * y0 - l21 * y1; + + // solve D*z = y + float z0 = y0 / d0; + float z1 = y1 / d1; + float z2 = y2 / d2; + + // substitute L^T*p = z + rz = z2; + ry = z1 - l21 * rz; + rx = z0 - l10 * ry - l20 * rz; + + return fabsf(d0) > eps && fabsf(d1) > eps && fabsf(d2) > eps; +} + +static float measureGrid(const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, int resolution, float* out_offset) +{ + size_t vertex_stride_float = vertex_positions_stride / sizeof(float); + + float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX}; + float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX}; + + for (size_t i = 0; i < vertex_count; ++i) + { + const float* v = vertex_positions_data + i * vertex_stride_float; + + for (int j = 0; j < 3; ++j) + { + float vj = v[j]; + + minv[j] = minv[j] > vj ? vj : minv[j]; + maxv[j] = maxv[j] < vj ? vj : maxv[j]; + } + } + + float extent = 0.f; + + extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]); + extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]); + extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]); + + // small subvoxel padding keeps vertices inside the last valid voxel to avoid clamping correctness issues + extent *= (resolution + 1e-2f) / float(resolution); + + // rescale extents to [0..resolution - 2], because the first and last voxels are used as empty padding + float scale = extent == 0 ? 0.f : (resolution - 2) / extent; + + // center mesh in the grid; this improves voxelization for symmetric inputs + for (int j = 0; j < 3; ++j) + out_offset[j] = minv[j] - (extent - (maxv[j] - minv[j])) * 0.5f; + + return scale; +} + +static void voxelize(unsigned char* grid, Voxel* voxels, const unsigned int* voxel_rows, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, float scale, const float offset[3], unsigned int options) +{ + (void)vertex_count; + + size_t vertex_stride_float = vertex_positions_stride / sizeof(float); + + for (size_t i = 0; i < index_count; i += 3) + { + unsigned int a = indices[i + 0], b = indices[i + 1], c = indices[i + 2]; + assert(a < vertex_count && b < vertex_count && c < vertex_count); + + const float* va = vertex_positions + a * vertex_stride_float; + const float* vb = vertex_positions + b * vertex_stride_float; + const float* vc = vertex_positions + c * vertex_stride_float; + + float ex = vb[0] - va[0], ey = vb[1] - va[1], ez = vb[2] - va[2]; + float fx = vc[0] - va[0], fy = vc[1] - va[1], fz = vc[2] - va[2]; + float gx = vc[0] - vb[0], gy = vc[1] - vb[1], gz = vc[2] - vb[2]; + + // use maximum edge length to establish sampling rate + // TODO: this is wasteful for thin triangles + float el = sqrtf(ex * ex + ey * ey + ez * ez); + float fl = sqrtf(fx * fx + fy * fy + fz * fz); + float gl = sqrtf(gx * gx + gy * gy + gz * gz); + + float max_edge = el > fl ? el : fl; + max_edge = max_edge > gl ? max_edge : gl; + + // we target 2 samples per voxel edge which should be enough to hit all voxels at any rotation + int samples = int(max_edge * scale * 2.f); + samples = samples > 1 ? samples : 1; + samples = samples < resolution * 2 ? samples : resolution * 2; + + // normal is used to compute area for weighting as well as for the quadric + float nx = ey * fz - ez * fy, ny = ez * fx - ex * fz, nz = ex * fy - ey * fx; + float area = sqrtf(nx * nx + ny * ny + nz * nz); + + float ns = area == 0.f ? 0.f : 1.f / area; + nx *= ns; + ny *= ns; + nz *= ns; + + float sx = va[0] - offset[0], sy = va[1] - offset[1], sz = va[2] - offset[2]; + float sr = 1.f / float(samples); + float weight = area / float((samples + 1) * (samples + 2)); + + // skip degenerate triangles; they don't contribute to voxelization and may result in NaN when computing voxel centroid + if (weight == 0.f) + continue; + + for (int u = 0; u <= samples; ++u) + for (int v = 0; v <= samples - u; ++v) + { + float su = float(u) * sr, sv = float(v) * sr; + float px = sx + su * ex + sv * fx; + float py = sy + su * ey + sv * fy; + float pz = sz + su * ez + sv * fz; + + // compute coordinates with an extra bit to use in octant mask + int hx = int(px * (scale * 2)); + int hy = int(py * (scale * 2)); + int hz = int(pz * (scale * 2)); + int x = hx >> 1, y = hy >> 1, z = hz >> 1; + + // safety: rounding errors and non-finite inputs may produce out of bounds coordinates, so we clamp them + int cutoff = resolution - 3; + x = unsigned(x) < unsigned(cutoff) ? x : cutoff; + y = unsigned(y) < unsigned(cutoff) ? y : cutoff; + z = unsigned(z) < unsigned(cutoff) ? z : cutoff; + + size_t row = (y + 1) + size_t(resolution) * (z + 1); + size_t idx = (x + 1) + size_t(resolution) * row; + + if (voxels) + { + assert(grid[idx] != 0 && grid[idx] != 0xff); + Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)]; + + vox.coord = (unsigned(x) << 20) | (unsigned(y) << 10) | unsigned(z); + vox.octants |= 1 << ((hx & 1) | ((hy & 1) << 1) | ((hz & 1) << 2)); + + voxelAccumulate(vox, px, py, pz, weight); + + if (options & meshopt_RemeshSolve) + { + float distance = nx * px + ny * py + nz * pz; + voxelAccumulateQ(vox, nx, ny, nz, -distance, weight); + } + } + else + { + grid[idx] = 1; + } + } + } +} + +static size_t rowpack(unsigned char* grid, unsigned int* voxel_rows, int resolution) +{ + size_t result = 0; + size_t slice = size_t(resolution) * size_t(resolution); + + for (size_t i = 0; i < slice; ++i) + { + unsigned char* data = grid + i * size_t(resolution); + + int count = 0; + + for (int x = 0; x < resolution; ++x) + { + assert(data[x] <= 1); // voxelize only produces 0/1 values + + count += data[x]; + data[x] = (data[x] != 0) ? (unsigned char)count : 0; + } + + assert(count < 255); // we store offsets in a single byte, with 0 reserved for empty voxels and 0xff reserved for interior voxels + + voxel_rows[i] = unsigned(result); + result += count; + } + + return result; +} + +static void solidifyQueue(unsigned int row, unsigned int* worklist, unsigned char* queued, size_t& pending) +{ + if (queued[row]) + return; + + queued[row] = 1; + worklist[pending++] = row; +} + +static void solidify(unsigned char* grid, unsigned int* worklist, unsigned char* queued, int resolution) +{ + size_t pending = 0; + memset(queued, 0, size_t(resolution) * size_t(resolution)); + + // mark the interior empty voxels as 'inside'; we will propagate 'empty' state from the boundary + for (int z = 1; z < resolution - 1; ++z) + for (int y = 1; y < resolution - 1; ++y) + { + unsigned char* data = grid + size_t(resolution) * (y + size_t(resolution) * z); + + for (int x = 1; x < resolution - 1; ++x) + data[x] = (data[x] == 0) ? 0xff : data[x]; + } + + // queue all rows for processing; note that boundary rows are always empty but we need them to propagate into neighboring rows + for (int z = 0; z < resolution; ++z) + for (int y = 0; y < resolution; ++y) + solidifyQueue(y + size_t(resolution) * z, worklist, queued, pending); + + while (pending) + { + unsigned int row = worklist[--pending]; + assert(queued[row]); + queued[row] = 0; + + unsigned char* data = grid + size_t(resolution) * row; + + // propagate outside state to the interior within row + for (int x = 1; x < resolution - 1; ++x) + data[x] = (data[x] == 0xff && data[x - 1] == 0) ? 0 : data[x]; + + for (int x = resolution - 2; x >= 1; --x) + data[x] = (data[x] == 0xff && data[x + 1] == 0) ? 0 : data[x]; + + // propagate outside state to the interior into neighboring rows + int y = row % resolution, z = row / resolution; + + for (int k = 0; k < 4; ++k) + { + int yn = y + (k == 0 ? -1 : (k == 1 ? 1 : 0)); + int zn = z + (k == 2 ? -1 : (k == 3 ? 1 : 0)); + + if (yn < 1 || yn >= resolution - 1 || zn < 1 || zn >= resolution - 1) + continue; + + unsigned char* datan = grid + size_t(resolution) * (yn + size_t(resolution) * zn); + unsigned char changed = 0; + + for (int x = 1; x < resolution - 1; ++x) + { + unsigned char rep = ((data[x] == 0) & (datan[x] == 0xff)) ? 0 : datan[x]; + changed |= rep ^ datan[x]; + datan[x] = rep; + } + + if (changed) + solidifyQueue(yn + size_t(resolution) * zn, worklist, queued, pending); + } + } +} + +static void solve(Voxel* voxels, size_t voxel_count, float scale, unsigned int options) +{ + // regularization factor is a tradeoff between solve stability/uniformity and adherence to the original surface + float factor = 3e-2f; + + // cutoff is a squared distance in model space; sqrt(3)/scale is the voxel diagonal + float cutoff = 3.f / (scale * scale); + float rscale = 1.f / scale; + + size_t solved = 0; + + for (size_t i = 0; i < voxel_count; ++i) + { + Voxel& vox = voxels[i]; + float px = vox.px / vox.w, py = vox.py / vox.w, pz = vox.pz / vox.w; + + // voxel corner, used for clamping and debug visualization + float cx = int((vox.coord >> 20) & 1023) * rscale; + float cy = int((vox.coord >> 10) & 1023) * rscale; + float cz = int((vox.coord >> 0) & 1023) * rscale; + + if (options & meshopt_RemeshSolve) + { + // compute minimizing point; this includes regularization to stabilize the solve + float sx, sy, sz; + if (voxelSolve(sx, sy, sz, vox, factor)) + { + float d = (sx - px) * (sx - px) + (sy - py) * (sy - py) + (sz - pz) * (sz - pz); + + // reject solutions that move the vertex too far from voxel centroid; this is a safety net in case regularization is insufficient + if (d < cutoff) + { + // use a clamped solution to reduce intersections + px = sx > cx ? sx : cx; + px = px < cx + rscale ? px : cx + rscale; + py = sy > cy ? sy : cy; + py = py < cy + rscale ? py : cy + rscale; + pz = sz > cz ? sz : cz; + pz = pz < cz + rscale ? pz : cz + rscale; + + solved++; + } + } + } + + if (options & meshopt_RemeshInternalDebug) + { + px = cx + 0.5f * rscale; + py = cy + 0.5f * rscale; + pz = cz + 0.5f * rscale; + } + + vox.px = px; + vox.py = py; + vox.pz = pz; + } + + (void)solved; + +#if TRACE + if (options & meshopt_RemeshSolve) + printf("remesher: solved %zu/%zu voxels\n", solved, voxel_count); +#endif +} + +static void emitVertex(float* result, int x, int y, int z, int corner, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution, const float offset[3]) +{ + int ox = corner & 1, oy = (corner >> 1) & 1, oz = (corner >> 2) & 1; + + size_t row = (y + oy) + size_t(resolution) * (z + oz); + size_t idx = (x + ox) + size_t(resolution) * row; + + assert(grid[idx] != 0 && grid[idx] != 0xff); + const Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)]; + + result[0] = vox.px + offset[0]; + result[1] = vox.py + offset[1]; + result[2] = vox.pz + offset[2]; +} + +static bool octantDecider(int x, int y, int z, int cube, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution) +{ + for (int c = 0; c < 8; ++c) + if (cube & (1 << c)) + { + int ox = c & 1, oy = (c >> 1) & 1, oz = (c >> 2) & 1; + + size_t row = (y + oy) + size_t(resolution) * (z + oz); + size_t idx = (x + ox) + size_t(resolution) * row; + + assert(grid[idx] != 0); + + // note: unlike vertex emission, which never consults interior occupied corners, here we need to skip corners without an occupied voxel + if (grid[idx] == 0xff) + continue; + + const Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)]; + + // test octant contained within the cell (opposite of the corner index) + if (vox.octants & (1 << (7 - c))) + return false; + } + + // select alternate configuration if *all* octants are empty + return true; +} + +static bool quadricDecider(int x, int y, int z, int cube, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution) +{ + // quads are encoded implicitly as 0xabc 0xcbd + unsigned int quad = (kTriangleTable[cube][0][0] << 4) | (kTriangleTable[cube][0][1] & 0xf); + + const Voxel* corner[4]; + + for (int i = 0; i < 4; ++i) + { + int c = (quad >> (12 - i * 4)) & 0xf; + size_t row = (y + ((c >> 1) & 1)) + size_t(resolution) * (z + ((c >> 2) & 1)); + size_t idx = (x + (c & 1)) + size_t(resolution) * row; + + assert(grid[idx] != 0 && grid[idx] != 0xff); + corner[i] = &voxels[voxel_rows[row] + (grid[idx] - 1)]; + } + + // evaluate error for midpoints of primary (bc) and alternate (ad) diagonal + float mx = (corner[1]->px + corner[2]->px) * 0.5f, my = (corner[1]->py + corner[2]->py) * 0.5f, mz = (corner[1]->pz + corner[2]->pz) * 0.5f; + float nx = (corner[0]->px + corner[3]->px) * 0.5f, ny = (corner[0]->py + corner[3]->py) * 0.5f, nz = (corner[0]->pz + corner[3]->pz) * 0.5f; + float error0 = 0, error1 = 0; + + for (int i = 0; i < 4; ++i) + { + error0 += voxelError(*corner[i], mx, my, mz); + error1 += voxelError(*corner[i], nx, ny, nz); + } + + // select alternate configuration if it is clearly better + return error1 < error0 * 0.9f; +} + +static size_t polygonize(float* destination, size_t max_triangle_count, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution, const float offset[3], unsigned int options) +{ + size_t result = 0; + size_t slice = size_t(resolution) * size_t(resolution); + + assert(gRemeshTablesInitialized); + + for (int z = 0; z < resolution - 1; ++z) + for (int y = 0; y < resolution - 1; ++y) + { + const unsigned char* data = grid + size_t(resolution) * (y + size_t(resolution) * z); + + // we track each slice as a 8-bit code (matching cube indexing) as we iterate through the row to avoid extra lookups + int last = (data[0] != 0) | ((data[resolution] != 0) << 2) | ((data[slice] != 0) << 4) | ((data[slice + resolution] != 0) << 6); + + for (int x = 0; x < resolution - 1; ++x) + { + int next = (data[x + 1] != 0) | ((data[x + 1 + resolution] != 0) << 2) | ((data[x + 1 + slice] != 0) << 4) | ((data[x + 1 + slice + resolution] != 0) << 6); + + // next slice bits are in position for the next iteration so we need to shift them to retrieve current cube code + int cube = last | (next << 1); + last = next; + + if (cube == 0 || cube == 0xff) + continue; + + if (!destination) + { + // fast path: we can statically determine the upper bound on triangles the loop below will output based on cube configuration + result += kTriangleCount[cube]; + continue; + } + + // deciders are only evaluated for cells with an alternate configuration to reduce overhead + int alt = 0; + + if (kTriangleAlt[cube] == 1) + alt = octantDecider(x, y, z, cube, grid, voxels, voxel_rows, resolution); + else if (kTriangleAlt[cube] == 2 && (options & meshopt_RemeshSolve)) + alt = quadricDecider(x, y, z, cube, grid, voxels, voxel_rows, resolution); + + const unsigned short* tris = kTriangleTable[cube][alt]; + + for (int i = 0; tris[i]; ++i) + { + // note: we only emit the triangle if we have space for it, but we count it regardless for consistent capacity estimation + if (result < max_triangle_count) + { + unsigned short tri = tris[i]; + + emitVertex(&destination[result * 9 + 0], x, y, z, (tri >> 8) & 0xf, grid, voxels, voxel_rows, resolution, offset); + emitVertex(&destination[result * 9 + 3], x, y, z, (tri >> 4) & 0xf, grid, voxels, voxel_rows, resolution, offset); + emitVertex(&destination[result * 9 + 6], x, y, z, (tri >> 0) & 0xf, grid, voxels, voxel_rows, resolution, offset); + } + + result++; + } + } + } + + return result; +} + +} // namespace meshopt + +size_t meshopt_remesh(float* destination, size_t max_triangle_count, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, unsigned int options) +{ + using namespace meshopt; + + assert(destination || max_triangle_count == 0); + assert(index_count % 3 == 0); + assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); + assert(vertex_positions_stride % sizeof(float) == 0); + assert(resolution >= 4 && resolution <= 256); + + meshopt_Allocator allocator; + + // measure voxel grid to compute position => voxel mapping + float offset[3] = {}; + float scale = measureGrid(vertex_positions, vertex_count, vertex_positions_stride, resolution, offset); + + // rasterize triangles into the voxel grid: in the first pass, this simply tags occupied voxels + unsigned char* grid = allocator.allocate(size_t(resolution) * size_t(resolution) * size_t(resolution)); + memset(grid, 0, size_t(resolution) * size_t(resolution) * size_t(resolution)); + + voxelize(grid, NULL, NULL, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride, resolution, scale, offset, options); + + // allocate additional voxel data for each occupied voxel; this can be filled in the second pass to compute positions + // note that we only do this if we need to compute output triangles; counting runs skip it for performance + Voxel* voxels = NULL; + unsigned int* voxel_rows = NULL; + size_t voxel_count = 0; + + if (destination) + { + voxel_rows = allocator.allocate(size_t(resolution) * size_t(resolution)); + voxel_count = rowpack(grid, voxel_rows, resolution); + + voxels = allocator.allocate(voxel_count); + memset(voxels, 0, voxel_count * sizeof(Voxel)); + +#if TRACE + printf("remesher: %zu voxels occupied\n", voxel_count); +#endif + } + + // fill in empty voxels that are not reachable from the grid boundary; the inside empty voxels are marked with 0xff + // note that these voxels do *not* have the associated Voxel data as they were never voxelized + if ((options & meshopt_RemeshShell) == 0) + { + unsigned int* worklist = allocator.allocate(size_t(resolution) * size_t(resolution)); + unsigned char* queued = allocator.allocate(size_t(resolution) * size_t(resolution)); + + solidify(grid, worklist, queued, resolution); + +#if TRACE + size_t inside_count = 0, occupied_count = 0; + for (size_t i = 0; i < size_t(resolution) * size_t(resolution) * size_t(resolution); ++i) + { + inside_count += (grid[i] == 0xff); + occupied_count += (grid[i] != 0 && grid[i] != 0xff); + } + + printf("remesher: %zu voxels occupied, %zu voxels inside\n", occupied_count, inside_count); +#endif + } + + // accumulate voxel positions: in the second pass, this computes enough data in each voxel to calculate positions + if (voxels) + voxelize(grid, voxels, voxel_rows, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride, resolution, scale, offset, options); + + // compute final voxel positions; each voxel has a single resulting position that will be emitted during polygonization + if (voxels) + solve(voxels, voxel_count, scale, options); + + // output triangles from the voxel grid; if destination is NULL, this still counts the number of triangles that would be generated + size_t result = polygonize(destination, max_triangle_count, grid, voxels, voxel_rows, resolution, offset, options); + +#if TRACE + printf("remesher: %zu triangles (%zu capacity)\n", result, max_triangle_count); +#endif + + return result; +} diff --git a/3rdparty/meshoptimizer/src/simplifier.cpp b/3rdparty/meshoptimizer/src/simplifier.cpp index efbdf91a6..e9774b237 100644 --- a/3rdparty/meshoptimizer/src/simplifier.cpp +++ b/3rdparty/meshoptimizer/src/simplifier.cpp @@ -367,8 +367,8 @@ static void classifyVertices(unsigned char* result, unsigned int* loop, unsigned memset(loopback, -1, vertex_count * sizeof(unsigned int)); // incoming & outgoing open edges: ~0u if no open edges, i if there are more than 1 - // note that this is the same data as required in loop[] arrays; loop[] data is only valid for border/seam - // but here it's okay to fill the data out for other types of vertices as well + // note that this is the same data as required in loop[] arrays; loop[] data is only used for border/seam by default + // in permissive mode we also use it to guide complex-complex collapses, so we fill it for all vertices unsigned int* openinc = loopback; unsigned int* openout = loop; @@ -941,7 +941,7 @@ static void quadricFromAttributes(Quadric& Q, QuadricGrad* G, const Vector3& p0, Q.c += w * (gw * gw); - // the only remaining sum components are ones that depend on attr; these will be addded during error evaluation, see quadricError + // the only remaining sum components are ones that depend on attr; these will be added during error evaluation, see quadricError G[k].gx = w * gx; G[k].gy = w * gy; G[k].gz = w * gz; @@ -1081,18 +1081,22 @@ static void fillFaceQuadrics(Quadric* vertex_quadrics, QuadricGrad* volume_gradi } } -static void fillVertexQuadrics(Quadric* vertex_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* remap, unsigned int options) +static void fillVertexQuadrics(Quadric* vertex_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* remap, const unsigned char* vertex_lock, const unsigned int* sparse_remap, unsigned int options) { // by default, we use a very small weight to improve triangulation and numerical stability without affecting the shape or error - float factor = (options & meshopt_SimplifyRegularize) ? 1e-1f : 1e-7f; + float factor = (options & meshopt_SimplifyRegularizeLight) ? 1e-2f : ((options & meshopt_SimplifyRegularize) ? 1e-1f : 1e-7f); for (size_t i = 0; i < vertex_count; ++i) { if (remap[i] != i) continue; + // increase regularization weight for vertices marked as priority; for now we only examine the primary vertex + unsigned int ri = sparse_remap ? sparse_remap[i] : unsigned(i); + bool priority = vertex_lock && (vertex_lock[ri] & meshopt_SimplifyVertex_Priority) != 0; + const Vector3& p = vertex_positions[i]; - float w = vertex_quadrics[i].w * factor; + float w = vertex_quadrics[i].w * (priority ? 1.0f : factor); Quadric Q; quadricFromPoint(Q, p.x, p.y, p.z, w); @@ -1152,6 +1156,68 @@ static void fillEdgeQuadrics(Quadric* vertex_quadrics, const unsigned int* indic } } +static unsigned int oppositeVertex(const EdgeAdjacency& adjacency, unsigned int a, unsigned int b, const unsigned int* remap, const unsigned int* wedge) +{ + unsigned int v = a; + unsigned int r = ~0u; + + do + { + unsigned int count = adjacency.offsets[v + 1] - adjacency.offsets[v]; + const EdgeAdjacency::Edge* edges = adjacency.data + adjacency.offsets[v]; + + for (size_t i = 0; i < count; ++i) + if (remap[edges[i].next] == remap[b]) + r = (r == ~0u) ? edges[i].prev : a; + + v = wedge[v]; + } while (v != a); + + // we only return the opposite vertex if it's unique + return r == a ? ~0u : r; +} + +static void fillFoldQuadrics(Quadric* vertex_quadrics, const EdgeAdjacency& adjacency, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* remap, const unsigned int* wedge) +{ + for (size_t i0 = 0; i0 < vertex_count; ++i0) + for (unsigned int i = adjacency.offsets[i0]; i < adjacency.offsets[i0 + 1]; ++i) + { + unsigned int i1 = adjacency.data[i].next; + unsigned int i2 = adjacency.data[i].prev; + + // since we only process paired edges each edge should occur twice; skip redundant edges + if (remap[i1] > remap[i0]) + continue; + + unsigned int i3 = oppositeVertex(adjacency, i1, unsigned(i0), remap, wedge); + + if (i3 == ~0u) + continue; + + Vector3 p10 = {vertex_positions[i1].x - vertex_positions[i0].x, vertex_positions[i1].y - vertex_positions[i0].y, vertex_positions[i1].z - vertex_positions[i0].z}; + Vector3 p20 = {vertex_positions[i2].x - vertex_positions[i0].x, vertex_positions[i2].y - vertex_positions[i0].y, vertex_positions[i2].z - vertex_positions[i0].z}; + Vector3 p30 = {vertex_positions[i3].x - vertex_positions[i0].x, vertex_positions[i3].y - vertex_positions[i0].y, vertex_positions[i3].z - vertex_positions[i0].z}; + + Vector3 normal = {p10.y * p20.z - p10.z * p20.y, p10.z * p20.x - p10.x * p20.z, p10.x * p20.y - p10.y * p20.x}; + Vector3 opposite = {p30.y * p10.z - p30.z * p10.y, p30.z * p10.x - p30.x * p10.z, p30.x * p10.y - p30.y * p10.x}; + float nl = sqrtf(normal.x * normal.x + normal.y * normal.y + normal.z * normal.z); + float ol = sqrtf(opposite.x * opposite.x + opposite.y * opposite.y + opposite.z * opposite.z); + + // identify folds by checking dihedral angle against a somewhat arbitrary ~30 degree cutoff + if (normal.x * opposite.x + normal.y * opposite.y + normal.z * opposite.z >= -0.85f * nl * ol) + continue; + + // use both triangles for symmetry since we won't visit this edge again + Quadric Q, QO; + quadricFromTriangleEdge(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f); + quadricFromTriangleEdge(QO, vertex_positions[i0], vertex_positions[i1], vertex_positions[i3], 1.f); + quadricAdd(Q, QO); + + quadricAdd(vertex_quadrics[remap[i0]], Q); + quadricAdd(vertex_quadrics[remap[i1]], Q); + } +} + static void fillAttributeQuadrics(Quadric* attribute_quadrics, QuadricGrad* attribute_gradients, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const float* vertex_attributes, size_t attribute_count) { for (size_t i = 0; i < index_count; i += 3) @@ -1271,6 +1337,20 @@ static float getNeighborhoodRadius(const EdgeAdjacency& adjacency, const Vector3 return sqrtf(result); } +static unsigned int getComplexTarget(unsigned int v, unsigned int target, const unsigned int* remap, const unsigned int* loop, const unsigned int* loopback) +{ + unsigned int r = remap[target]; + + // use loop metadata to guide complex collapses towards the correct wedge + // this works for edges on attribute discontinuities because loop/loopback track the single half-edge without a pair, similar to seams + if (loop[v] != ~0u && remap[loop[v]] == r) + return loop[v]; + else if (loopback[v] != ~0u && remap[loopback[v]] == r) + return loopback[v]; + else + return target; +} + static size_t boundEdgeCollapses(const EdgeAdjacency& adjacency, size_t vertex_count, size_t index_count, unsigned char* vertex_kind) { size_t dual_count = 0; @@ -1393,15 +1473,22 @@ static void rankEdgeCollapses(Collapse* collapses, size_t collapse_count, const } else { - // complex edges can have multiple wedges, so we need to aggregate errors for all wedges - // this is different from seams (where we aggregate pairwise) because all wedges collapse onto the same target + // complex edges can have multiple wedges, so we need to aggregate errors for all wedges based on the selected target if (vertex_kind[i0] == Kind_Complex) for (unsigned int v = wedge[i0]; v != i0; v = wedge[v]) - ei += quadricError(attribute_quadrics[v], &attribute_gradients[v * attribute_count], attribute_count, vertex_positions[i1], &vertex_attributes[i1 * attribute_count]); + { + unsigned int t = getComplexTarget(v, i1, remap, loop, loopback); + + ei += quadricError(attribute_quadrics[v], &attribute_gradients[v * attribute_count], attribute_count, vertex_positions[t], &vertex_attributes[t * attribute_count]); + } if (vertex_kind[i1] == Kind_Complex && bidi) for (unsigned int v = wedge[i1]; v != i1; v = wedge[v]) - ej += quadricError(attribute_quadrics[v], &attribute_gradients[v * attribute_count], attribute_count, vertex_positions[i0], &vertex_attributes[i0 * attribute_count]); + { + unsigned int t = getComplexTarget(v, i0, remap, loop, loopback); + + ej += quadricError(attribute_quadrics[v], &attribute_gradients[v * attribute_count], attribute_count, vertex_positions[t], &vertex_attributes[t * attribute_count]); + } } } @@ -1553,7 +1640,9 @@ static size_t performEdgeCollapses(unsigned int* collapse_remap, unsigned char* do { - collapse_remap[v] = i1; + unsigned int t = getComplexTarget(v, i1, remap, loop, loopback); + + collapse_remap[v] = t; v = wedge[v]; } while (v != i0); } @@ -2323,7 +2412,6 @@ size_t meshopt_simplifyEdge(unsigned int* destination, const unsigned int* indic assert(vertex_positions_stride % sizeof(float) == 0); assert(target_index_count <= index_count); assert(target_error >= 0); - assert((options & ~(meshopt_SimplifyLockBorder | meshopt_SimplifySparse | meshopt_SimplifyErrorAbsolute | meshopt_SimplifyPrune | meshopt_SimplifyRegularize | meshopt_SimplifyPermissive | meshopt_SimplifyInternalSolve | meshopt_SimplifyInternalDebug)) == 0); assert(vertex_attributes_stride >= attribute_count * sizeof(float) && vertex_attributes_stride <= 256); assert(vertex_attributes_stride % sizeof(float) == 0); assert(attribute_count <= kMaxAttributes); @@ -2417,9 +2505,12 @@ size_t meshopt_simplifyEdge(unsigned int* destination, const unsigned int* indic } fillFaceQuadrics(vertex_quadrics, volume_gradients, result, index_count, vertex_positions, remap); - fillVertexQuadrics(vertex_quadrics, vertex_positions, vertex_count, remap, options); + fillVertexQuadrics(vertex_quadrics, vertex_positions, vertex_count, remap, vertex_lock, sparse_remap, options); fillEdgeQuadrics(vertex_quadrics, result, index_count, vertex_positions, remap, vertex_kind, loop, loopback); + if (options & meshopt_SimplifyPreserveFolds) + fillFoldQuadrics(vertex_quadrics, adjacency, vertex_positions, vertex_count, remap, wedge); + if (attribute_count) fillAttributeQuadrics(attribute_quadrics, attribute_gradients, result, index_count, vertex_positions, vertex_attributes, attribute_count); @@ -2626,6 +2717,7 @@ size_t meshopt_simplifySloppy(unsigned int* destination, const unsigned int* ind assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); assert(vertex_positions_stride % sizeof(float) == 0); assert(target_index_count <= index_count); + assert(target_error >= 0); // we expect to get ~2 triangles/vertex in the output size_t target_cell_count = target_index_count / 6; @@ -2824,7 +2916,7 @@ size_t meshopt_simplifyPoints(unsigned int* destination, const float* vertex_pos size_t min_vertices = 0; size_t max_vertices = vertex_count; - // instead of starting in the middle, let's guess as to what the answer might be! triangle count usually grows as a square of grid size... + // instead of starting in the middle, let's guess as to what the answer might be! surface point count usually grows as a square of grid size... int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f); for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass) diff --git a/3rdparty/meshoptimizer/src/tangentspace.cpp b/3rdparty/meshoptimizer/src/tangentspace.cpp new file mode 100644 index 000000000..f9a60e97a --- /dev/null +++ b/3rdparty/meshoptimizer/src/tangentspace.cpp @@ -0,0 +1,730 @@ +// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details +#include "meshoptimizer.h" + +#include +#include +#include + +// This work is based on: +// Morten Mikkelsen. Simulation of wrinkled surfaces revisited. 2008 +// Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003 +// Cecil Hastings Jr. Approximations for digital computers. 1955 +namespace meshopt +{ + +struct Tangent +{ + float x, y, z, w; +}; + +static void computeFaceTangents(Tangent* result, size_t triangle_count, const unsigned int* indices, const float* vertex_positions, size_t vertex_positions_stride, const float* vertex_uvs, size_t vertex_uvs_stride) +{ + size_t vertex_position_stride_float = vertex_positions_stride / sizeof(float); + size_t vertex_uv_stride_float = vertex_uvs_stride / sizeof(float); + + for (size_t i = 0; i < triangle_count; ++i) + { + unsigned int a = indices ? indices[i * 3 + 0] : unsigned(i * 3 + 0); + unsigned int b = indices ? indices[i * 3 + 1] : unsigned(i * 3 + 1); + unsigned int c = indices ? indices[i * 3 + 2] : unsigned(i * 3 + 2); + + const float* pa = &vertex_positions[a * vertex_position_stride_float]; + const float* pb = &vertex_positions[b * vertex_position_stride_float]; + const float* pc = &vertex_positions[c * vertex_position_stride_float]; + + const float* ta = &vertex_uvs[a * vertex_uv_stride_float]; + const float* tb = &vertex_uvs[b * vertex_uv_stride_float]; + const float* tc = &vertex_uvs[c * vertex_uv_stride_float]; + + // compute face tangents using chain rule, see eq. 42 + float dp1x = pb[0] - pa[0], dp1y = pb[1] - pa[1], dp1z = pb[2] - pa[2]; + float dp2x = pc[0] - pa[0], dp2y = pc[1] - pa[1], dp2z = pc[2] - pa[2]; + float dt1x = tb[0] - ta[0], dt1y = tb[1] - ta[1]; + float dt2x = tc[0] - ta[0], dt2y = tc[1] - ta[1]; + + float rx = dt2y * dp1x - dt1y * dp2x; + float ry = dt2y * dp1y - dt1y * dp2y; + float rz = dt2y * dp1z - dt1y * dp2z; + + // .w stores orientation (+1/-1) or 0 for degenerate triangles, see eq. 48 + float det = dt1x * dt2y - dt1y * dt2x; + float rw = det == 0.f ? 0.f : (det > 0.f ? 1.f : -1.f); + + // treat zero area triangles as degenerate even if their UVs aren't for consistency + rw = (pa[0] == pb[0] && pa[1] == pb[1] && pa[2] == pb[2]) ? 0.f : rw; + rw = (pa[0] == pc[0] && pa[1] == pc[1] && pa[2] == pc[2]) ? 0.f : rw; + rw = (pb[0] == pc[0] && pb[1] == pc[1] && pb[2] == pc[2]) ? 0.f : rw; + + float rl = sqrtf(rx * rx + ry * ry + rz * rz); + float rs = rl != 0.f ? rw / rl : 0.f; + + result[i].x = rx * rs; + result[i].y = ry * rs; + result[i].z = rz * rs; + result[i].w = rw; + } +} + +struct Normal +{ + float x, y, z; + int id; +}; + +static void computeFaceNormals(Normal* result, size_t triangle_count, const unsigned int* indices, const float* vertex_positions, size_t vertex_positions_stride) +{ + size_t vertex_position_stride_float = vertex_positions_stride / sizeof(float); + + for (size_t i = 0; i < triangle_count; ++i) + { + unsigned int a = indices ? indices[i * 3 + 0] : unsigned(i * 3 + 0); + unsigned int b = indices ? indices[i * 3 + 1] : unsigned(i * 3 + 1); + unsigned int c = indices ? indices[i * 3 + 2] : unsigned(i * 3 + 2); + + const float* pa = &vertex_positions[a * vertex_position_stride_float]; + const float* pb = &vertex_positions[b * vertex_position_stride_float]; + const float* pc = &vertex_positions[c * vertex_position_stride_float]; + + float dp1x = pb[0] - pa[0], dp1y = pb[1] - pa[1], dp1z = pb[2] - pa[2]; + float dp2x = pc[0] - pa[0], dp2y = pc[1] - pa[1], dp2z = pc[2] - pa[2]; + + float nx = dp1y * dp2z - dp1z * dp2y; + float ny = dp1z * dp2x - dp1x * dp2z; + float nz = dp1x * dp2y - dp1y * dp2x; + + float nl = sqrtf(nx * nx + ny * ny + nz * nz); + float ns = nl != 0.f ? 1.f / nl : 0.f; + + result[i].x = nx * ns; + result[i].y = ny * ns; + result[i].z = nz * ns; + result[i].id = nl != 0.f; + } +} + +template +struct VertexHasherF +{ + const float* vertex_positions; + size_t vertex_positions_stride_float; + const float* vertex_normals; + size_t vertex_normals_stride_float; + const float* vertex_uvs; + size_t vertex_uvs_stride_float; + + size_t hash(unsigned int index) const + { + const unsigned int* p = reinterpret_cast(vertex_positions + index * vertex_positions_stride_float); + + unsigned int x = p[0], y = p[1], z = p[2]; + + // replace negative zero with zero + x = (x == 0x80000000) ? 0 : x; + y = (y == 0x80000000) ? 0 : y; + z = (z == 0x80000000) ? 0 : z; + + // scramble bits to make sure that integer coordinates have entropy in lower bits + x ^= x >> 17; + y ^= y >> 17; + z ^= z >> 17; + + if (PositionOnly) + return (x * 73856093) ^ (y * 19349663) ^ (z * 83492791); + + const unsigned int* n = reinterpret_cast(vertex_normals + index * vertex_normals_stride_float); + const unsigned int* t = reinterpret_cast(vertex_uvs + index * vertex_uvs_stride_float); + + // mix in normal bits + x ^= (n[0] == 0x80000000) ? 0 : n[0] >> 15; + y ^= (n[1] == 0x80000000) ? 0 : n[1] >> 15; + z ^= (n[2] == 0x80000000) ? 0 : n[2] >> 15; + + // collect texture coordinate bits (simplified -0 handling and scrambling) + unsigned int w = (t[0] ^ t[1]) & 0x7fffffff; + w ^= w >> 13; + + // Optimized Spatial Hashing for Collision Detection of Deformable Objects + return (x * 73856093) ^ (y * 19349663) ^ (z * 83492791) ^ (w * 50331653); + } + + bool equal(unsigned int lhs, unsigned int rhs) const + { + const float* lp = vertex_positions + lhs * vertex_positions_stride_float; + const float* rp = vertex_positions + rhs * vertex_positions_stride_float; + + if (PositionOnly) + return lp[0] == rp[0] && lp[1] == rp[1] && lp[2] == rp[2]; + + const float* ln = vertex_normals + lhs * vertex_normals_stride_float; + const float* lt = vertex_uvs + lhs * vertex_uvs_stride_float; + const float* rn = vertex_normals + rhs * vertex_normals_stride_float; + const float* rt = vertex_uvs + rhs * vertex_uvs_stride_float; + + return lp[0] == rp[0] && lp[1] == rp[1] && lp[2] == rp[2] && ln[0] == rn[0] && ln[1] == rn[1] && ln[2] == rn[2] && lt[0] == rt[0] && lt[1] == rt[1]; + } +}; + +static size_t hashBuckets4(size_t count) +{ + size_t buckets = 1; + while (buckets < count + count / 4) + buckets *= 2; + + return buckets; +} + +template +static T* hashLookup4(T* table, size_t buckets, const Hash& hash, const T& key, const T& empty) +{ + assert(buckets > 0); + assert((buckets & (buckets - 1)) == 0); + + size_t hashmod = buckets - 1; + size_t bucket = hash.hash(key) & hashmod; + + for (size_t probe = 0; probe <= hashmod; ++probe) + { + T& item = table[bucket]; + + if (item == empty) + return &item; + + if (hash.equal(item, key)) + return &item; + + // hash collision, quadratic probing + bucket = (bucket + probe + 1) & hashmod; + } + + assert(false && "Hash table is full"); // unreachable + return NULL; +} + +template +static void buildVertexRemap(unsigned int* remap, size_t vertex_count, const float* vertex_positions, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, const float* vertex_uvs, size_t vertex_uvs_stride, meshopt_Allocator& allocator) +{ + VertexHasherF vertex_hasher = {vertex_positions, vertex_positions_stride / sizeof(float), vertex_normals, vertex_normals_stride / sizeof(float), vertex_uvs, vertex_uvs_stride / sizeof(float)}; + + size_t vertex_table_size = hashBuckets4(vertex_count); + unsigned int* vertex_table = allocator.allocate(vertex_table_size); + memset(vertex_table, -1, vertex_table_size * sizeof(unsigned int)); + + for (size_t i = 0; i < vertex_count; ++i) + { + unsigned int index = unsigned(i); + unsigned int* entry = hashLookup4(vertex_table, vertex_table_size, vertex_hasher, index, ~0u); + + if (*entry == ~0u) + *entry = index; + + remap[index] = *entry; + } + + allocator.deallocate(vertex_table); +} + +struct CornerAdjacency +{ + unsigned int* offsets; + unsigned int* data; +}; + +static void buildCornerAdjacency(CornerAdjacency& adjacency, const unsigned int* indices, size_t index_count, const unsigned int* remap, size_t vertex_count, meshopt_Allocator& allocator) +{ + adjacency.offsets = allocator.allocate(vertex_count + 1); + adjacency.data = allocator.allocate(index_count); + + size_t face_count = index_count / 3; + unsigned int* offsets = adjacency.offsets + 1; + + // fill corner counts + memset(offsets, 0, vertex_count * sizeof(unsigned int)); + + for (size_t i = 0; i < index_count; ++i) + { + unsigned int v = indices ? remap[indices[i]] : remap[i]; + assert(v < vertex_count); + + offsets[v]++; + } + + // fill offset table + unsigned int offset = 0; + + for (size_t i = 0; i < vertex_count; ++i) + { + unsigned int count = offsets[i]; + offsets[i] = offset; + offset += count; + } + + assert(offset == index_count); + + // fill corner data + for (size_t i = 0; i < face_count; ++i) + { + unsigned int a = indices ? remap[indices[i * 3 + 0]] : remap[i * 3 + 0]; + unsigned int b = indices ? remap[indices[i * 3 + 1]] : remap[i * 3 + 1]; + unsigned int c = indices ? remap[indices[i * 3 + 2]] : remap[i * 3 + 2]; + + // encode corner index in the low 2 bits + adjacency.data[offsets[a]++] = unsigned(i << 2) | 0; + adjacency.data[offsets[b]++] = unsigned(i << 2) | 1; + adjacency.data[offsets[c]++] = unsigned(i << 2) | 2; + } + + // finalize offsets + adjacency.offsets[0] = 0; + assert(adjacency.offsets[vertex_count] == index_count); +} + +static unsigned int follow2(unsigned int* parents, unsigned int index) +{ + while (index != parents[index]) + { + unsigned int parent = parents[index]; + parents[index] = parents[parent]; + index = parent; + } + + return index; +} + +static void mergeTangentGroups(unsigned int* groups, unsigned int* facegroups, unsigned char* facesign, const unsigned int* data, size_t count, const unsigned int* indices, const unsigned int* remap) +{ + static const int next[4] = {1, 2, 0, 1}; + + for (size_t i = 0; i < count; ++i) + { + unsigned int ti = data[i] >> 2; + unsigned int cib = ti * 3 + next[(data[i] & 3) + 0]; + unsigned int cic = ti * 3 + next[(data[i] & 3) + 1]; + + unsigned int nib = indices ? remap[indices[cib]] : remap[cib]; + unsigned int nic = indices ? remap[indices[cic]] : remap[cic]; + + for (size_t j = i + 1; j < count; ++j) + { + unsigned int tj = data[j] >> 2; + unsigned int cjb = tj * 3 + next[(data[j] & 3) + 0]; + unsigned int cjc = tj * 3 + next[(data[j] & 3) + 1]; + + unsigned int njb = indices ? remap[indices[cjb]] : remap[cjb]; + unsigned int njc = indices ? remap[indices[cjc]] : remap[cjc]; + + // merge tangent groups if triangles are adjacent and orientations agree or either one is degenerate + if ((njb == nic || njc == nib) && (facesign[ti] | facesign[tj]) != 3) + { + // if either triangle is degenerate, it could have been assigned to a concrete orientation; we need to confirm that via facegroups + if ((facesign[ti] & facesign[tj]) == 0) + { + unsigned int gti = follow2(facegroups, ti); + unsigned int gtj = follow2(facegroups, tj); + + // skip merge if orientations diverge; otherwise, union the groups with gti as the root + if (gti != gtj) + { + if ((facesign[gti] | facesign[gtj]) == 3) + continue; + + facegroups[gtj] = gti; + facesign[gti] |= facesign[gtj]; + } + } + + // union tangent groups for individual corners with gi as the root; the orientations must match due to facegroups check above + unsigned int gi = follow2(groups, ti * 3 + (data[i] & 3)); + unsigned int gj = follow2(groups, tj * 3 + (data[j] & 3)); + + if (gi != gj) + groups[gj] = gi; + } + } + } +} + +static void mergeNormalGroups(unsigned int* groups, const unsigned int* data, size_t count, const unsigned int* indices, const unsigned int* remap, const Normal* face_normals, float crease_cutoff) +{ + static const int next[4] = {1, 2, 0, 1}; + + for (size_t i = 0; i < count; ++i) + { + unsigned int ti = data[i] >> 2; + unsigned int cib = ti * 3 + next[(data[i] & 3) + 0]; + unsigned int cic = ti * 3 + next[(data[i] & 3) + 1]; + + unsigned int nib = indices ? remap[indices[cib]] : remap[cib]; + unsigned int nic = indices ? remap[indices[cic]] : remap[cic]; + + const Normal& fi = face_normals[ti]; + + for (size_t j = i + 1; j < count; ++j) + { + unsigned int tj = data[j] >> 2; + unsigned int cjb = tj * 3 + next[(data[j] & 3) + 0]; + unsigned int cjc = tj * 3 + next[(data[j] & 3) + 1]; + + unsigned int njb = indices ? remap[indices[cjb]] : remap[cjb]; + unsigned int njc = indices ? remap[indices[cjc]] : remap[cjc]; + + const Normal& fj = face_normals[tj]; + + // merge normal groups if triangles are adjacent and normal angle is below crease threshold + if ((njb == nic || njc == nib) && fi.x * fj.x + fi.y * fj.y + fi.z * fj.z > crease_cutoff && fi.id == fj.id) + { + // union normal groups for individual corners with gi as the root + unsigned int gi = follow2(groups, ti * 3 + (data[i] & 3)); + unsigned int gj = follow2(groups, tj * 3 + (data[j] & 3)); + + if (gi != gj) + groups[gj] = gi; + } + } + } +} + +inline float optacos(float x) +{ + // approximation for acos(x) = sqrt(1 - abs(x)) * 2-degree polynomial, max error ~5e-4 + float ax = fabsf(x); + ax = ax > 1.f ? 1.f : ax; + float r = 1.570337f + ax * (-0.2053972f + ax * 0.05147786f); + r *= sqrtf(1.0f - ax); + return x < 0.f ? 3.1415926f - r : r; +} + +static void accumulateTangents(float* result, const unsigned int* groups, const unsigned int* indices, size_t index_count, const Tangent* face_tangents, const float* vertex_positions, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, unsigned int options) +{ + static const int next[4] = {1, 2, 0, 1}; + + size_t vertex_position_stride_float = vertex_positions_stride / sizeof(float); + size_t vertex_normal_stride_float = vertex_normals_stride / sizeof(float); + + size_t face_count = index_count / 3; + + for (size_t i = 0; i < face_count; ++i) + { + const Tangent& t = face_tangents[i]; + + if (t.w == 0.f) + continue; + + for (int k = 0; k < 3; ++k) + { + float* r = &result[size_t(groups[i * 3 + k]) * 4]; + + unsigned int a = indices ? indices[i * 3 + k] : unsigned(i * 3 + k); + unsigned int b = indices ? indices[i * 3 + next[k]] : unsigned(i * 3 + next[k]); + unsigned int c = indices ? indices[i * 3 + next[k + 1]] : unsigned(i * 3 + next[k + 1]); + + const float* pa = vertex_positions + a * vertex_position_stride_float; + const float* pb = vertex_positions + b * vertex_position_stride_float; + const float* pc = vertex_positions + c * vertex_position_stride_float; + + const float* n = vertex_normals + a * vertex_normal_stride_float; + + // reproject tangent vector onto vertex normal + float sd = t.x * n[0] + t.y * n[1] + t.z * n[2]; + float sx = t.x - n[0] * sd, sy = t.y - n[1] * sd, sz = t.z - n[2] * sd; + float sl = sqrtf(sx * sx + sy * sy + sz * sz); + + // compute incident angle for weighting, in tangent plane + // note: this step is absent from the paper, reference implementation computes angle after projecting edges to tangent plane... + float dp1x = pb[0] - pa[0], dp1y = pb[1] - pa[1], dp1z = pb[2] - pa[2]; + float dp2x = pc[0] - pa[0], dp2y = pc[1] - pa[1], dp2z = pc[2] - pa[2]; + + float dp1d = dp1x * n[0] + dp1y * n[1] + dp1z * n[2]; + float dp2d = dp2x * n[0] + dp2y * n[1] + dp2z * n[2]; + + dp1x -= n[0] * dp1d, dp1y -= n[1] * dp1d, dp1z -= n[2] * dp1d; + dp2x -= n[0] * dp2d, dp2y -= n[1] * dp2d, dp2z -= n[2] * dp2d; + + float dp1l = dp1x * dp1x + dp1y * dp1y + dp1z * dp1z; + float dp2l = dp2x * dp2x + dp2y * dp2y + dp2z * dp2z; + float dpl = sqrtf(dp1l * dp2l); + + float cosa = (dp1x * dp2x + dp1y * dp2y + dp1z * dp2z) * (dpl == 0.f ? 0.f : 1.f / dpl); + float angle = optacos(cosa); // optacos handles clamping to [-1..1] + + // accumulate renormalized tangent weighted by angle + float w = angle * (sl == 0.f ? 0.f : 1.f / sl); + + // weight larger adjacent triangles more to reduce interpolation artifacts from slivers; this deviates from reference implementation for higher quality + if ((options & meshopt_TangentCompatible) == 0) + w *= dpl; + + r[0] += sx * w; + r[1] += sy * w; + r[2] += sz * w; + } + } +} + +static void accumulateNormals(float* result, const unsigned int* groups, const unsigned int* indices, size_t index_count, const Normal* face_normals, const float* vertex_positions, size_t vertex_positions_stride) +{ + static const int next[4] = {1, 2, 0, 1}; + + size_t vertex_position_stride_float = vertex_positions_stride / sizeof(float); + + size_t face_count = index_count / 3; + + for (size_t i = 0; i < face_count; ++i) + { + const Normal& n = face_normals[i]; + + for (int k = 0; k < 3; ++k) + { + float* r = &result[size_t(groups[i * 3 + k]) * 3]; + + unsigned int a = indices ? indices[i * 3 + k] : unsigned(i * 3 + k); + unsigned int b = indices ? indices[i * 3 + next[k]] : unsigned(i * 3 + next[k]); + unsigned int c = indices ? indices[i * 3 + next[k + 1]] : unsigned(i * 3 + next[k + 1]); + + const float* pa = vertex_positions + a * vertex_position_stride_float; + const float* pb = vertex_positions + b * vertex_position_stride_float; + const float* pc = vertex_positions + c * vertex_position_stride_float; + + // compute incident angle for weighting + float dp1x = pb[0] - pa[0], dp1y = pb[1] - pa[1], dp1z = pb[2] - pa[2]; + float dp2x = pc[0] - pa[0], dp2y = pc[1] - pa[1], dp2z = pc[2] - pa[2]; + + float dp1l = dp1x * dp1x + dp1y * dp1y + dp1z * dp1z; + float dp2l = dp2x * dp2x + dp2y * dp2y + dp2z * dp2z; + float dpl = sqrtf(dp1l * dp2l); + + float cosa = (dp1x * dp2x + dp1y * dp2y + dp1z * dp2z) * (dpl == 0.f ? 0.f : 1.f / dpl); + float angle = optacos(cosa); // optacos handles clamping to [-1..1] + + // accumulate normal weighted by angle and edge length product; this weighting scheme matches the default scheme used for tangents + float w = angle * dpl; + + r[0] += n.x * w; + r[1] += n.y * w; + r[2] += n.z * w; + } + } +} + +static void smoothNormals(float* result, float* scratch, const unsigned int* groups, size_t index_count, float alpha) +{ + static const int next[4] = {1, 2, 0, 1}; + + size_t face_count = index_count / 3; + + memset(scratch, 0, index_count * 4 * sizeof(float)); + + // for each triangle, accumulate normal deltas alongside each edge in both directions + for (size_t i = 0; i < face_count; ++i) + { + for (int k = 0; k < 3; ++k) + { + unsigned int ga = groups[i * 3 + k]; + unsigned int gb = groups[i * 3 + next[k]]; + + const float* na = &result[ga * 3]; + const float* nb = &result[gb * 3]; + + // normal alignment is symmetric; we compute it once per edge and use dp^2 for stronger alignment + float dp = na[0] * nb[0] + na[1] * nb[1] + na[2] * nb[2]; + float w = (dp > 0.f ? dp * dp : 0.f); + + float nx = (nb[0] - na[0]) * w, ny = (nb[1] - na[1]) * w, nz = (nb[2] - na[2]) * w; + + // accumulate deltas; the last component is used to compute delta averages below + float* sa = &scratch[ga * 4]; + float* sb = &scratch[gb * 4]; + + sa[0] += nx; + sa[1] += ny; + sa[2] += nz; + sa[3] += 1.f; + + sb[0] -= nx; + sb[1] -= ny; + sb[2] -= nz; + sb[3] += 1.f; + } + } + + // apply average delta to group roots and renormalize the results + for (size_t i = 0; i < index_count; ++i) + if (groups[i] == i && scratch[i * 4 + 3] > 0.f) + { + float* r = &result[i * 3]; + const float* s = &scratch[i * 4]; + + float sw = alpha / s[3]; + float rx = r[0] + s[0] * sw, ry = r[1] + s[1] * sw, rz = r[2] + s[2] * sw; + + float rl = sqrtf(rx * rx + ry * ry + rz * rz); + float rs = rl == 0.f ? 0.f : 1.f / rl; + + r[0] = rx * rs; + r[1] = ry * rs; + r[2] = rz * rs; + } +} + +} // namespace meshopt + +void meshopt_generateTangents(float* result, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_normals, size_t vertex_normals_stride, const float* vertex_uvs, size_t vertex_uvs_stride, unsigned int options) +{ + using namespace meshopt; + + assert(indices || index_count == vertex_count); + assert(index_count % 3 == 0); + assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); + assert(vertex_normals_stride >= 12 && vertex_normals_stride <= 256); + assert(vertex_uvs_stride >= 8 && vertex_uvs_stride <= 256); + assert(vertex_positions_stride % sizeof(float) == 0 && vertex_normals_stride % sizeof(float) == 0 && vertex_uvs_stride % sizeof(float) == 0); + + meshopt_Allocator allocator; + + size_t face_count = index_count / 3; + + // compute vertex remap to unique vertex index + unsigned int* remap = allocator.allocate(vertex_count); + buildVertexRemap(remap, vertex_count, vertex_positions, vertex_positions_stride, vertex_normals, vertex_normals_stride, vertex_uvs, vertex_uvs_stride, allocator); + + // build adjacency information + CornerAdjacency adjacency = {}; + buildCornerAdjacency(adjacency, indices, index_count, remap, vertex_count, allocator); + + // compute per-triangle tangents + Tangent* face_tangents = allocator.allocate(face_count); + computeFaceTangents(face_tangents, face_count, indices, vertex_positions, vertex_positions_stride, vertex_uvs, vertex_uvs_stride); + + // compute per-corner tangent groups: triangles adjacent to the same vertex with the same orientation are merged in strips + unsigned int* groups = allocator.allocate(index_count); + for (size_t i = 0; i < index_count; ++i) + groups[i] = unsigned(i); + + // compute per-face orientation groups: degenerate triangles are merged with adjacent non-degenerate ones and inherit their orientation, stored as a bitmask + unsigned int* facegroups = allocator.allocate(face_count); + unsigned char* facesign = allocator.allocate(face_count); + + for (size_t i = 0; i < face_count; ++i) + { + facegroups[i] = unsigned(i); + facesign[i] = (face_tangents[i].w > 0.f ? 1 : 0) | (face_tangents[i].w < 0.f ? 2 : 0); + } + + for (size_t i = 0; i < vertex_count; ++i) + if (adjacency.offsets[i + 1] != adjacency.offsets[i]) + mergeTangentGroups(groups, facegroups, facesign, adjacency.data + adjacency.offsets[i], adjacency.offsets[i + 1] - adjacency.offsets[i], indices, remap); + + for (size_t i = 0; i < index_count; ++i) + groups[i] = follow2(groups, unsigned(i)); + + // accumulate tangents into their own respective groups + memset(result, 0, index_count * sizeof(float) * 4); + accumulateTangents(result, groups, indices, index_count, face_tangents, vertex_positions, vertex_positions_stride, vertex_normals, vertex_normals_stride, options); + + // finalize tangent signs for each face using facegroups + for (size_t i = 0; i < face_count; ++i) + { + unsigned int fg = follow2(facegroups, unsigned(i)); + float w = (facesign[fg] & 1) ? 1.f : -1.f; // note: if face is degenerate, we canonicalize to -1 for consistency + + result[(i * 3 + 0) * 4 + 3] = result[(i * 3 + 1) * 4 + 3] = result[(i * 3 + 2) * 4 + 3] = w; + } + + // finalize tangent vectors by normalizing roots and propagating the rest + for (size_t i = 0; i < index_count; ++i) + if (groups[i] == i) + { + float* r = &result[i * 4]; + float l = sqrtf(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]); + float s = l == 0.f ? 0.f : 1.f / l; + + r[0] *= s; + r[1] *= s; + r[2] *= s; + + // for isolated degenerate triangles, use (1, 0, 0) tangent for consistency + r[0] = (l == 0.f && (options & meshopt_TangentZeroFallback) == 0) ? 1.f : r[0]; + } + + for (size_t i = 0; i < index_count; ++i) + if (groups[i] != i) + memcpy(&result[i * 4], &result[size_t(groups[i]) * 4], sizeof(float) * 3); // .w was set per face earlier +} + +void meshopt_generateNormals(float* result, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float crease_angle, float smoothing) +{ + using namespace meshopt; + + assert(indices || index_count == vertex_count); + assert(index_count % 3 == 0); + assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256); + assert(vertex_positions_stride % sizeof(float) == 0); + assert(crease_angle >= 0 && crease_angle <= 3.1415927f); + assert(smoothing >= 0); + + meshopt_Allocator allocator; + + size_t face_count = index_count / 3; + float crease_cutoff = cosf(crease_angle); + + // compute vertex remap to unique vertex index + unsigned int* remap = allocator.allocate(vertex_count); + buildVertexRemap(remap, vertex_count, vertex_positions, vertex_positions_stride, NULL, 0, NULL, 0, allocator); + + // build adjacency information + CornerAdjacency adjacency = {}; + buildCornerAdjacency(adjacency, indices, index_count, remap, vertex_count, allocator); + + // compute per-triangle normals + Normal* face_normals = allocator.allocate(face_count); + computeFaceNormals(face_normals, face_count, indices, vertex_positions, vertex_positions_stride); + + // compute per-corner normal groups: triangles adjacent to the same vertex with similar normals are merged in strips + unsigned int* groups = allocator.allocate(index_count); + for (size_t i = 0; i < index_count; ++i) + groups[i] = unsigned(i); + + for (size_t i = 0; i < vertex_count; ++i) + if (adjacency.offsets[i + 1] != adjacency.offsets[i]) + mergeNormalGroups(groups, adjacency.data + adjacency.offsets[i], adjacency.offsets[i + 1] - adjacency.offsets[i], indices, remap, face_normals, crease_cutoff); + + for (size_t i = 0; i < index_count; ++i) + groups[i] = follow2(groups, unsigned(i)); + + // accumulate normals into their own respective groups + memset(result, 0, index_count * sizeof(float) * 3); + accumulateNormals(result, groups, indices, index_count, face_normals, vertex_positions, vertex_positions_stride); + + // normalize accumulated normals; we need to do this before smoothing as smoothing expects normalized inputs + for (size_t i = 0; i < index_count; ++i) + if (groups[i] == i) + { + float* r = &result[i * 3]; + float l = sqrtf(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]); + float s = l == 0.f ? 0.f : 1.f / l; + + r[0] *= s; + r[1] *= s; + r[2] *= s; + } + + if (smoothing > 0.f) + { + int passes = int(ceilf(smoothing)); + passes = passes < 10 ? passes : 10; // clamp to avoid pathological cases + + float* scratch = allocator.allocate(index_count * 4); + + for (int pass = 0; pass < passes; ++pass) + { + float rem = smoothing - float(pass); + float alpha = 0.5f * (rem < 1.f ? rem : 1.f); + + smoothNormals(result, scratch, groups, index_count, alpha); + } + } + + // finalize normal vectors by propagating the rest + for (size_t i = 0; i < index_count; ++i) + if (groups[i] != i) + memcpy(&result[i * 3], &result[size_t(groups[i]) * 3], sizeof(float) * 3); +} diff --git a/3rdparty/meshoptimizer/src/vertexcodec.cpp b/3rdparty/meshoptimizer/src/vertexcodec.cpp index 7085cce32..9aa131f63 100644 --- a/3rdparty/meshoptimizer/src/vertexcodec.cpp +++ b/3rdparty/meshoptimizer/src/vertexcodec.cpp @@ -7,19 +7,19 @@ // The block below auto-detects SIMD ISA that can be used on the target platform #ifndef MESHOPTIMIZER_NO_SIMD -// The SIMD implementation requires SSSE3, which can be enabled unconditionally through compiler settings -#if defined(__AVX__) || defined(__SSSE3__) +// The SIMD implementation requires SSSE3+POPCNT, which can be enabled unconditionally through compiler settings +#if defined(__AVX__) || (defined(__SSSE3__) && defined(__POPCNT__)) #define SIMD_SSE #endif // An experimental implementation using AVX512 instructions; it's only enabled when AVX512 is enabled through compiler settings -#if defined(__AVX512VBMI2__) && defined(__AVX512VBMI__) && defined(__AVX512VL__) && defined(__POPCNT__) +#if defined(__AVX512VBMI2__) && defined(__AVX512VBMI__) && defined(__AVX512VL__) #undef SIMD_SSE #define SIMD_AVX #endif -// MSVC supports compiling SSSE3 code regardless of compile options; we use a cpuid-based scalar fallback -#if !defined(SIMD_SSE) && !defined(SIMD_AVX) && defined(_MSC_VER) && !defined(__clang__) && (defined(_M_IX86) || defined(_M_X64)) +// MSVC supports compiling SSSE3+POPCNT code regardless of compile options; we use a cpuid-based scalar fallback +#if !defined(SIMD_SSE) && !defined(SIMD_AVX) && defined(_MSC_VER) && !defined(__clang__) && (defined(_M_IX86) || (defined(_M_X64) && !defined(_M_ARM64EC))) #define SIMD_SSE #define SIMD_FALLBACK #endif @@ -28,7 +28,7 @@ #if !defined(SIMD_SSE) && !defined(SIMD_AVX) && ((defined(__clang__) && __clang_major__ * 100 + __clang_minor__ >= 308) || (defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ >= 409)) && (defined(__i386__) || defined(__x86_64__)) #define SIMD_SSE #define SIMD_FALLBACK -#define SIMD_TARGET __attribute__((target("ssse3"))) +#define SIMD_TARGET __attribute__((target("ssse3,popcnt"))) #endif // GCC/clang define these when NEON support is available @@ -37,7 +37,7 @@ #endif // On MSVC, we assume that ARM builds always target NEON-capable devices -#if !defined(SIMD_NEON) && defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64)) +#if !defined(SIMD_NEON) && defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64) || defined(_M_ARM64EC)) #define SIMD_NEON #endif @@ -56,7 +56,7 @@ // When targeting AArch64/x64, optimize for latency to allow decoding of individual 16-byte groups to overlap // We don't do this for 32-bit systems because we need 64-bit math for this and this will hurt in-order CPUs -#if defined(__x86_64__) || defined(_M_X64) || defined(__aarch64__) || defined(_M_ARM64) +#if (defined(__x86_64__) || defined(_M_X64) || defined(__aarch64__) || defined(_M_ARM64) || defined(_M_ARM64EC)) && !defined(MESHOPTIMIZER_VERTEXCODEC_SIMDNOLOPT) #define SIMD_LATENCYOPT #endif @@ -72,7 +72,7 @@ #endif // !MESHOPTIMIZER_NO_SIMD #ifdef SIMD_SSE -#include +#include #endif #if defined(SIMD_SSE) && defined(SIMD_FALLBACK) @@ -708,6 +708,14 @@ static const unsigned char* decodeVertexBlock(const unsigned char* data, const u size_t vertex_count_aligned = (vertex_count + kByteGroupSize - 1) & ~(kByteGroupSize - 1); assert(vertex_count <= vertex_count_aligned); + // we can decode directly into the output buffer + // this uses strided writes and also reads the last vertex once, which is bad for performance for write-combined memory so we only enable this if configured +#ifdef MESHOPTIMIZER_VERTEXCODEC_ZEROCOPY + unsigned char* target = vertex_data; +#else + unsigned char* target = transposed; +#endif + size_t control_size = version == 0 ? 0 : vertex_size / 4; if (size_t(data_end - data) < control_size) return NULL; @@ -750,22 +758,23 @@ static const unsigned char* decodeVertexBlock(const unsigned char* data, const u switch (channel & 3) { case 0: - decodeDeltas1(buffer, transposed + k, vertex_count, vertex_size, last_vertex + k, 0); + decodeDeltas1(buffer, target + k, vertex_count, vertex_size, last_vertex + k, 0); break; case 1: - decodeDeltas1(buffer, transposed + k, vertex_count, vertex_size, last_vertex + k, 0); + decodeDeltas1(buffer, target + k, vertex_count, vertex_size, last_vertex + k, 0); break; case 2: - decodeDeltas1(buffer, transposed + k, vertex_count, vertex_size, last_vertex + k, (32 - (channel >> 4)) & 31); + decodeDeltas1(buffer, target + k, vertex_count, vertex_size, last_vertex + k, (32 - (channel >> 4)) & 31); break; default: return NULL; // invalid channel type } } - memcpy(vertex_data, transposed, vertex_count * vertex_size); + if (target == transposed) + memcpy(vertex_data, transposed, vertex_count * vertex_size); - memcpy(last_vertex, &transposed[vertex_size * (vertex_count - 1)], vertex_size); + memcpy(last_vertex, &target[vertex_size * (vertex_count - 1)], vertex_size); return data; } @@ -804,210 +813,145 @@ static bool gDecodeBytesGroupInitialized = decodeBytesGroupBuildTables(); #endif #ifdef SIMD_SSE +// sent mask, replicating shuffle, and two multipliers (even/odd) for multishift emulation +static const unsigned char kDecodeBytesGroupConfig[9][4][16] = { + {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128}, {0}, {0}}, + {{3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}, {0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3}, {4, 0, 64, 0, 4, 0, 64, 0, 4, 0, 64, 0, 4, 0, 64, 0}, {16, 0, 0, 1, 16, 0, 0, 1, 16, 0, 0, 1, 16, 0, 0, 1}}, + {{15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15}, {0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7}, {16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0}, {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}}, + {{0}, {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128}, {0}, {0}}, + {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128}, {0}, {0}}, + {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1}, {0, 1, 64, 0, 16, 0, 4, 0, 0, 1, 64, 0, 16, 0, 4, 0}, {128, 0, 32, 0, 8, 0, 2, 0, 128, 0, 32, 0, 8, 0, 2, 0}}, + {{3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}, {0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3}, {4, 0, 64, 0, 4, 0, 64, 0, 4, 0, 64, 0, 4, 0, 64, 0}, {16, 0, 0, 1, 16, 0, 0, 1, 16, 0, 0, 1, 16, 0, 0, 1}}, + {{15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15}, {0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7}, {16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0, 16, 0}, {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}}, + {{0}, {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128}, {0}, {0}}, +}; + SIMD_TARGET -inline __m128i decodeShuffleMask(unsigned char mask0, unsigned char mask1) +static inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) { + // 0 for 1-bit, 1 for 2-bit, 2 for 4-bit, 3 for 8-bit, and 4 for 0-bit as it makes some of the uses easier + static const int hbtn[9] = {4, 1, 2, 3, 4, 0, 1, 2, 3}; + + int n = hbtn[hbits]; + +#ifdef SIMD_LATENCYOPT + unsigned long long data64; + memcpy(&data64, data, 8); + data64 &= data64 >> n; + data64 &= data64 >> (n >> 1); + + // mask out one bit per group that is set if all group bits were 1 + static const unsigned long long lanes[9] = {0, 0x55555555, 0x1111111111111111ull, 0, 0, 0xffff, 0x55555555, 0x1111111111111111ull, 0}; + int datacnt = int(_mm_popcnt_u64(data64 & lanes[hbits])); +#endif + + // for 8-bit groups, instead of loading the bytes through 'data', we load them through 'skip' as they are easier to preserve + // for 0-bit groups, the load results get discarded because mask is always 0; in both cases the shift wraps to zero + const unsigned char* skip = data + ((2 << n) & 15); + + __m128i selb = _mm_loadl_epi64(reinterpret_cast(data)); + __m128i rest = _mm_loadu_si128(reinterpret_cast(skip)); + + // unpack 1, 2 or 4-bit values: shuffle replicates each source byte into both halves of a 16-bit lane + // mulhi extracts even and odd fields into the low byte; the results are interleaved back with shift/or + __m128i selw = _mm_shuffle_epi8(selb, _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][1]))); + __m128i sel0 = _mm_mulhi_epu16(selw, _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][2]))); + __m128i sel1 = _mm_mulhi_epu16(selw, _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][3]))); + __m128i seli = _mm_or_si128(sel0, _mm_slli_epi16(sel1, 8)); + + // the interleaved fields are masked by the bit count (special handling: for 0/8-bit values, mul produces 0) + __m128i sent = _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][0])); + __m128i sel = _mm_and_si128(seli, sent); + + // compare sel to sentinel; returns 0 for 0-bit (mul produces 0, sent is 1), 1 for 8-bit (mul produces 0, sent is 0) + __m128i mask = _mm_cmpeq_epi8(sel, sent); + int mask16 = _mm_movemask_epi8(mask); + unsigned char mask0 = (unsigned char)(mask16 & 255); + unsigned char mask1 = (unsigned char)(mask16 >> 8); + + // decode shuffle mask from two halves; second half needs to be shifted by popcount(mask0) __m128i sm0 = _mm_loadl_epi64(reinterpret_cast(&kDecodeBytesGroupShuffle[mask0])); __m128i sm1 = _mm_loadl_epi64(reinterpret_cast(&kDecodeBytesGroupShuffle[mask1])); - __m128i sm1off = _mm_set1_epi8(kDecodeBytesGroupCount[mask0]); - __m128i sm1r = _mm_add_epi8(sm1, sm1off); + // each lane of mask is 0x00 or 0xff; sad yields 255*popcount(mask0) in low word => low byte is -popcount(mask0) + __m128i npops = _mm_sad_epu8(mask, _mm_setzero_si128()); + __m128i sm1r = _mm_sub_epi8(sm1, _mm_shuffle_epi8(npops, _mm_setzero_si128())); + __m128i shuf = _mm_unpacklo_epi64(sm0, sm1r); - return _mm_unpacklo_epi64(sm0, sm1r); -} + // expand rest via shuffle mask and combine with sel; shuffle mask zeroes out bytes that are replaced by sel + __m128i result = _mm_or_si128(_mm_shuffle_epi8(rest, shuf), _mm_andnot_si128(mask, sel)); -SIMD_TARGET -inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) -{ - switch (hbits) - { - case 0: - case 4: - { - __m128i result = _mm_setzero_si128(); + _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - - return data; - } - - case 1: - case 6: - { -#ifdef __GNUC__ - typedef int __attribute__((aligned(1))) unaligned_int; +#ifdef SIMD_LATENCYOPT + // datacnt is 0 for 8-bit groups so we can't use skip to advance; 0-bit groups wrap the shift to zero + return data + ((2 << n) & 31) + datacnt; #else - typedef int unaligned_int; + return skip + _mm_popcnt_u32(mask16); #endif - -#ifdef SIMD_LATENCYOPT - unsigned int data32; - memcpy(&data32, data, 4); - data32 &= data32 >> 1; - - // arrange bits such that low bits of nibbles of data64 contain all 2-bit elements of data32 - unsigned long long data64 = ((unsigned long long)data32 << 30) | (data32 & 0x3fffffff); - - // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used mode 3 - int datacnt = int(((data64 & 0x1111111111111111ull) * 0x1111111111111111ull) >> 60); -#endif - - __m128i sel2 = _mm_cvtsi32_si128(*reinterpret_cast(data)); - __m128i rest = _mm_loadu_si128(reinterpret_cast(data + 4)); - - __m128i sel22 = _mm_unpacklo_epi8(_mm_srli_epi16(sel2, 4), sel2); - __m128i sel2222 = _mm_unpacklo_epi8(_mm_srli_epi16(sel22, 2), sel22); - __m128i sel = _mm_and_si128(sel2222, _mm_set1_epi8(3)); - - __m128i mask = _mm_cmpeq_epi8(sel, _mm_set1_epi8(3)); - int mask16 = _mm_movemask_epi8(mask); - unsigned char mask0 = (unsigned char)(mask16 & 255); - unsigned char mask1 = (unsigned char)(mask16 >> 8); - - __m128i shuf = decodeShuffleMask(mask0, mask1); - __m128i result = _mm_or_si128(_mm_shuffle_epi8(rest, shuf), _mm_andnot_si128(mask, sel)); - - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - -#ifdef SIMD_LATENCYOPT - return data + 4 + datacnt; -#else - return data + 4 + kDecodeBytesGroupCount[mask0] + kDecodeBytesGroupCount[mask1]; -#endif - } - - case 2: - case 7: - { -#ifdef SIMD_LATENCYOPT - unsigned long long data64; - memcpy(&data64, data, 8); - data64 &= data64 >> 1; - data64 &= data64 >> 2; - - // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used mode 3 - int datacnt = int(((data64 & 0x1111111111111111ull) * 0x1111111111111111ull) >> 60); -#endif - - __m128i sel4 = _mm_loadl_epi64(reinterpret_cast(data)); - __m128i rest = _mm_loadu_si128(reinterpret_cast(data + 8)); - - __m128i sel44 = _mm_unpacklo_epi8(_mm_srli_epi16(sel4, 4), sel4); - __m128i sel = _mm_and_si128(sel44, _mm_set1_epi8(15)); - - __m128i mask = _mm_cmpeq_epi8(sel, _mm_set1_epi8(15)); - int mask16 = _mm_movemask_epi8(mask); - unsigned char mask0 = (unsigned char)(mask16 & 255); - unsigned char mask1 = (unsigned char)(mask16 >> 8); - - __m128i shuf = decodeShuffleMask(mask0, mask1); - __m128i result = _mm_or_si128(_mm_shuffle_epi8(rest, shuf), _mm_andnot_si128(mask, sel)); - - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - -#ifdef SIMD_LATENCYOPT - return data + 8 + datacnt; -#else - return data + 8 + kDecodeBytesGroupCount[mask0] + kDecodeBytesGroupCount[mask1]; -#endif - } - - case 3: - case 8: - { - __m128i result = _mm_loadu_si128(reinterpret_cast(data)); - - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - - return data + 16; - } - - case 5: - { - __m128i rest = _mm_loadu_si128(reinterpret_cast(data + 2)); - - unsigned char mask0 = data[0]; - unsigned char mask1 = data[1]; - - __m128i shuf = decodeShuffleMask(mask0, mask1); - __m128i result = _mm_shuffle_epi8(rest, shuf); - - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - - return data + 2 + kDecodeBytesGroupCount[mask0] + kDecodeBytesGroupCount[mask1]; - } - - default: - SIMD_UNREACHABLE(); // unreachable - } } #endif #ifdef SIMD_AVX -static const __m128i kDecodeBytesGroupConfig[8][2] = { - {_mm_setzero_si128(), _mm_setzero_si128()}, - {_mm_set1_epi8(3), _mm_setr_epi8(6, 4, 2, 0, 14, 12, 10, 8, 22, 20, 18, 16, 30, 28, 26, 24)}, - {_mm_set1_epi8(15), _mm_setr_epi8(4, 0, 12, 8, 20, 16, 28, 24, 36, 32, 44, 40, 52, 48, 60, 56)}, - {_mm_setzero_si128(), _mm_setzero_si128()}, - {_mm_setzero_si128(), _mm_setzero_si128()}, - {_mm_set1_epi8(1), _mm_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)}, - {_mm_set1_epi8(3), _mm_setr_epi8(6, 4, 2, 0, 14, 12, 10, 8, 22, 20, 18, 16, 30, 28, 26, 24)}, - {_mm_set1_epi8(15), _mm_setr_epi8(4, 0, 12, 8, 20, 16, 28, 24, 36, 32, 44, 40, 52, 48, 60, 56)}, +// sent mask, multishift control +static const unsigned char kDecodeBytesGroupConfig[9][2][16] = { + {{0}, {0}}, + {{3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}, {6, 4, 2, 0, 14, 12, 10, 8, 22, 20, 18, 16, 30, 28, 26, 24}}, + {{15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15}, {4, 0, 12, 8, 20, 16, 28, 24, 36, 32, 44, 40, 52, 48, 60, 56}}, + {{0}, {0}}, + {{0}, {0}}, + {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}}, + {{3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}, {6, 4, 2, 0, 14, 12, 10, 8, 22, 20, 18, 16, 30, 28, 26, 24}}, + {{15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15}, {4, 0, 12, 8, 20, 16, 28, 24, 36, 32, 44, 40, 52, 48, 60, 56}}, + {{0}, {0}}, }; SIMD_TARGET -inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) +static inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) { - switch (hbits) - { - case 0: - case 4: - { - __m128i result = _mm_setzero_si128(); + // 0 for 1-bit, 1 for 2-bit, 2 for 4-bit, 3 for 8-bit, and 4 for 0-bit as it makes some of the uses easier + static const int hbtn[9] = {4, 1, 2, 3, 4, 0, 1, 2, 3}; - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); + int n = hbtn[hbits]; - return data; - } +#ifdef SIMD_LATENCYOPT + unsigned long long data64; + memcpy(&data64, data, 8); + data64 &= data64 >> n; + data64 &= data64 >> (n >> 1); - case 5: // 1-bit - case 1: // 2-bit - case 6: - case 2: // 4-bit - case 7: - { - const unsigned char* skip = data + (2 << (hbits < 3 ? hbits : hbits - 5)); + // mask out one bit per group that is set if all group bits were 1 + static const unsigned long long lanes[9] = {0, 0x55555555, 0x1111111111111111ull, 0, 0, 0xffff, 0x55555555, 0x1111111111111111ull, 0}; + int datacnt = int(_mm_popcnt_u64(data64 & lanes[hbits])); +#endif - __m128i selb = _mm_loadl_epi64(reinterpret_cast(data)); - __m128i rest = _mm_loadu_si128(reinterpret_cast(skip)); + // for 8-bit groups, instead of loading the bytes through 'data', we load them through 'skip' as they are easier to preserve + // for 0-bit groups, we use a masked load so that we load zero bytes; in both cases the shift wraps to zero + const unsigned char* skip = data + ((2 << n) & 15); - __m128i sent = kDecodeBytesGroupConfig[hbits][0]; - __m128i ctrl = kDecodeBytesGroupConfig[hbits][1]; + __m128i selb = _mm_loadl_epi64(reinterpret_cast(data)); + __m128i rest = _mm_maskz_loadu_epi8(__mmask16((n >> 2) - 1), skip); - __m128i selw = _mm_shuffle_epi32(selb, 0x44); - __m128i sel = _mm_and_si128(sent, _mm_multishift_epi64_epi8(ctrl, selw)); - __mmask16 mask16 = _mm_cmp_epi8_mask(sel, sent, _MM_CMPINT_EQ); + __m128i sent = _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][0])); + __m128i ctrl = _mm_loadu_si128(reinterpret_cast(kDecodeBytesGroupConfig[hbits][1])); - __m128i result = _mm_mask_expand_epi8(sel, mask16, rest); + // unpack 1, 2 or 4-bit values using multishift and mask the result; for 0/8-bit values, sel is always 0 + __m128i selw = _mm_shuffle_epi32(selb, 0x44); + __m128i sel = _mm_and_si128(sent, _mm_multishift_epi64_epi8(ctrl, selw)); - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); + // compare sel to sentinel and combine; mask is 0 for 0/8-bit as sent is 0 + __mmask16 mask16 = _mm_cmp_epi8_mask(sel, sent, _MM_CMPINT_EQ); + __m128i result = _mm_mask_expand_epi8(sel, mask16, rest); - return skip + _mm_popcnt_u32(mask16); - } + _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - case 3: - case 8: - { - __m128i result = _mm_loadu_si128(reinterpret_cast(data)); - - _mm_storeu_si128(reinterpret_cast<__m128i*>(buffer), result); - - return data + 16; - } - - default: - SIMD_UNREACHABLE(); // unreachable - } +#ifdef SIMD_LATENCYOPT + // datacnt is 0 for 8-bit groups so we can't use skip to advance; 0-bit groups wrap the shift to zero + return data + ((2 << n) & 31) + datacnt; +#else + // mask16 is all 1s for 0-bit groups and we need to use zero instead + return skip + _mm_popcnt_u32(n == 4 ? 0 : mask16); +#endif } #endif @@ -1037,7 +981,7 @@ inline void neonMoveMask(uint8x16_t mask, unsigned char& mask0, unsigned char& m } SIMD_TARGET -inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) +static inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) { switch (hbits) { @@ -1060,9 +1004,9 @@ inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsi data32 &= data32 >> 1; // arrange bits such that low bits of nibbles of data64 contain all 2-bit elements of data32 - unsigned long long data64 = ((unsigned long long)data32 << 30) | (data32 & 0x3fffffff); + unsigned long long data64 = ((unsigned long long)data32 << 30) | data32; - // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used mode 3 + // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used a cheaper mode int datacnt = int(((data64 & 0x1111111111111111ull) * 0x1111111111111111ull) >> 60); #endif @@ -1098,7 +1042,7 @@ inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsi data64 &= data64 >> 1; data64 &= data64 >> 2; - // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used mode 3 + // adds all 1-bit nibbles together; the sum fits in 4 bits because datacnt=16 would have used a cheaper mode int datacnt = int(((data64 & 0x1111111111111111ull) * 0x1111111111111111ull) >> 60); #endif @@ -1179,7 +1123,7 @@ inline void wasmMoveMask(v128_t mask, unsigned char& mask0, unsigned char& mask1 } SIMD_TARGET -inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) +static inline const unsigned char* decodeBytesGroupSimd(const unsigned char* data, unsigned char* buffer, int hbits) { switch (hbits) { @@ -1416,7 +1360,7 @@ inline v128_t unzigzag16(v128_t v) SIMD_TARGET inline v128_t rotate32(v128_t v, int r) { - return wasm_v128_or(wasm_i32x4_shl(v, r), wasm_i32x4_shr(v, 32 - r)); + return wasm_v128_or(wasm_i32x4_shl(v, r), wasm_u32x4_shr(v, 32 - r)); } #endif @@ -1443,6 +1387,22 @@ static const unsigned char* decodeBytesSimd(const unsigned char* data, const uns size_t header_offset = i / kByteGroupSize; unsigned char header_byte = header[header_offset / 4]; +#if defined(SIMD_SSE) || defined(SIMD_AVX) + // very-fast-path: for consecutive 4 groups that are all 0-bit (v0/0, v1/0/0000) or 8-bit (v0/3333, v1/1/3333), + // the branchless decoders are slower than branching over the decoding of 4 groups and issuing a few load/store ops + if (hshift != 5 && header_byte == 0) + { + memset(buffer + i, 0, kByteGroupSize * 4); + continue; + } + else if (hshift != 4 && header_byte == 255) + { + memcpy(buffer + i, data, kByteGroupSize * 4); + data += kByteGroupSize * 4; + continue; + } +#endif + data = decodeBytesGroupSimd(data, buffer + i + kByteGroupSize * 0, hshift + ((header_byte >> 0) & 3)); data = decodeBytesGroupSimd(data, buffer + i + kByteGroupSize * 1, hshift + ((header_byte >> 2) & 3)); data = decodeBytesGroupSimd(data, buffer + i + kByteGroupSize * 2, hshift + ((header_byte >> 4) & 3)); @@ -1560,6 +1520,14 @@ static const unsigned char* decodeVertexBlockSimd(const unsigned char* data, con size_t vertex_count_aligned = (vertex_count + kByteGroupSize - 1) & ~(kByteGroupSize - 1); + // we can decode directly into the output buffer if vertex count is aligned to 16 (delta decode works 16 vertices at a time) + // this uses strided writes and also reads the last vertex once, which is bad for performance for write-combined memory so we only enable this if configured +#ifdef MESHOPTIMIZER_VERTEXCODEC_ZEROCOPY + unsigned char* target = vertex_count == vertex_count_aligned ? vertex_data : transposed; +#else + unsigned char* target = transposed; +#endif + size_t control_size = version == 0 ? 0 : vertex_size / 4; if (size_t(data_end - data) < control_size) return NULL; @@ -1605,22 +1573,23 @@ static const unsigned char* decodeVertexBlockSimd(const unsigned char* data, con switch (channel & 3) { case 0: - decodeDeltas4Simd<0>(buffer, transposed + k, vertex_count_aligned, vertex_size, last_vertex + k, 0); + decodeDeltas4Simd<0>(buffer, target + k, vertex_count_aligned, vertex_size, last_vertex + k, 0); break; case 1: - decodeDeltas4Simd<1>(buffer, transposed + k, vertex_count_aligned, vertex_size, last_vertex + k, 0); + decodeDeltas4Simd<1>(buffer, target + k, vertex_count_aligned, vertex_size, last_vertex + k, 0); break; case 2: - decodeDeltas4Simd<2>(buffer, transposed + k, vertex_count_aligned, vertex_size, last_vertex + k, (32 - (channel >> 4)) & 31); + decodeDeltas4Simd<2>(buffer, target + k, vertex_count_aligned, vertex_size, last_vertex + k, (32 - (channel >> 4)) & 31); break; default: return NULL; // invalid channel type } } - memcpy(vertex_data, transposed, vertex_count * vertex_size); + if (target == transposed) + memcpy(vertex_data, transposed, vertex_count * vertex_size); - memcpy(last_vertex, &transposed[vertex_size * (vertex_count - 1)], vertex_size); + memcpy(last_vertex, &target[vertex_size * (vertex_count - 1)], vertex_size); return data; } @@ -1682,7 +1651,7 @@ size_t meshopt_encodeVertexBufferLevel(unsigned char* buffer, size_t buffer_size for (size_t k = 0; k < vertex_size; k += 4) { int rot = level >= 3 ? estimateRotate(vertex_data, vertex_count, vertex_size, k, /* group_size= */ 16) : 0; - int channel = estimateChannel(vertex_data, vertex_count, vertex_size, k, vertex_block_size, /* block_skip= */ 3, /* max_channels= */ level >= 3 ? 3 : 2, rot); + int channel = estimateChannel(vertex_data, vertex_count, vertex_size, k, vertex_block_size, /* block_skip= */ 3, /* max_channel= */ level >= 3 ? 3 : 2, rot); assert(unsigned(channel) < 2 || ((channel & 3) == 2 && unsigned(channel >> 4) < 8)); channels[k / 4] = (unsigned char)channel; @@ -1837,7 +1806,8 @@ int meshopt_decodeVertexBuffer(void* destination, size_t vertex_count, size_t ve const unsigned char* (*decode)(const unsigned char*, const unsigned char*, unsigned char*, size_t, size_t, unsigned char[256], const unsigned char*, int) = NULL; #if defined(SIMD_SSE) && defined(SIMD_FALLBACK) - decode = (cpuid & (1 << 9)) ? decodeVertexBlockSimd : decodeVertexBlock; + const unsigned int cpumask = (1 << 9) | (1 << 23); // SSSE3+POPCNT + decode = (cpuid & cpumask) == cpumask ? decodeVertexBlockSimd : decodeVertexBlock; #elif defined(SIMD_SSE) || defined(SIMD_AVX) || defined(SIMD_NEON) || defined(SIMD_WASM) decode = decodeVertexBlockSimd; #else diff --git a/3rdparty/meshoptimizer/src/vertexfilter.cpp b/3rdparty/meshoptimizer/src/vertexfilter.cpp index 7a7e67a39..a8d89d1aa 100644 --- a/3rdparty/meshoptimizer/src/vertexfilter.cpp +++ b/3rdparty/meshoptimizer/src/vertexfilter.cpp @@ -13,7 +13,7 @@ #endif // MSVC supports compiling SSE2 code regardless of compile options; we assume all 32-bit CPUs support SSE2 -#if !defined(SIMD_SSE) && defined(_MSC_VER) && !defined(__clang__) && (defined(_M_IX86) || defined(_M_X64)) +#if !defined(SIMD_SSE) && defined(_MSC_VER) && !defined(__clang__) && (defined(_M_IX86) || (defined(_M_X64) && !defined(_M_ARM64EC))) #define SIMD_SSE #endif @@ -23,7 +23,7 @@ #endif // On MSVC, we assume that ARM builds always target NEON-capable devices -#if !defined(SIMD_NEON) && defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64)) +#if !defined(SIMD_NEON) && defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64) || defined(_M_ARM64EC)) #define SIMD_NEON #endif @@ -113,7 +113,7 @@ static void decodeFilterQuat(short* data, size_t count) for (size_t i = 0; i < count; ++i) { - // recover scale from the high byte of the component + // recover scale from the high byte/word of the component int sf = data[i * 4 + 3] | 3; float s = float(sf); @@ -222,6 +222,11 @@ static void dispatchSimd(void (*process)(T*, size_t), T* data, size_t count, siz size_t count4 = count & ~size_t(3); process(data, count4); +#ifdef MESHOPTIMIZER_VERTEXFILTER_SIMDNOTAIL + // optionally omit tail processing to improve code size, expecting the caller to pass aligned counts + assert(count4 == count); + (void)stride; +#else if (count4 < count) { T tail[4 * 4] = {}; // max stride 4, max count 4 @@ -232,6 +237,7 @@ static void dispatchSimd(void (*process)(T*, size_t), T* data, size_t count, siz process(tail, count - count4); memcpy(data + count4 * stride, tail, tail_size); } +#endif } inline uint64_t rotateleft64(uint64_t v, int x) @@ -536,7 +542,7 @@ static void decodeFilterColorSimd16(unsigned short* data, size_t count) } #endif -#if defined(SIMD_NEON) && !defined(__aarch64__) && !defined(_M_ARM64) +#if defined(SIMD_NEON) && !defined(__aarch64__) && !(defined(_M_ARM64) || defined(_M_ARM64EC)) inline float32x4_t vsqrtq_f32(float32x4_t x) { float32x4_t r = vrsqrteq_f32(x); @@ -640,7 +646,7 @@ static void decodeFilterOctSimd16(short* data, size_t count) // compute normal length & scale float32x4_t ll = vfmaq_f32(vfmaq_f32(vmulq_f32(x, x), y, y), z, z); -#if !defined(__aarch64__) && !defined(_M_ARM64) +#if !defined(__aarch64__) && !(defined(_M_ARM64) || defined(_M_ARM64EC)) float32x4_t rl = vrsqrteq_f32(ll); rl = vmulq_f32(rl, vrsqrtsq_f32(vmulq_f32(rl, ll), rl)); // refine rsqrt estimate float32x4_t s = vmulq_f32(vdupq_n_f32(32767.f), rl); @@ -888,7 +894,7 @@ static void decodeFilterOctSimd8(signed char* data, size_t count) v128_t z = wasm_f32x4_sub(wasm_f32x4_convert_i32x4(zf), wasm_f32x4_add(wasm_f32x4_abs(x), wasm_f32x4_abs(y))); // fixup octahedral coordinates for z<0 - // note: i32x4_min with 0 is equvalent to f32x4_min + // note: i32x4_min with 0 is equivalent to f32x4_min v128_t t = wasm_i32x4_min(z, wasm_i32x4_splat(0)); x = wasm_f32x4_add(x, wasm_v128_xor(t, wasm_v128_and(x, sign))); @@ -919,8 +925,7 @@ static void decodeFilterOctSimd8(signed char* data, size_t count) static void decodeFilterOctSimd16(short* data, size_t count) { const v128_t sign = wasm_f32x4_splat(-0.f); - // TODO: volatile here works around LLVM mis-optimizing code; https://github.com/llvm/llvm-project/issues/149457 - volatile v128_t zmask = wasm_i32x4_splat(0x7fff); + const v128_t zmask = wasm_i32x4_splat(0x7fff); for (size_t i = 0; i < count; i += 4) { @@ -944,7 +949,7 @@ static void decodeFilterOctSimd16(short* data, size_t count) v128_t z = wasm_f32x4_sub(wasm_f32x4_convert_i32x4(zf), wasm_f32x4_add(wasm_f32x4_abs(x), wasm_f32x4_abs(y))); // fixup octahedral coordinates for z<0 - // note: i32x4_min with 0 is equvalent to f32x4_min + // note: i32x4_min with 0 is equivalent to f32x4_min v128_t t = wasm_i32x4_min(z, wasm_i32x4_splat(0)); x = wasm_f32x4_add(x, wasm_v128_xor(t, wasm_v128_and(x, sign))); @@ -1069,9 +1074,6 @@ static void decodeFilterExpSimd(unsigned int* data, size_t count) static void decodeFilterColorSimd8(unsigned char* data, size_t count) { - // TODO: volatile here works around LLVM mis-optimizing code; https://github.com/llvm/llvm-project/issues/149457 - volatile v128_t zero = wasm_i32x4_splat(0); - for (size_t i = 0; i < count; i += 4) { v128_t c4 = wasm_v128_load(&data[i * 4]); @@ -1080,7 +1082,7 @@ static void decodeFilterColorSimd8(unsigned char* data, size_t count) v128_t yf = wasm_v128_and(c4, wasm_i32x4_splat(0xff)); v128_t cof = wasm_i32x4_shr(wasm_i32x4_shl(c4, 16), 24); v128_t cgf = wasm_i32x4_shr(wasm_i32x4_shl(c4, 8), 24); - v128_t af = wasm_v128_or(zero, wasm_u32x4_shr(c4, 24)); + v128_t af = wasm_u32x4_shr(c4, 24); // recover scale from alpha high bit v128_t as = af; @@ -1120,9 +1122,6 @@ static void decodeFilterColorSimd8(unsigned char* data, size_t count) static void decodeFilterColorSimd16(unsigned short* data, size_t count) { - // TODO: volatile here works around LLVM mis-optimizing code; https://github.com/llvm/llvm-project/issues/149457 - volatile v128_t zero = wasm_i32x4_splat(0); - for (size_t i = 0; i < count; i += 4) { v128_t c4_0 = wasm_v128_load(&data[(i + 0) * 4]); @@ -1136,7 +1135,7 @@ static void decodeFilterColorSimd16(unsigned short* data, size_t count) v128_t yf = wasm_v128_and(c4_yco, wasm_i32x4_splat(0xffff)); v128_t cof = wasm_i32x4_shr(c4_yco, 16); v128_t cgf = wasm_i32x4_shr(wasm_i32x4_shl(c4_cga, 16), 16); - v128_t af = wasm_v128_or(zero, wasm_u32x4_shr(c4_cga, 16)); + v128_t af = wasm_u32x4_shr(c4_cga, 16); // recover scale from alpha high bit v128_t as = af; @@ -1189,8 +1188,9 @@ inline int optlog2(float v) } u; u.f = v; - // +1 accounts for implicit 1. in mantissa; denormalized numbers will end up clamped to min_exp by calling code - return v == 0 ? 0 : int((u.ui >> 23) & 0xff) - 127 + 1; + // +1 accounts for implicit 1. in mantissa; zero and denormalized numbers will end up clamped to min_exp by calling code + // this is important for shared exponent modes, as a zero component should not inflate the shared exponent + return int((u.ui >> 23) & 0xff) - 127 + 1; } // optimized variant of ldexp @@ -1276,7 +1276,8 @@ void meshopt_decodeFilterColor(void* buffer, size_t count, size_t stride) void meshopt_encodeFilterOct(void* destination, size_t count, size_t stride, int bits, const float* data) { assert(stride == 4 || stride == 8); - assert(bits >= 1 && bits <= 16); + assert(bits >= 2 && bits <= 16); + assert(stride == 8 || bits <= 8); signed char* d8 = static_cast(destination); short* d16 = static_cast(destination); @@ -1367,11 +1368,12 @@ void meshopt_encodeFilterExp(void* destination_, size_t count, size_t stride, in const int min_exp = -100; + // used to track maximum exponent for EncodeExpSharedComponent and zero exponent reuse for EncodeExpSeparate + for (size_t j = 0; j < stride_float; ++j) + component_exp[j] = (mode == meshopt_EncodeExpSeparate) ? 0 : min_exp; + if (mode == meshopt_EncodeExpSharedComponent) { - for (size_t j = 0; j < stride_float; ++j) - component_exp[j] = min_exp; - for (size_t i = 0; i < count; ++i) { const float* v = &data[i * stride_float]; @@ -1409,7 +1411,9 @@ void meshopt_encodeFilterExp(void* destination_, size_t count, size_t stride, in { int e = optlog2(v[j]); - component_exp[j] = (min_exp < e) ? e : min_exp; + // zero values inherit the last encoded exponent to keep the encoded data more compressible + if (v[j] != 0) + component_exp[j] = (min_exp < e) ? e : min_exp; } } else if (mode == meshopt_EncodeExpClamped) @@ -1435,9 +1439,14 @@ void meshopt_encodeFilterExp(void* destination_, size_t count, size_t stride, in exp -= (bits - 1); // compute renormalized rounded mantissa for each component - int mmask = (1 << 24) - 1; + const int mmask = (1 << 24) - 1; + const int mmax = mmask >> 1; + int m = int(v[j] * optexp2(-exp) + (v[j] >= 0 ? 0.5f : -0.5f)); + // clamp mantissa so that rounding away from zero doesn't overflow + m = m > mmax ? mmax : m; + d[j] = (m & mmask) | (unsigned(exp) << 24); } } @@ -1447,6 +1456,7 @@ void meshopt_encodeFilterColor(void* destination, size_t count, size_t stride, i { assert(stride == 4 || stride == 8); assert(bits >= 2 && bits <= 16); + assert(stride == 8 || bits <= 8); unsigned char* d8 = static_cast(destination); unsigned short* d16 = static_cast(destination); @@ -1469,7 +1479,7 @@ void meshopt_encodeFilterColor(void* destination, size_t count, size_t stride, i assert(unsigned((fy + fco - fcg) | (fy + fcg) | (fy - fco - fcg)) < (1u << bits)); // alpha: K-1-bit encoding with high bit set to 1 - int fa = meshopt_quantizeUnorm(c[3], bits - 1) | (1 << (bits - 1)); + int fa = (meshopt_quantizeUnorm(c[3], bits) >> 1) | (1 << (bits - 1)); if (stride == 4) { diff --git a/3rdparty/meshoptimizer/src/vfetchoptimizer.cpp b/3rdparty/meshoptimizer/src/vfetchoptimizer.cpp index 465d6df5c..4847465b6 100644 --- a/3rdparty/meshoptimizer/src/vfetchoptimizer.cpp +++ b/3rdparty/meshoptimizer/src/vfetchoptimizer.cpp @@ -6,8 +6,6 @@ size_t meshopt_optimizeVertexFetchRemap(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count) { - assert(index_count % 3 == 0); - memset(destination, -1, vertex_count * sizeof(unsigned int)); unsigned int next_vertex = 0; @@ -30,7 +28,6 @@ size_t meshopt_optimizeVertexFetchRemap(unsigned int* destination, const unsigne size_t meshopt_optimizeVertexFetch(void* destination, unsigned int* indices, size_t index_count, const void* vertices, size_t vertex_count, size_t vertex_size) { - assert(index_count % 3 == 0); assert(vertex_size > 0 && vertex_size <= 256); meshopt_Allocator allocator;