Convex Hull

Problemstellung

Eine konvexe Hülle einer Menge von Punkten ist das kleinste konvexe Polygon, welches alle Punkte enthält. Konvex bedeutet, dass es keinen Eckpunkt gibt, bei dem das Polygon nach innen gebogen ist.

Konvexe Fläche in der linken Abbildung, nicht konvexe Fläche in der rechten.

Die Fläche in der rechten Abbildung ist somit nicht konvex.

Die konvexe Hülle einer endlichen Menge von Punkten kann man sich mithilfe der Gummiband-Analogie vorstellen. Angenommen, die Punkte der Menge sind Nägel, die aus einem Brett hervorstehen, kann man ein Gummiband um alle Nägel spannen, welches versuchen wird, den kleinstmöglichen Umfang einzunehmen. Der Pfad des Gummibandes entspricht nun der Umrisslinie der konvexen Hülle.

Gummmiband-Analogie, Bild von `Wikipedia <https://en.wikipedia.org/wiki/Convex_hull>`__

Lösungsansatz

Um die konvexe Hülle zu finden bietet sich der Graham Scan an, der oft als einer der ersten Algorithmen der Computergeometrie bezeichnet wird.

Das gesuchte Polygon lässt sich als eine Abfolge einer Teilmenge der Punktmenge ausdrücken. Der Graham Scan lässt uns diese Abfolge von einem Anfangspunkt aus im Gegenuhrzeigersinn suchen.

Der Algorithmus beginnt damit, einen Punkt zu suchen, der garantiert auf der konvexen Hülle liegt. Dazu kann man zum Beispiel den Punkt mit der kleinsten y Koordinate verwenden, falls mehrere solche existieren, nehmen wir aus dieser Menge den Punkt mit der grössten x Koordinate. Es sind auch andere Wahlen möglich, zum beispiel der Punkt ganz links unten, wobei andere Optionen kleine Vor- und Nachteile haben. Dies ist der persönlichen Präferenz des Programmieres überlassen.

Anfangspunkt

Der gewählte Anfangspunkt ist hier in rot markiert.

Nun sortieren wir die anderen Punkte nach dem Winkel, den sie mit dem Anfangspunkt machen. Falls zwei Punkte den gleichen Winkel mit dem Anfangspunkt bilden, somit kollinear sind mit diesem, soll der Punkt der näher zum Anfangspunkt liegt früher im sortierten Array erscheinen.

Sortieren nach Winkel

Sortieren der anderen Punkte nach Winkel

Wenn wir nun den Punkten im Gegenuhrzeigersinn entlang gehen, können wir bei jedem Tripel entscheiden, ob der mittlere Punkt des Tripels einen Knick nach innen oder nach aussen repräsentiert. Falls der Knick nach innen ist, kann der mittlere Punkt nicht Teil der konvexen Hülle sein. Wir speichern die Punkte, die momentan Teil der Hülle sind auf einem Stack, somit können wir Objekte hinzufügen und entfernen. Dabei sind immer die ersten beiden Punkte des Tripels bereits auf dem Stack, den dritten betrachten wir gerade.

triplet test 1

Knick nach aussen.

Wir betrachten den Knick zwischen den roten Linien. Mithilfe des Kreuzproduktes können wir herausfinden, dass der Knick sich im Gegenuhrzeigersinn dreht, also nach aussen ist. Somit wird der dritte dieser Punkte zum Stack hinzugefügt.

Somit gehen wir also weiter zum nächsten Punkt, und überprüfen wieder, ob die Bedingung eingehalten wird.

triplet test 2

Knick nach innen, Bedingung wird also nicht erfüllt.

Hier dreht sich die Linie im Uhrzeigersinn, der Mittlere Punkt kann also nicht Teil der konvexen Hülle sein und wir entfernen diesen Punkt vom Stack.

triplet test 3

Beim nächsten Test wird der Punkt, den wir gerade aus dem Stack gelöscht haben ignoriert.

So können wir nun über alle Punkte loopen, bis wir die gesammte Hülle gefunden haben.

Korrektheit

Da wir bei jedem Punkt überprüft haben, ob der Winkel konvex ist, wissen wir, dass das resultierende Polygon konvex sein muss. Wir haben also ein konvexes Polygon, dessen Eckpunkte eine Teilmenge unserer Punkte sind. Wir müssen also nur noch beweisen, dass sich keine Punkte unserer Punktemenge ausserhalb dieses Polygons befinden.

(Zu beweisen, dass die ersten drei, so wie die letzten drei Punkte des Polygons nicht kollinear sind, ist als Übung dem Leser überlassen. Hint: betrachte die Funktion, die wir verwenden um die Punkte nach Winkel zu sortieren.)

Angenommen es existiert ein Punkt $P$ ausserhalb des Polygons, welches wir soeben gefunden haben. Wir wissen, dass dieser nicht der erste und auch nicht der letzte Punkt in der Sortierung sein kann, da diese sicher auf der Hülle liegen, wie einfach zu zeigen ist.

illegal point P

Der rote Punkt :math:`P` liegt ausserhalb des konvexen Polygons.

Da jeder Punkt einmal zum Stack hinzugefügt wird, wissen wir, dass $P$ einmal aus dem Stack entfernt wurde. Dies heisst, dass die beiden benachbarten Punkte zusammen mit $P$ einen Winkel im Uhrzeigersinn bilden. Die beiden benachbarten Punkte $O$ und $Q$ liegen auf der Grenze oder innerhalb des Polygons. Somit müsste $P$, da der Winkel nach rechts geht innerhalb des Polygons liegen, da das Polygon ja konvex ist. Konvex impliziert, dass alle Punkte auf der Strecke zwischen den Punkten $O$ und $Q$ im Polygon (oder auf dessen Grenze) liegen müssen, da $O$ und $Q$ im Polygon liegen.

Somit haben wir einen Wiederspruch erreicht, und der Punkt $P$ kann nicht existieren.

Implementation

Die hier implementierte Funktion convex_hull berechnet aus einer Menge von Punkten einen Array der Punkte, die auf der konvexen Hülle liegen. Die Punkte im zurückgegebenen Array sind im Gegenuhrzeigersinn sortiert.

In der Praxis ist es einfacher, nicht einen stack zu verwenden um die Punkte zu speichern, die momentan Teil der konvexen Hülle sind, sondern einen vector, denn damit können wir auch auf das zweitoberste Objekt zugreifen.

Punkte werden als pair zweier int Werte gespeichert, also als {x,y}.

#include<vector>
#include<cmath>
#include<algorithm>

#define point pair<int,int>

vector<point> convex_hull(vector<point> points){
  int n = points.size();

  // Triviale Eingaben gleich zu Begin ausschliessen.
  if(n<3) return points;

  // Wir suchen den Punkt unten rechts und speichern diesen in points[0]
  // so wie in der globalen Variable origin.
  int lowest = 0;
  for(int i = 1; i < n; ++i){
    if(make_pair(points[i].second, -points[i].first) <
       make_pair(points[lowest].second, -points[lowest].first)){
      lowest = i;
    }
  }
  swap(points[0], points[lowest]);
  origin = points[0];

  // Die restlichen Punkte (ohne origin) werden sortiert nach unserer
  // eigenen Sortierfunktion compareAngles
  sort(points.begin()+1, points.end(), compareAngles);

  // st agiert als stack um die Punkte der convex hull zu finden.
  vector<point> st;
  st.push_back(points[0]);
  st.push_back(points[1]);

  int pointer = 2;
  while(pointer < n){
    point a = st[st.size()-2];
    point b = st[st.size()-1];
    point c = points[pointer];
    if(convex(a,b,c)){
      // Die Punkte bilden einen Knick nach aussen. (ccw)
      st.push_back(c);
      ++pointer;
    }else{
      // Der mittlere Punkt ist nicht Teil der Hülle.
      st.pop_back();
      if(st.size()<2){
        st.push_back(c);
        ++pointer;
      }
    }
  }

  // Falls drei kollineare Punkte am Ende des Rundgangs um
  // die Hülle sind, entfernen wir hier den mittleren.
  point a = st[st.size()-2];
  point b = st[st.size()-1];
  point c = st[0];
  if(!convex(a,b,c)){
    st.pop_back();
  }

  return st;
}

Die Helferfunktionen die dazu verwendet werden, sind eine Funktion crossproduct um das Kreuzprodukt zweier Vektoren zu finden, so wie eine Funktion convex, die drei Punkte als Argument akzeptiert und angibt, ob der Pfad aus diesen Punkten sich im Gegenuhrzeigersinn biegt.

Die funktion compareAngles, die wir oben verwendet haben um den vector von Punkten damit zu sortieren, verwendet das Kreuzprodukt, um die Winkel zweier Punkte zum Punkt origin zu vergleichen, und muss so nie die Winkel explizit ausrechen. Betrachte auch, dass wir den Punkt origin deswegen als globale Variable speichern müssen.

Falls die beiden Punkte den gleichen Winkel zu origin bilden, vergleichen wir die Distanzen zu origin. Dazu reicht die Manhattan Distanz aus, dies macht die Berechnung einfacher.

point origin;

long long crossproduct(point a, point b){
  return a.first * b.second - a.second * b.first;
}

// Ist der Winkel zwischen den Vektoren ab und bc im Gegenuhrzeigersinn?
bool convex (point a, point b, point c){
  point vA = {b.first - a.first, b.second - a.second};
  point vB = {c.first - b.first, c.second - b.second};

  return crossproduct(vA, vB) > 0;
}

bool compareAngles(const point &a, const point &b){
  if(convex(origin, a, b)){
    return true;
  }
  if(convex(origin, b, a)){
    return false;
  }
  else{
    // die beiden Punkte sind kollinear mit dem origin, wir müssen
    // also die Distanze zum origin vergleichen.
    int distA = abs(a.first-origin.first) + abs(a.second-origin.second);
    int distB = abs(b.first-origin.first) + abs(b.second-origin.second);

    return distA < distB;
  }
}

Laufzeit- und Speicheranalyse