Minimum enclosing circle using Welzl’s algorithm

Prerequisites: Equation of circle when three points on the circle are given, Minimum Enclosing Circle.
Given an array arr[][] containing N points in a 2-D plane with integer coordinates. The task is to find the centre and the radius of the minimum enclosing circle(MEC). A minimum enclosing circle is a circle in which all the points lie either inside the circle or on its boundaries.
Examples:

Input: arr[][] = {{0, 0}, {0, 1}, {1, 0}} 
Output: Center = {0.5, 0.5}, Radius = 0.7071 
Explanation:
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: Center = {1.0, 1.0}, Radius = 5.000 

Approach: In the previous article, 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.
The idea is to use Welzl’s recursive algorithm. Using this algorithm, the MEC can be found in O(N). 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

Below is the implementation of the above approach:

CPP

// C++ program to find the minimum enclosing
// circle for N integer points in a 2-D plane
#include <algorithm>
#include <assert.h>
#include <iostream>
#include <math.h>
#include <vector>
using namespace std;

// Defining infinity
const double INF = 1e18;

// 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 is_inside(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 get_circle_center(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 circle_from(const Point& A, const Point& B,
                   const Point& C)
{
    Point I = get_circle_center(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 circle_from(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 is_valid_circle(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& p : P)
        if (!is_inside(c, p))
            return false;
    return true;
}

// Function to return the minimum enclosing
// circle for N <= 3
Circle min_circle_trivial(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 circle_from(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 = circle_from(P[i], P[j]);
            if (is_valid_circle(c, P))
                return c;
        }
    }
    return circle_from(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 welzl_helper(vector<Point>& P,
                    vector<Point> R, int n)
{
    // Base case when all points processed or |R| = 3
    if (n == 0 || R.size() == 3) {
        return min_circle_trivial(R);
    }

    // Pick a random point randomly
    int idx = rand() % n;
    Point p = 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 = welzl_helper(P, R, n - 1);

    // If d contains p, return d
    if (is_inside(d, p)) {
        return d;
    }

    // Otherwise, must be on the boundary of the MEC
    R.push_back(p);

    // Return the MEC for P - {p} and R U {p}
    return welzl_helper(P, R, n - 1);
}

Circle welzl(const vector<Point>& P)
{
    vector<Point> P_copy = P;
    random_shuffle(P_copy.begin(), P_copy.end());
    return welzl_helper(P_copy, {}, P_copy.size());
}

// Driver code
int main()
{
    Circle mec = welzl({ { 0, 0 },
                         { 0, 1 },
                         { 1, 0 } });
    cout << "Center = { " << mec.C.X << ", " << mec.C.Y
         << " } Radius = " << mec.R << endl;

    Circle mec2 = welzl({ { 5, -2 },
                          { -3, -2 },
                          { -2, 5 },
                          { 1, 6 },
                          { 0, 2 } });
    cout << "Center = { " << mec2.C.X << ", " << mec2.C.Y
         << " } Radius = " << mec2.R << endl;

    return 0;
}

Python3

# Python3 program to find the minimum enclosing
# circle for N integer points in a 2-D plane
from math import sqrt
from random import randint,shuffle

# Defining infinity
INF = 1e18

# Structure to represent a 2D point
class Point :
    def __init__(self,X=0,Y=0) -> None:
        self.X=X
        self.Y=Y

# Structure to represent a 2D circle
class Circle :
    def __init__(self,c=Point(),r=0) -> None:
        self.C=c
        self.R=r


# Function to return the euclidean distance
# between two points
def dist(a, 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
def is_inside(c, 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
def get_circle_center(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 the smallest circle
# that intersects 2 points
def circle_from1(A, B):
    # Set the center to be the midpoint of A and B
    C = Point((A.X + B.X) / 2.0, (A.Y + B.Y) / 2.0 )

    # Set the radius to be half the distance AB
    return Circle(C, dist(A, B) / 2.0 )

# Function to return a unique circle that
# intersects three points
def circle_from2(A, B, C):
    I = get_circle_center(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 check whether a circle
# encloses the given points
def is_valid_circle(c, P):

    # Iterating through all the points
    # to check  whether the points
    # lie inside the circle or not
    for p in P:
        if (not is_inside(c, p)):
            return False
    return True


# Function to return the minimum enclosing
# circle for N <= 3
def min_circle_trivial(P):
    assert(len(P) <= 3)
    if not P :
        return Circle()

    elif (len(P) == 1) :
        return Circle(P[0], 0)

    elif (len(P) == 2) :
        return circle_from1(P[0], P[1])


    # To check if MEC can be determined
    # by 2 points only
    for i in range(3):
        for j in range(i + 1,3):

            c = circle_from1(P[i], P[j])
            if (is_valid_circle(c, P)):
                return c


    return circle_from2(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.
def welzl_helper(P, R, n):
    # Base case when all points processed or |R| = 3
    if (n == 0 or len(R) == 3) :
        return min_circle_trivial(R)


    # Pick a random point randomly
    idx = randint(0,n-1)
    p = P[idx]

    # Put the picked point at the end of P
    # since it's more efficient than
    # deleting from the middle of the vector
    P[idx], P[n - 1]=P[n-1],P[idx]

    # Get the MEC circle d from the
    # set of points P - :p
    d = welzl_helper(P, R.copy(), n - 1)

    # If d contains p, return d
    if (is_inside(d, p)) :
        return d


    # Otherwise, must be on the boundary of the MEC
    R.append(p)

    # Return the MEC for P - :p and R U :p
    return welzl_helper(P, R.copy(), n - 1)


def welzl(P):
    P_copy = P.copy()
    shuffle(P_copy)
    return welzl_helper(P_copy, [], len(P_copy))


# Driver code
if __name__ == '__main__':
    mec = welzl([Point(0, 0) ,
                         Point(0, 1) ,
                         Point(1, 0)  ])
    print("Center = {",mec.C.X,",", mec.C.Y,"} Radius =",mec.R)

    mec2 = welzl([Point(5, -2) ,
                          Point(-3, -2) ,
                          Point(-2, 5) ,
                          Point(1, 6),
                          Point(0, 2)] )
    print("Center = {",mec2.C.X,",",mec2.C.Y,"} Radius =",mec2.R)

Javascript

// JS program to find the minimum enclosing
// circle for N integer points in a 2-D plane

// Function to get random int
function getRandomInt() {
  return Math.ceil(Math.random() * 32676);
}

let R = []

// Function to shuffle array
function shuffle(array) {
  let currentIndex = array.length,  randomIndex;

  // While there remain elements to shuffle.
  while (currentIndex != 0) {

    // Pick a remaining element.
    randomIndex = Math.floor(Math.random() * currentIndex);
    currentIndex--;

    // And swap it with the current element.
    [array[currentIndex], array[randomIndex]] = [
      array[randomIndex], array[currentIndex]];
  }

  return array;
}

// Defining infinity
let INF = 1e18;

// Structure to represent a 2D point
class Point {
    constructor(a = 0, b = 0)
    {
        this.X = a;
        this.Y = b;
    }
};

// Structure to represent a 2D circle
class Circle {
    constructor(p = new Point(0, 0), b = 0)
    {
        this.C = p;
        this.R = b;
    }
};



// Function to return the euclidean distance
// between two points
function dist(a, 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
function is_inside(c, 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
function get_circle_center( 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 to return a unique circle that
// intersects three points
function circle_from(A, B, C)
{
    let I = get_circle_center(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 to return the smallest circle
// that intersects 2 points
function circle_from(A, B)
{
    // Set the center to be the midpoint of A and B
    let C = new Point((A.X + B.X) / 2.0, (A.Y + B.Y) / 2.0);

    // Set the radius to be half the distance AB
    return new Circle(C, dist(A, B) / 2.0);
}

// Function to check whether a circle
// encloses the given points
function is_valid_circle(c, P)
{

    // Iterating through all the points
    // to check  whether the points
    // lie inside the circle or not
    for (var p in P)
        if (!is_inside(c, p))
            return false;
    return true;
}

// Function to return the minimum enclosing
// circle for N <= 3
function min_circle_trivial(P)
{

    if (P.length == 0) {
        return new Circle(new Point(0, 0), 0);
    }
    else if (P.length == 1) {
        return new Circle(P[0], 0);
    }
    else if (P.length == 2) {
        return circle_from(P[0], P[1]);
    }

    // To check if MEC can be determined
    // by 2 points only
    for (var i = 0; i < 3; i++) {
        for (var j = i + 1; j < 3; j++) {

            let c = circle_from(P[i], P[j]);
            if (is_valid_circle(c, P))
                return c;
        }
    }
    return circle_from(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.
function welzl_helper(P, n)
{
    // Base case when all points processed or |R| = 3
    if (n == 0 || R.length == 3) {
        return min_circle_trivial(R);
    }

    // Pick a random point randomly
    let idx = getRandomInt() % n;
    let p = P[idx];

    // Put the picked point at the end of P
    // since it's more efficient than
    // deleting from the middle of the vector
    let temp = P[idx]
    P[idx] = P[n - 1]
    P[n - 1] = temp

    // Get the MEC circle d from the
    // set of points P - {p}
    let d = welzl_helper(P, n - 1);

    // If d contains p, return d
    if (is_inside(d, p)) {
        return d;
    }

    // Otherwise, must be on the boundary of the MEC
    R.push(p);

    // Return the MEC for P - {p} and R U {p}
    return welzl_helper(P,  n - 1);
}

function welzl(P)
{
    let P_copy = [...P];

    R = []
    shuffle(P_copy);
    return welzl_helper(P_copy, P_copy.length)
}

// Driver code
let mec = welzl([new Point(0, 0), new Point(0, 1), new Point(1, 0)])

console.log("Center = {", mec.C.X, ",", mec.C.Y
         , "} Radius =", mec.R);

let mec2 = welzl([new Point(5, -2), new Point(-3, -2), new Point(-2, 5),
new Point(1, 6), new Point(0, 2)])

console.log("Center = {", mec2.C.X, ",", mec2.C.Y
         , "} Radius =", mec2.R);

// This code is contributed by phasing17

Output:

Center = { 0.5, 0.5 } Radius = 0.707107
Center = { 1, 1 } Radius = 5

Time Complexity: This algorithm has an expected time and space complexity of O(N) where N is the number of points. The space is due to the fact recursion is being used. To understand why the time complexity is linear, we can observe the number of different states to know how many calls can happen to the recursive function. 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).