perlのImager.plでマンデルブロを描いた

Math

深夜，唐突にフラクタルっぽいものを生成したくなったので勢いで描いてみました．

マンデルブロって何?

カオスな図形として良く取り上げられるので見たことある人は結構いるんじゃないかとは思いますがとりあえず．
Wikipediaによると

とあります．これだけだとなんなのかサッパリなので，↓に具体例を出してみます．


定義としては，複素数平面上における集合(複素数の集合!!)で，上のように数列の極限を用いて定義されるものです．
シンプルな定義なのですが，全体から見るとややこしい雪だるまっぽい図形になっているだとか，集合の境界部分の一部を拡大してみると再び雪だるまっぽい形状が現れるだとかといった面白い性質を持っています．
複素数の集合なので，画像に出力すると本来は白黒(集合に含まれるor Not)になるところなのですが，収束する場合の値や，発散に至った項数などで色を分けることでカラフルでサイケな感じの画像が生まれます．

どうやって計算するのか

上の定義にある通り，各点Cに対して，数列 z_0=0, z_{n+1}={z_n}^2+Cが発散するか否かを判定することになります．
一般に数列の収束を判定するには無限個の項を計算することになってそうそう簡単に出来るもんじゃないのですが，マンデルブロにおいてはある項で |z_k|>2になったとするとそこで発散すると言うことが知られているので，これを利用します．
具体的には， z_{50}とかくらいまで計算をして，途中で |z_k|>2になれば発散，そうでなければ収束と判定します．
当然，もっと先の項まで計算すれば精度の良い画像が得られますが，その分計算に時間がかかってしまいます．

ソース

perlにImager.pmという画像を処理できるライブラリがあるらしいのでそれを使いました．
集合の性質上，すぐ隣の点であっても各点ごとで振る舞いが全然変わってくるので，N*Nサイズの画像ならN*N個の点すべてに対して収束判定をしています．
複素数の計算についてはMath::Complex ライブラリを使うのがいいのかと思ったのですが，いざ使ってみると計算がとても重くてイライラしたので，実部と虚部で分けて普通の実数で計算しています．

#!/usr/local/bin/perl
use warnings;
use strict;
use Imager;

my ($bxl,$bxu,$byl,$byu) = (-1.8,0.6,-1.2,1.2);
my ($wx,$wy) = (450,450);
my $IterN = 100;
my $cband = 0.01;
my $CbN = int(2.0/$cband)+1;

sub convergence {
    my ($cx,$cy) = @_;
    my ($x,$y) = (0.0,0.0);
    for (my $r = 0; $r < $IterN; $r++) {
        my $xa = $x*$x - $y*$y + $cx;
        $y = 2*$x*$y + $cy;
        $x = $xa;
        if ($x*$x+$y*$y > 4) {
            return -$r;
        }
    }
    return int(($x*$x+$y*$y)/$cband);
}


my (@con_color, @div_color);
for (1..$CbN) {
    my $c = 255-int(255*$_/$CbN);
    my $clr = Imager::Color->new(0,$c,$c*0.3);
    push @con_color, $clr;
}
for (1..$IterN) {
    my $c = 255*$_/$IterN;
    my $clr = Imager::Color->new($c*0.6,0,$c);
    push @div_color, $clr;
}


my $img = Imager->new( xsize => $wx, ysize => $wy );
for (my $y=0; $y<$wy; $y++) {
    for (my $x=0; $x<$wx; $x++) {
        my $cx = $bxl + ($bxu-$bxl)*$x/$wx;
        my $cy = $byu - ($byu-$byl)*$y/$wy;
        my $conv = convergence($cx,$cy);
        if ($conv >= 0) { $img->setpixel( x=>$x, y=>$y, color=>$con_color[$conv]); }
        else { $img->setpixel( x=>$x, y=>$y, color=>$div_color[-$conv]); }
    }
    printf "$y\n";
}


$img->write( file=>'mandel_s.png' );

__END__

おまけ

漸化式では z_{n+1}={z_n}^2+Cとなっていますが，これを少し変えて z_{n+1}={z_n}^\alpha+Cに変えるとどうなるのかというのが気になったのでgifにして描画してみました．
α=1.0→11.0くらい?までを連続的に変化させています．
``こぶ''がどんどん分裂したり増えていっているのが見て取れると思います．

参考

マンデルブロ集合の不思議な世界 - マンデルブロ集合の作り方
Mandelbrot set - Wikipedia, the free encyclopedia