Mandelbrot with Rust
October 23, 2022
This write-up will explain how an image of the Mandelbrot set can be created with Rust. The image that the program generates is displayed below. It uses a smooth iteration counter that is used to get a smooth coloring. This is one of the first programs that I have written in Rust, and I have kept it pretty simple, meaning that it does not feature multi-threaded support.
Importing the image crate #
To generate images in Rust we can use the image crate. We can add this crate by adding a dependency in the Cargo.toml file.
1[dependencies]
2image = "0.24.4"And by importing the symbol in the top of our program.
1extern crate image;Defining Complex arithmetic #
To create the Mandelbrot set, we will first define a Complex type that will implement Add, Mul, and a method for the magnitude squared.
1#[derive(Clone, Copy)]
2struct Complex {
3 pub a: f32,
4 pub b: f32,
5}We also implement the Clone and Copy traits with the derive attribute.
Addition #
First we define the Add trait. Complex addition is defined as
$$
(a+bi) + (c+di) = a + b + (c+d)i
$$
1impl std::ops::Add for Complex {
2 type Output = Self;
3
4 fn add(self, rhs: Self) -> Self {
5 Complex {
6 a: self.a + rhs.a,
7 b: self.b + rhs.b,
8 }
9 }
10}Multiplication #
And then the Mul trait. Complex multiplication is defined as
$$
(a + bi) \cdot (c + di) = a^2 - b^2 + (ad + bc)i
$$
1impl std::ops::Mul for Complex {
2 type Output = Self;
3
4 fn mul(self, rhs: Self) -> Self {
5 Complex {
6 a: self.a * rhs.a - self.b * rhs.b,
7 b: self.a * rhs.b + self.b * rhs.a ,
8 }
9 }
10}Magnitude #
Finally, we will also add a method to calculate the magnitude. The magnitude is defined as $$ |\ a + bi\ | = \sqrt{a^2 + b^2} $$
However we will leave the square root out of it, since we do not need it. We define this function as arg_sq.
1impl Complex {
2 fn arg_sq(self) -> f32 {
3 self.a * self.a + self.b * self.b
4 }
5}Computational efficiency. While iterating the formula, the escape criteria is $\sqrt{a^2 + b^2} \gt 2$, which can be more efficiently computed by squaring both sides to get rid of the square root, giving $a^2 + b^2 \gt 4$.
The Mandelbrot set #
With that in place, we can define the mandelbrot(x, y) function. The Mandelbrot set is calculated by applying the recursive relationship:
$$
z_{n+1} = z_n^2 + c
$$
where $z$ is a complex number $(a + bi)$.
- If $|\ z\ | \rightarrow \infty$, then the complex number does not belong to the Mandelbrot set.
- If $|\ z\ | \leq 2$, then we say the the complex number does belong to the Mandelbrot set.
We can check if the complex numbers shoot off to infinity by evaluating the recursive relationship, say, for 256 times. If it hasn’t escaped by then, we assume that it is in the Mandelbrot set and we will color that pixel black. The more iterations that we use, the more detailed the image will look. The result of the function is a value $t \in [0, 1]$ that indicates the number of iterations it reached before it escaped to infinity. If the magnitude squared is greater than $4$, then the orbit will escape to infinity. We are checking with $32$ instead of $2^2$ to get rid of some artifacts in the image.
1fn mandelbrot(x: f32, y: f32) -> f32 {
2 let mut z = Complex { a: 0.0, b: 0.0 };
3 let c = Complex { a: x, b: y };
4 let max = 256;
5 let mut i = 0;
6 while i < max && z.arg_sq() < 32.0 {
7 z = z * z + c;
8 i += 1;
9 }
10 return (i as f32 - z.arg_sq().log2().log2()) / (max as f32);
11}It is also calculating the fractional part of $i$, which we use to get a continuous value for $t$. This is called a smooth iteration counter, and is used to create a nice transition between the colors for different discrete values of $i$.
Calculating the pixel color #
We use the value $t$ to loop over a color map function, which is going to return the RGB values of the pixel. The color map function is defined below.
1fn color(t: f32) -> [u8; 3] {
2 let a = (0.5, 0.5, 0.5);
3 let b = (0.5, 0.5, 0.5);
4 let c = (1.0, 1.0, 1.0);
5 let d = (0.0, 0.10, 0.20);
6 let r = b.0 * (6.28318 * (c.0 * t + d.0)).cos() + a.0;
7 let g = b.1 * (6.28318 * (c.1 * t + d.1)).cos() + a.1;
8 let b = b.2 * (6.28318 * (c.2 * t + d.2)).cos() + a.2;
9 [(255.0 * r) as u8, (255.0 * g) as u8, (255.0 * b) as u8]
10}The color function is created by Inigo Quilez, and is used a lot in the Shadertoy community.
Generating the image #
Finally we tie all of this together in the main method, where we are going to generate the image. Let $w$ be the width of the image in pixels, and $h$ the height of the image. Then the image has a domain of $x \in [0, \mathrm{w}]$ and $y \in [0, \mathrm{h}]$ which we need to normalize, because the interesting part of the Mandelbrot set resides in the rectangle $[-2, 2] \times [-2, 2]$.
We will also scale and translate the coordinates to get a better looking image.
1fn main() {
2 let image_width = 1920;
3 let image_height = 1080;
4
5 let mut image_buffer = image::ImageBuffer::new(
6 image_width, image_height);
7
8 for (x, y, pixel) in image_buffer.enumerate_pixels_mut() {
9 let u = x as f32 / image_height as f32;
10 let v = y as f32 / image_height as f32;
11 let t = mandelbrot(2.5 * (u - 0.5) - 1.4, 2.5 * (v - 0.5));
12 *pixel = image::Rgb(color((2.0 * t + 0.5) % 1.0));
13 }
14
15 image_buffer.save("mandelbrot.png").unwrap();
16}If we run this with cargo run, the result will be written to mandelbrot.png, which results in the image that has been included on the top of this document, and below here as well.
Julia sets #
Instead of assigning our pixel position to $c$, we are going to assign it to $z$. This leaves $c$ for another arbitrary value. If we generate the image with these parameters, we get something that is called a Julia set. All of the fractals have a corresponding Julia set, which we get by assigning the position to $z$ instead of $c$.
Note that to get an interesting image, we need to use a good value for $c$. I’ve changed the mandelbrot function to generate the Julia set, which now looks like this:
1fn julia(x: f32, y: f32) -> f32 {
2 let mut z = Complex { a: x, b: y };
3 let c = Complex { a: 0.38, b: 0.28 };
4 let max = 256;
5 let mut i = 0;
6 while i < max && z.arg_sq() < 32.0 {
7 z = z * z + c;
8 i += 1;
9 }
10 return (i as f32 - z.arg_sq().log2().log2()) / (max as f32);
11}We are using the value $(0.38, 0.28)$ for $c$ in this example. This will generate the following image:
Burning ship fractal #
If we take the absolute value of $z$ during the iteration, we can generate the Burning ship fractal. The absolute value of a complex number is defined as $$ \mathrm{abs}({z}) = |\ Re(z)\ | + i|\ Im(z)\ | $$ The absolute bars are already used to indicate the magnitude of $z$, so I have named the function $\mathrm{abs}$ just to be clear.
We will implement this as another trait:
1impl Complex {
2 fn abs(self) -> Self {
3 Complex {
4 a: self.a.abs(),
5 b: self.b.abs()
6 }
7 }
8}In the while loop in the Mandelbrot function, we take the absolute of $z$ at each iteration.
1while i < max && z.arg_sq() < 32.0 {
2 z = z.abs();
3 z = z * z + c;
4 i += 1;
5}If we now generate the image, we will get the following result:
The Burning ship fractal has a pretty cool Julia set, so make sure to check that out as well!
Final source code #
This is the final source code which contains everything that we have implemented above.
1extern crate image;
2
3#[derive(Clone, Copy)]
4struct Complex {
5 pub a: f32,
6 pub b: f32,
7}
8
9impl std::ops::Add for Complex {
10 type Output = Self;
11
12 fn add(self, rhs: Self) -> Self {
13 Complex {
14 a: self.a + rhs.a,
15 b: self.b + rhs.b,
16 }
17 }
18}
19
20impl std::ops::Mul for Complex {
21 type Output = Self;
22
23 fn mul(self, rhs: Self) -> Self {
24 Complex {
25 a: self.a * rhs.a - self.b * rhs.b,
26 b: self.a * rhs.b + self.b * rhs.a ,
27 }
28 }
29}
30
31impl Complex {
32 fn arg_sq(self) -> f32 {
33 self.a * self.a + self.b * self.b
34 }
35}
36
37impl Complex {
38 fn abs(self) -> Self {
39 Complex {
40 a: self.a.abs(),
41 b: self.b.abs()
42 }
43 }
44}
45
46fn mandelbrot(x: f32, y: f32) -> f32 {
47 let mut z = Complex { a: 0.0, b: 0.0 };
48 let c = Complex { a: x, b: y };
49 let max = 256;
50 let mut i = 0;
51 while i < max && z.arg_sq() < 32.0 {
52 z = z * z + c;
53 i += 1;
54 }
55 return (i as f32 - z.arg_sq().log2().log2()) / (max as f32);
56}
57
58fn julia(x: f32, y: f32) -> f32 {
59 let mut z = Complex { a: x, b: y };
60 let c = Complex { a: 0.38, b: 0.28 };
61 let max = 256;
62 let mut i = 0;
63 while i < max && z.arg_sq() < 32.0 {
64 z = z * z + c;
65 i += 1;
66 }
67 return (i as f32 - z.arg_sq().log2().log2()) / (max as f32);
68}
69
70fn burning_ship(x: f32, y: f32) -> f32 {
71 let mut z = Complex { a: 0.0, b: 0.0 };
72 let c = Complex { a: x, b: y };
73 let max = 256;
74 let mut i = 0;
75 while i < max && z.arg_sq() < 32.0 {
76 z = z.abs();
77 z = z * z + c;
78 i += 1;
79 }
80 return (i as f32 - z.arg_sq().log2().log2()) / (max as f32);
81}
82
83fn color(t: f32) -> [u8; 3] {
84 let a = (0.5, 0.5, 0.5);
85 let b = (0.5, 0.5, 0.5);
86 let c = (1.0, 1.0, 1.0);
87 let d = (0.0, 0.10, 0.20);
88 let r = b.0 * (6.28318 * (c.0 * t + d.0)).cos() + a.0;
89 let g = b.1 * (6.28318 * (c.1 * t + d.1)).cos() + a.1;
90 let b = b.2 * (6.28318 * (c.2 * t + d.2)).cos() + a.2;
91 [(255.0 * r) as u8, (255.0 * g) as u8, (255.0 * b) as u8]
92}
93
94fn main() {
95 let image_width = 1920;
96 let image_height = 1080;
97
98 let mut image_buffer = image::ImageBuffer::new(
99 image_width, image_height);
100
101 for (x, y, pixel) in image_buffer.enumerate_pixels_mut() {
102 let u = x as f32 / image_height as f32;
103 let v = y as f32 / image_height as f32;
104 let t = mandelbrot(2.5 * (u - 0.5) - 1.4, 2.5 * (v - 0.5));
105 *pixel = image::Rgb(color((2.0 * t + 0.5) % 1.0));
106 }
107
108 image_buffer.save("mandelbrot.png").unwrap();
109}