Quickhullアルゴリズム

　あるものを作りたいと思ったのですが，まずは凸包を作成する必要があったので，
今回は凸包を作成するのに使われるQuickhullアルゴリズムを実装してみました。

　QuickhullアルゴリズムはQuickという名の通り，クイックソートの考え方に似ています。
今回は2次元という風に限定して説明しますが，3次元の場合は同様にしてできるそうです。
まず下の図のように点の集合があると仮定します。


最初に，これらの点の集合の中からx座標が最大と最小となる2点を求めて，稜線を引きます。
次に，線が引かれることによって，線の上側と下側の2つの領域にわけることができます。下の図の場合はうまく2つの領域に分かれますが，場合によっては線の上側に点がない場合もあり，2つの領域にわかれない場合もあります。


つぎに，各領域に関して着目し，各点から稜線に向かって垂線を下ろし，垂線の長さが最大となる点を求めます。


垂直距離が最大となる点が見つからなかった場合はそこで処理を終了し，見つかった場合は新たに稜線を追加します。
このとき稜線を追加することによって三角形が形成されます。形成される三角形の内部に点がある場合は，凸包を求める処理からはずします。下の図でいうとオレンジ色の点が処理の対象から外れる点になります。


稜線を追加することによってまた線より上側と下側の領域にわけることができます。あとは，先ほどと同じように各点から稜線に向かって垂線を下ろし，垂線の長さが最大となる点を見つけます。見つからなかった場合はその領域での処理を終了します。


あとは，同じように処理を繰り返し行います。処理できる点がなくなれば，凸包の完成です。


quickhullアルゴリズムの基本的な考え方以上のようになります。

さて，プログラムの実装ですが，そんなに難しくはありません。
まず，最初の稜線を求める処理ですが，単純にx座標をif文を使って比較するだけなので問題ないでしょう。
つぎに稜線の上側か下側かを判定する処理ですが，うまい方法がすぐに浮かばなかったので，直線の方程式（1次関数）に，稜線を形成する2点と判定したい1点の計3点の座標を代入して，不等式を解くやり方にしました。これも高校まで行っている人ならばわかると思うので問題ないでしょう。
お次は，垂線の長さを求める処理ですが，「点と直線の距離の公式」とかいう公式もあるのですが，それを使わないで直接垂線の足の座標を求めて，2点間の長さを求めるやり方にしました。で，求め方なんですけど，下の図のように座標を適当に決めてやります。


上の図で，法線ベクトルはベクトルABを反時計回りに90度回転した点となるのですぐに求めることができます。
求めたい垂線の足Hは点Aを始点にして方向ベクトルABをt倍した点として表すことができます。またHは点Cを視点にして法線ベクトルをs倍した点としても表すことができます。
点Hを今述べたように2通りの考え方で式に表すことができます。


上の2つの式からx, yを代入法によって消去すると右上のような式になります。この式は未知数2つに式2つなので，連立方程式によって解くことができます。
まず，この式をsについて解いてみます。sについて解くと下のようになります。


次にsを再び式に代入し，tについて解きます。tについて解くと下のようになります。


このようして求まったtを2つ上の式に代入してやれば，垂線の足Hの座標値が求まります。Hの座標が求まったら，あとはいつも通りsqrt()関数を使って距離を求めるだけです。
最後に，三角形の内部にあるかどうかという判定処理ですが，三角形の面法線と三角形を形成する3点a, b, cと判定したい点dの外積( (b-a)×(d-a)，(c-b)×(d-b)，(a-c)×(d-c) )の内積を取って，3つの符号がすべて同じならば三角形内部，そうでないなら外部という風に判定するやり方を使ってみました。
imagireさんのページの「レイトレース：3角形と視線の交差判定」に説明が載っていますのでそちらを参照してください。
あとは，何回も繰り返して同じ処理をするので，再帰を使うと簡単に実装できるそうです。なので，再帰関数を使って実装してみました。
あんまし，きれいにコードは書けていませんが，その辺はご愛敬ということで…ｗ。

本ソースコードおよびプログラムを使用したことによる如何なる損害も製作者は責任を負いません。
本ソースコードおよびプログラムは自己責任でご使用ください。
プログラムの作成にはMicrosoft Visual Studio 2008 Professionalを用いています。

• QuickHull.cab (64.4KB)