フラクタル図形 ー バーンズリーのシダをPythonで描く

今年5月に掲載した、コッホ曲線、シェルピンスキーのギャスケットに続き、久しぶりにPythonでフラクタル図形を描画する。最終的にはジュリア集合やマンデルブロ集合の描画をやってみようと思っていた（実際に6月にジュリア集合をやりかけていたが途中かけになっていた）が、今回はバーンズリーのシダ（Barnsley fern）を描いてみた。

概要

バーンズリーのシダはマイケル・バーンズリーというイギリスの数学者が著書『Fractals Everywhere』の中で発表したもので、ジョージア工科大学の講義内容がもとになっている（Wikipediaより）らしい。

バーンズリーのシダの描き方

アフィン変換

• 平行移動
• \( T_{x} \)：　\( x軸 \)方向への移動
• \( T_{y} \)：　\( y軸 \)方向への移動
• 拡大・縮小
• \( a \)：　\( x軸 \)方向の拡大・縮小
• \( d \)：　\( y軸 \)方向の拡大・縮小
• 回転

• 反時計回りに\( \theta^\circ \)回転
• スキュー（変形）

• \( x \)軸から反時計回りに\( \alpha^\circ \)ねじり
• \( y \)軸から時計回りに\( \beta^\circ \)ねじり

バーンズリーのシダの変換係数

シダを茎・小さい葉・大きい葉（左右）の4パターンに分けて、アフィン変換の係数と実行する確率が設定されている。計算1回毎に乱数を生成して、これら4パターンのどれを用いて計算するかを決定し、ひたすら反復計算を行っていくとシダの模様になっていく。

茎（1%の確率で実行）

\( f_{1}(x,y)=\begin{bmatrix}0&0\\0&0.16\end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}0\\0\end{bmatrix} \)

\( \left\{ \begin{eqnarray}x_{n+1}&=&0\\ y_{n+1}&=&0.16y_{n} \end{eqnarray} \right. \)

\( f_{1} \) では、 \( x \)はゼロに \( y \) は0.16倍されることで、茎の付け根付近が描かれる。1%の確率で描画がリセットされるが原点 \( (0, 0) \) には戻らず、\( f_{2}, f_{3}, f_{4} \) の葉の描画の基となる点が生成される。

小さい葉（85%の確率で実行）

\( f_{2}(x,y)=\begin{bmatrix}0.85&0.04\\-0.04&0.85\end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}0\\1.60\end{bmatrix} \)  

\( \left\{ \begin{eqnarray}x_{n+1}&=&0.85x_{n}+0.04y_{n}\\ y_{n+1}&=&-0.04x_{n}+0.85y_{n}+1.60 \end{eqnarray} \right. \)

\( f_{2} \) では、時計回りに少し回転・少し縮小・y軸方向プラス側に移動させる動作をすることで、シダ全体を描いている。

左側の大きい葉（7%の確率で実行）

\( f_{3}(x,y)=\begin{bmatrix}0.20&-0.26\\0.23&0.22\end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}0\\1.60\end{bmatrix} \)

\( \left\{ \begin{eqnarray}x_{n+1}&=&0.20x_{n}-0.26y_{n}\\ y_{n+1}&=&0.23x_{n}+0.22y_{n}+1.60 \end{eqnarray} \right. \)

\( f_{3} \) では、反時計回りに90度くらい回転・30%程度に縮小・y軸方向プラス側に移動の動作により左側最下段の葉が描画される。

右側の大きい葉（7%の確率で実行）

\( f_{4}(x,y)=\begin{bmatrix}-0.15&0.28\\0.26&0.24\end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}0\\0.44\end{bmatrix} \)

\( \left\{ \begin{eqnarray}x_{n+1}&=&-0.15x_{n}+0.28y_{n}\\ y_{n+1}&=&0.26x_{n}+0.24y_{n}+0.44 \end{eqnarray} \right. \)

Pythonで描画

Matplotlibにて描画。単純な繰り返し計算ではあるが、計算時間はまずまずかかる。（100万回では約20分くらい要した）

import matplotlib.pyplot as plt
import random

# Init
p = 0
x = 0
y = 0
iteration = 100000     #計算回数

def BarnsleyFern(xn,yn,p):
    if p < 0.01:
        x = 0
        y = 0.16 * yn
    elif p < 0.86:
        x = 0.85 * xn + 0.04 * yn
        y = -0.04 * xn + 0.85 * yn + 1.6
    elif p < 0.93:
        x = 0.20 * xn - 0.26 * yn
        y = 0.23 * xn + 0.22 * yn + 1.6
    else:
        x = -0.15 * xn + 0.28 * yn
        y = 0.26 * xn + 0.24 * yn + 0.44
    return x,y

# Plotの設定
fig = plt.figure(figsize=(5, 5))

for n in range(iteration):
    p = random.random()
    (x,y) = BarnsleyFern(x, y, p)
    plt.plot(x,y,'.',markersize=0.1, markeredgecolor='g', markerfacecolor='g')

結果

100回計算 (markersize = 5)

500回計算 (markersize = 5)

1,000回計算 (markersize = 5)

5,000回計算 (markersize = 1)

10,000回計算 (markersize = 1)

50,000 (markersize = 0.5)

100,000回計算 (markersize = 0.1)

250,000回計算 (markersize = 0.1)

1,000,000回計算 (markersize = 0.01)

参考にしたサイト