Saturday, 20 June 2009

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

2 comments:

  1. Funny algorithm Jordi! And good implemetation, specially with those real number comparisions that are always pulling our leg! Are you reading any mathematics book or something like that? Lots of algorithms of that kind showing up lately! Or simply maybe they were stuck in your R-side up to now? ;)

    ReplyDelete
  2. Yes man! I'm improving myself with some of the most importants Computation geometry algorithms. And this book helps me a lot: http://www.amazon.com/Computational-Geometry-Applications-Mark-Berg/dp/3540779736/ref=sr_1_1?ie=UTF8&s=books&qid=1245534486&sr=8-1
    Now, I'm very relieved because I've solved these algorithms (Convex Hull and Smallest Enclosing Circle). Now I'm going to attack some of the computer vision algorithm's.

    ReplyDelete