C++で分割統治法を使って凸包を求める方法

このチュートリアルでは、与えられた点の集合に対して凸包(Convex Hull)を求めるプログラムについて解説します。

凸包とは、与えられたすべての点を「境界上」または「内部」に含む最小の凸多角形のことです。計算幾何学における基本的な問題の一つで、画像処理や衝突判定、地理情報システムなど幅広い分野で応用されています。

アルゴリズムの概要

本プログラムでは分割統治法(Divide and Conquer)を採用しています。全体の流れは以下のとおりです。

1. 点の集合を x 座標でソートした後、左右半分に分割する
2. 各半分について再帰的に凸包を求める
3. 点数が少ない場合(5点以下)は総当たり(ブルートフォース)で直接凸包を計算する
4. 2つの凸包を上下の共通接線(タンジェント)で結びつつマージし、最終的な凸包を構築する

この手法により、全体の計算量は O(n log n) に抑えられます。

C++での実装例

#include<bits/stdc++.h>
using namespace std;
// 多角形の中心点を格納
pair<int, int> mid;
// 特定の点が属する象限を計算
int quad(pair<int, int> p){
   if (p.first >= 0 && p.second >= 0)
   return 1;
   if (p.first <= 0 && p.second >= 0)
   return 2;
   if (p.first <= 0 && p.second <= 0)
   return 3;
   return 4;
}
// 線が多角形に対してどちら側にあるかを判定
int calc_line(pair<int, int> a, pair<int, int> b,
pair<int, int> c){
   int res = (b.second-a.second)*(c.first-b.first) -
   (c.second-b.second)*(b.first-a.first);
   if (res == 0)
   return 0;
   if (res > 0)
   return 1;
   return -1;
}
// ソート用の比較関数
bool compare(pair<int, int> p1, pair<int, int> q1){
   pair<int, int> p = make_pair(p1.first - mid.first,
   p1.second - mid.second);
   pair<int, int> q = make_pair(q1.first - mid.first,
   q1.second - mid.second);
   int one = quad(p);
   int two = quad(q);
   if (one != two)
   return (one < two);
   return (p.second*q.first < q.second*p.first);
}
// 2つの多角形の上部接線を求めてマージする
vector<pair<int, int>> merger(vector<pair<int, int> > a,
vector<pair<int, int> > b){
   int n1 = a.size(), n2 = b.size();
   int ia = 0, ib = 0;
   // a の右端の点を求める
   for (int i=1; i<n1; i++)
   if (a[i].first > a[ia].first)
   ia = i;
   // b の左端の点を求める
   for (int i=1; i<n2; i++)
   if (b[i].first < b[ib].first)
   ib=i;
   int inda = ia, indb = ib;
   bool done = 0;
   // 上部接線の計算
   while (!done){
      done = 1;
      while (calc_line(b[indb], a[inda], a[(inda+1)%n1]) >=0)
      inda = (inda + 1) % n1;
      while (calc_line(a[inda], b[indb], b[(n2+indb-1)%n2]) <=0){
         indb = (n2+indb-1)%n2;
         done = 0;
      }
   }
   int uppera = inda, upperb = indb;
   inda = ia, indb=ib;
   done = 0;
   int g = 0;
   // 下部接線の計算
   while (!done){
      done = 1;
      while (calc_line(a[inda], b[indb], b[(indb+1)%n2])>=0)
      indb=(indb+1)%n2;
      while (calc_line(b[indb], a[inda], a[(n1+inda-1)%n1])<=0){
      inda=(n1+inda-1)%n1;
      done=0;
   }
}
int lowera = inda, lowerb = indb;
vector<pair<int, int>> ret;
// 2つの多角形をマージして凸包を得る
int ind = uppera;
ret.push_back(a[uppera]);
while (ind != lowera){
   ind = (ind+1)%n1;
   ret.push_back(a[ind]);
}
ind = lowerb;
ret.push_back(b[lowerb]);
while (ind != upperb){
   ind = (ind+1)%n2;
   ret.push_back(b[ind]);
}
return ret;
}
// 総当たりで凸包を求める関数
vector<pair<int, int>> bruteHull(vector<pair<int, int>> a){
   set<pair<int, int> >s;
   for (int i=0; i<a.size(); i++){
      for (int j=i+1; j<a.size(); j++){
         int x1 = a[i].first, x2 = a[j].first;
         int y1 = a[i].second, y2 = a[j].second;
         int a1 = y1-y2;
         int b1 = x2-x1;
         int c1 = x1*y2-y1*x2;
         int pos = 0, neg = 0;
         for (int k=0; k<a.size(); k++){
            if (a1*a[k].first+b1*a[k].second+c1 <= 0)
            neg++;
            if (a1*a[k].first+b1*a[k].second+c1 >= 0)
            pos++;
         }
         // 全ての点が直線の同じ側にあれば辺として採用
         if (pos == a.size() || neg == a.size()){
            s.insert(a[i]);
            s.insert(a[j]);
         }
      }
   }
   vector<pair<int, int>>ret;
   for (auto e:s)
   ret.push_back(e);
   // 反時計回りに並べ替える
   mid = {0, 0};
   int n = ret.size();
   for (int i=0; i<n; i++){
      mid.first += ret[i].first;
      mid.second += ret[i].second;
      ret[i].first *= n;
      ret[i].second *= n;
   }
   sort(ret.begin(), ret.end(), compare);
   for (int i=0; i<n; i++)
   ret[i] = make_pair(ret[i].first/n, ret[i].second/n);
   return ret;
}
// 分割統治により凸包の値を返す
vector<pair<int, int>> divide(vector<pair<int, int>> a){
   // 点数が5以下なら総当たりで計算
   if (a.size() <= 5)
   return bruteHull(a);
   // left: 左半分の点 / right: 右半分の点
   vector<pair<int, int>>left, right;
   for (int i=0; i<a.size()/2; i++)
   left.push_back(a[i]);
   for (int i=a.size()/2; i<a.size(); i++)
   right.push_back(a[i]);
   // 左右それぞれの凸包を再帰的に求める
   vector<pair<int, int>>left_hull = divide(left);
   vector<pair<int, int>>right_hull = divide(right);
   // 2つの凸包をマージする
   return merger(left_hull, right_hull);
}
int main(){
   vector<pair<int, int> > a;
   a.push_back(make_pair(0, 0));
   a.push_back(make_pair(1, -4));
   a.push_back(make_pair(-1, -5));
   a.push_back(make_pair(-5, -3));
   a.push_back(make_pair(-3, -1));
   a.push_back(make_pair(-1, -3));
   a.push_back(make_pair(-2, -2));
   a.push_back(make_pair(-1, -1));
   a.push_back(make_pair(-2, -1));
   a.push_back(make_pair(-1, 1));
   int n = a.size();
   sort(a.begin(), a.end());
   vector<pair<int, int> >ans = divide(a);
   cout << "Convex Hull:\n";
   for (auto e:ans)
   cout << e.first << " "<< e.second << endl;
   return 0;
}

実行結果

Convex Hull:
-5 -3
-1 -5
1 -4
0 0
-1 1

コードのポイント

• quad 関数:中心点を原点としたときの各点の象限(1〜4)を判定し、角度順のソートに利用します。
• calc_line 関数:外積の符号によって点が直線のどちら側にあるかを判定します。戻り値が 0 の場合は同一直線上にあることを意味します。
• merger 関数:左側の凸包の右端点と右側の凸包の左端点からスタートし、上下の共通接線が見つかるまでインデックスを移動させます。
• bruteHull 関数:すべての点ペアについて、残りの点が全て直線の同一側に存在する場合、その2点は凸包の頂点であると判断します。計算量は O(n³) ですが、点数が少ない場合のみ使用されるため問題ありません。

このように分割統治法を用いることで、大量の点データでも効率的に凸包を計算できます。Andrew のモノトーンチェーン法など他のアルゴリズムと比較しながら、用途に応じて使い分けるとよいでしょう。