Update with experimental mesh optimizer.

Normals being optimized has better quality now.

Test simplify once and then use a slightly less tolerant 
error for the target error.

4 files changed, 418 insertions, 4 deletions

diff --git a/thirdparty/README.md b/thirdparty/README.md
index 7cedd1a2cb..d3599fd195 100644
--- a/thirdparty/README.md
+++ b/thirdparty/README.md
@@ -375,6 +375,9 @@ Files extracted from upstream repository:
 - All files in `src/`.
 - `LICENSE.md`.
 
+An [experimental upstream feature](https://github.com/zeux/meshoptimizer/tree/simplify-attr),
+has been backported, see patch in `patches` directory.
+
 
 ## miniupnpc
 
diff --git a/thirdparty/meshoptimizer/meshoptimizer.h b/thirdparty/meshoptimizer/meshoptimizer.h
index fe8d349731..e44b99ce52 100644
--- a/thirdparty/meshoptimizer/meshoptimizer.h
+++ b/thirdparty/meshoptimizer/meshoptimizer.h
@@ -299,6 +299,11 @@ MESHOPTIMIZER_EXPERIMENTAL void meshopt_decodeFilterExp(void* buffer, size_t ver
 MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error, float* result_error);
 
 /**
+ * Experimental: Mesh simplifier with attribute metric; attributes follow xyz position data atm (vertex data must contain 3 + attribute_count floats per vertex)
+ */
+MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_data, size_t vertex_count, size_t vertex_stride, size_t target_index_count, float target_error, float* result_error, const float* attributes, const float* attribute_weights, size_t attribute_count);
+
+/**
  * Experimental: Mesh simplifier (sloppy)
  * Reduces the number of triangles in the mesh, sacrificing mesh apperance for simplification performance
  * The algorithm doesn't preserve mesh topology but can stop short of the target goal based on target error.
diff --git a/thirdparty/meshoptimizer/patches/attribute-aware-simplify.patch b/thirdparty/meshoptimizer/patches/attribute-aware-simplify.patch
new file mode 100644
index 0000000000..cf648b0da3
--- /dev/null
+++ b/thirdparty/meshoptimizer/patches/attribute-aware-simplify.patch
@@ -0,0 +1,262 @@
+diff --git a/thirdparty/meshoptimizer/meshoptimizer.h b/thirdparty/meshoptimizer/meshoptimizer.h
+index fe8d349731..e44b99ce52 100644
+--- a/thirdparty/meshoptimizer/meshoptimizer.h
++++ b/thirdparty/meshoptimizer/meshoptimizer.h
+@@ -298,6 +298,11 @@ MESHOPTIMIZER_EXPERIMENTAL void meshopt_decodeFilterExp(void* buffer, size_t ver
+  */
+ MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error, float* result_error);
+ 
++/**
++ * Experimental: Mesh simplifier with attribute metric; attributes follow xyz position data atm (vertex data must contain 3 + attribute_count floats per vertex)
++ */
++MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_data, size_t vertex_count, size_t vertex_stride, size_t target_index_count, float target_error, float* result_error, const float* attributes, const float* attribute_weights, size_t attribute_count);
++
+ /**
+  * Experimental: Mesh simplifier (sloppy)
+  * Reduces the number of triangles in the mesh, sacrificing mesh apperance for simplification performance
+diff --git a/thirdparty/meshoptimizer/simplifier.cpp b/thirdparty/meshoptimizer/simplifier.cpp
+index b2cb589462..059cabb055 100644
+--- a/thirdparty/meshoptimizer/simplifier.cpp
++++ b/thirdparty/meshoptimizer/simplifier.cpp
+@@ -20,6 +20,8 @@
+ #define TRACESTATS(i) (void)0
+ #endif
+ 
++#define ATTRIBUTES 8
++
+ // 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
+@@ -358,6 +360,10 @@ static void classifyVertices(unsigned char* result, unsigned int* loop, unsigned
+ struct Vector3
+ {
+ 	float x, y, z;
++
++#if ATTRIBUTES
++	float a[ATTRIBUTES];
++#endif
+ };
+ 
+ static float rescalePositions(Vector3* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
+@@ -414,6 +420,13 @@ struct Quadric
+ 	float a10, a20, a21;
+ 	float b0, b1, b2, c;
+ 	float w;
++
++#if ATTRIBUTES
++	float gx[ATTRIBUTES];
++	float gy[ATTRIBUTES];
++	float gz[ATTRIBUTES];
++	float gw[ATTRIBUTES];
++#endif
+ };
+ 
+ struct Collapse
+@@ -456,6 +469,16 @@ static void quadricAdd(Quadric& Q, const Quadric& R)
+ 	Q.b2 += R.b2;
+ 	Q.c += R.c;
+ 	Q.w += R.w;
++
++#if ATTRIBUTES
++	for (int k = 0; k < ATTRIBUTES; ++k)
++	{
++		Q.gx[k] += R.gx[k];
++		Q.gy[k] += R.gy[k];
++		Q.gz[k] += R.gz[k];
++		Q.gw[k] += R.gw[k];
++	}
++#endif
+ }
+ 
+ static float quadricError(const Quadric& Q, const Vector3& v)
+@@ -481,6 +504,17 @@ static float quadricError(const Quadric& Q, const Vector3& v)
+ 	r += ry * v.y;
+ 	r += rz * v.z;
+ 
++#if ATTRIBUTES
++	// see quadricUpdateAttributes for general derivation; here we need to add the parts of (eval(pos) - attr)^2 that depend on attr
++	for (int k = 0; k < ATTRIBUTES; ++k)
++	{
++		float a = v.a[k];
++
++		r += a * a * Q.w;
++		r -= 2 * a * (v.x * Q.gx[k] + v.y * Q.gy[k] + v.z * Q.gz[k] + Q.gw[k]);
++	}
++#endif
++
+ 	float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
+ 
+ 	return fabsf(r) * s;
+@@ -504,6 +538,13 @@ static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d, flo
+ 	Q.b2 = c * dw;
+ 	Q.c = d * dw;
+ 	Q.w = w;
++
++#if ATTRIBUTES
++	memset(Q.gx, 0, sizeof(Q.gx));
++	memset(Q.gy, 0, sizeof(Q.gy));
++	memset(Q.gz, 0, sizeof(Q.gz));
++	memset(Q.gw, 0, sizeof(Q.gw));
++#endif
+ }
+ 
+ static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
+@@ -556,6 +597,84 @@ static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3
+ 	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
+ }
+ 
++#if ATTRIBUTES
++static void quadricUpdateAttributes(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float w)
++{
++	// for each attribute we want to encode the following function into the quadric:
++	// (eval(pos) - attr)^2
++	// where eval(pos) interpolates attribute across the triangle like so:
++	// eval(pos) = pos.x * gx + pos.y * gy + pos.z * gz + gw
++	// where gx/gy/gz/gw are gradients
++	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};
++
++	// we compute gradients using barycentric coordinates; barycentric coordinates can be computed as follows:
++	// v = (d11 * d20 - d01 * d21) / denom
++	// w = (d00 * d21 - d01 * d20) / denom
++	// u = 1 - v - w
++	// here v0, v1 are triangle edge vectors, v2 is a vector from point to triangle corner, and dij = dot(vi, vj)
++	const Vector3& v0 = p10;
++	const Vector3& v1 = p20;
++	float d00 = v0.x * v0.x + v0.y * v0.y + v0.z * v0.z;
++	float d01 = v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
++	float d11 = v1.x * v1.x + v1.y * v1.y + v1.z * v1.z;
++	float denom = d00 * d11 - d01 * d01;
++	float denomr = denom == 0 ? 0.f : 1.f / denom;
++
++	// precompute gradient factors
++	// these are derived by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w and factoring out common factors that are shared between attributes
++	float gx1 = (d11 * v0.x - d01 * v1.x) * denomr;
++	float gx2 = (d00 * v1.x - d01 * v0.x) * denomr;
++	float gy1 = (d11 * v0.y - d01 * v1.y) * denomr;
++	float gy2 = (d00 * v1.y - d01 * v0.y) * denomr;
++	float gz1 = (d11 * v0.z - d01 * v1.z) * denomr;
++	float gz2 = (d00 * v1.z - d01 * v0.z) * denomr;
++
++	for (int k = 0; k < ATTRIBUTES; ++k)
++	{
++		float a0 = p0.a[k], a1 = p1.a[k], a2 = p2.a[k];
++
++		// compute gradient of eval(pos) for x/y/z/w
++		// the formulas below are obtained by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w
++		float gx = gx1 * (a1 - a0) + gx2 * (a2 - a0);
++		float gy = gy1 * (a1 - a0) + gy2 * (a2 - a0);
++		float gz = gz1 * (a1 - a0) + gz2 * (a2 - a0);
++		float gw = a0 - p0.x * gx - p0.y * gy - p0.z * gz;
++
++		// quadric encodes (eval(pos)-attr)^2; this means that the resulting expansion needs to compute, for example, pos.x * pos.y * K
++		// since quadrics already encode factors for pos.x * pos.y, we can accumulate almost everything in basic quadric fields
++		Q.a00 += w * (gx * gx);
++		Q.a11 += w * (gy * gy);
++		Q.a22 += w * (gz * gz);
++
++		Q.a10 += w * (gy * gx);
++		Q.a20 += w * (gz * gx);
++		Q.a21 += w * (gz * gy);
++
++		Q.b0 += w * (gx * gw);
++		Q.b1 += w * (gy * gw);
++		Q.b2 += w * (gz * gw);
++
++		Q.c += w * (gw * gw);
++
++		// the only remaining sum components are ones that depend on attr; these will be addded during error evaluation, see quadricError
++		Q.gx[k] = w * gx;
++		Q.gy[k] = w * gy;
++		Q.gz[k] = w * gz;
++		Q.gw[k] = w * gw;
++
++#if TRACE > 2
++		printf("attr%d: %e %e %e\n",
++			k,
++			(gx * p0.x + gy * p0.y + gz * p0.z + gw - a0),
++			(gx * p1.x + gy * p1.y + gz * p1.z + gw - a1),
++			(gx * p2.x + gy * p2.y + gz * p2.z + gw - a2)
++			);
++#endif
++	}
++}
++#endif
++
+ 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)
+@@ -567,6 +686,9 @@ static void fillFaceQuadrics(Quadric* vertex_quadrics, const unsigned int* indic
+ 		Quadric Q;
+ 		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
+ 
++#if ATTRIBUTES
++		quadricUpdateAttributes(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], Q.w);
++#endif
+ 		quadricAdd(vertex_quadrics[remap[i0]], Q);
+ 		quadricAdd(vertex_quadrics[remap[i1]], Q);
+ 		quadricAdd(vertex_quadrics[remap[i2]], Q);
+@@ -1259,13 +1381,19 @@ unsigned int* meshopt_simplifyDebugLoopBack = 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, float* out_result_error)
++{
++	return meshopt_simplifyWithAttributes(destination, indices, index_count, vertex_positions_data, vertex_count, vertex_positions_stride, target_index_count, target_error, out_result_error, 0, 0, 0);
++}
++
++size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_data, size_t vertex_count, size_t vertex_stride, size_t target_index_count, float target_error, float* out_result_error, const float* attributes, const float* attribute_weights, size_t attribute_count)
+ {
+ 	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_stride > 0 && vertex_stride <= 256);
++	assert(vertex_stride % sizeof(float) == 0);
+ 	assert(target_index_count <= index_count);
++	assert(attribute_count <= ATTRIBUTES);
+ 
+ 	meshopt_Allocator allocator;
+ 
+@@ -1279,7 +1407,7 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
+ 	// 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);
++	buildPositionRemap(remap, wedge, vertex_data, vertex_count, vertex_stride, allocator);
+ 
+ 	// classify vertices; vertex kind determines collapse rules, see kCanCollapse
+ 	unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
+@@ -1303,7 +1431,21 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
+ #endif
+ 
+ 	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
+-	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
++	rescalePositions(vertex_positions, vertex_data, vertex_count, vertex_stride);
++
++#if ATTRIBUTES
++	for (size_t i = 0; i < vertex_count; ++i)
++	{
++		memset(vertex_positions[i].a, 0, sizeof(vertex_positions[i].a));
++
++		for (size_t k = 0; k < attribute_count; ++k)
++		{
++			float a = attributes[i * attribute_count + k];
++
++			vertex_positions[i].a[k] = a * attribute_weights[k];
++		}
++	}
++#endif
+ 
+ 	Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
+ 	memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
+@@ -1395,7 +1537,9 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
+ 
+ 	// result_error is quadratic; we need to remap it back to linear
+ 	if (out_result_error)
++	{
+ 		*out_result_error = sqrtf(result_error);
++	}
+ 
+ 	return result_count;
+ }
diff --git a/thirdparty/meshoptimizer/simplifier.cpp b/thirdparty/meshoptimizer/simplifier.cpp
index b2cb589462..059cabb055 100644
--- a/thirdparty/meshoptimizer/simplifier.cpp
+++ b/thirdparty/meshoptimizer/simplifier.cpp
@@ -20,6 +20,8 @@
 #define TRACESTATS(i) (void)0
 #endif
 
+#define ATTRIBUTES 8
+
 // 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
@@ -358,6 +360,10 @@ static void classifyVertices(unsigned char* result, unsigned int* loop, unsigned
 struct Vector3
 {
 	float x, y, z;
+
+#if ATTRIBUTES
+	float a[ATTRIBUTES];
+#endif
 };
 
 static float rescalePositions(Vector3* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
@@ -414,6 +420,13 @@ struct Quadric
 	float a10, a20, a21;
 	float b0, b1, b2, c;
 	float w;
+
+#if ATTRIBUTES
+	float gx[ATTRIBUTES];
+	float gy[ATTRIBUTES];
+	float gz[ATTRIBUTES];
+	float gw[ATTRIBUTES];
+#endif
 };
 
 struct Collapse
@@ -456,6 +469,16 @@ static void quadricAdd(Quadric& Q, const Quadric& R)
 	Q.b2 += R.b2;
 	Q.c += R.c;
 	Q.w += R.w;
+
+#if ATTRIBUTES
+	for (int k = 0; k < ATTRIBUTES; ++k)
+	{
+		Q.gx[k] += R.gx[k];
+		Q.gy[k] += R.gy[k];
+		Q.gz[k] += R.gz[k];
+		Q.gw[k] += R.gw[k];
+	}
+#endif
 }
 
 static float quadricError(const Quadric& Q, const Vector3& v)
@@ -481,6 +504,17 @@ static float quadricError(const Quadric& Q, const Vector3& v)
 	r += ry * v.y;
 	r += rz * v.z;
 
+#if ATTRIBUTES
+	// see quadricUpdateAttributes for general derivation; here we need to add the parts of (eval(pos) - attr)^2 that depend on attr
+	for (int k = 0; k < ATTRIBUTES; ++k)
+	{
+		float a = v.a[k];
+
+		r += a * a * Q.w;
+		r -= 2 * a * (v.x * Q.gx[k] + v.y * Q.gy[k] + v.z * Q.gz[k] + Q.gw[k]);
+	}
+#endif
+
 	float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
 
 	return fabsf(r) * s;
@@ -504,6 +538,13 @@ static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d, flo
 	Q.b2 = c * dw;
 	Q.c = d * dw;
 	Q.w = w;
+
+#if ATTRIBUTES
+	memset(Q.gx, 0, sizeof(Q.gx));
+	memset(Q.gy, 0, sizeof(Q.gy));
+	memset(Q.gz, 0, sizeof(Q.gz));
+	memset(Q.gw, 0, sizeof(Q.gw));
+#endif
 }
 
 static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
@@ -556,6 +597,84 @@ static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3
 	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
 }
 
+#if ATTRIBUTES
+static void quadricUpdateAttributes(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float w)
+{
+	// for each attribute we want to encode the following function into the quadric:
+	// (eval(pos) - attr)^2
+	// where eval(pos) interpolates attribute across the triangle like so:
+	// eval(pos) = pos.x * gx + pos.y * gy + pos.z * gz + gw
+	// where gx/gy/gz/gw are gradients
+	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};
+
+	// we compute gradients using barycentric coordinates; barycentric coordinates can be computed as follows:
+	// v = (d11 * d20 - d01 * d21) / denom
+	// w = (d00 * d21 - d01 * d20) / denom
+	// u = 1 - v - w
+	// here v0, v1 are triangle edge vectors, v2 is a vector from point to triangle corner, and dij = dot(vi, vj)
+	const Vector3& v0 = p10;
+	const Vector3& v1 = p20;
+	float d00 = v0.x * v0.x + v0.y * v0.y + v0.z * v0.z;
+	float d01 = v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
+	float d11 = v1.x * v1.x + v1.y * v1.y + v1.z * v1.z;
+	float denom = d00 * d11 - d01 * d01;
+	float denomr = denom == 0 ? 0.f : 1.f / denom;
+
+	// precompute gradient factors
+	// these are derived by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w and factoring out common factors that are shared between attributes
+	float gx1 = (d11 * v0.x - d01 * v1.x) * denomr;
+	float gx2 = (d00 * v1.x - d01 * v0.x) * denomr;
+	float gy1 = (d11 * v0.y - d01 * v1.y) * denomr;
+	float gy2 = (d00 * v1.y - d01 * v0.y) * denomr;
+	float gz1 = (d11 * v0.z - d01 * v1.z) * denomr;
+	float gz2 = (d00 * v1.z - d01 * v0.z) * denomr;
+
+	for (int k = 0; k < ATTRIBUTES; ++k)
+	{
+		float a0 = p0.a[k], a1 = p1.a[k], a2 = p2.a[k];
+
+		// compute gradient of eval(pos) for x/y/z/w
+		// the formulas below are obtained by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w
+		float gx = gx1 * (a1 - a0) + gx2 * (a2 - a0);
+		float gy = gy1 * (a1 - a0) + gy2 * (a2 - a0);
+		float gz = gz1 * (a1 - a0) + gz2 * (a2 - a0);
+		float gw = a0 - p0.x * gx - p0.y * gy - p0.z * gz;
+
+		// quadric encodes (eval(pos)-attr)^2; this means that the resulting expansion needs to compute, for example, pos.x * pos.y * K
+		// since quadrics already encode factors for pos.x * pos.y, we can accumulate almost everything in basic quadric fields
+		Q.a00 += w * (gx * gx);
+		Q.a11 += w * (gy * gy);
+		Q.a22 += w * (gz * gz);
+
+		Q.a10 += w * (gy * gx);
+		Q.a20 += w * (gz * gx);
+		Q.a21 += w * (gz * gy);
+
+		Q.b0 += w * (gx * gw);
+		Q.b1 += w * (gy * gw);
+		Q.b2 += w * (gz * gw);
+
+		Q.c += w * (gw * gw);
+
+		// the only remaining sum components are ones that depend on attr; these will be addded during error evaluation, see quadricError
+		Q.gx[k] = w * gx;
+		Q.gy[k] = w * gy;
+		Q.gz[k] = w * gz;
+		Q.gw[k] = w * gw;
+
+#if TRACE > 2
+		printf("attr%d: %e %e %e\n",
+			k,
+			(gx * p0.x + gy * p0.y + gz * p0.z + gw - a0),
+			(gx * p1.x + gy * p1.y + gz * p1.z + gw - a1),
+			(gx * p2.x + gy * p2.y + gz * p2.z + gw - a2)
+			);
+#endif
+	}
+}
+#endif
+
 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)
@@ -567,6 +686,9 @@ static void fillFaceQuadrics(Quadric* vertex_quadrics, const unsigned int* indic
 		Quadric Q;
 		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
 
+#if ATTRIBUTES
+		quadricUpdateAttributes(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], Q.w);
+#endif
 		quadricAdd(vertex_quadrics[remap[i0]], Q);
 		quadricAdd(vertex_quadrics[remap[i1]], Q);
 		quadricAdd(vertex_quadrics[remap[i2]], Q);
@@ -1260,12 +1382,18 @@ unsigned int* meshopt_simplifyDebugLoopBack = 0;
 
 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, float* out_result_error)
 {
+	return meshopt_simplifyWithAttributes(destination, indices, index_count, vertex_positions_data, vertex_count, vertex_positions_stride, target_index_count, target_error, out_result_error, 0, 0, 0);
+}
+
+size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_data, size_t vertex_count, size_t vertex_stride, size_t target_index_count, float target_error, float* out_result_error, const float* attributes, const float* attribute_weights, size_t attribute_count)
+{
 	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_stride > 0 && vertex_stride <= 256);
+	assert(vertex_stride % sizeof(float) == 0);
 	assert(target_index_count <= index_count);
+	assert(attribute_count <= ATTRIBUTES);
 
 	meshopt_Allocator allocator;
 
@@ -1279,7 +1407,7 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
 	// 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);
+	buildPositionRemap(remap, wedge, vertex_data, vertex_count, vertex_stride, allocator);
 
 	// classify vertices; vertex kind determines collapse rules, see kCanCollapse
 	unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
@@ -1303,7 +1431,21 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
 #endif
 
 	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
-	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
+	rescalePositions(vertex_positions, vertex_data, vertex_count, vertex_stride);
+
+#if ATTRIBUTES
+	for (size_t i = 0; i < vertex_count; ++i)
+	{
+		memset(vertex_positions[i].a, 0, sizeof(vertex_positions[i].a));
+
+		for (size_t k = 0; k < attribute_count; ++k)
+		{
+			float a = attributes[i * attribute_count + k];
+
+			vertex_positions[i].a[k] = a * attribute_weights[k];
+		}
+	}
+#endif
 
 	Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
 	memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
@@ -1395,7 +1537,9 @@ size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices,
 
 	// result_error is quadratic; we need to remap it back to linear
 	if (out_result_error)
+	{
 		*out_result_error = sqrtf(result_error);
+	}
 
 	return result_count;
 }