0x7o 9 июн в 13:12

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

Средний

14 мин

7.5K

Математика*Параллельное программирование*Визуализация данных*Rust*

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

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

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

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

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

Реализация на 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;

Напишем функцию для :

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. При первом запуске необходимо подождать, пока библиотеки скомпилируются.

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

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

Пусть имеется массив 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,
    );
}

Переносим код на 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,
    );
}

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

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

Сглаживание

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

Решить эту проблему можно путём добавления сэмплирования - это наиболее простой, но прожорливый метод. Его суть заключается в том, что мы раз случайным образом изменяем координату целевой точки, просчитываем 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, поэтому код далек от совершенства

Теги:

Хабы:

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

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

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

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

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

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

Сглаживание

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

Публикации

Истории

Работа

Ближайшие события