前回で基本画像PNGファイルを出力できるようになった。
1000x1000のPNG画像にマンデルブロ集合を描くとは、いったいどのような事なのだろうか？
そのためにはまず、マンデルブロ集合とはなにかを理解する必要がある。
上記のWikipediaにあるように以下の漸化式を満たす複素数の集合を一つ一つ色を付けてプロットしていけばマンデルブロ集合が画像として出来上がる。
\begin{align}
&z_0 = 0 \\
&z_{n+1} = z_n^2 + c
\end{align}
この場合のｃが複素数である。
たとえば以下の複素数を
\begin{align}
1+1i
\end{align}
1000x1000のPNG画像を縦軸に虚数、横軸に実数の複素平面としてプロットしてみる。
(便宜的に縦軸(虚数軸)は＋とーの方向を逆にしている）
またPNG画像の中心部が0,0地点、つまり 0+0i の複素数の値がある場所を表している。
\begin{align}
&z_0 = 0 \\
&z_{1} = 0^2 + (1+1i) \\
&z_{2} = (1+1i)^2 + (1+1i) \\
&. \\
&. \\
&. \\
\end{align}
1+1i の複素数の漸化式を３回まで計算してみると上記のようになる、これを1000回繰り返し、値の絶対値が２を超えるまでの回数をHSVのHとして360までの数値に変換し（具体的には360で割った余り）複素数のあった位置にカラーの点としてプロットする。
1000回までに発散しなかった場合は色0としてプロット。
入力された複素数を漸化式で1000回まで行って結果を返す関数を作成する。
for文を使っても良いのだかここではせっかくなのでパターンマッチと再帰を利用する。
// マンデルブロ集合を計算(漸化式)
fn mandelbrpt(c: &Complex<f64>, z: Complex<f64>, n: u32) -> u32 {
let z1 = z * z + *c;
match (c, 2.0, z1, n) {
//1000回で発散しなかった場合は 0 を返す
(_, _, _, 0) => 0,
//発散した場合はその回数を360で割って余りを返す
(_, _, z1, _) if z1.norm() > 2.0 => n % 360,
//それ以外は回数を一つ減らして再帰呼出し
(..) => mandelbrpt(c, z1, n - 1),
}
}
- 再帰回数によってはスタックオーバーフローすると思われるが、リリースビルドの場合末尾再帰は最適化されるため実際は発生しない。
plot関数を修正する
PNG画像の各点に対応する複素数の並びを作成し、mandelbrpt関数をマップして0,0から1000,1000までの100万個のプロットデータを作成する。
完成したコード
use image::{ImageBuffer, Rgb, RgbImage};
use irospace::{converter::*, ColorConverterBuilder, HsvColor, RgbColor};
use itertools::Itertools;
use num::{complex::Complex, integer::div_rem};
//画像のサイズを指定 (1000)
static SIZE: u32 = 1000;
//画像サイズに合わせる複素数の範囲を指定
static SPAN: f64 = 4.0;
//画像中心点の複素数の実数値を設定
static RE_CENTER: f64 = 0.0;
//画像中心点の複素数の虚数値を設定
static IM_CENTER: f64 = 0.0;
//画像サイズに合わせる複素数の左上の複素数をセット
static CPOINT1: Complex<f64> = Complex::new(RE_CENTER - (SPAN / 2.0), IM_CENTER - (SPAN / 2.0));
//画像サイズに合わせる複素数の右下の複素数をセット
static CPOINT2: Complex<f64> = Complex::new(RE_CENTER + (SPAN / 2.0), IM_CENTER + (SPAN / 2.0));
// 隣のピクセルとの複素数の差（X軸、Y軸）
static DX: Complex<f64> = Complex::new((CPOINT2.re - CPOINT1.re) / SIZE as f64, 0.0);
static DY: Complex<f64> = Complex::new(0.0, (CPOINT2.re - CPOINT1.re) / SIZE as f64);
fn main() {
// 1000x1000の画像データを作成
let plot_data: Vec<u32> = plot(SIZE);
// PNG画像を作成
let mut img: ImageBuffer<image::Rgb<u8>, Vec<u8>> = RgbImage::new(SIZE, SIZE);
// HSVからRGBへの変換器を作成
let converter = ColorConverterBuilder::new().from_hsv().to_rgb().build();
// 1000x1000の x,y座標の組み合わせを作成
let grid = (0..SIZE).cartesian_product(0..SIZE);
// x,y座標に色を割り当てる
grid.for_each(|x| {
// 色を取得
let color = plot_data[(x.0 * SIZE as u32 + x.1) as usize];
// HSVのhに色を割り当てる
let color_code = HsvColor::new(color as u16, 255, 255);
// HSVからRGBへ変換
let rgb: RgbColor = converter.convert(&color_code).unwrap();
// 画像に色を割り当てる
img.put_pixel(x.0, x.1, Rgb([rgb.r(), rgb.g(), rgb.b()]));
});
// 画像を保存
img.save("out_image.png").unwrap();
}
// プロットするデータを作成
fn plot(size: u32) -> Vec<u32> {
// 複素数のベクトルを作成
let cvec: Vec<Complex<f64>> = (0..(size * size))
.map(|c| {
let (q, r) = div_rem(c, SIZE);
CPOINT1 + DX.scale(q as f64) + DY.scale(r as f64)
})
.collect();
// 複素数の0を作成
let complex_num0: num::Complex<f64> = Complex::new(0.0, 0.0);
// マンデルブロ集合を計算
cvec.into_iter()
.map(|c: Complex<f64>| mandelbrpt(&c, complex_num0, 1000))
.collect()
}
// マンデルブロ集合を計算(漸化式)
fn mandelbrpt(c: &Complex<f64>, z: Complex<f64>, n: u32) -> u32 {
let z1 = z * z + *c;
match (c, 2.0, z1, n) {
//1000回で発散しなかった場合は 0 を返す
(_, _, _, 0) => 0,
//発散した場合はその回数を360で割って余りを返す
(_, _, z1, _) if z1.norm() > 2.0 => n % 360,
//それ以外は回数を一つ減らして再帰呼出し
(..) => mandelbrpt(c, z1, n - 1),
}
}
実行した画像
以下の定数をお好みで変化させると色々な場所が表示される。
//画像サイズに合わせる複素数の範囲を指定
static SPAN: f64 = 0.00002;
//画像中心点の複素数の実数値を設定
static RE_CENTER: f64 = -0.743643135;
//画像中心点の複素数の虚数値を設定
static IM_CENTER: f64 = 0.131825963;
