凸多角形の最適三角形分割

はいどうも！ この記事は、数学とコンピュータ・アドベントカレンダーの11日目の記事です。 凸多角形の最適三角形分割の話をします。

多角形の三角形分割 - Wikipedia

凸多角形の定義とは、文献１によれば、「多角形の内部のどの2点を結んでも、その2点を両端点とする線分がその多角形の内部に含まれるとき、その多角形は凸である」ということだそうです。「多角形Pが凸多角形であるための必要十分条件は、Pの全ての対角線がPの内部にあるときである」と言い換えられます。

こんな感じに三角分割する

凸多角形を分割する

凸多角形(単調な多角形)を任意の三角形で分割するだけなら、特別なアルゴリズムを考える必要もなさそうです。凸多角形の定義から、多角形の点を一つ選び、その線から多角形の他の点に対して線を引いていくだけでも凸多角形の三角形分割は可能です。というわけで、単に凸多角形を三角形に分解するだけなら簡単なので、凸多角形を三角形分割した時に、分割に使った対角線の総和が最小となるような分割を求めてみます。

この問題を分割統治法で解きます(後述しますが、この方法は、動的計画法によって解くのが一般的です)。

凸多角形を対角線で分割するため、解は必ず多角形の頂点のペアの列として表せます。つまり、頂点のペアの列を求めるわけですが、より正確には、合計の長さが最小となる対角線を表す頂点のインデックスのリストつまり、i番目とj番目の頂点のペアのリストが求めたい解となります。

凸多角形がどんな形で分割されたとしても、必ず、四角形以上であれば、必ず対角線が一本引かれます。つまり、$n$個の頂点を持つ多角形の対角線のうち、最初の1つの対角線を一本、どこかに引く必要があります。
そのような直線が一本決まると、残った部分の多角形の三角分割は、元の多角形の部分多角形の三角分割の問題となります。(三角分割なので、当然、対角線どうしが交差することはない為、決まった対角線によって分割された2部分多角形の三角分割の問題になります。)

まず始めに、凸多角形の全体を考えるのではなく、凸多角形の部分多角形の分割について考えてみたいと思います。話を簡単にするために、頂点を$v_i$(凸多角形の$i$番目の頂点)の$i$番目〜$j$番目間の分割を表す関数を$L[i, j]$とします。$L[i, j] (|i - j| < 2)$となる場合については、分割出来ないため、分割なし、分割を表す対角線の長さは0とします。また、インデックスですが、頂点$n$の多角形に対し$i+1$や$j-1$などは$\pmod{n}$で循環するものとします。

こんな風に分割

部分多角形の分割を表す$L[i, j]$に対して、部分多角形内の頂点$v_i, v_j$を含むような三角形を考えます。今求めようとしている部分多角形の性質上、頂点$v_i$と$v_j$間が別の対角線によって遮られる事はあり得ません。したがって、$v_i$, $v_j$をともうひとつの頂点から成る三角形によって分割される筈です。そのような三角形の頂点は、$v_i+1$cdots v_j-1$の中から必ず1つ選ばれる事になります。(選ばれなかった場合、三角分割として成り立たないため)
このような三角形の選び方は3通りあります。

1. $(v_i, v_j, v_{j-1})$で三角形が構成されるパターン
2. $(v_i, v_{i+1}, v_j)$で三角形が構成されるパターン
3. $(v_i, v_k, v_j)$で三角形が構成されるパターン(ただし、$i+1 < k < j-1$)

3通りの分割

1, 2, 3は3の形式でひとまとめに出来そうですが、1, 2は、引かれる対角線が1つであるのに対して、3では2つ対角線が引かれるという点が異なります。そして、それぞれ、以下のように3通りの部分多角形の三角分割を行った結果と足し合わせる形で表せます。

1. $L[i, j] = (i, j-1)$ 及び $L[i, j-1]$
2. $L[i, j] = (i+1, j)$ 及び $L[i+1, j]$
3. $L[i, j] = (i, k), (k, j)$ 及び $L[i, k]$ と $L[k, j]$

上記の条件をまとめると、三角分割した部分多角形$L[i, j]$の対角線の長さ$L'[i,j]$は、

$L'[i,j] = 0 (|i - j| < 2)$
$L'[i,j] = min \bigl\{$
    $L'[i + 1, j] + w[i + 1, j],$
    $L'[i, j - 1] + w[i, j - 1],$
    $\displaystyle \min_{k = (i+2) \dots (j-2)} \{ L'[i, k] + L'[k, j] + w[k, j] + w[k, j] \bigr\} \}$

であると言えます。$w[i, j]$は、頂点$v_i$, $v_j$間の重さ、$min { 〜 }$は、与えられた集合の中で最も小さい物を還す関数。添字の$k$も頂点を表すインデックスです。これにより、再帰的な関数によって線分$L'[i, j]$の長さが定義出来ることになりました。

ここでもう一度、元の問題、$n$個の頂点を持つ凸多角形の線分の長さが最小となる三角分割ですが、上記の部分多角形の三角分割$L[i, j]$の$L[0, n-1]$のケースであると見なせます。$L[0, n-1]$を解くことで、$n$個の頂点を持つ凸多角形の三角分割が可能になるのです。
ちなみに上記の$L'$ですが、コードに落とす時は、重さだけでなく、$L[i, j]$の要素のリストについても考慮する必要があります。

実装

というわけでClojureのコードです。
(実行結果を可視化するのが一番簡単だったため、ClojureとIncanterを使いました。)

(def weight st/euclidean-distance)

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [points i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min points i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight (points i+1) (points j)))]
          [lines w] (find-min points i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight (points i) (points j-1)))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min points i k)
                      [lines2 w2] (find-min points k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2
                      (weight (points i) (points k))
                      (weight (points k) (points j)))]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (find-min points 0 (- (count points) 1)))

weightは、重さ(長さ)を表す関数(incanterの関数)、前述の$w[i, j]$に相当するものです、+mは$\pmod{n}$を計算する関数、find-minが$L[i, j]$、triangulateは、$L[0, n-1]$を呼び出すだけの関数です。

find-minは、頂点のリストを引数に受け取り、頂点列のインデックスと対角線の長さの合計のペア[[$(i, j) \cdots$] $'L[i, j]$]を返します。skip-iとskip-jがそれぞれパターン1, 2の箇所で、(map search-inter 〜)で書かれている所がパターン3になります。最後にmin-pathで最も小さい対角線を求めています。ちなみに、紛らわしいかもしれませんが、スペースが入っていないi+1とj-1は変数です。

実行結果

こんな感じになります。一番上の画像と同じですが。

おわりに

さて、出来上がった関数ですが、計算の重複が多く、見事にメモ化再帰の実装対象になっています。文献１によれば、動的計画法でメモ化しながら、分割統治法的に問題を解くようです。(多角形の頂点の数を増やしていくと、めっちゃ重くなる。。。)つまり、メモ化再帰をテストするためのサンプルに使用できるということであります。メモ化再帰の実装サンプルに、凸多角形の最適三角分割を使用してみてはいかがでしょうか。
それでは！

全ソースコード

(require '[incanter.core :as in]
         '[incanter.stats :as st]
         '[incanter.charts :as ch])

;; ---- generate polygon: 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)))))))

(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 generate-sample []
  (let [ps (generate-sample-set)
        result (jarvis-march ps)]
    (vec (map (fn [[x y]] [(int (* x 100)) (int (* y 100))]) result))))

;; --- triangulation -----------------------
(def weight st/euclidean-distance)

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [points i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min points i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight (points i+1) (points j)))]
          [lines w] (find-min points i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight (points i) (points j-1)))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min points i k)
                      [lines2 w2] (find-min points k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2
                      (weight (points i) (points k))
                      (weight (points k) (points j)))]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (find-min points 0 (- (count points) 1)))

;; --- plot --------------------
(defn plot-triangles [points lines]
  (let [ys (concat points [(first points) (last points)])
        poly (-> (ch/scatter-plot (map first points) (map second points))
                 (ch/add-lines (map first ys) (map second ys) :auto-sort false))
        graph (reduce #(ch/add-lines %1 (map first %2) (map second %2) :auto-sort false)
                      poly
                      (partition 2 lines))]
    (in/view graph)))

(defn point-selector [points indices]
  (apply concat (map (fn [[start end]] [(points start) (points end)]) indices)))

(defn show [points]
  (let [[indices score] (triangulate points)
        xs (point-selector points indices)]
    (plot-triangles points xs)
    points))

参考文献

1. 浅野 哲夫, 計算幾何―理論の基礎から実装まで (アルゴリズム・サイエンスシリーズ―数理技法編), 共立出版, 2007