Implementando el Smallest Enclosing Circle con Delphi

En este artículo os muestro mi implementación para el problema del Smallest Enclosing Circle o minimum Enclosing Circle. Mi solución está basada en la utilización del algoritmo Elzinga-Hearn (1972) que es el de los más rápidos con un coste O(n^2), aunque luego en 1978 Shamos publicó un algoritmo con coste O(n logn). Más tarde en el año 1983 Nimrod Megiddo mostró un algoritmo con coste O(n) utilizando las técnicas prune-and-search Con éste, se mostró uno de los algoritmos más bonitos en geometría Computacional.
Podemos encontrar la fuente de la información en el siguiente enlace: minimax.pdf de D. Jack Elzinga, Donald W. Hearn, "Geometrical Solutions for some minimax location problems," Transportation Science, 6, (1972), pp 379 - 394. El D. Jack Elzinga publicó varios trabajos sobre diversos algoritmos y podemos encontrar más información sobre ellos en la Universidad de Florida. Por la red podemos encontrar muy buenas explicaciones de como funciona su algoritmo y applets en java demostrando su cálculo:

El algoritmo como se describe el los diferentes enlaces, consiste en conseguir el mínimo circulo 2D que engloba todos los puntos del plano mostrados. Cabe decir, que éste algoritmo se utiliza en muchos campos: estadística, Sistemas de información geográfica, militares, etc, Además su implementación es bastante sencilla. En éste artículo os dejo mi ejemplo en Delphi, utilizando diferentes objetos: Círculos, puntos y líneas puedo realizar los diferentes cálculos para conseguir el cálculo del smallest enclosing circle.

Mi aplicación llamada ThundaxEnclosingCircle permite la entrada de puntos sobre un TImage o se pueden generar puntos aleatoriamente. Luego se realiza el cálculo y se muestra el circulo que engloba todos los puntos.

• Implementando la solución

Para los tipos primitivos crearé los objetos TPointR, TLine y TCircle y como lista de puntos la TPointRList que la expliqué ayer en el post de aplicar la herencia con el TObjectList. La definición de los objetos es la siguiente:

uses
    Types, Contnrs, Lib_Math, SysUtils, Math;

type
    TPointR = class(TObject)
        X: double;
        Y: double;
        constructor Point(x, y: double); overload;
        constructor Point(point: TPointR); overload;
        function ToString(): string;
    end;

    TLine = class(TObject)
        p1, p2: TPointR;
    end;

    TCircle = class(TObject)
        Center: TPointR;
        radius: double;
        constructor Circle(); overload;
        constructor Circle(circ: TCircle); overload;
        constructor Circle(center: TPointR; radius: double); overload;
        procedure SetParam(x, y, radius : double);
    end;
  
{ TCircle }

constructor TCircle.Circle;
begin
    Self.Center := TPointR.Point(0, 0);
    Self.radius := 0.0;
end;

constructor TCircle.Circle(circ: TCircle);
begin
    Self.Center := circ.Center;
    Self.radius := circ.radius;
end;

constructor TCircle.Circle(center: TPointR; radius: double);
begin
    Self.Center := center;
    Self.radius := radius;
end;

procedure TCircle.SetParam(x, y, radius: double);
begin
    Self.center.X := x;
    Self.center.y := y;
    self.radius := radius;
end;

{ TPointR }

constructor TPointR.Point(x, y: double);
begin
    Self.x := x;
    Self.y := y;
end;

constructor TPointR.Point(point: TPointR);
begin
    Self.x := point.x;
    Self.Y := point.y;
end;

function TPointR.ToString: string;
begin
    result := FloatToStr(x) + ' - ' + FloatToStr(y);
end;

Luego, necesito una serie de funciones estándar para realizar cálculos geométricos, como la distancia entre 2 puntos, obtener el centro de 3 puntos, obtener un círculo a partir de 3 puntos, obtener el ángulo entre puntos, etc. Todas éstas están implementadas aquí:

function distancia(p1, p2: TpointR): double;
procedure swap(var p1, p2: TPointR);
function ObtenerCentro(p1, p2, p3: TPointR; var centro: TPointR): boolean;
function obtenerCirculo(p1, p2, p3: TPointR; var circulo: Tcircle): boolean;
function ObtenerAngulo(p1, p2, p3: TPointR): double;
function mismoLado(linea: TLine; p1, p2: TPointR): double;

function distancia(p1, p2: TPointR): double;
begin
    result := sqrt(sqr(p1.x - p2.x) + sqr(p1.y - p2.y));
end;

procedure swap(var p1, p2: TPointR);
var
    temp: TPointR;
begin
    temp := p1;
    p1 := p2;
    p2 := temp;
end;

function ObtenerCentro(p1, p2, p3: TPointR; var centro: TPointR): boolean;
var
    x, y: double;
    modulo1, modulo2: double;
    encontrado : boolean;
begin
    encontrado := true;
    modulo1 := 0; modulo2 := 0;

    if (Compare(p1.x,p2.x,'=') or Compare(p1.y,p2.y,'=')) then swap(p1, p3);
    if (Compare(p2.x,p3.x,'=')) then swap(p1, p2);

    if p1.x <> p2.x then modulo1 := (p2.y - p1.y) / (p2.x - p1.x)
    else encontrado := false;

    if p2.x <> p3.x then modulo2 := (p3.y - p2.y) / (p3.x - p2.x)
    else encontrado := false;

    if Compare(modulo1,0,'=') and Compare(modulo2,0,'=') then encontrado := false;
    if Compare(modulo1,modulo2,'=') then encontrado := false;

    if encontrado then
    begin
        x := (modulo1 * modulo2 * (p1.y - p3.y) + modulo2 * (p1.x + p2.x) - modulo1 * (p2.x + p3.x)) / (2 * (modulo2 - modulo1));
        if modulo1 <> 0 then y := -(x - (p1.x + p2.x) / 2) / modulo1 + (p1.y + p2.y) / 2
        else y := -(x - (p2.x + p3.x) / 2) / modulo2 + (p2.y + p3.y) / 2;
        centro.x := x;
        centro.y := y;
    end;
    result := encontrado;
end;

function obtenerCirculo(p1, p2, p3: TPointR; var circulo: Tcircle): boolean;
var
    encontrado : boolean;
begin
    encontrado := ObtenerCentro(p1, p2, p3, circulo.center);
    circulo.radius := 0;
    if encontrado then
        circulo.radius := sqrt(sqr(p1.x - circulo.center.x) + sqr(p1.y - circulo.center.y));
    result := encontrado;
end;

function ObtenerAngulo(p1, p2, p3: TPointR): double;
var
    distancia1, distancia2 : double;
begin
    distancia1 := distancia(p1,p2);
    distancia2 := distancia(p2,p3);
    result := arccos(ABS(((p1.x - p2.x) * (p2.x - p3.x) + (p1.y - p2.y) * (p2.y - p3.y)) /
        (distancia1 * distancia2)));
end;

function mismoLado(linea: TLine; p1, p2: TPointR): double;
var
    dx, dy, dx1, dy1, dx2, dy2: double;
begin
    dx := linea.p2.x - linea.p1.x;
    dy := linea.p2.y - linea.p1.y;
    dx1 := p1.x - linea.p1.x;
    dy1 := p1.y - linea.p1.y;
    dx2 := p2.x - linea.p2.x;
    dy2 := p2.y - linea.p2.y;
    result := (dx * dy1 - dy * dx1) * (dx * dy2 - dy * dx2);
end;

Por último solo necesito el algoritmo de Elzina-Hearn, que es el que tengo implementado aquí:

procedure EnclosingCircle(Sender: TObject);
    function EsAnguloAgudo(p1, p2, p3: integer; var q1, q2: integer): boolean;
    var
        angle: double;
        EsAgudo : boolean;
    begin
        EsAgudo := false;
        angle := ObtenerAngulo(p[p2], p[p1], p[p3]);
        if Compare(angle , pi/2, '<') then
        begin
            angle := ObtenerAngulo(p[p1], p[p2], p[p3]);
            if Compare(angle , pi/2, '<') then
            begin
                angle := ObtenerAngulo(p[p1], p[p3], p[p2]);
                if Compare(angle , pi/2, '<') then
                begin
                    q1 := 0; q2 := 0; EsAgudo := true;
                end
                else
                begin
                    q1 := p1; q2 := p2;
                end;
            end
            else
            begin
                q1 := p1; q2 := p3;
            end;
        end
        else
        begin
            q1 := p2; q2 := p3;
        end;
        result := EsAgudo;
    end;

    function ContieneTodosLosPuntos(circulo: TCircle): integer;
    var
        i, punto: integer;
        distanciaPunto: double;
    begin
        punto := -1;
        for i := 0 to p.count - 1 do
        begin
            distanciaPunto := distancia(p[i], circulo.center);
            if Compare(distanciaPunto, circulo.radius, '>') then
            begin
                punto := i; break;
            end;
        end;
        result := punto;
    end;

var
    Encontrado: boolean;
    distanciaMaxima, distPuntos: double;
    CentinelaI, CentinelaJ, P1, Q, i, j, PuntoPos, Seguridad, e1, e2: integer;
    circle: TCircle;
    linea: TLine;
begin
    distanciaMaxima := 0;
    CentinelaI := 0;
    CentinelaJ := 0;
    for i := 0 to p.count - 2 do
        for j := i + 1 to p.count - 1 do
        begin
            distPuntos := distancia(p[i], p[j]);
            if Compare(distPuntos, distanciaMaxima,'>') then
            begin
                CentinelaI := i;
                CentinelaJ := j;
                distanciaMaxima := distPuntos;
            end;
        end;

    circle := TCircle.Circle();
    linea := TLine.create();

    Encontrado := false;
    while not Encontrado do
    begin
        circle.center.x := (p[CentinelaI].x + p[CentinelaJ].x) / 2;
        circle.center.y := (p[CentinelaI].y + p[CentinelaJ].y) / 2;
        circle.radius := distancia(p[CentinelaI], Circle.center);
        PuntoPos := ContieneTodosLosPuntos(circle);
        if PuntoPos >= 0 then
        begin
            Seguridad := 0;
            while (not Encontrado) and (Seguridad < 1000) do
            begin
                inc(Seguridad);
                if EsAnguloAgudo(CentinelaI, CentinelaJ, PuntoPos, e1, e2) then
                begin
                    obtenerCirculo(p[CentinelaI], p[CentinelaJ], p[PuntoPos], circle);
                    P1 := ContieneTodosLosPuntos(circle);
                    if P1 >= 0 then
                    begin
                        distanciaMaxima := distancia(p[P1], p[CentinelaI]);
                        Q := CentinelaI;
                        distPuntos := distancia(p[P1], p[CentinelaJ]);
                        if distPuntos > distanciaMaxima then
                        begin
                            Q := CentinelaJ;
                            distanciaMaxima := distPuntos;
                        end;
                        if distancia(p[P1], p[PuntoPos]) > distanciaMaxima then
                            Q := PuntoPos;
                        linea.p1 := circle.center;
                        linea.p2 := TPointR(p[Q]);
                        if (Q <> CentinelaI) and (mismoLado(linea, p[P1], p[CentinelaI]) < 0) then
                        begin CentinelaJ := P1;
                            PuntoPos := q;
                        end
                        else if (Q <> CentinelaJ) and (mismoLado(linea, p[P1], p[CentinelaJ]) < 0) then
                        begin CentinelaI := P1;
                            PuntoPos := q;
                        end
                        else
                        begin CentinelaI := P1;
                            CentinelaJ := q;
                        end;
                    end
                    else
                        Encontrado := true;
                end
                else
                begin CentinelaI := e1;
                    CentinelaJ := e2;
                end;
            end;
            if Seguridad >= 1000 then
            begin
                showmessage('Bucles Infinitos (Sin Solución)');
                exit;
            end;
        end
        else
            Encontrado := true;
    end;
    Edit2.text := 'Center: x= ' + FloatToStr(circle.center.x) +
        ' y= ' + FloatToStr(circle.center.y) +
        ' Radius= ' + FloatToStr(circle.radius);
end;

Una vez pasamos el algoritmo sobre los puntos, éste nos entrega el centro del círculo mínimo y su radio. Espero que os sirva de ayuda.

• Enlaces de Interés:

Circle Covering Points
Geometrical Solution for some minimax location problem
Geometría Computacional - minimum Enclosing
Minimum Enclosing Solution in Ruby
Minimum Enclosing Solution in Matlab
Applet MinCircle Java