Dichtestes Punktpaar
# 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 ${\tbinom {n}{2}}={\tfrac {n\cdot (n-1)}{2}}={\tfrac {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

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:

* Das Punktpaar liegt in der linken oder der rechten Hälfte. In diesem Fall gibt es keine weiteren Probleme.
* 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. 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 ${\tfrac {\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({\tfrac {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 linear zur Anzahl der Punkte.

### 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({\tfrac {\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:

**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 ${\tfrac {\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**

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({\tfrac {\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 {\tfrac {\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 ${\tfrac {\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 ${\tfrac {\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 ${\tfrac {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}\longrightarrow \{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 {\tfrac {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
