Title: Find a minimal bounding circle of a set of points in C#

The example Find the convex hull of a set of points in C# finds the convex hull of a set of points. A convex hull is the smallest convex polygon that encloses the points. This example extends that result to find a minimal circle enclosing the points.

The key is to note that a minimal bounding circle passes through either two or three of the convex hull's points. The following picture shows the two possible scenarios. In the pictures, blue shows the convex hull, red shows a culling trapezoid, and orange shows a culling rectangle, all as described in the previous post.

The following code shows how this example finds minimal bounding circles.

// Find a minimal bounding circle. public static void FindMinimalBoundingCircle( List<PointF> points, out PointF center, out float radius) { // Find the convex hull. List<PointF> hull = MakeConvexHull(points); // The best solution so far. PointF best_center = points[0]; float best_radius2 = float.MaxValue; // Look at pairs of hull points. for (int i = 0; i < hull.Count - 1; i++) { for (int j = i + 1; j < hull.Count; j++) { // Find the circle through these two points. PointF test_center = new PointF( (hull[i].X + hull[j].X) / 2f, (hull[i].Y + hull[j].Y) / 2f); float dx = test_center.X - hull[i].X; float dy = test_center.Y - hull[i].Y; float test_radius2 = dx * dx + dy * dy; // See if this circle would be an improvement. if (test_radius2 < best_radius2) { // See if this circle encloses all of the points. if (CircleEnclosesPoints(test_center, test_radius2, hull, i, j, -1)) { // Save this solution. best_center = test_center; best_radius2 = test_radius2; } } } // for i } // for j // Look at triples of hull points. for (int i = 0; i < hull.Count - 2; i++) { for (int j = i + 1; j < hull.Count - 1; j++) { for (int k = j + 1; k < hull.Count; k++) { // Find the circle through these three points. PointF test_center; float test_radius2; FindCircle(hull[i], hull[j], hull[k], out test_center, out test_radius2); // See if this circle would be an improvement. if (test_radius2 < best_radius2) { // See if this circle encloses all the points. if (CircleEnclosesPoints(test_center, test_radius2, hull, i, j, k)) { // Save this solution. best_center = test_center; best_radius2 = test_radius2; } } } // for k } // for i } // for j center = best_center; if (best_radius2 == float.MaxValue) radius = 0; else radius = (float)Math.Sqrt(best_radius2); }

The code first uses the technique described in the previous post to find the convex hull.

It then loops through every pair of points on the hull to see if they lie on a bounding circle. For each pair of points, the program tests the circle with center exactly halfway between the two points and passing through the points. If the circle's radius squared is smaller than the best value found so far, the program calls the CircleEnclosesPoints method (described shortly) to see if the circle encloses all of the points. If the circle does enclose the points, the program updates its best circle center and radius.

After checking all pairs of points, the program loops through all triples of points. For each triple, the program uses the technique described in the post Draw a circle through three points in C# to get a circle passing through the three points. It compares the circle's radius squared to the best so far and calls CircleEnclosesPoints as before to see if it should update the best circle.

When it has finished checking all of the triples of points, the code compares best_radius2 to float.MaxValue to see if it found a circle. If the values are the same, that means the points array holds a single point. In that case, the program sets radius to 0 so it returns a circle centered at the single point with radius 0.

If best_radius2 doesn't equal float.MaxValue, the program sets the return radius result and ends.

The following code shows the CircleEnclosesPoints method.

// Return true if the indicated circle encloses all of the points. private static bool CircleEnclosesPoints(PointF center, float radius2, List<PointF> points, int skip1, int skip2, int skip3) { for (int i = 0; i < points.Count; i++) { if ((i != skip1) && (i != skip2) && (i != skip3)) { PointF point = points[i]; float dx = center.X - point.X; float dy = center.Y - point.Y; float test_radius2 = dx * dx + dy * dy; if (test_radius2 > radius2) return false; } } return true; }

This method takes as parameters a circle's center and radius squared, the list of points to examine, and three points that lie on the circle. It loops through the list of points, skipping the three on the circle, and determines whether the other points are all inside the circle. The code skips the three on the circle so rounding errors don't incorrectly make it seem like those points are not within the circle.

Because this method examines all triples of the points in the convex hull, it has run time O(N³) where N is the number of points in the convex hull. There are faster algorithms, but for most "typical" applications, the number of points in the convex hull isn't huge so this is fast enough. (It's also much simpler than the faster algorithms.)

Download the example to experiment with it and to see additional details.