曼德博集
曼德博集其实是一个“没什么用”的发现。
曼德博集（Mandelbrot Set）是一种在复平面上形成独特且复杂图案的点的集合。这个集合是以数学家本华·曼德博（Benoit Mandelbrot）的名字命名的，他在研究复杂结构和混沌理论时发现了这个集合。曼德博集是分形几何的一个经典例子，显示了一个简单的数学公式如何能产生无限复杂和美丽的图案。
曼德博集的定义相对简单。对于每一个复数
c
c
c，我们考虑以下迭代序列：
Z
n
+
1
=
z
n
2
+
c
其中
(
z
0
=
0
)
Z_{n+1} = z_n^2 + c \;\;\\ 其中 \;\; (z_0 = 0)
Zn+1=zn2+c其中(z0=0)
**曼德博集合由那些使得上述序列不趋于无限大的复数
c
c
c 组成**。在复平面上，这些点形成了一种独特的图案，通常以一种美丽且艺术的方式呈现。这个图案的边界非常复杂，包含了无限的细节和自相似的结构。这意味着无论你放大图案的哪一部分，你都会发现越来越精细的结构，这些结构在形式上与整体图案相似。
曼德博集合不仅在数学上有意义，也在艺术和科学中有广泛的应用，尤其是在研究混沌理论和复杂系统时。
具体可以看
- 维基百科-曼德博集
- bilibili - 2000亿倍放大曼德博集
目标功能
最终我们将实现一个命令行工具，它会根据我们输入的参数生成曼德博集图，使用如下：
./mandelbrot <FILE><PIXELS><UPPERLEFT><LOWERRIGHT>
FILE: 曼德博集图生成的图片路经。PIXELS: 图片分辨率，如4x3。UPPERLEFT: 指定在复平面中图片覆盖的左上角，如4.0,3.0。LOWERRIGHT: 制定在复平面中图片覆盖的右下角。
所以我们最终会根据指定的图片范围，截取
PIXELS
分辨率大小的曼德博集图：
基于以上目标，我们拆分成几个问题：
- 如何表示复数？
- 如何解析分辨率和坐标？
- 如何将图上像素映射到复数？
- 如何生成曼德博集图？即如何找到那些符合曼德博集的点，并将其进行着色标注？
- 如何写入图片文件？
- 如何渲染曼德博集？
- 如何解析命令行参数？
- 如何并发写入图片文件？
能学到什么
- 曼德博集是什么？
- Rust 中的复数的原理与应用。
- Rust 泛型初探。
- Rust 中的 Option 和 Result 初探。
- Rust 并发初探。
- Rust 中如何解析命令行参数？
- Rust 如何写入图像文件？
- Rust 如何写测试用例？
- Rust 实用 crate
num、image、crossbeam。
版本
[package]
name = "mandelbrot"
version = "0.1.0"
edition = "2021"
[dependencies]
image = {version = "0.13.0", features = ["default", "png"]}
num = "0.4.1"
crossbeam = "0.8"
rayon = "1.10.0"
完整代码：https://github.com/hedon-rust-road/mandelbrot/blob/master/src/main.rs
编码实现
0. 创建项目
cargo new mandelbrot
cd mandelbrot
1. 复数表示
使用复数，我们需要引入一个 crete：
num
：
cargoadd num
其中定义了一个复数类型
Complex
：
pubstructComplex<T>{/// 复数的实部pub re:T,/// 复数的虚部pub im:T,}
这里
T
是 Rust 中的泛型功能，表示任意类型
T
，确定好这个结构体的
T
的类型后，其中的属性
re
和
im
的类型也就随之确定了。
2. 解析分辨率和坐标
- 分辨率格式为：4000x3000
- 坐标格式为：-1.0,2.0
2.1 解析数对
我们要做的就是，将分辨率拆成 (4000,3000)，将坐标拆为 (-1.0, 2.0)。这里：
- 带解析的元素
s是一个字符串&str。 - 分隔符
separator是一个字符char。 - 返回值是一个元组
(T, T)，其中T这里可以是 u64/f32 等数字，它们都需要能从字符串转化而来，即<T:FromStr>。 - 因为解析可能出错，所以我们使用
Option来承载。
/// 把字符串 `s`（形如 `"400×600"` 或 ``"1.0,0.5"）解析成一个坐标对////// 具体来说，`s` 应该具有<left><sep><right>的格式，其中<sep>是由`separator`/// 参数给出的字符，而<left>和<right>是可以被 `T:from_str` 解析的字符串。/// `separator` 必须是 ASCII 字符////// 如果 `s` 具有正确的格式，就返回 `Some(x,y)`，否则返回 `None`fnparse_pair<T:FromStr>(s:&str, separator:char)->Option<(T,T)>{match s.find(separator){None=>None,Some(index)=>match(T::from_str(&s[..index]),T::from_str(&s[index +1..])){(Ok(l),Ok(r))=>Some((l, r)),
_ =>None,},}}
我们可以写几个测试用例来验证一下这个函数的正确性，这里我们用到
#[test]
和
assert_eq!
：
#[test]fntest_parse_pair(){assert_eq!(parse_pair::<i32>("",','),None);assert_eq!(parse_pair::<i32>("10,",','),None);assert_eq!(parse_pair::<i32>(",10",','),None);assert_eq!(parse_pair::<i32>("10,20",','),Some((10,20)));assert_eq!(parse_pair::<i32>("10,20xy",','),None);assert_eq!(parse_pair::<f64>("0.5x",'x'),None);assert_eq!(parse_pair::<f64>("0.5x1.5",'x'),Some((0.5,1.5)));}
2.2 转为复数
我们需要的参数
upper_left
和
lower_right
都是复平面中的一个点，所以从字符串中将数对解析完毕后，我们将其赋值到复数的实部和虚部，转为复数实例。
// 把一对用逗号隔开的浮点数解析为复数fnparse_complex(s:&str)->Option<Complex<f64>>{matchparse_pair(s,','){Some((re, im))=>Some(Complex{ re, im }),None=>None,}}
3. 将像素点映射成复数
第 2 步我们其实确定了两件事：
- 确定截取曼德博集的哪一部分。
- 要在这个部分中画多少个点。
这一步我们需要把
x
点转为复数，即确定它的横坐标和纵坐标。这部分可能需要发挥一下你的几何数学能力了（🤡🤡🤡）。
/// 给定输出图像重像素的行和列，返回复平面中对应的坐标////// `pixed` 是表示给图片中特定像素的 (column, row) 二元组。/// `upper_left` 参数和 `lower_right` 参数是在复平面中表示指定图像覆盖范围的点。fnpixed_to_point(/*
·--------------------> bounds.0 re
丨
丨
丨
bounds.1 im
*/
bounds:(usize,usize),
pixed:(usize,usize),
upper_left:Complex<f64>,
lower_right:Complex<f64>,)->Complex<f64>{let(width, height)=(
lower_right.re - upper_left.re,// 右-左
upper_left.im - lower_right.im,// 上-下);Complex{
re: upper_left.re + pixed.0asf64* width / bounds.0asf64,
im: upper_left.im - pixed.1asf64* height / bounds.1asf64,}}#[test]fntest_pixed_to_point(){assert_eq!(pixed_to_point((100,200),(25,175),Complex{ re:-1.0, im:1.0},Complex{ re:1.0, im:-1.0}),Complex{
re:-0.5,
im:-0.75,});}
4. 寻找曼德博集点
什么是曼德博集点？看看上面的定义：**曼德博集合由那些使得上述序列不趋于无限大的复数
c
c
c 组成**。
现在我们可以来表示上述的公式
Z
n
+
1
=
z
n
2
+
c
Z_{n+1} = z_n^2 + c
Zn+1=zn2+c 了：
fncomplex_square_add_loop(c:Complex<f64>){letmut z =Complex{ re:0.0, im:0.0};loop{
z = z * z + c
}}
其中我们将泛型结构体
Complex
的
T
确定为
f64
，并使用
loop
关键字进行无限循环。
所以我们的目标是什么？**找到令
z
不会“飞到”无穷远的
c
**。
由于复数
c
c
c 具有实部 re 和虚部 im，因此可以把它们视为笛卡尔平面上某个点的 x 坐标和 y 坐标，如果
c
c
c 在曼德博集中，就在其中用黑色着色，否则就用浅色。因此，对于图像中的每个像素，必须在复平面上相应点位运行前面的循环，看看它是否逃逸到无穷远还是永远绕着原点运行，并相应将其着色。
无限循环肯定是不现实的，我们总要找到退出循环的机会，有 2 个思路：
- 进行有限次数的迭代，这样可以获得该集合的一个不错的近似值，迭代的次数取决了精度的需要；
- 业界已证明，一旦
z离开了以原点为中心的半径 2 的圆，它最终一定会“飞到”无穷远。
所以我们最终确定的函数如下，其中
norm_sqr()
会返回
z
跟复平面原点的距离的平方：
/// 尝试测试 `c` 是否位于曼德博集中，使用最多 `limit` 次迭代来判定////// 如果 `c` 不是集合成员之一，则返回 `Some(i)`，其中 `i` 是 `c` 离开以原点/// 为中心的半径为 2 的圆时所需的迭代次数。如果 `c` 似乎是集群成员之一（确/// 切而言是达到了迭代次数限制但仍然无法证明 `c` 不是成员），则返回 `None`fnescape_time(c:Complex<f64>, limit:usize)->Option<usize>{letmut z =Complex{ re:0.0, im:0.0};for i in0..limit {if z.norm_sqr()>4.0{returnSome(i);}
z = z * z + c
}None}
5. 写入图片文件
我们可以使用
image
这个 crate 来写入图片文件，它支持多种格式图片的读写，并内置了多种颜色色值。
这里我们准备生成 png 图片，且需要对图片进行不同颜色的着色，所以我们引入
default
和
png
这两个 feature。
cargoadd image --features default,png
5.1 创建文件 File::create()
我们可以用标准库中的
File::create(filename)
来创建一个文件，成功的话会返回一个文件句柄：
let output =File::create(filename)?;
5.2 写入图片 PNGEncoder
image
中提供了
PNGEncoder
用于写入 png 图片，它有两个核心方法：
impl<W:Write>PNGEncoder<W>{/// Create a new encoder that writes its output to ```w```pubfnnew(w:W)->PNGEncoder<W>{PNGEncoder{
w: w
}}/// Encodes the image ```image```/// that has dimensions ```width```and ```height```/// and ```ColorType``````c```pubfnencode(self, data:&[u8], width:u32, height:u32, color:ColorType)->io::Result<()>{let(ct, bits)= color.into();letmut encoder =png::Encoder::new(self.w, width, height);
encoder.set(ct).set(bits);letmut writer =try!(encoder.write_header());
writer.write_image_data(data).map_err(|e| e.into())}}
new(w): 传进目标 writer，即我们上面创建的output。encode(): 写入图片信息，这里有几个参数： -width: u32: 图片宽度。-height: u32: 图片高度。-color: ColorType: 颜色类型，可以是 RGB, Gray(8) 等。-data: &[u8]: 像素色值列表，它的长度应该由上面 3 个字段共同决定，如果选取的颜色是 RGB，意味着需要 3 个 u8 才能表示一个像素点的颜色，所以长度为 width * height * 3，如果选取的颜色是 Gray(8)，那么我们用 1 个 u8 就可以表示一个像素点的灰度值，所以长度为 width * height * 1。本文中我们会采用 Gray(8) 来汇总曼德博集的黑白图。
6. 渲染曼德博集
这一步我们需要来确定上述
PNGEncoder::encode()
的 4 个参数：
width: u32: 图片宽度由命令行参数中指定即可。height: u32: 图片高度由命令行参数中指定即可。color: ColorType: 本文我们只绘制黑白图，这里使用ColorType::Gray(8)，它表示图像是一个灰度（单色）图像，每个像素用8位（即1个字节）来表示。在这种格式中，每个像素的灰度值范围是 0 到 255，其中 0 通常表示黑色，255 表示白色，中间值表示不同的灰度。data: &[u8]: 像素色值列表，我们需要确定 width * height 个像素的灰度值。首先我们根据第 3 步将像素点映射成复数 c c c，然后使用第 4 步中的escape_time()函数来判断复数 c c c 是否位于曼德博集中，如果是，则着黑色，即赋值0，如果不是，则看它迭代了多少次才失败，次数越多，则越接近曼德博集，颜色越深，即越靠近 0，所以赋值255-time。
最终我们实现的函数如下：
/// 将曼德博集对应的矩形渲染到像素缓冲区中////// `bounds` 参数会给缓冲区 `pixels` 的宽度和高度，此缓冲区的每个字节都/// 包含一个灰度像素。`upper_left` 和 `lower_right` 参数分别指定了/// 复平面中对应于像素缓冲区左上角和右上角的点。fnrender(
pixels:&mut[u8],
bounds:(usize,usize),
upper_left:Complex<f64>,
lower_right:Complex<f64>,){assert_eq!(pixels.len(), bounds.0* bounds.1);for raw in0..bounds.1{for column in0..bounds.0{let point =pixed_to_point(bounds,(column, raw), upper_left, lower_right);
pixels[raw * bounds.0+ column]=matchescape_time(point,255){None=>0,Some(count)=>255- count asu8,}}}}
7. 解析命令行参数
核心逻辑部分到这里其实就完成了，现在我们要做最后一步，就是解析命令行参数，让程序可以根据我们的要求绘制曼德博集图。
7.1 解析 std::env::args()
在 Rust 中解析命令行参数的一个常用方法是使用
std::env::args
函数，这个函数返回一个迭代器，它包含了命令行上传递给程序的所有参数。对于更复杂的命令行参数解析，可以使用像
clap
或
structopt
这样的第三方库，这些库提供了更高级的功能和更好的错误处理。
下面是一个使用
std::env::args
的基本例子：
usestd::env;fnmain(){let args:Vec<String>=env::args().collect();for arg in args.iter(){println!("{}", arg);}}
7.2 基础版程序
到这里，我们就可以实现完整的基础版程序了。
fnmain(){// 读取参数let args:Vec<String>=env::args().collect();// 参数个数 = 1 + 4，其中第 1 个是应用程序名if args.len()!=5{eprintln!("Usage: {} FILE PIXELS UPPERLEFT LOWERRIGHT", args[0]);eprintln!("Example: {} mandel.png 1000x700 -1.20,0.35 -1,0.20",
args[0]);std::process::exit(1);}// 解析参数let bounds =parse_pair(&args[2],'x').expect("error parsing image dimensions");let upper_left =parse_complex(&args[3]).expect("error parsing upper left corner point");let lower_right =parse_complex(&args[4]).expect("error parsing lower right corner point");letmut pixels =vec![0; bounds.0* bounds.1];// 渲染曼德博集render(&mut pixels, bounds, upper_left, lower_right);// 输出图片write_image(&args[1],&pixels, bounds).expect("error writing PNG file");}
我们在项目根目录下编译一下程序：
cargo build --release
会在 target/release 下生成可执行文件，执行：
./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20
执行后你应该可以看到我们生成的曼德博集图如下：
大概是处于这个位置：
8. 并发渲染
在 macOS 或 linux 系统下，我们可以使用
time
来输出程序的执行时间：
time ./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20
./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20 3.30s user 0.01s system 98% cpu 3.341 total
笔者使用的电脑为 macbook Pro m2 max 芯片 32 G 内存 12 核，可以看到在单核模式下，差不多需要 3~4s 的时间。
几乎所有的现代机器都有多个处理器核心，而当前这个程序只使用了一个。如果可以把此工作分派个机器提供的多个处理器核心，则应该可以更快地画完图像。
为此，我们可以将图像划分成多个部分，每个处理器负责其中的一个部分，并让每个处理器为分派给它的像素着色。为简单起见，可以将其分成一些水平条带，如下图所示：
crossbeam 是 Rust 中的一个并发编程工具箱，它广泛用于提供各种并发和多线程编程的组件。
crossbeam::scope
是 crossbeam 提供的一个非常有用的功能，它允许你安全地创建临时的线程，并确保这些线程在离开作用域之前结束。
这里我们引入 crossbeam：
cargoadd crossbeam
我们将
fn main()
中的：
render(&mut pixels, bounds, upper_left, lower_right);
替换成：
// 使用 8 个线程来并发执行let threads =8;// 计算每个线程负责渲染的高度，向上取整let rows_per_band = bounds.1/ threads +1;{// chunks_mut() 会返回一个迭代器，该迭代器会生成此缓冲区的可变且不可迭代的切片let bands:Vec<&mut[u8]>= pixels.chunks_mut(rows_per_band * bounds.0).collect();// crossbeam::scope 确保所有子线程在作用域结束之前完成，// 这防止了悬垂指针和其他数据竞争问题。crossbeam::scope(|spawner|{// 遍历像素缓冲区的各个条带，// 这里 into_iter() 迭代器会为循环体的每次迭代赋予独占一个条带的所有权，// 确保一次只有一个线程可以写入它。for(i, band)in bands.into_iter().enumerate(){// 确定每个条带的参数let top = rows_per_band * i;let height = band.len()/ bounds.0;let band_bounds =(bounds.0, height);let band_upper_left =pixed_to_point(bounds,(0, top), upper_left, lower_right);let band_lower_right =pixed_to_point(bounds,(bounds.0, top + height), upper_left, lower_right);// 创建一个线程，渲染图像// move 表示这个闭包会接手它所用遍历的所有权，// 所以只有此闭关，即只有此线程可以使用可变切片 band。
spawner.spawn(move|_|{render(band, band_bounds, band_upper_left, band_lower_right);});}}).unwrap();}
再次执行：
time ./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20
./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20 3.57s user 0.01s system 335% cpu 1.067 total
可以看到虽然总共使用的 CPU 时间还是 3~4s，但是整个程序的执行时间只缩短到 1s 左右了。
9. rayon 工作窃取
前面我们使用 8 个工作线程优化了曼德博集的绘制速度，大概是 4 倍的速度提升。其实这还不够快。
问题的根源在于我们没有平均分配工作量。计算图像的一个像素相当于运行一个循环。事实上，图像的浅灰色部分（循环会快速退出的地方）比黑色部分（循环会运行整整 255 次迭代的地方）渲染速度要快得多。因此，虽然我们将整个区域划分成了大小相等的水平条带，但创建了不均等的工作负载，
使用 rayon 很容易解决这个问题。我们可以为输出中的每一行像素启动一个并行任务。这会创建数百个任务，而 rayon 可以在其线程中分配这些任务。有了工作窃取机制，任务的规模是无关紧要的。rayon 会对这些工作进行平衡。
我们先引入
rayon
：
cargoadd rayon
在
main.rs
中引入
rayon
:
userayon::prelude::*;
然后
main
中并发绘制的部分替换为下面的代码：
let bands:Vec<(usize,&mut[u8])>= pixels.chunks_mut(bounds.0).enumerate().collect();
bands.into_par_iter().for_each(|(i, band)|{let top = i;let band_bounds =(bounds.0,1);let band_upper_left =pixed_to_point(bounds,(0, top), upper_left, lower_right);let band_lower_right =pixed_to_point(bounds,(bounds.0, top +1), upper_left, lower_right);render(band, band_bounds, band_upper_left, band_lower_right);});
首先，创建 bands，也就是要传给 rayon 的任务集合。每个任务只是一个元组类型 (usize, &mut [u8])：第一个是计算所需的行号，第二个是要填充的 pixels 切片。我们使用 chunks_mut 方法将图像缓冲区分成一些行，enumerate 则会给每一行添加行号，然后 collect 会将所有数值切片对放入一个向量中。（这里需要一个向量，因为 rayon 只能从数组和向量中创建并行迭代器。）
编译：
cargo build --release
再次执行：
time ./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20
./target/release/mandelbrot mandel.png 4000x3000 -1.20,0.35 -1,0.20 3.96s user 0.01s system 973% cpu 0.408 total
可以看到，这次速度提升更加明显，总共只用了 0.4s 左右的时间。
以上就是实用 Rust 绘制曼德博集实战的全部内容，enjoy，happy coding~
版权归原作者 后端工程师孔乙己 所有， 如有侵权，请联系我们删除。