Konvexe Hülle um eine Punktemenge

Die konvexe Hülle (engl. convex hull) ist ein Polygon, das eine Punktemenge einschließt. Konvex ist ein Ausdruck aus der Mathematik. In Verbindung mit der Hülle lässt er sich sehr gut veranschaulichen (s. auch Bild unten)

Man stelle sich vor, die Punkte seien Nägel, die in ein Brett geschlagen wurden. Spannt man nun ein Gummiband um die Gesamtheit der Nägel, so ist die Linie des Bandes, die sich daraus ergibt, die konvexe Hülle
Eine andere Beschreibung, die auch bei der Ermittlung der konvexen Hülle genutzt wird: Startet man bei irgendeinem Punkt der Hülle und folgt der Linie gegen den Uhrzeigersinn, so darf man auf diesem Weg niemals rechts abbiegen.

Die Graham-Methode zur Ermittlung der konvexen Hülle

Algorithmen zur Bestimmung der konvexen Hülle gibt es viele. Die bekanntesten sind die Graham-Methode (Graham Scan) und Jarvis‘ March. Hier wird der Graham Scan benutzt, ein sehr eleganter Algorithmus, der einen Stapel (stack) benutzt, um die Punkte der Hülle zusammen zu stellen. Außerdem wird eine Sortierfunktion gebraucht.

Das Vorgehen in sehr grober Form (in Anlehnung an Cormen et al.: Introduction to Algorithms):

1. Bestimme den untersten Punkt p(0) (falls nicht eindeutig, dann von den untersten den am weitesten links liegenden).

2. Ordne die restlichen Punkte P(i) nach dem Winkel zwischen der Horizontalen durch P(0) und der Verbindung zwischen P(0) und P(i) (sog. Polarwinkel). Gibt es mehrere Punkte mit demselben Winkel, so nimm nur den am weitesten von P(0) entfernten.

3. Durchlaufe die Punkte P(0), P(1), …, P(n) in dieser Reihenfolge und sammle die Punkte der Hülle nach der folgenden Regel: wird auf dem Weg von P(i) nach P(k) über P(j) in P(j) nicht links abgebogen, so gehört P(j) nicht zur Hülle.

Die vorbereitende Sortierung (Schritt 2) ist im Vergleich zum eigentlichen Scan (Schritt 3) aufwändiger. Schritt 3 ist mit Hilfe eines Stapels ohne großen Aufwand durchzuführen. Hier der detaillierte Ablauf in Pseudocode:

Beginne mit einem leeren Stapel S

Push (S, P(0))

Push (S, P(1))

Push (S, P(2))

For i = 3 to n

      Do While (der von NextToTop(S), Top(S) und P(i) gebildete Winkel
                                                                               ist nicht nach links gerichtet)
             Pop(S)

Loop

Push (S, P(i))

Next i

Die Funktion convexHull_GS

Die Ermittlung der konvexen Hülle wurde als Funktion convexHull_GS programmiert, die sich in einem Modul convexHull befindet. Die Funktion bezieht alle notwendigen Informationen über ihre Parameter und benutzt keine globalen Variablen. Sie ist also sehr gut wiederverwendbar. Ihre Kopfzeile lautet:

Public Function convexHull_GS (ByRef p() as Double) As Double()

Das Array p ist zweidimensional und 1-basiert. Es enthält die Koordinaten der Punkte, die umhüllt werden sollen. Die Funktion liefert die Folge der Hüllenpunkte ebenfalls in Form eines zweidimensionalen, 1-basierten Arrays. Dabei ist der Startpunkt zweimal enthalten, damit der Kreis geschlossen wird.

Aus der Funktion werden zwei kleine Hilfsfunktionen aufgerufen, die sich in demselben Modul (convexHull) befinden:

Funktion dist ermittelt die Entfernung zwischen zwei Punkten nach dem Satz von Pythagoras (euklidische Distanz)
Funktion isLeft prüft für drei Punkte p1, p2 und p3, ob auf dem Weg p1-p2-p3 links abgebogen wird.

Außerdem bedient sich convexHull_GS eines Stapels (Klasse stack), der neben den gängigen Stapel-Methoden die Methoden nextToTop und stackArray bietet. Mit nextToTop wird das zweite Element von oben gelesen, ohne dass ein Element aus dem Stapel entfernt würde. Mit der Methode stackArray wird der gesamte Stapel in einem Array geliefert, wobei alle Elemente auf dem Stapel verbleiben.

Die vorbereitende Sortierung der Punkte (nach steigenden Polarwinkeln und, innerhalb gleicher Polarwinkel nach steigender Entfernung von P(0)) wird an eine Funktion matrSort2 im Modul sort delegiert.

Code des Moduls convexHull

Da die Funktion ConvexHull_GS keine globalen Variablen benötigt, ist es egal, ob man sie in einem Standardmodul unterbringt oder in einem Klassenmodul. Der Code der Funktion ist in beiden Fällen derselbe. Man kann die Wahl der Modulart davon abhängig machen, wie die gesamte Anwendung aufgebaut ist.

Bei der Programmierung der Funktion convexHull_GS wurde darauf verzichtet, das letzte Quäntchen an Effizienz heraus zu pressen. Stattdessen wurde mehr darauf geachtet, die einzelnen Schritte deutlich zu zeigen. Hier ist der Code:

'produces a list of points representing the convex hull for the

'set of points p by use of the Graham scan method

'-------------------------------------------------------------------------------

Public Function convexHull_GS(ByRef p() As Double) As Double()

Dim i As Long, j As Integer, k As Long

Dim numPts As Long, numNpts As Long, minPt As Long

Dim h() As Double 'points, without p0, with direction cosines & distance from p0

Dim s() As Double 'points sorted according to direction cosine & distance from p0

Dim t() As Double 'points, unnecessary ones omitted

Dim st As stack

numPts = UBound(p, 1)

    'search starting point (point having min x within min y
    minPt = 1
    For i = 2 To numPts
        If p(i, 2) < p(minPt, 2) Then
            minPt = i
        Else
            If p(i, 2) = p(minPt, 2) And p(i, 1) < p(minPt, 1) Then minPt = i
        End If
    Next i

'copy points (without starting point) to h;

'add direction cosines and distances from starting point

ReDim h(1 To numPts - 1, 1 To 4)

j = 1

For i = 1 To numPts

If i <> minPt Then

h(j, 1) = p(i, 1)

h(j, 2) = p(i, 2)

'calculating negative of direction cosine

h(j, 3) = -((p(i, 1) - p(minPt, 1)) / _

((p(i, 1) - p(minPt, 1)) ^ 2 + (p(i, 2) - p(minPt, 2)) ^ 2) ^ 0.5)

'add distance from starting point

h(j, 4) = Dist(p(minPt, 1), p(minPt, 2), p(i, 1), p(i, 2))

j = j + 1

End If

Next i

'sort according to negative of direction cosine and distance

'from starting point and transfer to s

s = sort.matrSort2(h, 3, 4)

'eliminate all but one point with identical cosines

'by copying selectively to t

ReDim t(0 To numPts - 1, 1 To 3) 'index 0 is for starting point

Dim copy As Boolean

numNpts = 0

i = 1

Do While i <= numPts - 1

If i = numPts - 1 Then

copy = True

Else

If s(i, 3) <> s(i + 1, 3) Then copy = True Else copy = False

End If

If copy Then

numNpts = numNpts + 1

For j = 1 To 3

t(numNpts, j) = s(i, j)

Next j

End If

i = i + 1

Loop

t(0, 1) = p(minPt, 1) 'insert starting point

t(0, 2) = p(minPt, 2)

'scan and hold back the points of the hull in the stack

'the elements of the stack are indices referring to array t

Set st = New stack

st.initStack

st.push (0)

st.push (1)

st.push (2)

For i = 3 To numNpts

'remove from stack points that are not part of the hull

Do While Not isLeft(t(st.nextToTop, 1), t(st.nextToTop, 2), _

t(st.top, 1), t(st.top, 2), t(i, 1), t(i, 2))

st.pop

Loop

'push the current point on the stack

st.push (i)

Next i

'get stack contents in total; copy its elements from t to d

Dim sa() As Variant

sa = st.stackArray 'indices of the points in the hull

Dim d() As Double 'coordinates of the points in the hull

ReDim d(1 To UBound(sa) + 1, 1 To 2)

For i = 1 To UBound(d, 1) - 1

For j = 1 To UBound(d, 2)

d(i, j) = CDbl(t(CLng(sa(i)), j))

Next j

Next i

d(UBound(d, 1), 1) = t(0, 1) 'add starting point a second time

d(UBound(d, 1), 2) = t(0, 2) 'in order to close the circle

convexHull_GS = d

End Function

Die beiden Hilfsfunktionen Dist und isLeft sind, wie bereits erwähnt, ebenfalls im Modul enthalten und konnten deshalb aus convexHull_GS ohne Nennung des Moduls aufgerufen werden. Ihr Code lautet:

'takes the coordinates of three points p1, p2, p3, and checks if

'we turn to the left in p2 when following the path p1p2p3

‘---------------------------------------------------------------------------------

Private Function isLeft(ByVal x1 As Double, ByVal y1 As Double, _

ByVal x2 As Double, ByVal y2 As Double, _

ByVal x3 As Double, ByVal y3 As Double) As Boolean

isLeft = (x3 - x2) * (y2 - y1) - (y3 - y2) * (x2 - x1) < 0

End Function

'calculates euclidian distance between (x1,y1) and (x2,y2)

'---------------------------------------------------------------------

Private Function Dist(ByVal x1 As Double, ByVal y1 As Double, _

ByVal x2 As Double, ByVal y2 As Double) As Double

Dist = ((x2 - x1) ^ 2 + (y2 - y1) ^ 2) ^ 0.5

End Function

Der Code der Funktion matrSort2

Der oben dargestellte Code der Funktion convexHull_GS beruht auf der Annahme, dass die Funktion matrSort2 in einem Standardmodul namens sort untergebracht ist. Dies ist natürlich nicht zwingend, man könnte auch eine Klasse dafür verwenden. In diesem Fall müsste convexHull_GS entsprechend angepasst werden.

Die Funktion matrSort2 sortiert die Zeilen einer Matrix aufsteigend nach den Werten zweier Spalten. Beim Aufruf der Funktion werden die Nummern dieser Spalten übergeben, wobei die Spaltenzählung bei 1 beginnt. Die Zeilen der Matrix werden zunächst nach den Werten der zuerst genannten Spalte c1 sortiert (primäres Sortierkriterium), innerhalb gleicher Werte von c1 dann nach Werten der zweiten genannten Spalte c2.

'sorts the rows of matrix m according to the values in columns c1

'and c2 (for identical values of c1) using a variation of selection sort;

'smallest values first; indices begin with 1

'-----------------------------------------------------------------------

Public Function matrSort2(ByRef m() As Double, _

ByVal c1 As Long, _

ByVal c2 As Long) As Double()

Dim i As Long, j As Long, minPos As Long

Dim numr As Long, numc As Long 'number of rows, columns

Dim s() As Double

Dim h() As Boolean

numr = UBound(m, 1)

numc = UBound(m, 2)

ReDim s(1 To numr, 1 To numc) 'the result

ReDim h(1 To numr) 'tells whether a certain row is still in m

For i = 1 To numr 'in the beginning, all rows are still in m

h(i) = True

Next i

For i = 1 To numr 'for all rows in the result array

'search initial value for minPos

Dim found As Boolean

found = False

minPos = 0

Do While Not found

minPos = minPos + 1

If h(minPos) Then found = True

Loop

        'search position of smallest row in m
        For j = 1 To numr
          If h(j) Then
              If m(j, c1) < m(minPos, c1) Then
                  minPos = j
              Else
                  If m(j, c1) = m(minPos, c1) And m(j, c2) < m(minPos, c2) Then minPos = j
              End If
          End If
        Next j

'transfer smallest row from m to s

For j = 1 To numc

s(i, j) = m(minPos, j)

h(minPos) = False

Next j

Next i

matrSort2 = s

End Function

Die Klasse stack

wurde auf der Basis eines Variant-Arrays implementiert. Wie bereits erwähnt, wurde die Klasse um zwei Methoden erweitert, die nicht zur Standardausstattung eines Stapels gehören, nämlich nextToTop und stackArray. Die beiden Methoden wären auch für den Graham Scan nicht unbedingt nötig gewesen, erleichtern aber die Programmierung sehr.

Private s() As Variant

Private tos As Integer

Public Sub initStack()

tos = 0

ReDim s(1 To 1)

End Sub

Public Sub push(ByVal e As Variant)

tos = tos + 1

ReDim Preserve s(1 To tos)

s(tos) = e

End Sub

Public Function pop() As Variant

pop = s(tos)

tos = tos - 1

End Function

Public Function top() As Variant

top = s(tos)

End Function

Public Function nextToTop() As Variant

nextToTop = s(tos - 1)

End Function

Public Function isEmpty() As Boolean

isEmpty = tos = 0

End Function

Public Function stackArray() As Variant

Dim a() As Variant

Dim i As Integer

ReDim a(1 To tos)

For i = 1 To tos

a(i) = s(i)

Next i

stackArray = a

End Function

Public Function getTos() As Integer

getTos = tos

End Function