I would hold off on replacing it. The implementation I currently replaced it with is fast (O(n)) but again it does not always generate the optimal circle. For example, when generating a circle for the triangle-'sh shape in the testbed, it generates an oversized circle.

It looks like the current code is either Welzl's algorithm or a variant of it. I'm pretty sure I'll be able to figure out what's going on and either fix it or replace it with another implementation of it. I just need time to really look at it (which I haven't had). If nobody else is running into problems but me, I'd just hold off.

Though, you are completely welcome to the existing code if you want, it's a direct port of what's on that page I posted above:

Just replace the body of calculate_minimum_enclosing_disc(...) with this:

Code:

// The following is a port of the code listed at http://softsurfer.com/Archive/algorithm_0107/algorithm_0107.htm#fastBall()
// which falls under the below copyright.
//
// Copyright 2001, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
CL_Pointf C;
float rad;
float rad22;
float xmin, xmax, ymin, ymax;
int Pxmin, Pxmax, Pymin, Pymax;
xmin = xmax = points[0].x;
ymin = ymax = points[0].y;
Pxmin = Pxmax = Pymin = Pymax = 0;
for(int i=start+1; i<=end; i++) {
if(points[i].x < xmin) {
xmin = points[i].x;
Pxmin = i;
}
else if(points[i].x > xmax) {
xmax = points[i].x;
Pxmax = i;
}
if(points[i].y < ymin) {
ymin = points[i].y;
Pymin = i;
}
else if(points[i].y > ymax) {
ymax = points[i].y;
Pymax = i;
}
}
CL_Pointf dVx = points[Pxmax] - points[Pxmin];
CL_Pointf dVy = points[Pymax] - points[Pymin];
float dx2 = dVx.dot(dVx);
float dy2 = dVy.dot(dVy);
if(dx2>=dy2) {
C = points[Pxmin] + (dVx / 2.0f);
CL_Pointf tmp = points[Pxmax] - C;
rad22 = tmp.dot(tmp);
}
else {
C = points[Pymin] + (dVy / 2.0f);
CL_Pointf tmp = points[Pymax] - C;
rad22 = tmp.dot(tmp);
}
rad = sqrtf(rad22);
CL_Pointf dV;
float dist, dist2;
for(int i=start; i<=end; i++) {
dV = points[i] - C;
dist2 = dV.dot(dV);
if(dist2<=rad22)
continue;
dist = sqrtf(dist2);
rad = (rad + dist) / 2.0f;
rad22 = rad * rad;
C = C + (dV * ((dist-rad)/dist));
}
smalldisc.position = C;
smalldisc.radius = rad;

## Bookmarks