どちらもアルゴリズム的には, シンプルですが, Quickhullの方は, 理解するのに少し時間がかかりました. Jarvis's Marchは, ギフト包装法という呼び方の方が一般的なようです. どちらも基本的には幾何学の問題で, 高校数学の基礎的な内容と外積を駆使します.
画像としては以下のような感じ. 点の集合を囲む凹みの無いつつみが凸包です. 凸包は, 一部の線形計画法の問題などで使われますが, 抽象解釈などにも凸多面体などとして, 応用されているという点が面白いですね.
以下では, 2次元について考えます.
Jarvis's Marchから.
求めるのは, 凸包を構成する点列です. そしてその点列に含まれるのが自明な点は, 最左, 最右と, 最上, 最下の点でしょう. Jarvis's Marchのアルゴリズムは, この４つのどの点でもいいのですが, とりあえずとっかかりとして, 最左の点(x座標が一番小さい点)から始めます.
1. 現在の点から他の全ての点へのベクトルの角度を求め, 角度が最小の点を凸包の点列に追加し, これを次の点とする.
2. 始めた点へ戻ってくるまで1を続ける.
というのがアルゴリズムで, 非常に簡単でした. 1で角度が同じ場合は, 現在の点から距離が近い方を選びます. 角度を求める手法は何でも良いらしいですが, Wikipediaに従い, 外積を用いることにしました. (外積の意味については, 参考 : 基礎の基礎その1 内積と外積の使い方 - ○×つくろーどっとコム)
Wikipediaの擬似コードを無理やり関数型言語へ押し込むと以下のようになりました. 多分, 手続き型言語よりも整理出来ていると思います. (incanterを使っています.)
(require '[incanter.core :as in]
'[incanter.stats :as st]
'[incanter.charts :as ch])
;; -- Jarvis's March to find convex hull
(defn left-most-point [[x1 y1 :as p1] [x2 y2 :as p2]]
(if (< x1 x2) p1 p2))
(defn ex-product
"exterior product"
[[x1 y1 :as p1] [x2 y2 :as p2]]
(- (* x1 y2) (* x2 y1)))
(defn -vec [[x1 y1 :as p1] [x2 y2 :as p2]]
[(- x1 x2) (- y2 y1)])
(defn min-angle [current-p point-b point-c]
(let [v (ex-product (-vec current-p point-b) (-vec current-p point-c))]
(if-not (or (= point-b current-p)
(< 0 v)
(and (< -0.001 v 0.001)
(< (st/euclidean-distance current-p point-b)
(st/euclidean-distance current-p point-c))))
point-b point-c)))
(defn jarvis-march
"find convex hull, points = [[x1, y1], ...]"
[points]
(let [start-p (reduce left-most-point points)]
(loop [current-p start-p hull-ps [start-p]]
(let [next-p (reduce (partial min-angle current-p) points)]
(if (and (= next-p start-p))
hull-ps
(recur next-p (cons next-p hull-ps)))))))
;; --- visualization --- (use incanter)
(defn plot-convex-hull [points hull-ps]
(let [hull-ps (cons (last hull-ps) hull-ps)]
(ch/add-lines
(ch/scatter-plot (map first points) (map second points))
(map first hull-ps)
(map second hull-ps)
:auto-sort false)))
;; --- test ---
(defn generate-sample-set []
(letfn [(generate-by [m]
(st/sample-normal
(+ 100 (rand-int 100)) :mean m :sd (+ 0.15 (rand 0.1))))]
(map list
(generate-by (dec (rand 2)))
(generate-by (dec (rand 2))))))
(defn trial []
(let [ps (generate-sample-set)
result (jarvis-march ps)]
(in/view (plot-convex-hull ps result))))
このアルゴリズムは, 計算量が大きそうですが, 自明な凸包の輪郭とならない点を削ることで, 計算量を減らすことは可能なようです. 例えば, 前述の上下左右4点の内側の点は全て凸包の輪郭には含まれないのでこれらの点は無視できますね.意外に気づかずにハマったのは, incanterのchart/add-linesが引数として渡されたxy各座標のリストを勝手にソートしてくれることです. :auto-sortをfalseに設定しています.
次に, Quickhullのアルゴリズム. これは, Quicksort風の分割統治法によるアルゴリズムで, やはり, このアルゴリズムも自明な点から始めます.
最左と最右の点を初期値にします. 上記のWikipediaの記事とは書き方が少し異なりますが, 内容としては同じ(はず)です.
1. 最左/最右の点を求める.
2. 前述の二点(p-left-most, p-right-most)から線分を求め, これで点の集合を二分し, 上半分についてのみ考え, 3.に上記二点(p-left-most, p-right-most)と上半分の集合を渡す.
3. (関数. 引数は 二点(p1, p2)と点の集合)
3-1. 線分から最も遠い点fpを求める
3-2. これを凸包の点の一つとし, 線分とfpからなる三角形の内側の点を除外.
3-3. 三角形の外側に残る点のうち, fpより左側の点の集合と, fpより右側の点の集合に分ける.
3-4. 3.へp1, fpと左半分を渡し, p1〜pf間の凸包の部分列(輪郭の点列)を求める.
3-5. 3.へfp, p2と右半分を渡し, pf〜p2間の凸包の部分列を求める.
3-6. 3-4, pf, 3-5の結果を合わせて, p1〜p2間の凸包の部分列として返す.
4. 下半分も3.へ渡す.
5. 上半分, 下半分について求めた結果と, 最左最右の点を纏める => 凸包を構成する点列.
Jarvis's Marchと比較して, 少し手間がかかりますが, こんなイメージ.
全体として分割統治法で, 再帰的に求めていくと, 答えが求まるという流れになります. 線分からの距離は, y=ax+b(直線)からの距離を求めればOK. 三角形に含まれるかどうかも, 外積から求まります.
(require '[incanter.core :as in]
'[incanter.stats :as st]
'[incanter.charts :as ch])
(defn get-ab
"calc a,b for y = ax + b from p1, p2"
[[x1 y1] [x2 y2]]
(let [a (/ (- y1 y2) (- x1 x2))]
[a (- y2 (* a x2))]))
(defn ex-product
"exterior product"
[[x1 y1 :as p1] [x2 y2 :as p2]]
(- (* x1 y2) (* x2 y1)))
(defn -vec
"p1 - p2 (p1 & p2 is vector)"
[[x1 y1 :as p1] [x2 y2 :as p2]]
[(- x1 x2) (- y2 y1)])
(defn in-triangle?
[[x1 y1 :as p1] [x2 y2 :as p2] [x3 y3 :as p3] [tx ty :as tp]]
(let [from-a (ex-product (-vec p3 p1) (-vec p1 tp))
from-b (ex-product (-vec p1 p2) (-vec p2 tp))
from-c (ex-product (-vec p2 p3) (-vec p3 tp))
abc [from-a from-b from-c]]
(or (reduce #(and %1 %2) (map pos? abc))
(reduce #(and %1 %2) (map neg? abc)))))
(defn furthest-point [[x1 y1 :as xy1] [x2 y2 :as xy2] points]
(let [[a b] (get-ab xy1 xy2)
sqa2b2 (. Math sqrt (+ (. Math pow a 2) (. Math pow b 2)))]
(letfn [(get-dist [[x0 y0]] (/ (. Math abs (+ (* a x0) (* -1 y0) b)) sqa2b2))]
(reduce #(if (< (get-dist %1) (get-dist %2)) %2 %1) points))))
(defn find-points [p1 p2 points] ;; find points between p1 & p2
(if (empty? points)
[]
(let [[xf yf :as pf] (furthest-point p1 p2 points) ;; furthest point
rest-points (remove #(or (in-triangle? p1 p2 pf %) (= pf %)) points)]
(concat (find-points p1 pf (filter #(< (first %) xf) rest-points))
[pf]
(find-points pf p2 (filter #(<= xf (first %)) rest-points))))))
(defn quick-hull [points]
(let [p-left-most (reduce #(if (< (first %1) (first %2)) %1 %2) points)
p-right-most (reduce #(if (< (first %1) (first %2)) %2 %1) points)
[a b] (get-ab p-left-most p-right-most)]
(letfn [(upper? [[x y]] (< (+ (* a x) b) y))]
(concat [p-left-most]
(find-points p-left-most p-right-most (filter upper? points))
[p-right-most]
(->> (filter #(not (upper? %)) points)
(find-points p-left-most p-right-most)
reverse)))))
;; --- visualization --- (use incanter)
(defn plot-convex-hull [points hull-ps]
(let [hull-ps (cons (last hull-ps) hull-ps)]
(ch/add-lines
(ch/scatter-plot (map first points) (map second points))
(map first hull-ps)
(map second hull-ps)
:auto-sort false)))
;; --- test ---
(defn generate-sample-set []
(letfn [(generate-by [m]
(st/sample-normal
(+ 1000 (rand-int 500)) :mean m :sd (+ 0.15 (rand 0.1))))]
(map list
(generate-by (dec (rand 2)))
(generate-by (dec (rand 2))))))
(defn trial []
(let [ps (generate-sample-set)
result (quick-hull ps)]
(in/view (plot-convex-hull ps result))))
filterやconcatの使い方が, Haskell系の偽Quicksortにそっくりですね.参考文献
Quickhullアルゴリズム Project ASURA