Merge pull request #853 from zeux/bsphere

clusterizer: Expose sphere bounds computation

Author: zeux

demo/nanite.cpp

@@ -86,35 +86,17 @@ static LODBounds bounds(const std::vector<Vertex>& vertices, const std::vector<u
 
 static LODBounds boundsMerge(const std::vector<Cluster>& clusters, const std::vector<int>& group)
 {
-	LODBounds result = {};
-
-	// we approximate merged bounds center as weighted average of cluster centers
-	// (could also use bounds() center, but we can't use bounds() radius so might as well just merge manually)
-	float weight = 0.f;
+	std::vector<LODBounds> bounds(group.size());
 	for (size_t j = 0; j < group.size(); ++j)
-	{
-		result.center[0] += clusters[group[j]].self.center[0] * clusters[group[j]].self.radius;
-		result.center[1] += clusters[group[j]].self.center[1] * clusters[group[j]].self.radius;
-		result.center[2] += clusters[group[j]].self.center[2] * clusters[group[j]].self.radius;
-		weight += clusters[group[j]].self.radius;
-	}
+		bounds[j] = clusters[group[j]].self;
 
-	if (weight > 0)
-	{
-		result.center[0] /= weight;
-		result.center[1] /= weight;
-		result.center[2] /= weight;
-	}
+	meshopt_Bounds merged = meshopt_computeSphereBounds(&bounds[0].center[0], bounds.size(), sizeof(LODBounds), &bounds[0].radius, sizeof(LODBounds));
 
-	// merged bounds must strictly contain all cluster bounds
-	result.radius = 0.f;
-	for (size_t j = 0; j < group.size(); ++j)
-	{
-		float dx = clusters[group[j]].self.center[0] - result.center[0];
-		float dy = clusters[group[j]].self.center[1] - result.center[1];
-		float dz = clusters[group[j]].self.center[2] - result.center[2];
-		result.radius = std::max(result.radius, clusters[group[j]].self.radius + sqrtf(dx * dx + dy * dy + dz * dz));
-	}
+	LODBounds result = {};
+	result.center[0] = merged.center[0];
+	result.center[1] = merged.center[1];
+	result.center[2] = merged.center[2];
+	result.radius = merged.radius;
 
 	// merged bounds error must be conservative wrt cluster errors
 	result.error = 0.f;

demo/tests.cpp

@@ -1059,6 +1059,36 @@ static void clusterBoundsDegenerate()
 	assert(bounds2.center[2] - bounds2.radius <= 0 && bounds2.center[2] + bounds2.radius >= 1);
 }
 
+static void sphereBounds()
+{
+	const float vbr[] = {
+	    0, 0, 0, 0,
+	    0, 1, 0, 1,
+	    0, 0, 1, 2,
+	    1, 0, 1, 3, // clang-format
+	};
+
+	// without the radius, the center is somewhere inside the tetrahedron
+	// note that we currently compute a somewhat suboptimal sphere here due to the tetrahedron being perfecly axis aligned
+	meshopt_Bounds bounds = meshopt_computeSphereBounds(vbr, 4, sizeof(float) * 4, NULL, 0);
+	assert(bounds.radius < 0.97f);
+
+	// forcing a better initial guess for the center by using different radii produces a close to optimal sphere
+	float eps[4] = {1e-3f, 2e-3f, 3e-3f, 4e-3f};
+	meshopt_Bounds boundse = meshopt_computeSphereBounds(vbr, 4, sizeof(float) * 4, eps, sizeof(float));
+	assert(fabsf(boundse.center[0] - 0.5f) < 1e-2f);
+	assert(fabsf(boundse.center[1] - 0.5f) < 1e-2f);
+	assert(fabsf(boundse.center[2] - 0.5f) < 1e-2f);
+	assert(boundse.radius < 0.87f);
+
+	// when using the radius, the last sphere envelops the entire set
+	meshopt_Bounds boundsr = meshopt_computeSphereBounds(vbr, 4, sizeof(float) * 4, vbr + 3, sizeof(float) * 4);
+	assert(fabsf(boundsr.center[0] - 1.f) < 1e-2f);
+	assert(fabsf(boundsr.center[1] - 0.f) < 1e-2f);
+	assert(fabsf(boundsr.center[2] - 1.f) < 1e-2f);
+	assert(fabsf(boundsr.radius - 3.f) < 1e-2f);
+}
+
 static void meshletsEmpty()
 {
 	const float vbd[4 * 3] = {};
@@ -2256,6 +2286,7 @@ void runTests()
 	encodeFilterExpClamp();
 
 	clusterBoundsDegenerate();
+	sphereBounds();
 
 	meshletsEmpty();
 	meshletsDense();

src/clusterizer.cpp

@@ -145,67 +145,82 @@ static void buildTriangleAdjacencySparse(TriangleAdjacency2& adjacency, const un
 	}
 }
 
-static void computeBoundingSphere(float result[4], const float points[][3], size_t count)
+static void computeBoundingSphere(float result[4], const float* points, size_t count, size_t points_stride, const float* radii, size_t radii_stride)
 {
 	assert(count > 0);
 
+	size_t points_stride_float = points_stride / sizeof(float);
+	size_t radii_stride_float = radii_stride / sizeof(float);
+
 	// find extremum points along all 3 axes; for each axis we get a pair of points with min/max coordinates
 	size_t pmin[3] = {0, 0, 0};
 	size_t pmax[3] = {0, 0, 0};
 
 	for (size_t i = 0; i < count; ++i)
 	{
-		const float* p = points[i];
+		const float* p = points + i * points_stride_float;
+		float r = radii[i * radii_stride_float];
 
 		for (int axis = 0; axis < 3; ++axis)
 		{
-			pmin[axis] = (p[axis] < points[pmin[axis]][axis]) ? i : pmin[axis];
-			pmax[axis] = (p[axis] > points[pmax[axis]][axis]) ? i : pmax[axis];
+			float bmin = points[pmin[axis] * points_stride_float + axis] - radii[pmin[axis] * radii_stride_float];
+			float bmax = points[pmax[axis] * points_stride_float + axis] + radii[pmax[axis] * radii_stride_float];
+
+			pmin[axis] = (p[axis] - r < bmin) ? i : pmin[axis];
+			pmax[axis] = (p[axis] + r > bmax) ? i : pmax[axis];
 		}
 	}
 
 	// find the pair of points with largest distance
-	float paxisd2 = 0;
 	int paxis = 0;
+	float paxisdr = 0;
 
 	for (int axis = 0; axis < 3; ++axis)
 	{
-		const float* p1 = points[pmin[axis]];
-		const float* p2 = points[pmax[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 (d2 > paxisd2)
+		if (dr > paxisdr)
 		{
-			paxisd2 = d2;
+			paxisdr = dr;
 			paxis = axis;
 		}
 	}
 
 	// use the longest segment as the initial sphere diameter
-	const float* p1 = points[pmin[paxis]];
-	const float* p2 = points[pmax[paxis]];
+	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]) / 2, (p1[1] + p2[1]) / 2, (p1[2] + p2[2]) / 2};
-	float radius = sqrtf(paxisd2) / 2;
+	float center[3] = {p1[0] + (p2[0] - p1[0]) * paxisk, p1[1] + (p2[1] - p1[1]) * paxisk, p1[2] + (p2[2] - p1[2]) * paxisk};
+	float radius = paxisdr / 2;
 
 	// iteratively adjust the sphere up until all points fit
 	for (size_t i = 0; i < count; ++i)
 	{
-		const float* p = points[i];
+		const float* p = points + i * points_stride_float;
+		float r = radii[i * radii_stride_float];
+
 		float d2 = (p[0] - center[0]) * (p[0] - center[0]) + (p[1] - center[1]) * (p[1] - center[1]) + (p[2] - center[2]) * (p[2] - center[2]);
+		float d = sqrtf(d2);
 
-		if (d2 > radius * radius)
+		if (d + r > radius)
 		{
-			float d = sqrtf(d2);
-			assert(d > 0);
-
-			float k = 0.5f + (radius / d) / 2;
+			float k = d > 0 ? (d + r - radius) / (2 * d) : 0.f;
 
-			center[0] = center[0] * k + p[0] * (1 - k);
-			center[1] = center[1] * k + p[1] * (1 - k);
-			center[2] = center[2] * k + p[2] * (1 - k);
-			radius = (radius + d) / 2;
+			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;
 		}
 	}
 
@@ -1004,15 +1019,17 @@ meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t
 	if (triangles == 0)
 		return bounds;
 
+	const float rzero = 0.f;
+
 	// compute cluster bounding sphere; we'll use the center to determine normal cone apex as well
 	float psphere[4] = {};
-	computeBoundingSphere(psphere, corners[0], triangles * 3);
+	computeBoundingSphere(psphere, corners[0][0], triangles * 3, sizeof(float) * 3, &rzero, 0);
 
 	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, triangles);
+	computeBoundingSphere(nsphere, normals[0], triangles, sizeof(float) * 3, &rzero, 0);
 
 	float axis[3] = {nsphere[0], nsphere[1], nsphere[2]};
 	float axislength = sqrtf(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]);
@@ -1122,6 +1139,33 @@ meshopt_Bounds meshopt_computeMeshletBounds(const unsigned int* meshlet_vertices
 	return meshopt_computeClusterBounds(indices, triangle_count * 3, vertex_positions, vertex_count, vertex_positions_stride);
 }
 
+meshopt_Bounds meshopt_computeSphereBounds(const float* positions, size_t count, size_t positions_stride, const float* radii, size_t radii_stride)
+{
+	using namespace meshopt;
+
+	assert(positions_stride >= 12 && positions_stride <= 256);
+	assert(positions_stride % sizeof(float) == 0);
+	assert((radii_stride >= 4 && radii_stride <= 256) || radii == NULL);
+	assert(radii_stride % sizeof(float) == 0);
+
+	meshopt_Bounds bounds = {};
+
+	if (count == 0)
+		return bounds;
+
+	const float rzero = 0.f;
+
+	float psphere[4] = {};
+	computeBoundingSphere(psphere, positions, count, positions_stride, radii ? radii : &rzero, radii ? radii_stride : 0);
+
+	bounds.center[0] = psphere[0];
+	bounds.center[1] = psphere[1];
+	bounds.center[2] = psphere[2];
+	bounds.radius = psphere[3];
+
+	return bounds;
+}
+
 void meshopt_optimizeMeshlet(unsigned int* meshlet_vertices, unsigned char* meshlet_triangles, size_t triangle_count, size_t vertex_count)
 {
 	using namespace meshopt;

src/meshoptimizer.h

@@ -434,7 +434,7 @@ MESHOPTIMIZER_EXPERIMENTAL size_t meshopt_simplifySloppy(unsigned int* destinati
  *
  * destination must contain enough space for the target index buffer (target_vertex_count elements)
  * vertex_positions should have float3 position in the first 12 bytes of each vertex
- * vertex_colors should can be NULL; when it's not NULL, it should have float3 color in the first 12 bytes of each vertex
+ * vertex_colors can be NULL; when it's not NULL, it should have float3 color in the first 12 bytes of each vertex
  * color_weight determines relative priority of color wrt position; 1.0 is a safe default
  */
 MESHOPTIMIZER_API size_t meshopt_simplifyPoints(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, const float* vertex_colors, size_t vertex_colors_stride, float color_weight, size_t target_vertex_count);
@@ -619,6 +619,15 @@ struct meshopt_Bounds
 MESHOPTIMIZER_API struct meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride);
 MESHOPTIMIZER_API struct 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);
 
+/**
+ * Experimental: Sphere bounds generator
+ * Creates bounding sphere around a set of points or a set of spheres; returns the center and radius of the sphere, with other fields of the result set to 0.
+ *
+ * positions should have float3 position in the first 12 bytes of each element
+ * radii can be NULL; when it's not NULL, it should have a non-negative float radius in the first 4 bytes of each element
+ */
+MESHOPTIMIZER_EXPERIMENTAL struct meshopt_Bounds meshopt_computeSphereBounds(const float* positions, size_t count, size_t positions_stride, const float* radii, size_t radii_stride);
+
 /**
  * Experimental: Cluster partitioner
  * Partitions clusters into groups of similar size, prioritizing grouping clusters that share vertices.