Generate Mandelbrot Fractal Image in C++

Question

This is a follow-up question for An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. I am trying to implement Mandelbrot Fractal image generator in C++ in this post.

The experimental implementation

• generate_mandelbrot Template Function Implementation

  namespace TinyDIP
  {
      /**
       * @brief Generates a Mandelbrot set fractal image.
       *
       * @tparam ExecutionPolicy The execution policy (e.g., std::execution::seq, std::execution::par).
       * @tparam FloatingPoint The floating-point type for calculations (e.g., double, long double).
       * @param policy The execution policy instance.
       * @param image_width The width of the output image in pixels.
       * @param image_height The height of the output image in pixels.
       * @param x_min The minimum value of the real component (complex plane).
       * @param x_max The maximum value of the real component (complex plane).
       * @param y_min The minimum value of the imaginary component (complex plane).
       * @param y_max The maximum value of the imaginary component (complex plane).
       * @param max_iterations The maximum number of iterations for each point.
       * @return An Image<TinyDIP::RGB> containing the Mandelbrot set.
       */
      template<class ExecutionPolicy, std::floating_point FloatingPoint = double>
      requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
      [[nodiscard]] Image<RGB> generate_mandelbrot(
          ExecutionPolicy&& policy,
          const std::size_t image_width,
          const std::size_t image_height,
          const FloatingPoint x_min,
          const FloatingPoint x_max,
          const FloatingPoint y_min,
          const FloatingPoint y_max,
          const std::size_t max_iterations)
      {
          Image<RGB> image(image_width, image_height);
          const std::size_t total_pixels = image.count();
  
          // Create a view of all pixel indices from 0 to total_pixels - 1.
          auto pixel_indices = std::ranges::views::iota(std::size_t{0}, total_pixels);
  
          // Process all pixels, potentially in parallel.
          std::for_each(
              std::forward<ExecutionPolicy>(policy),
              std::ranges::begin(pixel_indices),
              std::ranges::end(pixel_indices),
              [&](const std::size_t pixel_index)
              {
                  // 1. Map 1D pixel index to 2D image coordinates
                  const std::size_t px = pixel_index % image_width;
                  const std::size_t py = pixel_index / image_width;
  
                  // 2. Map image coordinates to a point in the complex plane
                  const auto x0 = static_cast<FloatingPoint>(px) / (image_width - 1) * (x_max - x_min) + x_min;
                  const auto y0 = static_cast<FloatingPoint>(py) / (image_height - 1) * (y_max - y_min) + y_min;
  
                  // 3. Perform Mandelbrot iteration
                  auto x = static_cast<FloatingPoint>(0);
                  auto y = static_cast<FloatingPoint>(0);
                  std::size_t iteration = 0;
  
                  while (x * x + y * y <= static_cast<FloatingPoint>(4) && iteration < max_iterations)
                  {
                      const auto xtemp = x * x - y * y + x0;
                      y = static_cast<FloatingPoint>(2) * x * y + y0;
                      x = xtemp;
                      iteration++;
                  }
  
                  // 4. Map the result to a color and set the pixel
                  image.set(pixel_index) = map_iterations_to_color(iteration, max_iterations);
              }
          );
  
          return image;
      }
  }
  

• map_iterations_to_color Function Implementation

  namespace TinyDIP
  {
      /**
       * @brief Maps the number of iterations from a fractal calculation to an RGB color.
       * @param iterations The number of iterations completed.
       * @param max_iterations The maximum number of iterations allowed.
       * @return A TinyDIP::RGB struct representing the calculated color.
       */
      [[nodiscard]] constexpr RGB map_iterations_to_color(const std::size_t iterations, const std::size_t max_iterations) noexcept
      {
          if (iterations >= max_iterations)
          {
              return RGB{ 0, 0, 0 }; // Black for points inside the set
          }
  
          // Use a smooth coloring formula based on sine waves for a psychedelic effect.
          // These constants can be tweaked to produce different color palettes.
          constexpr double freq_r = 0.1;
          constexpr double freq_g = 0.15;
          constexpr double freq_b = 0.2;
          constexpr double phase_r = 3.0;
          constexpr double phase_g = 2.5;
          constexpr double phase_b = 1.0;
  
          const double t = static_cast<double>(iterations);
  
          const auto r = static_cast<std::uint8_t>(sin(freq_r * t + phase_r) * 127.5 + 127.5);
          const auto g = static_cast<std::uint8_t>(sin(freq_g * t + phase_g) * 127.5 + 127.5);
          const auto b = static_cast<std::uint8_t>(sin(freq_b * t + phase_b) * 127.5 + 127.5);
  
          return RGB{ r, g, b };
      }
  }
  

The usage of generate_mandelbrot template function:

int main(int argc, char* argv[])
{
    TinyDIP::Timer timer1;
    // --- Fractal Parameters ---
    constexpr std::size_t width = 1200;
    constexpr std::size_t height = 800;
    constexpr std::size_t max_iterations = 255;

    // Region of the complex plane to render
    constexpr double x_min = -2.0;
    constexpr double x_max = 1.0;
    constexpr double y_min = -1.0;
    constexpr double y_max = 1.0;

    // --- Sequential Execution ---
    std::cout << "Generating Mandelbrot set (Sequential)...\n";
    auto start_seq = std::chrono::high_resolution_clock::now();
    
    auto mandelbrot_image_seq = TinyDIP::generate_mandelbrot(
        std::execution::seq,
        width, height,
        x_min, x_max,
        y_min, y_max,
        max_iterations
    );
    
    auto end_seq = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_seq = end_seq - start_seq;
    std::cout << "Sequential generation took: " << diff_seq.count() << " s\n";

    TinyDIP::bmp_write("mandelbrot_sequential", mandelbrot_image_seq);

    // --- Parallel Execution ---
    std::cout << "Generating Mandelbrot set (Parallel)...\n";
    auto start_par = std::chrono::high_resolution_clock::now();
    
    auto mandelbrot_image_par = TinyDIP::generate_mandelbrot(
        std::execution::par,
        width, height,
        x_min, x_max,
        y_min, y_max,
        max_iterations
    );

    auto end_par = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_par = end_par - start_par;
    std::cout << "Parallel generation took: " << diff_par.count() << " s\n";

    TinyDIP::bmp_write("mandelbrot_parallel", mandelbrot_image_par);

    std::cout << "Performance improvement: " << (diff_seq.count() / diff_par.count()) << "x\n";

    return EXIT_SUCCESS;
}

The output image:

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++

• What changes has been made in the code since last question?

  I implemented generate_mandelbrot template function for generating Mandelbrot Fractal image in this post.

• Why a new review is being asked for?

  Please review the implementation of generate_mandelbrot template function, map_iterations_to_color function, and its tests.

Answer 1

This should not be part of TinyDIP

You are writing a digital image processing library. I don't see how generating a Mandelbrot fractal would be part of that. I recommend you remove generate_mandelbrot() from TinyDIP. The same goes for map_iterations_to_color().

If you do add this kind of functionality, you also have to wonder: where does it stop? Do you add Julia fractals as well? What about other fractals? What about other image generating algorithms? Are you going to maintain everything?

Keep it simple

The way you generate pixel indices, then std::for_each() over all indices, then converting those to x and y coordinates is way more complicated than necessary. I'd make it look more like a good old nested for-loop over the x and y coordinates:

Image<RGB> image(image_width, image_height);

for (std::size_t py = 0; py != image_height; ++py) {
    for (std::size_t px = 0; px != image_width; ++px) {
        const auto x0 = (x_max - x_min) * px / image_width + x_min;
        const auto y0 = (y_max - y_min) * py / image_width + y_min;
        …
        image.at(px, py) = …;
    }
}

You could slap a #pragma omp parallel on it. But perhaps it's even better to create a helper function that returns a range of pixels with coordinates:

for (auto& [pixel_value, px, py]: image.pixels_with_coordinates()) {
    const auto x0 = (x_max - x_min) * px / image_width + x_min;
    const auto y0 = (y_max - y_min) * py / image_width + y_min;
    …
    pixel_value = …;
}

That could be used in combination with std::for_each():

Image<RGB> image(image_width, image_height);
auto pixels = image.pixels_with_coordinates();

std::for_each(policy, std::ranges::begin(pixels), std::ranges::end(pixels),
    [&](auto& value, auto px, auto py) {
        const auto x0 = (x_max - x_min) * px / image_width + x_min;
        const auto y0 = (y_max - y_min) * py / image_width + y_min;
        …
        pixel_value = …;
    }
);

Use std::execution::par_unseq

Your implementation of the Mandelbrot fractal doesn't require that pixels are processed in any particular order, so you can use std::execution::par_unseq. I do not expect this to make any noticable difference though.

Answer 2

map_iterations_to_color() could be generalized to mapping a gray-scale value to a color from a color map. Your Mandelbrot function would then first create the color map, and then use the generalized color mapping function to index into the color map.

Computing the color for each pixel is quite expensive: you're using three calls to sin(). Computing a color map of 256 values would be more than enough to create a seamless output (actually, you only need max_iterations values), and would speed up your code significantly.

The generalized color mapping function would then also be useful in other applications.