乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

フラクタル図形

昨日に引き続きフラクタル図形の話です。

Pythonには便利なライブラリが有るので
それを使ってみましょう。

解説動画はこちら




Pythonには
タートル(turtle)と言う描画用の
ライブラリがあります。

これは亀さんを操作して
亀さんが色々作図してくれるという
ライブラリです(笑)。


ライブラリの呼び出しは
from turtle import *

主なメソッドとしては
fd(距離):タートルを距離分前に進ませ線を描く
rt(角度):タートルの向きを右向きに角度分変える
circle(大きさ):円を描く
done():最後につける

などがあります。

早速作図してみましょう。

四角形を描画するなら
from turtle import *

for _ in range(4):
    fd(100)
    rt(90)

done()
スクリーンショット 2020-12-13 16.53.22
多角形を描くなら
角度を変えてあげます。
from turtle import *

N=12
for _ in range(N):
    fd(50)
    rt(360//N)

done()
スクリーンショット 2020-12-13 16.53.35
円形に円を並べて描いてみます。
circle()で円を小さな描いています。
from turtle import *

speed(0)
pu()
setpos(0, 50)
pd()
pensize(4)

n = 30
angle = 360//n

for i in range(n):
    pencolor(((i/n), 0, 1.0-(i/n)))
    fd(15)
    rt(angle)
    circle(30)

done()
スクリーンショット 2020-12-13 16.53.48
色々組み合わせると面白いものが描けます。
六角形を並べてオウムガイっぽいのを
描いてみます。
from turtle import *

def hexagon(length):
    for _ in range(6):
        fd(length)
        rt(60)

def spindle(n, length, angle):
    for _ in range(n):
        hexagon(length)
        rt(angle)
        length = length*1.05

n = 100
length = 3
angle = 10
speed(0)

spindle(n, length, angle)
done()
スクリーンショット 2020-12-13 16.54.20
次は森っぽいものです。
三角形の木を描きますが
その途中で三角形を再帰を用いて
付け足していくと
木から森っぽくなります。
from turtle import *

def forest(n, length=1500):
    if n <= 0:
        fd(length)
        return
    forest(n-1, length * 0.5)
    rt(80)
    forest(n-1, length / 3)
    lt(160)
    forest(n-1, length / 3)
    rt(80)
    forest(n-1, length * 0.3)

pu()
setpos(500, -200)
pd()
lt(180)
speed(0)
forest(5)
done()
スクリーンショット 2020-12-13 16.55.13
お次はヒルベルト曲線です。
from turtle import *

def hilbert(n, length):
    l_turn(n, length)

def l_turn(n, length):
    if n <= 0:
        return
    rt(90)
    r_turn(n-1, length)
    fd(length)
    lt(90)
    l_turn(n-1, length)
    fd(length)
    l_turn(n-1, length)
    lt(90)
    fd(length)
    r_turn(n-1, length)
    rt(90)

def r_turn(n, length):
    if n <= 0:
        return
    lt(90)
    l_turn(n-1, length)
    fd(length)
    rt(90)
    r_turn(n-1, length)
    fd(length)
    r_turn(n-1, length)
    rt(90)
    fd(length)
    l_turn(n-1, length)
    lt(90)

pu()
setpos(-200, 200)
pd()
speed(0)
pensize(10)
hilbert(4, 20)
done()
スクリーンショット 2020-12-13 16.55.39
再帰を用いて描いています。

お次はシェルピンスキー曲線です。
from turtle import *

def sierpinski(n, length=8, angle=45):
    if n <= 0:
        fd(length)
        return
    rt(angle)
    sierpinski(n-1, length, -angle)
    lt(angle)
    fd(length)
    lt(angle)
    sierpinski(n-1, length, -angle)
    rt(angle)

lt(90)
speed(0)
length=10
pensize(3)
for _ in range(4):
    sierpinski(7, length)
    rt(45)
    fd(length)
    rt(45)
done()
スクリーンショット 2020-12-13 16.56.44
なんて表現すれば良いか分かりませんが
だんだん複雑になってきました。

次はドラゴン曲線です。
ラーメン丼にアリそうなやつです。
from turtle import *

def dragon(n, length=500):
    dragon_r(n, length)

def dragon_r(n, length):
    length /= 1.414
    if n <= 0:
        fd(length)
        return
    rt(45)
    dragon_r(n-1, length)
    lt(90)
    doragon_l(n-1, length)
    rt(45)

def doragon_l(n, length):
    length /= 1.414
    if n <= 0:
        fd(length)
        return
    lt(45)
    dragon_r(n-1, length)
    rt(90)
    doragon_l(n-1, length)
    lt(45)

pu()
setpos(-100, 0)
pd()
speed(0)
pensize(3)
dragon(10)
done()

複雑なやつは描くのに時間がかかりますね。

タートルにはデモスクリプトがついているので
簡単に試して結果を見ることができます。

スクリーンショット 2020-12-13 16.53.02
複雑なフラクタル図形も
有るみたいなので
色々遊んでみると面白いです。

興味がある人は
是非遊んでみてください。

それでは。

最近見ているドラマでフラクタル図形が
取り上げられていたので
気になっちゃいました。

今回は色々なフラクタル図形を作図します。

解説動画はこちら




さてフラクタル図形と言うものは
どう言う図形でしょうか?

フランスの数学者ブノワ・マンデルブロが
導入した幾何学の概念で
図形の一部が全体と
自己相似な構造を持っている
図形のことです。

結構たくさんの種類があります。

Python言語では作図の方法が沢山あります。

最もメジャーな方法としては
matplotlibを用いた作図です。

x,y座標を計算してプロットするものです。


マンデルブロ集合
まず最初はマンデルブロ集合です。
数学者ブノワ・マンデルブロの
名に因む集合体の図形です。

次のようなコードで試すことができます。
from ipywidgets import interact,IntSlider
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

h,w = 1000,1000 
y,x = np.ogrid[ -1.4:1.4:h*1j, -2:0.8:w*1j ]
c = x+y*1j

def mandelbro_plt(rate=30):
    z = c
    times = rate + np.zeros(z.shape, dtype=int)
    for i in range(rate):
        z = z**2 + c
        diverge = z*np.conj(z) > 2**2
        div = diverge & (times==rate)
        times[div] = i
        z[diverge] = 2
    return times

mdb_rate = IntSlider(min=1, max=30, step=1, value=1)
@interact(mdb_rate=mdb_rate)
def plot(mdb_rate):
    plt.figure(figsize=(10,10))
    plt.imshow(mandelbro_plt(mdb_rate))
    plt.show()
download-2

スライダーを上げていくと
より細かい感じになっていきます。

一部を見ると貝のような感じにも
見受けられる図形です。


シェルピンスキーのギャスケット
ポーランドの数学者
ヴァツワフ・シェルピンスキに
因んで名づけられた

自己相似的な無数の
三角形からなる図形です。

matplotlib.patchesの
Polygonで三角形を描画できるので
これを使っていきます。

patches.Polygon(xy = [座標1, 座標2, 座標2])
少し複雑なプログラムですが
スライダーで三角形の
深さを変えることができます。
from ipywidgets import interact,IntSlider
import matplotlib.pyplot as plt
import matplotlib.patches as pat
import numpy as np
%matplotlib inline

def dist(p1,p2):
    return np.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)

def points(p, tri):
    d1,d2,d3 = dist(p, tri[0]),dist(p, tri[1]),dist(p, tri[2])
    if d1 > d2:
        return [p, tri[1], tri[2]] if d1 > d3 else [p, tri[0], tri[1]]
    else:
        return [p, tri[0],tri[2]] if d2 > d3 else [p, tri[0],tri[1]]

def make_tri(tri):
    x1,y1 = (tri[0][0] + tri[1][0])/2,(tri[0][1] + tri[1][1])/2
    x2,y2 = (tri[1][0] + tri[2][0])/2,(tri[1][1] + tri[2][1])/2
    x3,y3 = (tri[2][0] + tri[0][0])/2,(tri[2][1] + tri[0][1])/2
    return [(x1,y1),(x2,y2),(x3,y3)]

def make_fractal(ax,bef, ite):
    if ite == 0: return 0
    p1,p2,p3 = bef[0],bef[1],bef[2]
    tri = make_tri(bef)
    ax.add_patch(pat.Polygon(xy = tri,fc = "white", ec = "black"))
    make_fractal(ax,points(p1,tri), ite-1)
    make_fractal(ax,points(p2,tri), ite-1)
    make_fractal(ax,points(p3,tri), ite-1)

tri_rate = IntSlider(min=0, max=7, step=1, value=0)
@interact(tri_rate=tri_rate)
def plot(tri_rate):
    tri = [(0.2, 0.2), (0.8, 0.2), (0.5, 0.8)]
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(1,1,1)
    ax.add_patch(pat.Polygon(xy = tri,fc = "white", ec = "black"))
    make_fractal(ax,tri,tri_rate)
    plt.show()
download-1


コッホ曲線
スウェーデンの数学者
ヘルゲ・フォン・コッホが考案した
線分を3等分し分割した2点を
頂点とする正三角形の作図を
無限に繰り返すことによって
得られる図形です。

PythonのPILのImage.drawを用いて
座標を結ぶ直線を描くことができるので
これを使っていきます。

最初は直線を引き
それを分割するような点を求め
三角形の頂点を求め・・・

と言うのを繰り返すことで
描画する座標を求めて
その点同士を直線で結んでいきます。

from PIL import Image,ImageDraw
from math import cos,sin,radians
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

def divine(num,tp):
     if num>0:
        (Ax,Ay),ind=tp[0],1
        for (Bx,By) in tp[1:]:
            ABx,ABy= Bx - Ax , By - Ay
            ACx,ACy= ABx / 3,ABy / 3
            Cx,Cy = Ax + ACx,Ay + ACy
            Dx,Dy = Ax + ACx*2,Ay + ACy*2
            Ex,Ey = cos(r)*ACx - sin(r)*ACy + Cx,sin(r)*ACx + cos(r)*ACy + Cy
            tp.insert(ind,  (Cx,Cy))
            tp.insert(ind+1,(Ex,Ey))
            tp.insert(ind+2,(Dx,Dy))
            ind,Ax,Ay=ind+4,Bx,By
        divine(num-1,tp)

def calc_cell(n,t):
    divine(n,t)
    return t

n1_rate = IntSlider(min=0, max=5, step=1, value=0)
n2_rate = IntSlider(min=3, max=12, step=1, value=3)
@interact(n1_rate=n1_rate,n2_rate=n2_rate)
def plot(n1_rate,n2_rate):
    white,black=(255,255,255),(0,0,0)
    w,h = 2000,2000
    img=Image.new('RGB',(h,w), white)
    draw=ImageDraw.Draw(img)
    x,y,r=[],[],w//2
    for _x in np.linspace(-180,180,361):
        if _x%(360//n2_rate)==0:
            x.append(int(sin(radians(_x))*r*3//4)+r)
            y.append(int(cos(radians(_x))*r*3//4)+r)
    x.append(x[0])
    y.append(y[0])
    poslist = []
    xys = [(xx,yy) for xx,yy in zip(x,y)]
    tmps = [[xys[b],xys[b+1]]for b in range(len(xys)-1)]
    for tmp in tmps:
        poslist += calc_cell(n1_rate,tmp)
    (Ax,Ay)=poslist[0]
    for (Bx,By) in poslist[1:]:
        draw.line((Ax,Ay, Bx,By),fill=black,width=5)
        Ax,Ay = Bx,By
    plt.figure(figsize=(10,10))
    plt.imshow(img)
    plt.show()

download

スライダーで分割数と
多角形の辺の数を定義できます。

より深くしていくと
曲線に近くなりますね。

ドラマ
「危険なビーナス」では
このフラクタル図形と言うのが
物語の鍵になっているようです。

人間が描こうとすると
かなり至難の技ですが
コンピューターであれば
描く点の座標さえ計算できれば
どんな複雑な図形であっても
描くことができます。

ドラマのフラクタル図形は
もっと複雑でしたね。

これを作るのは大変そう・・・
興味がある方は
作ってみてください。

それでは。

このページのトップヘ