Minimum enclosing circle using Welzl’s algorithm

Last Updated : 17 Feb, 2025

Given an array arr[][] containing n points in a 2-D plane with integer coordinates. The task is to find the center and the radius of the minimum enclosing circle (MEC).

Note: A minimum enclosing circle is a circle in which all the points lie either inside the circle or on its boundaries.
Prerequisite: Equation of circle when three points on the circle are given

Examples:

Input: arr[][] = [[0, 0], [0, 1], [1, 0]]
Output: 0.5 0.5 0.7071
Explanation: The first two values are the co-ordinates of the center of circle, and the last value is radius of the circle.
On plotting the above circle with radius 0.707 and center (0.5, 0.5), it can be observed clearly that all the mentioned points lie either inside or on the circle.

Input: arr[][] = [[5, -2], [-3, -2], [-2, 5], [1, 6], [0, 2]]
Output: 1.0 1.0 5.0

Note: In the article Minimum Enclosing Circle, a naive approach and an optimized approach by getting the convex full of the set of points first and then performing a naive approach has been discussed. Although the optimized solution would work very well for certain inputs, the worst-case time complexity after that optimization was still O(n⁴). In this article, an optimized approach has been discussed.

Minimum Enclosing Circle using Welzl’s Algorithm

The idea is to use Welzl’s recursive algorithm. Using this algorithm, the MEC can be found in linear time. The working of the algorithm depends on the observations and conclusions drawn from the previous article. The idea of the algorithm is to randomly remove a point from the given input set to form a circle equation. Once the equation is formed, check if the point which was removed is enclosed by the equation or not. If it doesn’t, then the point must lie on the boundary of the MEC. Therefore, this point is considered as a boundary point and the function is recursively called.

The detailed working of the algorithm is as follows:

The algorithm takes a set of points P and a set R that’s initially empty and used to represent the points on the boundary of the MEC as the input.
The base case of the algorithm is when P becomes empty or the size of the set R is equal to 3:

If P is empty, then all the points have been processed.
If |R| = 3, then 3 points have already been found that lie on the circle boundary, and since a circle can be uniquely determined using 3 points only, the recursion can be stopped.

When the algorithm reaches the base case above, it returns the trivial solution for R, being:

If |R| = 1, we return a circle centered at R[0] with radius = 0
If |R| = 2, we return the MEC for R[0] and R[2]
If |R| = 3, we return the MEC by trying the 3 pairs (R[0], R[1]), (R[0], R[2]), (R[1], R[2])
If none of these pairs is valid, we return the circle defined by the 3 points in R

If the base case is not reached yet, we do the following:

Pick a random point p from P and remove it from P
Call the algorithm on P and R to get circle d
If p is enclosed by d, then we return d
otherwise, p must lie on the boundary of the MEC
- Add p to R
- Return the output of the algorithm on P and R

CPP

// c++ program to find the minimum enclosing // circle for N integer points in a 2-d plane #include <bits/stdc++.h> using namespace std;  // Structure to represent a 2D point struct Point {     double x, y; };  // Structure to represent a 2D circle struct Circle {     Point c;     double r; };  // Function to return the euclidean distance // between two points double dist(const Point& a, const Point& b) {     return sqrt(pow(a.x - b.x, 2)                 + pow(a.y - b.y, 2)); }  // Function to check whether a point lies inside // or on the boundaries of the circle bool isInside(const Circle& c, const Point& p) {     return dist(c.c, p) <= c.r; }  // The following two functions are used // To find the equation of the circle when // three points are given.  // Helper method to get a circle defined by 3 points Point getCircleCenter(double bx, double by,                         double cx, double cy) {     double b = bx * bx + by * by;     double c = cx * cx + cy * cy;     double d = bx * cy - by * cx;     return { (cy * b - by * c) / (2 * d),              (bx * c - cx * b) / (2 * d) }; }  // Function to return a unique circle that // intersects three points Circle circleFrom(const Point& a, const Point& b,                    const Point& c) {     Point i = getCircleCenter(b.x - a.x, b.y - a.y,                                 c.x - a.x, c.y - a.y);      i.x += a.x;     i.y += a.y;     return { i, dist(i, a) }; }  // Function to return the smallest circle // that intersects 2 points Circle circleFrom(const Point& a, const Point& b) {      // Set the center to be the midpoint of a and b     Point c = { (a.x + b.x) / 2.0, (a.y + b.y) / 2.0 };      // Set the radius to be half the distance AB     return { c, dist(a, b) / 2.0 }; }  // Function to check whether a circle // encloses the given points bool isValidCircle(const Circle& c,                      const vector<Point>& p) {      // Iterating through all the points     // to check  whether the points     // lie inside the circle or not     for (const Point& i : p)         if (!isInside(c, i))             return false;     return true; }  // Function to return the minimum enclosing // circle for N <= 3 Circle minCircleTrivial(vector<Point>& p) {     assert(p.size() <= 3);     if (p.empty()) {         return { { 0, 0 }, 0 };     }     else if (p.size() == 1) {         return { p[0], 0 };     }     else if (p.size() == 2) {         return circleFrom(p[0], p[1]);     }      // To check if MEC can be determined     // by 2 points only     for (int i = 0; i < 3; i++) {         for (int j = i + 1; j < 3; j++) {              Circle c = circleFrom(p[i], p[j]);             if (isValidCircle(c, p))                 return c;         }     }     return circleFrom(p[0], p[1], p[2]); }  // Returns the MEC using Welzl's algorithm // Takes a set of input points p and a set r // points on the circle boundary. // n represents the number of points in p // that are not yet processed. Circle welzlHelper(vector<Point>& p,                     vector<Point> r, int n) {      // Base case when all points processed or |r| = 3     if (n == 0 || r.size() == 3) {         return minCircleTrivial(r);     }      // Pick a random point randomly     int idx = rand() % n;     Point pnt = p[idx];      // Put the picked point at the end of p     // since it's more efficient than     // deleting from the middle of the vector     swap(p[idx], p[n - 1]);      // Get the MEC circle d from the     // set of points p - {p}     Circle d = welzlHelper(p, r, n - 1);      // If d contains pnt, return d     if (isInside(d, pnt)) {         return d;     }      // Otherwise, must be on the boundary of the MEC     r.push_back(pnt);      // Return the MEC for p - {p} and r U {p}     return welzlHelper(p, r, n - 1); }  Circle welzl(const vector<Point>& p) {      vector<Point> pCopy = p;     random_shuffle(pCopy.begin(), pCopy.end());     return welzlHelper(pCopy, {}, pCopy.size()); }  int main() {      Circle mec = welzl({ { 5, -2 },                           { -3, -2 },                           { -2, 5 },                           { 1, 6 },                           { 0, 2 } });     cout << mec.c.x << " " << mec.c.y <<" "<< mec.r;     return 0; }

Java

import java.util.*;  class GFG {      static class Point {         double x, y;         Point(double x, double y) { this.x = x; this.y = y; }     }      static class Circle {         Point c;         double r;         Circle(Point c, double r) { this.c = c; this.r = r; }     }      static double dist(Point a, Point b) {         return Math.sqrt(Math.pow(a.x - b.x, 2) + Math.pow(a.y - b.y, 2));     }      static boolean isInside(Circle c, Point p) {         return dist(c.c, p) <= c.r;     }      static Point getCircleCenter(double bx, double by, double cx, double cy) {         double b = bx * bx + by * by;         double c = cx * cx + cy * cy;         double d = bx * cy - by * cx;         return new Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d));     }      static Circle circleFrom(Point a, Point b, Point c) {         Point i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);         i.x += a.x; i.y += a.y;         return new Circle(i, dist(i, a));     }      static Circle circleFrom(Point a, Point b) {         Point c = new Point((a.x + b.x) / 2.0, (a.y + b.y) / 2.0);         return new Circle(c, dist(a, b) / 2.0);     }      static boolean isValidCircle(Circle c, List<Point> p) {         for (Point i : p) if (!isInside(c, i)) return false;         return true;     }      static Circle minCircleTrivial(List<Point> p) {         assert p.size() <= 3;         if (p.isEmpty()) return new Circle(new Point(0, 0), 0);         if (p.size() == 1) return new Circle(p.get(0), 0);         if (p.size() == 2) return circleFrom(p.get(0), p.get(1));          for (int i = 0; i < 3; i++) {             for (int j = i + 1; j < 3; j++) {                 Circle c = circleFrom(p.get(i), p.get(j));                 if (isValidCircle(c, p)) return c;             }         }         return circleFrom(p.get(0), p.get(1), p.get(2));     }      static Circle welzlHelper(List<Point> p, List<Point> r, int n) {         if (n == 0 || r.size() == 3) return minCircleTrivial(new ArrayList<>(r));          int idx = new Random().nextInt(n);         Point pnt = p.get(idx);         Collections.swap(p, idx, n - 1);          Circle d = welzlHelper(p, r, n - 1);         if (isInside(d, pnt)) return d;          List<Point> newR = new ArrayList<>(r);         newR.add(pnt);         return welzlHelper(p, newR, n - 1);     }      static Circle welzl(List<Point> p) {         List<Point> pCopy = new ArrayList<>(p);         Collections.shuffle(pCopy);         return welzlHelper(pCopy, new ArrayList<>(), pCopy.size());     }      public static void main(String[] args) {         Circle mec = welzl(Arrays.asList(             new Point(5, -2), new Point(-3, -2), new Point(-2, 5),             new Point(1, 6), new Point(0, 2)         ));          System.out.println(mec.c.x + " " + mec.c.y + " " + mec.r);     } }

Python

import math import random  # Structure to represent a 2D point class Point:     def __init__(self, x, y):         self.x = x         self.y = y  # Structure to represent a 2D circle class Circle:     def __init__(self, c, r):         self.c = c         self.r = r  # Function to return the euclidean distance between two points def dist(a, b):     return math.sqrt((a.x - b.x) ** 2 + (a.y - b.y) ** 2)  # Function to check whether a point lies inside or on the boundaries of the circle def isInside(c, p):     return dist(c.c, p) <= c.r  # Helper method to get a circle defined by 3 points def getCircleCenter(bx, by, cx, cy):     b = bx * bx + by * by     c = cx * cx + cy * cy     d = bx * cy - by * cx     return Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d))  # Function to return a unique circle that intersects three points def circleFrom(a, b, c):     i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)     i.x += a.x     i.y += a.y     return Circle(i, dist(i, a))  # Function to return the smallest circle that intersects 2 points def circleFromTwo(a, b):     c = Point((a.x + b.x) / 2.0, (a.y + b.y) / 2.0)     return Circle(c, dist(a, b) / 2.0)  # Function to check whether a circle encloses the given points def isValidCircle(c, p):     return all(isInside(c, point) for point in p)  # Function to return the minimum enclosing circle for N <= 3 def minCircleTrivial(p):     assert len(p) <= 3     if not p:         return Circle(Point(0, 0), 0)     elif len(p) == 1:         return Circle(p[0], 0)     elif len(p) == 2:         return circleFromTwo(p[0], p[1])      for i in range(3):         for j in range(i + 1, 3):             c = circleFromTwo(p[i], p[j])             if isValidCircle(c, p):                 return c     return circleFrom(p[0], p[1], p[2])  # Returns the MEC using Welzl's algorithm def welzlHelper(p, r, n):     if n == 0 or len(r) == 3:         return minCircleTrivial(r[:])  # Ensure we pass a copy      idx = random.randint(0, n - 1)     pnt = p[idx]     p[idx], p[n - 1] = p[n - 1], p[idx]      d = welzlHelper(p, r, n - 1)      if isInside(d, pnt):         return d      return welzlHelper(p, r + [pnt], n - 1)  # Append without modifying original r  def welzl(p):     pCopy = list(p)     random.shuffle(pCopy)     return welzlHelper(pCopy, [], len(pCopy))  # Main function def main():     test_points = [         Point(5, -2),         Point(-3, -2),         Point(-2, 5),         Point(1, 6),         Point(0, 2)     ]      mec = welzl(test_points)     print(mec.c.x, mec.c.y, mec.r)  # Run the main function if __name__ == "__main__":     main()

// C# program to find the minimum enclosing  // circle for N integer points in a 2-d plane using System; using System.Collections.Generic;  class GFG {      // Structure to represent a 2D point     class Point {         public double x, y;          public Point(double x, double y) {             this.x = x;             this.y = y;         }     }      // Structure to represent a 2D circle     class Circle {         public Point c;         public double r;          public Circle(Point c, double r) {             this.c = c;             this.r = r;         }     }      // Function to return the euclidean distance     // between two points     static double Dist(Point a, Point b) {         return Math.Sqrt(Math.Pow(a.x - b.x, 2)                         + Math.Pow(a.y - b.y, 2));     }      // Function to check whether a point lies inside     // or on the boundaries of the circle     static bool IsInside(Circle c, Point p) {         return Dist(c.c, p) <= c.r;     }      // The following two functions are used     // To find the equation of the circle when     // three points are given.      // Helper method to get a circle defined by 3 points     static Point GetCircleCenter(double bx, double by, double cx, double cy) {         double b = bx * bx + by * by;         double c = cx * cx + cy * cy;         double d = bx * cy - by * cx;         return new Point((cy * b - by * c) / (2 * d),                          (bx * c - cx * b) / (2 * d));     }      // Function to return a unique circle that     // intersects three points     static Circle CircleFrom(Point a, Point b, Point c) {         Point i = GetCircleCenter(b.x - a.x, b.y - a.y,                                   c.x - a.x, c.y - a.y);          i.x += a.x;         i.y += a.y;         return new Circle(i, Dist(i, a));     }      public static void Main() {         Circle mec = CircleFrom(new Point(5, -2),                                 new Point(-3, -2),                                 new Point(-2, 5));          Console.WriteLine(mec.c.x + " " + mec.c.y + " " + mec.r);     } }

JavaScript

class Point {     constructor(x, y) { this.x = x; this.y = y; } }  class Circle {     constructor(c, r) { this.c = c; this.r = r; } }  function dist(a, b) {     return Math.sqrt((a.x - b.x) ** 2 + (a.y - b.y) ** 2); }  function isInside(c, p) {     return dist(c.c, p) <= c.r; }  function getCircleCenter(bx, by, cx, cy) {     let b = bx * bx + by * by;     let c = cx * cx + cy * cy;     let d = bx * cy - by * cx;     return new Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d)); }  function circleFrom(a, b, c) {     let i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);     i.x += a.x; i.y += a.y;     return new Circle(i, dist(i, a)); }  function circleFromTwo(a, b) {     let c = new Point((a.x + b.x) / 2.0, (a.y + b.y) / 2.0);     return new Circle(c, dist(a, b) / 2.0); }  function isValidCircle(c, p) {     return p.every(point => isInside(c, point)); }  function minCircleTrivial(p) {     if (p.length === 0) return new Circle(new Point(0, 0), 0);     if (p.length === 1) return new Circle(p[0], 0);     if (p.length === 2) return circleFromTwo(p[0], p[1]);      for (let i = 0; i < 3; i++) {         for (let j = i + 1; j < 3; j++) {             let c = circleFromTwo(p[i], p[j]);             if (isValidCircle(c, p)) return c;         }     }     return circleFrom(p[0], p[1], p[2]); }  function welzlHelper(p, r, n) {     if (n === 0 || r.length === 3) return minCircleTrivial([...r]);      let idx = Math.floor(Math.random() * n);     let pnt = p[idx];     [p[idx], p[n - 1]] = [p[n - 1], p[idx]];      let d = welzlHelper(p, r, n - 1);     if (isInside(d, pnt)) return d;      return welzlHelper(p, [...r, pnt], n - 1); }  function welzl(p) {     let pCopy = [...p];     pCopy.sort(() => Math.random() - 0.5);     return welzlHelper(pCopy, [], pCopy.length); }  let mec = welzl([     new Point(5, -2), new Point(-3, -2), new Point(-2, 5),     new Point(1, 6), new Point(0, 2) ]);  console.log(mec.c.x, mec.c.y, mec.r);

Output

1 1 5

Time Complexity: O(n), where n is the number of points. With every call, the size of P gets reduced by one. Also, the size of R can remain the same or can be increased by one. Since |R| cannot exceed 3, then the number of different states would be 3N. Therefore, this makes the time complexity to be O(N).
Auxiliary Space: O(n), considering the recursive call stack.

Minimum rotations to unlock a circular lock

marwan tammam

Improve

Article Tags :

Practice Tags :

Minimum enclosing circle using Welzl’s algorithm

Minimum Enclosing Circle using Welzl’s Algorithm

Similar Reads