C#3.0とLINQでマンデルブロ集合を描いてみた

ただいまLINQ勉強中。無理やり（？）LINQを使ってマンデルブロ集合を描いてみたよ、という話。LINQ面白いよLINQ。

マンデルブロ集合の定義をWikipediaから引用しておく

マンデルブロ集合（まんでるぶろしゅうごう、Mandelbrot set）とは、次の漸化式
$z_{n+1} = z_n^2 + c$
$z_0=0$
で定義される複素数列 $\{z_n\} n \in N$ が n → ∞ の極限で無限大に発散しないという条件を満たす複素数 c 全体が作る集合のことである。
マンデルブロ集合 - Wikipedia

LINQをどう使うのか

LINQは基本的には検索用の構文だ。マンデルブロ集合の描画にどう使うの？って疑問に思うかもしれない。
でも考えてみよう。マンデルブロ集合ってのは複素数全体の部分集合なのだ。そして、検索というのはある集合からある条件に適う部分集合を抽出することを言うのだ。ならば複素数全体からマンデルブロ集合を「検索」するロジックもLINQで書けるに違いない。
つまり、発想としてはこんな感じ。

from c in 複素数全体 where c ∈ マンデルブロ集合 select c;

もちろん、現実には複素数全体から検索してくるのは不可能だ。今回は次の複素平面上の矩形領域から検索することにした。
$D = \left{z = x + iy | -1.5 \le x < 1.5, -1.5 \le y < 1.5 \right}$
従って、検索文はこうなる。

from c in 矩形領域D where c ∈ マンデルブロ集合 select c;

まずは複素数型 Complex を作る

ここは今回の主題じゃないのでさらっと、次のような複素数型を定義したよ、ってだけで説明終わり。

struct Complex
{
    public readonly double X, Y;

    public Complex(double x, double y) { X = x; Y = y; }
    public Complex(double x) : this(x, 0) { }

    public double Abs() { return Math.Sqrt(X * X + Y * Y); }

    public static Complex operator +(Complex z1, Complex z2)
    {
        return new Complex(z1.X + z2.X, z1.Y + z2.Y);
    }

    public static Complex operator *(Complex z1, Complex z2)
    {
        double x = z1.X * z2.X - z1.Y * z2.Y;
        double y = z1.X * z2.Y + z1.Y * z2.X;
        return new Complex(x, y);
    }
}

漸化式を定義する

冒頭のWikipedia引用部分にある漸化式を、複素数のIEnumerableとして定義した。下記のループには終了判定がないところがミソ。これで無限の長さの数列が定義できたことになる。遅延評価を利用すればコンピュータにも「無限」が扱えるんですなぁ。

static IEnumerable<Complex> MandelbrotSequence(Complex c)
{
    for (var z = new Complex(); true; z = z * z + c) yield return z;
}

マンデルブロ集合に含まれるか判定する

ある複素数 c がマンデルブロ集合に含まれるかどうかを、次の関数で判定する。

static bool IsMandelbrot(Complex c)
{
    return MandelbrotSequence(c).Take(100).All(z => z.Abs() < 100.0);
}

この関数では、漸化式を最大100番目まで評価していて、それらの複素数 z の値がすべて複素平面上の半径100の円に収まっているのなら発散しなかった、つまりマンデルブロ集合の点である、と判定している。
これが一行で書けちゃうなんて！超クール！

マンデルブロ集合のBitmapを作成する。

ではいよいよ、次のようなBitmapを生成する関数を定義することにしたい。

static Bitmap GenerateMandelbrotImage(Size size) { ... }

その実装に入る前に、まずはピクセル座標を複素平面にマッピングする関数を用意しておこう。下記の関数を使えば、各ピクセルを複素平面上の矩形領域 D にマッピングさせることが出来る。

static Complex PixelToComplex(int x, int y, Size size)
{
    return new Complex(3.0 * x / size.Width - 1.5, 3.0 * y / size.Height - 1.5);
}

準備が整ったので、GenerateMandelbrotImage() の実装に突入だ。まず、イメージに含まれるピクセル全体は、次のLINQ構文で取得することが出来る。

var seq =
    from x in Enumerable.Range(0, size.Width)
    from y in Enumerable.Range(0, size.Height)
    select new Point(x, y);

ここではマンデルブロ集合に含まれるピクセルだけを取り出したいので、次のように where 文を追加すれば良い。

var mandelbrots =
    from x in Enumerable.Range(0, size.Width)
    from y in Enumerable.Range(0, size.Height)
    where IsMandelbrot(PixelToComplex(x, y, size))
    select new Point(x, y);

あとは、得られたピクセルを黒く塗るだけでOKだ。LINQ万歳！

Bitmap bitmap = new Bitmap(size.Width, size.Height);
foreach (var p in mandelbrots) bitmap.SetPixel(p.X, p.Y, Color.Black);
return bitmap;

フォームに描画する

もう後は、得られたビットマップをフォームに描画するだけだ。簡単簡単。

public partial class Form1 : Form
{
    Bitmap _bitmap;

    public Form1()
    {
        InitializeComponent();
        _bitmap = GenerateMandelbrotImage(this.ClientSize);
    }

    protected override void OnSizeChanged(EventArgs e)
    {
        _bitmap = GenerateMandelbrotImage(this.ClientSize);
        this.Invalidate();
        base.OnSizeChanged(e);
    }

    protected override void OnPaint(PaintEventArgs e)
    {
        e.Graphics.DrawImage(_bitmap, new Point());
        base.OnPaint(e);
    }
 
    ...
}

最後に

LINQには色んな使い道がありそうだ。数値計算にも使えたら面白いな。ぼちぼち模索してみることにしよう。