Rust 实战丨绘制曼德博集

最新推荐文章于 2024-07-20 08:56:52 发布

后端工程师孔乙己

最新推荐文章于 2024-07-20 08:56:52 发布

阅读量1.4k

点赞数 26

分类专栏： Rust Rust实战文章标签： rust 开发语言后端

本文链接：https://blog.csdn.net/Hedon954/article/details/139676109

版权

Rust 同时被 2 个专栏收录

10 篇文章 0 订阅

订阅专栏

Rust实战

6 篇文章 0 订阅

订阅专栏

曼德博集

曼德博集其实是一个“没什么用”的发现。

曼德博集（Mandelbrot Set）是一种在复平面上形成独特且复杂图案的点的集合。这个集合是以数学家本华·曼德博（Benoit Mandelbrot）的名字命名的，他在研究复杂结构和混沌理论时发现了这个集合。曼德博集是分形几何的一个经典例子，显示了一个简单的数学公式如何能产生无限复杂和美丽的图案。

曼德博集的定义相对简单。对于每一个复数 $c$ ，我们考虑以下迭代序列：
$Z_{n+1} = z_n^2 + c \;\;\\ 其中 \;\; (z_0 = 0)$
曼德博集合由那些使得上述序列不趋于无限大的复数 $c$ 组成。在复平面上，这些点形成了一种独特的图案，通常以一种美丽且艺术的方式呈现。这个图案的边界非常复杂，包含了无限的细节和自相似的结构。这意味着无论你放大图案的哪一部分，你都会发现越来越精细的结构，这些结构在形式上与整体图案相似。

曼德博集合不仅在数学上有意义，也在艺术和科学中有广泛的应用，尤其是在研究混沌理论和复杂系统时。

具体可以看

目标功能

最终我们将实现一个命令行工具，它会根据我们输入的参数生成曼德博集图，使用如下：

./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：

cargo add num

其中定义了一个复数类型 Complex：

pub struct Complex<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`
fn parse_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]
fn test_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 都是复平面中的一个点，所以从字符串中将数对解析完毕后，我们将其赋值到复数的实部和虚部，转为复数实例。

// 把一对用逗号隔开的浮点数解析为复数
fn parse_complex(s: &str) -> Option<Complex<f64>> {
    match parse_pair(s, ',') {
        Some((re, im)) => Some(Complex { re, im }),
        None => None,
    }
}

3. 将像素点映射成复数

第 2 步我们其实确定了两件事：

确定截取曼德博集的哪一部分。
要在这个部分中画多少个点。

目标区域中的像素点

这一步我们需要把 x 点转为复数，即确定它的横坐标和纵坐标。这部分可能需要发挥一下你的几何数学能力了（🤡🤡🤡）。

/// 给定输出图像重像素的行和列，返回复平面中对应的坐标
///
/// `pixed` 是表示给图片中特定像素的 (column, row) 二元组。
/// `upper_left` 参数和 `lower_right` 参数是在复平面中表示指定图像覆盖范围的点。
fn pixed_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.0 as f64 * width / bounds.0 as f64,
        im: upper_left.im - pixed.1 as f64 * height / bounds.1 as f64,
    }
}

#[test]
fn test_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$ 组成。

现在我们可以来表示上述的公式 $Z_{n+1} = z_n^2 + c$ 了：

fn complex_square_add_loop(c: Complex<f64>) {
    let mut z = Complex { re: 0.0, im: 0.0 };
    loop {
        z = z * z + c
    }
}

其中我们将泛型结构体 Complex 的 T 确定为 f64，并使用 loop 关键字进行无限循环。

所以我们的目标是什么？找到令 z 不会“飞到”无穷远的 c。

由于复数 $c$ 具有实部 re 和虚部 im，因此可以把它们视为笛卡尔平面上某个点的 x 坐标和 y 坐标，如果 $c$ 在曼德博集中，就在其中用黑色着色，否则就用浅色。因此，对于图像中的每个像素，必须在复平面上相应点位运行前面的循环，看看它是否逃逸到无穷远还是永远绕着原点运行，并相应将其着色。

无限循环肯定是不现实的，我们总要找到退出循环的机会，有 2 个思路：

进行有限次数的迭代，这样可以获得该集合的一个不错的近似值，迭代的次数取决了精度的需要；
业界已证明，一旦 z 离开了以原点为中心的半径 2 的圆，它最终一定会“飞到”无穷远。

所以我们最终确定的函数如下，其中 norm_sqr() 会返回 z 跟复平面原点的距离的平方：

/// 尝试测试 `c` 是否位于曼德博集中，使用最多 `limit` 次迭代来判定
///
/// 如果 `c` 不是集合成员之一，则返回 `Some(i)`，其中 `i` 是 `c` 离开以原点
/// 为中心的半径为 2 的圆时所需的迭代次数。如果 `c` 似乎是集群成员之一（确
/// 切而言是达到了迭代次数限制但仍然无法证明 `c` 不是成员），则返回 `None`
fn escape_time(c: Complex<f64>, limit: usize) -> Option<usize> {
    let mut z = Complex { re: 0.0, im: 0.0 };
    for i in 0..limit {
        if z.norm_sqr() > 4.0 {
            return Some(i);
        }
        z = z * z + c
    }
    None
}

5. 写入图片文件

我们可以使用 image 这个 crate 来写入图片文件，它支持多种格式图片的读写，并内置了多种颜色色值。

这里我们准备生成 png 图片，且需要对图片进行不同颜色的着色，所以我们引入 default 和 png 这两个 feature。

cargo add 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```
    pub fn new(w: W) -> PNGEncoder<W> {
        PNGEncoder {
            w: w
        }
    }

    /// Encodes the image ```image```
    /// that has dimensions ```width```and ```height```
    /// and ```ColorType``````c```
    pub fn encode(self, data: &[u8], width: u32, height: u32, color: ColorType) -> io::Result<()> {
        let (ct, bits) = color.into();
        let mut encoder = png::Encoder::new(self.w, width, height);
        encoder.set(ct).set(bits);
        let mut 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$ ，然后使用第 4 步中的 escape_time() 函数来判断复数 $c$ 是否位于曼德博集中，如果是，则着黑色，即赋值 0，如果不是，则看它迭代了多少次才失败，次数越多，则越接近曼德博集，颜色越深，即越靠近 0，所以赋值 255-time。

最终我们实现的函数如下：

/// 将曼德博集对应的矩形渲染到像素缓冲区中
///
/// `bounds` 参数会给缓冲区 `pixels` 的宽度和高度，此缓冲区的每个字节都
/// 包含一个灰度像素。`upper_left` 和 `lower_right` 参数分别指定了
/// 复平面中对应于像素缓冲区左上角和右上角的点。
fn render(
    pixels: &mut [u8],
    bounds: (usize, usize),
    upper_left: Complex<f64>,
    lower_right: Complex<f64>,
) {
    assert_eq!(pixels.len(), bounds.0 * bounds.1);

    for raw in 0..bounds.1 {
        for column in 0..bounds.0 {
            let point = pixed_to_point(bounds, (column, raw), upper_left, lower_right);
            pixels[raw * bounds.0 + column] = match escape_time(point, 255) {
                None => 0,
                Some(count) => 255 - count as u8,
            }
        }
    }
}

7. 解析命令行参数

核心逻辑部分到这里其实就完成了，现在我们要做最后一步，就是解析命令行参数，让程序可以根据我们的要求绘制曼德博集图。

7.1 解析 std::env::args()

在 Rust 中解析命令行参数的一个常用方法是使用std::env::args函数，这个函数返回一个迭代器，它包含了命令行上传递给程序的所有参数。对于更复杂的命令行参数解析，可以使用像clap或structopt这样的第三方库，这些库提供了更高级的功能和更好的错误处理。

下面是一个使用std::env::args的基本例子：

use std::env;

fn main() {
    let args: Vec<String> = env::args().collect();

    for arg in args.iter() {
        println!("{}", arg);
    }
}

7.2 基础版程序

到这里，我们就可以实现完整的基础版程序了。

fn main() {
  	// 读取参数
    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");
    let mut 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：

cargo add 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：

cargo add rayon

在 main.rs 中引入 rayon:

use rayon::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~