This page is in response to a message that was posted on the sci.math.num-analysis news group. The post asked if there was an algorithm that would calculate a minimum bounding circle which would enclose an arbitrary set of circles of different radii.
I thought this was an interesting problem, so I spent an afternoon coming up with an algorithm, and writing a small Mac program that would create a set of random circles, and then calculate such a minimum bounding circle for that set. The results of that program are displayed below for various circle sets.
An analytic solution for this problem is pretty difficult. It is fairly easy to create a minimum bounding circle for a point set, because there are an infinite number of circles which intersect 2 non-coincident points, and such circles all lie on the straight line between the 2 circle centers. Finding a circle which intersects more than 2 points thus involves finding the intersection of several straight lines.
There are also an infinite number of circles which enclose (and, in this case, osculate (i.e. are tangent) with) any two non-coincident circles, but this family lies on a straight line only if the two circles have equal radii. If the radii are not equal, the family of osculating circles lies on a curve (a parabola, I think), so that finding a circle which encloses several circles of different radii involves solving for the intersection(s) of several parabolas. Not a very inviting task!
The algorithm I use involves an interative search, somewhat like a gradient descent method for solving minimization problems. This was much easier to write, and does not take much time to run, even for a large number of circles. The algorithm begins with a "guess" circle which is centered on one of the input circles, and which has a minimum radius to enclose all of the circles.
At each step in the iteration, the center of the bounding circle is moved in a
direction which will allow its radius to decrease. A grid is used to determine
the step direction and size, and this grid size is dynamically made finer and
finer whenever possible. The code for determining the direction
and step size is listed below. Eventually, the circle settles on the minimum
bound. In the case of just 2 circles, the center is along the line connecting the
centers of the two circles.
For three circles, the bounding circle can be made "tight" against each of them
(or it may osculate only two, with the third well inside).
The algorithm also works well for a larger number of circles, and does not
really need any more steps to accomplish its task.
Here is the C++ code for performing this calculation. This code was compiled and run on a 68040-Mac with the Symantec C++ 7.0 compiler and project manager.
typedef struct circle_type
{
long double x; // x coordinate of center
long double y; // y coordinate of center
long double r; // radius
} circle_type;
class mainwin_type: public win_type
{
public:
// ...other fields
int circles; // number of circles
circle_type **circle; // handle to circle array
circle_type guess; // current bounding circle
long double grid_size; // calculation grid (initially 100.0)
long double step_size; // percentage of grid to step each iteration (0.5)
long double stop_size; // stop iterating when grid gets small (1.0)
int count; // number of iterations
};
void mainwin_type::new_set(void)
{
// create circles with random ceters and radii
for (int i = 0; i < circles; i ++)
{
(*circle)[i].x = 75 + 150.0 * rand() / 32767.0;
(*circle)[i].y = 75 + 150.0 * rand() / 32767.0;
(*circle)[i].r = 50.0 * rand() / 32767.0;
}
// start bound with first circle
guess.x = (*circle)[0].x;
guess.y = (*circle)[0].y;
guess.r = min_rad(guess.x, guess.y);
// initialize calculation params
grid_size = GRID_SIZE; // 100.0
step_size = STEP_SIZE; // 0.5
stop_size = STOP_SIZE; // 1.0
count = 0;
refresh(); // redraw window
}
void mainwin_type::recalc(void)
{
count = 0; // initialize count to zero
// stop iterating if grid size gets small enough, or if too many steps
// have gone by
while (count < 100 && grid_size > stop_size)
{
// find minimum bounding circles for 4 points near current guess:
// (1) x + grid_size
// (2) x - grid_size
// (3) y + grid_size
// (4) y - grid_size
long double r1 = min_rad(guess.x + grid_size, guess.y);
long double r2 = min_rad(guess.x - grid_size, guess.y);
long double r3 = min_rad(guess.x, guess.y + grid_size);
long double r4 = min_rad(guess.x, guess.y - grid_size);
// if any new circle radius is smaller than current one, take a
// step in that direction, proportional to current grid size
if (r1 < guess.r) guess.x += step_size * grid_size;
else if (r2 < guess.r) guess.x -= step_size * grid_size;
else if (r3 < guess.r) guess.y += step_size * grid_size;
else if (r4 < guess.r) guess.y -= step_size * grid_size;
// if none of the 4 directions gives a smaller circle (which will happen),
// look for alternate directions to go. in this case, potential alternate
// directions are perpendicular to a line between any given circle center
// and the center of the current bounding circle. (actually, this is
// wasteful, since only those circles which are currently on the periphery
// of the set need to be checked.)
else
{
int found = 0; // reset flag
// look for a new direction tangent to all circles
for (int i = 0; i < circles; i ++)
{
// compute vector between circle center and bounding circle center
long double dx = (*circle)[i].x - guess.x;
long double dy = (*circle)[i].y - guess.y;
long double len = sqrt(dx * dx + dy * dy);
// check only circles not coincident with bounding circle
if (len > 0)
{
// compute vector perpendicular to connecting line
long double dx2 = dy / len;
long double dy2 = -dx / len;
// compute minimum bounding radius for points in both directions along
// tangent
long double r5 = min_rad(guess.x + dx2 * grid_size, guess.y + dy2 * grid_size);
long double r6 = min_rad(guess.x - dx2 * grid_size, guess.y - dy2 * grid_size);
// if either direction produces a decrease in bounding radius, take a
// step in that direction
if (r5 < guess.r)
{
guess.x += dx2 * step_size * grid_size;
guess.y += dy2 * step_size * grid_size;
found = 1;
break;
}
else if (r6 < guess.r)
{
guess.x -= dx2 * step_size * grid_size;
guess.y -= dy2 * step_size * grid_size;
found = 1;
break;
}
}
}
// if no direction decreases bounding radius, decrease grid size and try
// again next time
if (!found) grid_size *= step_size;
}
// update bounding radius based on new circle coords
guess.r = min_rad(guess.x, guess.y);
count ++; // update count
refresh(); // redraw window
update(); // redraw window
}
}
long double mainwin_type::min_rad(long double x, long double y)
{
long double rad = 0; // initialize radius to zero
// check distance to all circles
for (int i = 0; i < circles; i ++)
{
// compute distance between point and circle center
long double dx = x - (*circle)[i].x;
long double dy = y - (*circle)[i].y;
// add circle radius to distance
long double rad2 = sqrt(dx * dx + dy * dy) + (*circle)[i].r;
// compute maximum value
if (rad2 > rad) rad = rad2;
}
return rad;
}
I make no claim that this is the best strategy to use for this problem, or that this particular implementation is an efficient one. However, the algorithm seems to work for all cases, and is fairly fast for a large (~100) number of circles.
Download the program file: circles.sea.hqx; it is about 88K.
This is a binhex of a self-extracting archive for the Macintosh. You will need to unbinhex it once you download it (or your web client may do this for you). Just double-click on the resulting application to unstuff the program. There is a "READ-ME" file in the folder; please read it.