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.

 

 

konvexe Hülle um eine Punktemenge
konvexe Hülle um eine Punktemenge

 

 

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