Python + matplotlibを使ってバーンズリーのシダを描画してみます。
バーンズリーのシダとは
バーンズリーのシダ（Barnsley's Fern）とはイギリスの数学者のマイケル・バーンズリーさん（Barnsley Michael）にちなんだ1988年ごろに発表された数学的に生成されるフラクタルのパターンです。シダの植物のようなパターンが描画されます。
※写真はロイヤルティフリーの木陰で共存する白くて小さな花とシダ植物の写真を拡大・トリミング（JPG画像）より。
どういった計算なのか
フラクタルの計算の一種で、前の座標点から確率に応じた変化処理が反映され新しい座標が算出され、結果的にシダの植物のような形状のパターンを生成します。
※フラクタルに関しては以下のWikipediaの記事などを必要に応じてご確認ください。
シダの生成に使われる変換処理は以下の4つが存在します。それぞれ確率が設けられており、確率に応じてランダムに対象の変換が4つの中から選択されます。$x$と$y$は現在のX座標とY座標、$x_1$と$y_1$は変換処理反映後のXとYの座標を表しています。
1つ目の変換（85%の確率で選択）:
x_1 = 0.85x + 0.04y\\
y_1 = -0.04x + 0.85y + 1.6
2つ目の変換（7％の確率で選択）:
x_1 = 0.2x - 0.26y\\
y_1 = 0.23x + 0.22y + 1.6
3つ目の変換（7％の確率で選択）:
x_1 = -0.15x + 0.28y\\
y_1 = 0.26x + 0.24y + 0.44
4つ目の変換（1％の確率で選択）:
x_1 = 0\\
y_1 = 0.16y
上記変換処理を数千回、ケースによっては数万回繰り返してパターンを生成します。それぞれの変換がどんな雰囲気なのかは後の節で触れていきます。
使うもの
- Python 3.9.0
- Jupyter (VS Code上のものを利用)
- matplotlib==3.4.3
検証で1つ目の変換のみのコードを書いてみる
大雑把に1つ目の変換処理がどんな座標の推移をするのか、1つ目の変換だけ反映する形でコードを書いてみます。
import random
import matplotlib.pyplot as plt
from typing import TypedDict, List
class Point(TypedDict):
"""XとY座標単体を扱うための辞書。
"""
x: float
y: float
def apply_transformation_1(point: Point) -> Point:
"""
1つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.85 * x + 0.04 * y
y1: float = -0.04 * x + 0.85 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
def draw_fern(n: int) -> None:
"""
バーンズリーのシダの計算を行い結果をmatplotlibで描画する。
Parameters
----------
n : int
座標変換処理の実行回数。
"""
# 変換処理が反映された座標を格納するリストを作成しています。
# 最初の座標（最初のインデックス）はXと両方とも0からスタートしています。
x_list: List[float] = [0]
y_list: List[float] = [0]
current_point: Point = {'x': x_list[0], 'y': y_list[0]}
for _ in range(n):
current_point = apply_transformation_1(point=current_point)
x_list.append(current_point['x'])
y_list.append(current_point['y'])
plt.figure(figsize=(6, 6))
plt.plot(x_list, y_list, 'o')
plt.show()
N: int = 3000
draw_fern(n=N)
まずはXとY座標単体の組み合わせを保持するためのPointというTypedDictの辞書のクラスを作っています。
参考 :
また、座標変換用にapply_transformation_1という関数を設けています。この関数では以下の数式部分（1つ目の変換式）を担当します。
x_1 = 0.85x + 0.04y\\
y_1 = -0.04x + 0.85y + 1.6
draw_fern関数では渡された引数の回数分ループを行い変換処理（今回はapply_transformation_1）を反映していきます。この関数内ではループが終わったら結果の各X座標とY座標のリスト（x_listとy_list）でmatplotlibを使って可視化しています。今回は単一点を丸で描画したかったのでplotの第三引数は'o'としています。
N: int = 3000とし、3000回変換処理を行うようにしています。
実行してみると以下のようにプロットされます。
XとY共に0の座標からスタートしているので、この変換式では右上の方向に座標が移動していっていることが分かります。
検証で2つ目の変換の処理を追加してみる
続いて2つ目の変換処理を加えていきましょう。本来の確率ではありませんが、大体の確率の割り振りとして1つ目の変換は89％、2つ目の変換は11％で設定しています。
import random
import matplotlib.pyplot as plt
from typing import TypedDict, List, Protocol
class Point(TypedDict):
"""XとY座標単体を扱うための辞書。
"""
x: float
y: float
class Transformation(Protocol):
def __call__(self, point: Point) -> Point:
...
def apply_transformation_1(point: Point) -> Point:
"""
1つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.85 * x + 0.04 * y
y1: float = -0.04 * x + 0.85 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_2(point: Point) -> Point:
"""
2つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.2 * x - 0.26 * y
y1: float = 0.23 * x + 0.22 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_randomly(point: Point) -> Point:
"""
確率に応じたいずれかの座標変換処理を反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
probabilities: List[float] = [0.89, 0.11]
transformations: List[Transformation] = [
apply_transformation_1, apply_transformation_2]
random_value: float = random.random()
current_probability: float = 0.0
for probability, transformation in zip(probabilities, transformations):
current_probability += probability
if random_value <= current_probability:
point = transformation(point=point)
return point
point = transformations[-1](point=point)
return point
def draw_fern(n: int) -> None:
"""
バーンズリーのシダの計算を行い結果をmatplotlibで描画する。
Parameters
----------
n : int
座標変換処理の実行回数。
"""
# 変換処理が反映された座標を格納するリストを作成しています。
# 最初の座標（最初のインデックス）はXと両方とも0からスタートしています。
x_list: List[float] = [0]
y_list: List[float] = [0]
current_point: Point = {'x': x_list[0], 'y': y_list[0]}
for _ in range(n):
current_point = apply_transformation_randomly(point=current_point)
x_list.append(current_point['x'])
y_list.append(current_point['y'])
plt.figure(figsize=(6, 6))
plt.plot(x_list, y_list, 'o')
plt.show()
N: int = 3000
draw_fern(n=N)
コードとしては1つ目の時点と同じところも多いのですが、以下の変更を加えています。
まずはProtocolクラスを使った変換の各関数の型アノテーション用のクラスを設けています。ミスを減らすために使っているので、型アノテーションをしない場合などには不要です。
class Transformation(Protocol):
def __call__(self, point: Point) -> Point:
...
参考 :
また、2つ目の変換用の関数を追加してあります。
def apply_transformation_2(point: Point) -> Point:
"""
2つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.2 * x - 0.26 * y
y1: float = 0.23 * x + 0.22 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
この関数は以下の数式部分を担当しています。
x_1 = 0.2x - 0.26y\\
y_1 = 0.23x + 0.22y + 1.6
そしてランダムに変換の関数を反映する関数も追加しています。
def apply_transformation_randomly(point: Point) -> Point:
"""
確率に応じたいずれかの座標変換処理を反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
probabilities: List[float] = [0.89, 0.11]
transformations: List[Transformation] = [
apply_transformation_1, apply_transformation_2]
random_value: float = random.random()
current_probability: float = 0.0
for probability, transformation in zip(probabilities, transformations):
current_probability += probability
if random_value <= current_probability:
point = transformation(point=point)
return point
point = transformations[-1](point=point)
return point
2つの変換における確率の割り振りは[0.89, 0.11]という箇所で指定しています。あとはビルトインのrandomパッケージを使って乱数を計算し得られた値に応じた変換を行っています。
あとはdraw_fern関数内で追加したapply_transformation_randomly関数を呼び出して変換をかけています。
...
for _ in range(n):
current_point = apply_transformation_randomly(point=current_point)
x_list.append(current_point['x'])
y_list.append(current_point['y'])
...
他は同じです。実行してプロットしてみます。
若干右側にも座標が増えたりしていますが、主に左側に座標点が増えた・・・という所感です。
検証で3つ目の変換の処理を追加してみる
続いて3つ目の変換処理も加えていきます。コードはほとんど2つ目と一緒で、以下の点のみ対応します。
- 3つ目の変換の関数を追加します。
- 各関数での確率割り振りに3つ目の変換を加えます。
まずは3つ目の変換の関数を追加します。
def apply_transformation_3(point: Point) -> Point:
"""
3つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = -0.15 * x + 0.28 * y
y1: float = 0.26 * x + 0.24 * y + 0.44
point['x'] = x1
point['y'] = y1
return point
上記の変換の関数は以下の数式部分に該当します。
x_1 = -0.15x + 0.28y\\
y_1 = 0.26x + 0.24y + 0.44
あとは各変換の確率の割り振りを設定するだけです。この時点では1つ目が大雑把に85.5％、2つ目が7.25％、3つ目も7.25％としています。
def apply_transformation_randomly(point: Point) -> Point:
"""
確率に応じたいずれかの座標変換処理を反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
probabilities: List[float] = [0.855, 0.0725, 0.0725]
transformations: List[Transformation] = [
apply_transformation_1, apply_transformation_2,
apply_transformation_3]
実行してプロットしてみると右側の方のの座標が充実したような挙動になりました。
4つ目の変換の処理を追加してみる
最後に4つ目の変換処理を加えていきます。やることは3つ目の追加時とほぼ同じです。
まずは4つ目の変換用の関数を追加します。
def apply_transformation_4(point: Point) -> Point:
"""
4つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
y: float = point['y']
x1: float = 0
y1 = 0.16 * y
point['x'] = x1
point['y'] = y1
return point
この関数は以下の数式部分に該当します。
x_1 = 0
y_1 = 0.16y
あとは確率の割り振りを本来のものに調整します。本来の確率は1つ目が85％、2つ目が7％、3つ目が7％、4つ目が1％となります。
def apply_transformation_randomly(point: Point) -> Point:
"""
確率に応じたいずれかの座標変換処理を反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
probabilities: List[float] = [0.85, 0.07, 0.07, 0.01]
transformations: List[Transformation] = [
apply_transformation_1, apply_transformation_2,
apply_transformation_3, apply_transformation_4]
...
x = 0の箇所に座標が増えて茎のような雰囲気が追加になりました。これで完成です！
コード全体
以下最終的な全体のコードです。
import random
import matplotlib.pyplot as plt
from typing import TypedDict, List, Protocol
class Point(TypedDict):
"""XとY座標単体を扱うための辞書。
"""
x: float
y: float
class Transformation(Protocol):
def __call__(self, point: Point) -> Point:
...
def apply_transformation_1(point: Point) -> Point:
"""
1つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.85 * x + 0.04 * y
y1: float = -0.04 * x + 0.85 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_2(point: Point) -> Point:
"""
2つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = 0.2 * x - 0.26 * y
y1: float = 0.23 * x + 0.22 * y + 1.6
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_3(point: Point) -> Point:
"""
3つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
x: float = point['x']
y: float = point['y']
x1: float = -0.15 * x + 0.28 * y
y1: float = 0.26 * x + 0.24 * y + 0.44
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_4(point: Point) -> Point:
"""
4つ目の変換処理を指定された座標に反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
y: float = point['y']
x1: float = 0
y1 = 0.16 * y
point['x'] = x1
point['y'] = y1
return point
def apply_transformation_randomly(point: Point) -> Point:
"""
確率に応じたいずれかの座標変換処理を反映する。
Parameters
----------
point : Point
反映対象の座標。
Returns
-------
point : Point
反映処理後の座標。
"""
probabilities: List[float] = [0.85, 0.07, 0.07, 0.01]
transformations: List[Transformation] = [
apply_transformation_1, apply_transformation_2,
apply_transformation_3, apply_transformation_4]
random_value: float = random.random()
current_probability: float = 0.0
for probability, transformation in zip(probabilities, transformations):
current_probability += probability
if random_value <= current_probability:
point = transformation(point=point)
return point
point = transformations[-1](point=point)
return point
def draw_fern(n: int) -> None:
"""
バーンズリーのシダの計算を行い結果をmatplotlibで描画する。
Parameters
----------
n : int
座標変換処理の実行回数。
"""
# 変換処理が反映された座標を格納するリストを作成しています。
# 最初の座標（最初のインデックス）はXと両方とも0からスタートしています。
x_list: List[float] = [0]
y_list: List[float] = [0]
current_point: Point = {'x': x_list[0], 'y': y_list[0]}
for _ in range(n):
current_point = apply_transformation_randomly(point=current_point)
x_list.append(current_point['x'])
y_list.append(current_point['y'])
plt.figure(figsize=(6, 6))
plt.plot(x_list, y_list, 'o')
plt.show()
N: int = 3000
draw_fern(n=N)
参考文献・参考サイトまとめ