Generate Julia Fractal Image in C++

Question

This is a follow-up question for Generate Mandelbrot Fractal Image in C++ and An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. Besides Mandelbrot Fractal image, I am trying to implement a Julia Fractal image generator function in C++ as an example of using TinyDIP library.

The experimental implementation

• generate_fractal_color_map Function Implementation

  As Cris Luengo's suggestion, I implemented generate_fractal_color_map function to pre-compute a color map of 256 values.

  /**
      * @brief Generates a pre-computed color map for fractal rendering.
      * @return An std::array of 256 RGB colors.
      */
  [[nodiscard]] constexpr std::array<TinyDIP::RGB, 256> generate_fractal_color_map() noexcept
  {
      std::array<TinyDIP::RGB, 256> color_map{};
      for (std::size_t i = 0; i < 256; ++i)
      {
          // Use a smooth coloring formula based on sine waves for a psychedelic effect.
          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>(i);
  
          const auto r = static_cast<std::uint8_t>(std::sin(freq_r * t + phase_r) * 127.5 + 127.5);
          const auto g = static_cast<std::uint8_t>(std::sin(freq_g * t + phase_g) * 127.5 + 127.5);
          const auto b = static_cast<std::uint8_t>(std::sin(freq_b * t + phase_b) * 127.5 + 127.5);
  
          color_map[i] = TinyDIP::RGB{ r, g, b };
      }
      return color_map;
  }
  

• map_iterations_to_color Function Implementation

  /**
      * @brief Maps the number of iterations to an RGB color using a pre-computed color map.
      * @param iterations The number of iterations completed.
      * @param max_iterations The maximum number of iterations allowed.
      * @param color_map A pre-computed table of colors.
      * @return A TinyDIP::RGB struct representing the calculated color.
      */
  [[nodiscard]] constexpr TinyDIP::RGB map_iterations_to_color(const std::size_t iterations, const std::size_t max_iterations, const std::array<TinyDIP::RGB, 256>& color_map) noexcept
  {
      if (iterations >= max_iterations)
      {
          return TinyDIP::RGB{ 0, 0, 0 }; // Black for points inside the set
      }
  
      // Look up the color from the pre-computed map.
      // The modulo wraps the iteration count to the size of the map.
      return color_map[iterations % 256];
  }
  

• generate_julia Template Function Implementation

  /**
      * @brief Generates a Julia set fractal image for a given constant 'c'.
      *
      * @tparam ExecutionPolicy The execution policy (e.g., std::execution::seq, std::execution::par).
      * @tparam FloatingPoint The floating-point type for calculations (e.g., float, double, long double).
      * @param policy The execution policy instance.
      * @param image_width The width of the output image.
      * @param image_height The height of the output image.
      * @param x_min The minimum value of the real component.
      * @param x_max The maximum value of the real component.
      * @param y_min The minimum value of the imaginary component.
      * @param y_max The maximum value of the imaginary component.
      * @param c The complex constant 'c'.
      * @param max_iterations The maximum number of iterations.
      * @return An Image<TinyDIP::RGB> containing the Julia set.
      */
  template<class ExecutionPolicy, std::floating_point FloatingPoint = double>
  requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
  [[nodiscard]] TinyDIP::Image<TinyDIP::RGB> generate_julia(
      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::complex<FloatingPoint> c,
      const std::size_t max_iterations)
  {
      TinyDIP::Image<TinyDIP::RGB> image(image_width, image_height);
      const FloatingPoint width_float = image_width > 1 ? static_cast<FloatingPoint>(image_width - 1) : 1.0;
      const FloatingPoint height_float = image_height > 1 ? static_cast<FloatingPoint>(image_height - 1) : 1.0;
      const FloatingPoint x_range = x_max - x_min;
      const FloatingPoint y_range = y_max - y_min;
  
      auto proxy = image.pixels_with_coordinates();
  
      static const auto color_map = generate_fractal_color_map();
  
      std::for_each(
          std::forward<ExecutionPolicy>(policy),
          std::ranges::begin(proxy),
          std::ranges::end(proxy),
          [&](auto&& pixel_tuple)
          {
              auto& [pixel_value, px, py] = pixel_tuple;
  
              auto z_real = x_range * static_cast<FloatingPoint>(px) / width_float + x_min;
              auto z_imag = y_range * static_cast<FloatingPoint>(py) / height_float + y_min;
  
              std::size_t iteration = 0;
              while (z_real * z_real + z_imag * z_imag <= static_cast<FloatingPoint>(4) && iteration < max_iterations)
              {
                  const auto z_real_temp = z_real * z_real - z_imag * z_imag + c.real();
                  z_imag = static_cast<FloatingPoint>(2) * z_real * z_imag + c.imag();
                  z_real = z_real_temp;
                  iteration++;
              }
  
              pixel_value = map_iterations_to_color(iteration, max_iterations, color_map);
          }
      );
  
      return image;
  }
  

• The Image class implementation with pixels_with_coordinates member function:

  template <typename ElementT>
  class Image
  {
  public:
      Image() = default;
  
      template<std::same_as<std::size_t>... Sizes>
      Image(Sizes... sizes): size{sizes...}, image_data((1 * ... * sizes))
      {}
  
      template<std::same_as<int>... Sizes>
      Image(Sizes... sizes)
      {
          size.reserve(sizeof...(sizes));
          (size.emplace_back(sizes), ...);
          image_data.resize(
              std::reduce(
                  std::ranges::cbegin(size),
                  std::ranges::cend(size),
                  std::size_t{1},
                  std::multiplies<>()
                  )
          );
      }
  
      //  Image constructor
      template<std::ranges::input_range Sizes>
      requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
      Image(const Sizes& sizes)
      {
          if (sizes.empty())
          {
              throw std::runtime_error("Image size vector is empty!");
          }
          size.resize(sizes.size());
          std::transform(std::ranges::cbegin(sizes), std::ranges::cend(sizes), std::ranges::begin(size), [&](auto&& element) { return element; });
          image_data.resize(
              std::reduce(
                  std::ranges::cbegin(sizes),
                  std::ranges::cend(sizes),
                  std::size_t{1},
                  std::multiplies<>()
                  ));
      }
  
      //  Image constructor
      #ifdef __cpp_lib_containers_ranges
          template<std::ranges::input_range Range,
                   std::same_as<std::size_t>... Sizes>
          Image(Range&& input, Sizes... sizes):
              size{sizes...}, image_data(std::from_range, std::forward<Range>(input))
          {
              if (image_data.size() != (1 * ... * sizes)) {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
          }
      #else
          template<std::ranges::input_range Range,
                   std::same_as<std::size_t>... Sizes>
          Image(const Range&& input, Sizes... sizes):
              size{sizes...}, image_data(input.begin(), input.end())
          {
              if (image_data.size() != (1 * ... * sizes)) {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
          }
      #endif
  
      //  Image constructor
      #ifdef __cpp_lib_containers_ranges
          template<std::ranges::input_range Range, std::same_as<std::size_t>... Sizes>
          Image(const Range& input, Sizes... sizes) :
              size{ sizes... }
          {
              if (input.empty())
              {
                  throw std::runtime_error("Input vector is empty!");
              }
              image_data = std::vector(std::from_range, input);
              if (image_data.size() != (1 * ... * sizes)) {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
          }
      #else
          template<std::ranges::input_range Range, std::same_as<std::size_t>... Sizes>
          Image(const Range& input, Sizes... sizes):
              size{sizes...}
          {
              if (input.empty())
              {
                  throw std::runtime_error("Input vector is empty!");
              }
              image_data = std::vector(input.begin(), input.end());
              if (image_data.size() != (1 * ... * sizes)) {
                  throw std::runtime_error("Image data input and the given size are mismatched!");
              }
          }
      #endif
  
      //  Image constructor
      template<std::ranges::input_range Range, std::ranges::input_range Sizes>
      requires(std::same_as<std::ranges::range_value_t<Sizes>, std::size_t>)
      Image(const Range& input, const Sizes& sizes)
      {
          if (input.empty())
          {
              throw std::runtime_error("Input vector is empty!");
          }
          size.resize(sizes.size());
          std::transform(std::ranges::cbegin(sizes), std::ranges::cend(sizes), std::ranges::begin(size), [&](auto&& element) { return element; });
          image_data = std::vector(std::ranges::cbegin(input), std::ranges::cend(input));
          auto count = std::reduce(std::ranges::cbegin(sizes), std::ranges::cend(sizes), 1, std::multiplies());
          if (image_data.size() != count) {
              throw std::runtime_error("Image data input and the given size are mismatched!");
          }
      }
  
      Image(const std::vector<std::vector<ElementT>>& input)
      {
          if (input.empty())
          {
              throw std::runtime_error("Input vector is empty!");
          }
          size.reserve(2);
          size.emplace_back(input[0].size());
          size.emplace_back(input.size());
          for (auto& rows : input)
          {
              image_data.insert(image_data.end(), std::ranges::begin(rows), std::ranges::end(rows));    //  flatten
          }
          return;
      }
  
      //  at template function implementation
      template<std::same_as<std::size_t>... Args>
      constexpr ElementT& at(const Args... indexInput)
      {
          return const_cast<ElementT&>(static_cast<const Image &>(*this).at(indexInput...));
      }
  
      //  at template function implementation
      //  Reference: https://codereview.stackexchange.com/a/288736/231235
      template<std::same_as<std::size_t>... Args>
      constexpr ElementT const& at(const Args... indexInput) const
      {
          checkBoundary(indexInput...);
          return at_without_boundary_check(indexInput...);
      }
  
      //  at template function implementation
      template<std::same_as<int>... Args>
      constexpr ElementT& at(const Args... indexInput)
      {
          return at(static_cast<std::size_t>(indexInput)...);
      }
  
      //  at template function implementation
      template<std::same_as<int>... Args>
      constexpr ElementT const& at(const Args... indexInput) const
      {
          return at(static_cast<std::size_t>(indexInput)...);
      }
  
      //  at template function implementation (std::ranges::input_range case)
      template<std::ranges::input_range Indices>
      requires(std::same_as<std::ranges::range_value_t<Indices>, std::size_t> ||
               std::same_as<std::ranges::range_value_t<Indices>, int>)
      constexpr ElementT const& at(const Indices indexInput) const
      {
          for (std::size_t i = 0; i < indexInput.size(); ++i)
          {
              if (indexInput[i] > size[i])
              {
                  throw std::out_of_range("Given index out of range!");
              }
          }
          return at_without_boundary_check(indexInput);
      }
  
      //  at template function implementation (std::ranges::input_range case)
      template<std::ranges::input_range Indices>
      requires(std::same_as<std::ranges::range_value_t<Indices>, std::size_t> ||
               std::same_as<std::ranges::range_value_t<Indices>, int>)
      constexpr ElementT& at(const Indices indexInput)
      {
          return const_cast<ElementT&>(static_cast<const Image&>(*this).at(indexInput));
      }
  
      //  at_without_boundary_check template function implementation
      template<std::same_as<std::size_t>... Args>
      constexpr ElementT& at_without_boundary_check(const Args... indexInput)
      {
          return const_cast<ElementT&>(static_cast<const Image &>(*this).at_without_boundary_check(indexInput...));
      }
  
      //  at_without_boundary_check template function implementation
      template<std::same_as<std::size_t>... Args>
      constexpr ElementT const& at_without_boundary_check(const Args... indexInput) const
      {
          constexpr std::size_t n = sizeof...(Args);
          if (n != size.size()) {
              throw std::runtime_error("Dimensionality mismatched!");
          }
          std::size_t index = calculateIndex(indexInput...);
          return image_data[index];
      }
  
      //  at_without_boundary_check template function implementation
      template<std::same_as<int>... Args>
      constexpr ElementT& at_without_boundary_check(const Args... indexInput)
      {
          return at_without_boundary_check(static_cast<std::size_t>(indexInput)...);
      }
  
      //  at_without_boundary_check template function implementation
      template<std::same_as<int>... Args>
      constexpr ElementT const& at_without_boundary_check(const Args... indexInput) const
      {
          return at_without_boundary_check(static_cast<std::size_t>(indexInput)...);
      }
  
      //  at_without_boundary_check template function implementation (std::ranges::input_range case)
      template<std::ranges::input_range Indices>
      requires(std::same_as<std::ranges::range_value_t<Indices>, std::size_t> ||
               std::same_as<std::ranges::range_value_t<Indices>, int>)
      constexpr ElementT const& at_without_boundary_check(const Indices indexInput) const
      {
          std::vector<std::size_t> index_size_t(indexInput.size());
          std::transform(
              std::ranges::cbegin(indexInput),
              std::ranges::cend(indexInput),
              std::ranges::begin(index_size_t),
              [&](auto&& element) {return static_cast<std::size_t>(element); }        //  casting to std::size_t
          );
          return image_data[calculateIndex(index_size_t)];
      }
  
      //  at_without_boundary_check template function implementation (std::ranges::input_range case)
      template<std::ranges::input_range Indices>
      requires(std::same_as<std::ranges::range_value_t<Indices>, std::size_t> ||
               std::same_as<std::ranges::range_value_t<Indices>, int>)
      constexpr ElementT& at_without_boundary_check(const Indices indexInput)
      {
          return const_cast<ElementT&>(static_cast<const Image&>(*this).at_without_boundary_check(indexInput));
      }
  
      //  get function implementation
      constexpr ElementT get(std::size_t index) const noexcept
      {
          return image_data[index];
      }
  
      //  set function implementation
      constexpr ElementT& set(const std::size_t index)
      {
          if (index >= count())
          {
              std::cout << "index = " << index << ", count = " << count() << '\n';
              throw std::out_of_range("Given index out of range!");
          }
          return image_data[index];
      }
  
      //  set template function implementation
      template<class TupleT>
      requires(is_tuple<TupleT>::value and
               check_tuple_element_type<std::size_t, TupleT>::value)
      constexpr bool set(const TupleT location, const ElementT draw_value)
      {
          if (checkBoundaryTuple(location))
          {
              image_data[tuple_location_to_index(location)] = draw_value;
              return true;
          }
          return false;
      }
  
      //  cast template function implementation
      template<typename TargetT>
      constexpr Image<TargetT> cast()
      {
          std::vector<TargetT> output_data;
          output_data.resize(image_data.size());
          std::transform(
              std::ranges::cbegin(image_data),
              std::ranges::cend(image_data),
              std::ranges::begin(output_data),
              [&](const auto& input){ return static_cast<TargetT>(input); }
              );
          Image<TargetT> output(output_data, size);
          return output;
      }
  
      constexpr std::size_t count() const noexcept
      {
          return std::reduce(std::ranges::cbegin(size), std::ranges::cend(size), std::size_t{ 1 }, std::multiplies());
      }
  
      //  count member function implementation
      template<class ExPo>
      requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
      constexpr std::size_t count(ExPo&& execution_policy) const
      {
          if (size.empty()) return 0;
          return std::reduce(std::forward<ExPo>(execution_policy), std::ranges::cbegin(size), std::ranges::cend(size), std::size_t{ 1 }, std::multiplies());
      }
  
      constexpr std::size_t getDimensionality() const noexcept
      {
          return size.size();
      }
  
      constexpr std::size_t getWidth() const noexcept
      {
          return size[0];
      }
  
      constexpr std::size_t getHeight() const noexcept
      {
          return (getDimensionality() > 1) ? size[1] : 0;
      }
  
      //  getSize function implementation
      constexpr auto getSize() const noexcept
      {
          return size;
      }
  
      //  getSize function implementation
      constexpr auto getSize(std::size_t index) const noexcept
      {
          return size[index];
      }
  
      //  getStride function implementation
      constexpr std::size_t getStride(std::size_t index) const noexcept
      {
          if(index == 0)
          {
              return std::size_t{1};
          }
          std::size_t output = std::size_t{1};
          for(std::size_t i = 0; i < index; ++i)
          {
              output *= size[i];
          }
          return output;
      }
  
      std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data
  
      /**
       * print function implementation
       * @brief Prints the image content to an output stream.
       * This function is generic and supports printing N-dimensional images.
       * @param separator The separator to use between elements.
       * @param os The output stream to write to.
       */
      void print(std::string_view separator = "\t", std::ostream& os = std::cout) const
      {
          if (getDimensionality() == 0)
          {
              return;
          }
  
          // A lambda to define how a single element should be printed.
          auto element_printer = [&](const ElementT& value) {
              if constexpr (is_MultiChannel<ElementT>::value)
              {
                   os << "( ";
                   // Assumes .channels is an array-like member of the multi-channel type
                  for (std::size_t i = 0; i < std::size(value.channels); ++i) {
                      os << +value.channels[i] << (i == std::size(value.channels) - 1 ? "" : " ");
                  }
                  os << ") ";
              }
              else if constexpr (is_streamable<ElementT> && !std::is_fundamental_v<ElementT>)
              {
                  os << value;
              }
              else
              {
                  // Use unary '+' to ensure char types are printed as numbers
                  os << +value;
              }
          };
  
          std::vector<std::size_t> indices(getDimensionality(), 0);
          print_recursive_helper(indices, getDimensionality() - 1, element_printer, separator, os);
      }
  
      Image<ElementT>& setAllValue(const ElementT& input)
      {
          std::fill(std::ranges::begin(image_data), std::ranges::end(image_data), input);
          return *this;
      }
  
      //  setAllValue template function implementation (with Execution Policy)
      template<class ExecutionPolicy>
      requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)
      Image<ElementT>& setAllValue(ExecutionPolicy&& execution_policy, const ElementT& input)
      {
          std::fill(std::forward<ExecutionPolicy>(execution_policy), std::ranges::begin(image_data), std::ranges::end(image_data), input);
          return *this;
      }
  
      /**
       * @brief Sets all elements, automatically choosing between sequential and parallel execution.
       * @param The auto_exec tag to invoke this specific overload.
       * @param input The value to set all elements to.
       * @return A reference to the modified image.
       */
      Image<ElementT>& setAllValue(const TinyDIP::auto_execution_policy&, const ElementT& input)
      {
          // This threshold is a heuristic. It represents the number of *bytes* to process
          // before parallel overhead is considered worthwhile. This value often requires
          // empirical tuning for best results on target hardware.
          // Let's set it to ~4MB as a reasonable starting point.
          constexpr std::size_t PARALLEL_FILL_COST_THRESHOLD = 4 * 1024 * 1024;
  
          const std::size_t cost = count() * sizeof(ElementT);
  
          if (cost > PARALLEL_FILL_COST_THRESHOLD)
          {
              // The job is large enough, go parallel!
              std::cout << "[Auto Dispatcher]: Workload is large (" << cost << " bytes). Choosing PARALLEL execution.\n";
              setAllValue(std::execution::par, input);
          }
          else
          {
              // The job is small, stick to sequential to avoid overhead.
              std::cout << "[Auto Dispatcher]: Workload is small (" << cost << " bytes). Choosing SEQUENTIAL execution.\n";
              setAllValue(input); // Call the base sequential version
          }
  
          return *this;
      }
  
      friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
      {
          const std::string separator = "\t";
          rhs.print(separator, os);
          return os;
      }
  
      Image<ElementT>& operator+=(const Image<ElementT>& rhs)
      {
          check_size_same(rhs, *this);
          std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                 std::ranges::begin(image_data), std::plus<>{});
          return *this;
      }
  
      Image<ElementT>& operator-=(const Image<ElementT>& rhs)
      {
          check_size_same(rhs, *this);
          std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                 std::ranges::begin(image_data), std::minus<>{});
          return *this;
      }
  
      Image<ElementT>& operator*=(const Image<ElementT>& rhs)
      {
          check_size_same(rhs, *this);
          std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                 std::ranges::begin(image_data), std::multiplies<>{});
          return *this;
      }
  
      Image<ElementT>& operator/=(const Image<ElementT>& rhs)
      {
          check_size_same(rhs, *this);
          std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                 std::ranges::begin(image_data), std::divides<>{});
          return *this;
      }
  
      friend bool operator==(Image<ElementT> const&, Image<ElementT> const&) = default;
  
      friend bool operator!=(Image<ElementT> const&, Image<ElementT> const&) = default;
  
      friend Image<ElementT> operator+(Image<ElementT> input1, const Image<ElementT>& input2)
      {
          return input1 += input2;
      }
  
      friend Image<ElementT> operator-(Image<ElementT> input1, const Image<ElementT>& input2)
      {
          return input1 -= input2;
      }
  
      friend Image<ElementT> operator*(Image<ElementT> input1, ElementT input2)
      {
          return multiplies(input1, input2);
      }
  
      friend Image<ElementT> operator*(ElementT input1, Image<ElementT> input2)
      {
          return multiplies(input2, input1);
      }
  
      // Nested class for the custom iterator
      class pixel_iterator
      {
      public:
          using iterator_category = std::forward_iterator_tag;
          using value_type = std::tuple<ElementT&, std::size_t, std::size_t>;
          using difference_type = std::ptrdiff_t;
          using pointer = value_type*;
          using reference = value_type;
  
          pixel_iterator() = default;
  
          pixel_iterator(ElementT* ptr, std::size_t x, std::size_t y, std::size_t width)
              : current_ptr(ptr), px(x), py(y), image_width(width)
          {
          }
  
          reference operator*() const
          {
              return {*current_ptr, px, py};
          }
  
          pixel_iterator& operator++()
          {
              ++current_ptr;
              ++px;
              if (px == image_width)
              {
                  px = 0;
                  ++py;
              }
              return *this;
          }
  
          pixel_iterator operator++(int)
          {
              pixel_iterator tmp = *this;
              ++(*this);
              return tmp;
          }
  
          bool operator==(const pixel_iterator& other) const = default;
  
      private:
          ElementT* current_ptr = nullptr;
          std::size_t px = 0;
          std::size_t py = 0;
          std::size_t image_width = 0;
      };
  
      // Nested class for the custom const iterator
      class const_pixel_iterator
      {
      public:
          using iterator_category = std::forward_iterator_tag;
          using value_type = std::tuple<const ElementT&, std::size_t, std::size_t>;
          using difference_type = std::ptrdiff_t;
          using pointer = value_type*;
          using reference = value_type;
  
          const_pixel_iterator() = default;
  
          const_pixel_iterator(const ElementT* ptr, std::size_t x, std::size_t y, std::size_t width)
              : current_ptr(ptr), px(x), py(y), image_width(width)
          {
          }
  
          reference operator*() const
          {
              return {*current_ptr, px, py};
          }
  
          const_pixel_iterator& operator++()
          {
              ++current_ptr;
              ++px;
              if (px == image_width)
              {
                  px = 0;
                  ++py;
              }
              return *this;
          }
  
          const_pixel_iterator operator++(int)
          {
              const_pixel_iterator tmp = *this;
              ++(*this);
              return tmp;
          }
  
          bool operator==(const const_pixel_iterator& other) const = default;
  
      private:
          const ElementT* current_ptr = nullptr;
          std::size_t px = 0;
          std::size_t py = 0;
          std::size_t image_width = 0;
      };
  
      // Nested class for the range proxy object
      class pixel_proxy
      {
      public:
          pixel_proxy(Image<ElementT>& image) : img(image)
          {
              if (img.getDimensionality() != 2)
              {
                  throw std::logic_error("pixels_with_coordinates is only supported for 2D images.");
              }
          }
  
          [[nodiscard]] auto begin()
          {
              return pixel_iterator(img.image_data.data(), 0, 0, img.getWidth());
          }
  
          [[nodiscard]] auto end()
          {
              return pixel_iterator(img.image_data.data() + img.count(), 0, img.getHeight(), img.getWidth());
          }
  
      private:
          Image<ElementT>& img;
      };
  
      // Nested class for the const range proxy object
      class const_pixel_proxy
      {
      public:
          const_pixel_proxy(const Image<ElementT>& image) : img(image)
          {
              if (img.getDimensionality() != 2)
              {
                  throw std::logic_error("pixels_with_coordinates is only supported for 2D images.");
              }
          }
  
          [[nodiscard]] auto begin() const
          {
              return const_pixel_iterator(img.image_data.data(), 0, 0, img.getWidth());
          }
  
          [[nodiscard]] auto end() const
          {
              return const_pixel_iterator(img.image_data.data() + img.count(), 0, img.getHeight(), img.getWidth());
          }
  
      private:
          const Image<ElementT>& img;
      };
  
      /**
       * @brief Returns a range-like object for iterating over pixels with their 2D coordinates.
       * @details This is intended for use in range-based for loops with structured bindings:
       * for (auto& [value, x, y] : image.pixels_with_coordinates()) { ... }
       * @note This function is only valid for 2D images and will throw an exception otherwise.
       * @return A proxy object with begin() and end() methods.
       */
      [[nodiscard]] auto pixels_with_coordinates()
      {
          return pixel_proxy(*this);
      }
  
      /**
       * @brief Returns a const range-like object for iterating over pixels with their 2D coordinates.
       * @details This is intended for use in range-based for loops with structured bindings:
       * for (const auto& [value, x, y] : image.pixels_with_coordinates()) { ... }
       * @note This function is only valid for 2D images and will throw an exception otherwise.
       * @return A proxy object with begin() and end() methods.
       */
      [[nodiscard]] auto pixels_with_coordinates() const
      {
          return const_pixel_proxy(*this);
      }
  
  #ifdef USE_BOOST_SERIALIZATION
  
      void Save(std::string filename)
      {
          const std::string filename_with_extension = filename + ".dat";
          //    Reference: https://stackoverflow.com/questions/523872/how-do-you-serialize-an-object-in-c
          std::ofstream ofs(filename_with_extension, std::ios::binary);
          boost::archive::binary_oarchive ArchiveOut(ofs);
          //    write class instance to archive
          ArchiveOut << *this;
          //    archive and stream closed when destructors are called
          ofs.close();
      }
  
  #endif
  private:
      std::vector<std::size_t> size;
      std::vector<ElementT> image_data;
  
  
      /**
       * print_recursive_helper template function implementation
       * @brief Recursively iterates through dimensions to print the image.
       *
       * @tparam PrintFunc A callable type for printing a single element.
       * @param indices A vector to hold the current N-dimensional index.
       * @param current_dim The current dimension being iterated over.
       * @param print_element The function that prints a single element.
       * @param separator The separator string.
       * @param os The output stream.
       */
      template <std::invocable<const ElementT&> PrintFunc>
      void print_recursive_helper(
          std::vector<std::size_t>& indices,
          std::size_t current_dim,
          const PrintFunc& print_element,
          std::string_view separator,
          std::ostream& os) const
      {
          if (current_dim == 0)                                   // Base case: The innermost dimension
          {
              for (std::size_t i = 0; i < getSize(current_dim); ++i)
              {
                  indices[current_dim] = i;
                  print_element(at_without_boundary_check(indices));
                  os << separator;
              }
          }
          else                                                    // Recursive step: Outer dimensions
          {
              for (std::size_t i = 0; i < getSize(current_dim); ++i)
              {
                  indices[current_dim] = i;
                  print_recursive_helper(indices, current_dim - 1, print_element, separator, os);
                  // After a full "row" of the lower dimension, add a newline.
                  os << '\n';
              }
              // Add an extra newline to separate higher-dimensional "planes".
               os << '\n';
          }
      }
  
      // calculateIndex template function implementation
      template<std::same_as<std::size_t>... Args>
      constexpr auto calculateIndex(const Args... indices) const
      {
          std::size_t index = 0;
          std::size_t stride = 1;
          std::size_t i = 0;
          auto update_index = [&](std::size_t idx) {
              index += idx * stride;
              stride *= size[i++];
              };
          (update_index(indices), ...);
          return index;
      }
  
      //  calculateIndex template function implementation
      template<std::ranges::input_range Indices>
      requires (std::same_as<std::ranges::range_value_t<Indices>, std::size_t>)
      constexpr std::size_t calculateIndex(const Indices& indices) const
      {
          std::size_t index = 0;
          std::size_t stride = 1;
          std::size_t i = 0;
          for (const auto& idx : indices) {
              index += idx * stride;
              stride *= size[i++];
          }
          return index;
      }
  
      template<typename... Args>
      void checkBoundary(const Args... indexInput) const
      {
          constexpr std::size_t n = sizeof...(Args);
          if(n != size.size())
          {
              throw std::runtime_error("Dimensionality mismatched!");
          }
          std::size_t parameter_pack_index = 0;
          auto function = [&](auto index) {
              if (index >= size[parameter_pack_index])
                  throw std::out_of_range("Given index out of range!");
              parameter_pack_index = parameter_pack_index + 1;
          };
  
          (function(indexInput), ...);
      }
  
      //  checkBoundaryTuple template function implementation
      template<class TupleT>
      requires(TinyDIP::is_tuple<TupleT>::value)
      constexpr bool checkBoundaryTuple(const TupleT location)
      {
          constexpr std::size_t n = std::tuple_size<TupleT>{};
          if(n != size.size())
          {
              throw std::runtime_error("Dimensionality mismatched!");
          }
          std::size_t parameter_pack_index = 0;
          auto function = [&](auto index) {
              if (std::cmp_greater_equal(index, size[parameter_pack_index]))
                  return false;
              parameter_pack_index = parameter_pack_index + 1;
              return true;
          };
          return std::apply([&](auto&&... args) { return ((function(args))&& ...);}, location);
      }
  
      //  tuple_location_to_index template function implementation
      template<class TupleT>
      requires(TinyDIP::is_tuple<TupleT>::value)
      constexpr std::size_t tuple_location_to_index(TupleT location)
      {
          std::size_t i = 0;
          std::size_t stride = 1;
          std::size_t position = 0;
          auto update_position = [&](auto index) {
                  position += index * stride;
                  stride *= size[i++];
              };
          std::apply([&](auto&&... args) {((update_position(args)), ...);}, location);
          return position;
      }
  
  #ifdef USE_BOOST_SERIALIZATION
      friend class boost::serialization::access;
      template<class Archive>
      void serialize(Archive& ar, const unsigned int version)
      {
          ar& size;
          ar& image_data;
      }
  #endif
  
  };
  

• is_streamable concept

  //  is_streamable concept
  // Concept to check if a type supports the << operator with std::ostream
  template<typename T>
  concept is_streamable = requires(std::ostream & os, const T & val)
  {
      { os << val } -> std::same_as<std::ostream&>;
  };
  

The usage of generate_julia template function:

void generate_julia_set()
{
    // --- Fractal Parameters ---
    constexpr std::size_t width = 2400;
    constexpr std::size_t height = 1600;
    constexpr std::size_t max_iterations = 255;
    constexpr double x_min = -1.5;
    constexpr double x_max = 1.5;
    constexpr double y_min = -1.0;
    constexpr double y_max = 1.0;
    // An interesting constant for the Julia set
    constexpr std::complex<double> c{ -0.7269 , 0.1889 };

    std::cout << "--- Generating Julia Set (c = " << c << ") ---\n";

    // --- Sequential Execution ---
    std::cout << "Executing sequentially...\n";
    auto start_seq = std::chrono::high_resolution_clock::now();
    auto julia_image_seq = generate_julia(
        std::execution::seq, width, height, x_min, x_max, y_min, y_max, c, 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("julia_sequential", julia_image_seq);

    // --- Parallel Execution ---
    std::cout << "Executing in parallel...\n";
    auto start_par = std::chrono::high_resolution_clock::now();
    auto julia_image_par = generate_julia(
        std::execution::par_unseq, width, height, x_min, x_max, y_min, y_max, c, 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("julia_parallel", julia_image_par);

    std::cout << "Julia performance improvement: " << (diff_seq.count() / diff_par.count()) << "x\n";
    std::cout << "-------------------------------------------\n\n";
}

int main(int argc, char* argv[])
{
    TinyDIP::Timer timer1;
    generate_julia_set();
    return EXIT_SUCCESS;
}

The output image:

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Generate Mandelbrot Fractal Image in C++ and

  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_julia and generate_fractal_color_map functions for generating Julia Fractal image in this post.

• Why a new review is being asked for?

  Please review the implementation of generate_julia template function, map_iterations_to_color function, and its tests. Moreover, please check the modified Image class with pixels_with_coordinates member function.