マンデルブロ集合をコードに実装する方法

42 Tokyoの「fract-ol」という課題で、フラクタル図形の描写に取り組んだので、

前日はマンデルブロ集合の数式について説明したので、本日は実際のコードに落とし込む方法について見ていきたい。

コードに実装するには？

実際にコードに落とし込む際に複素数は使えないので、マンデルブロ集合を複素数を使わずに書き直すことが必要になる。

書き直しを行うには、Zn を点 (xn, yn) に、C を点 (a, b) にそれぞれ置き代える。

この置き換えを$${Z_{n + 1}  = Z_{n}^{2}  + C}$$の式に当てはめると、

$$
x_{n+1}+y_{n+1}=\left( x_{n}+y_{n}\right) ^{2}+a+b
$$

と表せる。

2乗の箇所の計算をすると、以下のように展開される。

$$
x_{n+1}+y_{n+1}=x_{n}^2+2x_{n}y_{n}-y_{n}^2+a+b
$$


ここで注意しなければいけないことは

$$
-y_{n}^2
$$

とマイナスになっていることだ。

これは、yは虚数iを含んでおり、虚数iの特徴は「2乗すると-1」となることだ。なので、

$$
-y_{n}^2
$$

はマイナスになり、完全なる実数となった。

ここで、実数と虚数部分を分けると、以下の2つの式に分けることができる。

$$
①   x_{n+1}=x_{n}^2-y_{n}^2+a
$$


$$
②   y_{n+1}=2x_{n}y_{n}+b
$$


横軸（ｘ軸）の座標を求めるには上の式を、縦軸（ｙ軸）の座標を求めるのには下の式を使用することになる。

実際のコード

次に実際に使用したコードを見ていく。

前回、マンデルブロ集合は以下の漸化式（前に出した数字を利用して次の値を出す）で定義され、nを大きくしていった時に、「数が大きくならずに安定する」と述べたが、安定する時の条件は、「2乗した際に、4以下」であれば安定しているとみなす。

$$
\begin{cases}
Z_{0} =\ 0\\
Z_{n\ +\ 1} \ =\ Z_{n}^{2} \ +\ C\
\end{cases}
$$

//この関数の中で、Znが安定するのか、発散するのかを確認し、結果によって色付けを行っている。
unsigned int	ft_pick_color(t_data *data)
{
	int		i;
	int		color;
	double	tmp_r;

	i = 0;
	
	//data->z_rは、本文中の式でいうとxn+1にあたる。rはreal number、つまり実数を表している。
	//data->z_iは、本文中の式でいうとyn+1にあたる。iはimaginary number、つまり虚数を表している。
	//このwhile文ではまず、実数部分と虚数部分のそれぞれを2乗したものが4以内、かつ
	//data->max-itは何回計算を行うかを設定しており、初期値では100に設定してあるので、ここでは100回while分が回ることになる。
	while (pow(data->z_r, 2.0) + pow(data->z_i, 2.0) <= 4 && i < data->max_it)
	{
		//以下が、上記の①の式に当たる。tmp_rというテンポラリーの変数に収めている理由は、ここで求めているはdata->z_rだが、data->z_rは次の式でも使用するので、内容を変更しないため。
		//data->c_rはCの実数部分に当たり、①ではaに当たる。
		tmp_r = pow(data->z_r, 2.0) - pow(data->z_i, 2.0) + data->c_r;
		
		//以下が、上記の②の式に当たる。data->c_iはCの虚数部分、②の式だとbに当たる
		data->z_i = 2 * data->z_r * data->z_i + data->c_i;
		data->z_r = tmp_r;
		i++;
	}
	
	//以下は、計算回数が上限に達した時、白（RGBで255,255,255）を塗るようにしている。
	if (i == data->max_it)
		color = ft_rgb_to_hex(255, 255, 255);
	
	//以下では、発散した場合の色つけを行っている。
	//色のつけ方はネット上でたくさんでているが、ここでは
	//計算回数（i）と実数部分（data->z_r）の大きさによって、色の塗り分けを行う関数を呼び出している。
	else
		color = ft_gradation(data, (double)i
				/ (double)data->max_it, pow(data->z_r, 2.0));
	return (color);
}

int	ft_draw_mandelbrot(t_data *data)
{
	int	x;
	int	y;

	//ここでは、描画するwindow上で1ピクセルだけ進んだ際の、複素平面乗で移動する差分を求めている。
	//data->delta_rは実数部分の差分を求めている。
	//data->max_rは複素平面の実数部分の最大値で、data->min_rは複素平面の実数部分の最小値。
	//横幅（WIDTH）は描画するwindowの横軸のピクセル数。
	data->delta_r = (data->max_r - data->min_r) / (double)WIDTH;
	//data->delta_iは虚数部分の差分を求めている。
	//data->max_iは複素平面の虚数部分の最大値で、data->min_iは複素平面の虚数部分の最小値。
	//縦幅（HEIGHT）は描画するwindowの縦軸のピクセル数。
	data->delta_i = (data->max_i - data->min_i) / (double)HEIGHT;
	y = 0;
	while (y < HEIGHT)
	{
		x = 0;
		while (x < WIDTH)
		{
			//ここでは描画するwindow上での1ピクセルずつ移動し、色を決めておくことになる。
			//data->c_rは、Cの実数部分。①のaにあたる。
			//data->c_iは、Cの虚数部分。①のbにあたる。
			//data->z_rとdata->z_iはZnの初期値だが、両方ともOとなる。
			data->c_r = data->min_r + x * data->delta_r;
			data->c_i = data->min_i + y * data->delta_i;
			data->z_r = 0;
			data->z_i = 0;

			//ここでは紹介しないが、ft_my_mlx_pixel_putは実際に色を塗る関数。
			//ft_pick_color関数は上に紹介してある通り、漸化式が安定するか発散するかによって色を選ぶ関数となる。
			ft_my_mlx_pixel_put(data, x, y, ft_pick_color(data));
			x++;
		}
		y++;
	}
	return (0);
}