feat: LOD 网格分级——meshopt 解析期简化 + 按距离/级联选级

vendored meshoptimizer 1.2（submodule 环境故障改直接 vendor）。Mesh
多 LOD：各级索引连续拼进同一 index buffer（共享顶点，LodRange 区间
选择），解析期 meshopt_simplify 生成链（LOD0=原始，每级目标减半，
相对误差 1e-2，小网格复用上级）；都在 worker 线程，加载零额外卡顿。

SceneRenderer 按 (mesh, LOD) 聚合 batch（lookup key = 指针<<8|lod）：
主视图按屏占比近似 radius/dist 选级（阈值 10/20 倍半径）；阴影级联
按级联序降级（越远越简）。Draw 按 batch.lod 选索引区间。

A/B 对比（100 背包 10×10，RTX 5060 Laptop，vsync off）：
- 纯背包：LOD 开 93-95 fps / 关 90-107 fps —— 无显著差异
- 加 SunTemple：LOD 开 52-56 / 关 51-55 —— 无显著差异
结论：当前场景瓶颈不在顶点数（5060 顶点吞吐富余），LOD 基础设施
保留（几何复杂资产/低端 GPU/更大场景生效），帧率问题转向 Tracy
GPU 时间线定位真实瓶颈。

坑：LodRange 记录"LOD 不足复用最后一级"时 first 不得继续累加
（buffer 无重复数据，累加必越界——VUID-08798 抓的）。

Lib/meshoptimizer/src/allocator.cpp

@@ -0,0 +1,21 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+
+#ifdef MESHOPTIMIZER_ALLOC_EXPORT
+meshopt_Allocator::Storage& meshopt_Allocator::storage()
+{
+	static Storage s = {::operator new, ::operator delete };
+	return s;
+}
+#endif
+
+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;
+}

Lib/meshoptimizer/src/indexanalyzer.cpp

@@ -0,0 +1,125 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <string.h>
+
+meshopt_VertexCacheStatistics meshopt_analyzeVertexCache(const unsigned int* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int warp_size, unsigned int primgroup_size)
+{
+	assert(index_count % 3 == 0);
+	assert(cache_size >= 3);
+	assert(warp_size == 0 || warp_size >= 3);
+
+	meshopt_Allocator allocator;
+
+	meshopt_VertexCacheStatistics result = {};
+
+	unsigned int warp_offset = 0;
+	unsigned int primgroup_offset = 0;
+
+	unsigned int* cache_timestamps = allocator.allocate<unsigned int>(vertex_count);
+	memset(cache_timestamps, 0, vertex_count * sizeof(unsigned int));
+
+	unsigned int timestamp = cache_size + 1;
+
+	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);
+
+		bool ac = (timestamp - cache_timestamps[a]) > cache_size;
+		bool bc = (timestamp - cache_timestamps[b]) > cache_size;
+		bool cc = (timestamp - cache_timestamps[c]) > cache_size;
+
+		// flush cache if triangle doesn't fit into warp or into the primitive buffer
+		if ((primgroup_size && primgroup_offset == primgroup_size) || (warp_size && warp_offset + ac + bc + cc > warp_size))
+		{
+			result.warps_executed += warp_offset > 0;
+
+			warp_offset = 0;
+			primgroup_offset = 0;
+
+			// reset cache
+			timestamp += cache_size + 1;
+		}
+
+		// update cache and add vertices to warp
+		for (int j = 0; j < 3; ++j)
+		{
+			unsigned int index = indices[i + j];
+
+			if (timestamp - cache_timestamps[index] > cache_size)
+			{
+				cache_timestamps[index] = timestamp++;
+				result.vertices_transformed++;
+				warp_offset++;
+			}
+		}
+
+		primgroup_offset++;
+	}
+
+	size_t unique_vertex_count = 0;
+
+	for (size_t i = 0; i < vertex_count; ++i)
+		unique_vertex_count += cache_timestamps[i] > 0;
+
+	result.warps_executed += warp_offset > 0;
+
+	result.acmr = index_count == 0 ? 0 : float(result.vertices_transformed) / float(index_count / 3);
+	result.atvr = unique_vertex_count == 0 ? 0 : float(result.vertices_transformed) / float(unique_vertex_count);
+
+	return result;
+}
+
+meshopt_VertexFetchStatistics meshopt_analyzeVertexFetch(const unsigned int* indices, size_t index_count, size_t vertex_count, size_t vertex_size)
+{
+	assert(vertex_size > 0 && vertex_size <= 256);
+
+	meshopt_Allocator allocator;
+
+	meshopt_VertexFetchStatistics result = {};
+
+	unsigned char* vertex_visited = allocator.allocate<unsigned char>(vertex_count);
+	memset(vertex_visited, 0, vertex_count);
+
+	const size_t kCacheLine = 64;
+	const size_t kCacheSize = 128 * 1024;
+
+	// simple direct mapped cache; on typical mesh data this is close to 4-way cache, and this model is a gross approximation anyway
+	size_t cache[kCacheSize / kCacheLine] = {};
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		vertex_visited[index] = 1;
+
+		size_t start_address = index * vertex_size;
+		size_t end_address = start_address + vertex_size;
+
+		size_t start_tag = start_address / kCacheLine;
+		size_t end_tag = (end_address + kCacheLine - 1) / kCacheLine;
+
+		assert(start_tag < end_tag);
+
+		for (size_t tag = start_tag; tag < end_tag; ++tag)
+		{
+			size_t line = tag % (sizeof(cache) / sizeof(cache[0]));
+
+			// we store +1 since cache is filled with 0 by default
+			result.bytes_fetched += (cache[line] != tag + 1) * kCacheLine;
+			cache[line] = tag + 1;
+		}
+	}
+
+	size_t unique_vertex_count = 0;
+
+	for (size_t i = 0; i < vertex_count; ++i)
+		unique_vertex_count += vertex_visited[i];
+
+	result.overfetch = unique_vertex_count == 0 ? 0 : float(result.bytes_fetched) / float(unique_vertex_count * vertex_size);
+
+	return result;
+}

Lib/meshoptimizer/src/indexcodec.cpp

@@ -0,0 +1,703 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <string.h>
+
+// This work is based on:
+// Fabian Giesen. Simple lossless index buffer compression & follow-up. 2013
+// Conor Stokes. Vertex Cache Optimised Index Buffer Compression. 2014
+namespace meshopt
+{
+
+const unsigned char kIndexHeader = 0xe0;
+const unsigned char kSequenceHeader = 0xd0;
+
+static int gEncodeIndexVersion = 1;
+const int kDecodeIndexVersion = 1;
+
+typedef unsigned int VertexFifo[16];
+typedef unsigned int EdgeFifo[16][2];
+
+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
+};
+
+static int rotateTriangle(unsigned int a, unsigned int b, unsigned int c, unsigned int next)
+{
+	(void)a;
+
+	return (b == next) ? 1 : (c == next ? 2 : 0);
+}
+
+static int getEdgeFifo(EdgeFifo fifo, unsigned int a, unsigned int b, unsigned int c, size_t offset)
+{
+	for (int i = 0; i < 16; ++i)
+	{
+		size_t index = (offset - 1 - i) & 15;
+
+		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 pushEdgeFifo(EdgeFifo fifo, unsigned int a, unsigned int b, size_t& offset)
+{
+	fifo[offset][0] = a;
+	fifo[offset][1] = b;
+	offset = (offset + 1) & 15;
+}
+
+static int getVertexFifo(VertexFifo fifo, unsigned int v, size_t offset)
+{
+	for (int i = 0; i < 16; ++i)
+	{
+		size_t index = (offset - 1 - i) & 15;
+
+		if (fifo[index] == v)
+			return i;
+	}
+
+	return -1;
+}
+
+static void pushVertexFifo(VertexFifo fifo, unsigned int v, size_t& offset, int cond = 1)
+{
+	fifo[offset] = v;
+	offset = (offset + cond) & 15;
+}
+
+static void encodeVByte(unsigned char*& data, unsigned int v)
+{
+	// encode 32-bit value in up to 5 7-bit groups
+	do
+	{
+		*data++ = (v & 127) | (v > 127 ? 128 : 0);
+		v >>= 7;
+	} while (v);
+}
+
+static unsigned int decodeVByte(const unsigned char*& data)
+{
+	unsigned char lead = *data++;
+
+	// fast path: single byte
+	if (lead < 128)
+		return lead;
+
+	// slow path: up to 4 extra bytes
+	// note that this loop always terminates, which is important for malformed data
+	unsigned int result = lead & 127;
+	unsigned int shift = 7;
+
+	for (int i = 0; i < 4; ++i)
+	{
+		unsigned char group = *data++;
+		result |= unsigned(group & 127) << shift;
+		shift += 7;
+
+		if (group < 128)
+			break;
+	}
+
+	return result;
+}
+
+static void encodeIndex(unsigned char*& data, unsigned int index, unsigned int last)
+{
+	unsigned int d = index - last;
+	unsigned int v = (d << 1) ^ (int(d) >> 31);
+
+	encodeVByte(data, v);
+}
+
+static unsigned int decodeIndex(const unsigned char*& data, unsigned int last)
+{
+	unsigned int v = decodeVByte(data);
+	unsigned int d = (v >> 1) ^ -int(v & 1);
+
+	return last + d;
+}
+
+static int getCodeAuxIndex(unsigned char v, const unsigned char* table)
+{
+	for (int i = 0; i < 16; ++i)
+		if (table[i] == v)
+			return i;
+
+	return -1;
+}
+
+static void* writeTriangle(void* destination, size_t index_size, unsigned int a, unsigned int b, unsigned int c)
+{
+	if (index_size == 2)
+	{
+		unsigned short* tri = static_cast<unsigned short*>(destination);
+		tri[0] = (unsigned short)(a);
+		tri[1] = (unsigned short)(b);
+		tri[2] = (unsigned short)(c);
+
+		return tri + 3;
+	}
+	else
+	{
+		unsigned int* tri = static_cast<unsigned int*>(destination);
+		tri[0] = a;
+		tri[1] = b;
+		tri[2] = c;
+
+		return tri + 3;
+	}
+}
+
+} // namespace meshopt
+
+size_t meshopt_encodeIndexBuffer(unsigned char* buffer, size_t buffer_size, const unsigned int* indices, size_t index_count)
+{
+	using namespace meshopt;
+
+	assert(index_count % 3 == 0);
+
+	// the minimum valid encoding is header, 1 byte per triangle and a 16-byte codeaux table
+	if (buffer_size < 1 + index_count / 3 + 16)
+		return 0;
+
+	int version = gEncodeIndexVersion;
+
+	buffer[0] = (unsigned char)(kIndexHeader | version);
+
+	EdgeFifo edgefifo;
+	memset(edgefifo, -1, sizeof(edgefifo));
+
+	VertexFifo vertexfifo;
+	memset(vertexfifo, -1, sizeof(vertexfifo));
+
+	size_t edgefifooffset = 0;
+	size_t vertexfifooffset = 0;
+
+	unsigned int next = 0;
+	unsigned int last = 0;
+
+	unsigned char* code = buffer + 1;
+	unsigned char* data = code + index_count / 3;
+	unsigned char* data_safe_end = buffer + buffer_size - 16;
+
+	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;
+
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		// make sure we have enough space to write a triangle
+		// each triangle writes at most 16 bytes: 1b for codeaux and 5b for each free index
+		// after this we can be sure we can write without extra bounds checks
+		if (data > data_safe_end)
+			return 0;
+
+		int fer = getEdgeFifo(edgefifo, indices[i + 0], indices[i + 1], indices[i + 2], edgefifooffset);
+
+		if (fer >= 0 && (fer >> 2) < 15)
+		{
+			// note: getEdgeFifo implicitly rotates triangles by matching a/b to existing edge
+			const int* order = rotations + (fer & 3);
+
+			unsigned int a = indices[i + order[0]], b = indices[i + order[1]], c = indices[i + order[2]];
+
+			// encode edge index and vertex fifo index, next or free index
+			int fe = fer >> 2;
+			int fc = getVertexFifo(vertexfifo, c, vertexfifooffset);
+
+			int fec = (fc >= 1 && fc < fecmax) ? fc : (c == next ? (next++, 0) : 15);
+
+			if (fec == 15 && version >= 1)
+			{
+				// encode last-1 and last+1 to optimize strip-like sequences
+				if (c + 1 == last)
+					fec = 13, last = c;
+				if (c == last + 1)
+					fec = 14, last = c;
+			}
+
+			*code++ = (unsigned char)((fe << 4) | fec);
+
+			// note that we need to update the last index since free indices are delta-encoded
+			if (fec == 15)
+				encodeIndex(data, c, last), last = c;
+
+			// we only need to push third vertex since first two are likely already in the vertex fifo
+			if (fec == 0 || fec >= fecmax)
+				pushVertexFifo(vertexfifo, c, vertexfifooffset);
+
+			// we only need to push two new edges to edge fifo since the third one is already there
+			pushEdgeFifo(edgefifo, c, b, edgefifooffset);
+			pushEdgeFifo(edgefifo, a, c, edgefifooffset);
+		}
+		else
+		{
+			int rotation = rotateTriangle(indices[i + 0], indices[i + 1], indices[i + 2], next);
+			const int* order = rotations + rotation;
+
+			unsigned int a = indices[i + order[0]], b = indices[i + order[1]], c = indices[i + order[2]];
+
+			// if a/b/c are 0/1/2, we emit a reset code
+			bool reset = false;
+
+			if (a == 0 && b == 1 && c == 2 && next > 0 && version >= 1)
+			{
+				reset = true;
+				next = 0;
+
+				// reset vertex fifo to make sure we don't accidentally reference vertices from that in the future
+				// this makes sure next continues to get incremented instead of being stuck
+				memset(vertexfifo, -1, sizeof(vertexfifo));
+			}
+
+			int fb = getVertexFifo(vertexfifo, b, vertexfifooffset);
+			int fc = getVertexFifo(vertexfifo, c, vertexfifooffset);
+
+			// after rotation, a is almost always equal to next, so we don't waste bits on FIFO encoding for a
+			// note: decoder implicitly assumes that if feb=fec=0, then fea=0 (reset code); this is enforced by rotation
+			int fea = (a == next) ? (next++, 0) : 15;
+			int feb = (fb >= 0 && fb < 14) ? fb + 1 : (b == next ? (next++, 0) : 15);
+			int fec = (fc >= 0 && fc < 14) ? fc + 1 : (c == next ? (next++, 0) : 15);
+
+			// we encode feb & fec in 4 bits using a table if possible, and as a full byte otherwise
+			unsigned char codeaux = (unsigned char)((feb << 4) | fec);
+			int codeauxindex = getCodeAuxIndex(codeaux, codeaux_table);
+
+			// <14 encodes an index into codeaux table, 14 encodes fea=0, 15 encodes fea=15
+			if (fea == 0 && codeauxindex >= 0 && codeauxindex < 14 && !reset)
+			{
+				*code++ = (unsigned char)((15 << 4) | codeauxindex);
+			}
+			else
+			{
+				*code++ = (unsigned char)((15 << 4) | 14 | fea);
+				*data++ = codeaux;
+			}
+
+			// note that we need to update the last index since free indices are delta-encoded
+			if (fea == 15)
+				encodeIndex(data, a, last), last = a;
+
+			if (feb == 15)
+				encodeIndex(data, b, last), last = b;
+
+			if (fec == 15)
+				encodeIndex(data, c, last), last = c;
+
+			// only push vertices that weren't already in fifo
+			if (fea == 0 || fea == 15)
+				pushVertexFifo(vertexfifo, a, vertexfifooffset);
+
+			if (feb == 0 || feb == 15)
+				pushVertexFifo(vertexfifo, b, vertexfifooffset);
+
+			if (fec == 0 || fec == 15)
+				pushVertexFifo(vertexfifo, c, vertexfifooffset);
+
+			// all three edges aren't in the fifo; pushing all of them is important so that we can match them for later triangles
+			pushEdgeFifo(edgefifo, b, a, edgefifooffset);
+			pushEdgeFifo(edgefifo, c, b, edgefifooffset);
+			pushEdgeFifo(edgefifo, a, c, edgefifooffset);
+		}
+	}
+
+	// make sure we have enough space to write codeaux table
+	if (data > data_safe_end)
+		return 0;
+
+	// add codeaux encoding table to the end of the stream; this is used for decoding codeaux *and* as padding
+	// we need padding for decoding to be able to assume that each triangle is encoded as <= 16 bytes of extra data
+	// this is enough space for aux byte + 5 bytes per varint index which is the absolute worst case for any input
+	for (size_t i = 0; i < 16; ++i)
+	{
+		// decoder assumes that table entries never refer to separately encoded indices
+		assert((codeaux_table[i] & 0xf) != 0xf && (codeaux_table[i] >> 4) != 0xf);
+
+		*data++ = codeaux_table[i];
+	}
+
+	// since we encode restarts as codeaux without a table reference, we need to make sure 00 is encoded as a table reference
+	assert(codeaux_table[0] == 0);
+
+	assert(data >= buffer + index_count / 3 + 16);
+	assert(data <= buffer + buffer_size);
+
+	return data - buffer;
+}
+
+size_t meshopt_encodeIndexBufferBound(size_t index_count, size_t vertex_count)
+{
+	assert(index_count % 3 == 0);
+
+	// compute number of bits required for each index
+	unsigned int vertex_bits = 1;
+
+	while (vertex_bits < 32 && vertex_count > size_t(1) << vertex_bits)
+		vertex_bits++;
+
+	// worst-case encoding is 2 header bytes + 3 varint-7 encoded index deltas
+	unsigned int vertex_groups = (vertex_bits + 1 + 6) / 7;
+
+	return 1 + (index_count / 3) * (2 + 3 * vertex_groups) + 16;
+}
+
+void meshopt_encodeIndexVersion(int version)
+{
+	assert(unsigned(version) <= unsigned(meshopt::kDecodeIndexVersion));
+
+	meshopt::gEncodeIndexVersion = version;
+}
+
+int meshopt_decodeIndexVersion(const unsigned char* buffer, size_t buffer_size)
+{
+	if (buffer_size < 1)
+		return -1;
+
+	unsigned char header = buffer[0];
+
+	if ((header & 0xf0) != meshopt::kIndexHeader && (header & 0xf0) != meshopt::kSequenceHeader)
+		return -1;
+
+	int version = header & 0x0f;
+	if (version > meshopt::kDecodeIndexVersion)
+		return -1;
+
+	return version;
+}
+
+int meshopt_decodeIndexBuffer(void* destination, size_t index_count, size_t index_size, const unsigned char* buffer, size_t buffer_size)
+{
+	using namespace meshopt;
+
+	assert(index_count % 3 == 0);
+	assert(index_size == 2 || index_size == 4);
+
+	// the minimum valid encoding is header, 1 byte per triangle and a 16-byte codeaux table
+	if (buffer_size < 1 + index_count / 3 + 16)
+		return -2;
+
+	if ((buffer[0] & 0xf0) != kIndexHeader)
+		return -1;
+
+	int version = buffer[0] & 0x0f;
+	if (version > kDecodeIndexVersion)
+		return -1;
+
+	EdgeFifo edgefifo;
+	memset(edgefifo, -1, sizeof(edgefifo));
+
+	VertexFifo vertexfifo;
+	memset(vertexfifo, -1, sizeof(vertexfifo));
+
+	size_t edgefifooffset = 0;
+	size_t vertexfifooffset = 0;
+
+	unsigned int next = 0;
+	unsigned int last = 0;
+
+	int fecmax = version >= 1 ? 13 : 15;
+
+	// 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* 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;
+
+	while (code < code_end)
+	{
+		unsigned char codetri = *code++;
+
+		if (codetri < 0xf0)
+		{
+			int fe = codetri >> 4;
+
+			// fifo reads are wrapped around 16 entry buffer
+			unsigned int a = edgefifo[(edgefifooffset - 1 - fe) & 15][0];
+			unsigned int b = edgefifo[(edgefifooffset - 1 - fe) & 15][1];
+			unsigned int c = 0;
+
+			int fec = codetri & 15;
+
+			// note: this is the most common path in the entire decoder
+			// inside this if we try to stay branchless (by using cmov/etc.) since these aren't predictable
+			if (fec < fecmax)
+			{
+				// 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;
+				next += fec0;
+
+				// push vertex fifo must match the encoding step *exactly* otherwise the data will not be decoded correctly
+				pushVertexFifo(vertexfifo, c, vertexfifooffset, fec0);
+			}
+			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;
+
+				// 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 * 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);
+			}
+
+			// push edge fifo must match the encoding step *exactly* otherwise the data will not be decoded correctly
+			pushEdgeFifo(edgefifo, c, b, edgefifooffset);
+			pushEdgeFifo(edgefifo, a, c, edgefifooffset);
+
+			// output triangle
+			destination = writeTriangle(destination, index_size, a, b, c);
+		}
+		else
+		{
+			// fast path: read codeaux from the table
+			if (codetri < 0xfe)
+			{
+				unsigned char codeaux = codeaux_table[codetri & 15];
+
+				// note: table can't contain feb/fec=15
+				int feb = codeaux >> 4;
+				int fec = codeaux & 15;
+
+				// fifo reads are wrapped around 16 entry buffer
+				// also note that we increment next for all three vertices before decoding indices - this matches encoder behavior
+				unsigned int a = next++;
+
+				unsigned int bf = vertexfifo[(vertexfifooffset - feb) & 15];
+				unsigned int b = (feb == 0) ? next : bf;
+
+				int feb0 = feb == 0;
+				next += feb0;
+
+				unsigned int cf = vertexfifo[(vertexfifooffset - fec) & 15];
+				unsigned int c = (fec == 0) ? next : cf;
+
+				int fec0 = fec == 0;
+				next += fec0;
+
+				// output triangle
+				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);
+				pushVertexFifo(vertexfifo, b, vertexfifooffset, feb0);
+				pushVertexFifo(vertexfifo, c, vertexfifooffset, fec0);
+
+				pushEdgeFifo(edgefifo, b, a, edgefifooffset);
+				pushEdgeFifo(edgefifo, c, b, edgefifooffset);
+				pushEdgeFifo(edgefifo, a, c, edgefifooffset);
+			}
+			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++;
+
+				int fea = codetri == 0xfe ? 0 : 15;
+				int feb = codeaux >> 4;
+				int fec = codeaux & 15;
+
+				// reset: codeaux is 0 but encoded as not-a-table
+				if (codeaux == 0)
+					next = 0;
+
+				// fifo reads are wrapped around 16 entry buffer
+				// also note that we increment next for all three vertices before decoding indices - this matches encoder behavior
+				unsigned int a = (fea == 0) ? next++ : 0;
+				unsigned int b = (feb == 0) ? next++ : vertexfifo[(vertexfifooffset - feb) & 15];
+				unsigned int c = (fec == 0) ? next++ : vertexfifo[(vertexfifooffset - fec) & 15];
+
+				// note that we need to update the last index since free indices are delta-encoded
+				if (fea == 15)
+					last = a = decodeIndex(data, last);
+
+				if (feb == 15)
+					last = b = decodeIndex(data, last);
+
+				if (fec == 15)
+					last = c = decodeIndex(data, last);
+
+				// output triangle
+				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);
+				pushVertexFifo(vertexfifo, b, vertexfifooffset, (feb == 0) | (feb == 15));
+				pushVertexFifo(vertexfifo, c, vertexfifooffset, (fec == 0) | (fec == 15));
+
+				pushEdgeFifo(edgefifo, b, a, edgefifooffset);
+				pushEdgeFifo(edgefifo, c, b, edgefifooffset);
+				pushEdgeFifo(edgefifo, a, c, edgefifooffset);
+			}
+		}
+	}
+
+	// we should've read all data bytes and stopped at the boundary between data and codeaux table
+	if (data != data_safe_end)
+		return -3;
+
+	return 0;
+}
+
+size_t meshopt_encodeIndexSequence(unsigned char* buffer, size_t buffer_size, const unsigned int* indices, size_t index_count)
+{
+	using namespace meshopt;
+
+	// the minimum valid encoding is header, 1 byte per index and a 4-byte tail
+	if (buffer_size < 1 + index_count + 4)
+		return 0;
+
+	int version = gEncodeIndexVersion;
+
+	buffer[0] = (unsigned char)(kSequenceHeader | version);
+
+	unsigned int last[2] = {};
+	unsigned int current = 0;
+
+	unsigned char* data = buffer + 1;
+	unsigned char* data_safe_end = buffer + buffer_size - 4;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		// make sure we have enough data to write
+		// each index writes at most 5 bytes of data; there's a 4 byte tail after data_safe_end
+		// after this we can be sure we can write without extra bounds checks
+		if (data >= data_safe_end)
+			return 0;
+
+		unsigned int index = indices[i];
+
+		// this is a heuristic that switches between baselines when the delta grows too large
+		// we want the encoded delta to fit into one byte (7 bits), but 2 bits are used for sign and baseline index
+		// for now we immediately switch the baseline when delta grows too large - this can be adjusted arbitrarily
+		int cd = int(index - last[current]);
+		current ^= ((cd < 0 ? -cd : cd) >= 30);
+
+		// encode delta from the last index
+		unsigned int d = index - last[current];
+		unsigned int v = (d << 1) ^ (int(d) >> 31);
+
+		// note: low bit encodes the index of the last baseline which will be used for reconstruction
+		encodeVByte(data, (v << 1) | current);
+
+		// update last for the next iteration that uses it
+		last[current] = index;
+	}
+
+	// make sure we have enough space to write tail
+	if (data > data_safe_end)
+		return 0;
+
+	for (int k = 0; k < 4; ++k)
+		*data++ = 0;
+
+	return data - buffer;
+}
+
+size_t meshopt_encodeIndexSequenceBound(size_t index_count, size_t vertex_count)
+{
+	// compute number of bits required for each index
+	unsigned int vertex_bits = 1;
+
+	while (vertex_bits < 32 && vertex_count > size_t(1) << vertex_bits)
+		vertex_bits++;
+
+	// worst-case encoding is 1 varint-7 encoded index delta for a K bit value and an extra bit
+	unsigned int vertex_groups = (vertex_bits + 1 + 1 + 6) / 7;
+
+	return 1 + index_count * vertex_groups + 4;
+}
+
+int meshopt_decodeIndexSequence(void* destination, size_t index_count, size_t index_size, const unsigned char* buffer, size_t buffer_size)
+{
+	using namespace meshopt;
+
+	// the minimum valid encoding is header, 1 byte per index and a 4-byte tail
+	if (buffer_size < 1 + index_count + 4)
+		return -2;
+
+	if ((buffer[0] & 0xf0) != kSequenceHeader)
+		return -1;
+
+	int version = buffer[0] & 0x0f;
+	if (version > kDecodeIndexVersion)
+		return -1;
+
+	const unsigned char* data = buffer + 1;
+	const unsigned char* data_safe_end = buffer + buffer_size - 4;
+
+	unsigned int last[2] = {};
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		// make sure we have enough data to read
+		// each index reads at most 5 bytes of data; there's a 4 byte tail after data_safe_end
+		// after this we can be sure we can read without extra bounds checks
+		if (data >= data_safe_end)
+			return -2;
+
+		unsigned int v = decodeVByte(data);
+
+		// decode the index of the last baseline
+		unsigned int current = v & 1;
+		v >>= 1;
+
+		// reconstruct index as a delta
+		unsigned int d = (v >> 1) ^ -int(v & 1);
+		unsigned int index = last[current] + d;
+
+		// update last for the next iteration that uses it
+		last[current] = index;
+
+		if (index_size == 2)
+		{
+			static_cast<unsigned short*>(destination)[i] = (unsigned short)(index);
+		}
+		else
+		{
+			static_cast<unsigned int*>(destination)[i] = index;
+		}
+	}
+
+	// we should've read all data bytes and stopped at the boundary between data and tail
+	if (data != data_safe_end)
+		return -3;
+
+	return 0;
+}

Lib/meshoptimizer/src/indexgenerator.cpp

@@ -0,0 +1,825 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <string.h>
+
+// This work is based on:
+// Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003
+// John McDonald, Mark Kilgard. Crack-Free Point-Normal Triangles using Adjacent Edge Normals. 2010
+// John Hable. Variable Rate Shading with Visibility Buffer Rendering. 2024
+namespace meshopt
+{
+
+static unsigned int hashUpdate(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;
+	}
+
+	if (len > 0)
+	{
+		h ^= key[0] | (key[len / 2] << 8) | (key[len - 1] << 16);
+		h *= m;
+		h ^= h >> r;
+	}
+
+	return h;
+}
+
+struct VertexHasher
+{
+	const unsigned char* vertices;
+	size_t vertex_size;
+	size_t vertex_stride;
+
+	size_t hash(unsigned int index) const
+	{
+		return hashUpdate(0, vertices + index * vertex_stride, vertex_size);
+	}
+
+	bool equal(unsigned int lhs, unsigned int rhs) const
+	{
+		return memcmp(vertices + lhs * vertex_stride, vertices + rhs * vertex_stride, vertex_size) == 0;
+	}
+};
+
+struct VertexStreamHasher
+{
+	const meshopt_Stream* streams;
+	size_t stream_count;
+
+	size_t hash(unsigned int index) const
+	{
+		unsigned int h = 0;
+
+		for (size_t i = 0; i < stream_count; ++i)
+		{
+			const meshopt_Stream& s = streams[i];
+			const unsigned char* data = static_cast<const unsigned char*>(s.data);
+
+			h = hashUpdate(h, data + index * s.stride, s.size);
+		}
+
+		return h;
+	}
+
+	bool equal(unsigned int lhs, unsigned int rhs) const
+	{
+		for (size_t i = 0; i < stream_count; ++i)
+		{
+			const meshopt_Stream& s = streams[i];
+			const unsigned char* data = static_cast<const unsigned char*>(s.data);
+
+			if (memcmp(data + lhs * s.stride, data + rhs * s.stride, s.size) != 0)
+				return false;
+		}
+
+		return true;
+	}
+};
+
+struct VertexCustomHasher
+{
+	const float* vertex_positions;
+	size_t vertex_stride_float;
+
+	int (*callback)(void*, unsigned int, unsigned int);
+	void* context;
+
+	size_t hash(unsigned int index) const
+	{
+		const unsigned int* key = reinterpret_cast<const unsigned int*>(vertex_positions + index * vertex_stride_float);
+
+		unsigned int x = key[0], y = key[1], z = key[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;
+
+		// Optimized Spatial Hashing for Collision Detection of Deformable Objects
+		return (x * 73856093) ^ (y * 19349663) ^ (z * 83492791);
+	}
+
+	bool equal(unsigned int lhs, unsigned int rhs) const
+	{
+		const float* lp = vertex_positions + lhs * vertex_stride_float;
+		const float* rp = vertex_positions + rhs * vertex_stride_float;
+
+		if (lp[0] != rp[0] || lp[1] != rp[1] || lp[2] != rp[2])
+			return false;
+
+		return callback ? callback(context, lhs, rhs) : true;
+	}
+};
+
+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;
+
+	size_t hash(unsigned long long edge) const
+	{
+		unsigned int e0 = unsigned(edge >> 32);
+		unsigned int e1 = unsigned(edge);
+
+		unsigned int h1 = remap[e0];
+		unsigned int h2 = remap[e1];
+
+		const unsigned int m = 0x5bd1e995;
+
+		// MurmurHash64B finalizer
+		h1 ^= h2 >> 18;
+		h1 *= m;
+		h2 ^= h1 >> 22;
+		h2 *= m;
+		h1 ^= h2 >> 17;
+		h1 *= m;
+		h2 ^= h1 >> 19;
+		h2 *= m;
+
+		return h2;
+	}
+
+	bool equal(unsigned long long lhs, unsigned long long rhs) const
+	{
+		unsigned int l0 = unsigned(lhs >> 32);
+		unsigned int l1 = unsigned(lhs);
+
+		unsigned int r0 = unsigned(rhs >> 32);
+		unsigned int r1 = unsigned(rhs);
+
+		return remap[l0] == remap[r0] && remap[l1] == remap[r1];
+	}
+};
+
+static size_t hashBuckets(size_t count)
+{
+	size_t buckets = 1;
+	while (buckets < count + count / 4)
+		buckets *= 2;
+
+	return buckets;
+}
+
+template <typename T, typename Hash>
+static T* hashLookup(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;
+}
+
+static void buildPositionRemap(unsigned int* remap, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, meshopt_Allocator& allocator)
+{
+	VertexHasher vertex_hasher = {reinterpret_cast<const unsigned char*>(vertex_positions), 3 * sizeof(float), vertex_positions_stride};
+
+	size_t vertex_table_size = hashBuckets(vertex_count);
+	unsigned int* vertex_table = allocator.allocate<unsigned int>(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 = hashLookup(vertex_table, vertex_table_size, vertex_hasher, index, ~0u);
+
+		if (*entry == ~0u)
+			*entry = index;
+
+		remap[index] = *entry;
+	}
+
+	allocator.deallocate(vertex_table);
+}
+
+template <typename Hash>
+static size_t generateVertexRemap(unsigned int* remap, const unsigned int* indices, size_t index_count, size_t vertex_count, const Hash& hash, meshopt_Allocator& allocator)
+{
+	memset(remap, -1, vertex_count * sizeof(unsigned int));
+
+	size_t table_size = hashBuckets(vertex_count);
+	unsigned int* table = allocator.allocate<unsigned int>(table_size);
+	memset(table, -1, table_size * sizeof(unsigned int));
+
+	unsigned int next_vertex = 0;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices ? indices[i] : unsigned(i);
+		assert(index < vertex_count);
+
+		if (remap[index] != ~0u)
+			continue;
+
+		unsigned int* entry = hashLookup(table, table_size, hash, index, ~0u);
+
+		if (*entry == ~0u)
+		{
+			*entry = index;
+			remap[index] = next_vertex++;
+		}
+		else
+		{
+			assert(remap[*entry] != ~0u);
+			remap[index] = remap[*entry];
+		}
+	}
+
+	assert(next_vertex <= vertex_count);
+	return next_vertex;
+}
+
+template <size_t BlockSize>
+static void remapVertices(void* destination, const void* vertices, size_t vertex_count, size_t vertex_size, const unsigned int* remap)
+{
+	size_t block_size = BlockSize == 0 ? vertex_size : BlockSize;
+	assert(block_size == vertex_size);
+
+	for (size_t i = 0; i < vertex_count; ++i)
+		if (remap[i] != ~0u)
+		{
+			assert(remap[i] < vertex_count);
+			memcpy(static_cast<unsigned char*>(destination) + remap[i] * block_size, static_cast<const unsigned char*>(vertices) + i * block_size, block_size);
+		}
+}
+
+template <typename Hash>
+static void generateShadowBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, const Hash& hash, meshopt_Allocator& allocator)
+{
+	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
+	memset(remap, -1, vertex_count * sizeof(unsigned int));
+
+	size_t table_size = hashBuckets(vertex_count);
+	unsigned int* table = allocator.allocate<unsigned int>(table_size);
+	memset(table, -1, table_size * sizeof(unsigned int));
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		if (remap[index] == ~0u)
+		{
+			unsigned int* entry = hashLookup(table, table_size, hash, index, ~0u);
+
+			if (*entry == ~0u)
+				*entry = index;
+
+			remap[index] = *entry;
+		}
+
+		destination[i] = remap[index];
+	}
+}
+
+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<TriangleKey>(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)
+{
+	using namespace meshopt;
+
+	assert(indices || index_count == vertex_count);
+	assert(vertex_size > 0 && vertex_size <= 256);
+
+	meshopt_Allocator allocator;
+	VertexHasher hasher = {static_cast<const unsigned char*>(vertices), vertex_size, vertex_size};
+
+	return generateVertexRemap(destination, indices, index_count, vertex_count, hasher, allocator);
+}
+
+size_t meshopt_generateVertexRemapMulti(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(indices || index_count == vertex_count);
+	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};
+
+	return generateVertexRemap(destination, indices, index_count, vertex_count, hasher, allocator);
+}
+
+size_t meshopt_generateVertexRemapCustom(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int (*callback)(void*, unsigned int, unsigned int), void* context)
+{
+	using namespace meshopt;
+
+	assert(indices || index_count == vertex_count);
+	assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+	VertexCustomHasher hasher = {vertex_positions, vertex_positions_stride / sizeof(float), callback, context};
+
+	return generateVertexRemap(destination, indices, index_count, vertex_count, hasher, allocator);
+}
+
+void meshopt_remapVertexBuffer(void* destination, const void* vertices, size_t vertex_count, size_t vertex_size, const unsigned int* remap)
+{
+	using namespace meshopt;
+
+	assert(vertex_size > 0 && vertex_size <= 256);
+
+	meshopt_Allocator allocator;
+
+	// support in-place remap
+	if (destination == vertices)
+	{
+		unsigned char* vertices_copy = allocator.allocate<unsigned char>(vertex_count * vertex_size);
+		memcpy(vertices_copy, vertices, vertex_count * vertex_size);
+		vertices = vertices_copy;
+	}
+
+	// specialize the loop for common vertex sizes to ensure memcpy is compiled as an inlined intrinsic
+	switch (vertex_size)
+	{
+	case 4:
+		return remapVertices<4>(destination, vertices, vertex_count, vertex_size, remap);
+
+	case 8:
+		return remapVertices<8>(destination, vertices, vertex_count, vertex_size, remap);
+
+	case 12:
+		return remapVertices<12>(destination, vertices, vertex_count, vertex_size, remap);
+
+	case 16:
+		return remapVertices<16>(destination, vertices, vertex_count, vertex_size, remap);
+
+	default:
+		return remapVertices<0>(destination, vertices, vertex_count, vertex_size, remap);
+	}
+}
+
+void meshopt_remapIndexBuffer(unsigned int* destination, const unsigned int* indices, size_t index_count, const unsigned int* remap)
+{
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices ? indices[i] : unsigned(i);
+		assert(remap[index] != ~0u);
+
+		destination[i] = remap[index];
+	}
+}
+
+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<const unsigned char*>(vertices), vertex_size, vertex_stride};
+
+	unsigned int* remap = allocator.allocate<unsigned int>(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<unsigned int>(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(index_count % 3 == 0);
+	assert(vertex_size > 0 && vertex_size <= 256);
+	assert(vertex_size <= vertex_stride);
+
+	meshopt_Allocator allocator;
+	VertexHasher hasher = {static_cast<const unsigned char*>(vertices), vertex_size, vertex_stride};
+
+	generateShadowBuffer(destination, indices, index_count, vertex_count, hasher, allocator);
+}
+
+void meshopt_generateShadowIndexBufferMulti(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};
+
+	generateShadowBuffer(destination, indices, index_count, vertex_count, hasher, allocator);
+}
+
+void meshopt_generatePositionRemap(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
+{
+	using namespace meshopt;
+
+	assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+	VertexCustomHasher hasher = {vertex_positions, vertex_positions_stride / sizeof(float), NULL, NULL};
+
+	size_t table_size = hashBuckets(vertex_count);
+	unsigned int* table = allocator.allocate<unsigned int>(table_size);
+	memset(table, -1, table_size * sizeof(unsigned int));
+
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		unsigned int* entry = hashLookup(table, table_size, hasher, unsigned(i), ~0u);
+
+		if (*entry == ~0u)
+			*entry = unsigned(i);
+
+		destination[i] = *entry;
+	}
+}
+
+void meshopt_generateAdjacencyIndexBuffer(unsigned int* destination, 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(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	static const int next[4] = {1, 2, 0, 1};
+
+	// build position remap: for each vertex, which other (canonical) vertex does it map to?
+	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
+	buildPositionRemap(remap, vertex_positions, vertex_count, vertex_positions_stride, allocator);
+
+	// build edge set; this stores all triangle edges but we can look these up by any other wedge
+	EdgeHasher edge_hasher = {remap};
+
+	size_t edge_table_size = hashBuckets(index_count);
+	unsigned long long* edge_table = allocator.allocate<unsigned long long>(edge_table_size);
+	unsigned int* edge_vertex_table = allocator.allocate<unsigned int>(edge_table_size);
+
+	memset(edge_table, -1, edge_table_size * sizeof(unsigned long long));
+	memset(edge_vertex_table, -1, edge_table_size * sizeof(unsigned int));
+
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		for (int e = 0; e < 3; ++e)
+		{
+			unsigned int i0 = indices[i + e];
+			unsigned int i1 = indices[i + next[e]];
+			unsigned int i2 = indices[i + next[e + 1]];
+			assert(i0 < vertex_count && i1 < vertex_count && i2 < vertex_count);
+
+			unsigned long long edge = ((unsigned long long)i0 << 32) | i1;
+			unsigned long long* entry = hashLookup(edge_table, edge_table_size, edge_hasher, edge, ~0ull);
+
+			if (*entry == ~0ull)
+			{
+				*entry = edge;
+
+				// store vertex opposite to the edge
+				edge_vertex_table[entry - edge_table] = i2;
+			}
+		}
+	}
+
+	// build resulting index buffer: 6 indices for each input triangle
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		unsigned int patch[6];
+
+		for (int e = 0; e < 3; ++e)
+		{
+			unsigned int i0 = indices[i + e];
+			unsigned int i1 = indices[i + next[e]];
+			assert(i0 < vertex_count && i1 < vertex_count);
+
+			// note: this refers to the opposite edge!
+			unsigned long long edge = ((unsigned long long)i1 << 32) | i0;
+			unsigned long long* oppe = hashLookup(edge_table, edge_table_size, edge_hasher, edge, ~0ull);
+
+			patch[e * 2 + 0] = i0;
+			patch[e * 2 + 1] = (*oppe == ~0ull) ? i0 : edge_vertex_table[oppe - edge_table];
+		}
+
+		memcpy(destination + i * 2, patch, sizeof(patch));
+	}
+}
+
+void meshopt_generateTessellationIndexBuffer(unsigned int* destination, 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(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	static const int next[3] = {1, 2, 0};
+
+	// build position remap: for each vertex, which other (canonical) vertex does it map to?
+	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
+	buildPositionRemap(remap, vertex_positions, vertex_count, vertex_positions_stride, allocator);
+
+	// build edge set; this stores all triangle edges but we can look these up by any other wedge
+	EdgeHasher edge_hasher = {remap};
+
+	size_t edge_table_size = hashBuckets(index_count);
+	unsigned long long* edge_table = allocator.allocate<unsigned long long>(edge_table_size);
+	memset(edge_table, -1, edge_table_size * sizeof(unsigned long long));
+
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		for (int e = 0; e < 3; ++e)
+		{
+			unsigned int i0 = indices[i + e];
+			unsigned int i1 = indices[i + next[e]];
+			assert(i0 < vertex_count && i1 < vertex_count);
+
+			unsigned long long edge = ((unsigned long long)i0 << 32) | i1;
+			unsigned long long* entry = hashLookup(edge_table, edge_table_size, edge_hasher, edge, ~0ull);
+
+			if (*entry == ~0ull)
+				*entry = edge;
+		}
+	}
+
+	// build resulting index buffer: 12 indices for each input triangle
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		unsigned int patch[12];
+
+		for (int e = 0; e < 3; ++e)
+		{
+			unsigned int i0 = indices[i + e];
+			unsigned int i1 = indices[i + next[e]];
+			assert(i0 < vertex_count && i1 < vertex_count);
+
+			// note: this refers to the opposite edge!
+			unsigned long long edge = ((unsigned long long)i1 << 32) | i0;
+			unsigned long long oppe = *hashLookup(edge_table, edge_table_size, edge_hasher, edge, ~0ull);
+
+			// use the same edge if opposite edge doesn't exist (border)
+			oppe = (oppe == ~0ull) ? edge : oppe;
+
+			// triangle index (0, 1, 2)
+			patch[e] = i0;
+
+			// opposite edge (3, 4; 5, 6; 7, 8)
+			patch[3 + e * 2 + 0] = unsigned(oppe);
+			patch[3 + e * 2 + 1] = unsigned(oppe >> 32);
+
+			// dominant vertex (9, 10, 11)
+			patch[9 + e] = remap[i0];
+		}
+
+		memcpy(destination + i * 4, patch, sizeof(patch));
+	}
+}
+
+size_t meshopt_generateProvokingIndexBuffer(unsigned int* destination, unsigned int* reorder, const unsigned int* indices, size_t index_count, size_t vertex_count)
+{
+	assert(index_count % 3 == 0);
+
+	meshopt_Allocator allocator;
+
+	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
+	memset(remap, -1, vertex_count * sizeof(unsigned int));
+
+	// compute vertex valence; this is used to prioritize least used corner
+	// note: we use 8-bit counters for performance; for outlier vertices the valence is incorrect but that just affects the heuristic
+	unsigned char* valence = allocator.allocate<unsigned char>(vertex_count);
+	memset(valence, 0, vertex_count);
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		valence[index]++;
+	}
+
+	unsigned int reorder_offset = 0;
+
+	// assign provoking vertices; leave the rest for the next pass
+	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);
+
+		// try to rotate triangle such that provoking vertex hasn't been seen before
+		// if multiple vertices are new, prioritize the one with least valence
+		// this reduces the risk that a future triangle will have all three vertices seen
+		unsigned int va = remap[a] == ~0u ? valence[a] : ~0u;
+		unsigned int vb = remap[b] == ~0u ? valence[b] : ~0u;
+		unsigned int vc = remap[c] == ~0u ? valence[c] : ~0u;
+
+		if (vb != ~0u && vb <= va && vb <= vc)
+		{
+			// abc -> bca
+			unsigned int t = a;
+			a = b, b = c, c = t;
+		}
+		else if (vc != ~0u && vc <= va && vc <= vb)
+		{
+			// abc -> cab
+			unsigned int t = c;
+			c = b, b = a, a = t;
+		}
+
+		unsigned int newidx = reorder_offset;
+
+		// now remap[a] = ~0u or all three vertices are old
+		// recording remap[a] makes it possible to remap future references to the same index, conserving space
+		if (remap[a] == ~0u)
+			remap[a] = newidx;
+
+		// we need to clone the provoking vertex to get a unique index
+		// if all three are used the choice is arbitrary since no future triangle will be able to reuse any of these
+		reorder[reorder_offset++] = a;
+
+		// note: first vertex is final, the other two will be fixed up in next pass
+		destination[i + 0] = newidx;
+		destination[i + 1] = b;
+		destination[i + 2] = c;
+
+		// update vertex valences for corner heuristic
+		valence[a]--;
+		valence[b]--;
+		valence[c]--;
+	}
+
+	// remap or clone non-provoking vertices (iterating to skip provoking vertices)
+	int step = 1;
+
+	for (size_t i = 1; i < index_count; i += step, step ^= 3)
+	{
+		unsigned int index = destination[i];
+
+		if (remap[index] == ~0u)
+		{
+			// we haven't seen the vertex before as a provoking vertex
+			// to maintain the reference to the original vertex we need to clone it
+			unsigned int newidx = reorder_offset;
+
+			remap[index] = newidx;
+			reorder[reorder_offset++] = index;
+		}
+
+		destination[i] = remap[index];
+	}
+
+	assert(reorder_offset <= vertex_count + index_count / 3);
+	return reorder_offset;
+}

Lib/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 <assert.h>
+#include <float.h>
+#include <math.h>
+#include <string.h>
+
+// 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;
+}

Lib/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 <assert.h>
+#include <math.h>
+#include <string.h>
+
+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<const unsigned char*>(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 <typename T, typename Hash>
+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 <int States>
+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 <int States>
+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<States>(result, index * 4 + i, p0[2], p1[2], p2[2], ac, em[tri[3]], em[tri[4]], em[tri[5]], edgeres);
+	}
+}
+
+template <int States>
+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<States>(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<States>(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<States>(result, index * 4 + 0, level - 1, edgeres, c0, c01, c20, texture);
+	rasterizeOpacityRec<States>(result, index * 4 + 1, level - 1, edgeres, c20, c12, c01, texture);
+	rasterizeOpacityRec<States>(result, index * 4 + 2, level - 1, edgeres, c01, c1, c12, texture);
+	rasterizeOpacityRec<States>(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<unsigned int>(table_size);
+	memset(table, -1, table_size * sizeof(unsigned int));
+
+	TriangleOMM* triangles = allocator.allocate<TriangleOMM>(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<unsigned char>(data_size);
+	memcpy(data_old, data, data_size);
+
+	size_t table_size = hashBuckets3(omm_count);
+	unsigned int* table = allocator.allocate<unsigned int>(table_size);
+	memset(table, -1, table_size * sizeof(unsigned int));
+
+	OMMHasher hasher = {data, offsets, levels, states};
+
+	int* remap = allocator.allocate<int>(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;
+}

Lib/meshoptimizer/src/overdrawoptimizer.cpp

@@ -0,0 +1,333 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+
+// This work is based on:
+// Pedro Sander, Diego Nehab and Joshua Barczak. Fast Triangle Reordering for Vertex Locality and Reduced Overdraw. 2007
+namespace meshopt
+{
+
+static void calculateSortData(float* sort_data, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const unsigned int* clusters, size_t cluster_count)
+{
+	size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
+
+	float mesh_centroid[3] = {};
+
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		const float* p = vertex_positions + vertex_stride_float * i;
+
+		mesh_centroid[0] += p[0];
+		mesh_centroid[1] += p[1];
+		mesh_centroid[2] += p[2];
+	}
+
+	mesh_centroid[0] /= float(vertex_count);
+	mesh_centroid[1] /= float(vertex_count);
+	mesh_centroid[2] /= float(vertex_count);
+
+	for (size_t cluster = 0; cluster < cluster_count; ++cluster)
+	{
+		size_t cluster_begin = clusters[cluster] * 3;
+		size_t cluster_end = (cluster + 1 < cluster_count) ? clusters[cluster + 1] * 3 : index_count;
+		assert(cluster_begin < cluster_end);
+
+		float cluster_area = 0;
+		float cluster_centroid[3] = {};
+		float cluster_normal[3] = {};
+
+		for (size_t i = cluster_begin; i < cluster_end; i += 3)
+		{
+			const float* p0 = vertex_positions + vertex_stride_float * indices[i + 0];
+			const float* p1 = vertex_positions + vertex_stride_float * indices[i + 1];
+			const float* p2 = vertex_positions + vertex_stride_float * indices[i + 2];
+
+			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);
+
+			cluster_centroid[0] += (p0[0] + p1[0] + p2[0]) * (area / 3);
+			cluster_centroid[1] += (p0[1] + p1[1] + p2[1]) * (area / 3);
+			cluster_centroid[2] += (p0[2] + p1[2] + p2[2]) * (area / 3);
+			cluster_normal[0] += normalx;
+			cluster_normal[1] += normaly;
+			cluster_normal[2] += normalz;
+			cluster_area += area;
+		}
+
+		float inv_cluster_area = cluster_area == 0 ? 0 : 1 / cluster_area;
+
+		cluster_centroid[0] *= inv_cluster_area;
+		cluster_centroid[1] *= inv_cluster_area;
+		cluster_centroid[2] *= inv_cluster_area;
+
+		float cluster_normal_length = sqrtf(cluster_normal[0] * cluster_normal[0] + cluster_normal[1] * cluster_normal[1] + cluster_normal[2] * cluster_normal[2]);
+		float inv_cluster_normal_length = cluster_normal_length == 0 ? 0 : 1 / cluster_normal_length;
+
+		cluster_normal[0] *= inv_cluster_normal_length;
+		cluster_normal[1] *= inv_cluster_normal_length;
+		cluster_normal[2] *= inv_cluster_normal_length;
+
+		float centroid_vector[3] = {cluster_centroid[0] - mesh_centroid[0], cluster_centroid[1] - mesh_centroid[1], cluster_centroid[2] - mesh_centroid[2]};
+
+		sort_data[cluster] = centroid_vector[0] * cluster_normal[0] + centroid_vector[1] * cluster_normal[1] + centroid_vector[2] * cluster_normal[2];
+	}
+}
+
+static void calculateSortOrderRadix(unsigned int* sort_order, const float* sort_data, unsigned short* sort_keys, size_t cluster_count)
+{
+	// compute sort data bounds and renormalize, using fixed point snorm
+	float sort_data_max = 1e-3f;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		float dpa = fabsf(sort_data[i]);
+
+		sort_data_max = (sort_data_max < dpa) ? dpa : sort_data_max;
+	}
+
+	const int sort_bits = 11;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		// note that we flip distribution since high dot product should come first
+		float sort_key = 0.5f - 0.5f * (sort_data[i] / sort_data_max);
+
+		sort_keys[i] = meshopt_quantizeUnorm(sort_key, sort_bits) & ((1 << sort_bits) - 1);
+	}
+
+	// fill histogram for counting sort
+	unsigned int histogram[1 << sort_bits];
+	memset(histogram, 0, sizeof(histogram));
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		histogram[sort_keys[i]]++;
+	}
+
+	// compute offsets based on histogram data
+	size_t histogram_sum = 0;
+
+	for (size_t i = 0; i < 1 << sort_bits; ++i)
+	{
+		size_t count = histogram[i];
+		histogram[i] = unsigned(histogram_sum);
+		histogram_sum += count;
+	}
+
+	assert(histogram_sum == cluster_count);
+
+	// compute sort order based on offsets
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		sort_order[histogram[sort_keys[i]]++] = unsigned(i);
+	}
+}
+
+static unsigned int updateCache(unsigned int a, unsigned int b, unsigned int c, unsigned int cache_size, unsigned int* cache_timestamps, unsigned int& timestamp)
+{
+	unsigned int cache_misses = 0;
+
+	// if vertex is not in cache, put it in cache
+	if (timestamp - cache_timestamps[a] > cache_size)
+	{
+		cache_timestamps[a] = timestamp++;
+		cache_misses++;
+	}
+
+	if (timestamp - cache_timestamps[b] > cache_size)
+	{
+		cache_timestamps[b] = timestamp++;
+		cache_misses++;
+	}
+
+	if (timestamp - cache_timestamps[c] > cache_size)
+	{
+		cache_timestamps[c] = timestamp++;
+		cache_misses++;
+	}
+
+	return cache_misses;
+}
+
+static size_t generateHardBoundaries(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, unsigned int cache_size, unsigned int* cache_timestamps)
+{
+	memset(cache_timestamps, 0, vertex_count * sizeof(unsigned int));
+
+	unsigned int timestamp = cache_size + 1;
+
+	size_t face_count = index_count / 3;
+
+	size_t result = 0;
+
+	for (size_t i = 0; i < face_count; ++i)
+	{
+		unsigned int m = updateCache(indices[i * 3 + 0], indices[i * 3 + 1], indices[i * 3 + 2], cache_size, &cache_timestamps[0], timestamp);
+
+		// when all three vertices are not in the cache it's usually relatively safe to assume that this is a new patch in the mesh
+		// that is disjoint from previous vertices; sometimes it might come back to reference existing vertices but that frequently
+		// suggests an inefficiency in the vertex cache optimization algorithm
+		// usually the first triangle has 3 misses unless it's degenerate - thus we make sure the first cluster always starts with 0
+		if (i == 0 || m == 3)
+		{
+			destination[result++] = unsigned(i);
+		}
+	}
+
+	assert(result <= index_count / 3);
+
+	return result;
+}
+
+static size_t generateSoftBoundaries(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, const unsigned int* clusters, size_t cluster_count, unsigned int cache_size, float threshold, unsigned int* cache_timestamps)
+{
+	memset(cache_timestamps, 0, vertex_count * sizeof(unsigned int));
+
+	unsigned int timestamp = 0;
+
+	size_t result = 0;
+
+	for (size_t it = 0; it < cluster_count; ++it)
+	{
+		size_t start = clusters[it];
+		size_t end = (it + 1 < cluster_count) ? clusters[it + 1] : index_count / 3;
+		assert(start < end);
+
+		// reset cache
+		timestamp += cache_size + 1;
+
+		// measure cluster ACMR
+		unsigned int cluster_misses = 0;
+
+		for (size_t i = start; i < end; ++i)
+		{
+			unsigned int m = updateCache(indices[i * 3 + 0], indices[i * 3 + 1], indices[i * 3 + 2], cache_size, &cache_timestamps[0], timestamp);
+
+			cluster_misses += m;
+		}
+
+		float cluster_threshold = threshold * (float(cluster_misses) / float(end - start));
+
+		// first cluster always starts from the hard cluster boundary
+		destination[result++] = unsigned(start);
+
+		// reset cache
+		timestamp += cache_size + 1;
+
+		unsigned int running_misses = 0;
+		unsigned int running_faces = 0;
+
+		for (size_t i = start; i < end; ++i)
+		{
+			unsigned int m = updateCache(indices[i * 3 + 0], indices[i * 3 + 1], indices[i * 3 + 2], cache_size, &cache_timestamps[0], timestamp);
+
+			running_misses += m;
+			running_faces += 1;
+
+			if (float(running_misses) / float(running_faces) <= cluster_threshold)
+			{
+				// we have reached the target ACMR with the current triangle so we need to start a new cluster on the next one
+				// note that this may mean that we add 'end` to destination for the last triangle, which will imply that the last
+				// cluster is empty; however, the 'pop_back' after the loop will clean it up
+				destination[result++] = unsigned(i + 1);
+
+				// reset cache
+				timestamp += cache_size + 1;
+
+				running_misses = 0;
+				running_faces = 0;
+			}
+		}
+
+		// each time we reach the target ACMR we flush the cluster
+		// this means that the last cluster is by definition not very good - there are frequent cases where we are left with a few triangles
+		// in the last cluster, producing a very bad ACMR and significantly penalizing the overall results
+		// thus we remove the last cluster boundary, merging the last complete cluster with the last incomplete one
+		// there are sometimes cases when the last cluster is actually good enough - in which case the code above would have added 'end'
+		// to the cluster boundary array which we need to remove anyway - this code will do that automatically
+		if (destination[result - 1] != start)
+		{
+			result--;
+		}
+	}
+
+	assert(result >= cluster_count);
+	assert(result <= index_count / 3);
+
+	return result;
+}
+
+} // namespace meshopt
+
+void meshopt_optimizeOverdraw(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, float threshold)
+{
+	using namespace meshopt;
+
+	assert(index_count % 3 == 0);
+	assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	// guard for empty meshes
+	if (index_count == 0 || vertex_count == 0)
+		return;
+
+	// support in-place optimization
+	if (destination == indices)
+	{
+		unsigned int* indices_copy = allocator.allocate<unsigned int>(index_count);
+		memcpy(indices_copy, indices, index_count * sizeof(unsigned int));
+		indices = indices_copy;
+	}
+
+	unsigned int cache_size = 16;
+
+	unsigned int* cache_timestamps = allocator.allocate<unsigned int>(vertex_count);
+
+	// generate hard boundaries from full-triangle cache misses
+	unsigned int* hard_clusters = allocator.allocate<unsigned int>(index_count / 3);
+	size_t hard_cluster_count = generateHardBoundaries(hard_clusters, indices, index_count, vertex_count, cache_size, cache_timestamps);
+
+	// generate soft boundaries
+	unsigned int* soft_clusters = allocator.allocate<unsigned int>(index_count / 3 + 1);
+	size_t soft_cluster_count = generateSoftBoundaries(soft_clusters, indices, index_count, vertex_count, hard_clusters, hard_cluster_count, cache_size, threshold, cache_timestamps);
+
+	const unsigned int* clusters = soft_clusters;
+	size_t cluster_count = soft_cluster_count;
+
+	// fill sort data
+	float* sort_data = allocator.allocate<float>(cluster_count);
+	calculateSortData(sort_data, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride, clusters, cluster_count);
+
+	// sort clusters using sort data
+	unsigned short* sort_keys = allocator.allocate<unsigned short>(cluster_count);
+	unsigned int* sort_order = allocator.allocate<unsigned int>(cluster_count);
+	calculateSortOrderRadix(sort_order, sort_data, sort_keys, cluster_count);
+
+	// fill output buffer
+	size_t offset = 0;
+
+	for (size_t it = 0; it < cluster_count; ++it)
+	{
+		unsigned int cluster = sort_order[it];
+		assert(cluster < cluster_count);
+
+		size_t cluster_begin = clusters[cluster] * 3;
+		size_t cluster_end = (cluster + 1 < cluster_count) ? clusters[cluster + 1] * 3 : index_count;
+		assert(cluster_begin < cluster_end);
+
+		memcpy(destination + offset, indices + cluster_begin, (cluster_end - cluster_begin) * sizeof(unsigned int));
+		offset += cluster_end - cluster_begin;
+	}
+
+	assert(offset == index_count);
+}

Lib/meshoptimizer/src/partition.cpp

@@ -0,0 +1,621 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+
+// This work is based on:
+// Takio Kurita. An efficient agglomerative clustering algorithm using a heap. 1991
+namespace meshopt
+{
+
+// To avoid excessive recursion for malformed inputs, we switch to bisection after some depth
+const int kMergeDepthCutoff = 40;
+
+struct ClusterAdjacency
+{
+	unsigned int* offsets;
+	unsigned int* clusters;
+	unsigned int* shared;
+};
+
+static void filterClusterIndices(unsigned int* data, unsigned int* offsets, const unsigned int* cluster_indices, const unsigned int* cluster_index_counts, size_t cluster_count, unsigned char* used, size_t vertex_count, size_t total_index_count)
+{
+	(void)vertex_count;
+	(void)total_index_count;
+
+	size_t cluster_start = 0;
+	size_t cluster_write = 0;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		offsets[i] = unsigned(cluster_write);
+
+		// copy cluster indices, skipping duplicates
+		for (size_t j = 0; j < cluster_index_counts[i]; ++j)
+		{
+			unsigned int v = cluster_indices[cluster_start + j];
+			assert(v < vertex_count);
+
+			data[cluster_write] = v;
+			cluster_write += 1 - used[v];
+			used[v] = 1;
+		}
+
+		// reset used flags for the next cluster
+		for (size_t j = offsets[i]; j < cluster_write; ++j)
+			used[data[j]] = 0;
+
+		cluster_start += cluster_index_counts[i];
+	}
+
+	assert(cluster_start == total_index_count);
+	assert(cluster_write <= total_index_count);
+	offsets[cluster_count] = unsigned(cluster_write);
+}
+
+static float computeClusterBounds(const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_positions_stride, float* out_center)
+{
+	size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
+
+	float center[3] = {0, 0, 0};
+
+	// approximate center of the cluster by averaging all vertex positions
+	for (size_t j = 0; j < index_count; ++j)
+	{
+		const float* p = vertex_positions + indices[j] * vertex_stride_float;
+
+		center[0] += p[0];
+		center[1] += p[1];
+		center[2] += p[2];
+	}
+
+	// note: technically clusters can't be empty per meshopt_partitionCluster but we check for a division by zero in case that changes
+	if (index_count)
+	{
+		center[0] /= float(index_count);
+		center[1] /= float(index_count);
+		center[2] /= float(index_count);
+	}
+
+	// compute radius of the bounding sphere for each cluster
+	float radiussq = 0;
+
+	for (size_t j = 0; j < index_count; ++j)
+	{
+		const float* p = vertex_positions + indices[j] * vertex_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]);
+
+		radiussq = radiussq < d2 ? d2 : radiussq;
+	}
+
+	memcpy(out_center, center, sizeof(center));
+	return sqrtf(radiussq);
+}
+
+static void buildClusterAdjacency(ClusterAdjacency& adjacency, const unsigned int* cluster_indices, const unsigned int* cluster_offsets, size_t cluster_count, size_t vertex_count, meshopt_Allocator& allocator)
+{
+	unsigned int* ref_offsets = allocator.allocate<unsigned int>(vertex_count + 1);
+
+	// compute number of clusters referenced by each vertex
+	memset(ref_offsets, 0, vertex_count * sizeof(unsigned int));
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		for (size_t j = cluster_offsets[i]; j < cluster_offsets[i + 1]; ++j)
+			ref_offsets[cluster_indices[j]]++;
+	}
+
+	// compute (worst-case) number of adjacent clusters for each cluster
+	size_t total_adjacency = 0;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		size_t count = 0;
+
+		// worst case is every vertex has a disjoint cluster list
+		for (size_t j = cluster_offsets[i]; j < cluster_offsets[i + 1]; ++j)
+			count += ref_offsets[cluster_indices[j]] - 1;
+
+		// ... but only every other cluster can be adjacent in the end
+		total_adjacency += count < cluster_count - 1 ? count : cluster_count - 1;
+	}
+
+	// we can now allocate adjacency buffers
+	adjacency.offsets = allocator.allocate<unsigned int>(cluster_count + 1);
+	adjacency.clusters = allocator.allocate<unsigned int>(total_adjacency);
+	adjacency.shared = allocator.allocate<unsigned int>(total_adjacency);
+
+	// convert ref counts to offsets
+	size_t total_refs = 0;
+
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		size_t count = ref_offsets[i];
+		ref_offsets[i] = unsigned(total_refs);
+		total_refs += count;
+	}
+
+	unsigned int* ref_data = allocator.allocate<unsigned int>(total_refs);
+
+	// fill cluster refs for each vertex
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		for (size_t j = cluster_offsets[i]; j < cluster_offsets[i + 1]; ++j)
+			ref_data[ref_offsets[cluster_indices[j]]++] = unsigned(i);
+	}
+
+	// after the previous pass, ref_offsets contain the end of the data for each vertex; shift it forward to get the start
+	memmove(ref_offsets + 1, ref_offsets, vertex_count * sizeof(unsigned int));
+	ref_offsets[0] = 0;
+
+	unsigned int* slots = allocator.allocate<unsigned int>(cluster_count);
+	memset(slots, -1, cluster_count * sizeof(unsigned int));
+
+	// fill cluster adjacency for each cluster...
+	adjacency.offsets[0] = 0;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		unsigned int* adj = adjacency.clusters + adjacency.offsets[i];
+		unsigned int* shd = adjacency.shared + adjacency.offsets[i];
+		size_t count = 0;
+
+		for (size_t j = cluster_offsets[i]; j < cluster_offsets[i + 1]; ++j)
+		{
+			unsigned int v = cluster_indices[j];
+
+			// merge the entire cluster list of each vertex into current list
+			for (size_t k = ref_offsets[v]; k < ref_offsets[v + 1]; ++k)
+			{
+				unsigned int c = ref_data[k];
+				assert(c < cluster_count);
+
+				if (c == unsigned(i))
+					continue;
+
+				// if the cluster is already in the list, increment the shared count
+				// we do not need to clear slots across clusters due to adj[slot] confirmation check
+				unsigned int slot = slots[c];
+				if (slot < count && adj[slot] == c)
+					shd[slot]++;
+				else
+				{
+					slots[c] = unsigned(count);
+					adj[count] = c;
+					shd[count] = 1;
+					count++;
+				}
+			}
+		}
+
+		// mark the end of the adjacency list; the next cluster will start there as well
+		adjacency.offsets[i + 1] = adjacency.offsets[i] + unsigned(count);
+	}
+
+	assert(adjacency.offsets[cluster_count] <= total_adjacency);
+
+	// ref_offsets can't be deallocated as it was allocated before adjacency
+	allocator.deallocate(slots);
+	allocator.deallocate(ref_data);
+}
+
+struct ClusterGroup
+{
+	int group;
+	int next;
+	unsigned int size; // 0 unless root
+	unsigned int vertices;
+
+	float center[3];
+	float radius;
+};
+
+struct GroupOrder
+{
+	unsigned int id;
+	int order;
+};
+
+static void heapPush(GroupOrder* heap, size_t size, GroupOrder item)
+{
+	// insert a new element at the end (breaks heap invariant)
+	heap[size++] = item;
+
+	// bubble up the new element to its correct position
+	size_t i = size - 1;
+	while (i > 0 && heap[i].order < heap[(i - 1) / 2].order)
+	{
+		size_t p = (i - 1) / 2;
+
+		GroupOrder temp = heap[i];
+		heap[i] = heap[p];
+		heap[p] = temp;
+		i = p;
+	}
+}
+
+static GroupOrder heapPop(GroupOrder* heap, size_t size)
+{
+	assert(size > 0);
+	GroupOrder top = heap[0];
+
+	// move the last element to the top (breaks heap invariant)
+	heap[0] = heap[--size];
+
+	// bubble down the new top element to its correct position
+	size_t i = 0;
+	while (i * 2 + 1 < size)
+	{
+		// find the smallest child
+		size_t j = i * 2 + 1;
+		j += (j + 1 < size && heap[j + 1].order < heap[j].order);
+
+		// if the parent is already smaller than both children, we're done
+		if (heap[j].order >= heap[i].order)
+			break;
+
+		// otherwise, swap the parent and child and continue
+		GroupOrder temp = heap[i];
+		heap[i] = heap[j];
+		heap[j] = temp;
+		i = j;
+	}
+
+	return top;
+}
+
+static void mergeBounds(ClusterGroup& target, const ClusterGroup& source)
+{
+	float r1 = target.radius, r2 = source.radius;
+	float dx = source.center[0] - target.center[0], dy = source.center[1] - target.center[1], dz = source.center[2] - target.center[2];
+	float d = sqrtf(dx * dx + dy * dy + dz * dz);
+
+	if (d + r1 < r2)
+	{
+		target.center[0] = source.center[0];
+		target.center[1] = source.center[1];
+		target.center[2] = source.center[2];
+		target.radius = source.radius;
+		return;
+	}
+
+	if (d + r2 > r1)
+	{
+		float k = d > 0 ? (d + r2 - r1) / (2 * d) : 0.f;
+
+		target.center[0] += dx * k;
+		target.center[1] += dy * k;
+		target.center[2] += dz * k;
+		target.radius = (d + r2 + r1) / 2;
+	}
+}
+
+static float boundsScore(const ClusterGroup& target, const ClusterGroup& source)
+{
+	float r1 = target.radius, r2 = source.radius;
+	float dx = source.center[0] - target.center[0], dy = source.center[1] - target.center[1], dz = source.center[2] - target.center[2];
+	float d = sqrtf(dx * dx + dy * dy + dz * dz);
+
+	float mr = d + r1 < r2 ? r2 : (d + r2 < r1 ? r1 : (d + r2 + r1) / 2);
+
+	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, unsigned int* shared_acc, unsigned int& out_shared)
+{
+	assert(groups[id].size > 0);
+
+	float group_rsqrt = 1.f / sqrtf(float(int(groups[id].vertices)));
+
+	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)
+				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;
+
+			float other_rsqrt = 1.f / sqrtf(float(int(groups[other].vertices)));
+
+			// normalize shared count by the expected boundary of each group (+ keeps scoring symmetric)
+			float score = float(int(shared)) * (group_rsqrt + other_rsqrt);
+
+			// incorporate spatial score to favor merging nearby groups
+			if (use_bounds)
+				score *= 1.f + 0.4f * boundsScore(groups[id], groups[other]);
+
+			if (score > best_score)
+			{
+				best_group = other;
+				best_score = score;
+				best_shared = shared;
+			}
+		}
+
+	out_shared = best_shared;
+	return best_group;
+}
+
+static void mergeLeaf(ClusterGroup* groups, unsigned int* order, size_t count, size_t target_partition_size, size_t max_partition_size)
+{
+	for (size_t i = 0; i < count; ++i)
+	{
+		unsigned int id = order[i];
+		if (groups[id].size == 0 || groups[id].size >= target_partition_size)
+			continue;
+
+		float best_score = -1.f;
+		int best_group = -1;
+
+		for (size_t j = 0; j < count; ++j)
+		{
+			unsigned int other = order[j];
+			if (id == other || groups[other].size == 0)
+				continue;
+
+			if (groups[id].size + groups[other].size > max_partition_size)
+				continue;
+
+			// favor merging nearby groups
+			float score = boundsScore(groups[id], groups[other]);
+
+			if (score > best_score)
+			{
+				best_score = score;
+				best_group = other;
+			}
+		}
+
+		// merge id *into* best_group; that way, we may merge more groups into the same best_group, maximizing the chance of reaching target
+		if (best_group != -1)
+		{
+			// combine groups by linking them together
+			unsigned int tail = best_group;
+			while (groups[tail].next >= 0)
+				tail = groups[tail].next;
+
+			groups[tail].next = id;
+
+			// update group sizes; note, we omit vertices update for simplicity as it's not used for spatial merge
+			groups[best_group].size += groups[id].size;
+			groups[id].size = 0;
+
+			// merge bounding spheres
+			mergeBounds(groups[best_group], groups[id]);
+			groups[id].radius = 0.f;
+		}
+	}
+}
+
+static size_t mergePartition(unsigned int* order, size_t count, const ClusterGroup* groups, int axis, float pivot)
+{
+	size_t m = 0;
+
+	// invariant: elements in range [0, m) are < pivot, elements in range [m, i) are >= pivot
+	for (size_t i = 0; i < count; ++i)
+	{
+		float v = groups[order[i]].center[axis];
+
+		// swap(m, i) unconditionally
+		unsigned int t = order[m];
+		order[m] = order[i];
+		order[i] = t;
+
+		// when v >= pivot, we swap i with m without advancing it, preserving invariants
+		m += v < pivot;
+	}
+
+	return m;
+}
+
+static void mergeSpatial(ClusterGroup* groups, unsigned int* order, size_t count, size_t target_partition_size, size_t max_partition_size, size_t leaf_size, int depth)
+{
+	size_t total = 0;
+	for (size_t i = 0; i < count; ++i)
+		total += groups[order[i]].size;
+
+	if (total <= max_partition_size || count <= leaf_size)
+		return mergeLeaf(groups, order, count, target_partition_size, max_partition_size);
+
+	float mean[3] = {};
+	float vars[3] = {};
+	float runc = 1, runs = 1;
+
+	// gather statistics on the points in the subtree using Welford's algorithm
+	for (size_t i = 0; i < count; ++i, runc += 1.f, runs = 1.f / runc)
+	{
+		const float* point = groups[order[i]].center;
+
+		for (int k = 0; k < 3; ++k)
+		{
+			float delta = point[k] - mean[k];
+			mean[k] += delta * runs;
+			vars[k] += delta * (point[k] - mean[k]);
+		}
+	}
+
+	// split axis is one where the variance is largest
+	int axis = (vars[0] >= vars[1] && vars[0] >= vars[2]) ? 0 : (vars[1] >= vars[2] ? 1 : 2);
+
+	float split = mean[axis];
+	size_t middle = mergePartition(order, count, groups, axis, split);
+
+	// enforce balance for degenerate partitions
+	// this also ensures recursion depth is bounded on pathological inputs
+	if (middle <= leaf_size / 2 || count - middle <= leaf_size / 2 || depth >= kMergeDepthCutoff)
+		middle = count / 2;
+
+	// recursion depth is logarithmic and bounded due to max depth check above
+	mergeSpatial(groups, order, middle, target_partition_size, max_partition_size, leaf_size, depth + 1);
+	mergeSpatial(groups, order + middle, count - middle, target_partition_size, max_partition_size, leaf_size, depth + 1);
+}
+
+} // namespace meshopt
+
+size_t meshopt_partitionClusters(unsigned int* destination, const unsigned int* 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)
+{
+	using namespace meshopt;
+
+	assert((vertex_positions == NULL || vertex_positions_stride >= 12) && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+	assert(target_partition_size > 0);
+
+	size_t max_partition_size = target_partition_size + target_partition_size / 3;
+
+	meshopt_Allocator allocator;
+
+	unsigned char* used = allocator.allocate<unsigned char>(vertex_count);
+	memset(used, 0, vertex_count);
+
+	unsigned int* cluster_newindices = allocator.allocate<unsigned int>(total_index_count);
+	unsigned int* cluster_offsets = allocator.allocate<unsigned int>(cluster_count + 1);
+
+	// make new cluster index list that filters out duplicate indices
+	filterClusterIndices(cluster_newindices, cluster_offsets, cluster_indices, cluster_index_counts, cluster_count, used, vertex_count, total_index_count);
+	cluster_indices = cluster_newindices;
+
+	// build cluster adjacency along with edge weights (shared vertex count)
+	ClusterAdjacency adjacency = {};
+	buildClusterAdjacency(adjacency, cluster_indices, cluster_offsets, cluster_count, vertex_count, allocator);
+
+	ClusterGroup* groups = allocator.allocate<ClusterGroup>(cluster_count);
+	memset(groups, 0, sizeof(ClusterGroup) * cluster_count);
+
+	GroupOrder* order = allocator.allocate<GroupOrder>(cluster_count);
+	size_t pending = 0;
+
+	unsigned int* shared_acc = allocator.allocate<unsigned int>(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)
+	{
+		groups[i].group = int(i);
+		groups[i].next = -1;
+		groups[i].size = 1;
+		groups[i].vertices = cluster_offsets[i + 1] - cluster_offsets[i];
+		assert(groups[i].vertices > 0);
+
+		// compute bounding sphere for each cluster if positions are provided
+		if (vertex_positions)
+			groups[i].radius = computeClusterBounds(cluster_indices + cluster_offsets[i], cluster_offsets[i + 1] - cluster_offsets[i], vertex_positions, vertex_positions_stride, groups[i].center);
+
+		GroupOrder item = {};
+		item.id = unsigned(i);
+		item.order = groups[i].vertices;
+
+		heapPush(order, pending++, item);
+	}
+
+	// iteratively merge the smallest group with the best group
+	while (pending)
+	{
+		GroupOrder top = heapPop(order, pending--);
+
+		// this group was merged into another group earlier
+		if (groups[top.id].size == 0)
+			continue;
+
+		// disassociate clusters from the group to prevent them from being merged again; we will re-associate them if the group is reinserted
+		for (int i = top.id; i >= 0; i = groups[i].next)
+		{
+			assert(groups[i].group == int(top.id));
+			groups[i].group = -1;
+		}
+
+		// the group is large enough, emit as is
+		if (groups[top.id].size >= target_partition_size)
+			continue;
+
+		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;
+
+		// combine groups by linking them together
+		unsigned int tail = top.id;
+		while (groups[tail].next >= 0)
+			tail = groups[tail].next;
+
+		groups[tail].next = best_group;
+
+		// 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 > best_shared) ? groups[top.id].vertices - best_shared : 1;
+
+		groups[best_group].size = 0;
+		groups[best_group].vertices = 0;
+
+		// merge bounding spheres if bounds are available
+		if (vertex_positions)
+		{
+			mergeBounds(groups[top.id], groups[best_group]);
+			groups[best_group].radius = 0;
+		}
+
+		// re-associate all clusters back to the merged group
+		for (int i = top.id; i >= 0; i = groups[i].next)
+			groups[i].group = int(top.id);
+
+		top.order = groups[top.id].vertices;
+		heapPush(order, pending++, top);
+	}
+
+	// if vertex positions are provided, we do a final pass to see if we can merge small groups based on spatial locality alone
+	if (vertex_positions)
+	{
+		unsigned int* merge_order = reinterpret_cast<unsigned int*>(order);
+		size_t merge_offset = 0;
+
+		for (size_t i = 0; i < cluster_count; ++i)
+			if (groups[i].size)
+				merge_order[merge_offset++] = unsigned(i);
+
+		mergeSpatial(groups, merge_order, merge_offset, target_partition_size, max_partition_size, /* leaf_size= */ 8, 0);
+	}
+
+	// output each remaining group
+	size_t next_group = 0;
+
+	for (size_t i = 0; i < cluster_count; ++i)
+	{
+		if (groups[i].size == 0)
+			continue;
+
+		for (int j = int(i); j >= 0; j = groups[j].next)
+			destination[j] = unsigned(next_group);
+
+		next_group++;
+	}
+
+	assert(next_group <= cluster_count);
+	return next_group;
+}

Lib/meshoptimizer/src/quantization.cpp

@@ -0,0 +1,127 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <float.h>
+#include <math.h>
+
+union FloatBits
+{
+	float f;
+	unsigned int ui;
+};
+
+unsigned short meshopt_quantizeHalf(float v)
+{
+	FloatBits u = {v};
+	unsigned int ui = u.ui;
+
+	int s = (ui >> 16) & 0x8000;
+	int em = ui & 0x7fffffff;
+
+	// bias exponent and round to nearest; 112 is relative exponent bias (127-15)
+	int h = (em - (112 << 23) + (1 << 12)) >> 13;
+
+	// underflow: flush to zero; 113 encodes exponent -14
+	h = (em < (113 << 23)) ? 0 : h;
+
+	// overflow: infinity; 143 encodes exponent 16
+	h = (em >= (143 << 23)) ? 0x7c00 : h;
+
+	// NaN; note that we convert all types of NaN to qNaN
+	h = (em > (255 << 23)) ? 0x7e00 : h;
+
+	return (unsigned short)(s | h);
+}
+
+float meshopt_quantizeFloat(float v, int N)
+{
+	assert(N >= 0 && N <= 23);
+
+	FloatBits u = {v};
+	unsigned int ui = u.ui;
+
+	const int mask = (1 << (23 - N)) - 1;
+	const int round = (1 << (23 - N)) >> 1;
+
+	int e = ui & 0x7f800000;
+	unsigned int rui = (ui + round) & ~mask;
+
+	// round all numbers except inf/nan; this is important to make sure nan doesn't overflow into -0
+	ui = e == 0x7f800000 ? ui : rui;
+
+	// flush denormals to zero
+	ui = e == 0 ? 0 : ui;
+
+	u.ui = ui;
+	return u.f;
+}
+
+float meshopt_dequantizeHalf(unsigned short h)
+{
+	unsigned int s = unsigned(h & 0x8000) << 16;
+	int em = h & 0x7fff;
+
+	// bias exponent and pad mantissa with 0; 112 is relative exponent bias (127-15)
+	int r = (em + (112 << 10)) << 13;
+
+	// denormal: flush to zero
+	r = (em < (1 << 10)) ? 0 : r;
+
+	// infinity/NaN; note that we preserve NaN payload as a byproduct of unifying inf/nan cases
+	// 112 is an exponent bias fixup; since we already applied it once, applying it twice converts 31 to 255
+	r += (em >= (31 << 10)) ? (112 << 23) : 0;
+
+	FloatBits u;
+	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;
+}

Lib/meshoptimizer/src/rasterizer.cpp

@@ -0,0 +1,289 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <float.h>
+#include <string.h>
+
+// This work is based on:
+// Nicolas Capens. Advanced Rasterization. 2004
+namespace meshopt
+{
+
+const int kViewport = 256;
+
+struct OverdrawBuffer
+{
+	float z[kViewport][kViewport][2];
+	unsigned int overdraw[kViewport][kViewport][2];
+};
+
+static float computeDepthGradients(float& dzdx, float& dzdy, float x1, float y1, float z1, float x2, float y2, float z2, float x3, float y3, float z3)
+{
+	// z2 = z1 + dzdx * (x2 - x1) + dzdy * (y2 - y1)
+	// z3 = z1 + dzdx * (x3 - x1) + dzdy * (y3 - y1)
+	// (x2-x1 y2-y1)(dzdx) = (z2-z1)
+	// (x3-x1 y3-y1)(dzdy)   (z3-z1)
+	// we'll solve it with Cramer's rule
+	float det = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1);
+	float invdet = (det == 0) ? 0 : 1 / det;
+
+	dzdx = ((z2 - z1) * (y3 - y1) - (y2 - y1) * (z3 - z1)) * invdet;
+	dzdy = ((x2 - x1) * (z3 - z1) - (z2 - z1) * (x3 - x1)) * invdet;
+
+	return det;
+}
+
+// half-space fixed point triangle rasterizer
+static void rasterize(OverdrawBuffer* buffer, float v1x, float v1y, float v1z, float v2x, float v2y, float v2z, float v3x, float v3y, float v3z)
+{
+	// compute depth gradients
+	float DZx, DZy;
+	float det = computeDepthGradients(DZx, DZy, v1x, v1y, v1z, v2x, v2y, v2z, v3x, v3y, v3z);
+	int sign = det > 0;
+
+	// flip backfacing triangles to simplify rasterization logic
+	if (sign)
+	{
+		// flipping v2 & v3 preserves depth gradients since they're based on v1; only v1z is used below
+		float t;
+		t = v2x, v2x = v3x, v3x = t;
+		t = v2y, v2y = v3y, v3y = t;
+
+		// flip depth since we rasterize backfacing triangles to second buffer with reverse Z; only v1z is used below
+		v1z = kViewport - v1z;
+		DZx = -DZx;
+		DZy = -DZy;
+	}
+
+	// coordinates, 28.4 fixed point
+	int X1 = int(16.0f * v1x + 0.5f);
+	int X2 = int(16.0f * v2x + 0.5f);
+	int X3 = int(16.0f * v3x + 0.5f);
+
+	int Y1 = int(16.0f * v1y + 0.5f);
+	int Y2 = int(16.0f * v2y + 0.5f);
+	int Y3 = int(16.0f * v3y + 0.5f);
+
+	// bounding rectangle, clipped against viewport
+	// since we rasterize pixels with covered centers, min >0.5 should round up
+	// as for max, due to top-left filling convention we will never rasterize right/bottom edges
+	// so max >= 0.5 should round down for inclusive bounds, and up for exclusive (in our case)
+	int minx = X1 < X2 ? X1 : X2;
+	minx = minx < X3 ? minx : X3;
+	minx = (minx + 7) >> 4;
+	minx = minx < 0 ? 0 : minx;
+
+	int miny = Y1 < Y2 ? Y1 : Y2;
+	miny = miny < Y3 ? miny : Y3;
+	miny = (miny + 7) >> 4;
+	miny = miny < 0 ? 0 : miny;
+
+	int maxx = X1 > X2 ? X1 : X2;
+	maxx = maxx > X3 ? maxx : X3;
+	maxx = (maxx + 7) >> 4;
+	maxx = maxx > kViewport ? kViewport : maxx;
+
+	int maxy = Y1 > Y2 ? Y1 : Y2;
+	maxy = maxy > Y3 ? maxy : Y3;
+	maxy = (maxy + 7) >> 4;
+	maxy = maxy > kViewport ? kViewport : maxy;
+
+	// deltas, 28.4 fixed point
+	int DX12 = X1 - X2;
+	int DX23 = X2 - X3;
+	int DX31 = X3 - X1;
+
+	int DY12 = Y1 - Y2;
+	int DY23 = Y2 - Y3;
+	int DY31 = Y3 - Y1;
+
+	// fill convention correction
+	int TL1 = DY12 < 0 || (DY12 == 0 && DX12 > 0);
+	int TL2 = DY23 < 0 || (DY23 == 0 && DX23 > 0);
+	int TL3 = DY31 < 0 || (DY31 == 0 && DX31 > 0);
+
+	// half edge equations, 24.8 fixed point
+	// note that we offset minx/miny by half pixel since we want to rasterize pixels with covered centers
+	int FX = (minx << 4) + 8;
+	int FY = (miny << 4) + 8;
+	int CY1 = DX12 * (FY - Y1) - DY12 * (FX - X1) + TL1 - 1;
+	int CY2 = DX23 * (FY - Y2) - DY23 * (FX - X2) + TL2 - 1;
+	int CY3 = DX31 * (FY - Y3) - DY31 * (FX - X3) + TL3 - 1;
+	float ZY = v1z + (DZx * float(FX - X1) + DZy * float(FY - Y1)) * (1 / 16.f);
+
+	for (int y = miny; y < maxy; y++)
+	{
+		int CX1 = CY1;
+		int CX2 = CY2;
+		int CX3 = CY3;
+		float ZX = ZY;
+
+		for (int x = minx; x < maxx; x++)
+		{
+			// check if all CXn are non-negative
+			if ((CX1 | CX2 | CX3) >= 0)
+			{
+				if (ZX >= buffer->z[y][x][sign])
+				{
+					buffer->z[y][x][sign] = ZX;
+					buffer->overdraw[y][x][sign]++;
+				}
+			}
+
+			// signed left shift is UB for negative numbers so use unsigned-signed casts
+			CX1 -= int(unsigned(DY12) << 4);
+			CX2 -= int(unsigned(DY23) << 4);
+			CX3 -= int(unsigned(DY31) << 4);
+			ZX += DZx;
+		}
+
+		// signed left shift is UB for negative numbers so use unsigned-signed casts
+		CY1 += int(unsigned(DX12) << 4);
+		CY2 += int(unsigned(DX23) << 4);
+		CY3 += int(unsigned(DX31) << 4);
+		ZY += DZy;
+	}
+}
+
+static float transformTriangles(float* triangles, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
+{
+	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 + 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]);
+
+	float scale = extent == 0.f ? 0.f : kViewport / extent;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		const float* v = vertex_positions + index * vertex_stride_float;
+
+		triangles[i * 3 + 0] = (v[0] - minv[0]) * scale;
+		triangles[i * 3 + 1] = (v[1] - minv[1]) * scale;
+		triangles[i * 3 + 2] = (v[2] - minv[2]) * scale;
+	}
+
+	return extent;
+}
+
+static void rasterizeTriangles(OverdrawBuffer* buffer, const float* triangles, size_t index_count, int axis)
+{
+	for (size_t i = 0; i < index_count; i += 3)
+	{
+		const float* vn0 = &triangles[3 * (i + 0)];
+		const float* vn1 = &triangles[3 * (i + 1)];
+		const float* vn2 = &triangles[3 * (i + 2)];
+
+		switch (axis)
+		{
+		case 0:
+			rasterize(buffer, vn0[2], vn0[1], vn0[0], vn1[2], vn1[1], vn1[0], vn2[2], vn2[1], vn2[0]);
+			break;
+		case 1:
+			rasterize(buffer, vn0[0], vn0[2], vn0[1], vn1[0], vn1[2], vn1[1], vn2[0], vn2[2], vn2[1]);
+			break;
+		case 2:
+			rasterize(buffer, vn0[1], vn0[0], vn0[2], vn1[1], vn1[0], vn1[2], vn2[1], vn2[0], vn2[2]);
+			break;
+		}
+	}
+}
+
+} // namespace meshopt
+
+meshopt_OverdrawStatistics meshopt_analyzeOverdraw(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(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	meshopt_OverdrawStatistics result = {};
+
+	float* triangles = allocator.allocate<float>(index_count * 3);
+	transformTriangles(triangles, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride);
+
+	OverdrawBuffer* buffer = allocator.allocate<OverdrawBuffer>(1);
+
+	for (int axis = 0; axis < 3; ++axis)
+	{
+		memset(buffer, 0, sizeof(OverdrawBuffer));
+		rasterizeTriangles(buffer, triangles, index_count, axis);
+
+		for (int y = 0; y < kViewport; ++y)
+			for (int x = 0; x < kViewport; ++x)
+				for (int s = 0; s < 2; ++s)
+				{
+					unsigned int overdraw = buffer->overdraw[y][x][s];
+
+					result.pixels_covered += overdraw > 0;
+					result.pixels_shaded += overdraw;
+				}
+	}
+
+	result.overdraw = result.pixels_covered ? float(result.pixels_shaded) / float(result.pixels_covered) : 0.f;
+
+	return result;
+}
+
+meshopt_CoverageStatistics meshopt_analyzeCoverage(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(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	meshopt_CoverageStatistics result = {};
+
+	float* triangles = allocator.allocate<float>(index_count * 3);
+	float extent = transformTriangles(triangles, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride);
+
+	OverdrawBuffer* buffer = allocator.allocate<OverdrawBuffer>(1);
+
+	for (int axis = 0; axis < 3; ++axis)
+	{
+		memset(buffer, 0, sizeof(OverdrawBuffer));
+		rasterizeTriangles(buffer, triangles, index_count, axis);
+
+		unsigned int covered = 0;
+
+		for (int y = 0; y < kViewport; ++y)
+			for (int x = 0; x < kViewport; ++x)
+				covered += (buffer->overdraw[y][x][0] | buffer->overdraw[y][x][1]) > 0;
+
+		result.coverage[axis] = float(covered) / float(kViewport * kViewport);
+	}
+
+	result.extent = extent;
+
+	return result;
+}

Lib/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 <assert.h>
+#include <float.h>
+#include <math.h>
+#include <string.h>
+
+#ifndef TRACE
+#define TRACE 0
+#endif
+
+#if TRACE
+#include <stdio.h>
+#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<unsigned char>(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<unsigned int>(size_t(resolution) * size_t(resolution));
+		voxel_count = rowpack(grid, voxel_rows, resolution);
+
+		voxels = allocator.allocate<Voxel>(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<unsigned int>(size_t(resolution) * size_t(resolution));
+		unsigned char* queued = allocator.allocate<unsigned char>(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;
+}

Lib/meshoptimizer/src/spatialorder.cpp

@@ -0,0 +1,341 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <float.h>
+#include <string.h>
+
+// This work is based on:
+// Fabian Giesen. Decoding Morton codes. 2009
+namespace meshopt
+{
+
+// "Insert" two 0 bits after each of the 20 low bits of x
+inline unsigned long long part1By2(unsigned long long x)
+{
+	x &= 0x000fffffull;                          // x = ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- jihg fedc ba98 7654 3210
+	x = (x ^ (x << 32)) & 0x000f00000000ffffull; // x = ---- ---- ---- jihg ---- ---- ---- ---- ---- ---- ---- ---- fedc ba98 7654 3210
+	x = (x ^ (x << 16)) & 0x000f0000ff0000ffull; // x = ---- ---- ---- jihg ---- ---- ---- ---- fedc ba98 ---- ---- ---- ---- 7654 3210
+	x = (x ^ (x << 8)) & 0x000f00f00f00f00full;  // x = ---- ---- ---- jihg ---- ---- fedc ---- ---- ba98 ---- ---- 7654 ---- ---- 3210
+	x = (x ^ (x << 4)) & 0x00c30c30c30c30c3ull;  // x = ---- ---- ji-- --hg ---- fe-- --dc ---- ba-- --98 ---- 76-- --54 ---- 32-- --10
+	x = (x ^ (x << 2)) & 0x0249249249249249ull;  // x = ---- --j- -i-- h--g --f- -e-- d--c --b- -a-- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
+	return x;
+}
+
+static void computeOrder(unsigned long long* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, bool morton)
+{
+	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]);
+
+	// rescale each axis to 16 bits to get 48-bit Morton codes
+	float scale = extent == 0 ? 0.f : 65535.f / extent;
+
+	// generate Morton order based on the position inside a unit cube
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		const float* v = vertex_positions_data + i * vertex_stride_float;
+
+		int x = int((v[0] - minv[0]) * scale + 0.5f);
+		int y = int((v[1] - minv[1]) * scale + 0.5f);
+		int z = int((v[2] - minv[2]) * scale + 0.5f);
+
+		if (morton)
+			result[i] = part1By2(x) | (part1By2(y) << 1) | (part1By2(z) << 2);
+		else
+			result[i] = ((unsigned long long)x << 0) | ((unsigned long long)y << 20) | ((unsigned long long)z << 40);
+	}
+}
+
+static void radixSort10(unsigned int* destination, const unsigned int* source, const unsigned short* keys, size_t count)
+{
+	unsigned int hist[1024];
+	memset(hist, 0, sizeof(hist));
+
+	// compute histogram (assume keys are 10-bit)
+	for (size_t i = 0; i < count; ++i)
+		hist[keys[i]]++;
+
+	unsigned int sum = 0;
+
+	// replace histogram data with prefix histogram sums in-place
+	for (int i = 0; i < 1024; ++i)
+	{
+		unsigned int h = hist[i];
+		hist[i] = sum;
+		sum += h;
+	}
+
+	assert(sum == count);
+
+	// reorder values
+	for (size_t i = 0; i < count; ++i)
+	{
+		unsigned int id = keys[source[i]];
+
+		destination[hist[id]++] = source[i];
+	}
+}
+
+static void computeHistogram(unsigned int (&hist)[256][2], const unsigned short* data, size_t count)
+{
+	memset(hist, 0, sizeof(hist));
+
+	// compute 2 8-bit histograms in parallel
+	for (size_t i = 0; i < count; ++i)
+	{
+		unsigned long long id = data[i];
+
+		hist[(id >> 0) & 255][0]++;
+		hist[(id >> 8) & 255][1]++;
+	}
+
+	unsigned int sum0 = 0, sum1 = 0;
+
+	// replace histogram data with prefix histogram sums in-place
+	for (int i = 0; i < 256; ++i)
+	{
+		unsigned int h0 = hist[i][0], h1 = hist[i][1];
+
+		hist[i][0] = sum0;
+		hist[i][1] = sum1;
+
+		sum0 += h0;
+		sum1 += h1;
+	}
+
+	assert(sum0 == count && sum1 == count);
+}
+
+static void radixPass(unsigned int* destination, const unsigned int* source, const unsigned short* keys, size_t count, unsigned int (&hist)[256][2], int pass)
+{
+	int bitoff = pass * 8;
+
+	for (size_t i = 0; i < count; ++i)
+	{
+		unsigned int id = unsigned(keys[source[i]] >> bitoff) & 255;
+
+		destination[hist[id][pass]++] = source[i];
+	}
+}
+
+static void partitionPoints(unsigned int* target, const unsigned int* order, const unsigned char* sides, size_t split, size_t count)
+{
+	size_t l = 0, r = split;
+
+	for (size_t i = 0; i < count; ++i)
+	{
+		unsigned char side = sides[order[i]];
+		target[side ? r : l] = order[i];
+		l += 1;
+		l -= side;
+		r += side;
+	}
+
+	assert(l == split && r == count);
+}
+
+static void splitPoints(unsigned int* destination, unsigned int* orderx, unsigned int* ordery, unsigned int* orderz, const unsigned long long* keys, size_t count, void* scratch, size_t cluster_size)
+{
+	if (count <= cluster_size)
+	{
+		memcpy(destination, orderx, count * sizeof(unsigned int));
+		return;
+	}
+
+	unsigned int* axes[3] = {orderx, ordery, orderz};
+
+	int bestk = -1;
+	unsigned int bestdim = 0;
+
+	for (int k = 0; k < 3; ++k)
+	{
+		const unsigned int mask = (1 << 20) - 1;
+		unsigned int dim = (unsigned(keys[axes[k][count - 1]] >> (k * 20)) & mask) - (unsigned(keys[axes[k][0]] >> (k * 20)) & mask);
+
+		if (dim >= bestdim)
+		{
+			bestk = k;
+			bestdim = dim;
+		}
+	}
+
+	assert(bestk >= 0);
+
+	// split roughly in half, with the left split always being aligned to cluster size
+	size_t split = ((count / 2) + cluster_size - 1) / cluster_size * cluster_size;
+	assert(split > 0 && split < count);
+
+	// mark sides of split for partitioning
+	unsigned char* sides = static_cast<unsigned char*>(scratch) + count * sizeof(unsigned int);
+
+	for (size_t i = 0; i < split; ++i)
+		sides[axes[bestk][i]] = 0;
+
+	for (size_t i = split; i < count; ++i)
+		sides[axes[bestk][i]] = 1;
+
+	// partition all axes into two sides, maintaining order
+	unsigned int* temp = static_cast<unsigned int*>(scratch);
+
+	for (int k = 0; k < 3; ++k)
+	{
+		if (k == bestk)
+			continue;
+
+		unsigned int* axis = axes[k];
+		memcpy(temp, axis, sizeof(unsigned int) * count);
+		partitionPoints(axis, temp, sides, split, count);
+	}
+
+	// recursion depth is logarithmic and bounded as we always split in approximately half
+	splitPoints(destination, orderx, ordery, orderz, keys, split, scratch, cluster_size);
+	splitPoints(destination + split, orderx + split, ordery + split, orderz + split, keys, count - split, scratch, cluster_size);
+}
+
+} // namespace meshopt
+
+void meshopt_spatialSortRemap(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
+{
+	using namespace meshopt;
+
+	assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	meshopt_Allocator allocator;
+
+	unsigned long long* keys = allocator.allocate<unsigned long long>(vertex_count);
+	computeOrder(keys, vertex_positions, vertex_count, vertex_positions_stride, /* morton= */ true);
+
+	unsigned int* scratch = allocator.allocate<unsigned int>(vertex_count * 2); // 4b for order + 2b for keys
+	unsigned short* keyk = (unsigned short*)(scratch + vertex_count);
+
+	for (size_t i = 0; i < vertex_count; ++i)
+		destination[i] = unsigned(i);
+
+	unsigned int* order[] = {scratch, destination};
+
+	// 5-pass radix sort computes the resulting order into scratch
+	for (int k = 0; k < 5; ++k)
+	{
+		// copy 10-bit key segments into keyk to reduce cache pressure during radix pass
+		for (size_t i = 0; i < vertex_count; ++i)
+			keyk[i] = (unsigned short)((keys[i] >> (k * 10)) & 1023);
+
+		radixSort10(order[k % 2], order[(k + 1) % 2], keyk, vertex_count);
+	}
+
+	// since our remap table is mapping old=>new, we need to reverse it
+	for (size_t i = 0; i < vertex_count; ++i)
+		destination[scratch[i]] = unsigned(i);
+}
+
+void meshopt_spatialSortTriangles(unsigned int* destination, 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(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+
+	(void)vertex_count;
+
+	size_t face_count = index_count / 3;
+	size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
+
+	meshopt_Allocator allocator;
+
+	float* centroids = allocator.allocate<float>(face_count * 3);
+
+	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);
+
+		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;
+
+		centroids[i * 3 + 0] = (va[0] + vb[0] + vc[0]) / 3.f;
+		centroids[i * 3 + 1] = (va[1] + vb[1] + vc[1]) / 3.f;
+		centroids[i * 3 + 2] = (va[2] + vb[2] + vc[2]) / 3.f;
+	}
+
+	unsigned int* remap = allocator.allocate<unsigned int>(face_count);
+
+	meshopt_spatialSortRemap(remap, centroids, face_count, sizeof(float) * 3);
+
+	// support in-order remap
+	if (destination == indices)
+	{
+		unsigned int* indices_copy = allocator.allocate<unsigned int>(index_count);
+		memcpy(indices_copy, indices, index_count * sizeof(unsigned int));
+		indices = indices_copy;
+	}
+
+	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];
+		unsigned int r = remap[i];
+
+		destination[r * 3 + 0] = a;
+		destination[r * 3 + 1] = b;
+		destination[r * 3 + 2] = c;
+	}
+}
+
+void meshopt_spatialClusterPoints(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t cluster_size)
+{
+	using namespace meshopt;
+
+	assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
+	assert(vertex_positions_stride % sizeof(float) == 0);
+	assert(cluster_size > 0);
+
+	meshopt_Allocator allocator;
+
+	unsigned long long* keys = allocator.allocate<unsigned long long>(vertex_count);
+	computeOrder(keys, vertex_positions, vertex_count, vertex_positions_stride, /* morton= */ false);
+
+	unsigned int* order = allocator.allocate<unsigned int>(vertex_count * 3);
+	unsigned int* scratch = allocator.allocate<unsigned int>(vertex_count * 2); // 4b for order + 1b for side or 2b for keys
+	unsigned short* keyk = reinterpret_cast<unsigned short*>(scratch + vertex_count);
+
+	for (int k = 0; k < 3; ++k)
+	{
+		// copy 16-bit key segments into keyk to reduce cache pressure during radix pass
+		for (size_t i = 0; i < vertex_count; ++i)
+			keyk[i] = (unsigned short)(keys[i] >> (k * 20));
+
+		unsigned int hist[256][2];
+		computeHistogram(hist, keyk, vertex_count);
+
+		for (size_t i = 0; i < vertex_count; ++i)
+			order[k * vertex_count + i] = unsigned(i);
+
+		radixPass(scratch, order + k * vertex_count, keyk, vertex_count, hist, 0);
+		radixPass(order + k * vertex_count, scratch, keyk, vertex_count, hist, 1);
+	}
+
+	splitPoints(destination, order, order + vertex_count, order + 2 * vertex_count, keys, vertex_count, scratch, cluster_size);
+}

Lib/meshoptimizer/src/stripifier.cpp

@@ -0,0 +1,296 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <limits.h>
+#include <string.h>
+
+// This work is based on:
+// Francine Evans, Steven Skiena and Amitabh Varshney. Optimizing Triangle Strips for Fast Rendering. 1996
+namespace meshopt
+{
+
+static unsigned int findStripFirst(const unsigned int buffer[][3], unsigned int buffer_size, const unsigned char* valence)
+{
+	unsigned int index = 0;
+	unsigned int iv = ~0u;
+
+	for (size_t i = 0; i < buffer_size; ++i)
+	{
+		unsigned char va = valence[buffer[i][0]], vb = valence[buffer[i][1]], vc = valence[buffer[i][2]];
+		unsigned int v = (va < vb && va < vc) ? va : (vb < vc ? vb : vc);
+
+		if (v < iv)
+		{
+			index = unsigned(i);
+			iv = v;
+		}
+	}
+
+	return index;
+}
+
+static int findStripNext(const unsigned int buffer[][3], unsigned int buffer_size, unsigned int e0, unsigned int e1)
+{
+	for (size_t i = 0; i < buffer_size; ++i)
+	{
+		unsigned int a = buffer[i][0], b = buffer[i][1], c = buffer[i][2];
+
+		if (e0 == a && e1 == b)
+			return (int(i) << 2) | 2;
+		else if (e0 == b && e1 == c)
+			return (int(i) << 2) | 0;
+		else if (e0 == c && e1 == a)
+			return (int(i) << 2) | 1;
+	}
+
+	return -1;
+}
+
+} // namespace meshopt
+
+size_t meshopt_stripify(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, unsigned int restart_index)
+{
+	assert(destination != indices);
+	assert(index_count % 3 == 0);
+
+	using namespace meshopt;
+
+	meshopt_Allocator allocator;
+
+	const size_t buffer_capacity = 8;
+
+	unsigned int buffer[buffer_capacity][3] = {};
+	unsigned int buffer_size = 0;
+
+	size_t index_offset = 0;
+
+	unsigned int strip[2] = {};
+	unsigned int parity = 0;
+
+	size_t strip_size = 0;
+
+	// compute vertex valence; this is used to prioritize starting triangle for strips
+	// note: we use 8-bit counters for performance; for outlier vertices the valence is incorrect but that just affects the heuristic
+	unsigned char* valence = allocator.allocate<unsigned char>(vertex_count);
+	memset(valence, 0, vertex_count);
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		valence[index]++;
+	}
+
+	int next = -1;
+
+	while (buffer_size > 0 || index_offset < index_count)
+	{
+		assert(next < 0 || (size_t(next >> 2) < buffer_size && (next & 3) < 3));
+
+		// fill triangle buffer
+		while (buffer_size < buffer_capacity && index_offset < index_count)
+		{
+			buffer[buffer_size][0] = indices[index_offset + 0];
+			buffer[buffer_size][1] = indices[index_offset + 1];
+			buffer[buffer_size][2] = indices[index_offset + 2];
+
+			buffer_size++;
+			index_offset += 3;
+		}
+
+		assert(buffer_size > 0);
+
+		if (next >= 0)
+		{
+			unsigned int i = next >> 2;
+			unsigned int a = buffer[i][0], b = buffer[i][1], c = buffer[i][2];
+			unsigned int v = buffer[i][next & 3];
+
+			// ordered removal from the buffer
+			memmove(buffer[i], buffer[i + 1], (buffer_size - i - 1) * sizeof(buffer[0]));
+			buffer_size--;
+
+			// update vertex valences for strip start heuristic
+			valence[a]--;
+			valence[b]--;
+			valence[c]--;
+
+			// find next triangle (note that edge order flips on every iteration)
+			// in some cases we need to perform a swap to pick a different outgoing triangle edge
+			// for [a b c], the default strip edge is [b c], but we might want to use [a c]
+			int cont = findStripNext(buffer, buffer_size, parity ? strip[1] : v, parity ? v : strip[1]);
+			int swap = cont < 0 ? findStripNext(buffer, buffer_size, parity ? v : strip[0], parity ? strip[0] : v) : -1;
+
+			if (cont < 0 && swap >= 0)
+			{
+				// [a b c] => [a b a c]
+				destination[strip_size++] = strip[0];
+				destination[strip_size++] = v;
+
+				// next strip has same winding
+				// ? a b => b a v
+				strip[1] = v;
+
+				next = swap;
+			}
+			else
+			{
+				// emit the next vertex in the strip
+				destination[strip_size++] = v;
+
+				// next strip has flipped winding
+				strip[0] = strip[1];
+				strip[1] = v;
+				parity ^= 1;
+
+				next = cont;
+			}
+		}
+		else
+		{
+			// if we didn't find anything, we need to find the next new triangle
+			// we use a heuristic to maximize the strip length
+			unsigned int i = findStripFirst(buffer, buffer_size, valence);
+			unsigned int a = buffer[i][0], b = buffer[i][1], c = buffer[i][2];
+
+			// ordered removal from the buffer
+			memmove(buffer[i], buffer[i + 1], (buffer_size - i - 1) * sizeof(buffer[0]));
+			buffer_size--;
+
+			// update vertex valences for strip start heuristic
+			valence[a]--;
+			valence[b]--;
+			valence[c]--;
+
+			// we need to pre-rotate the triangle so that we will find a match in the existing buffer on the next iteration
+			int ea = findStripNext(buffer, buffer_size, c, b);
+			int eb = findStripNext(buffer, buffer_size, a, c);
+			int ec = findStripNext(buffer, buffer_size, b, a);
+
+			// in some cases we can have several matching edges; since we can pick any edge, we pick the one with the smallest
+			// triangle index in the buffer. this reduces the effect of stripification on ACMR and additionally - for unclear
+			// reasons - slightly improves the stripification efficiency
+			int mine = INT_MAX;
+			mine = (ea >= 0 && mine > ea) ? ea : mine;
+			mine = (eb >= 0 && mine > eb) ? eb : mine;
+			mine = (ec >= 0 && mine > ec) ? ec : mine;
+
+			if (ea == mine)
+			{
+				// keep abc
+				next = ea;
+			}
+			else if (eb == mine)
+			{
+				// abc -> bca
+				unsigned int t = a;
+				a = b, b = c, c = t;
+
+				next = eb;
+			}
+			else if (ec == mine)
+			{
+				// abc -> cab
+				unsigned int t = c;
+				c = b, b = a, a = t;
+
+				next = ec;
+			}
+
+			if (restart_index)
+			{
+				if (strip_size)
+					destination[strip_size++] = restart_index;
+
+				destination[strip_size++] = a;
+				destination[strip_size++] = b;
+				destination[strip_size++] = c;
+
+				// new strip always starts with the same edge winding
+				strip[0] = b;
+				strip[1] = c;
+				parity = 1;
+			}
+			else
+			{
+				if (strip_size)
+				{
+					// connect last strip using degenerate triangles
+					destination[strip_size++] = strip[1];
+					destination[strip_size++] = a;
+				}
+
+				// note that we may need to flip the emitted triangle based on parity
+				// we always end up with outgoing edge "cb" in the end
+				unsigned int e0 = parity ? c : b;
+				unsigned int e1 = parity ? b : c;
+
+				destination[strip_size++] = a;
+				destination[strip_size++] = e0;
+				destination[strip_size++] = e1;
+
+				strip[0] = e0;
+				strip[1] = e1;
+				parity ^= 1;
+			}
+		}
+	}
+
+	return strip_size;
+}
+
+size_t meshopt_stripifyBound(size_t index_count)
+{
+	assert(index_count % 3 == 0);
+
+	// worst case without restarts is 2 degenerate indices and 3 indices per triangle
+	// worst case with restarts is 1 restart index and 3 indices per triangle
+	return (index_count / 3) * 5;
+}
+
+size_t meshopt_unstripify(unsigned int* destination, const unsigned int* indices, size_t index_count, unsigned int restart_index)
+{
+	assert(destination != indices);
+
+	size_t offset = 0;
+	size_t start = 0;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		if (restart_index && indices[i] == restart_index)
+		{
+			start = i + 1;
+		}
+		else if (i - start >= 2)
+		{
+			unsigned int a = indices[i - 2], b = indices[i - 1], c = indices[i];
+
+			// flip winding for odd triangles
+			if ((i - start) & 1)
+			{
+				unsigned int t = a;
+				a = b, b = t;
+			}
+
+			// although we use restart indices, strip swaps still produce degenerate triangles, so skip them
+			if (a != b && a != c && b != c)
+			{
+				destination[offset + 0] = a;
+				destination[offset + 1] = b;
+				destination[offset + 2] = c;
+				offset += 3;
+			}
+		}
+	}
+
+	return offset;
+}
+
+size_t meshopt_unstripifyBound(size_t index_count)
+{
+	assert(index_count == 0 || index_count >= 3);
+
+	return (index_count == 0) ? 0 : (index_count - 2) * 3;
+}

Lib/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 <assert.h>
+#include <math.h>
+#include <string.h>
+
+// 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 <bool PositionOnly>
+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<const unsigned int*>(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<const unsigned int*>(vertex_normals + index * vertex_normals_stride_float);
+		const unsigned int* t = reinterpret_cast<const unsigned int*>(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 <typename T, typename Hash>
+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 <bool PositionOnly>
+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<PositionOnly> 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<unsigned int>(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<unsigned int>(vertex_count + 1);
+	adjacency.data = allocator.allocate<unsigned int>(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<unsigned int>(vertex_count);
+	buildVertexRemap<false>(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<Tangent>(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<unsigned int>(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<unsigned int>(face_count);
+	unsigned char* facesign = allocator.allocate<unsigned char>(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<unsigned int>(vertex_count);
+	buildVertexRemap<true>(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<Normal>(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<unsigned int>(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<float>(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);
+}

Lib/meshoptimizer/src/vcacheoptimizer.cpp

@@ -0,0 +1,467 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <string.h>
+
+// This work is based on:
+// Tom Forsyth. Linear-Speed Vertex Cache Optimisation. 2006
+// Pedro Sander, Diego Nehab and Joshua Barczak. Fast Triangle Reordering for Vertex Locality and Reduced Overdraw. 2007
+namespace meshopt
+{
+
+const size_t kCacheSizeMax = 16;
+const size_t kValenceMax = 8;
+
+struct VertexScoreTable
+{
+	float cache[1 + kCacheSizeMax];
+	float live[1 + kValenceMax];
+};
+
+// Tuned to minimize the ACMR of a GPU that has a cache profile similar to NVidia and AMD
+static const VertexScoreTable kVertexScoreTable = {
+    {0.f, 0.779f, 0.791f, 0.789f, 0.981f, 0.843f, 0.726f, 0.847f, 0.882f, 0.867f, 0.799f, 0.642f, 0.613f, 0.600f, 0.568f, 0.372f, 0.234f},
+    {0.f, 0.995f, 0.713f, 0.450f, 0.404f, 0.059f, 0.005f, 0.147f, 0.006f},
+};
+
+// Tuned to minimize the encoded index buffer size
+static const VertexScoreTable kVertexScoreTableStrip = {
+    {0.f, 1.000f, 1.000f, 1.000f, 0.453f, 0.561f, 0.490f, 0.459f, 0.179f, 0.526f, 0.000f, 0.227f, 0.184f, 0.490f, 0.112f, 0.050f, 0.131f},
+    {0.f, 0.956f, 0.786f, 0.577f, 0.558f, 0.618f, 0.549f, 0.499f, 0.489f},
+};
+
+struct TriangleAdjacency
+{
+	unsigned int* counts;
+	unsigned int* offsets;
+	unsigned int* data;
+};
+
+static void buildTriangleAdjacency(TriangleAdjacency& adjacency, const unsigned int* indices, size_t index_count, size_t vertex_count, meshopt_Allocator& allocator)
+{
+	size_t face_count = index_count / 3;
+
+	// allocate arrays
+	adjacency.counts = allocator.allocate<unsigned int>(vertex_count);
+	adjacency.offsets = allocator.allocate<unsigned int>(vertex_count);
+	adjacency.data = allocator.allocate<unsigned int>(index_count);
+
+	// fill triangle counts
+	memset(adjacency.counts, 0, vertex_count * sizeof(unsigned int));
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		assert(indices[i] < vertex_count);
+
+		adjacency.counts[indices[i]]++;
+	}
+
+	// fill offset table
+	unsigned int offset = 0;
+
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		adjacency.offsets[i] = offset;
+		offset += adjacency.counts[i];
+	}
+
+	assert(offset == index_count);
+
+	// fill triangle data
+	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];
+
+		adjacency.data[adjacency.offsets[a]++] = unsigned(i);
+		adjacency.data[adjacency.offsets[b]++] = unsigned(i);
+		adjacency.data[adjacency.offsets[c]++] = unsigned(i);
+	}
+
+	// fix offsets that have been disturbed by the previous pass
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		assert(adjacency.offsets[i] >= adjacency.counts[i]);
+
+		adjacency.offsets[i] -= adjacency.counts[i];
+	}
+}
+
+static unsigned int getNextVertexDeadEnd(const unsigned int* dead_end, unsigned int& dead_end_top, unsigned int& input_cursor, const unsigned int* live_triangles, size_t vertex_count)
+{
+	// check dead-end stack
+	while (dead_end_top)
+	{
+		unsigned int vertex = dead_end[--dead_end_top];
+
+		if (live_triangles[vertex] > 0)
+			return vertex;
+	}
+
+	// input order
+	while (input_cursor < vertex_count)
+	{
+		if (live_triangles[input_cursor] > 0)
+			return input_cursor;
+
+		++input_cursor;
+	}
+
+	return ~0u;
+}
+
+static unsigned int getNextVertexNeighbor(const unsigned int* next_candidates_begin, const unsigned int* next_candidates_end, const unsigned int* live_triangles, const unsigned int* cache_timestamps, unsigned int timestamp, unsigned int cache_size)
+{
+	unsigned int best_candidate = ~0u;
+	int best_priority = -1;
+
+	for (const unsigned int* next_candidate = next_candidates_begin; next_candidate != next_candidates_end; ++next_candidate)
+	{
+		unsigned int vertex = *next_candidate;
+
+		// otherwise we don't need to process it
+		if (live_triangles[vertex] > 0)
+		{
+			int priority = 0;
+
+			// will it be in cache after fanning?
+			if (2 * live_triangles[vertex] + timestamp - cache_timestamps[vertex] <= cache_size)
+			{
+				priority = timestamp - cache_timestamps[vertex]; // position in cache
+			}
+
+			if (priority > best_priority)
+			{
+				best_candidate = vertex;
+				best_priority = priority;
+			}
+		}
+	}
+
+	return best_candidate;
+}
+
+static float vertexScore(const VertexScoreTable* table, int cache_position, unsigned int live_triangles)
+{
+	assert(cache_position >= -1 && cache_position < int(kCacheSizeMax));
+
+	unsigned int live_triangles_clamped = live_triangles < kValenceMax ? live_triangles : kValenceMax;
+
+	return table->cache[1 + cache_position] + table->live[live_triangles_clamped];
+}
+
+static unsigned int getNextTriangleDeadEnd(unsigned int& input_cursor, const unsigned char* emitted_flags, size_t face_count)
+{
+	// input order
+	while (input_cursor < face_count)
+	{
+		if (!emitted_flags[input_cursor])
+			return input_cursor;
+
+		++input_cursor;
+	}
+
+	return ~0u;
+}
+
+} // namespace meshopt
+
+void meshopt_optimizeVertexCacheTable(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, const meshopt::VertexScoreTable* table)
+{
+	using namespace meshopt;
+
+	assert(index_count % 3 == 0);
+
+	meshopt_Allocator allocator;
+
+	// guard for empty meshes
+	if (index_count == 0 || vertex_count == 0)
+		return;
+
+	// support in-place optimization
+	if (destination == indices)
+	{
+		unsigned int* indices_copy = allocator.allocate<unsigned int>(index_count);
+		memcpy(indices_copy, indices, index_count * sizeof(unsigned int));
+		indices = indices_copy;
+	}
+
+	unsigned int cache_size = 16;
+	assert(cache_size <= kCacheSizeMax);
+
+	size_t face_count = index_count / 3;
+
+	// build adjacency information
+	TriangleAdjacency adjacency = {};
+	buildTriangleAdjacency(adjacency, indices, index_count, vertex_count, allocator);
+
+	// live triangle counts; note, we alias adjacency.counts as we remove triangles after emitting them so the counts always match
+	unsigned int* live_triangles = adjacency.counts;
+
+	// emitted flags
+	unsigned char* emitted_flags = allocator.allocate<unsigned char>(face_count);
+	memset(emitted_flags, 0, face_count);
+
+	// compute initial vertex scores
+	float* vertex_scores = allocator.allocate<float>(vertex_count);
+
+	for (size_t i = 0; i < vertex_count; ++i)
+		vertex_scores[i] = vertexScore(table, -1, live_triangles[i]);
+
+	// compute triangle scores
+	float* triangle_scores = allocator.allocate<float>(face_count);
+
+	for (size_t i = 0; i < face_count; ++i)
+	{
+		unsigned int a = indices[i * 3 + 0];
+		unsigned int b = indices[i * 3 + 1];
+		unsigned int c = indices[i * 3 + 2];
+
+		triangle_scores[i] = vertex_scores[a] + vertex_scores[b] + vertex_scores[c];
+	}
+
+	unsigned int cache_holder[2 * (kCacheSizeMax + 4)];
+	unsigned int* cache = cache_holder;
+	unsigned int* cache_new = cache_holder + kCacheSizeMax + 4;
+	size_t cache_count = 0;
+
+	unsigned int current_triangle = 0;
+	unsigned int input_cursor = 1;
+
+	unsigned int output_triangle = 0;
+
+	while (current_triangle != ~0u)
+	{
+		assert(output_triangle < face_count);
+
+		unsigned int a = indices[current_triangle * 3 + 0];
+		unsigned int b = indices[current_triangle * 3 + 1];
+		unsigned int c = indices[current_triangle * 3 + 2];
+
+		// output indices
+		destination[output_triangle * 3 + 0] = a;
+		destination[output_triangle * 3 + 1] = b;
+		destination[output_triangle * 3 + 2] = c;
+		output_triangle++;
+
+		// update emitted flags
+		emitted_flags[current_triangle] = true;
+		triangle_scores[current_triangle] = 0;
+
+		// new triangle
+		size_t cache_write = 0;
+		cache_new[cache_write++] = a;
+		cache_new[cache_write++] = b;
+		cache_new[cache_write++] = c;
+
+		// old triangles
+		for (size_t i = 0; i < cache_count; ++i)
+		{
+			unsigned int index = cache[i];
+
+			cache_new[cache_write] = index;
+			cache_write += (index != a) & (index != b) & (index != c);
+		}
+
+		unsigned int* cache_temp = cache;
+		cache = cache_new, cache_new = cache_temp;
+		cache_count = cache_write > cache_size ? cache_size : cache_write;
+
+		// remove emitted triangle from adjacency data
+		// this makes sure that we spend less time traversing these lists on subsequent iterations
+		// live triangle counts are updated as a byproduct of these adjustments
+		for (size_t k = 0; k < 3; ++k)
+		{
+			unsigned int index = indices[current_triangle * 3 + k];
+
+			unsigned int* neighbors = &adjacency.data[0] + adjacency.offsets[index];
+			size_t neighbors_size = adjacency.counts[index];
+
+			for (size_t i = 0; i < neighbors_size; ++i)
+			{
+				unsigned int tri = neighbors[i];
+
+				if (tri == current_triangle)
+				{
+					neighbors[i] = neighbors[neighbors_size - 1];
+					adjacency.counts[index]--;
+					break;
+				}
+			}
+		}
+
+		unsigned int best_triangle = ~0u;
+		float best_score = 0;
+
+		// update cache positions, vertex scores and triangle scores, and find next best triangle
+		for (size_t i = 0; i < cache_write; ++i)
+		{
+			unsigned int index = cache[i];
+
+			// no need to update scores if we are never going to use this vertex
+			if (adjacency.counts[index] == 0)
+				continue;
+
+			int cache_position = i >= cache_size ? -1 : int(i);
+
+			// update vertex score
+			float score = vertexScore(table, cache_position, live_triangles[index]);
+			float score_diff = score - vertex_scores[index];
+
+			vertex_scores[index] = score;
+
+			// update scores of vertex triangles
+			const unsigned int* neighbors_begin = &adjacency.data[0] + adjacency.offsets[index];
+			const unsigned int* neighbors_end = neighbors_begin + adjacency.counts[index];
+
+			for (const unsigned int* it = neighbors_begin; it != neighbors_end; ++it)
+			{
+				unsigned int tri = *it;
+				assert(!emitted_flags[tri]);
+
+				float tri_score = triangle_scores[tri] + score_diff;
+				assert(tri_score > 0);
+
+				best_triangle = best_score < tri_score ? tri : best_triangle;
+				best_score = best_score < tri_score ? tri_score : best_score;
+
+				triangle_scores[tri] = tri_score;
+			}
+		}
+
+		// step through input triangles in order if we hit a dead-end
+		current_triangle = best_triangle;
+
+		if (current_triangle == ~0u)
+		{
+			current_triangle = getNextTriangleDeadEnd(input_cursor, &emitted_flags[0], face_count);
+		}
+	}
+
+	assert(input_cursor == face_count);
+	assert(output_triangle == face_count);
+}
+
+void meshopt_optimizeVertexCache(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count)
+{
+	meshopt_optimizeVertexCacheTable(destination, indices, index_count, vertex_count, &meshopt::kVertexScoreTable);
+}
+
+void meshopt_optimizeVertexCacheStrip(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count)
+{
+	meshopt_optimizeVertexCacheTable(destination, indices, index_count, vertex_count, &meshopt::kVertexScoreTableStrip);
+}
+
+void meshopt_optimizeVertexCacheFifo(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count, unsigned int cache_size)
+{
+	using namespace meshopt;
+
+	assert(index_count % 3 == 0);
+	assert(cache_size >= 3);
+
+	meshopt_Allocator allocator;
+
+	// guard for empty meshes
+	if (index_count == 0 || vertex_count == 0)
+		return;
+
+	// support in-place optimization
+	if (destination == indices)
+	{
+		unsigned int* indices_copy = allocator.allocate<unsigned int>(index_count);
+		memcpy(indices_copy, indices, index_count * sizeof(unsigned int));
+		indices = indices_copy;
+	}
+
+	size_t face_count = index_count / 3;
+
+	// build adjacency information
+	TriangleAdjacency adjacency = {};
+	buildTriangleAdjacency(adjacency, indices, index_count, vertex_count, allocator);
+
+	// live triangle counts
+	unsigned int* live_triangles = allocator.allocate<unsigned int>(vertex_count);
+	memcpy(live_triangles, adjacency.counts, vertex_count * sizeof(unsigned int));
+
+	// cache time stamps
+	unsigned int* cache_timestamps = allocator.allocate<unsigned int>(vertex_count);
+	memset(cache_timestamps, 0, vertex_count * sizeof(unsigned int));
+
+	// dead-end stack
+	unsigned int* dead_end = allocator.allocate<unsigned int>(index_count);
+	unsigned int dead_end_top = 0;
+
+	// emitted flags
+	unsigned char* emitted_flags = allocator.allocate<unsigned char>(face_count);
+	memset(emitted_flags, 0, face_count);
+
+	unsigned int current_vertex = 0;
+
+	unsigned int timestamp = cache_size + 1;
+	unsigned int input_cursor = 1; // vertex to restart from in case of dead-end
+
+	unsigned int output_triangle = 0;
+
+	while (current_vertex != ~0u)
+	{
+		const unsigned int* next_candidates_begin = &dead_end[0] + dead_end_top;
+
+		// emit all vertex neighbors
+		const unsigned int* neighbors_begin = &adjacency.data[0] + adjacency.offsets[current_vertex];
+		const unsigned int* neighbors_end = neighbors_begin + adjacency.counts[current_vertex];
+
+		for (const unsigned int* it = neighbors_begin; it != neighbors_end; ++it)
+		{
+			unsigned int triangle = *it;
+
+			if (!emitted_flags[triangle])
+			{
+				unsigned int a = indices[triangle * 3 + 0], b = indices[triangle * 3 + 1], c = indices[triangle * 3 + 2];
+
+				// output indices
+				destination[output_triangle * 3 + 0] = a;
+				destination[output_triangle * 3 + 1] = b;
+				destination[output_triangle * 3 + 2] = c;
+				output_triangle++;
+
+				// update dead-end stack
+				dead_end[dead_end_top + 0] = a;
+				dead_end[dead_end_top + 1] = b;
+				dead_end[dead_end_top + 2] = c;
+				dead_end_top += 3;
+
+				// update live triangle counts
+				live_triangles[a]--;
+				live_triangles[b]--;
+				live_triangles[c]--;
+
+				// update cache info
+				// if vertex is not in cache, put it in cache
+				if (timestamp - cache_timestamps[a] > cache_size)
+					cache_timestamps[a] = timestamp++;
+
+				if (timestamp - cache_timestamps[b] > cache_size)
+					cache_timestamps[b] = timestamp++;
+
+				if (timestamp - cache_timestamps[c] > cache_size)
+					cache_timestamps[c] = timestamp++;
+
+				// update emitted flags
+				emitted_flags[triangle] = true;
+			}
+		}
+
+		// next candidates are the ones we pushed to dead-end stack just now
+		const unsigned int* next_candidates_end = &dead_end[0] + dead_end_top;
+
+		// get next vertex
+		current_vertex = getNextVertexNeighbor(next_candidates_begin, next_candidates_end, &live_triangles[0], &cache_timestamps[0], timestamp, cache_size);
+
+		if (current_vertex == ~0u)
+		{
+			current_vertex = getNextVertexDeadEnd(&dead_end[0], dead_end_top, input_cursor, &live_triangles[0], vertex_count);
+		}
+	}
+
+	assert(output_triangle == face_count);
+}

Lib/meshoptimizer/src/vfetchoptimizer.cpp

@@ -0,0 +1,71 @@
+// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
+#include "meshoptimizer.h"
+
+#include <assert.h>
+#include <string.h>
+
+size_t meshopt_optimizeVertexFetchRemap(unsigned int* destination, const unsigned int* indices, size_t index_count, size_t vertex_count)
+{
+	memset(destination, -1, vertex_count * sizeof(unsigned int));
+
+	unsigned int next_vertex = 0;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		if (destination[index] == ~0u)
+		{
+			destination[index] = next_vertex++;
+		}
+	}
+
+	assert(next_vertex <= vertex_count);
+
+	return next_vertex;
+}
+
+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(vertex_size > 0 && vertex_size <= 256);
+
+	meshopt_Allocator allocator;
+
+	// support in-place optimization
+	if (destination == vertices)
+	{
+		unsigned char* vertices_copy = allocator.allocate<unsigned char>(vertex_count * vertex_size);
+		memcpy(vertices_copy, vertices, vertex_count * vertex_size);
+		vertices = vertices_copy;
+	}
+
+	// build vertex remap table
+	unsigned int* vertex_remap = allocator.allocate<unsigned int>(vertex_count);
+	memset(vertex_remap, -1, vertex_count * sizeof(unsigned int));
+
+	unsigned int next_vertex = 0;
+
+	for (size_t i = 0; i < index_count; ++i)
+	{
+		unsigned int index = indices[i];
+		assert(index < vertex_count);
+
+		unsigned int& remap = vertex_remap[index];
+
+		if (remap == ~0u) // vertex was not added to destination VB
+		{
+			// add vertex
+			memcpy(static_cast<unsigned char*>(destination) + next_vertex * vertex_size, static_cast<const unsigned char*>(vertices) + index * vertex_size, vertex_size);
+
+			remap = next_vertex++;
+		}
+
+		// modify indices in place
+		indices[i] = remap;
+	}
+
+	assert(next_vertex <= vertex_count);
+
+	return next_vertex;
+}