Círculo envolvente mínimo | Conjunto 2 – Algoritmo de Welzl

Prerrequisitos: Ecuación de la circunferencia cuando se dan tres puntos de la circunferencia , Circunferencia mínima envolvente .
Dada una array arr[][] que contiene N puntos en un plano 2-D con coordenadas enteras. La tarea es encontrar el centro y el radio del círculo envolvente mínimo (MEC). Un círculo envolvente mínimo es un círculo en el que todos los puntos se encuentran dentro del círculo o en sus límites.
Ejemplos: 
 

Entrada: arr[][] = {{0, 0}, {0, 1}, {1, 0}} Salida 
: Centro = {0.5, 0.5}, Radio = 0.7071 
Explicación: 
Al trazar el círculo anterior con radio 0.707 y centro (0.5, 0.5), se puede observar claramente que todos los puntos mencionados se encuentran dentro o sobre el círculo. 
 

Entrada: arr[][] = {{5, -2}, {-3, -2}, {-2, 5}, {1, 6}, {0, 2}} Salida: Centro = { 
1.0 , 1.0}, Radio = 5.000 
 

Enfoque: En el artículo anterior, se discutió un enfoque ingenuo y un enfoque optimizado al llenar primero el convexo del conjunto de puntos y luego realizar un enfoque ingenuo. Aunque la solución optimizada funcionaría muy bien para ciertas entradas, la complejidad de tiempo en el peor de los casos después de esa optimización seguía siendo O(N 4 ) . En este artículo, se ha discutido un enfoque optimizado.
La idea es utilizar el algoritmo recursivo de Welzl . Usando este algoritmo, el MEC se puede encontrar en O (N). El funcionamiento del algoritmo depende de las observaciones y conclusiones extraídas del artículo anterior. La idea del algoritmo es eliminar aleatoriamente un punto del conjunto de entrada dado para formar una ecuación circular. Una vez formada la ecuación, compruebe si el punto que se eliminó está encerrado en la ecuación o no. Si no es así, entonces el punto debe estar en el límite del MEC. Por lo tanto, este punto se considera como un punto límite y la función se llama recursivamente . El funcionamiento detallado del algoritmo es el siguiente:
el algoritmo toma un conjunto de puntos P y un conjunto R que inicialmente está vacío y se usa para representar los puntos en el límite del MEC como entrada.
El caso base del algoritmo es cuando Pqueda vacío o el tamaño del conjunto R es igual a 3
 

  • Si P está vacío, entonces todos los puntos han sido procesados.
  • Si |R| = 3 , entonces ya se han encontrado 3 puntos que se encuentran en el límite del círculo, y dado que un círculo se puede determinar de manera única usando solo 3 puntos, la recursividad se puede detener.

Cuando el algoritmo alcanza el caso base anterior, devuelve la solución trivial para R, siendo: 
 

  • Si |R| = 1 , devolvemos un círculo centrado en R[0] con radio = 0
  • Si |R| = 2 , devolvemos el MEC para R[0] y R[2]
  • Si |R| = 3 , devolvemos el MEC probando los 3 pares (R[0], R[1]), (R[0], R[2]), (R[1], R[2]) 
    • Si ninguno de estos pares es válido, devolvemos el círculo definido por los 3 puntos en R

Si aún no se alcanza el caso base, hacemos lo siguiente: 
 

  • Elija un punto aleatorio p de P y elimínelo de P
  • Llame al algoritmo en P y R para obtener el círculo d
  • Si p está encerrado por d, entonces devolvemos d
  • de lo contrario, p debe estar en el límite del MEC 
    • Añadir p a R
    • Devolver la salida del algoritmo en P y R

A continuación se muestra la implementación del enfoque anterior:
 

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)
Producción: 

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

 

Complejidad de tiempo: este algoritmo tiene una complejidad de tiempo y espacio esperada de O(N) donde N es el número de puntos. El espacio se debe al hecho de que se está utilizando la recursividad . Para entender por qué la complejidad del tiempo es lineal, podemos observar la cantidad de estados diferentes para saber cuántas llamadas pueden ocurrirle a la función recursiva. Con cada llamada, el tamaño de P se reduce en uno . Además, el tamaño de R puede permanecer igual o puede aumentarse en uno. Desde |R| no puede exceder de 3 , entonces el número de estados diferentes sería 3N . Por lo tanto, esto hace que la complejidad del tiempo sea O(N) .
 

Publicación traducida automáticamente

Artículo escrito por marwan tammam y traducido por Barcelona Geeks. The original can be accessed here. Licence: CCBY-SA

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *