Рисуем фракталы на Rust и CUDA

от автора

Фракталы — это бесконечные самоподобные фигуры. Они определяются простыми математическими формулами, которые создают удивительную красоту!

В этой статье мы рассмотрим алгоритм визуализации одного из самых известных фракталов на языке Rust с аппаратным ускорением NVIDIA, масштабированием, сглаживанием и многопоточностью.

Множество Мандельброта

Пожалуй, это самый известный фрактал, описанный в 1905 году Пьером Фату. Множество задаётся следующей формулой, где z и C — комплексные числа:

f(z) = z^2 + C

Фрактал можно визуализировать последовательным вычислением функции f(z), начиная с z_0. Доказано, что множество ограничено окружностью с радиусом 2. Если за n итераций точка z осталась внутри окружности (|z| < 2), то считаем, что z принадлежит множеству и закрашиваем её в чёрный цвет.

Реализация на CPU

Для начала давайте напишем программу для визуализации набора без использования NVIDIA CUDA. Создадим пустой проект с помощью cargo init:

$ cargo init mandelbrot_set & cd mandelbrot_set $ tree ├── Cargo.lock ├── Cargo.toml ├── README.md └── src     └── main.rs

В Cargo.toml добавим следующие зависимости:

[package] name = "mandelbrot_set" version = "0.1.0" edition = "2021"  [dependencies] +num = "0.4.3" +image = "0.25.1"
  • image для работы с изображениями

  • num для поддержки комплексных чисел

Импортируем эти библиотеки:

use image::{Rgb, RgbImage}; use num::complex::Complex;

Напишем функцию для f(z) = z^2 + C:

fn num_iters(cx: f64, cy: f64, max_iters: u32) -> u32 {     let mut z = Complex::new(0.0, 0.0);     let c = Complex::new(cx, cy);      for i in 0..=max_iters {         if z.norm() > 2.0 {             return i;         }         z = z * z + c;     }      max_iters }

Обозначим функцию generate_set, которая сохранит визуализацию множества в файл:

fn generate_set(     file_name: String,     max_iters: u32,     x_min: f64,     y_min: f64,     x_max: f64,     y_max: f64,     resolution: u32, ) { }

Создадим и сохраним пустое RGB изображение:

let mut buffer = RgbImage::new(resolution, resolution); buffer.save(&file_name).unwrap();

Рассчитаем количество итераций из num_iters для каждого пикселя изображения:

for x in 0..resolution {     for y in 0..resolution {         let x_percent = x as f64 / resolution as f64;         let y_percent = y as f64 / resolution as f64;         let cx = x_min + (x_max - x_min) * x_percent;         let cy = y_min + (y_max - y_min) * y_percent;         let iters = num_iters(cx, cy, max_iters);     } }

Если предел итераций достигнут, мы оставляем пиксель по умолчанию черным. В противном случае закрашиваем его белым цветом:

let pixel = buffer.get_pixel_mut(x, y); if iters < max_iters {     *pixel = Rgb([255, 255, 255]); }
Полный код функции generate_set
fn generate_set(     file_name: String,     max_iters: u32,     x_min: f64,     y_min: f64,     x_max: f64,     y_max: f64,     resolution: u32, ) {     let mut buffer = RgbImage::new(resolution, resolution);      for x in 0..resolution {         for y in 0..resolution {             let x_percent = x as f64 / resolution as f64;             let y_percent = y as f64 / resolution as f64;             let cx = x_min + (x_max - x_min) * x_percent;             let cy = y_min + (y_max - y_min) * y_percent;             let iters = num_iters(cx, cy, max_iters);             let pixel = buffer.get_pixel_mut(x, y);              if iters < max_iters {                 *pixel = Rgb([255, 255, 255]);             }         }     }      buffer.save(&file_name).unwrap(); }

Вызовем функцию из main:

fn main() {     let resolution = 1024;;     let max_iters = 100;      generate_set(         "fractal.png".to_string(),         max_iters,         -1.5,         -1.0,         0.5,         1.0,         resolution,     ); }
Полный код main.rs
use image::{Rgb, RgbImage}; use num::complex::Complex;   fn generate_set(     file_name: String,     max_iters: u32,     x_min: f64,     y_min: f64,     x_max: f64,     y_max: f64,     resolution: u32, ) {     let mut buffer = RgbImage::new(resolution, resolution);      for x in 0..resolution {         for y in 0..resolution {             let x_percent = x as f64 / resolution as f64;             let y_percent = y as f64 / resolution as f64;             let cx = x_min + (x_max - x_min) * x_percent;             let cy = y_min + (y_max - y_min) * y_percent;             let iters = num_iters(cx, cy, max_iters);             let pixel = buffer.get_pixel_mut(x, y);              if iters < max_iters {                 *pixel = Rgb([255, 255, 255]);             }         }     }      buffer.save(&file_name).unwrap(); }  fn num_iters(cx: f64, cy: f64, max_iters: u32) -> u32 {     let mut z = Complex::new(0.0, 0.0);     let c = Complex::new(cx, cy);      for i in 0..=max_iters {         if z.norm() > 2.0 {             return i;         }         z = z * z + c;     }      max_iters }  fn main() {     let resolution = 1024;;     let max_iters = 100;      generate_set(         "fractal.png".to_string(),         max_iters,         -1.5,         -1.0,         0.5,         1.0,         resolution,     ); }

Запустим программу — cargo run --release. При первом запуске необходимо подождать, пока библиотеки скомпилируются.

fractal.png

fractal.png

В результате получается такой фрактал. Строго математически, это правильная визуализация — точка либо принадлежит множеству (черный цвет), либо нет (белый цвет), но никто не мешает нам добавить немного цвета! Один из способов раскрасить фрактал — рассмотреть, как быстро становится понятно, что точка не принадлежит множеству.

Добавим цвета

Пусть имеется массив gradient с max_iters элементов, в котором будут заключены различные цвета. Если iters < max_iters, то мы берём gradient[iters] и окрашиваем пиксель в этот цвет. В противном случае оставляем точку черной:

fn generate_set(     ...     colors: Vec<&str>,     ... ) {   ...   let gradient = get_gradient(colors, max_iters);   ...   for x in 0..resolution {       for y in 0..resolution {           ...           let pixel = buffer.get_pixel_mut(x, y);           let color = gradient.get(iters as usize).unwrap_or(&[0, 0, 0]);           *pixel = Rgb(*color);       }   } }

Здесь реализованы функции для создания линейного градиента из нескольких цветов HEX:

hex2rgb, lerp_color и get_gradient
fn hex2rgb(hex: &str) -> Result<Vec<u8>, String> {     let hex = hex.trim_start_matches('#');     if hex.len() != 6 {         return Err("Invalid HEX color length".to_string());     }      let r = u8::from_str_radix(&hex[0..2], 16).map_err(|_| "Invalid HEX color")?;     let g = u8::from_str_radix(&hex[2..4], 16).map_err(|_| "Invalid HEX color")?;     let b = u8::from_str_radix(&hex[4..6], 16).map_err(|_| "Invalid HEX color")?;      Ok(vec![r, g, b]) }  fn lerp_color(color1: &[u8; 3], color2: &[u8; 3], value: f64) -> [u8; 3] {     [         (color1[0] as f64 + (color2[0] as f64 - color1[0] as f64) * value) as u8,         (color1[1] as f64 + (color2[1] as f64 - color1[1] as f64) * value) as u8,         (color1[2] as f64 + (color2[2] as f64 - color1[2] as f64) * value) as u8,     ] }  fn get_gradient(gradient_colors: Vec<&str>, max_iters: u32) -> Vec<[u8; 3]> {     let mut colors = vec![];     let mut gradient_colors_rgb = vec![];     for color in &gradient_colors {         let rgb = hex2rgb(color).unwrap();         gradient_colors_rgb.push([rgb[0], rgb[1], rgb[2]]);     }      for i in 0..max_iters {         let color_index = (i as usize * (gradient_colors.len() - 1)) / max_iters as usize;         let color_value = (i as f64 * (gradient_colors.len() as f64 - 1.0)) / max_iters as f64;         let value = color_value % 1.0;         colors.push(lerp_color(             &gradient_colors_rgb[color_index],             &gradient_colors_rgb[color_index + 1],             value,         ));     }      colors }

Создавать цветовые палитры можно на этом сайте — https://color.adobe.com/

Полный код main.rs
use image::{Rgb, RgbImage}; use num::complex::Complex;  fn hex2rgb(hex: &str) -> Result<Vec<u8>, String> {     let hex = hex.trim_start_matches('#');     if hex.len() != 6 {         return Err("Invalid HEX color length".to_string());     }      let r = u8::from_str_radix(&hex[0..2], 16).map_err(|_| "Invalid HEX color")?;     let g = u8::from_str_radix(&hex[2..4], 16).map_err(|_| "Invalid HEX color")?;     let b = u8::from_str_radix(&hex[4..6], 16).map_err(|_| "Invalid HEX color")?;      Ok(vec![r, g, b]) }  fn lerp_color(color1: &[u8; 3], color2: &[u8; 3], value: f64) -> [u8; 3] {     [         (color1[0] as f64 + (color2[0] as f64 - color1[0] as f64) * value) as u8,         (color1[1] as f64 + (color2[1] as f64 - color1[1] as f64) * value) as u8,         (color1[2] as f64 + (color2[2] as f64 - color1[2] as f64) * value) as u8,     ] }  fn get_gradient(gradient_colors: Vec<&str>, max_iters: u32) -> Vec<[u8; 3]> {     let mut colors = vec![];     let mut gradient_colors_rgb = vec![];     for color in &gradient_colors {         let rgb = hex2rgb(color).unwrap();         gradient_colors_rgb.push([rgb[0], rgb[1], rgb[2]]);     }      for i in 0..max_iters {         let color_index = (i as usize * (gradient_colors.len() - 1)) / max_iters as usize;         let color_value = (i as f64 * (gradient_colors.len() as f64 - 1.0)) / max_iters as f64;         let value = color_value % 1.0;         colors.push(lerp_color(             &gradient_colors_rgb[color_index],             &gradient_colors_rgb[color_index + 1],             value,         ));     }      colors }  fn generate_set(     file_name: String,     max_iters: u32,     colors: Vec<&str>,     x_min: f64,     y_min: f64,     x_max: f64,     y_max: f64,     resolution: u32, ) {     let mut buffer = RgbImage::new(resolution, resolution);     let gradient = get_gradient(colors, max_iters);      for x in 0..resolution {         for y in 0..resolution {             let x_percent = x as f64 / resolution as f64;             let y_percent = y as f64 / resolution as f64;             let cx = x_min + (x_max - x_min) * x_percent;             let cy = y_min + (y_max - y_min) * y_percent;             let iters = num_iters(cx, cy, max_iters);             let pixel = buffer.get_pixel_mut(x, y);             let color = gradient.get(iters as usize).unwrap_or(&[0, 0, 0]);             *pixel = Rgb(*color);         }     }      buffer.save(&file_name).unwrap(); }  fn num_iters(cx: f64, cy: f64, max_iters: u32) -> u32 {     let mut z = Complex::new(0.0, 0.0);     let c = Complex::new(cx, cy);      for i in 0..=max_iters {         if z.norm() > 2.0 {             return i;         }         z = z * z + c;     }      max_iters }  fn main() {     let resolution = 1024;     let max_iters = 100;      generate_set(         "fractal.png".to_string(),         max_iters,         vec!["#1C448E", "#6F8695", "#CEC288", "#FFE381", "#DBFE87"],         -1.5,         -1.0,         0.5,         1.0,         resolution,     ); }

fractal.png

fractal.png

Переносим код на NVIDIA CUDA

Связать C++ с Rust можно через библиотеки libc и cc (источник). Добавим их в Cargo.toml:

[package] name = "mandelbrot_set" version = "0.1.0" edition = "2021"  [dependencies] num = "0.4.3" image = "0.25.1" +libc = "0.2.155"  +[build-dependencies] +cc = "1.0.98"

Создадим в папке src файлы mandelbrot_gpu.h, mandelbrot.cpp и mandelbrot_gpu.cu. В новых версиях CUDA встроена библиотека thrust/complex.h для нативной поддержки комплексных чисел. Перепишем функцию num_iters на CUDA:

#include <thrust/complex.h>  _device__ unsigned int num_iters(double cx, double cy, unsigned int max_iters) {     thrust::complex<double> z(0.0, 0.0);     thrust::complex<double> c(cx, cy);      for (unsigned int i = 0; i <= max_iters; ++i) {         if (thrust::abs(z) > 2.0) {             return i;         }         z = z * z + c;     }      return max_iters; }

Создадим функцию gpu_mandelbrot, которая в многопоточном режиме просчитает количество итераций для каждого пикселя изображения:

__global__ void mandelbrot_kernel(unsigned int* results, double* cx, double* cy, unsigned int max_iters, int num_points) {     int idx = blockDim.x * blockIdx.x + threadIdx.x;     if (idx < num_points) {         results[idx] = num_iters(cx[idx], cy[idx], max_iters);     } }  unsigned int* gpu_mandelbrot(double* cx, double* cy, size_t num_points, unsigned int max_iters) {     unsigned int *dev_results, *results;     double *dev_cx, *dev_cy;      results = new unsigned int[num_points];      cudaMalloc((void**)&dev_cx, num_points * sizeof(double));     cudaMalloc((void**)&dev_cy, num_points * sizeof(double));     cudaMalloc((void**)&dev_results, num_points * sizeof(unsigned int));      cudaMemcpy(dev_cx, cx, num_points * sizeof(double), cudaMemcpyHostToDevice);     cudaMemcpy(dev_cy, cy, num_points * sizeof(double), cudaMemcpyHostToDevice);      int threads_per_block = 256;     int blocks_per_grid = (num_points + threads_per_block - 1) / threads_per_block;      mandelbrot_kernel<<<blocks_per_grid, threads_per_block>>>(dev_results, dev_cx, dev_cy, max_iters, num_points);      cudaMemcpy(results, dev_results, num_points * sizeof(unsigned int), cudaMemcpyDeviceToHost);      cudaFree(dev_cx);     cudaFree(dev_cy);     cudaFree(dev_results);      return results; }

В файле mandelbrot.cpp создадим основную функцию для вызова из Rust:

extern "C" {     void calculate_mandelbrot(const double* cx, const double* cy, int num_points, unsigned int max_iters, unsigned int* output) {         double* non_const_cx = const_cast<double*>(cx);         double* non_const_cy = const_cast<double*>(cy);          unsigned int* gpu_results = gpu_mandelbrot(non_const_cx, non_const_cy, num_points, max_iters);          for (int i = 0; i < num_points; ++i) {             output[i] = gpu_results[i];         }          delete[] gpu_results;     } }
Полный код

mandelbrot_gpu.h

#ifndef MANDELBROT_GPU_H #define MANDELBROT_GPU_H  unsigned int* gpu_mandelbrot(double* cx, double* cy, size_t num_points, unsigned int max_iters);  #endif // MANDELBROT_GPU_H

mandelbrot_gpu.cu

#include <thrust/complex.h> #include <cuda_runtime.h> #include "mandelbrot_gpu.h"  __device__ unsigned int num_iters(double cx, double cy, unsigned int max_iters) {     thrust::complex<double> z(0.0, 0.0);     thrust::complex<double> c(cx, cy);      for (unsigned int i = 0; i <= max_iters; ++i) {         if (thrust::abs(z) > 2.0) {             return i;         }         z = z * z + c;     }      return max_iters; }  __global__ void mandelbrot_kernel(unsigned int* results, double* cx, double* cy, unsigned int max_iters, int num_points) {     int idx = blockDim.x * blockIdx.x + threadIdx.x;     if (idx < num_points) {         results[idx] = num_iters(cx[idx], cy[idx], max_iters);     } }  unsigned int* gpu_mandelbrot(double* cx, double* cy, size_t num_points, unsigned int max_iters) {     unsigned int *dev_results, *results;     double *dev_cx, *dev_cy;      results = new unsigned int[num_points];      cudaMalloc((void**)&dev_cx, num_points * sizeof(double));     cudaMalloc((void**)&dev_cy, num_points * sizeof(double));     cudaMalloc((void**)&dev_results, num_points * sizeof(unsigned int));      cudaMemcpy(dev_cx, cx, num_points * sizeof(double), cudaMemcpyHostToDevice);     cudaMemcpy(dev_cy, cy, num_points * sizeof(double), cudaMemcpyHostToDevice);      int threads_per_block = 256;     int blocks_per_grid = (num_points + threads_per_block - 1) / threads_per_block;      mandelbrot_kernel<<<blocks_per_grid, threads_per_block>>>(dev_results, dev_cx, dev_cy, max_iters, num_points);      cudaMemcpy(results, dev_results, num_points * sizeof(unsigned int), cudaMemcpyDeviceToHost);      cudaFree(dev_cx);     cudaFree(dev_cy);     cudaFree(dev_results);      return results; }

mandelbrot.cpp

#include <iostream> #include "mandelbrot_gpu.h"  extern "C" {     void calculate_mandelbrot(const double* cx, const double* cy, int num_points, unsigned int max_iters, unsigned int* output) {         double* non_const_cx = const_cast<double*>(cx);         double* non_const_cy = const_cast<double*>(cy);          unsigned int* gpu_results = gpu_mandelbrot(non_const_cx, non_const_cy, num_points, max_iters);          for (int i = 0; i < num_points; ++i) {             output[i] = gpu_results[i];         }          delete[] gpu_results;     } }

Создадим build.rs — файл, который соберёт и скомпилирует CUDA‑часть:

use cc;  fn main() {     cc::Build::new()         .cuda(true)         .cpp(true)         .flag("-cudart=shared")         .flag("-ccbin=gcc")         .files(&["./src/mandelbrot.cpp", "./src/mandelbrot_gpu.cu"])         .compile("mandelbrot.a"); }

Убедитесь, что вы имеете установленный компилятор CUDA. Выполним cargo build для сборки всего проекта:

$ cargo build     Compiling mandelbrot_set v0.1.0 (/home/danil/home/danil/RustroverProjects/mandelbrot_set) warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used warning: mandelbrot_set@0.1.0: nvcc warning : incompatible redefinition for option 'compiler-bindir', the last value of this option was used     Finished `dev` profile [unoptimized + debuginfo] target(s) in 3.87s

Если вы хотите использовать другой компилятор или получаете ошибки, замените gcc в build.rs на другое значение. Например:

...         .flag("-ccbin=clang") ...

Определим функцию calculate_mandelbrot в main.rs:

... use libc::{c_double, c_int, c_uint};  extern "C" {     fn calculate_mandelbrot(         cx: *mut c_double,         cy: *mut c_double,         num_points: c_int,         max_iters: c_uint,         output: *mut c_uint,     ); } ...

Изменим generate_set для работы с C++:

fn generate_set(     ... ) {     ...     let mut h_cx = vec![];     let mut h_cy = vec![];     let mut output = vec![0; (resolution * resolution) as usize];      for x in 0..resolution {         for y in 0..resolution {             ...             h_cx.push(cx);             h_cy.push(cy);         }     }      unsafe {         calculate_mandelbrot(             h_cx.as_mut_ptr(),             h_cy.as_mut_ptr(),             (resolution * resolution) as c_int,             max_iters,             output.as_mut_ptr(),         );     }      for (x, row) in output.chunks(resolution as usize).enumerate() {         for (y, column) in row.iter().enumerate() {             let pixel = buffer.get_pixel_mut(x as u32, y as u32);             let color = gradient.get(*column as usize).unwrap_or(&[0, 0, 0]);             *pixel = Rgb(*color);         }     }      buffer.save(&file_name).unwrap(); }
Полный код main.rs
use image::{Rgb, RgbImage}; use num::complex::Complex; use libc::{c_double, c_int, c_uint};  extern "C" {     fn calculate_mandelbrot(         cx: *mut c_double,         cy: *mut c_double,         num_points: c_int,         max_iters: c_uint,         output: *mut c_uint,     ); }  fn hex2rgb(hex: &str) -> Result<Vec<u8>, String> {     let hex = hex.trim_start_matches('#');     if hex.len() != 6 {         return Err("Invalid HEX color length".to_string());     }      let r = u8::from_str_radix(&hex[0..2], 16).map_err(|_| "Invalid HEX color")?;     let g = u8::from_str_radix(&hex[2..4], 16).map_err(|_| "Invalid HEX color")?;     let b = u8::from_str_radix(&hex[4..6], 16).map_err(|_| "Invalid HEX color")?;      Ok(vec![r, g, b]) }  fn lerp_color(color1: &[u8; 3], color2: &[u8; 3], value: f64) -> [u8; 3] {     [         (color1[0] as f64 + (color2[0] as f64 - color1[0] as f64) * value) as u8,         (color1[1] as f64 + (color2[1] as f64 - color1[1] as f64) * value) as u8,         (color1[2] as f64 + (color2[2] as f64 - color1[2] as f64) * value) as u8,     ] }  fn get_gradient(gradient_colors: Vec<&str>, max_iters: u32) -> Vec<[u8; 3]> {     let mut colors = vec![];     let mut gradient_colors_rgb = vec![];     for color in &gradient_colors {         let rgb = hex2rgb(color).unwrap();         gradient_colors_rgb.push([rgb[0], rgb[1], rgb[2]]);     }      for i in 0..max_iters {         let color_index = (i as usize * (gradient_colors.len() - 1)) / max_iters as usize;         let color_value = (i as f64 * (gradient_colors.len() as f64 - 1.0)) / max_iters as f64;         let value = color_value % 1.0;         colors.push(lerp_color(             &gradient_colors_rgb[color_index],             &gradient_colors_rgb[color_index + 1],             value,         ));     }      colors }  fn generate_set(     file_name: String,     max_iters: u32,     colors: Vec<&str>,     x_min: f64,     y_min: f64,     x_max: f64,     y_max: f64,     resolution: u32, ) {     let mut buffer = RgbImage::new(resolution, resolution);     let gradient = get_gradient(colors, max_iters);     let mut h_cx = vec![];     let mut h_cy = vec![];     let mut output = vec![0; (resolution * resolution) as usize];      for x in 0..resolution {         for y in 0..resolution {             let x_percent = x as f64 / resolution as f64;             let y_percent = y as f64 / resolution as f64;             let cx = x_min + (x_max - x_min) * x_percent;             let cy = y_min + (y_max - y_min) * y_percent;             h_cx.push(cx);             h_cy.push(cy);         }     }      unsafe {         calculate_mandelbrot(             h_cx.as_mut_ptr(),             h_cy.as_mut_ptr(),             (resolution * resolution) as c_int,             max_iters,             output.as_mut_ptr(),         );     }      for (x, row) in output.chunks(resolution as usize).enumerate() {         for (y, column) in row.iter().enumerate() {             let pixel = buffer.get_pixel_mut(x as u32, y as u32);             let color = gradient.get(*column as usize).unwrap_or(&[0, 0, 0]);             *pixel = Rgb(*color);         }     }      buffer.save(&file_name).unwrap(); }  fn num_iters(cx: f64, cy: f64, max_iters: u32) -> u32 {     let mut z = Complex::new(0.0, 0.0);     let c = Complex::new(cx, cy);      for i in 0..=max_iters {         if z.norm() > 2.0 {             return i;         }         z = z * z + c;     }      max_iters }  fn main() {     let resolution = 1024;     let max_iters = 100;      generate_set(         "fractal.png".to_string(),         max_iters,         vec!["#1C448E", "#6F8695", "#CEC288", "#FFE381", "#DBFE87"],         -1.5,         -1.0,         0.5,         1.0,         resolution,     ); }

Масштабирование

Добавим переменную scale, которая будет отвечать за приближения к определённой точке фрактала:

fn main() {     let resolution = 1024;     let target_x = -0.6582034218739634;     let target_y = 0.44967917993930356;     let max_iters = 5000;     let scale = 500_000.0;      let x_min = target_x - (1.0 / scale);     let x_max = target_x + (1.0 / scale);     let y_min = target_y - (1.0 / scale);     let y_max = target_y + (1.0 / scale);      generate_set(         "fractal.png".to_string(),         ...         x_min,         y_min,         x_max,         y_max,         resolution,     ); }
fractal.png

fractal.png

К сожалению, из-за ограниченной точности float64 качественный масштаб может быть выполнен до 10^{15} раз:

let resolution = 256; ... let scale = 10f64.powf(15.0); ...
fractal.png

fractal.png

Сглаживание

Как видно из последней визуализации, наше изображение имеет большое количество «случайных точек». Это не является ошибкой, но создаёт неприятный шум.

Решить эту проблему можно путём добавления сэмплирования — это наиболее простой, но прожорливый метод. Его суть заключается в том, что мы n раз случайным образом изменяем координату целевой точки, просчитываем iters, берём среднее значение итераций и красим целевую точку. Вот, как это выглядит.

Для генерации случайных значений, добавим соответствующую библиотеку:

[package] name = "mandelbrot_set" version = "0.1.0" edition = "2021"  [dependencies] num = "0.4.3" image = "0.25.1" libc = "0.2.155" +rand = "0.9.0-alpha.1"  [build-dependencies] cc = "1.0.98"

Изменим generate_set:

use rand::Rng;   fn generate_set( ...   samples: u32, ... ) {   let mut rng = rand::thread_rng();   ...   for _ in 0..samples {       ...       let x_percent = (x as f64 + rng.gen::<f64>()) / resolution as f64;       let y_percent = (y as f64 + rng.gen::<f64>()) / resolution as f64;       ...   }   unsafe {     calculate_mandelbrot(       ...       (resolution * resolution * samples) as c_int,       ...     );   }   for (x, row) in output.chunks((resolution * samples) as usize).enumerate() {       for (y, column) in row.chunks(samples as usize).enumerate() {           let mut sum = 0;           for iteration in column {               sum += *iteration as usize;           }           let pixel = buffer.get_pixel_mut(x as u32, y as u32);           let color = gradient.get(sum / column.len()).unwrap_or(&[0, 0, 0]);           *pixel = Rgb(*color);       }   } }
fn main() {     ...     let samples = 16;     ...     generate_set(         ...         samples     ); }

Выглядит гораздо лучше!

Финальные штрихи

Воспользуемся экспоненциальным ростом масштаба, чтобы анимировать масштабирование фрактала:

let scale = 10.0_f64.powf((i as f64 / frames as f64) * args.max_scale.ilog10() as f64);

Добавим библиотеку rayon для рендера в многопоточном режиме и clap для поддержки аргументов командной строки. Итоговая версия проекта лежит на GitHub:

$ git clone https://github.com/0x7o/mandelbrot_set $ cd mandelbrot_set $ cargo build --release
$ ./target/release/mandelbrot_set \     --resolution 1024 \     --colors "#00A3BC, #8B00BD, #81BD00, #BD5400" \     --x -0.6582034218739634 \     --y 0.44967917993930356 \     --iters 3000 \     --max-scale 1000000000000000 \     --fps 24 \     --seconds 60 \     --n-samples 4 \     --output ./output
$ ffmpeg -framerate 24 -i output/frame_%09d.png -c:v libx264 -pix_fmt yuv420p -crf 18 -y video.mp4 

Спасибо за прочтение! Я только начинаю изучать Rust, поэтому код далек от совершенства


ссылка на оригинал статьи https://habr.com/ru/articles/820489/