第23章
凸包の内側
（やさしい版）
Pearls of Functional Algorithm Design（関数プログラミングによるアルゴリズム設計の真珠）
どんな問題？
平面や空間にいくつかの点が散らばっているとき、それらをぐるっと包み込む一番外側の輪郭を凸包（とつほう、convex hull）と呼びます。平面ならゴム紐で点を囲んだときの形、空間なら点をシュリンクフィルムで包んだときの形をイメージしてください。
この章のテーマは、d 次元（2次元でも3次元でもよい一般の次元）の点集合について、次の二つを扱うことです。
- 凸包の仕様（「点 p が凸包の中にあるか？」をどう定義するか）
- 凸包を計算する逐次アルゴリズム（点を1つずつ足していく素朴な方法）
この章の面白いポイント
- アルゴリズムの導出はしない。その代わり、QuickCheck というテストライブラリで正しさを確かめる。
- 実際にテストしてみたらバグが見つかった。読者が気づけるかどうかを試すため、そのバグはあえて残してある。
準備 ― 点と単体を定義する
点は「有理数の座標」で表す
幾何アルゴリズムは、浮動小数点の丸め誤差が入ると簡単に壊れます。そこでこの章では、整数の演算だけで済ませるために、座標が有理数である点だけを扱います。d 次元の有理座標を持つ空間を Q(d) と書きます。
点は d+1 個の整数のリスト [x0, x1, ..., xd]（ただし xd ≠ 0）で表します。実際の座標は最後の値で割った [x0/xd, ..., xd−1/xd] になります（分母を明示的に持つイメージ）。
dimension :: Point → Int
dimension ps = length ps − 1
Dart
// 点は分母を末尾に持つ整数リスト。次元は「要素数 − 1」。
typedef Point = List<int>;
int dimension(Point ps) => ps.length - 1;
単体（simplex）とは何か
d-単体とは、Q(d) 上の d+1 個の点と、その並び順から決まる「向き」（+1 か −1）を組にしたものです。次元によって呼び名が変わります。
- 1-単体：辺（2点）
- 2-単体：三角形（3点）
- 3-単体：四面体（4点）
type Simplex = ([Point], Int)
Dart
// 単体 = 点のリスト + 向き（+1 / −1）
typedef Simplex = (List<Point> us, int b);
typedef Facet = (List<Point> us, int b);
図 23.1 正の向きを持つ四面体の例。a, b, c が (x, y) 平面にあり、d がその上に浮かんでいる。
一般の位置（general position）
頂点が「一般の位置にある」とは、
都合の悪い一致がないことです。つまり、
- 同じ点がない
- 3点が同一直線上にない
- 4点が同一平面上にない、…と次元まで続く
これは、頂点の座標を並べた行列
X = (
xij) の
行列式が 0 でないことと同値です。
向き（orientation）を行列式の符号で決める
行列式の符号を、そのまま単体の「向き」とします。1
orientation :: [Point] → Int
orientation = fromIntegral · signum · det
Dart
// 行列式：ライプニッツ展開（順列に沿って符号付きで加算）。次元は小さい想定。
int det(List<List<int>> m) {
final n = m.length;
final idx = [for (var i = 0; i < n; i++) i];
var total = 0;
void perm(int k, int sign) {
if (k == n) {
var t = sign;
for (var i = 0; i < n; i++) t *= m[i][idx[i]];
total += t;
return;
}
for (var i = k; i < n; i++) {
final t = idx[k]; idx[k] = idx[i]; idx[i] = t;
perm(k + 1, (i == k) ? sign : -sign);
final u = idx[k]; idx[k] = idx[i]; idx[i] = u;
}
}
perm(0, 1);
return total;
}
// 単体の向き = 行列式の符号（+1, 0, −1）。
int orientation(List<Point> us) => det(us).sign;
向きの直観的な意味は次のとおりです。
- Q(1) の辺 [a, b]：a > b なら +1、a < b なら −1、a = b なら 0。
- Q(2) の三角形 [a, b, c]：反時計回りなら +1、時計回りなら −1、同一直線上なら 0。
- Q(3) の四面体 [a, b, c, d]：a から三角形 [b, c, d] を見たとき、それが正の向きなら +1。
頂点を偶数回入れ替えても向きは変わらず、奇数回入れ替えると符号が反転します。
単体の面（facet）と、単体の中にあるかの判定
単体 smp の 面（facet）とは、頂点を1つ抜いてできる「一段低い次元の単体」に、元の向きから決まる向きを添えたものです。
- 辺の面 → その両端の頂点（0-単体）
- 三角形の面 → その3本の辺
- 四面体の面 → その4つの三角形
facets :: Simplex → [Facet]
facets (us, b) = zip (minors us) (cycle [b, −b])
minors :: [a] → [[a]]
minors [] = []
minors (x : xs) = xs : map (x :) (minors xs)
Dart
// 1要素ずつ抜いた部分列。minors([a,b,c,d]) = [[b,c,d],[a,c,d],[a,b,d],[a,b,c]]
List<List<T>> minors<T>(List<T> xs) {
if (xs.isEmpty) return [];
final x = xs.first;
final rest = xs.sublist(1);
return [rest, ...minors(rest).map((m) => [x, ...m])];
}
// 各面の向きは +b, −b, +b, −b, ... と交互。
List<Facet> facets(Simplex smp) {
final (us, b) = smp;
final ms = minors(us);
return [for (var i = 0; i < ms.length; i++) (ms[i], i.isEven ? b : -b)];
}
ここで Facet = ([Point], Int)。minors は1要素ずつ抜いた部分列を並べる関数で、たとえば minors "abcd" = ["bcd", "acd", "abd", "abc"] となります。
各面には、元の向きに対して +1, −1, +1, −1, ... と交互の符号が付きます。たとえば正の向きの三角形 [a, b, c] の3つの面は、
([b, c], +1), ([a, c], −1), ([a, b], +1)
「中にあるか」の判定のアイデア
点 p が単体の中にあるかを判定するには、各面の向きと同じ側に p がいるかを、面ごとにチェックします。具体的には、面 (us, b) について、p を us の先頭に付けた新しい単体の向きが b と一致すれば、p はその面の「内側」にあるとみなします。
これを述語として書くと次のようになります。
insideCS :: Simplex → Point → Bool
insideCS smp p
= and [0 ≤ b ∗ orientation (p : us) | (us, b) ← facets smp]
Dart
// 点 p が単体 smp の内部か境界上にあれば true。
bool insideCS(Simplex smp, Point p) {
for (final (us, b) in facets(smp)) {
if (b * orientation([p, ...us]) < 0) return false;
}
return true;
}
すべての面について同符号（積が 0 以上）ならば、p は単体の内部か境界上にいます。積がちょうど 0 になる面があれば境界上、すべての積が正なら厳密に内部です。
凸包の仕様
「あらゆる単体の和集合」として定義する
Q(d) 上の点集合 vs の凸包 CH(vs) は、次のように定義できます。
凸包の定義
vs の点から作れるすべての d-単体 smp」について、CS(smp)（単体の中の点全体）を集めた和集合。
たとえば平面（Q(2)）なら、vs の3点で作れるすべての三角形の内部を寄せ集めたものが凸包です。これをそのまま述語にすると：
insideCH :: [Point] → Point → Bool
insideCH vs p = or [insideCS smp p | smp ← simplexes vs]
simplexes :: [Point] → [Simplex]
simplexes vs = [(us, b) | us ← tuples (d + 1) vs,
let b = orientation us, b ≠ 0]
where d = dimension (head vs)
Dart
// vs から n 個を選ぶすべての組合せ（順序保持）。
Iterable<List<T>> tuples<T>(int n, List<T> vs) sync* {
if (n == 0) { yield []; return; }
for (var i = 0; i <= vs.length - n; i++) {
for (final rest in tuples(n - 1, vs.sublist(i + 1))) {
yield [vs[i], ...rest];
}
}
}
// d+1 点で作れる、退化していないすべての単体。
List<Simplex> simplexes(List<Point> vs) {
final d = dimension(vs.first);
final result = <Simplex>[];
for (final us in tuples(d + 1, vs)) {
final b = orientation(us);
if (b != 0) result.add((us, b));
}
return result;
}
// 仕様版：p がいずれかの単体に入っていれば凸包内。
bool insideCH(List<Point> vs, Point p) {
for (final smp in simplexes(vs)) {
if (insideCS(smp, p)) return true;
}
return false;
}
tuples n vs は「vs から n 個選んだ部分列すべて」を返します。定義は演習として省略されています。
仕様の限界
この定義では、点がすべて同一直線上（Q(2)）や同一平面上（Q(3)）にあるなど、d-単体がひとつも作れない場合、凸包は空集合になります。次元を下げれば凸包を定義できますが、ここでは仕様をそのままにしておきます（あとでこの点がバグの引き金になります）。
逐次アルゴリズム
まずは効率を上げる方向を決める
n 点から d-単体は Ω(nd+1) 個も作れてしまうので、insideCH の素朴な実装は非常に遅くなります。多くの単体は重なっているので、重なりなく凸包全体を覆う「分割」が得られれば十分です。そこで simplexes を partition（分割を返す関数）に取り替えます。
insideCH′ :: [Point] → Point → Bool
insideCH′ vs p = or [insideCS smp p | smp ← partition vs]
Dart
// 分割版：重ならない単体のリストで凸包を覆う。
bool insideCHp(List<Point> vs, Point p) {
for (final smp in partition(vs)) {
if (insideCS(smp, p)) return true;
}
return false;
}
分割を作る戦略：1点ずつ足す
凸包の分割は、次の手順で作れます。
- 最初に、単体をひとつ見つけて出発点とする（findSimplex）。
- その頂点を除いた残りの点を1つずつ順に取り出し、現在の分割を更新していく（update）。
partition :: [Point] → [Simplex]
partition vs
= case findSimplex vs of
Nothing → [ ]
Just [smp] → foldl update [smp] (vs \\ vertices smp)
vertices :: Simplex → [Point]
vertices = sort · fst
Dart
// 点リストの辞書順比較。
int cmpPoint(Point a, Point b) {
for (var i = 0; i < a.length && i < b.length; i++) {
final c = a[i].compareTo(b[i]);
if (c != 0) return c;
}
return a.length.compareTo(b.length);
}
// 単体の頂点集合（ソート済み）。
List<Point> vertices(Simplex smp) {
final vs = [...smp.$1];
vs.sort(cmpPoint);
return vs;
}
// ソート済みリストの差集合 xs \ ys。
List<Point> diff(List<Point> xs, List<Point> ys) {
return [for (final x in xs) if (!ys.any((y) => cmpPoint(x, y) == 0)) x];
}
List<Simplex> partition(List<Point> vs) {
final smp = findSimplex(vs);
if (smp == null) return [];
var acc = [smp];
for (final v in diff(vs, vertices(smp))) {
acc = update(acc, v);
}
return acc;
}
vs をソート済みにしておけば、\\（差集合）は順序付きリストの差として効率よく行えます。あとは findSimplex と update を作れば完成です。
単体を見つける ― findSimplex
素朴な定義とその弱点
まず一番素直な定義は次のようになります。
findSimplex vs = if null smps then Nothing else Just (head smps)
where smps = simplexes vs
Dart
// 素朴版：すべての単体を列挙して先頭を返す。
Simplex? findSimplexNaive(List<Point> vs) {
final smps = simplexes(vs);
return smps.isEmpty ? null : smps.first;
}
ところが、最悪の場合 Ω(nd+1) ステップかかってしまい、せっかく効率化しようとした意味が薄れます。実用上は問題ないことが多いですが、別のやり方も紹介されています。
「点を1つずつ足しながら一般の位置になるか調べる」方法
方針はこうです。
- 先頭の v0 から始める。
- 次に v0 と重ならない v1 を探す。
- その次に v0, v1 と同一直線上にない v2 を探す。
- これを繰り返し、一般の位置にある d+1 点が揃うまで続ける。
難しいところ
Q(d) 上の k+1 点（k < d）が一般の位置にあるかを判定するとき、点の座標行列は (k+1) × (d+1) の正方でない行列になるので、単純に行列式では判定できません。
正方部分行列を全部見る「退化判定」
そこで、最初の d 列から k 列を選ぶあらゆる組合せに、最終列（分母の列）を合わせて、(k+1) × (k+1) の正方部分行列を作ります。これらの行列式がすべて 0 なら退化している、と判定します。
degenerate k = all (== 0) · map det · submatrices k · transpose
submatrices k vs = map (++ [last vs]) (tuples k (init vs))
Dart
// 行列の転置。
List<List<int>> transpose(List<List<int>> m) {
if (m.isEmpty) return [];
return [for (var j = 0; j < m.first.length; j++)
[for (var i = 0; i < m.length; i++) m[i][j]]];
}
// 最初の要素から k 個を選び、最終列（分母列）を必ず付ける。
List<List<List<int>>> submatrices(int k, List<List<int>> vs) {
final init = vs.sublist(0, vs.length - 1);
final last = vs.last;
return [for (final t in tuples(k, init)) [...t, last]];
}
// k+1 点（それを転置した行列）が退化しているか。
bool degenerate(int k, List<Point> us) {
for (final sm in submatrices(k, transpose(us))) {
if (det(sm) != 0) return false;
}
return true;
}
転置してから列ではなく行を選ぶ形にしていますが、行列式は転置しても変わらないので同じことです。部分行列は O(dk) 個、行列式の計算に O(k3) かかるので、退化判定は k ≤ d のとき O(dd+3) ステップで行えます。
findSimplex :: [Point] → Maybe Simplex
findSimplex [] = Nothing
findSimplex (v : vs) = search (length v − 1) 1 [v] vs
search d k us vs
| k == d + 1 = Just (us, orientation us)
| null vs = Nothing
| degenerate k (v : us) = search d k us (tail vs)
| otherwise = search d (k + 1) (v : us) (tail vs)
where v = head vs
Dart
// 1点ずつ足しながら、一般の位置にある d+1 点が揃うまで探す。
Simplex? findSimplex(List<Point> vs) {
if (vs.isEmpty) return null;
final v0 = vs.first;
return _search(v0.length - 1, 1, [v0], vs.sublist(1));
}
Simplex? _search(int d, int k, List<Point> us, List<Point> vs) {
while (true) {
if (k == d + 1) return (us, orientation(us));
if (vs.isEmpty) return null;
final v = vs.first;
if (degenerate(k, [v, ...us])) {
vs = vs.sublist(1); // v を使わず次へ
} else {
us = [v, ...us]; // v を採用して次元を上げる
vs = vs.sublist(1);
k++;
}
}
}
計算量のまとめ
n 点に対する findSimplex は O(dd+3n) ステップ。n に関して線形だが、次元 d に関して非常に急速に増える定数を持つ。
更新 ― update
面には「内部」と「外部」の2種類がある
いま扱っている点に対する分割 smps を考えます。すべての単体の面を並べると、面は次の2種類に分けられます。
- 内部面：ちょうど 2 回、互いに反対の向きで現れる。
- 外部面：ちょうど 1 回だけ現れる。
例：正方形 [a, b, c, d] の三角形分割
- [a, b, c] と [c, d, a] に分けると、辺 [a, c] が内部。
- [a, b, d] と [b, c, d] に分けると、辺 [b, d] が内部。
それ以外の辺（外周の4辺）は外部面です。
外部面を取り出すには、すべての面を並べてソートし、隣り合う同じ頂点集合の面を打ち消し合わせます。
external :: [Simplex] → [Facet]
external = foldr op [ ] · sort · concatMap facets
op smp [ ] = [ ]
op smp (smp′ : smps)
= if vertices smp == vertices smp′ then smps
else smp : smp′ : smps
Dart
// 面の頂点集合（ソート済み）で比較。
int cmpFacetVerts(Facet a, Facet b) {
final av = [...a.$1]..sort(cmpPoint);
final bv = [...b.$1]..sort(cmpPoint);
for (var i = 0; i < av.length && i < bv.length; i++) {
final c = cmpPoint(av[i], bv[i]);
if (c != 0) return c;
}
return av.length.compareTo(bv.length);
}
// すべての面を並べ、同じ頂点集合の面（内部面）を打ち消し、外部面だけ残す。
List<Facet> external(List<Simplex> smps) {
final all = <Facet>[for (final s in smps) ...facets(s)];
all.sort(cmpFacetVerts);
final result = <Facet>[];
for (final f in all) {
if (result.isNotEmpty && cmpFacetVerts(result.last, f) == 0) {
result.removeLast(); // ペアで消える → 内部面
} else {
result.add(f);
}
}
return result;
}
コストはソートが支配的で O(dS log dS) です（S は単体の個数）。
電球のたとえ ― 可視な面を見つける
新しい点 v が凸包の外側に落ちたとき、外部面は「v から見える面（可視）」と「見えない面（不可視）」に分かれます。
電球のたとえ
点 v の位置に電球があると想像しましょう。この電球が照らせるのが可視な面、影になるのが不可視な面です。式で言うと、面 (us, b) について orientation(v : us) の符号が b と反対のとき可視です。
visible :: Point → [Facet] → [Facet]
visible v fs = [(us, b) | (us, b) ← fs, b ∗ orientation (v : us) < 0]
Dart
// v から「見える」面：v を面の前に付けた向きが元の向きと反対のもの。
List<Facet> visible(Point v, List<Facet> fs) {
return [for (final (us, b) in fs) if (b * orientation([v, ...us]) < 0) (us, b)];
}
v が現在の凸包の中や境界にあるときは可視な面がありません。ですから、入力に重複点があっても構いません（凸包は変わらない）。
可視な面ごとに新しい単体を追加
更新は「可視な面それぞれに、新しい点 v を頂点として付け加えた単体を作る」だけです。
newSimplex :: Point → Facet → Simplex
newSimplex v (us, b) = (v : us, −b)
update :: [Simplex] → Point → [Simplex]
update smps v
= smps ++ map (newSimplex v) (visible v (external smps))
Dart
// 可視な面それぞれに v を頂点として付けて新しい単体を作る。向きは反転。
Simplex newSimplex(Point v, Facet f) {
final (us, b) = f;
return ([v, ...us], -b);
}
List<Simplex> update(List<Simplex> smps, Point v) {
final vis = visible(v, external(smps));
return [...smps, for (final f in vis) newSimplex(v, f)];
}
付ける向きが −b になる理由は簡単です。可視の定義から orientation(v : us) = −b となっているからです。
計算量
update のコストは可視面の計算が支配的で O(dS log dS)。全体としては O(dnS log dS)。単体の最大個数 S = O(ne)（e = ⌊d/2⌋）が知られているので、insideCH′ は O(ne+1 log n) となります。insideCH よりずっと速いですが、まだ改善できます。
さらなる改良 ― 単体ではなく「面」を保持する
アイデア：単体は捨てて、外部面だけを持ち続ける
アルゴリズムを見直すと、毎回外部面だけを使って可視性を判定していることに気づきます。だったら、最初から単体ではなく外部面だけを保持すればいいはず。faces = external · partition と置けば、次の insideCH″ ができます。
insideCH″ vs p = and [0 ≤ b ∗ orientation (p : us) | (us, b) ← faces vs]
Dart
// 面リスト版：外部面だけで凸包を判定する（注意：バグあり版）。
bool insideCHppBuggy(List<Point> vs, Point p) {
for (final (us, b) in faces(vs)) {
if (b * orientation([p, ...us]) < 0) return false;
}
return true; // 面が空でも true になってしまう！
}
これは計算幾何学で凸包を「外部面のリスト」として表す一般的なやり方に対応しています。
foldl の融合則で導出
faces の効率的な計算は、foldl の融合則から次を満たす update′ を求めれば得られます。
external (update smps v) = update′ (external smps) v
faces vs
= case findSimplex vs of
Nothing → [ ]
Just [smp] → foldl update′ (facets smp) (vs \\ vertices vs)
update′ fs v = (fs \\ fs′) ++ map (newFacet v) (external fs′)
where fs′ = visible v fs
newFacet v (us, b) = (v : us, b)
Dart
// 新しい面は「可視面の外周 + 新点 v」で張る。向きはそのまま。
Facet newFacet(Point v, Facet f) {
final (us, b) = f;
return ([v, ...us], b);
}
// external を Facet に対して適用（内部面のペアを打ち消す）。
List<Facet> externalF(List<Facet> fs) {
final all = [...fs]..sort(cmpFacetVerts);
final result = <Facet>[];
for (final f in all) {
if (result.isNotEmpty && cmpFacetVerts(result.last, f) == 0) {
result.removeLast();
} else {
result.add(f);
}
}
return result;
}
List<Facet> updateP(List<Facet> fs, Point v) {
final vis = visible(v, fs);
final kept = [
for (final f in fs)
if (!vis.any((g) => cmpFacetVerts(f, g) == 0)) f
];
return [...kept, for (final f in externalF(vis)) newFacet(v, f)];
}
List<Facet> faces(List<Point> vs) {
final smp = findSimplex(vs);
if (smp == null) return [];
var fs = facets(smp);
for (final v in diff(vs, vertices(smp))) {
fs = updateP(fs, v);
}
return fs;
}
言葉で言うと
- 新しい点から見える面を全部消す。
- そのすき間を埋めるように、可視な面の集合の外周（(d−2)-単体の集合）を境に、新しい点を頂点として面を張り直す。
- 3次元なら「見える三角形の集合をシールみたいに剥がして、その外周の辺と新しい点で新しい三角形を作る」イメージ。
面全体を毎回スキャンして可視面を探しているので、まだ非効率です。可視面は局所的にまとまっているので、隣接情報を使えばもっと速くできます。Blelloch et al. (2001) の Bulldozer アルゴリズムなどがその方向の代表例ですが、本章では立ち入りません。
QuickCheck で試す ― そしてバグ発見
ランダムに点を作って比較する
Claessen と Hughes による QuickCheck を使うと、ランダムなテストケースで性質を確かめられます。まず、点とその集合を生成するジェネレータを用意します。
point :: Int → Gen [Integer]
point d = do {xs ← vector d; return (xs ++ [1])}
points :: Int → Int → Gen [[Integer]]
points d 0 = return [ ]
points d (n + 1) = do {p ← point d; ps ← points d n;
return (p : ps)}
Dart
import 'dart:math';
// d 次元のランダムな点（分母は 1 で固定）。座標範囲は −range..range。
Point randomPoint(Random rng, int d, {int range = 5}) {
return [for (var i = 0; i < d; i++) rng.nextInt(2 * range + 1) - range, 1];
}
// n 個の点。
List<Point> randomPoints(Random rng, int d, int n, {int range = 5}) {
return [for (var i = 0; i < n; i++) randomPoint(rng, d, range: range)];
}
そして、仕様 insideCH と実装 insideCH′ が一致することを性質として書きます。
prop_Hull :: Int → Int → Property
prop_Hull d n = forAll (points d n) $ λvs →
forAll (point d) $ λv →
insideCH vs v == insideCH′ vs v
OK, passed 100 tests.
Dart
// 仕様 insideCH と実装 insideCHp が一致するかをランダムに検証。
bool propHull(int d, int n, {int trials = 100, int seed = 42}) {
final rng = Random(seed);
for (var t = 0; t < trials; t++) {
final vs = randomPoints(rng, d, n);
final v = randomPoint(rng, d);
if (insideCH(vs, v) != insideCHp(vs, v)) return false;
}
return true;
}
insideCH″ に置き換えると…失敗！
ところが、prop_Hull の中で insideCH′ を insideCH″ に置き換えるとテストが落ちます。
Main> quickCheck (prop_Hull 2 4)
Falsifiable, after 2 tests:
[[0,0,1],[0,0,1],[−1,−1,1]]
[1,0,1]
バグの原因
反例の4点はすべて同一直線上にあります。したがって分割も面も作れず、faces vs は空リスト。
- insideCH vs は正しく False を返す。
- ところが insideCH″ vs は、空リストに対する
and が定義上 True を返してしまう。
つまり、
「面がひとつもない」場合を特別扱いする必要があるのに、それを忘れていたわけです。
修正
insideCH″ vs v = if null fs then False else
and [0 ≤ b ∗ orientation (v : us) | (us, b) ← fs]
where fs = faces vs
Dart
// 修正版：面がひとつも作れなければ False を返す。
bool insideCHpp(List<Point> vs, Point v) {
final fs = faces(vs);
if (fs.isEmpty) return false;
for (final (us, b) in fs) {
if (b * orientation([v, ...us]) < 0) return false;
}
return true;
}
Dart
void main() {
// 平面 (d=2) の正方形の4頂点で凸包を作り、点の内外を判定する。
final square = <Point>[
[0, 0, 1],
[2, 0, 1],
[2, 2, 1],
[0, 2, 1],
];
final inside = [1, 1, 1]; // 正方形の中央
final outside = [3, 1, 1]; // 右外
print(insideCH(square, inside)); // true
print(insideCH(square, outside)); // false
print(insideCHp(square, inside)); // true
print(insideCHpp(square, outside)); // false
// 反例：4点すべて同一直線上 → faces は空 → 修正版は false を返す。
final degen = <Point>[
[0, 0, 1], [0, 0, 1], [-1, -1, 1],
];
print(insideCH(degen, [1, 0, 1])); // false
print(insideCHpp(degen, [1, 0, 1])); // false（修正版で正しく）
}
面が空なら False を返すように直せば QuickCheck も通ります。この間違い、あなたは気づけたでしょうか？
おわりに
凸包アルゴリズムを扱う計算幾何学の教科書は数多くあり、O'Rourke (1998) や Preparata and Shamos (1985) はそのうちの2冊です。特に O'Rourke の本には凸包に丁寧な2章が割かれ、参考文献も充実しています。ただし Dijkstra (1976) の3次元の扱いは載っていません。この真珠は Karimipour and Frank (2009) の詳細に取り組む過程で生まれましたが、細部はかなり異なります。凸包について有益な議論を交わしてくれた Irina Voiculescu 氏に感謝します。
参考文献
Claessen, K. and Hughes, J. (2000). QuickCheck: a lightweight tool for random testing of Haskell programs. ACM SIGPLAN International Conference of Functional Programming, Montreal, Canada, pp. 268–79. http://www.cs.chalmers.se/~rjmh/QuickCheck/ も参照。
Blelloch, G., Burch, H., Crary, K., et al. (2001). Persistent triangulations. Journal of Functional Programming 11 (5), 441–66.
Dijkstra, E. W. (1976). A Discipline of Programming. Englewood Cliffs, NJ: Prentice-Hall.
Karimipour, F. and Frank, A. U. (2009). A dimension independent convex hull algorithm. Unpublished.
O'Rourke, J. (1998). Computational Geometry, second edition. Cambridge, UK: Cambridge University Press.
Preparata, F. P. and Shamos, M. I. (1985). Computational Geometry. New York, NY: Springer-Verlag.