2017/12/11

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

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

凸多角形の定義とは、文献1によれば、「多角形の内部のどの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は変数です。

実行結果

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

おわりに

さて、出来上がった関数ですが、計算の重複が多く、見事にメモ化再帰の実装対象になっています。文献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

2017/11/05

HylomorphismとMetamorphism

はいどうも! Recursion SchemeのHylomorphismとMetamorphismの話です。 今年のアドベントカレンダーに向けたブログ投稿リハビリ用の記事になります。 

hylomorphismは、catamorphismとanamorphismの合成、 metamorphismは、hylmorphismの双対、anamorphismとcatamorphismの合成に当たる射になります。

catamorphismやanamorphismは、(単方向連結)Listで言う所のfoldrとunfoldに当たるものです。

(以下、長いため、catamorphismをcataと書くなど、適宜morphismを省いて書きます。)

hylomorphismは、cataとanaの合成射であり、
hylo = cata . ana
という合成です。 

プログラムにおけるhyloは、foldr/unfoldの合成、つまりデータの生成処理(unfold)と、 データの消費処理(foldr)が融合するような関数になります。

hyloを直接再帰の形式で書くことにより、anaで生成されるデータをcataが直接消費することで、 anaが生成する(代数的データ型による)データを実際に生成することなくcata/anaの合成処理が記述することが可能になります。

このような合成は、プログラムの融合変換に応用できます。
実際の所、recursion schemeは、cata、ana、hyloをプログラムで直接手で記述するというよりは、コンパイル時の中間処理として登場するイメージに近いと思います。

例えば、プログラム中に登場するmap、filter、foldr、unfoldなどをcata、ana、hyloに展開し、展開後のコードの合成部分を見て、cata/ana融合が登場すると、直接再帰の形式のhyloに書き換えるといった最適化が可能になります。 (例えば、参考文献3では、mapやsum等、リストを扱う関数をhylo形式に書き換え、hylo同士を融合するといった最適化を行っています。)

この種の最適化は、実行時に生成される中間のデータの生成、消費処理を除去し、実行時のヒープ使用量を減らすなどの効果を目当てとして、行われます。

次に、metamorphismですが、これは、hyloの双対にあたり、cata(foldr)が生成した計算結果を元に、 ana(unfoldr)によりデータを展開するという処理の記述の一般化された形式で、 
meta = ana . cata
となる合成関数です。 

hylo同様、metaもana/cataの融合であることから、同じように説明することが出来ます。 つまり、展開されているデータを畳み込み、その結果から再度、データを展開する。 ところで、そんな関数は、リスト操作上だとどのような関数でしょうか? 

その一つの答えは、リスト操作関数におけるmap関数はcataでありanaであるという事です。 例えば、anaやcataにconsを生成させる関数を渡す事で、cataとana、いずれでもmapの役割を果たすことが出来ます。

この事は、map関数でありがちな、map f . map g = map (f . g)の融合をcata/ana融合、言い換えると hyloで表現できることを意味するわけですが、それと同時に、metaでも実現出来るわけです。 ana/cata融合、つまりmetaもまた、map/map融合を表す射になりうるのです。

つまり、 map f . map g = map (f . g) は、hyloで表すと、
(cata f') . (ana g') = hylo f' g'
であり、 これと同様の意味をmetaで書くと、
(ana f'') . (cata g'') = meta g'' f'' 
と書くことができます。
f'、f''、g'、g''は、それぞれ、元のf、gをhylo/meta用に修正(wrapping)した関数になります。 cata/ana共に、そもそも、foldr/unfoldなので、mapしたい関数を直接渡す事はできません。。。

実際に、Haskellのコードで、リスト上のcata/anaを書くと以下のようになります。 まず、標準のリスト上のcataとanaです。
lcata :: (alpha -> beta -> beta) -> beta -> [alpha] -> beta
lcata f b []     = b
lcata f b (a:as) = f a (lcata f b as)

lana :: (beta -> Maybe(alpha, beta)) -> beta -> [alpha]
lana f b = case f b of
  Just (a, b') -> a:lana f b'
  Nothing      -> []
anaでは、JustとNothingでデータを生成し続けるのか、停止するのか、切り替えています。 ちなみに、これは基本的にフラグなので、代わりに1, 2とそのタプルを使う文献も多いです。(下の参考文献2、3あたりを参照。というかそっちが主流かも?) 今回は、後で紹介するコードに合わせて、Maybeで書きました。

次に直接再帰形式のhylo/metaです。cata/ana同様にHaskellの標準のリストを対象にしています。
lhylo :: beta -> (alpha -> Maybe(gam, alpha)) -> (gam -> beta -> beta) -> alpha -> beta
lhylo init f g a = case f a of
  Nothing     -> init
  Just (x, y) -> g x (lhylo init f g y)

lmeta :: (gam -> Maybe (beta, gam)) -> (alpha -> gam -> gam) -> gam -> [alpha] -> [beta]
lmeta f g c x = case f c of
  Just(b, c') -> b:(lmeta f g c' x)
  Nothing     -> case x of
    []   -> []
    a:x' -> lmeta f g (g a c) x'
hyloの直接再帰の形式はよく見るのですが、metaについては、ネットではあまり見かけません。。。上記のmetaの形式は、参考文献1の10pの直接再帰形式のコードからHaskellのリスト向けに復元したものです。

このcata/anaで書いたコードに対するmap関数は以下のようになります。元はfoldr/unfoldだけあって、mapしたい関数を直接渡す事はできません。。。
mapC f = lcata func []
  where
    func x xs = (f x):xs

mapA f = lana func
  where
    func []     = Nothing
    func (x:xs) = Just(f x, xs)
次に、hylo/metaによる(map f) . (map g) = map (f . g)の合成。それぞれ、hylo版がmapH、meta版がmapEになります。
mapH f g = lcata func1 [] . lana func2
  where
    func1 x xs = (f x):xs
    func2 []     = Nothing
    func2 (x:xs) = Just(g x, xs)

mapE f g = lana func2 . lcata func1 []
  where
    func1 x xs = (g x):xs
    func2 []     = Nothing
    func2 (x:xs) = Just(f x, xs)
しかし、この例はあまりに面白くありませんね。。。

mapをrecursion schemeで表すなら前述したようにcataかana十分なはずで、 そもそもmap自体、簡単に融合出来てしまいます。 致命的なのが、metaの特徴でもある畳込み後の値が使えるという性質を上手く活かしきれていません。

しかし、上記の例から分かることは、foldで畳み込んだ結果をリストとすることで、metamorphismもそれなりに使いみちが出てくるかも知れません。

というわけで、やってみたのが、metamorphismによるfilterとmapの融合です。
mapfilterM f g = lana func2 . lcata func1 []
  where
    func1 x xs   = if f x then x:xs else xs
    func2 []     = Nothing
    func2 (x:xs) = Just(g x, xs)
これを直接再帰の形式のmetaに差し替えると、こうなります。
mapfilterMF f g = lmeta func2 func1 []
  where
    func1 x xs   = if f x then x:xs else xs
    func2 []     = Nothing
    func2 (x:xs) = Just(g x, xs)
fにfilter用の関数を渡し、gにmapに渡す関数を渡します。
map after filter(filter関数の実結果にmap関数を実行)の融合変換は、map関数を展開があるため、metamorphismを使うことで綺麗に表現出来ています。。。(多分

以下実行結果です。
*Main> mapfilterMF (\x-> (mod x 5) == 0) (3 +) [1..20]
[8,13,18,23]
*Main> map (3 +) $ filter (\x-> (mod x 5) == 0) [1..20]
[8,13,18,23]
*Main> mapfilterM (\x-> (mod x 5) == 0) (3 +) [1..20]
[8,13,18,23]
このようにして、hylomorphismやmetamorphismは融合変換に使用することができます。

そして、最後に! ググったら出てきた、面白いhylomorphism/metamorphismの例を紹介します。(参考文献1より)
それが、hylomorphismとmetamorphismによるQuicksort & Heapsortです。 例によって、in-placeではないので、それぞれ偽のQuicksort、Heapsortになります。

まず、hylomorphismによるQuicksort(参考文献1からの引用)!
data Tree alpha = Node(Maybe(alpha, Tree(alpha), Tree(alpha))) deriving Show

foldt :: (Maybe(alpha, beta, beta) -> beta) -> Tree(alpha) -> beta
foldt f (Node Nothing)         = f Nothing
foldt f (Node (Just(a, t, u))) = f $ Just(a, foldt f t, foldt f u)

unfoldt :: (beta -> Maybe(alpha, beta, beta)) -> beta -> Tree(alpha)
unfoldt f b = case f b of
  Nothing         -> Node Nothing
  Just(a, b1, b2) -> Node(Just(a, unfoldt f b1, unfoldt f b2))

partition :: (Ord alpha) => [alpha] -> Maybe(alpha, [alpha], [alpha])
partition []          = Nothing
partition (a:as)      = Just(a, filter (< a) as, filter (> a) as)

join :: Maybe (alpha, [alpha], [alpha]) -> [alpha]
join Nothing          = []
join (Just(a, x, y))  = x ++ [a] ++ y

quicksort :: (Ord alpha) => [alpha] -> [alpha]
quicksort = foldt join . unfoldt partition
主に、partitionとjsonそして、unfoldt/foldtは木構造を生成して、それをリストに引き戻している処理になります。 Haskellのfake Quicksortのコードがjoinとpartitionの箇所にそのまま現れていることがわかります。

同様にして、metamorphismによるHeapsort(Quicksort同様、参考文献1からのほぼ引用ですが、一部修正)。
insert :: (Ord alpha) =>  alpha -> Tree(alpha) -> Tree(alpha)
insert a t = merge (Node(Just(a, e, e)), t)
  where
    e = Node Nothing

splitMin ::  (Ord alpha) => Tree(alpha) -> Maybe(alpha, Tree(alpha))
splitMin (Node t) = case t of
  Nothing       -> Nothing
  Just(a, u, v) -> Just(a, merge(v, u))

merge :: (Ord alpha) => (Tree(alpha), Tree(alpha)) -> Tree(alpha)
merge(t, Node Nothing) = t
merge(Node Nothing, u) = u
merge(Node x, Node y)  = if a < b
  then Node(Just(a, t2, merge(t1, Node y)))
  else Node(Just(b, u2, merge(u1, Node x)))
  where
    Just(a, t1, t2) = x
    Just(b, u1, u2) = y

heapsort :: (Ord alpha) => [alpha] -> [alpha]
heapsort = lana splitMin . lcata insert (Node Nothing)
こちらは少し特殊ですが、補助関数mergeを作っています。 mergeが別に再帰的な処理を行っていますが、これは本体の再帰とは別の処理です。

まあ、ここまで来ると、素直に再帰を書けば良い気もしますが、再帰関数がこのように分離出来るのは面白いですね。

Haskellのcataやana、hyloなどで書ける関数は、The pointless-haskell package(のExample)に色々紹介してあるので、それらを使って 遊んで見るのも良いかも知れません。

今回作成したコードの全体はここ(Gist)に置いてます。

ps. 急ぎで書いてしまったので、後で修正が入るかも知れません。

参考文献

  1. J.Gibbsons, Metamorphisms: Streaming Representation-Changers, 2005(PDF)
  2. F. Domínguez, Alberto Pardo, Exploiting algebra/coalgebra duality for program fusion extensions, 2011 (PDF)
  3. Yoshiyuki Onoue, Zhenjiang Hu, Hideya Iwasaki, Masato Takeichi, A Calculational Fusion System HYLO, 1997 (PDF)

2017/05/14

Emacsのギリシャ文字/キリル文字のフォント表示

長らくブログを書いてなかったので、復帰がてらに、最近やったEmacsの設定について、メモします。

Emacsを使う場合、日本語の設定を
;; 日本語
(set-fontset-font
 'nil 'japanese-jisx0208 (font-spec :family "Takaoゴシック" :height 90))
のように、何かしらフォントの設定を行うと思いますが、その時、ギリシャ文字やキリル文字のフォントも全角化されてしまいます。以下のように...


しかし、全角のギリシャ/キリル文字は読みづらいので、できれば、別のフォントで表示させたくなります。というわけで、こんな感じで別のコードポイントを割り当てます。
;; ギリシャ文字
(set-fontset-font
 'nil '(#x0370 . #x03FF) (font-spec :family "Ubuntu Mono" :height 100))

;; キリル文字
(set-fontset-font
 'nil '(#x0400 . #x04FF) (font-spec :family "Ubuntu Mono" :height 100))

ギリシャ文字のコードポイントは、#x0370〜#x03FF、キリル文字のコードポイントは、#x0400〜#x04FFなので、上記のように設定します。
そして、再表示した結果が以下。

しっかり、半角化されて、綺麗な表示になりました。

ちなみに、私はこんなフォント設定をしています。
(defvar default-font-family "Ubuntu Mono")
(defvar default-font-family-jp "Takaoゴシック")

;; デフォルトフォント設定
(set-face-attribute
 'default nil :family default-font-family :height 100)

;; 日本語のフォントセット : あいうえお ... 日本語
(set-fontset-font
 'nil 'japanese-jisx0208 (font-spec :family default-font-family-jp :height 90))

;; ギリシャ文字のフォントセット : αβγκλ ... ΛΩ
(set-fontset-font
 'nil '(#x0370 . #x03FF) (font-spec :family default-font-family :height 100))

;; キリル文字のフォントセット : Эта статья ... Русский
(set-fontset-font
 'nil '(#x0400 . #x04FF) (font-spec :family default-font-family :height 100))