LoginSignup
2
1

More than 5 years have passed since last update.

@sage-git

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

Last updated at Posted at 2018-06-04

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

目的としては、

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

を考えています。

Brownian tree

系の設定

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

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になったら終了。

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

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

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を定義したのはこのため。

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を追加で書き足す。
       # 成長の１ステップごとに画像を作る場合、毎回１から作るよりメモリ・処理の節約を目指した。

初期化

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

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 = []

実行

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

経路探索

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

ブラウン樹の準備

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()

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

グラフの準備

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

    while not T.is_grown():
        T.grow()

ついでに、左上に入り口と右下の出口用の枝を追加します。

    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からデータを作ります。

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$ への経路を計算します。

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

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

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

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

可視化

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

    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

その他１

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

その他２

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

2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1