Brownian treeのシミュレーションしてからグラフ理論で最短経路を求める

ブラウン樹の結晶成長シミュレーションを眺めていたら、すべてのノードが連結されているので、難易度はともかく迷路として取り扱えるかなと。
とりあえず樹を目一杯成長させてから、ある端から別の端までの経路をグラフ理論を使って経路を求めてみます。

目的としては、

* 格子上のブラウン運動程度の簡単なシミュレータでpythonの練習
* networkxの使い方

を考えています。

# Brownian tree

## 系の設定

高さ`H`マス、幅`W`マスの盤面を考えます。

```python
import random
import numpy as np

class BTree(object):
    def __init__(self, H, W):
        self.N = H*W
        self.H = H
        self.W = W
        self.tree = np.zeros((H, W), dtype=int)
        self.root = (random.randint(0, W - 1), random.randint(0, H - 1))
        self.tree[self.root[1], self.root[0]] = 1
        self.edges = []
```

* `tree`は結晶がいるマスが`1`、空きマスが`0`になる行列
* `root`は種になるマスの座標の(x, y)タプル、とりあえず乱数で決定
* `edges`は結晶が成長した方向を保存、(x1, y1, x2, y2)のタプルをリストで持つ

## 結晶の成長

結晶の成長ステップは以下の通り。

* 空きマスにランダムに自由粒子が置かれる。
* 自由粒子はランダムに上下左右に動く。システムの端より外には行かない。
* `tree`が`1`になっているマス(x1, y1)の隣に来たら運動をやめる。
  - 粒子がいたところ(x2, y2)について、`tree[y2, x2]`を`1`にする。
  - 複数の`tree = 1`と隣接していたら、その中からランダムで(x1, y1)を選択。
  - `edges`に((x1, y1), (x2, y2))を追加。
* `tree`がすべて`1`になったら終了。

通常は見たい対象に応じて終了条件を設定すると思います。例えばブラウン樹が系の端に達したら終了、など。今回は迷路を作りたいというだけなので、すべて埋まるまで粒子を投入し続けます。

また、今回は自由粒子の移動方向と隣接マスの方向は同じ上下左右とします。

```python
move = [
    np.array([ 1,  0], dtype=int), 
    np.array([ 0,  1], dtype=int), 
    np.array([-1,  0], dtype=int), 
    np.array([ 0, -1], dtype=int), 
]

neighbor_grid = move

class BTree(object):
    def grow(self):
        free_place = np.array(np.where(self.tree == 0)).T
        p = random.choice(free_place)[-1::-1]
        while not self._check_touched_at(p):
            p = p + random.choice(move)
            p[p < 0] = 0
            p[0] = min(p[0], self.W - 1)
            p[1] = min(p[1], self.H - 1)
        x, y = p
        self.tree[y, x] = 1

    def _check_touched_at(self, p):
        qs = []
        for dx in neighbor_grid:
            q = p + dx
            if q[0] < 0 or q[1] < 0 or q[1] >= self.H or q[0] >= self.W:
                continue
            if self.tree[q[1], q[0]] == 1:
                qs.append(q)
        if len(qs) > 0:
            self.edges.append((p, random.choice(qs)))
            return True

        return False

    def is_grown(self):
        return np.all(self.tree == 1)
```
ざっとポイントを確認します。

`free_place = np.array(np.where(self.tree == 0)).T` は、`tree`の中で`0`である座標がN行２列(つまり[[y0, x0], [y1, x1], ... ])という形で得られます。
`_check_touched_at()`は、自由粒子が座標`p`に来た時に、結晶に触れているかどうかをチェックし、触れていたら`True`を返します。
同時に`self.edges`に新しく生えたedgeを追加します。本当はそれ用の関数を用意したかったのですが、せっかく見つけた隣接ノードデータ`qs`を捨てたくなかったので。設計をミスりました。
それはともかく、自由粒子が成長した結晶の枝に触れた時、複数のノードに同時に触れる可能性もあります。そこで全ての触れたノードに枝を生やさず、ランダムに１つ選んでedgeを生やすことにします。それが`self.edges.append((p, random.choice(qs)))`です。
これで、`grow()`を実行すると、１つの粒子が結晶に追加される処理が行われます。

`is_grown()`は結晶が成長しきったかどうかをチェックします。今回は結晶がすべてのマスを埋め尽くしたかどうかとしました。


## 可視化

成長の方向を踏まえて、枝っぽく伸ばしていこうかなと。わざわざ`edges`を定義したのはこのため。

```python:BTree
import cv2

...

    def make_image(self, grid_px=11, edge_width=5):
        w = self.W*grid_px
        h = self.H*grid_px
        img = np.zeros((h, w, 3), dtype=np.uint8)
        center = int(grid_px*0.5 + 0.5)

        cx = self.root[0]*grid_px + center
        cy = self.root[1]*grid_px + center
        cv2.circle(img, (cx, cy), grid_px, (0, 120, 120), -1)

        for e in self.edges:
            p1, p2 = e
            c1 = p1*grid_px + center
            c2 = p2*grid_px + center
            cv2.line(img, tuple(c1), tuple(c2), (0, 200, 0), edge_width)
        return img

    def add_edge_to_image(self, img, last_n_edge=1, grid_px=11, edge_width=5):
       # 略。last_n_edgeまで描写した画像に対し、それ以降のedgeを追加で書き足す。
       # 成長の１ステップごとに画像を作る場合、毎回１から作るよりメモリ・処理の節約を目指した。
```

## 初期化

系の初期化を行う関数もとりあえず作っておきます。

```python:BTree
    def reset(self):
        self.tree = np.zeros((self.H, self.W), dtype=int)
        self.root = (random.randint(0, self.W - 1), random.randint(0, self.H - 1))
        self.tree[self.root[1], self.root[0]] = 1
        del self.edges
        self.edges = []
```

## 実行

```python:brownian_tree.py

import sys
import random
import numpy as np
import cv2

# 略。以下を定義する
# move
# neighbor_grid
# BTree(object)


def main():
    T = BTree(24, 32)
    while True:
        img = T.make_image()
        while not T.is_grown():
            img = T.add_edge_to_image(img)
            cv2.imshow("test", img)
            cv2.waitKey(5)
            T.grow()
        print("grown!")
        img = T.add_edge_to_image(img)
        cv2.imshow("test", img)
        cv2.waitKey(1000)
        T.reset()
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

```

これで以下のようなアニメーションを延々と繰り返します。

![test.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/12122134-8fe1-623e-353b-3b2122753974.gif)


# 経路探索

まあ一言で言えばDijkstraのアルゴリズムを使うだけなのですが、networkxの使い方の確認も兼ねてソースコードをざっと眺めます。

## ブラウン樹の準備

```python
import networkx as nx
import numpy as np
import cv2
from brownian_tree import BTree

def main():
    H = 48 
    W = 64
    T = BTree(H, W)

    img = T.make_image()
    while not T.is_grown():
        img = T.add_edge_to_image(img)
        cv2.imshow("test", img)
        cv2.waitKey(1)
        T.grow()
```

先の節のような結晶が育つアニメーションを眺めて待ちます。

## グラフの準備

```python
    H = 48 
    W = 64
    T = BTree(H, W)

    while not T.is_grown():
        T.grow()
```
ついでに、左上に入り口と右下の出口用の枝を追加します。

```python
    T.edges.append((np.array([0, -1], dtype=int), np.array([0, 0], dtype=int)))
    T.edges.append((np.array([W - 1, H - 1], dtype=int), np.array([W - 1, H], dtype=int)))
```

この`T.edges`からデータを作ります。

```python
import networkx as nx

（中略）

    G = nx.Graph()
    for e in T.edges:
        u1, u2 = e
        index1 = u1[1]*W + u1[0]
        index2 = u2[1]*W + u2[0]
        G.add_edge(index1, index2)
```
`T.edges`の要素を１つ取ってみると`((x1, y1), (x2, y2))`となって、
１つのエッジの端から端になります。
エッジの端はノードです。ここでノードを`Graph`で管理できるように、ノードの２次元座標`(x1, y1)`を通し番号に変換します。それが`index1 = u1[1]*W + u1[0]`で、数式にすると $i=Wy + x $ となります。
マップの上から横方向に通し番号を振ったものですね。

なお、先ほどわざわざ追加した入口と出口は、入り口の通し番号が$-W$、出口の番号が$H(W + 1) - 1$ となっています。


## 経路の検出

左上のマス $0$ から、右下のマス $WH - 1$ への経路を計算します。

```python
   route = nx.dijkstra_path(G, 0, W*H - 1)
```

`route`には`[i1, i2, i3, ... , (ゴール手前, ゴール)]`という感じでデータが入っています。

ついでに、入口、出口のノードも経路に追加しておきます。

```python
    route = [-W] + route + [W*H + W - 1]
```

## 可視化

`route` から $n$ 番目と $n + 1$ 番目を取り出した $r_n$ と $r_{n+1}$ を結んだedgeは経路となっています。
これを順番に描写させると迷路を通っている映像になります。

```python
    img = T.make_image()
    cv2.imshow("test", img)
    cv2.waitKey(1)
    count = 0
    grid_px = 11
    edge_width = 5
    center = int(grid_px*0.5 + 0.5)
    for u1, u2 in zip(route[:-1], route[1:]):
        p1 = ((u1%W)*grid_px + center, int(u1/W)*grid_px + center)
        p2 = ((u2%W)*grid_px + center, int(u2/W)*grid_px + center)
        cv2.line(img, p1, p2, (200, 0, 0), edge_width - 2)
        cv2.imshow("test", img)
        cv2.waitKey(20)
    cv2.imshow("test", img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
```

こうすると、成長したbrownian treeを迷路としたときの経路を、青い線が通って行く様子が見られます。

![maze_test.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/0ee5cba7-8b35-f348-73c1-4b80708f7c40.gif)

# その他１

かなり前に、Javascriptでサンゴっぽいものを成長させる遊びをしたことがあります。Three.js使用。
http://jsdo.it/sage_k/0yKG

# その他２

迷路を作るなら迷路用のアルゴリズムを使いましょう。