text: // 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
  
  // This work is based on:
  // Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997
  // Michael Garland. Quadric-based polygonal surface simplification. 1999
  // Peter Lindstrom. Out-of-Core Simplification of Large Polygonal Models. 2000
  // Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003
  // Peter Van Sandt, Yannis Chronis, Jignesh M. Patel. Efficiently Searching In-Memory Sorted Arrays: Revenge of the Interpolation Search? 2019
  namespace meshopt
  {
  
  typedef unsigned long long VertexID;
  
  const unsigned long long POS_GRID_SIZE_BITS = 15;
  const unsigned long long UV_ID_BITS = 13;
  const unsigned long long NRM_GRID_SIZE_BITS = 2;
  static_assert(POS_GRID_SIZE_BITS * 3 + UV_ID_BITS + NRM_GRID_SIZE_BITS * 3 <= 64, "VertexID overflow");
  
  struct EdgeAdjacency
  {
  	unsigned int* counts;
  	unsigned int* offsets;
  	unsigned int* data;
  };
  
  static void buildEdgeAdjacency(EdgeAdjacency& 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 edge 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 edge 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]++] = b;
  		adjacency.data[adjacency.offsets[b]++] = c;
  		adjacency.data[adjacency.offsets[c]++] = a;
  	}
  
  	// 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];
  	}
  }
  
  struct PositionHasher
  {
  	const float* vertex_positions;
  	size_t vertex_stride_float;
  
  	size_t hash(unsigned int index) const
  	{
  		// MurmurHash2
  		const unsigned int m = 0x5bd1e995;
  		const int r = 24;
  
  		unsigned int h = 0;
  		const unsigned int* key = reinterpret_cast<const unsigned int*>(vertex_positions + index * vertex_stride_float);
  
  		for (size_t i = 0; i < 3; ++i)
  		{
  			unsigned int k = key[i];
  
  			k *= m;
  			k ^= k >> r;
  			k *= m;
  
  			h *= m;
  			h ^= k;
  		}
  
  		return h;
  	}
  
  	bool equal(unsigned int lhs, unsigned int rhs) const
  	{
  		return memcmp(vertex_positions + lhs * vertex_stride_float, vertex_positions + rhs * vertex_stride_float, sizeof(float) * 3) == 0;
  	}
  };
  
  static size_t hashBuckets2(size_t count)
  {
  	size_t buckets = 1;
  	while (buckets < count)
  		buckets *= 2;
  
  	return buckets;
  }
  
  template <typename T, typename Hash>
  static T* hashLookup2(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 0;
  }
  
  static void buildPositionRemap(unsigned int* remap, unsigned int* wedge, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, meshopt_Allocator& allocator)
  {
  	PositionHasher hasher = {vertex_positions_data, vertex_positions_stride / sizeof(float)};
  
  	size_t table_size = hashBuckets2(vertex_count);
  	unsigned int* table = allocator.allocate<unsigned int>(table_size);
  	memset(table, -1, table_size * sizeof(unsigned int));
  
  	// build forward remap: for each vertex, which other (canonical) vertex does it map to?
  	// we use position equivalence for this, and remap vertices to other existing vertices
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		unsigned int index = unsigned(i);
  		unsigned int* entry = hashLookup2(table, table_size, hasher, index, ~0u);
  
  		if (*entry == ~0u)
  			*entry = index;
  
  		remap[index] = *entry;
  	}
  
  	// build wedge table: for each vertex, which other vertex is the next wedge that also maps to the same vertex?
  	// entries in table form a (cyclic) wedge loop per vertex; for manifold vertices, wedge[i] == remap[i] == i
  	for (size_t i = 0; i < vertex_count; ++i)
  		wedge[i] = unsigned(i);
  
  	for (size_t i = 0; i < vertex_count; ++i)
  		if (remap[i] != i)
  		{
  			unsigned int r = remap[i];
  
  			wedge[i] = wedge[r];
  			wedge[r] = unsigned(i);
  		}
  }
  
  enum VertexKind
  {
  	Kind_Manifold, // not on an attribute seam, not on any boundary
  	Kind_Border,   // not on an attribute seam, has exactly two open edges
  	Kind_Seam,     // on an attribute seam with exactly two attribute seam edges
  	Kind_Complex,  // none of the above; these vertices can move as long as all wedges move to the target vertex
  	Kind_Locked,   // none of the above; these vertices can't move
  
  	Kind_Count
  };
  
  // manifold vertices can collapse onto anything
  // border/seam vertices can only be collapsed onto border/seam respectively
  // complex vertices can collapse onto complex/locked
  // a rule of thumb is that collapsing kind A into kind B preserves the kind B in the target vertex
  // for example, while we could collapse Complex into Manifold, this would mean the target vertex isn't Manifold anymore
  const unsigned char kCanCollapse[Kind_Count][Kind_Count] = {
      {1, 1, 1, 1, 1},
      {0, 1, 0, 0, 0},
      {0, 0, 1, 0, 0},
      {0, 0, 0, 1, 1},
      {0, 0, 0, 0, 0},
  };
  
  // if a vertex is manifold or seam, adjoining edges are guaranteed to have an opposite edge
  // note that for seam edges, the opposite edge isn't present in the attribute-based topology
  // but is present if you consider a position-only mesh variant
  const unsigned char kHasOpposite[Kind_Count][Kind_Count] = {
      {1, 1, 1, 0, 1},
      {1, 0, 1, 0, 0},
      {1, 1, 1, 0, 1},
      {0, 0, 0, 0, 0},
      {1, 0, 1, 0, 0},
  };
  
  static bool hasEdge(const EdgeAdjacency& adjacency, unsigned int a, unsigned int b)
  {
  	unsigned int count = adjacency.counts[a];
  	const unsigned int* data = adjacency.data + adjacency.offsets[a];
  
  	for (size_t i = 0; i < count; ++i)
  		if (data[i] == b)
  			return true;
  
  	return false;
  }
  
  static unsigned int findWedgeEdge(const EdgeAdjacency& adjacency, const unsigned int* wedge, unsigned int a, unsigned int b)
  {
  	unsigned int v = a;
  
  	do
  	{
  		if (hasEdge(adjacency, v, b))
  			return v;
  
  		v = wedge[v];
  	} while (v != a);
  
  	return ~0u;
  }
  
  static size_t countOpenEdges(const EdgeAdjacency& adjacency, unsigned int vertex, unsigned int* last = 0)
  {
  	size_t result = 0;
  
  	unsigned int count = adjacency.counts[vertex];
  	const unsigned int* data = adjacency.data + adjacency.offsets[vertex];
  
  	for (size_t i = 0; i < count; ++i)
  		if (!hasEdge(adjacency, data[i], vertex))
  		{
  			result++;
  
  			if (last)
  				*last = data[i];
  		}
  
  	return result;
  }
  
  static void classifyVertices(unsigned char* result, unsigned int* loop, size_t vertex_count, const EdgeAdjacency& adjacency, const unsigned int* remap, const unsigned int* wedge)
  {
  	for (size_t i = 0; i < vertex_count; ++i)
  		loop[i] = ~0u;
  
  #if TRACE
  	size_t lockedstats[4] = {};
  #define TRACELOCKED(i) lockedstats[i]++;
  #else
  #define TRACELOCKED(i) (void)0
  #endif
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		if (remap[i] == i)
  		{
  			if (wedge[i] == i)
  			{
  				// no attribute seam, need to check if it's manifold
  				unsigned int v = 0;
  				size_t edges = countOpenEdges(adjacency, unsigned(i), &v);
  
  				// note: we classify any vertices with no open edges as manifold
  				// this is technically incorrect - if 4 triangles share an edge, we'll classify vertices as manifold
  				// it's unclear if this is a problem in practice
  				// also note that we classify vertices as border if they have *one* open edge, not two
  				// this is because we only have half-edges - so a border vertex would have one incoming and one outgoing edge
  				if (edges == 0)
  				{
  					result[i] = Kind_Manifold;
  				}
  				else if (edges == 1)
  				{
  					result[i] = Kind_Border;
  					loop[i] = v;
  				}
  				else
  				{
  					result[i] = Kind_Locked;
  					TRACELOCKED(0);
  				}
  			}
  			else if (wedge[wedge[i]] == i)
  			{
  				// attribute seam; need to distinguish between Seam and Locked
  				unsigned int a = 0;
  				size_t a_count = countOpenEdges(adjacency, unsigned(i), &a);
  				unsigned int b = 0;
  				size_t b_count = countOpenEdges(adjacency, wedge[i], &b);
  
  				// seam should have one open half-edge for each vertex, and the edges need to "connect" - point to the same vertex post-remap
  				if (a_count == 1 && b_count == 1)
  				{
  					unsigned int ao = findWedgeEdge(adjacency, wedge, a, wedge[i]);
  					unsigned int bo = findWedgeEdge(adjacency, wedge, b, unsigned(i));
  
  					if (ao != ~0u && bo != ~0u)
  					{
  						result[i] = Kind_Seam;
  
  						loop[i] = a;
  						loop[wedge[i]] = b;
  					}
  					else
  					{
  						result[i] = Kind_Locked;
  						TRACELOCKED(1);
  					}
  				}
  				else
  				{
  					result[i] = Kind_Locked;
  					TRACELOCKED(2);
  				}
  			}
  			else
  			{
  				// more than one vertex maps to this one; we don't have classification available
  				result[i] = Kind_Locked;
  				TRACELOCKED(3);
  			}
  		}
  		else
  		{
  			assert(remap[i] < i);
  
  			result[i] = result[remap[i]];
  		}
  	}
  
  #if TRACE
  	printf("locked: many open edges %d, disconnected seam %d, many seam edges %d, many wedges %d\n",
  	       int(lockedstats[0]), int(lockedstats[1]), int(lockedstats[2]), int(lockedstats[3]));
  #endif
  }
  
  struct Vector3
  {
  	float x, y, z;
  };
  
  static void rescalePositions(Vector3* result, const float* vertex_positions_data, 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_data + i * vertex_stride_float;
  
  		result[i].x = v[0];
  		result[i].y = v[1];
  		result[i].z = v[2];
  
  		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 ? 0.f : 1.f / extent;
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		result[i].x = (result[i].x - minv[0]) * scale;
  		result[i].y = (result[i].y - minv[1]) * scale;
  		result[i].z = (result[i].z - minv[2]) * scale;
  	}
  }
  
  struct Quadric
  {
  	float a00, a11, a22;
  	float a10, a20, a21;
  	float b0, b1, b2, c;
  	float w;
  };
  
  struct Collapse
  {
  	unsigned int v0;
  	unsigned int v1;
  
  	union {
  		unsigned int bidi;
  		float error;
  		unsigned int errorui;
  	};
  };
  
  static float normalize(Vector3& v)
  {
  	float length = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
  
  	if (length > 0)
  	{
  		v.x /= length;
  		v.y /= length;
  		v.z /= length;
  	}
  
  	return length;
  }
  
  static void quadricAdd(Quadric& Q, const Quadric& R)
  {
  	Q.a00 += R.a00;
  	Q.a11 += R.a11;
  	Q.a22 += R.a22;
  	Q.a10 += R.a10;
  	Q.a20 += R.a20;
  	Q.a21 += R.a21;
  	Q.b0 += R.b0;
  	Q.b1 += R.b1;
  	Q.b2 += R.b2;
  	Q.c += R.c;
  	Q.w += R.w;
  }
  
  static float quadricError(const Quadric& Q, const Vector3& v)
  {
  	float rx = Q.b0;
  	float ry = Q.b1;
  	float rz = Q.b2;
  
  	rx += Q.a10 * v.y;
  	ry += Q.a21 * v.z;
  	rz += Q.a20 * v.x;
  
  	rx *= 2;
  	ry *= 2;
  	rz *= 2;
  
  	rx += Q.a00 * v.x;
  	ry += Q.a11 * v.y;
  	rz += Q.a22 * v.z;
  
  	float r = Q.c;
  	r += rx * v.x;
  	r += ry * v.y;
  	r += rz * v.z;
  
  	float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
  
  	return fabsf(r) * s;
  }
  
  static void quadricFromPlane(Quadric& Q, 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;
  
  	Q.a00 = a * aw;
  	Q.a11 = b * bw;
  	Q.a22 = c * cw;
  	Q.a10 = a * bw;
  	Q.a20 = a * cw;
  	Q.a21 = b * cw;
  	Q.b0 = a * dw;
  	Q.b1 = b * dw;
  	Q.b2 = c * dw;
  	Q.c = d * dw;
  	Q.w = w;
  }
  
  static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
  {
  	// we need to encode (x - X) ^ 2 + (y - Y)^2 + (z - Z)^2 into the quadric
  	Q.a00 = w;
  	Q.a11 = w;
  	Q.a22 = w;
  	Q.a10 = 0.f;
  	Q.a20 = 0.f;
  	Q.a21 = 0.f;
  	Q.b0 = -2.f * x * w;
  	Q.b1 = -2.f * y * w;
  	Q.b2 = -2.f * z * w;
  	Q.c = (x * x + y * y + z * z) * w;
  	Q.w = w;
  }
  
  static void quadricFromTriangle(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
  {
  	Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
  	Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
  
  	// normal = cross(p1 - p0, p2 - p0)
  	Vector3 normal = {p10.y * p20.z - p10.z * p20.y, p10.z * p20.x - p10.x * p20.z, p10.x * p20.y - p10.y * p20.x};
  	float area = normalize(normal);
  
  	float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
  
  	// we use sqrtf(area) so that the error is scaled linearly; this tends to improve silhouettes
  	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, sqrtf(area) * weight);
  }
  
  static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
  {
  	Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
  	float length = normalize(p10);
  
  	// p20p = length of projection of p2-p0 onto normalize(p1 - p0)
  	Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
  	float p20p = p20.x * p10.x + p20.y * p10.y + p20.z * p10.z;
  
  	// normal = altitude of triangle from point p2 onto edge p1-p0
  	Vector3 normal = {p20.x - p10.x * p20p, p20.y - p10.y * p20p, p20.z - p10.z * p20p};
  	normalize(normal);
  
  	float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
  
  	// note: the weight is scaled linearly with edge length; this has to match the triangle weight
  	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
  }
  
  static void fillFaceQuadrics(Quadric* vertex_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap)
  {
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		unsigned int i0 = indices[i + 0];
  		unsigned int i1 = indices[i + 1];
  		unsigned int i2 = indices[i + 2];
  
  		Quadric Q;
  		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
  
  		quadricAdd(vertex_quadrics[remap[i0]], Q);
  		quadricAdd(vertex_quadrics[remap[i1]], Q);
  		quadricAdd(vertex_quadrics[remap[i2]], Q);
  	}
  }
  
  static void fillEdgeQuadrics(Quadric* vertex_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap, const unsigned char* vertex_kind, const unsigned int* loop)
  {
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		static const int next[3] = {1, 2, 0};
  
  		for (int e = 0; e < 3; ++e)
  		{
  			unsigned int i0 = indices[i + e];
  			unsigned int i1 = indices[i + next[e]];
  
  			unsigned char k0 = vertex_kind[i0];
  			unsigned char k1 = vertex_kind[i1];
  
  			// check that i0 and i1 are border/seam and are on the same edge loop
  			// loop[] tracks half edges so we only need to check i0->i1
  			if (k0 != k1 || (k0 != Kind_Border && k0 != Kind_Seam) || loop[i0] != i1)
  				continue;
  
  			unsigned int i2 = indices[i + next[next[e]]];
  
  			// we try hard to maintain border edge geometry; seam edges can move more freely
  			// due to topological restrictions on collapses, seam quadrics slightly improves collapse structure but aren't critical
  			const float kEdgeWeightSeam = 1.f;
  			const float kEdgeWeightBorder = 10.f;
  
  			float edgeWeight = (k0 == Kind_Seam) ? kEdgeWeightSeam : kEdgeWeightBorder;
  
  			Quadric Q;
  			quadricFromTriangleEdge(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], edgeWeight);
  
  			quadricAdd(vertex_quadrics[remap[i0]], Q);
  			quadricAdd(vertex_quadrics[remap[i1]], Q);
  		}
  	}
  }
  
  static size_t pickEdgeCollapses(Collapse* collapses, const unsigned int* indices, size_t index_count, const unsigned int* remap, const unsigned char* vertex_kind, const unsigned int* loop)
  {
  	size_t collapse_count = 0;
  
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		static const int next[3] = {1, 2, 0};
  
  		for (int e = 0; e < 3; ++e)
  		{
  			unsigned int i0 = indices[i + e];
  			unsigned int i1 = indices[i + next[e]];
  
  			// this can happen either when input has a zero-length edge, or when we perform collapses for complex
  			// topology w/seams and collapse a manifold vertex that connects to both wedges onto one of them
  			// we leave edges like this alone since they may be important for preserving mesh integrity
  			if (remap[i0] == remap[i1])
  				continue;
  
  			unsigned char k0 = vertex_kind[i0];
  			unsigned char k1 = vertex_kind[i1];
  
  			// the edge has to be collapsible in at least one direction
  			if (!(kCanCollapse[k0][k1] | kCanCollapse[k1][k0]))
  				continue;
  
  			// manifold and seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
  			if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
  				continue;
  
  			// two vertices are on a border or a seam, but there's no direct edge between them
  			// this indicates that they belong to two different edge loops and we should not collapse this edge
  			// loop[] tracks half edges so we only need to check i0->i1
  			if (k0 == k1 && (k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
  				continue;
  
  			// edge can be collapsed in either direction - we will pick the one with minimum error
  			// note: we evaluate error later during collapse ranking, here we just tag the edge as bidirectional
  			if (kCanCollapse[k0][k1] & kCanCollapse[k1][k0])
  			{
  				Collapse c = {i0, i1, {/* bidi= */ 1}};
  				collapses[collapse_count++] = c;
  			}
  			else
  			{
  				// edge can only be collapsed in one direction
  				unsigned int e0 = kCanCollapse[k0][k1] ? i0 : i1;
  				unsigned int e1 = kCanCollapse[k0][k1] ? i1 : i0;
  
  				Collapse c = {e0, e1, {/* bidi= */ 0}};
  				collapses[collapse_count++] = c;
  			}
  		}
  	}
  
  	return collapse_count;
  }
  
  static void rankEdgeCollapses(Collapse* collapses, size_t collapse_count, const Vector3* vertex_positions, const Quadric* vertex_quadrics, const unsigned int* remap)
  {
  	for (size_t i = 0; i < collapse_count; ++i)
  	{
  		Collapse& c = collapses[i];
  
  		unsigned int i0 = c.v0;
  		unsigned int i1 = c.v1;
  
  		// most edges are bidirectional which means we need to evaluate errors for two collapses
  		// to keep this code branchless we just use the same edge for unidirectional edges
  		unsigned int j0 = c.bidi ? i1 : i0;
  		unsigned int j1 = c.bidi ? i0 : i1;
  
  		const Quadric& qi = vertex_quadrics[remap[i0]];
  		const Quadric& qj = vertex_quadrics[remap[j0]];
  
  		float ei = quadricError(qi, vertex_positions[i1]);
  		float ej = quadricError(qj, vertex_positions[j1]);
  
  		// pick edge direction with minimal error
  		c.v0 = ei <= ej ? i0 : j0;
  		c.v1 = ei <= ej ? i1 : j1;
  		c.error = ei <= ej ? ei : ej;
  	}
  }
  
  #if TRACE > 1
  static void dumpEdgeCollapses(const Collapse* collapses, size_t collapse_count, const unsigned char* vertex_kind)
  {
  	size_t ckinds[Kind_Count][Kind_Count] = {};
  	float cerrors[Kind_Count][Kind_Count] = {};
  
  	for (int k0 = 0; k0 < Kind_Count; ++k0)
  		for (int k1 = 0; k1 < Kind_Count; ++k1)
  			cerrors[k0][k1] = FLT_MAX;
  
  	for (size_t i = 0; i < collapse_count; ++i)
  	{
  		unsigned int i0 = collapses[i].v0;
  		unsigned int i1 = collapses[i].v1;
  
  		unsigned char k0 = vertex_kind[i0];
  		unsigned char k1 = vertex_kind[i1];
  
  		ckinds[k0][k1]++;
  		cerrors[k0][k1] = (collapses[i].error < cerrors[k0][k1]) ? collapses[i].error : cerrors[k0][k1];
  	}
  
  	for (int k0 = 0; k0 < Kind_Count; ++k0)
  		for (int k1 = 0; k1 < Kind_Count; ++k1)
  			if (ckinds[k0][k1])
  				printf("collapses %d -> %d: %d, min error %e\n", k0, k1, int(ckinds[k0][k1]), cerrors[k0][k1]);
  }
  
  static void dumpLockedCollapses(const unsigned int* indices, size_t index_count, const unsigned char* vertex_kind)
  {
  	size_t locked_collapses[Kind_Count][Kind_Count] = {};
  
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		static const int next[3] = {1, 2, 0};
  
  		for (int e = 0; e < 3; ++e)
  		{
  			unsigned int i0 = indices[i + e];
  			unsigned int i1 = indices[i + next[e]];
  
  			unsigned char k0 = vertex_kind[i0];
  			unsigned char k1 = vertex_kind[i1];
  
  			locked_collapses[k0][k1] += !kCanCollapse[k0][k1] && !kCanCollapse[k1][k0];
  		}
  	}
  
  	for (int k0 = 0; k0 < Kind_Count; ++k0)
  		for (int k1 = 0; k1 < Kind_Count; ++k1)
  			if (locked_collapses[k0][k1])
  				printf("locked collapses %d -> %d: %d\n", k0, k1, int(locked_collapses[k0][k1]));
  }
  #endif
  
  static void sortEdgeCollapses(unsigned int* sort_order, const Collapse* collapses, size_t collapse_count)
  {
  	const int sort_bits = 11;
  
  	// fill histogram for counting sort
  	unsigned int histogram[1 << sort_bits];
  	memset(histogram, 0, sizeof(histogram));
  
  	for (size_t i = 0; i < collapse_count; ++i)
  	{
  		// skip sign bit since error is non-negative
  		unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
  
  		histogram[key]++;
  	}
  
  	// 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 == collapse_count);
  
  	// compute sort order based on offsets
  	for (size_t i = 0; i < collapse_count; ++i)
  	{
  		// skip sign bit since error is non-negative
  		unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
  
  		sort_order[histogram[key]++] = unsigned(i);
  	}
  }
  
  static size_t performEdgeCollapses(unsigned int* collapse_remap, unsigned char* collapse_locked, Quadric* vertex_quadrics, const Collapse* collapses, size_t collapse_count, const unsigned int* collapse_order, const unsigned int* remap, const unsigned int* wedge, const unsigned char* vertex_kind, size_t triangle_collapse_goal, float error_goal, float error_limit)
  {
  	size_t edge_collapses = 0;
  	size_t triangle_collapses = 0;
  
  	for (size_t i = 0; i < collapse_count; ++i)
  	{
  		const Collapse& c = collapses[collapse_order[i]];
  
  		if (c.error > error_limit)
  			break;
  
  		if (c.error > error_goal && triangle_collapses > triangle_collapse_goal / 10)
  			break;
  
  		if (triangle_collapses >= triangle_collapse_goal)
  			break;
  
  		unsigned int i0 = c.v0;
  		unsigned int i1 = c.v1;
  
  		unsigned int r0 = remap[i0];
  		unsigned int r1 = remap[i1];
  
  		// we don't collapse vertices that had source or target vertex involved in a collapse
  		// it's important to not move the vertices twice since it complicates the tracking/remapping logic
  		// it's important to not move other vertices towards a moved vertex to preserve error since we don't re-rank collapses mid-pass
  		if (collapse_locked[r0] | collapse_locked[r1])
  			continue;
  
  		assert(collapse_remap[r0] == r0);
  		assert(collapse_remap[r1] == r1);
  
  		quadricAdd(vertex_quadrics[r1], vertex_quadrics[r0]);
  
  		if (vertex_kind[i0] == Kind_Complex)
  		{
  			unsigned int v = i0;
  
  			do
  			{
  				collapse_remap[v] = r1;
  				v = wedge[v];
  			} while (v != i0);
  		}
  		else if (vertex_kind[i0] == Kind_Seam)
  		{
  			// remap v0 to v1 and seam pair of v0 to seam pair of v1
  			unsigned int s0 = wedge[i0];
  			unsigned int s1 = wedge[i1];
  
  			assert(s0 != i0 && s1 != i1);
  			assert(wedge[s0] == i0 && wedge[s1] == i1);
  
  			collapse_remap[i0] = i1;
  			collapse_remap[s0] = s1;
  		}
  		else
  		{
  			assert(wedge[i0] == i0);
  
  			collapse_remap[i0] = i1;
  		}
  
  		collapse_locked[r0] = 1;
  		collapse_locked[r1] = 1;
  
  		// border edges collapse 1 triangle, other edges collapse 2 or more
  		triangle_collapses += (vertex_kind[i0] == Kind_Border) ? 1 : 2;
  		edge_collapses++;
  	}
  
  	return edge_collapses;
  }
  
  static size_t remapIndexBuffer(unsigned int* indices, size_t index_count, const unsigned int* collapse_remap)
  {
  	size_t write = 0;
  
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		unsigned int v0 = collapse_remap[indices[i + 0]];
  		unsigned int v1 = collapse_remap[indices[i + 1]];
  		unsigned int v2 = collapse_remap[indices[i + 2]];
  
  		// we never move the vertex twice during a single pass
  		assert(collapse_remap[v0] == v0);
  		assert(collapse_remap[v1] == v1);
  		assert(collapse_remap[v2] == v2);
  
  		if (v0 != v1 && v0 != v2 && v1 != v2)
  		{
  			indices[write + 0] = v0;
  			indices[write + 1] = v1;
  			indices[write + 2] = v2;
  			write += 3;
  		}
  	}
  
  	return write;
  }
  
  static void remapEdgeLoops(unsigned int* loop, size_t vertex_count, const unsigned int* collapse_remap)
  {
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		if (loop[i] != ~0u)
  		{
  			unsigned int l = loop[i];
  			unsigned int r = collapse_remap[l];
  
  			// i == r is a special case when the seam edge is collapsed in a direction opposite to where loop goes
  			loop[i] = (i == r) ? loop[l] : r;
  		}
  	}
  }
  
  struct CellHasher
  {
  	const VertexID* vertex_ids;
  
  	size_t hash(unsigned int i) const
  	{
  		VertexID h = vertex_ids[i];
  
  		// MurmurHash2 finalizer
  		h ^= h >> 47;
  		h *= 0xc6a4a7935bd1e995;
  		h ^= h >> 47;
  		return h;
  	}
  
  	bool equal(unsigned int lhs, unsigned int rhs) const
  	{
  		return vertex_ids[lhs] == vertex_ids[rhs];
  	}
  };
  
  struct IdHasher
  {
  	size_t hash(VertexID id) const
  	{
  		VertexID h = id;
  
  		// MurmurHash2 finalizer
  		h ^= h >> 47;
  		h *= 0xc6a4a7935bd1e995;
  		h ^= h >> 47;
  		return h;
  	}
  
  	bool equal(VertexID lhs, VertexID rhs) const
  	{
  		return lhs == rhs;
  	}
  };
  
  struct TriangleHasher
  {
  	unsigned int* indices;
  
  	size_t hash(unsigned int i) const
  	{
  		const unsigned int* tri = indices + i * 3;
  
  		// Optimized Spatial Hashing for Collision Detection of Deformable Objects
  		return (tri[0] * 73856093) ^ (tri[1] * 19349663) ^ (tri[2] * 83492791);
  	}
  
  	bool equal(unsigned int lhs, unsigned int rhs) const
  	{
  		const unsigned int* lt = indices + lhs * 3;
  		const unsigned int* rt = indices + rhs * 3;
  
  		return lt[0] == rt[0] && lt[1] == rt[1] && lt[2] == rt[2];
  	}
  };
  
  struct Triangle
  {
  	VertexID a;
  	VertexID b;
  	VertexID c;
  
  	bool operator==(const Triangle& rhs) const
  	{
  		return a == rhs.a && b == rhs.b && c == rhs.c;
  	}
  };
  
  struct MemorylessTriangleHasher
  {
  	size_t hash(const Triangle& tri) const
  	{
  		// Optimized Spatial Hashing for Collision Detection of Deformable Objects (with some random 64-bit primes)
  		return (tri.a * 14372678844667432909ull) ^ (tri.b * 3832277515459471913ull) ^ (tri.c * 10359333880808892121ull);
  	}
  
  	bool equal(const Triangle& lhs, const Triangle& rhs) const
  	{
  		return lhs == rhs;
  	}
  };
  
  static void computeVertexIds(VertexID* vertex_ids, const Vector3* vertex_positions, const float* vertex_normals, const unsigned int* uv_islands, size_t vertex_normals_stride, size_t vertex_count, int grid_size)
  {
  	assert(grid_size >= 1 && grid_size <= (1 << POS_GRID_SIZE_BITS));
  	float cell_scale = float(grid_size - 1);
  	size_t normals_stride_floats = vertex_normals_stride / sizeof(float);
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		const Vector3& pos = vertex_positions[i];
  
  		VertexID xi = VertexID(pos.x * cell_scale + 0.5f);
  		VertexID yi = VertexID(pos.y * cell_scale + 0.5f);
  		VertexID zi = VertexID(pos.z * cell_scale + 0.5f);
  
  		VertexID xni = 0;
  		VertexID yni = 0;
  		VertexID zni = 0;
  		if (vertex_normals)
  		{
  			// Split [-1, +1] into 3 regions, which seems to be enough to prevent issues.
  			// It's important that there isn't a region boundary at 0, because this would prevent
  			// nearly-parallel normals from being merged in more cases (assuming a roughly axis-aligned mesh).
  			const float* np = &vertex_normals[i * normals_stride_floats];
  			Vector3 n = {np[0], np[1], np[2]};
  			assert(fabsf(n.x * n.x + n.y * n.y + n.z * n.z - 1.f) < 1e-3f);
  			xni = VertexID((n.x + 1.f) * 1.5f - 1e-5f);
  			yni = VertexID((n.y + 1.f) * 1.5f - 1e-5f);
  			zni = VertexID((n.z + 1.f) * 1.5f - 1e-5f);
  		}
  
  		VertexID uv_island = 0;
  		if (uv_islands)
  		{
  			const VertexID max_island = (1 << UV_ID_BITS);
  			assert(uv_islands[i] < max_island);
  			uv_island = (uv_islands[i] < max_island) ? uv_islands[i] : 0;
  		}
  		
  		vertex_ids[i] = (xni << 62ull) | (yni << 60ull) | (zni << 58ull) 
  			| (uv_island << (3 * POS_GRID_SIZE_BITS)) 
  			| (xi << (2 * POS_GRID_SIZE_BITS)) | (yi << POS_GRID_SIZE_BITS) | zi;
  	}
  }
  
  static size_t countTriangles(Triangle* tritable, size_t tritable_size, const VertexID* vertex_ids, const VertexID* vertex_ids_pos_only, const unsigned int* indices, size_t index_count, const bool* locked_vertices)
  {
  	size_t result = 0;
  
  	memset(tritable, -1, tritable_size * sizeof(Triangle));
  	const Triangle invalid_triangle = {~0ull, ~0ull, ~0ull};
  
  	MemorylessTriangleHasher hasher;
  
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		// Use IDs based on position only to check for degeneracy.
  		VertexID a2 = vertex_ids_pos_only[indices[i + 0]];
  		VertexID b2 = vertex_ids_pos_only[indices[i + 1]];
  		VertexID c2 = vertex_ids_pos_only[indices[i + 2]];
  
  		// If one of the vertices in this triangle is locked, add 1 and continue to the next.
  		if (locked_vertices)
  		{
  			if (locked_vertices[indices[i + 0]] || locked_vertices[indices[i + 1]] || locked_vertices[indices[i + 2]])
  			{
  				++result;
  			}
  		}
  
  		if (a2 != b2 && a2 != c2 && b2 != c2)
  		{
  			VertexID a = vertex_ids[indices[i + 0]];
  			VertexID b = vertex_ids[indices[i + 1]];
  			VertexID c = vertex_ids[indices[i + 2]];
  			assert(a != b && a != c && b != c);
  
  			// Sort IDs to ensure equivalent triangles compare equal.
  			if (b < a && b < c)
  			{
  				VertexID t = a;
  				a = b, b = c, c = t;
  			}
  			else if (c < a && c < b)
  			{
  				VertexID t = c;
  				c = b, b = a, a = t;
  			}
  
  			// Check if triangle has already been seen.
  			Triangle tri = {a, b, c};
  			Triangle* entry = hashLookup2(tritable, tritable_size, hasher, tri, invalid_triangle);
  
  			if (*entry == invalid_triangle)
  			{
  				*entry = tri;
  				++result;
  			}
  		}
  	}
  
  	return result;
  }
  
  static size_t fillVertexCells(unsigned int* table, size_t table_size, unsigned int* vertex_cells, const VertexID* vertex_ids, size_t vertex_count, const bool* locked_vertices)
  {
  	CellHasher hasher = {vertex_ids};
  
  	memset(table, -1, table_size * sizeof(unsigned int));
  
  	size_t result = 0;
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		unsigned int* entry = hashLookup2(table, table_size, hasher, unsigned(i), ~0u);
  
  		// If this vertex was marked as "locked", ignore the hash and vertex id and whatnot, give it it's own cell iod.
  		if (locked_vertices)
  		{
  			if (locked_vertices[i])
  			{
  				vertex_cells[i] = unsigned(result++);
  				continue;
  			}
  		}
  
  		if (*entry == ~0u)
  		{
  			*entry = unsigned(i);
  			vertex_cells[i] = unsigned(result++);
  		}
  		else
  		{
  			vertex_cells[i] = vertex_cells[*entry];
  		}
  	}
  
  	return result;
  }
  
  static size_t countVertexCells(VertexID* table, size_t table_size, const VertexID* vertex_ids, size_t vertex_count)
  {
  	IdHasher hasher;
  
  	memset(table, -1, table_size * sizeof(VertexID));
  
  	size_t result = 0;
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		VertexID id = vertex_ids[i];
  		VertexID* entry = hashLookup2(table, table_size, hasher, id, ~0ull);
  
  		result += (*entry == ~0ull);
  		*entry = id;
  	}
  
  	return result;
  }
  
  static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells)
  {
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		unsigned int i0 = indices[i + 0];
  		unsigned int i1 = indices[i + 1];
  		unsigned int i2 = indices[i + 2];
  
  		unsigned int c0 = vertex_cells[i0];
  		unsigned int c1 = vertex_cells[i1];
  		unsigned int c2 = vertex_cells[i2];
  
  		bool single_cell = (c0 == c1) & (c0 == c2);
  
  		Quadric Q;
  		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], single_cell ? 3.f : 1.f);
  
  		if (single_cell)
  		{
  			quadricAdd(cell_quadrics[c0], Q);
  		}
  		else
  		{
  			quadricAdd(cell_quadrics[c0], Q);
  			quadricAdd(cell_quadrics[c1], Q);
  			quadricAdd(cell_quadrics[c2], Q);
  		}
  	}
  }
  
  static void fillCellQuadrics(Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* vertex_cells)
  {
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		unsigned int c = vertex_cells[i];
  		const Vector3& v = vertex_positions[i];
  
  		Quadric Q;
  		quadricFromPoint(Q, v.x, v.y, v.z, 1.f);
  
  		quadricAdd(cell_quadrics[c], Q);
  	}
  }
  
  static void fillCellRemap(unsigned int* cell_remap, float* cell_errors, size_t cell_count, const unsigned int* vertex_cells, const Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count)
  {
  	memset(cell_remap, -1, cell_count * sizeof(unsigned int));
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		unsigned int cell = vertex_cells[i];
  		float error = quadricError(cell_quadrics[cell], vertex_positions[i]);
  
  		if (cell_remap[cell] == ~0u || cell_errors[cell] > error)
  		{
  			cell_remap[cell] = unsigned(i);
  			cell_errors[cell] = error;
  		}
  	}
  }
  
  static void fillBasicCellRemap(unsigned int* cell_remap, size_t cell_count, const unsigned int* vertex_cells, size_t vertex_count)
  {
  	memset(cell_remap, -1, cell_count * sizeof(unsigned int));
  
  	for (size_t i = 0; i < vertex_count; ++i)
  	{
  		unsigned int cell = vertex_cells[i];
  		if (cell_remap[cell] == ~0u)
  		{
  			cell_remap[cell] = unsigned(i);
  		}
  	}
  }
  
  static size_t filterTriangles(unsigned int* destination, unsigned int* tritable, size_t tritable_size, const unsigned int* indices, size_t index_count, const unsigned int* vertex_cells, const unsigned int* cell_remap, const unsigned int* vertex_cells_pos_only)
  {
  	TriangleHasher hasher = {destination};
  
  	memset(tritable, -1, tritable_size * sizeof(unsigned int));
  
  	size_t result = 0;
  
  	for (size_t i = 0; i < index_count; i += 3)
  	{
  		unsigned int id0 = vertex_cells_pos_only[indices[i + 0]];
  		unsigned int id1 = vertex_cells_pos_only[indices[i + 1]];
  		unsigned int id2 = vertex_cells_pos_only[indices[i + 2]];
  
  		if (id0 != id1 && id0 != id2 && id1 != id2)
  		{
  			unsigned int c0 = vertex_cells[indices[i + 0]];
  			unsigned int c1 = vertex_cells[indices[i + 1]];
  			unsigned int c2 = vertex_cells[indices[i + 2]];
  			assert(c0 != c1 && c0 != c2 && c1 != c2);
  
  			unsigned int a = cell_remap[c0];
  			unsigned int b = cell_remap[c1];
  			unsigned int c = cell_remap[c2];
  
  			if (b < a && b < c)
  			{
  				unsigned int t = a;
  				a = b, b = c, c = t;
  			}
  			else if (c < a && c < b)
  			{
  				unsigned int t = c;
  				c = b, b = a, a = t;
  			}
  
  			destination[result * 3 + 0] = a;
  			destination[result * 3 + 1] = b;
  			destination[result * 3 + 2] = c;
  			
  			unsigned int* entry = hashLookup2(tritable, tritable_size, hasher, unsigned(result), ~0u);
  
  			if (*entry == ~0u)
  				*entry = unsigned(result++);
  		}
  	}
  
  	return result * 3;
  }
  
  static void reconcileVertexAttribute(float* vertex_data, size_t vertex_size, size_t vertex_stride, const unsigned int* indices, size_t index_count, const unsigned int* vertex_cells, const unsigned int* cell_remap)
  {
  	unsigned int vertex_stride_floats = (unsigned int)(vertex_stride / sizeof(float));
  	for (size_t i = 0; i < index_count; ++i)
  	{
  		unsigned int vertex = indices[i];
  		unsigned int best_vertex = cell_remap[vertex_cells[vertex]];
  		memcpy(&vertex_data[vertex_stride_floats * vertex], &vertex_data[vertex_stride_floats * best_vertex], vertex_size);
  	}
  }
  
  static void writeAverageNormals(float* vertex_normals, size_t vertex_normals_stride, Vector3* accumulation_buffer, const unsigned int* new_indices, size_t new_index_count, const unsigned int* original_indices, size_t original_index_count, const unsigned int* vertex_cells, size_t cell_count)
  {
  	// calculate new normals by taking the average of the normals that have been collapsed into one,
  	// weighted by the number of triangles using each one
  	unsigned int normals_stride_floats = (unsigned int)(vertex_normals_stride / sizeof(float));
  	memset(accumulation_buffer, 0, cell_count * sizeof(Vector3));
  	for (size_t i = 0; i < original_index_count; ++i)
  	{
  		unsigned int vertex = original_indices[i];
  		unsigned int cell = vertex_cells[vertex];
  		Vector3& n = accumulation_buffer[cell];
  		n.x += vertex_normals[vertex * normals_stride_floats + 0];
  		n.y += vertex_normals[vertex * normals_stride_floats + 1];
  		n.z += vertex_normals[vertex * normals_stride_floats + 2];
  	}
  
  	// assign new normals back to vertex_normals
  	for (size_t i = 0; i < new_index_count; ++i)
  	{
  		unsigned int vertex = new_indices[i];
  		unsigned int cell = vertex_cells[vertex];
  		Vector3& src = accumulation_buffer[cell];
  		normalize(src);
  		float* dst = &vertex_normals[vertex * normals_stride_floats];
  		dst[0] = src.x;
  		dst[1] = src.y;
  		dst[2] = src.z;
  	}
  }
  
  static float interpolate(float y, float x0, float y0, float x1, float y1, float x2, float y2)
  {
  	// three point interpolation from "revenge of interpolation search" paper
  	float num = (y1 - y) * (x1 - x2) * (x1 - x0) * (y2 - y0);
  	float den = (y2 - y) * (x1 - x2) * (y0 - y1) + (y0 - y) * (x1 - x0) * (y1 - y2);
  	return x1 + num / den;
  }
  
  } // namespace meshopt
  
  #if TRACE
  unsigned char* meshopt_simplifyDebugKind = 0;
  unsigned int* meshopt_simplifyDebugLoop = 0;
  #endif
  
  size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error)
  {
  	using namespace meshopt;
  
  	assert(index_count % 3 == 0);
  	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
  	assert(vertex_positions_stride % sizeof(float) == 0);
  	assert(target_index_count <= index_count);
  
  	meshopt_Allocator allocator;
  
  	unsigned int* result = destination;
  
  	// build adjacency information
  	EdgeAdjacency adjacency = {};
  	buildEdgeAdjacency(adjacency, indices, index_count, vertex_count, allocator);
  
  	// build position remap that maps each vertex to the one with identical position
  	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
  	unsigned int* wedge = allocator.allocate<unsigned int>(vertex_count);
  	buildPositionRemap(remap, wedge, vertex_positions_data, vertex_count, vertex_positions_stride, allocator);
  
  	// classify vertices; vertex kind determines collapse rules, see kCanCollapse
  	unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
  	unsigned int* loop = allocator.allocate<unsigned int>(vertex_count);
  	classifyVertices(vertex_kind, loop, vertex_count, adjacency, remap, wedge);
  
  #if TRACE
  	size_t unique_positions = 0;
  	for (size_t i = 0; i < vertex_count; ++i)
  		unique_positions += remap[i] == i;
  
  	printf("position remap: %d vertices => %d positions\n", int(vertex_count), int(unique_positions));
  
  	size_t kinds[Kind_Count] = {};
  	for (size_t i = 0; i < vertex_count; ++i)
  		kinds[vertex_kind[i]] += remap[i] == i;
  
  	printf("kinds: manifold %d, border %d, seam %d, complex %d, locked %d\n",
  	       int(kinds[Kind_Manifold]), int(kinds[Kind_Border]), int(kinds[Kind_Seam]), int(kinds[Kind_Complex]), int(kinds[Kind_Locked]));
  #endif
  
  	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
  	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
  
  	Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
  	memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
  
  	fillFaceQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap);
  	fillEdgeQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap, vertex_kind, loop);
  
  	if (result != indices)
  		memcpy(result, indices, index_count * sizeof(unsigned int));
  
  #if TRACE
  	size_t pass_count = 0;
  	float worst_error = 0;
  #endif
  
  	Collapse* edge_collapses = allocator.allocate<Collapse>(index_count);
  	unsigned int* collapse_order = allocator.allocate<unsigned int>(index_count);
  	unsigned int* collapse_remap = allocator.allocate<unsigned int>(vertex_count);
  	unsigned char* collapse_locked = allocator.allocate<unsigned char>(vertex_count);
  
  	size_t result_count = index_count;
  
  	// target_error input is linear; we need to adjust it to match quadricError units
  	float error_limit = target_error * target_error;
  
  	while (result_count > target_index_count)
  	{
  		size_t edge_collapse_count = pickEdgeCollapses(edge_collapses, result, result_count, remap, vertex_kind, loop);
  
  		// no edges can be collapsed any more due to topology restrictions
  		if (edge_collapse_count == 0)
  			break;
  
  		rankEdgeCollapses(edge_collapses, edge_collapse_count, vertex_positions, vertex_quadrics, remap);
  
  #if TRACE > 1
  		dumpEdgeCollapses(edge_collapses, edge_collapse_count, vertex_kind);
  #endif
  
  		sortEdgeCollapses(collapse_order, edge_collapses, edge_collapse_count);
  
  		// most collapses remove 2 triangles; use this to establish a bound on the pass in terms of error limit
  		// note that edge_collapse_goal is an estimate; triangle_collapse_goal will be used to actually limit collapses
  		size_t triangle_collapse_goal = (result_count - target_index_count) / 3;
  		size_t edge_collapse_goal = triangle_collapse_goal / 2;
  
  		// we limit the error in each pass based on the error of optimal last collapse; since many collapses will be locked
  		// as they will share vertices with other successfull collapses, we need to increase the acceptable error by this factor
  		const float kPassErrorBound = 1.5f;
  
  		float error_goal = edge_collapse_goal < edge_collapse_count ? edge_collapses[collapse_order[edge_collapse_goal]].error * kPassErrorBound : FLT_MAX;
  
  		for (size_t i = 0; i < vertex_count; ++i)
  			collapse_remap[i] = unsigned(i);
  
  		memset(collapse_locked, 0, vertex_count);
  
  		size_t collapses = performEdgeCollapses(collapse_remap, collapse_locked, vertex_quadrics, edge_collapses, edge_collapse_count, collapse_order, remap, wedge, vertex_kind, triangle_collapse_goal, error_goal, error_limit);
  
  		// no edges can be collapsed any more due to hitting the error limit or triangle collapse limit
  		if (collapses == 0)
  			break;
  
  		remapEdgeLoops(loop, vertex_count, collapse_remap);
  
  		size_t new_count = remapIndexBuffer(result, result_count, collapse_remap);
  		assert(new_count < result_count);
  
  #if TRACE
  		float pass_error = 0.f;
  		for (size_t i = 0; i < edge_collapse_count; ++i)
  		{
  			Collapse& c = edge_collapses[collapse_order[i]];
  
  			if (collapse_remap[c.v0] == c.v1)
  				pass_error = c.error;
  		}
  
  		pass_count++;
  		worst_error = (worst_error < pass_error) ? pass_error : worst_error;
  
  		printf("pass %d: triangles: %d -> %d, collapses: %d/%d (goal: %d), error: %e (limit %e goal %e)\n", int(pass_count), int(result_count / 3), int(new_count / 3), int(collapses), int(edge_collapse_count), int(edge_collapse_goal), pass_error, error_limit, error_goal);
  #endif
  
  		result_count = new_count;
  	}
  
  #if TRACE
  	printf("passes: %d, worst error: %e\n", int(pass_count), worst_error);
  #endif
  
  #if TRACE > 1
  	dumpLockedCollapses(result, result_count, vertex_kind);
  #endif
  
  #if TRACE
  	if (meshopt_simplifyDebugKind)
  		memcpy(meshopt_simplifyDebugKind, vertex_kind, vertex_count);
  
  	if (meshopt_simplifyDebugLoop)
  		memcpy(meshopt_simplifyDebugLoop, loop, vertex_count * sizeof(unsigned int));
  #endif
  
  	return result_count;
  }
  
  size_t meshopt_simplifySloppy(unsigned int* destination, const unsigned int* indices, size_t index_count, 
  	float* vertex_positions_data, float* vertex_normals, float* vertex_uvs, 
  	size_t vertex_count, size_t vertex_positions_stride, size_t vertex_normals_stride, size_t vertex_uvs_stride, size_t target_index_count, const unsigned int* uv_islands, const bool* locked_vertices)
  {
  	using namespace meshopt;
  	assert(index_count % 3 == 0);
  	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
  	assert(vertex_positions_stride % sizeof(float) == 0);
  	assert(vertex_normals_stride % sizeof(float) == 0);
  	assert(vertex_uvs_stride % sizeof(float) == 0);
  	assert(target_index_count <= index_count);
  
  	// we expect to get ~2 triangles/vertex in the output
  	size_t target_cell_count = target_index_count / 6;
  
  	if (target_cell_count == 0)
  		return 0;
  
  	meshopt_Allocator allocator;
  
  	// If true, we will make one full grid and one grid that only considers positions. See big comment below.
  	const bool extra_grid_attributes = (vertex_normals || uv_islands);
  
  	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
  	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
  
  	// find the optimal grid size using guided binary search
  #if TRACE
  	printf("source: %d vertices, %d triangles\n", int(vertex_count), int(index_count / 3));
  	printf("target: %d cells, %d triangles\n", int(target_cell_count), int(target_index_count / 3));
  #endif
  
  	VertexID* vertex_ids = allocator.allocate<VertexID>(vertex_count);
  	VertexID* vertex_ids_pos_only = extra_grid_attributes ? allocator.allocate<VertexID>(vertex_count) : vertex_ids;
  
  	// allocate tritable to be used both for countTriangles (as Triangle*) and filterTriangles (as unsigned int*)
  	size_t tritable_mem_size = hashBuckets2(index_count / 3) * sizeof(Triangle);
  	unsigned char* tritable_mem = allocator.allocate<unsigned char>(tritable_mem_size);
  
  	const int kInterpolationPasses = 5;
  
  	// invariant: # of triangles in min_grid <= target_count
  	int min_grid = 0;
  	int max_grid = (1 << POS_GRID_SIZE_BITS) + 1;
  	size_t min_triangles = 0;
  	size_t max_triangles = index_count / 3;
  
  	// instead of starting in the middle, let's guess as to what the answer might be! triangle count usually grows as a square of grid size...
  	int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
  
  	for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
  	{
  		assert(min_triangles < target_index_count / 3);
  		assert(max_grid - min_grid > 1);
  
  		// we clamp the prediction of the grid size to make sure that the search converges
  		int grid_size = next_grid_size;
  		grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
  
  		// assign IDs to each vertex based on position only, and use this to count triangles.
  		computeVertexIds(vertex_ids, vertex_positions, vertex_normals, uv_islands, vertex_normals_stride, vertex_count, grid_size);
  		if (vertex_normals)
  		{
  			computeVertexIds(vertex_ids_pos_only, vertex_positions, NULL, NULL, 0, vertex_count, grid_size);
  		}
  		size_t triangles = countTriangles((Triangle*)tritable_mem, tritable_mem_size / sizeof(Triangle), vertex_ids, vertex_ids_pos_only, indices, index_count, locked_vertices);
  		assert(triangles <= index_count / 3);
  
  #if TRACE
  		printf("pass %d (%s): grid size %d, triangles %d, %s\n",
  		       pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
  		       grid_size, int(triangles),
  		       (triangles <= target_index_count / 3) ? "under" : "over");
  #endif
  
  		float tip = interpolate(float(target_index_count / 3), float(min_grid), float(min_triangles), float(grid_size), float(triangles), float(max_grid), float(max_triangles));
  
  		if (triangles <= target_index_count / 3)
  		{
  			min_grid = grid_size;
  			min_triangles = triangles;
  		}
  		else
  		{
  			max_grid = grid_size;
  			max_triangles = triangles;
  		}
  
  		if (triangles == target_index_count / 3 || max_grid - min_grid <= 1)
  			break;
  
  		// we start by using interpolation search - it usually converges faster
  		// however, interpolation search has a worst case of O(N) so we switch to binary search after a few iterations which converges in O(logN)
  		next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
  	}
  
  	if (min_triangles == 0)
  		return 0;
  
  	// build vertex->cell association by mapping all vertices with the same quantized position to the same cell
  	size_t table_size = hashBuckets2(vertex_count);
  	unsigned int* table = allocator.allocate<unsigned int>(table_size);
  
  	unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
  	unsigned int* vertex_cells_pos_only = vertex_normals ? allocator.allocate<unsigned int>(vertex_count) : vertex_cells;
  
  	computeVertexIds(vertex_ids, vertex_positions, vertex_normals, uv_islands, vertex_normals_stride, vertex_count, min_grid);
  	size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count, locked_vertices);
  	size_t cell_count_pos_only = cell_count;
  
  	// branching code paths below based on whether or not we have normals:
  	//  - if we do have normals or UV islands, then:
  	//    -- calculate a basic remap first that maps each vertex in a cell
  	//       to the first vertex in that cell. 
  	//    -- ignore normals and UVs when choosing the most representative vertex for each spatial cell
  	//    -- copy that vertex' position onto every normal that is in the same spatial cell.
  	//       This allows us to remove triangles between vertices that share a spatial cell but have
  	//       different normals.
  	//  - if we don't have normals or UV islands, then:
  	//    -- calculate a remapping that maps each vertex in a cell to the representative vertex
  	//       for that cell
  	//    -- do not modify vertex_positions
  
  	unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
  
  	if (extra_grid_attributes)
  	{
  		// get IDs with no normals, so we can ensure vertices in the same spatial cell use the same position
  		// (even if they have different normals).
  		// re-use table...
  		computeVertexIds(vertex_ids_pos_only, vertex_positions, NULL, NULL, 0, vertex_count, min_grid);
  		cell_count_pos_only = fillVertexCells(table, table_size, vertex_cells_pos_only, vertex_ids_pos_only, vertex_count, locked_vertices);
  
  		// build the actual remap table by simply mapping vertices that share a cell (accounting for normals/UVs) to the first vertex in the cell
  		fillBasicCellRemap(cell_remap, cell_count, vertex_cells, vertex_count);
  	}
  
  	// build a quadric for each target spatial cell (ignoring normals/UVs)
  	Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count_pos_only);
  	memset(cell_quadrics, 0, cell_count_pos_only * sizeof(Quadric));
  	fillCellQuadrics(cell_quadrics, indices, index_count, vertex_positions, vertex_cells_pos_only);
  
  	// for each target spatial cell, find the vertex with the minimal error
  	unsigned int* cell_remap_pos_only = extra_grid_attributes ? allocator.allocate<unsigned int>(cell_count_pos_only) : cell_remap;
  	float* cell_errors = allocator.allocate<float>(cell_count_pos_only);
  	fillCellRemap(cell_remap_pos_only, cell_errors, cell_count_pos_only, vertex_cells_pos_only, cell_quadrics, vertex_positions, vertex_count);
  
  	// collapse triangles!
  	// note that we need to filter out triangles that we've already output because we very frequently generate redundant triangles between cells :(
  	// filter triangles with two or more verts in the same spatial cell, regardless of normals
  	size_t filter_tritable_size = hashBuckets2(min_triangles);
  	assert(filter_tritable_size * sizeof(unsigned int) <= tritable_mem_size);
  
  	size_t write = filterTriangles(destination, (unsigned int*)tritable_mem, filter_tritable_size, indices, index_count, vertex_cells, cell_remap, vertex_cells_pos_only);
  	assert(write == min_triangles * 3);
  
  	if (extra_grid_attributes)
  	{
  		// copy positions for unfiltered vertices based on the above remap
  		reconcileVertexAttribute(vertex_positions_data, 3 * sizeof(float), vertex_positions_stride, destination, write, vertex_cells_pos_only, cell_remap_pos_only);
  	}
  
  	if (vertex_normals)
  	{
  		// average out normals that share a cell
  		unsigned char* accumulation_buffer = allocator.allocate<unsigned char>(cell_count * 3 * sizeof(float));
  		writeAverageNormals(vertex_normals, vertex_normals_stride, (Vector3*)accumulation_buffer, destination, write, indices, index_count, vertex_cells, cell_count);
  	}
  
  	if (vertex_uvs && extra_grid_attributes)
  	{
  		// Copy UV coordinates around to match the chosen vertex for each cell, which improves texturing quality.
  
  		// If there is a single connected UV island, then these are already equivalent.
  		unsigned int* vertex_cells_pos_uv = vertex_cells_pos_only;
  		unsigned int* cell_remap_pos_uv = cell_remap_pos_only;
  
  		if (uv_islands)
  		{
  			// Otherwise, generate a new cell remapping specifically for improving UV quality, using only positions and UV island info.
  			// This achieves a much better quality than if we just used the existing remappings that account for normals.
  			// We can just overwrite a bunch of these buffers at this point rather than allocating yet more memory.
  			VertexID* vertex_ids_pos_uv = vertex_ids;
  			computeVertexIds(vertex_ids_pos_uv, vertex_positions, NULL, uv_islands, 0, vertex_count, min_grid);
  
  			size_t cell_count_pos_uv = fillVertexCells(table, table_size, vertex_cells_pos_uv, vertex_ids_pos_uv, vertex_count, locked_vertices);
  				
  			// build a quadric for each target cell
  			Quadric* cell_quadrics_pos_uv = allocator.allocate<Quadric>(cell_count_pos_uv);
  			float* cell_errors_pos_uv = allocator.allocate<float>(cell_count_pos_uv);
  			memset(cell_quadrics_pos_uv, 0, cell_count_pos_uv * sizeof(Quadric));
  			fillCellQuadrics(cell_quadrics_pos_uv, indices, index_count, vertex_positions, vertex_cells_pos_uv);
  	
  			// for each target cell, find the vertex with the minimal error
  			assert(cell_count_pos_uv <= cell_count);
  			cell_remap_pos_uv = cell_remap; // This needs to be a bigger buffer.
  			fillCellRemap(cell_remap_pos_uv, cell_errors_pos_uv, cell_count_pos_uv, vertex_cells_pos_uv, cell_quadrics_pos_uv, vertex_positions, vertex_count);
  		}
  
  		// Do the actual copy.
  		reconcileVertexAttribute(vertex_uvs, 2 * sizeof(float), vertex_uvs_stride, destination, write, vertex_cells_pos_uv, cell_remap_pos_uv);
  	}
  
  #if TRACE
  	printf("result: %d cells, %d triangles (%d unfiltered)\n", int(cell_count), int(write / 3), int(min_triangles));
  #endif
  
  	return write;
  }
  
  size_t meshopt_simplifyPoints(unsigned int* destination, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, size_t target_vertex_count)
  {
  	using namespace meshopt;
  
  	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
  	assert(vertex_positions_stride % sizeof(float) == 0);
  	assert(target_vertex_count <= vertex_count);
  
  	size_t target_cell_count = target_vertex_count;
  
  	if (target_cell_count == 0)
  		return 0;
  
  	meshopt_Allocator allocator;
  
  	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
  	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
  
  	// find the optimal grid size using guided binary search
  #if TRACE
  	printf("source: %d vertices\n", int(vertex_count));
  	printf("target: %d cells\n", int(target_cell_count));
  #endif
  
  	VertexID* vertex_ids = allocator.allocate<VertexID>(vertex_count);
  
  	size_t table_size = hashBuckets2(vertex_count);
  	unsigned char* table_mem = allocator.allocate<unsigned char>(table_size * sizeof(VertexID));	
  
  	const int kInterpolationPasses = 5;
  
  	// invariant: # of vertices in min_grid <= target_count
  	int min_grid = 0;
  	int max_grid = 1025;
  	size_t min_vertices = 0;
  	size_t max_vertices = vertex_count;
  
  	// instead of starting in the middle, let's guess as to what the answer might be! triangle count usually grows as a square of grid size...
  	int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
  
  	for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
  	{
  		assert(min_vertices < target_vertex_count);
  		assert(max_grid - min_grid > 1);
  
  		// we clamp the prediction of the grid size to make sure that the search converges
  		int grid_size = next_grid_size;
  		grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
  
  		computeVertexIds(vertex_ids, vertex_positions, NULL, NULL, 0, vertex_count, grid_size);
  		size_t vertices = countVertexCells((VertexID*)table_mem, table_size, vertex_ids, vertex_count);
  
  #if TRACE
  		printf("pass %d (%s): grid size %d, vertices %d, %s\n",
  		       pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
  		       grid_size, int(vertices),
  		       (vertices <= target_vertex_count) ? "under" : "over");
  #endif
  
  		float tip = interpolate(float(target_vertex_count), float(min_grid), float(min_vertices), float(grid_size), float(vertices), float(max_grid), float(max_vertices));
  
  		if (vertices <= target_vertex_count)
  		{
  			min_grid = grid_size;
  			min_vertices = vertices;
  		}
  		else
  		{
  			max_grid = grid_size;
  			max_vertices = vertices;
  		}
  
  		if (vertices == target_vertex_count || max_grid - min_grid <= 1)
  			break;
  
  		// we start by using interpolation search - it usually converges faster
  		// however, interpolation search has a worst case of O(N) so we switch to binary search after a few iterations which converges in O(logN)
  		next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
  	}
  
  	if (min_vertices == 0)
  		return 0;
  
  	// build vertex->cell association by mapping all vertices with the same quantized position to the same cell
  	unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
  
  	computeVertexIds(vertex_ids, vertex_positions, NULL, NULL, 0, vertex_count, min_grid);
  	size_t cell_count = fillVertexCells((unsigned int*)table_mem, table_size, vertex_cells, vertex_ids, vertex_count, NULL);
  
  	// build a quadric for each target cell
  	Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
  	memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
  
  	fillCellQuadrics(cell_quadrics, vertex_positions, vertex_count, vertex_cells);
  
  	// for each target cell, find the vertex with the minimal error
  	unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
  	float* cell_errors = allocator.allocate<float>(cell_count);
  
  	fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
  
  	// copy results to the output
  	assert(cell_count <= target_vertex_count);
  	memcpy(destination, cell_remap, sizeof(unsigned int) * cell_count);
  
  #if TRACE
  	printf("result: %d cells\n", int(cell_count));
  #endif
  
  	return cell_count;
  }
  