viernes, 29 de octubre de 2010

El problema de los puntos más cercanos

Ahora que estoy con el proyecto todo el día, estoy redescubriendo problemas de cuando estaba en la carrera. El caso es que en uno de los pasos de la implementación, necesitaba obtener los dos puntos más cercanos.

El problema de los puntos más cercanos, es un problema geométrico, en el que dado n puntos (en mi caso en 2D, es decir en un plano), quiero encontrar un par de puntos cuya distancia entre ellos sea la menor.

Inicialmente pensé, uy que problema más fácil, con dos bucles anidados resuelvo el problema:


int p1, p2;
double dist, dist_min = 100000000;

for (int i = 0; i < num_puntos; i++){

    for (int i = 0; i < num_puntos; i++){

        if (i != j){
            dist = sqrt((puntos[i].X - puntos[j].X)^2 + (puntos[i].Y - puntos[j].Y)^2);

            if (dist < dist_min){

                p1 = i;
                p2 = j;
                dist = dist_min;
            }
        }
    }
}



Este caso serie por la fuerza bruta, sin pensar en ningún tipo de eficiencia ni nada por el estilo, Por lo que si tienes que repetir esta búsqueda muchas veces, el programa puede tirarse un ratillo en acabar ... Y ya que estamos os hablo de eficiencia. En el caso de la fuerza bruta, la eficiencia sería de O(n2), donde n es el número de puntos.

Hay varias formas de mejorar esta eficiencia, y resolver el problema en un tiempo menor:

  • Usar un algorimo divide y vencerás.Donde se consiguie una eficiencia de O(n log2n)
  • Usar programación dinámica y crear una serie de estructuras de datos para facilitar el problema. En este caso se podría llegar a alcanzar una eficiencia de O(log n)
En mi caso he implementado la primera opción. El algoritmo divide y vencerás realizará lo siguiente:
  1. Ordenamos los puntos según su componente X
  2. El caso base será cuando queden 2 ó 3 puntos, en ese caso calcular la distancia mínima.
  3. Si el número de puntos es mayor, calcular la distancia mínima como el mínimo entra la distancia en el lado derecho e izquierdo.
  4. Además, tenemos que tener en cuenta los puntos cerca de la zona de división, en este caso se calcularán los puntos que estén a menos de dist_min de la línea L. Esta nueva lista de puntos se ordena con respecto a la componente Y, y obtenemos la distancia mínima. Si esta es menor que la distancia mínimal actual se modifica el valor.
No sé si está muy claro el ejemplo pero de todas formas os dejo el trozo de código (en c++) que he implementado.


punto calcular_distancia_minima(int ini, int fin){
    punto p1, p2, p, paux;
    double dist1, dist2, dist3;

    if (fin - ini + 1 > 3){
        p1 = calcular_distancia_minima(ini, ini + (fin - ini + 1)/2 - 1);
        p2 = calcular_distancia_minima(ini + (fin - ini + 1)/2, fin);

        if (p1.get_peso() < p2.get_peso())

            p = p1;
        else
            p = p2;

        vector c1;

        for (int i = ini; i <= fin; i++){

            if (puntos[i].X > (puntos[ini + (fin - ini + 1)/2].X - p.distancia)){
                paux.X = puntos[i].X;
                paux.Y = puntos[i].Y;
                paux.posicion = i;
                c1.aniadir(paux);
            }

            else {

                if (puntos[i].X < (puntos[ini + (fin - ini + 1)/2].X + p.distancia)){

                    paux.X = puntos[i].X;
                    paux.Y = puntos[i].Y;
                    paux.posicion = i;
                    c1.aniadir(paux);
                }
            }
        }

        c1.mergesortY(0, c1.tamanio()-1);

        double dist;
        int j;

        for (int i = 0; i < c1.tamanio() - 1; i++){

            j = i + 1;

            while ((j < c1.tamanio()) && (c1[j].Y - c1[i].Y < p.distancia)){                 dist = sqrt((c1[i].X - c1[j].X)^2 + (c1[i].Y - c1[j].Y)^2);

                if (dist < p.distancia){

                    p.X = c1[i].posicion;
                    p.Y = c1[j].posicion;
                    p.distancia = dist;
                }

                j++;
            }
        }
    }

    else {

        if (fin - ini + 1 == 2){
            p.X = ini;
            p.Y = fin;
            p.distancia = sqrt((puntos[ini].X - puntos[fin].X)^2 + (puntos[ini].Y - puntos[fin].Y)^2);
        }

        else{
            dist1 = sqrt((puntos[ini].X - puntos[fin].X)^2 + (puntos[ini].Y - puntos[fin].Y)^2);
            dist2 = sqrt((puntos[ini].X - puntos[ini + 1].X)^2 + (puntos[ini].Y - puntos[ini + 1].Y)^2);
            dist3 = sqrt((puntos[ini + 1].X - puntos[fin].X)^2 + (puntos[ini + 1].Y - puntos[fin].Y)^2);

            if (dist1 <= dist2 && dist1 <= dist3){

                p.X = ini;
                p.Y = fin;
                p.distancia = dist1;
            }

            else{

                if (dist2 <= dist1 && dist2 <= dist3){

                    p.X = ini;
                    p.Y = ini + 1;
                    p.distancia = dist2;
                }

                else{
                    p.X = ini + 1;
                    p.Y = fin;
                    p.distancia = dist3;
                }
            }
        }
    }

    return p;
}


Programa main que llama a la función de calcular distancia:


// creamos un vector de puntos
vector puntos(2000);
puntos.inicializar();

int main(){
    puntos.mergesortX(0, puntos.tamanio() - 1);
    punto p_distancia = calcular_distancia_minima(0, puntos.tamanio() - 1);
}


Fuentes:
Closest Pair of Points Problem - Wikipedia
Closest Pair Problem - Subhash Suri
losest Pair of Points - Rahul Gokhale

2 comentarios:

  1. Bastante interesante el ejemplo , yo lo estaba implementando para una practica de la facultad y me ha venido de perlas como guia

    Una duda , que hace la funcion get_peso() de la clase punto ?

    ResponderEliminar
  2. Muy buenas, el ejemplo que he puesto es sacado de algo que he estado implementando. Y se me ha colado, en realidad get_peso() sería igual que distancia; es decir, que p1.get_peso() = p1.distancia.

    Lo que me pasa siempre cuando paso de código a pseudocódigo :P.

    Y me alegro que te sirviera Francisco Javier :).

    ResponderEliminar