2008-11-26 11 views
11

Dado un conjunto ordenado de ubicaciones de píxeles 2D (adyacentes o adyacentes a la diagonal) que forman una ruta completa sin repeticiones, ¿cómo determino la dimensión lineal más grande del polígono cuyo perímetro es ese conjunto de píxeles? (donde el GLD es la mayor distancia lineal de cualquier par de puntos en el conjunto)Dimensión lineal más grande 2d conjunto de puntos

Para mis propósitos, la obvia solución O (n^2) probablemente no sea lo suficientemente rápida para cifras de miles de puntos. ¿Hay buenos métodos heurísticos o de búsqueda que acerquen la complejidad del tiempo a O (n) u O (log (n))?

Respuesta

18

Una manera fácil es encontrar primero el casco convexo de los puntos, que se puede hacer en el tiempo O (n log n) de muchas maneras. [Me gusta Graham scan (ver animation), pero el algoritmo incremental también es popular, como lo son others, aunque algunos toman more time.]

entonces usted puede encontrar el par más alejado (el diámetro) partiendo de dos puntos (por ejemplo xey) en el casco convexo, moviéndose y en el sentido de las agujas del reloj hasta que esté más alejado de x, luego moviendo x, moviendo y de nuevo, etc. Puedes probar que todo esto solo toma O (n) tiempo (amortizado). Entonces, es O (n log n) + O (n) = O (n log n) en total, y posiblemente O (nh) si usa envoltura de regalo como su algoritmo de casco convexo. Esta idea se llama rotating calipers, como mencionaste.

Aquí está code by David Eppstein (investigador de geometría computacional, consulte también su Python Algorithms and Data Structures para referencia futura).

Todo esto no es muy difícil de codificar (debe ser un centenar de líneas como máximo, es menos de 50 en el código Python anterior), pero antes de hacerlo, primero debe considerar si realmente lo necesita. Si, como dices, solo tienes "miles de puntos", entonces el algoritmo trivial O (n^2) (que compara todos los pares) se ejecutará en menos de un segundo en cualquier lenguaje de programación razonable. Incluso con un millón de puntos, no debería tomar más de una hora. :-)

Debe elegir el algoritmo más simple que funciona.

2

En esta página:

Esto demuestra que se puede determinar el diámetro máximo de un polígono convexo en O (n). Solo necesito convertir mi conjunto de puntos en un polígono convexo primero (probablemente usando el escaneo de Graham).

Aquí hay un código C# me encontré para el cálculo de la envolvente convexa:

0

Se podría quizá dibujar un círculo que era más grande que el polígono y reducirlo lentamente, comprobando si ya se ha cruzado algún punto. Entonces tu diámetro es el número que estás buscando. No estoy seguro si esto es un buen método, suena en algún lugar entre O (n) y O (n^2)

+0

¿Dónde coloca el centro? Probablemente cerca del centroide, pero apuesto a que podría llegar a situaciones en las que el centro de ese círculo tenga un impacto importante en si encuentras el GLD correcto o no. –

+2

esto es conceptualmente defectuoso - el diámetro de la circunferencia circunscrita de un triángulo equilátero es sqrt (3) veces el GLD, que es igual a la longitud lateral, – Jimmy

0

Mi solución fuera de la manga es intentar un enfoque de partición binaria, donde se dibuja una somwwhere línea en el medio y verifique las distancias de todos los puntos desde el medio de esa línea. Eso le proporcionaría 2 puntos presumiblemente muy lejanos. Luego verifique la distancia de esos dos y repita la verificación de distancia anterior. Repita este proceso por un tiempo.

Mi instinto dice que esto es una heurística n log n que te acercará bastante.

2

I porté el código de Python a C#. Parece funcionar.

using System; 
using System.Collections.Generic; 
using System.Drawing; 

// Based on code here: 
// http://code.activestate.com/recipes/117225/ 
// Jared Updike ported it to C# 3 December 2008 

public class Convexhull 
{ 
    // given a polygon formed by pts, return the subset of those points 
    // that form the convex hull of the polygon 
    // for integer Point structs, not float/PointF 
    public static Point[] ConvexHull(Point[] pts) 
    { 
     PointF[] mpts = FromPoints(pts); 
     PointF[] result = ConvexHull(mpts); 
     int n = result.Length; 
     Point[] ret = new Point[n]; 
     for (int i = 0; i < n; i++) 
      ret[i] = new Point((int)result[i].X, (int)result[i].Y); 
     return ret; 
    } 

    // given a polygon formed by pts, return the subset of those points 
    // that form the convex hull of the polygon 
    public static PointF[] ConvexHull(PointF[] pts) 
    { 
     PointF[][] l_u = ConvexHull_LU(pts); 
     PointF[] lower = l_u[0]; 
     PointF[] upper = l_u[1]; 
     // Join the lower and upper hull 
     int nl = lower.Length; 
     int nu = upper.Length; 
     PointF[] result = new PointF[nl + nu]; 
     for (int i = 0; i < nl; i++) 
      result[i] = lower[i]; 
     for (int i = 0; i < nu; i++) 
      result[i + nl] = upper[i]; 
     return result; 
    } 

    // returns the two points that form the diameter of the polygon formed by points pts 
    // takes and returns integer Point structs, not PointF 
    public static Point[] Diameter(Point[] pts) 
    { 
     PointF[] fpts = FromPoints(pts); 
     PointF[] maxPair = Diameter(fpts); 
     return new Point[] { new Point((int)maxPair[0].X, (int)maxPair[0].Y), new Point((int)maxPair[1].X, (int)maxPair[1].Y) }; 
    } 

    // returns the two points that form the diameter of the polygon formed by points pts 
    public static PointF[] Diameter(PointF[] pts) 
    { 
     IEnumerable<Pair> pairs = RotatingCalipers(pts); 
     double max2 = Double.NegativeInfinity; 
     Pair maxPair = null; 
     foreach (Pair pair in pairs) 
     { 
      PointF p = pair.a; 
      PointF q = pair.b; 
      double dx = p.X - q.X; 
      double dy = p.Y - q.Y; 
      double dist2 = dx * dx + dy * dy; 
      if (dist2 > max2) 
      { 
       maxPair = pair; 
       max2 = dist2; 
      } 
     } 

     // return Math.Sqrt(max2); 
     return new PointF[] { maxPair.a, maxPair.b }; 
    } 

    private static PointF[] FromPoints(Point[] pts) 
    { 
     int n = pts.Length; 
     PointF[] mpts = new PointF[n]; 
     for (int i = 0; i < n; i++) 
      mpts[i] = new PointF(pts[i].X, pts[i].Y); 
     return mpts; 
    } 

    private static double Orientation(PointF p, PointF q, PointF r) 
    { 
     return (q.Y - p.Y) * (r.X - p.X) - (q.X - p.X) * (r.Y - p.Y); 
    } 

    private static void Pop<T>(List<T> l) 
    { 
     int n = l.Count; 
     l.RemoveAt(n - 1); 
    } 

    private static T At<T>(List<T> l, int index) 
    { 
     int n = l.Count; 
     if (index < 0) 
      return l[n + index]; 
     return l[index]; 
    } 

    private static PointF[][] ConvexHull_LU(PointF[] arr_pts) 
    { 
     List<PointF> u = new List<PointF>(); 
     List<PointF> l = new List<PointF>(); 
     List<PointF> pts = new List<PointF>(arr_pts.Length); 
     pts.AddRange(arr_pts); 
     pts.Sort(Compare); 
     foreach (PointF p in pts) 
     { 
      while (u.Count > 1 && Orientation(At(u, -2), At(u, -1), p) <= 0) Pop(u); 
      while (l.Count > 1 && Orientation(At(l, -2), At(l, -1), p) >= 0) Pop(l); 
      u.Add(p); 
      l.Add(p); 
     } 
     return new PointF[][] { l.ToArray(), u.ToArray() }; 
    } 

    private class Pair 
    { 
     public PointF a, b; 
     public Pair(PointF a, PointF b) 
     { 
      this.a = a; 
      this.b = b; 
     } 
    } 

    private static IEnumerable<Pair> RotatingCalipers(PointF[] pts) 
    { 
     PointF[][] l_u = ConvexHull_LU(pts); 
     PointF[] lower = l_u[0]; 
     PointF[] upper = l_u[1]; 
     int i = 0; 
     int j = lower.Length - 1; 
     while (i < upper.Length - 1 || j > 0) 
     { 
      yield return new Pair(upper[i], lower[j]); 
      if (i == upper.Length - 1) j--; 
      else if (j == 0) i += 1; 
      else if ((upper[i + 1].Y - upper[i].Y) * (lower[j].X - lower[j - 1].X) > 
       (lower[j].Y - lower[j - 1].Y) * (upper[i + 1].X - upper[i].X)) 
       i++; 
      else 
       j--; 
     } 
    } 

    private static int Compare(PointF a, PointF b) 
    { 
     if (a.X < b.X) 
     { 
      return -1; 
     } 
     else if (a.X == b.X) 
     { 
      if (a.Y < b.Y) 
       return -1; 
      else if (a.Y == b.Y) 
       return 0; 
     } 
     return 1; 
    } 
} 
Cuestiones relacionadas