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.

[trifonov-borislav]

In some cases I was getting radius up to 20% larger than the optimal using Ritter, but the average case is not much worse. I implemented EPOS-6 before and with and works well (within a couple %) while it's still fast enough, but it relies on Welzl as a subroutine: basically you run Welzl on the projected by EPOS point subset. The trick is that Welzl has numerical issues and can explode, so you need to run Ritter beforehand and pass that in as an upper bound to detect such a situation (simply checking determinant is not sufficient). Definitely adds complexity but it's not too bad

[zeux]

FWIW I continue to have reservations about using Welzl but I've adjusted the code to use the axes EPOS uses (EPOS-14 variant to be more precise; I've tested EPOS-26 axis selection and the extra axes introduced there don't seem to meaningfully improve the radius on clusters I've tested) in #883.

On some meshes this results in 1-2% tighter average cluster bounds - I haven't measured the impact on worst case cluster bounds but I'd expect it to be more significant obviously. As a result this improves the cluster occlusion culling efficiency a little bit; the time to compute the bounds is a little longer but it's mitigated by small code tweaks and only using this for positions, as I didn't find the effect on cone angles to be significant.