Minimum bounding sphere computation for meshes/meshlets

[expenses]

Hi, I've been using your library to compute meshlets and their bounds and it's been working out well so far. The current code for this uses a linear scan:

meshoptimizer/src/clusterizer.cpp

Lines 78 to 146 in 28bd987

	static void computeBoundingSphere(float result[4], const float points[][3], size_t count)
	{
	assert(count > 0);
	// 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];
	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];
	}
	}
	// find the pair of points with largest distance
	float paxisd2 = 0;
	int paxis = 0;
	for (int axis = 0; axis < 3; ++axis)
	{
	const float* p1 = points[pmin[axis]];
	const float* p2 = points[pmax[axis]];
	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]);
	if (d2 > paxisd2)
	{
	paxisd2 = d2;
	paxis = axis;
	}
	}
	// use the longest segment as the initial sphere diameter
	const float* p1 = points[pmin[paxis]];
	const float* p2 = points[pmax[paxis]];
	float center[3] = {(p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2, (p1[2] + p2[2]) / 2};
	float radius = sqrtf(paxisd2) / 2;
	// iteratively adjust the sphere up until all points fit
	for (size_t i = 0; i < count; ++i)
	{
	const float* p = points[i];
	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]);
	if (d2 > radius * radius)
	{
	float d = sqrtf(d2);
	assert(d > 0);
	float k = 0.5f + (radius / d) / 2;
	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;
	}
	}
	result[0] = center[0];
	result[1] = center[1];
	result[2] = center[2];
	result[3] = radius;
	}

. This is good enough for just-in-time uses, but for offline computation it'd be great to use an implementation of Welzl's algorithm that calculates the exact minimum bounding sphere.

Something like gltfpack could compute and store the minimum bounding-sphere for meshes as well.

[zeux]

I'd need to go back to look at my notes; my recollection is that I considered Welzl and rejected it because of stability and performance/determinism concerns (it's amortized linear in theory but theory and practice tend to diverge); note that the currently used algorithm doesn't guarantee minimal results but it provides, generally, results that are fairly good. I think I experimented with other guaranteed-linear-time algorithms and this was the best among the ones that are based on heuristics.

Out of curiosity, do you have specific examples where the spheres are significantly suboptimal? I'm very open to improving this but I'd prefer to explore options that aren't Welzl.

[castano]

There's also Ritter's and EPOS algorithms: https://ep.liu.se/ecp/034/009/ecp083409.pdf
Here's an implementation of all of them: https://gist.github.com/castano/893695b3291d949242f2aee3943dfd55
Note there are some robustness issues in the Sphere constructor when the points are coplanar. I've been meaning to revisit and robustify that code, but haven't had a good excuse to do so.

[zeux]

Yeah the clusterizer code uses Ritter right now.

[expenses]

Out of curiosity, do you have specific examples where the spheres are significantly suboptimal?

No, I'm just a little wary that some performance is left on the table. @castano thanks for much for posting these! I was having trouble finding drop-in implementations. I'll do some testing with my own meshes and report back.