ボックスカウント法でフラクタル次元を計算 - でかいチーズをベーグルする

いろいろ調べてたら楽しくなってきたから実装してみた。

やっぱりフラクタルには夢が詰まってると思う。

実装はGithubにあげてみた。
yamaguchiyuto/fractal · GitHub

まとめ

ボックスカウント法でコッホ曲線、シェルピンスキーのギャスケット、直線、点のフラクタル次元を計算した。
コッホ曲線は自分で描いた
結果は理論値と一致した（シェルピンスキーのギャスケットは少しずれたけどそれは画像の問題っぽい）

フラクタル次元

フラクタル次元 - Wikipedia

次元の値が整数にならないすごいやつ。

説明はWikipediaに完全に丸投げ。

ボックスカウント法

画像をボックス（格子状とか）に区切って、フラクタル図形が含まれるボックスの数を数えていろいろやってフラクタル次元を計算する方法。

フラクタル次元をD、ボックスの幅をδ、フラクタル図形を含むボックスの個数をN(δ)とすると、

${ \displaystyle D = - \lim_{\delta \rightarrow 0} \frac{\log N(\delta)}{\log \delta} }$

となる。

実際には極限とか取れないので、いろいろなボックスの幅についてデータ（σ、N(σ)）をとって、log-logプロットをとって、線形回帰をすると、フィットした直線の傾きがボックスカウント法で計算したフラクタル次元になる。

詳しい説明はこのへんが詳しいかも。

Lecture on Fractals: １．フラクタル
 http://www.ae.keio.ac.jp/lab/soc/takeuchi/japla/works/symposium/03/boxcount.pdf

実装は細かいところを気にしなければ割りと簡単。

（直線をフィットして傾きを求めるとかめんどかったからデータをとる部分だけ）

画像を入力するといろいろな大きさのボックスについてデータを取って、（log(δ), log(N(δ))）を出力する。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys
import math
import Image

def count(img, x, y, b):
    """ img: 対象画像、x,y: 画像のサイズ、b: ボックスのサイズ """
    i = 0
    j = 0
    c = 0
    while i < x and j < y:
        flag = False
        for k in range(0, b):
            for l in range(0, b):
                if i+k < x and j+l < y:
                    if img.getpixel((i+k,j+l)) == 0:
                        """ ボックスに図形が含まれていたらカウントして次の図形へ """
                        c += 1
                        flag = True
                        break
            if flag:
                break
        i += b
        if i >= x:
            """ ボックスが右端に達したら左端に戻す """
            i = 0
            j += b

    return c # 図形が含まれていたボックスの数を返す


img = Image.open(sys.argv[1]).convert('1') # 画像を読み込んで二値化
x, y = img.size # 画像サイズを取得

i = x/10
while i > x/1000:
    n = count(img, x, y, i)
    print math.log(i), math.log(n)
    i = i / 2

実験

ちゃんと実装できてるか実験。

ここでは有名ドコロ二つと、フラクタル図形じゃないけど直線と点に対してやってみた。

コッホ曲線

コッホ曲線 - Wikipedia

f:id:yamaguchiyuto:20140428084709p:plain

こんなやつ。

フラクタル次元は大体1.26。

よさ気な画像がなかったので自分で描いた。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys
import Image,ImageDraw
import math

def koch(xs,ys,xt,yt,l):
    """
左から各点は
(xs,ys), (x1,y1), (x2,y2), (x3,y3), (xt,yt)
"""

    pl = l/3.0 # 1/3の長さ

    x1 = (2*xs+xt)/3.0
    y1 = (2*ys+yt)/3.0
    x3 = (2*xt+xs)/3.0
    y3 = (2*yt+ys)/3.0

    a = abs(x3-x1)/pl
    if a > 1.0:
        a = 1.0
    if a < -1.0:
        a = -1.0
    theta = math.acos(a)

    if xt > xs:
        if yt > ys:
            x2 = x3 - math.cos(math.pi/3 + theta) * pl
            y2 = y3 - math.sin(math.pi/3 + theta) * pl
        else:
            x2 = x1 + math.cos(math.pi/3 + theta) * pl
            y2 = y1 - math.sin(math.pi/3 + theta) * pl
    else:
        if yt > ys:
            x2 = x1 - math.cos(math.pi/3 + theta) * pl
            y2 = y1 + math.sin(math.pi/3 + theta) * pl
        else:
            x2 = x3 + math.cos(math.pi/3 + theta) * pl
            y2 = y3 + math.sin(math.pi/3 + theta) * pl

    # 長さが閾値以下になったら描画
    if pl < 10:
        draw.line((xs,ys)+(x1,y1), 0, width=1)
        draw.line((x1,y1)+(x2,y2), 0, width=1)
        draw.line((x2,y2)+(x3,y3), 0, width=1)
        draw.line((x3,y3)+(xt,yt), 0, width=1)
        return

    # 次のレベル
    koch(xs,ys,x1,y1,pl)
    koch(x1,y1,x2,y2,pl)
    koch(x2,y2,x3,y3,pl)
    koch(x3,y3,xt,yt,pl)

if __name__ == '__main__' :
    width = height = int(sys.argv[1])
    im = Image.new("L",(width,height),255)
    draw = ImageDraw.Draw(im)
    xs = 0
    ys = height/2
    xt = xs+width
    yt = ys
    l = math.sqrt((xs-xt)**2+(ys-yt)**2)
    koch(xs,ys,xt,yt,l)
    im.save("koch.png","PNG")

結果

> python koch.py 10000
> python boxcount.py koch.png > koch.box
> cat koch.box
6.90775527898 3.17805383035
6.21460809842 4.06044301055
5.52146091786 4.91265488574
4.8283137373 5.74300318781
4.12713438505 6.6066501862
3.43398720449 7.48493028329
2.7080502011 8.38640090117

これを直線でフィットしたらこんな感じ。
横軸はlog(δ)、縦軸はlog(N(δ))。

f:id:yamaguchiyuto:20140428093907p:plain