Minimum Bounding Circles

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.

### Here is the starting configuration for 2 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.

### Here is the final configuration for 2 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).

### Here is the final configuration for 3 circles: The algorithm also works well for a larger number of circles, and does not really need any more steps to accomplish its task.

### Here are the final configurations for 10 and 100 circles:  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.

### First, a struct for a circle type is defined:

```typedef struct circle_type
{
long double  x;	// x coordinate of center
long double  y;	// y coordinate of center
} circle_type;
```

### Here is part of the window class which contains the random circle set and the final bounding circle:

```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
};
```

### This method creates a new random circle set and initializes the calculation values:

```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).x;
guess.y = (*circle).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
}
```

### Here is the method that does the actual calculations:

```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

count ++;	// update count

refresh();	// redraw window
update(); 	// redraw window
}
}
```

### Here is the method that computes a minimum bounding circle centered on given x and y coordinates:

```long double mainwin_type::min_rad(long double x, long double y)
{

// 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;

long double rad2 = sqrt(dx * dx + dy * dy) + (*circle)[i].r;

// compute maximum value

}

}

```

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.