C++で実装するグラハムスキャンアルゴリズムによる凸包の求め方

凸包(とつほう)とは、平面上に与えられたすべてのデータ点を包含できる最小の閉領域のことです。

グラハムスキャン(Graham Scan)アルゴリズムは、この凸包の境界となる頂点を効率的に求める古典的な手法です。まず、最も下にある点(y座標が最小の点)を選び、それを凸包の開始点とします。続いて、残りのn-1個の頂点を、開始点から見た反時計回りの角度に基づいてソートします。同じ角度を持つ点が複数ある場合は、開始点から最も遠い点を除いてすべて取り除きます。

その後、残りの点を順にスタックへプッシュしていきます。スタックの最上位の点・その下の点・新しく選択した点points[i]の3点が反時計回りになっていない場合は、スタックから要素を1つずつポップして除去し、方向が確認できた時点でpoints[i]をスタックに挿入します。

なお、グラハムスキャンの計算量はO(n log n)であり、支配的なのは点群のソート処理です。

入力: 点の集合 {(-7,8), (-4,6), (2,6), (6,4), (8,6), (7,-2), (4,-6), (8,-7), (0,0), (3,-2), (6,-10), (0,-6), (-9,-5), (-8,-2), (-8,0), (-10,3), (-2,2), (-10,4)}
出力: 凸包の境界点は (-9, -5) (-10, 3) (-10, 4) (-7, 8) (8, 6) (8, -7) (6, -10)

アルゴリズム

findConvexHull(points, n)

入力: 点の集合と点の個数

出力: 凸包の境界上の点

Begin
    minY := points[0].y
    min := 0
    for i := 1 to n-1 do
        y := points[i].y
        if y < minY または (minY == y かつ points[i].x < points[min].x) then
            minY := points[i].y
            min := i
    done
    points[0] と points[min] を交換
    p0 := points[0]
    points[1] 以降を反時計回りの角度でソート
    arrSize := 1
    for i := 1 to n do
        i < n-1 かつ (p0, points[i], points[i+1]) が同一直線上にある間
            i := i + 1
        done
        points[arrSize] := points[i]
        arrSize := arrSize + 1
    done
    if arrSize < 3 then
        return cHullPoints
    points[0]、points[1]、points[2] をスタックにプッシュ
    for i := 3 to arrSize do
        スタック最上位の点・その下の点・points[i] が反時計回りになっていない間
            スタックの最上位要素を削除
        done
        points[i] をスタックにプッシュ
    done
    スタックが空になるまで
        cHullPoints にスタック最上位の要素を追加
        スタックからポップ
    done
End

処理手順のポイント

• y座標が最小の点(同値の場合はx座標も最小)を基準点p0として選択する
• 残りの点をp0からの反時計回りの角度でソートする
• 同一角度の点は、最も遠い点のみを残して重複を除去する
• 有効な点が3未満の場合は凸包が成立しないため、空の結果を返す
• スタックを用いて常に「左回り(反時計回り)」を維持しながら頂点を追加・除去する

サンプルコード(C++)

#include<iostream>
#include<stack>
#include<algorithm>
#include<vector>
using namespace std;

struct point { // 2次元平面上の点を定義
    int x, y;
};

point p0; // 他の2点との比較に使用する基準点

// スタックの上位から2番目の要素を取得する
point secondTop(stack<point> &stk) {
    point tempPoint = stk.top();
    stk.pop();
    point res = stk.top(); // 2番目の要素を取得
    stk.push(tempPoint);   // 元の最上位要素を戻す
    return res;
}

// 2点間の距離の2乗を返す
int squaredDist(point p1, point p2) {
    return ((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y));
}

// 3点の回転方向を判定する
int direction(point a, point b, point c) {
    int val = (b.y-a.y)*(c.x-b.x)-(b.x-a.x)*(c.y-b.y);
    if (val == 0)
        return 0; // 同一直線上(コリニア)
    else if (val < 0)
        return 2; // 反時計回り
    return 1;     // 時計回り
}

// qsort用の比較関数
int comp(const void *point1, const void *point2) {
    point *p1 = (point*)point1;
    point *p2 = (point*)point2;
    int dir = direction(p0, *p1, *p2);
    if (dir == 0)
        return (squaredDist(p0, *p2) >= squaredDist(p0, *p1)) ? -1 : 1;
    return (dir == 2) ? -1 : 1;
}

// 凸包を求める本体
vector<point> findConvexHull(point points[], int n) {
    vector<point> convexHullPoints;
    int minY = points[0].y, minIndex = 0;
    for (int i = 1; i < n; i++) {
        int y = points[i].y;
        // 最も下にある点(同値の場合は最も左の点)を探す
        if ((y < minY) || (minY == y && points[i].x < points[minIndex].x)) {
            minY = points[i].y;
            minIndex = i;
        }
    }
    swap(points[0], points[minIndex]); // 最小点を先頭へ移動
    p0 = points[0];
    qsort(&points[1], n-1, sizeof(point), comp); // 残りの点を極角でソート

    int arrSize = 1;
    for (int i = 1; i < n; i++) {
        // 角度が同じ点は除去し、最も遠い点だけを残す
        while (i < n-1 && direction(p0, points[i], points[i+1]) == 0)
            i++;
        points[arrSize] = points[i];
        arrSize++;
    }

    if (arrSize < 3)
        return convexHullPoints; // 点が3未満の場合は空のリストを返す

    // スタックを作成し、最初の3点を追加
    stack<point> stk;
    stk.push(points[0]);
    stk.push(points[1]);
    stk.push(points[2]);

    for (int i = 3; i < arrSize; i++) { // 残りの頂点を処理
        // 左回りにならない限り、スタックの要素を除去
        while (direction(secondTop(stk), stk.top(), points[i]) != 2)
            stk.pop();
        stk.push(points[i]);
    }

    while (!stk.empty()) {
        convexHullPoints.push_back(stk.top()); // スタックから点を取り出す
        stk.pop();
    }
    return convexHullPoints;
}

int main() {
    point points[] = {{-7,8},{-4,6},{2,6},{6,4},{8,6},{7,-2},{4,-6},{8,-7},{0,0},
                      {3,-2},{6,-10},{0,-6},{-9,-5},{-8,-2},{-8,0},{-10,3},{-2,2},{-10,4}};
    int n = 18;
    vector<point> result;
    result = findConvexHull(points, n);
    cout << "Boundary points of convex hull are: " << endl;
    vector<point>::iterator it;
    for (it = result.begin(); it != result.end(); it++)
        cout << "(" << it->x << ", " << it->y << ") ";
}

実行結果

Boundary points of convex hull are:
(-9, -5) (-10, 3) (-10, 4) (-7, 8) (8, 6) (8, -7) (6, -10)