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^{3}) 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.)

Dear Rod,

first, many thanks for prividing this very useful algorithms!

In your example program the minimal bounding circle works well. But after I used it in another context, I get wrong results. This was, because in the Geometry.cs in the downloadable file, there is the loop

if (CircleEnclosesPoints(test_center, test_radius2, points, i, j, -1))

in the FindMinimalBoundingCircle method. It uses “points” instead of “hull”. In your text above it’s correct!

Best regards

Reinhard

Hi Reinhard,

It shouldn’t matter whether you call CircleEnclosesPoints for all of the points or just those in the hull. If the circle encloses one, then it should enclose the other. But if you have a set of points that the program can’t handle, please let me know. If you tell me the points, I’ll look into it and see what’s wrong.

But checking the hull will be faster than checking every point, so I’ve updated the example to do that. Excellent suggestion!

Rod

You claim that your minimal bounding circle algorithm runs in O(N^3) time, but actually it runs in O(N^4) time. This is because you have a triple for-loop containing a call to CircleEnclosesPoints(), which has its own loop.

If you’re interested in reading, I implemented a solution to this problem in O(N) linear time, with about 3 times as much code as yours: https://www.nayuki.io/page/smallest-enclosing-circle

Oops. I think you’re right about the O(N^4). It does get some speed up by culling, but not if all of the points lie on the convex hull.