Dichtestes Punktpaar

Das Problem des dichtesten Punktpaares (englisch closest pair of points problem) ist die Suche nach den zwei am dichtesten beieinander liegenden Punkten in einer Ebene. Gegeben ist eine beliebige Menge von Punkten in der Ebene und gesucht sind zwei dieser Punkte, sodass der euklidische Abstand minimal ist. Ein ähnliches Problem ist die Suche nach den zwei am weitesten voneinander entfernten Punkten in der Ebene, also den zwei Punkten mit dem maximalen euklidischen Abstand.

Kleinster Abstand von Punkten in der Ebene

Brute-force-Algorithmus

Der Brute-force-Algorithmus berechnet die Abstände zwischen allen möglichen Punktpaaren und wählt das Punktpaar mit dem kleinsten Abstand aus. Wenn $n$ die Anzahl der Punkte ist, gibt es $(\genfrac{}{}{0ex}{}{n}{2})=\frac{n\cdot (n-1)}{2}=\frac{1}{2}\cdot (n^2-n)$ mögliche Punktpaare, für die die Abstände berechnet werden müssen. Die Laufzeit des Algorithmus ist also quadratisch und liegt in $O(n^2)$.

Divide-and-conquer-Algorithmus

Vorlage:Siehe auch Zunächst wird die gegebene Menge der Punkte einmal nach x-Koordinaten und einmal nach y-Koordinaten sortiert. Man erhält zwei sortierte Listen $L_x$ und $L_y$.

Der Divide-Schritt teilt die nach x-Koordinaten sortierte Liste $L_x$ in zwei Hälften $P_L$ und $P_R$ links und rechts des Medians auf.

Der rekursive Aufruf geschieht nun jeweils für die beiden Hälften. Man erhält das jeweils dichteste Punktpaar beider Hälften, ohne dass eventuelle Punktpaare zwischen beiden Hälften berücksichtigt werden.

Der Conquer-Schritt kombiniert nun diese beiden Hälften. Es wird das Punktpaar mit dem kleinsten gefundenen Abstand ausgewählt. Hierbei sind zwei Fälle zu beachten:

1. Das Punktpaar liegt in der linken oder der rechten Hälfte. In diesem Fall gibt es keine weiteren Probleme.
2. Es gibt zwei Punkte in unterschiedlichen Hälften, deren Abstand kleiner ist, als der bisher in den Hälften gefundene. Dieser Fall ist nun gesondert zu berücksichtigen.
   

   Datei:Punkteinfuegen d and c.jpg

   Prüfung der Nachbarn von P

Es ist nicht nötig, alle solchen Punktpaare einzeln durchzuprüfen. Ist $\delta$ der kleinste gefundene Abstand in einer der beiden Hälften, dann ist dies eine obere Grenze für den kleinsten Abstand. Daher müssen nur noch Punkte betrachtet werden, deren Abstand zum Median in x-Richtung höchstens $\delta$ beträgt. Außerdem müssen für diese Punkte nur noch Partner betrachtet werden, deren Abstand in y-Richtung kleiner als $\delta$ ist.
Daraus ergibt sich für jeden dieser Punkte, dass lediglich der Abstand zu einer konstanten Anzahl von anderen Punkten zu prüfen ist: Denkt man sich ein Gitter der Weite $\frac{\delta }{2}$ in der Umgebung des Medians, dann kann in jeder Gitterzelle nur einer dieser Punkte liegen, weil sonst bereits ein Abstand kleiner als $\delta$ gefunden worden wäre. Weil nur diejenigen Gitterzellen zu prüfen sind, die von P aus höchstens Abstand $\delta$ haben, sind maximal 24 weitere Abstände für jeden Punkt zu berechnen, nämlich jeweils maximal 12 Abstände nach oben und nach unten, womit statt quadratischer Laufzeit nur noch lineare Laufzeit notwendig ist. Insgesamt ergibt sich für die Anzahl der berechneten Abstände die Rekursionsgleichung $T(n)=2\cdot T(\frac{n}{2})+24\cdot n$, sodass die Laufzeit in $O(n\cdot \log n)$ liegt.

Programmierung

Das folgende Computerprogramm in der Programmiersprache C++ zeigt eine Implementierung des Divide-and-conquer-Algorithmus. Bei der Ausführung des Programms wird die Methode main verwendet, die das Punktpaar und den kleinsten Abstand zwischen dem Punktpaar das auf der Konsole ausgibt.

Die Funktion getSmallestDistanceBruteForce verwendet für Teilbereiche aus 2 oder 3 Punkten den Brute-force-Algorithmus. Der Conquer-Schritt wird von der rekursiven Funktion getSmallestDistanceRecursive durchgeführt, die sich selbst für die linke Hälfte und für die rechte Hälfte aufruft. Die Funktion getSmallestDistanceOfStrip bestimmt das Punktpaar und den kleinsten Abstand für den Streifen, der Punkte aus beiden Hälften enthält. Die Laufzeit dieser Funktion ist nicht linear zur Anzahl der Punkte, da sortiert wird. Dadurch ist die Laufzeit dieser Implementierung

$𝒪(n\cdot {\log }^2(n))$

.

#include <iostream>
#include <sstream>
using namespace std;

// Deklariert den Datentyp für Punkte mit x-Koordinate und y-Koordinate in der Ebene
struct Point
{
    int x, y;
};

// Vergleichsfunktion für das Sortieren der Punkte nach der x-Koordinate
int compareX(const void* a, const void* b)
{
    Point* point1 = (Point*)a, * point2 = (Point*)b;
    return point1->x - point2->x;
}
// Vergleichsfunktion für das Sortieren der Punkte nach der y-Koordinate
int compareY(const void* a, const void* b)
{
    Point* point1 = (Point*)a, * point2 = (Point*)b;
    return (point1->y - point2->y);
}

// Diese Funktion berechnet den euklidischen Abstand zwischen zwei Punkten
double getDistance(Point point1, Point point2)
{
    return sqrt((point1.x - point2.x) * (point1.x - point2.x) + (point1.y - point2.y) * (point1.y - point2.y));
}

// Diese Funktion gibt den kleinsten Abstand der Punkte bis zum gegebenen Index zurück, indem alle Paare von Punkten durchlaufen und geprüft werden
double getSmallestDistanceBruteForce(Point points[], Point& point1, Point& point2, int index)
{
    double minimumDistance = FLT_MAX;
    for (int i = 0; i < index; i++)
    {
        for (int j = i + 1; j < index; j++)
        {
            // Die for-Schleifen durchlaufen alle Paare von Punkten
            double currentDistance = getDistance(points[i], points[j]); // Ruft die Funktion für die Berechnung des euklidischen Abstands auf
            if (currentDistance < minimumDistance) // Prüft, ob der Abstand zwischen den aktuellen Punkten kleiner als der bisher kleinste Abstand ist
            {
                point1 = points[i];
                point2 = points[j];
                minimumDistance = currentDistance;
            }
        }
    }
    return minimumDistance;
}

// Diese Funktion gibt den kleinsten Abstand der Punkte im gegebenen Streifen zurück, wenn der Abstand kleiner als der vorgegebene Abstand ist
double getSmallestDistanceOfStrip(Point points[], Point& point1, Point& point2, int index, double distance)
{
    double minimumDistance = distance; // Initialisiert den kleinsten Abstand
    qsort(points, index, sizeof(Point), compareY); // Sortiert die Punkte nach der y-Koordinate mithilfe der Vergleichsfunktion compareY
    for (int i = 0; i < index; i++)
    {
        for (int j = i + 1; j < index && points[j].y - points[i].y < minimumDistance; j++) // Durchläuft die Punkte, so lange die Differenz der y-Koordinaten kleiner als der vorgegebene Abstand ist
        {
            double currentDistance = getDistance(points[i], points[j]); // Ruft die Funktion für die Berechnung des euklidischen Abstands auf
            if (currentDistance < minimumDistance) // Prüft, ob der Abstand zwischen den aktuellen Punkten kleiner als der bisher kleinste Abstand ist
            {
                point1 = points[i];
                point2 = points[j];
                minimumDistance = currentDistance;
            }
        }
    }
    return minimumDistance;
}

// Diese rekursive Funktion gibt den kleinsten Abstand der Punkte zurück. Sie verwendet das Teile-und-herrsche-Verfahren, indem sie die Menge der Punkte in zwei Hälften teilt und den kleinsten Abstand für die beiden Hälften und den Streifen zwischen den beiden Hälften berechnet.
double getSmallestDistanceRecursive(Point points[], Point& point1, Point& point2, int index)
{
    if (index <= 3)
    {
        return getSmallestDistanceBruteForce(points, point1, point2, index); // Wenn das Array aus 2 oder 3 Punkten besteht, wird die Brute Force Funktion aufgerufen
    }
    int medianIndex = index / 2;
    Point medianPoint = points[medianIndex];
    Point leftPoint1, leftPoint2, rightPoint1, rightPoint2; // Deklariert Referenzen für die Punkte mit dem kleinsten Abstand
    double leftDistance = getSmallestDistanceRecursive(points, leftPoint1, leftPoint2, medianIndex); // Rekursiver Aufruf der Funktion für die linke Hälfte der Punkte
    double rightDistance = getSmallestDistanceRecursive(points + medianIndex, rightPoint1, rightPoint2, index - medianIndex); // Rekursiver Aufruf der Funktion für die rechte Hälfte der Punkte
    // Bestimmt den kleineren der beiden Abstände und weist die Referenzen auf die Punkte zu
    double distance;
    if (leftDistance < rightDistance) // Wenn der Abstand in der linken Hälfte kleiner ist
    {
        point1 = leftPoint1;
        point2 = leftPoint2;
        distance = leftDistance;
    }
    else
    {
        point1 = rightPoint1;
        point2 = rightPoint2;
        distance = rightDistance;
    }
    // Erzeugt ein Array mit einem Streifen von Punkten auf, deren x-Koordinate einen kleineren Abstand hat
    Point* stripPoints = new Point[index];
    int j = 0;
    for (int i = 0; i < index; i++)
    {
        if (abs(points[i].x - medianPoint.x) < distance) // Prüft, ob der Abstand der x-Koordinate zwischen den aktuellen Punkten kleiner als der bisher kleinste Abstand ist
        {
            stripPoints[j] = points[i];
            j++;
        }
    }
    getSmallestDistanceOfStrip(stripPoints, point1, point2, j, distance); // Aufruf der Funktion für den Streifen von Punkten aus beiden Hälften. Wenn der Abstand in diesem Streifen kleiner ist, wird dieser zurückgegeben.
}

// Diese Funktion gibt den kleinsten Abstand der Punkte und Referenzen auf die zwei Punkte mit dem kleinsten Abstand zurück
double getSmallestDistance(Point points[], Point& point1, Point& point2, int numberOfPoints)
{
    qsort(points, numberOfPoints, sizeof(Point), compareX); // Sortiert die Punkte nach der x-Koordinate mithilfe der Vergleichsfunktion compareX
    return getSmallestDistanceRecursive(points, point1, point2, numberOfPoints); // Aufruf der rekursiven Funktion
}

// Diese Funktion gibt den Punkt in der Form (x, y) zurück
string pointToString(Point* point)
{
    stringstream text;
    text << "(" << point->x << ", " << point->y << ")" << endl;
    return text.str(); // Typumwandlung von stringstream nach string
}

// Hauptfunktion die das Programm ausführt
int main()
{
    Point points[] = { {2, 3}, {12, 30}, {40, 50}, {5, 1}, {12, 10}, {3, 4} }; // Deklariert und initialisiert ein Array mit 6 Punkten
    int numberOfPoints = sizeof(points) / sizeof(points[0]);
    Point point1, point2; // Deklariert Referenzen für die zwei Punkte mit dem kleinsten Abstand
    cout << "Der kleinste Abstand zwischen zwei Punkten ist " << getSmallestDistance(points, point1, point2, numberOfPoints) << endl; // Aufruf der Funktion, Ausgabe auf der Konsole
    cout << "Erster Punkt: " << pointToString(&point1); // Ausgabe auf der Konsole
    cout << "Zweiter Punkt: " << pointToString(&point2); // Ausgabe auf der Konsole
}

Folgende Implementierung hat die Laufzeit

$\mathrm{\Theta }(n\cdot \log (n))$

des Ansatzes Divide and Conquer.

Code-Schnipsel

#include <iostream> #include <vector> #include <algorithm> #include <cmath> #include <cfloat> using namespace std; struct Point { int x, y; }; // Überschreibt den << Operator, um einen Punkt im Format (x, y) auszugeben ostream& operator<<(ostream &os, const Point &p) { os << "(" << p.x << ", " << p.y << ")"; return os; } // Berechnet die euklidische Distanz zwischen zwei Punkten double getDistance(const Point &p1, const Point &p2) { return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)); } // Eine Brute-Force Methode für kleine Eingabegrößen double bruteForce(const vector<Point>& points, Point &best1, Point &best2) { double minDist = FLT_MAX; int n = points.size(); for (int i = 0; i < n; i++) { for (int j = i + 1; j < n; j++) { double d = getDistance(points[i], points[j]); if (d < minDist) { minDist = d; best1 = points[i]; best2 = points[j]; } } } return minDist; } // Gegeben ist ein Streifen von Punkten (nach y sortiert). // Finde das nächstgelegene Punktpaar im Streifen. double stripClosest(const vector<Point>& strip, double d, Point &best1, Point &best2) { double minDist = d; int n = strip.size(); for (int i = 0; i < n; i++) { // Da der Streifen nach y sortiert ist, werden nur Punkte verglichen, deren y-Differenz kleiner als minDist ist. // Hier sollte nochmals hervorgehoben werden, dass die innere Schleife abbricht, sobald die Distanz überschritten wird. // Genau dies macht die Anzahl an Iterationen in dieser Schleife konstant. for (int j = i + 1; j < n && (strip[j].y - strip[i].y) < minDist; j++) { double dCurrent = getDistance(strip[i], strip[j]); if (dCurrent < minDist) { minDist = dCurrent; best1 = strip[i]; best2 = strip[j]; } } } return minDist; } // Rekursive Funktion, die die kleinste Distanz berechnet. // Die Funktion erhält sowohl die Punkte sortiert nach x (Px) als auch nach y (Py). double closestPairRec(vector<Point>& Px, vector<Point>& Py, Point &best1, Point &best2) { int n = Px.size(); if (n <= 3) { return bruteForce(Px, best1, best2); } int mid = n / 2; Point midPoint = Px[mid]; // Teile die Punkte in Px in linke und rechte Hälften. vector<Point> Qx(Px.begin(), Px.begin() + mid); vector<Point> Rx(Px.begin() + mid, Px.end()); // Teile Py in Qy und Ry, sodass Qy die Punkte aus Qx enthält und Ry die Punkte aus Rx. // Dadurch bleibt die Sortierung nach y erhalten. vector<Point> Qy, Ry; for (const auto &p : Py) { if (p.x <= midPoint.x) { Qy.push_back(p); } else { Ry.push_back(p); } } Point best1Left, best2Left, best1Right, best2Right; double dLeft = closestPairRec(Qx, Qy, best1Left, best2Left); double dRight = closestPairRec(Rx, Ry, best1Right, best2Right); double d; if (dLeft < dRight) { d = dLeft; best1 = best1Left; best2 = best2Left; } else { d = dRight; best1 = best1Right; best2 = best2Right; } // Erzeuge den Streifen: Punkte in Py, deren x-Koordinate innerhalb von d um midPoint.x liegt. vector<Point> strip; for (const auto &p : Py) { if (abs(p.x - midPoint.x) < d) { strip.push_back(p); } } // Überprüfe den Streifen auf ein näherliegendes Paar. double dStrip = stripClosest(strip, d, best1, best2); return min(d, dStrip); } // Die Hauptfunktion, die die sortierten Arrays vorbereitet und die rekursive Funktion aufruft. double closestPair(vector<Point>& points, Point &best1, Point &best2) { vector<Point> Px = points; // Punkte sortiert nach x vector<Point> Py = points; // Punkte sortiert nach y sort(Px.begin(), Px.end(), [](const Point &a, const Point &b) { return a.x < b.x; }); sort(Py.begin(), Py.end(), [](const Point &a, const Point &b) { return a.y < b.y; }); return closestPairRec(Px, Py, best1, best2); } int main() { vector<Point> points = { {2, 3}, {12, 30}, {40, 50}, {5, 1}, {12, 10}, {3, 4} }; Point best1, best2; double dist = closestPair(points, best1, best2); cout << "Kleinste Distanz: " << dist << "\n"; cout << "Nächstgelegenes Paar: " << best1 << " und " << best2 << "\n"; return 0; }

Folgende Implementierung hat $𝒪(n)$ Speicherverbrauch.

Code-Schnipsel

#include <iostream> #include <vector> #include <algorithm> #include <cmath> #include <cfloat> #include <iterator> using namespace std; struct Point { int x, y; }; // Überschreibt den << Operator, um einen Punkt im Format (x, y) auszugeben ostream& operator<<(ostream &os, const Point &p) { os << "(" << p.x << ", " << p.y << ")"; return os; } // Berechnet die euklidische Distanz zwischen zwei Punkten double getDistance(const Point &p1, const Point &p2) { return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)); } // Brute-Force-Lösung für sehr kleine Bereiche (weniger als oder gleich 3 Punkte) // Diese Version arbeitet mit Iteratoren für den Bereich [xBegin, xEnd) template<typename It> double bruteForce(It xBegin, It xEnd, Point &best1, Point &best2) { double minDist = FLT_MAX; for (auto it1 = xBegin; it1 != xEnd; ++it1) { for (auto it2 = next(it1); it2 != xEnd; ++it2) { double d = getDistance(*it1, *it2); if (d < minDist) { minDist = d; best1 = *it1; best2 = *it2; } } } return minDist; } // Überprüft den Streifen (bereits nach y sortiert) auf ein näherliegendes Paar. // Arbeitet mit Iteratoren, um keine zusätzlichen Kopien anzulegen. template<typename It> double stripClosest(It begin, It end, double d, Point &best1, Point &best2) { double minDist = d; for (auto it1 = begin; it1 != end; ++it1) { auto it2 = next(it1); for (; it2 != end && ((it2->y - it1->y) < minDist); ++it2) { double dCurrent = getDistance(*it1, *it2); if (dCurrent < minDist) { minDist = dCurrent; best1 = *it1; best2 = *it2; } } } return minDist; } // Rekursive Funktion, die das nächstgelegene Punktpaar findet. // Es werden Iteratoren für den x‑sortierten Bereich [xBegin, xEnd) // und den y‑sortierten Bereich [PyBegin, PyEnd) übergeben. // 'aux' ist ein Hilfsvektor, der einmalig O(n) Speicher reserviert und wiederverwendet wird. double closestPairRec(vector<Point>::iterator xBegin, vector<Point>::iterator xEnd, vector<Point>::iterator PyBegin, vector<Point>::iterator PyEnd, vector<Point>& aux, Point &best1, Point &best2) { int n = distance(xBegin, xEnd); if (n <= 3) { return bruteForce(xBegin, xEnd, best1, best2); } // Bestimme den mittleren Iterator im x‑sortierten Bereich. auto xMid = xBegin + n / 2; Point midPoint = *xMid; // In-Place Partitionierung von Py: // Teile den Bereich [PyBegin, PyEnd) so, dass alle Punkte mit p.x <= midPoint.x am Anfang stehen. auto midIt = stable_partition(PyBegin, PyEnd, [&](const Point &p) { return p.x <= midPoint.x; }); Point best1Left, best2Left, best1Right, best2Right; double dLeft = closestPairRec(xBegin, xMid, PyBegin, midIt, aux, best1Left, best2Left); double dRight = closestPairRec(xMid, xEnd, midIt, PyEnd, aux, best1Right, best2Right); double d; if (dLeft < dRight) { d = dLeft; best1 = best1Left; best2 = best2Left; } else { d = dRight; best1 = best1Right; best2 = best2Right; } // Erzeuge den Streifen: Durchlaufe den gesamten Py-Bereich und wähle Punkte, // deren Abstand zur vertikalen Linie (midPoint.x) kleiner als d ist. aux.clear(); for (auto it = PyBegin; it != PyEnd; ++it) { if (abs(it->x - midPoint.x) < d) { aux.push_back(*it); } } Point best1Strip, best2Strip; double dStrip = stripClosest(aux.begin(), aux.end(), d, best1Strip, best2Strip); if (dStrip < d) { d = dStrip; best1 = best1Strip; best2 = best2Strip; } return d; } // Diese Funktion bereitet die sortierten Vektoren vor und ruft den rekursiven Algorithmus auf. double closestPair(vector<Point>& points, Point &best1, Point &best2) { int n = points.size(); vector<Point> Px = points; // Punkte sortiert nach x vector<Point> Py = points; // Punkte sortiert nach y sort(Px.begin(), Px.end(), [](const Point &a, const Point &b) { return a.x < b.x; }); sort(Py.begin(), Py.end(), [](const Point &a, const Point &b) { return a.y < b.y; }); vector<Point> aux; aux.reserve(n); // Reserviere O(n) Speicher für den Hilfsvektor return closestPairRec(Px.begin(), Px.end(), Py.begin(), Py.end(), aux, best1, best2); } int main() { vector<Point> points = { {2, 3}, {12, 30}, {40, 50}, {5, 1}, {12, 10}, {3, 4} }; Point best1, best2; double dist = closestPair(points, best1, best2); cout << "Kleinste Distanz: " << dist << "\n"; cout << "Nächstgelegenes Paar: " << best1 << " und " << best2 << "\n"; return 0; }

Randomisierter Algorithmus

Funktionsweise

Der randomisierte Algorithmus beruht darauf, dass die Punkte in zufälliger Reihenfolge in eine Hashtabelle eingefügt werden, wobei das Einfügen eines Punktes ebenso wie das Testen, ob er den vorher bekannten minimalen Abstand $\delta$ unterbietet, konstante Laufzeit hat. Die Hashtabelle bildet alle Gitterzellen der Größe ${\left(\frac{\delta }{2}\right)}^2$ auf den eventuell darin liegenden Punkt ab. Es muss zwar bei jeder solchen Aktualisierung von $\delta$ die Hashtabelle neu aufgebaut werden, die Gesamtzahl der Einfügeoperationen in sämtliche dieser Hashtabellen liegt aber trotzdem in $O(n)$.
Ohne Beachtung, wie die Hashtabelle konkret genutzt wird, ergibt sich folgender Algorithmus:

Bilde eine zufällige Reihenfolge der Punkte P1, ..., Pn
δ = Abstand zwischen P1 und P2
Füge P1 und P2 in die Hashtabelle ein (Gittergröße δ/2)
Für i = 3, ..., n
    Falls Pi einen Abstand δ < δ' zu einem der Punkte P1, ..., Pi-1 hat:
        δ = δ'
        Baue die Hashtabelle neu auf (Gittergröße δ'/2)
    Füge Pi in die Hashtabelle ein

Gitterstruktur und Hashfunktion

Man kann vereinfachend annehmen, dass alle Punkte Koordinaten zwischen (0,0) und (1,1) haben. Gegebenenfalls kann hierfür eine Koordinatentransformation vorgenommen werden. Die Ebene, in der die Punkte liegen, wird sodann durch ein Gitter der Weite $\frac{\delta }{2}$ überdeckt. Es wird nun eine universelle Hashfunktion benötigt; sie ermöglicht es einerseits, von einer gegebenen Gitterzelle in konstanter Laufzeit festzustellen, ob sich in dieser Gitterzelle einer der Punkte befindet, und andererseits neue Punkte in konstanter Laufzeit in die bestehende Hashtabelle einzufügen.

Einfügen neuer Punkte

Datei:Engstes punktpaar pkt einfuegen.jpg

Einfügen von P

Bei jedem Einfügen eines neuen Punktes P ist zu prüfen, ob dadurch der bereits bekannte minimale Abstand $\delta$ unterschritten wird. Anhand der Koordinaten von P lässt sich sofort bestimmen, in welcher der Gitterzellen der Punkt P liegt. Aufgrund der Aufteilung in Gitterzellen der Größe ${\left(\frac{\delta }{2}\right)}^2$ muss nur festgestellt werden, ob in einer der umliegenden Gitterzellen schon ein Punkt liegt. Umliegend sind dabei alle Gitterzellen, die einen solchen Abstand begründen könnten, also nur das umgebende 5x5-Rechteck. Alle Punkte, die außerhalb dieses Bereichs liegen, können keinen kleineren Abstand zu P haben, weil sie aufgrund der Gittergröße mindestens den Abstand $2\cdot \frac{\delta }{2}=\delta$ haben. Weil in jeder Gitterzelle nur ein einziger Punkt liegen kann, denn sonst wäre vorher bereits ein kleinerer Abstand gefunden worden, müssen also auch nur höchstens 25 Punkte geprüft werden. Sofern in diesem Bereich keine Punkte liegen, kann P bedenkenlos in die bestehende Hashtabelle eingefügt werden.

Sind jedoch in diesem Bereich bereits weitere Punkte vorhanden, so wird der Abstand $\delta$ auf den neu gefundenen minimalen Abstand gesetzt. Dies hat zur Folge, dass die bisherige Hashtabelle nutzlos geworden ist, weil ihre Gitterweite nicht mehr mit $\frac{\delta }{2}$ übereinstimmt. Wir benötigen also eine neue Hashtabelle mit korrektem $\delta$. Wir erstellen einfach eine solche Tabelle für alle Punkte bis zu P von Grund auf neu. Der Effizienz halber kann man sich beim Neuerstellen die Abstandsprüfung sparen, weil der minimale Abstand aller Punkte zu P bereits bekannt ist.

Schließlich hat man nach dem Einfügen eines neuen Punktes immer eine korrekte Hashtabelle mit Gitterweite $\frac{\delta }{2}$ und in jedem Gitterquadrat liegt höchstens ein Punkt.

Laufzeit

Relevant für die Laufzeit ist lediglich die Anzahl der Einfüge-Operationen neuer Punkte. Trivialerweise muss jeder Punkt mindestens einmal eingefügt werden, also ist die Laufzeit mindestens linear.

Fraglich ist jedoch, wie sich der gelegentlich vorkommende Neuaufbau der Hashtabelle auswirkt, weil hierfür alle bereits bekannten Punkte erneut eingefügt werden müssen. Hierfür ist es notwendig, die Wahrscheinlichkeit anzugeben, mit der beim Einfügen eines Punkts ein Neuaufbau erforderlich wird, und die dafür notwendigen Kosten einzuberechnen. Intuitiv sieht man, dass mit zunehmender Anzahl der Punkte der Aufwand für die Reorganisation immer größer wird, die Wahrscheinlichkeit dafür aber immer kleiner.

Der Erwartungswert für die Laufzeit kann wie folgt berechnet werden:

Die Wahrscheinlichkeit dafür, dass der Punkt $P_i$ beim Einfügen einen kleineren Abstand erzeugt, ist gleich $\frac{2}{i}$, denn dafür müsste $P_i$ einer der 2 Punkte sein, die diesen Abstand zueinander haben.

Definieren wir nun die Zufallsvariable $X_i⟶\{0,1\}$. Sie sei 1, falls beim Einfügen des Punkts $P_i$ ein kleinerer Abstand entsteht, und sonst 0.

Dann gilt: Die Gesamtanzahl der Einfügeoperationen ist $n+\sum _{i=2}^{n}i\cdot X_i$, also die Anzahl $n$ der eingefügten Punkte plus die Anzahl der Einfügeoperationen durch die Neuorganisationen der Hashtabelle.

Für den Erwartungswert gilt aufgrund dessen Linearität: $E(X)=n+\sum _{i=2}^{n}i\cdot E(X_i)=n+\sum _{i=2}^{n}i\cdot \frac{2}{i}=n+2\cdot (n-1)=3\cdot n-2$.

Das heißt, die Gesamtanzahl der erwarteten Einfügeoperationen ist lediglich gleich $3\cdot n-2$. Insgesamt ist die erwartete Laufzeit somit linear, liegt also in $O(n)$.

Größter Abstand von Punkten in der Ebene

Rotating calipers Algorithmus

Der Rotating calipers Algorithmus ist danach benannt, dass ein Messschieber um die Außenseite eines konvexen Polygons gedreht wird. Jedes Mal, wenn ein Messschenkel flach an einer Seite des Polygons liegt, bildet es ein antipodales Punktpaar, wobei der Punkt oder die Seite den entgegengesetzte Messschenkel berührt. Die vollständige Drehung des Messschiebers um das Polygon erkennt alle antipodalen Punktpaare und bestimmt das Punktepaar mit dem größten Abstand.

Um den Algorithmus auf beliebige Punkte in der Ebene anwenden zu können, muss erst die konvexe Hülle der Punkte bestimmt werden. Dafür wird die Laufzeit $O(n\cdot \log n)$ benötigt.

Zwei Punkte mit maximalem Abstand müssen auf den Rand des konvexen Polygons liegen, das die konvexe Hülle bildet. Der Algorithmus beginnt mit den Punkten $P_n$ und $P_1$ und berechnet den Flächeninhalt der Dreiecke $P_n{P}_1{P}_k$. Zunächst wird $k=2$ gesetzt und so lange erhöht, wie der Flächeninhalt zunimmt. So wird der antipodalen Punkt für $P_1$ bestimmt. Auf ähnliche Weise wird der antipodale Punkt für $P_2$ bestimmt wird, indem der Flächeninhalt der Dreiecke $P_1{P}_2{P}_k$ berechnet wird und der Index $k$ weiter erhöht wird, so lange der Flächeninhalt zunimmt. Diese Schritte werden mit den anderen Punkten $P_i$ wiederholt, bis der Index $k$ erreicht ist, also $i=k$ ist. Die Abstände der antipodalen Punktpaare werden berechnet und mit dem größten bisher gefundenen Abstand verglichen. Die Laufzeit für die Berechnung der Flächeninhalte und Abstände liegt in $O(n)$.

Insgesamt ergibt sich die Laufzeit $O(n\cdot \log n)+O(n)=O(n\cdot \log n)$.

Programmierung

Das folgende Computerprogramm in der Programmiersprache C++ zeigt eine Implementierung des Rotating calipers Algorithmus. Bei der Ausführung des Programms wird die Methode main verwendet, die das Punktepaar und den kleinsten Abstand zwischen dem Punktpaar das auf der Konsole ausgibt.

Die Funktion getConvexHull bestimmt die konvexe Hülle der Punkte mit dem Graham-Scan-Algorithmus. Die Funktion getOrientation bestimmt die Orientierung der Dreiecke mithilfe des Flächeninhalts.

Code-Schnipsel

#include <iostream> #include <sstream> #include <stack> #include <vector> using namespace std; // Deklariert den Datentyp für Punkte mit x-Koordinate und y-Koordinate in der Ebene struct Point { int x, y; }; Point lowestPoint; // Globale Variable für den Punkt mit der kleinsten y-Koordinate // Diese Funktion berechnet das Doppelte des Flächeninhalts des Dreiecks aus den drei Punkten. Wenn die Punkte im Uhrzeigersinn sind, ist das Vorzeichen positiv. Wenn die Punkte im Gegenuhrzeigersinn sind, ist das Vorzeichen negativ. int getArea(Point point1, Point point2, Point point3) { return ((point1.y - point2.y) * (point2.x - point3.x) - (point2.y - point3.y) * (point1.x - point2.x)); } // Diese Funktion bestimmt die Orientierung des drei Punkte. Wenn die Punkte im Uhrzeigersinn sind, wird der Wert 1 zurückgegeben. Wenn die Punkte im Gegenuhrzeigersinn sind, wird der Wert -1 zurückgegeben. Wenn die Punkte kollinear sind, wird der Wert 0 zurückgegeben. int getOrientation(Point point1, Point point2, Point point3) { int area = getArea(point1, point2, point3); // Aufruf der Funktion, die das Doppelte des Flächeninhalts mit Vorzeichen berechnet if (area > 0) // Wenn das Vorzeichen positiv ist, sind die Punkte im Uhrzeigersinn { return 1; } if (area < 0) // Wenn das Vorzeichen negativ ist, sind die Punkte im Gegenuhrzeigersinn { return -1; } return 0; // Wenn der Flächeninhalt gleich 0 ist, sind die Punkte kollinear } // Diese Funktion berechnet das Doppelte des Flächeninhalts des Dreiecks aus den drei Punkten int absoluteArea(Point point1, Point point2, Point point3) { return abs(getArea(point1, point2, point3)); // Aufruf der Funktion, die das Doppelte des Flächeninhalts mit Vorzeichen berechnet } // Diese Funktion berechnet den euklidischen Abstand zwischen zwei Punkten double getDistance(Point point1, Point point2) { return sqrt((point1.x - point2.x) * (point1.x - point2.x) + (point1.y - point2.y) * (point1.y - point2.y)); } // Vergleichsfunktion für das Sortieren der Punkte nach Orientierung im Uhrzeigersinn int compareByOrientation(const void* a, const void* b) { Point* point1 = (Point*)a, * point2 = (Point*)b; int orientation = getOrientation(lowestPoint, *point1, *point2); // Aufruf der Funktion, die die Orientierung der Punkte bestimmt if (orientation == 0) // Wenn die Punkte kollinear sind, werden die Abstände zum Referenzpunkt verglichen { return (int)(getDistance(*point1, lowestPoint) - getDistance(*point2, lowestPoint)); } if (orientation == 1) // Wenn die Punkte im Uhrzeigersinn sind, ist der erste Punkt größer { return 1; } return -1; // Wenn die Punkte im Gegenuhrzeigersinn sind, ist der zweite Punkt größer } // Diese Funktion bestimmt die konvexe Hülle der gegebenen Punkte mit dem Algorithmus Graham Scan vector<Point> getConvexHull(Point points[], int numberOfPoints) { // Bestimmt den Punkt mit der kleinsten y-Koordinate und seinen Index. Wenn die y-Koordinate gleich ist, wird der Punkt mit der Punkt mit der kleinsten x-Koordinate bestimmt. lowestPoint = { INT_MAX,INT_MAX }; int index = -1; for (int i = 0; i < numberOfPoints; i++) // for-Schleife, die alle Punkte durchläuft { if (points[i].y < lowestPoint.y) // Vergleicht die y-Koordinate der Punkte { lowestPoint.y = points[i].y; lowestPoint.x = points[i].x; index = i; } else { if (lowestPoint.y == points[i].y && points[i].x < lowestPoint.x) // Vergleicht die x-Koordinate der Punkte, wenn die y-Koordinate gleich ist { lowestPoint.x = points[i].x; index = i; } } } points[index] = points[0]; points[0] = lowestPoint; // Setzt den Punkt an den Anfang des Arrays qsort(points, numberOfPoints, sizeof(Point), compareByOrientation); // Sortiert die Punkte nach Orientierung im Uhrzeigersinn stack<Point> pointStack; // Stack für die konvexe Hülle pointStack.push(points[0]); // Fügt den ersten Punkt der konvexen Hülle hinzu int k = 1; while (k < numberOfPoints - 1 && getOrientation(points[0], points[k], points[k + 1]) == 0) // Bestimmt den nächsten Punkt, der nicht kollinear ist { k++; } pointStack.push(points[k]); // Fügt den zweiten Punkt der konvexen Hülle hinzu for (int i = k + 1; i < numberOfPoints; i++) // for-Schleife, die alle weiteren Punkte des Stack durchläuft und die konvexe Hülle erzeugt { Point point = pointStack.top(); pointStack.pop(); // Entfernt das oberste Element vom Stack while (getOrientation(pointStack.top(), point, points[i]) >= 0) // Bestimmt den nächsten Punkt, der mit den beiden obersten Elementen des Stack im Uhrzeigersinn ist { point = pointStack.top(); pointStack.pop(); // Entfernt das oberste Element vom Stack } pointStack.push(point); pointStack.push(points[i]); // Fügt den nächsten Punkt dem Stack hinzu } vector<Point> convexHull; // Vektor für die konvexe Hülle int n = pointStack.size(); for (int i = 0; i < n; i++) // for-Schleife, die alle Punkte des Stack durchläuft und dem Vektor hinzufügt { convexHull.push_back(pointStack.top()); // Fügt das oberste Element des Stack dem Vektor hinzu pointStack.pop(); } return convexHull; } // Diese Funktion gibt den größten Abstand der Punkte und Referenzen auf die zwei Punkte mit dem kleinsten Abstand zurück und verwendet den Rotating calipers Algorithmus double getLargestDistance(Point points[], Point& point1, Point& point2, int numberOfPoints) { vector<Point> convexHull = getConvexHull(points, numberOfPoints); // Aufruf der Funktion für die konvexe Hülle Point* hull = new Point[convexHull.size()]; // Referenz auf das Array für die konvexe Hülle int n = convexHull.size(); // Variable für die Anzahl der Punkte der konvexen Hülle for (int i = 0; i < n; i++) // for-Schleife, die das Array für die konvexe Hülle füllt { hull[i] = convexHull[i]; } if (n == 1) // Wenn die konvexe Hülle aus 1 Punkt besteht, ist der größte Abstand gleich 0 { return 0; } if (n == 2) // Wenn die konvexe Hülle aus 2 Punkten besteht, ist der größte Abstand der Abstand dieser 2 Punkte { return sqrt(getDistance(hull[0], hull[1])); } int k = 1; // Variable für den Index des ersten Punkts while (absoluteArea(hull[n - 1], hull[0], hull[(k + 1) % n]) > absoluteArea(hull[n - 1], hull[0], hull[k])) // while-Schleife, die den Punkt mit dem größte Abstand zum ersten und letzten Punkt der konvexe Hülle bestimmt { k++; } double maximumDistance = 0; int i = 0; int j = k; // Bestimmt die zwei Punkte der konvexe Hülle mit dem größten Abstand mit dem Rotating calipers Algorithmus while (i <= k && j < n) // while-Schleife, die den Index des ersten Punkts durchläuft { double currentDistance = getDistance(hull[i], hull[j]); // Ruft die Funktion für die Berechnung des euklidischen Abstands auf if (currentDistance > maximumDistance) // Prüft, ob der Abstand zwischen den aktuellen Punkten größer als der bisher größte Abstand ist { point1 = hull[i]; point2 = hull[j]; maximumDistance = currentDistance; } while (j < n && absoluteArea(hull[i], hull[(i + 1) % n], hull[(j + 1) % n]) > absoluteArea(hull[i], hull[(i + 1) % n], hull[j])) // while-Schleife, die den Index des zweiten Punkts durchläuft { double currentDistance = getDistance(hull[i], hull[(j + 1) % n]); // Ruft die Funktion für die Berechnung des euklidischen Abstands auf if (currentDistance > maximumDistance) // Prüft, ob der Abstand zwischen den aktuellen Punkten größer als der bisher größte Abstand ist { point1 = hull[i]; point2 = hull[(j + 1) % n]; maximumDistance = currentDistance; } j++; } i++; } return maximumDistance; } // Diese Funktion gibt den Punkt in der Form (x, y) zurück string pointToString(Point* point) { stringstream text; text << "(" << point->x << ", " << point->y << ")" << endl; return text.str(); // Typumwandlung von stringstream nach string } // Hauptfunktion die das Programm ausführt int main() { Point points[] = { {4, 0}, {0, 2}, {-1, -7}, {1, 10}, {2, -3} }; // Deklariert und initialisiert ein Array mit 5 Punkten int numberOfPoints = sizeof(points) / sizeof(points[0]); Point point1, point2; // Deklariert Referenzen für die zwei Punkte mit dem größten Abstand cout << "Der größte Abstand zwischen zwei Punkten ist " << getLargestDistance(points, point1, point2, numberOfPoints) << endl; // Aufruf der Funktion, Ausgabe auf der Konsole cout << "Erster Punkt: " << pointToString(&point1); // Ausgabe auf der Konsole cout << "Zweiter Punkt: " << pointToString(&point2); // Ausgabe auf der Konsole }

Literatur

• Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, und Clifford Stein. Algorithmen – Eine Einführung. Oldenbourg-Verlag, 2004. ISBN 3-486-27515-1. Seiten 959–963 (Kapitel 33.4: Berechnung des dichtesten Punktepaares).
• Jon Kleinberg, Éva Tardos. Algorithm Design. Pearson International Edition, 2006. ISBN 0-321-37291-3. Seiten 225ff (Divide & Conquer); Seiten 741ff (Randomisierter Algorithmus).