algorithm - applications - Dimensión lineal más grande 2d conjunto de puntos

cgal windows

Question

Dado un conjunto ordenado de ubicaciones de píxeles 2D (adyacentes o diagonales adyacentes) 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)

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))?

Answer 1

En esta página:

• http://en.wikipedia.org/wiki/Rotating_calipers
• http://cgm.cs.mcgill.ca/~orm/diam.html

muestra que puedes 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).

• http://en.wikipedia.org/wiki/Convex_hull_algorithms

Este es un código C # que encontré para calcular el casco convexo:

• http://www.itu.dk/~sestoft/gcsharp/index.html#hull

Answer 2

Mi solución original es probar un enfoque de partición binaria, donde dibujas una línea en el medio y verificas 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 intuición dice que esto es una heurística n log n que te acercará bastante.

Answer 3

Tal vez podrías dibujar un círculo que fuera más grande que el polígono y reducirlo lentamente, verificando si ya has cruzado algún punto. Entonces tu diámetro es el número que estás buscando. No estoy seguro si este es un buen método, suena entre O (n) y O (n ^ 2)

Answer 4

Transmití 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;  
    }  
}  

Answer 5

Una forma fácil es encontrar primero el casco convexo de los puntos, lo que se puede hacer en el tiempo O (n log n) de muchas maneras. [Me gusta Graham scan (ver animación ), pero el algoritmo incremental también es popular, al igual que otros , aunque algunos toman más tiempo .]

Luego puede encontrar el par más lejano (el diámetro) comenzando con dos puntos (por ejemplo, xey) en el casco convexo, moviéndose y en el sentido de las agujas del reloj hasta el más alejado de x, luego moviendo x, moviendo y nuevamente, etc. Puede 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 calibradores giratorios , como mencionaste.

Aquí está el código de David Eppstein (investigador de geometría computacional, ver también sus algoritmos de Python y estructuras de datos 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 de 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 funcione.