読者です 読者をやめる 読者になる 読者になる

gnuplotでフラクタル図形を描こう

gnuplot

gnuplot は, プログラムちっくなお絵かきするのに優れています.
あまりご存じないかもしれませんが, gnuplot では条件分岐, ループ, 関数等のプログラミングをするには十分な機能を備えています.
再帰呼び出しをすることもできます.


再帰呼び出しと言えば, フラクタル図形ですね.
gnuplotフラクタル図形を描くことはとても簡単です.


シェルピンスキーのギャスケット

とても簡単に描けます.

reset

number = 8
sizex = 900
sizey = 900
wid = 0.95
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / 2 * sqrt(3) / sizey) / 2

set terminal jpeg size sizex,sizey
set output "sierpinski.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

sierpinski(n, x, y, dx) \
  = n >= number ?  \
    sprintf("set object polygon from %f,%f to %f,%f to %f,%f to %f,%f fillstyle solid;", x, y, x + dx, y, x + dx / 2, y + dx / 2 * sqrt(3), x, y) : \
    sierpinski(n + 1, x, y, dx / 2). \
    sierpinski(n + 1, x + dx / 2, y, dx / 2). \
    sierpinski(n + 1, x + dx / 4, y + dx / 4 * sqrt(3), dx / 2);

eval(sierpinski(0, 0.0, 0.0, 1.0))

plot -100


再帰的に関数を呼び, sprintf で生成した大量の arrow を文字列結合演算子「.」で結合しまくり, 最後に eval すれば描画できます.
条件分岐には, 三項演算子を使います.

シェルピンスキーのカーペット

こちらは, 直ぐに再帰爆発するから, あまりnumberを増やせない.

reset

number = 5
sizex = 800
sizey = 800
wid = 0.95
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / sizey) / 2

set terminal jpeg size sizex,sizey
set output "carpet.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

carpet(n, x, y, d) \
  = n >= number ?  \
    sprintf("set object rect from %f,%f to %f,%f fc rgb '#000000' fs solid;", x, y, x + d, y + d) : \
    carpet(n + 1, x, y, d / 3) . \
    carpet(n + 1, x + d / 3, y, d / 3) . \
    carpet(n + 1, x + 2 * d / 3, y, d / 3) . \
    carpet(n + 1, x, y + d / 3, d / 3) . \
    carpet(n + 1, x + 2 * d / 3, y + d / 3, d / 3) . \
    carpet(n + 1, x, y + 2 * d / 3, d / 3) . \
    carpet(n + 1, x + d / 3, y + 2 * d / 3, d / 3) . \
    carpet(n + 1, x + 2 * d / 3, y + 2 * d / 3, d / 3);

eval(carpet(0, 0.0, 0.0, 1.0))

plot -100


コッホ曲線

回転行列が計算できたら簡単に描けます.

reset

number = 6
sizex = 1200
sizey = 400
wid = 0.9
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / 6 * sqrt(3) / sizey) / 2

set terminal jpeg size sizex,sizey
set output "koch.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

rotx(x,y,t) = cos(t) * x - sin(t) * y
roty(x,y,t) = sin(t) * x + cos(t) * y
koch(n, x, y, dx, dy) \
  = n >= number ?  \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx, y+dy) : \
    koch(n + 1 \
        , x \
        , y \
        , dx / 3 \
        , dy / 3). \
    koch(n + 1 \
        , x + dx / 3 \
        , y + dy / 3 \
        , rotx(dx / 3, dy / 3, pi / 3) \
        , roty(dx / 3, dy / 3, pi / 3)). \
    koch(n + 1 \
        , x + dx / 3 + rotx(dx / 3, dy / 3, pi / 3) \
        , y + dy / 3 + roty(dx / 3, dy / 3, pi / 3) \
        , rotx(dx / 3, dy / 3, - pi / 3) \
        , roty(dx / 3, dy / 3, - pi / 3)). \
    koch(n + 1 \
        , x + 2 * dx / 3 \
        , y + 2 * dy / 3 \
        , dx / 3 \
        , dy / 3);

eval(koch(0, 0.0, 0.0, 1.0, 0.0))

plot -100


number = 6 の行を変化させたりして遊んでみて下さい.

3つ描いたらコッホ雪片.

高木曲線

難しそうに思える高木曲線も, 意外と簡単なプログラムで描けます.
これも, number を変化させて遊んでみて下さい.

reset

number = 9
sizex = 400
sizey = 400
wid = 0.95
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / 2 * sqrt(3) / sizey) / 2

set terminal jpeg size sizex,sizey
set output "takagi.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

takagi(n, x, y, dx, dy) \
  = n >= number ?  \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x + dx, y + dy) : \
    takagi(n + 1, x, y, dx / 2, dy / 2 + dx / 2). \
    takagi(n + 1, x + dx / 2, y + dy / 2 + dx / 2, dx / 2, dy / 2 - dx / 2);

eval(takagi(0, 0.0, 0.0, 1.0, 0.0))

plot -100


コッホ曲線なんかよりも, 逆に簡単ですね...

ドラゴン曲線

こちらもシンプルなシステム.

reset

number = 12
sizex = 600
sizey = 600
wid = 0.95
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / sizey) / 2

set terminal jpeg size sizex,sizey
set output "dragon.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

hilbert(n, x, y, dx, dy) \
  = n >= number ?  \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x + dx, y + dy) : \
    hilbert(n + 1, x, y, (dx - dy) / 2, (dy + dx) / 2) . \
    hilbert(n + 1, x + dx, y + dy, - (dx + dy) / 2, (dx - dy) / 2);

eval(hilbert(0, 0.2, 0.4, 0.7, 0.0))

plot -100

ヒルベルト曲線

これが, 一番アタマを捻ったし, 時間がかかった.
隣と接続するのが結構難しい.
でも, きちんと動いたプログラムを見なおしたら大したことなかった.
手で図をきちんと描けば分かる.

reset

number = 4
sizex = 600
sizey = 600
wid = 0.95
pad = (1 - wid) / 2 * 100;
padbottom = (1 - sizex * wid / sizey) / 2

set terminal jpeg size sizex,sizey
set output "hilbert.jpeg"
set lmargin pad
set rmargin pad
set tmargin 0
set bmargin 0
set origin 0,padbottom
set border 0
set notitle
set nokey
unset xtics
unset ytics
set xrange [0:1]
set yrange [0:1.0/sizex*sizey/wid]
set style arrow 1 nohead linewidth 1

hilbert(n, x, y, dx, dy, dx2, dy2, dx3, dy3) \
  = n >= number ?  \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - dy) / 4, y + (dy + dx) / 4, x + (3 * dx - dy) / 4, y + (3 * dy + dx) / 4) . \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x + (3 * dx - dy) / 4, y + (3 * dy + dx) / 4, x + (3 * dx - 3 * dy) / 4, y + (3 * dy + 3 * dx) / 4) . \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - 3 * dy) / 4, y + (dy + 3 * dx) / 4, x + (3 * dx - 3 * dy) / 4, y + (3 * dy + 3 * dx) / 4) . \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - 3 * dy) / 4, y + (dy + 3 * dx) / 4, x + (dx - 3 * dy) / 4 + dx2 / 2, y + (dy + 3 * dx) / 4 + dy2 / 2) . \
    sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - dy) / 4, y + (dy + dx) / 4,  x + (dx - dy) / 4 + dx3 / 2, y + (dy + dx) / 4 + dy3 / 2) : \
    hilbert(n + 1, x + dx / 2, y + dy / 2, - dy / 2, dx / 2, dx3 / 2, dy3 / 2, dx / 2, dy / 2). \
    hilbert(n + 1, x + dx / 2, y + dy / 2, dx / 2, dy / 2, - dy / 2, dx / 2, 0.0, 0.0). \
    hilbert(n + 1, x + dx / 2 - dy / 2, y + dy / 2 + dx / 2, dx / 2, dy / 2, -dx / 2, -dy / 2, 0.0, 0.0). \
    hilbert(n + 1, x - dy, y + dx, dy / 2, - dx / 2, 0.0, 0.0, dx2 / 2, dy2 / 2);

eval(hilbert(0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0))

plot -100


number = 5にした時はこちら

迷路みたいー♡

まとめ

gnuplot は, お絵かきするのに使える
フラクタル図形はおもしろい