Dictionary based non-local mean implementation in C++

Question

This is a follow-up question for Manhattan distance calculation between two images in C++ and Dictionary based non-local mean implementation in Matlab. For learning C++20 and researching purposes, I am attempting to implement function create_dictionary and dictionaryBasedNonlocalMean which purpose are similar to the linked follow-up question, i.e. calculate dictionary based non-local mean. In mathematical form, output of function dictionaryBasedNonlocalMean is calculated with the following way.

$$output = K\sum_{{i = 1}}^{{N}_{D}} \left[ G_{\mu = 0, \sigma} (\left\|input - X(i)\right\|_{1}) \cdot Y(i) \right] $$

where N_D is the count of X-Y pairs (the cardinality of set X and set Y) in the dictionary and the Gaussian function

$$G_{\mu = 0, \sigma} (\left\|input - X(i)\right\|_{1}) = \left.e^{\frac{-(\left\|input - X(i)\right\|_{1} - \mu)^2}{2 {\sigma}^{2}}}\right|_{\mu=0} $$

Moreover, K is a normalization factor

$$K = \frac{1}{\sum_{{i = 1}}^{{N}_{D}} G_{\mu = 0, \sigma} (\left\|input - X(i)\right\|_{1})} $$

Illustration

The experimental implementation

• create_dictionary template function implementation: The purpose of this function is similar to Matlab version CreateDictionary in the linked follow-up question. The N_D_ paired codewords are created in std::tuple<std::vector<...>, std::vector<...>> structure.

  template<class ElementT = double>
  constexpr static auto create_dictionary(const std::size_t ND, const std::size_t xsize, const std::size_t ysize, const std::size_t zsize)
  {
   auto code_words_x = TinyDIP::n_dim_vector_generator<1>(std::vector<TinyDIP::Image<ElementT>>(), ND);
   auto code_words_y = TinyDIP::n_dim_vector_generator<1>(std::vector<TinyDIP::Image<ElementT>>(), ND);
   for (std::size_t i = 0; i < ND; ++i)
   {
   code_words_x[i] = TinyDIP::n_dim_vector_generator<1>(
   TinyDIP::Image<ElementT>(xsize, ysize, static_cast<ElementT>(i) / static_cast<ElementT>(ND)), zsize);
   code_words_y[i] = TinyDIP::n_dim_vector_generator<1>(
   TinyDIP::Image<ElementT>(xsize, ysize, 1.0 + static_cast<ElementT>(i) / static_cast<ElementT>(ND)), zsize);
   }
   return std::make_tuple(code_words_x, code_words_y);
  }
  

• dictionaryBasedNonlocalMean template function implementation:

  template<class ExPo, class ElementT = double>
  requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
  constexpr static auto dictionaryBasedNonlocalMean( ExPo execution_policy,
   const std::tuple<
   std::vector<std::vector<TinyDIP::Image<ElementT>>>,
   std::vector<std::vector<TinyDIP::Image<ElementT>>>
   >& dictionary,
   const std::vector<TinyDIP::Image<ElementT>>& input,
   const double gaussian_sigma = 3.0,
   const double gaussian_mean = 0,
   const double threshold = 1e-160) noexcept
  {
   std::vector<TinyDIP::Image<ElementT>> output = 
   TinyDIP::n_dim_vector_generator<1>(
   TinyDIP::Image(input[0].getWidth(), input[0].getHeight(), 0.0), input.size());
   auto code_words_x = std::get<0>(dictionary);
   auto code_words_y = std::get<1>(dictionary);
   if (code_words_x.size() != code_words_y.size())
   {
   throw std::runtime_error("Size of data in dictionary incorrect.");
   }
   auto weights = TinyDIP::recursive_transform<1>(
   execution_policy,
   [&](auto&& element)
   {
   return TinyDIP::normalDistribution1D(
   TinyDIP::recursive_reduce(
   TinyDIP::recursive_transform<1>(
   [&](auto&& each_plane_x, auto&& each_plane_input) { return TinyDIP::manhattan_distance(each_plane_x, each_plane_input); },
   element, input),
   ElementT{}) + gaussian_mean,
   gaussian_sigma);
   }, code_words_x);
   auto sum_of_weights = TinyDIP::recursive_reduce(weights, ElementT{});
   std::cout << "sum_of_weights: " << std::to_string(sum_of_weights) << '\n';
   if (sum_of_weights < threshold)
   {
   return output;
   }
   auto outputs = TinyDIP::recursive_transform<1>(
   [&](auto&& input1, auto&& input2)
   {
   return TinyDIP::multiplies(input1, TinyDIP::n_dim_vector_generator<1>(TinyDIP::Image(input1[0].getWidth(), input1[0].getHeight(), input2), input1.size()));
   }, code_words_y, weights);
   for (std::size_t i = 0; i < outputs.size(); ++i)
   {
   output = TinyDIP::plus(output, outputs[i]);
   }
   output = TinyDIP::divides(output, TinyDIP::n_dim_vector_generator<1>(TinyDIP::Image(output[0].getWidth(), output[0].getHeight(), sum_of_weights), output.size()));
   return output;
  }
  

  • manhattan_distance template function implementation:

  template<arithmetic ElementT = double>
  constexpr static ElementT manhattan_distance(const Image<ElementT>& input1, const Image<ElementT>& input2)
  {
   is_size_same(input1, input2);
   return recursive_reduce(difference(input1, input2).getImageData(), ElementT{});
  }
  

• difference template function implementation:

  template<arithmetic ElementT = double>
  constexpr static auto difference(const Image<ElementT>& input1, const Image<ElementT>& input2)
  {
   return pixelwiseOperation([](auto&& element1, auto&& element2) { return std::abs(element1 - element2); }, input1, input2);
  }
  

• pixelwiseOperation template function implementation:

  template<std::size_t unwrap_level = 1, class... Args>
  constexpr static auto pixelwiseOperation(auto op, const Args&... inputs)
  {
   auto output = Image(
   recursive_transform<unwrap_level>(
   [&](auto&& element1, auto&&... elements) 
   {
   return op(element1, elements...);
   },
   inputs.getImageData()...),
   first_of(inputs...).getWidth(),
   first_of(inputs...).getHeight());
   return output;
  }
  template<std::size_t unwrap_level = 1, class ExPo, class InputT>
  requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
  constexpr static auto pixelwiseOperation(ExPo execution_policy, auto op, const Image<InputT>& input1)
  {
   auto output = Image(
   recursive_transform<unwrap_level>(
   execution_policy,
   [&](auto&& element1) 
   {
   return op(element1);
   },
   (input1.getImageData())),
   input1.getWidth(),
   input1.getHeight());
   return output;
  }
  

Full Testing Code

Tests for dictionaryBasedNonlocalMean template function:

#include <algorithm>
#include <array>
#include <cassert>
#include <chrono>
#include <cmath>
#include <concepts>
#include <exception>
#include <execution>
#include <fstream>
#include <functional>
#include <iostream>
#include <iterator>
#include <mutex>
#include <numeric>
#include <ranges>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
using BYTE = unsigned char;
struct RGB
{
 BYTE channels[3];
};
using GrayScale = BYTE;
namespace TinyDIP
{
 #define is_size_same(x, y) {assert(x.getWidth() == y.getWidth()); assert(x.getHeight() == y.getHeight());}
 // Reference: https://stackoverflow.com/a/58067611/6667035
 template <typename T>
 concept arithmetic = std::is_arithmetic_v<T>;
 // recursive_depth function implementation
 template<typename T>
 constexpr std::size_t recursive_depth()
 {
 return 0;
 }
 template<std::ranges::input_range Range>
 constexpr std::size_t recursive_depth()
 {
 return recursive_depth<std::ranges::range_value_t<Range>>() + 1;
 }
 template<std::size_t index = 1, typename Arg, typename... Args>
 constexpr static auto& get_from_variadic_template(const Arg& first, const Args&... inputs)
 {
 if constexpr (index > 1)
 return get_from_variadic_template<index - 1>(inputs...);
 else
 return first;
 }
 template<typename... Args>
 constexpr static auto& first_of(Args&... inputs) {
 return get_from_variadic_template<1>(inputs...);
 }
 // recursive_reduce implementation
 // Reference: https://codereview.stackexchange.com/a/251310/231235
 template<class T, class ValueType, class Function = std::plus<ValueType>>
 constexpr auto recursive_reduce(const T& input, ValueType init, const Function& f)
 {
 return f(init, input);
 }
 template<std::ranges::range Container, class ValueType, class Function = std::plus<ValueType>>
 constexpr auto recursive_reduce(const Container& input, ValueType init, const Function& f = std::plus<ValueType>())
 {
 for (const auto& element : input) {
 auto result = recursive_reduce(element, ValueType{}, f);
 init = f(init, result);
 }
 return init;
 }
 // recursive_invoke_result_t implementation
 template<std::size_t, typename, typename>
 struct recursive_invoke_result { };
 template<typename T, typename F>
 struct recursive_invoke_result<0, F, T> { using type = std::invoke_result_t<F, T>; };
 template<std::size_t unwrap_level, typename F, template<typename...> typename Container, typename... Ts>
 requires (std::ranges::input_range<Container<Ts...>> &&
 requires { typename recursive_invoke_result<unwrap_level - 1, F, std::ranges::range_value_t<Container<Ts...>>>::type; })
 struct recursive_invoke_result<unwrap_level, F, Container<Ts...>>
 {
 using type = Container<typename recursive_invoke_result<unwrap_level - 1, F, std::ranges::range_value_t<Container<Ts...>>>::type>;
 };
 template<std::size_t unwrap_level, typename F, typename T>
 using recursive_invoke_result_t = typename recursive_invoke_result<unwrap_level, F, T>::type;
 // recursive_variadic_invoke_result_t implementation
 template<std::size_t, typename, typename, typename...>
 struct recursive_variadic_invoke_result { };
 template<typename F, class...Ts1, template<class...>class Container1, typename... Ts>
 struct recursive_variadic_invoke_result<1, F, Container1<Ts1...>, Ts...>
 {
 using type = Container1<std::invoke_result_t<F,
 std::ranges::range_value_t<Container1<Ts1...>>,
 std::ranges::range_value_t<Ts>...>>;
 };
 template<std::size_t unwrap_level, typename F, class...Ts1, template<class...>class Container1, typename... Ts>
 requires ( std::ranges::input_range<Container1<Ts1...>> &&
 requires { typename recursive_variadic_invoke_result<
 unwrap_level - 1,
 F,
 std::ranges::range_value_t<Container1<Ts1...>>,
 std::ranges::range_value_t<Ts>...>::type; }) // The rest arguments are ranges
 struct recursive_variadic_invoke_result<unwrap_level, F, Container1<Ts1...>, Ts...>
 {
 using type = Container1<
 typename recursive_variadic_invoke_result<
 unwrap_level - 1,
 F,
 std::ranges::range_value_t<Container1<Ts1...>>,
 std::ranges::range_value_t<Ts>...
 >::type>;
 };
 template<std::size_t unwrap_level, typename F, typename T1, typename... Ts>
 using recursive_variadic_invoke_result_t = typename recursive_variadic_invoke_result<unwrap_level, F, T1, Ts...>::type;
 template<typename OutputIt, typename NAryOperation, typename InputIt, typename... InputIts>
 OutputIt transform(OutputIt d_first, NAryOperation op, InputIt first, InputIt last, InputIts... rest) {
 while (first != last) {
 *d_first++ = op(*first++, (*rest++)...);
 }
 return d_first;
 }
 // recursive_transform for the multiple parameters cases (the version with unwrap_level)
 template<std::size_t unwrap_level = 1, class F, class Arg1, class... Args>
 constexpr auto recursive_transform(const F& f, const Arg1& arg1, const Args&... args)
 {
 if constexpr (unwrap_level > 0)
 {
 static_assert(unwrap_level <= recursive_depth<Arg1>(),
 "unwrap level higher than recursion depth of input");
 recursive_variadic_invoke_result_t<unwrap_level, F, Arg1, Args...> output{};
 transform(
 std::inserter(output, std::ranges::end(output)),
 [&f](auto&& element1, auto&&... elements) { return recursive_transform<unwrap_level - 1>(f, element1, elements...); },
 std::ranges::cbegin(arg1),
 std::ranges::cend(arg1),
 std::ranges::cbegin(args)...
 );
 return output;
 }
 else
 {
 return f(arg1, args...);
 }
 }
 // recursive_transform implementation (the version with unwrap_level, with execution policy)
 template<std::size_t unwrap_level = 1, class ExPo, class T, class F>
 requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
 constexpr auto recursive_transform(ExPo execution_policy, const F& f, const T& input)
 {
 if constexpr (unwrap_level > 0)
 {
 static_assert(unwrap_level <= recursive_depth<T>(),
 "unwrap level higher than recursion depth of input");
 recursive_invoke_result_t<unwrap_level, F, T> output{};
 output.resize(input.size());
 std::mutex mutex;
 std::transform(execution_policy, std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
 [&](auto&& element)
 {
 std::lock_guard lock(mutex);
 return recursive_transform<unwrap_level - 1>(execution_policy, f, element);
 });
 return output;
 }
 else
 {
 return f(input);
 }
 }
 template<std::size_t dim, class T>
 constexpr auto n_dim_vector_generator(T input, std::size_t times)
 {
 if constexpr (dim == 0)
 {
 return input;
 }
 else
 {
 auto element = n_dim_vector_generator<dim - 1>(input, times);
 std::vector<decltype(element)> output(times, element);
 return output;
 }
 }
 template <typename ElementT>
 class Image
 {
 public:
 Image() = default;
 Image(const std::size_t width, const std::size_t height):
 width(width),
 height(height),
 image_data(width * height) { }
 Image(const std::size_t width, const std::size_t height, const ElementT initVal):
 width(width),
 height(height),
 image_data(width * height, initVal) {}
 Image(const std::vector<ElementT> input, std::size_t newWidth, std::size_t newHeight):
 width(newWidth),
 height(newHeight)
 {
 if (input.size() != newWidth * newHeight)
 {
 throw std::runtime_error("Image data input and the given size are mismatched!");
 }
 image_data = std::move(input);
 }
 constexpr ElementT& at(const unsigned int x, const unsigned int y)
 { 
 checkBoundary(x, y);
 return image_data[y * width + x];
 }
 constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
 {
 checkBoundary(x, y);
 return image_data[y * width + x];
 }
 constexpr std::size_t getWidth() const
 {
 return width;
 }
 constexpr std::size_t getHeight() const
 {
 return height;
 }
 constexpr auto getSize()
 {
 return std::make_tuple(width, height);
 }
 std::vector<ElementT> const& getImageData() const { return image_data; } // expose the internal data
 void print(std::string separator = "\t", std::ostream& os = std::cout) const
 {
 for (std::size_t y = 0; y < height; ++y)
 {
 for (std::size_t x = 0; x < width; ++x)
 {
 // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
 os << +at(x, y) << separator;
 }
 os << "\n";
 }
 os << "\n";
 return;
 }
 // Enable this function if ElementT = RGB
 void print(std::string separator = "\t", std::ostream& os = std::cout) const
 requires(std::same_as<ElementT, RGB>)
 {
 for (std::size_t y = 0; y < height; ++y)
 {
 for (std::size_t x = 0; x < width; ++x)
 {
 os << "( ";
 for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
 {
 // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
 os << +at(x, y).channels[channel_index] << separator;
 }
 os << ")" << separator;
 }
 os << "\n";
 }
 os << "\n";
 return;
 }
 friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
 {
 const std::string separator = "\t";
 for (std::size_t y = 0; y < rhs.height; ++y)
 {
 for (std::size_t x = 0; x < rhs.width; ++x)
 {
 // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
 os << +rhs.at(x, y) << separator;
 }
 os << "\n";
 }
 os << "\n";
 return os;
 }
 Image<ElementT>& operator+=(const Image<ElementT>& rhs)
 {
 assert(rhs.width == this->width);
 assert(rhs.height == this->height);
 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)
 {
 assert(rhs.width == this->width);
 assert(rhs.height == this->height);
 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)
 {
 assert(rhs.width == this->width);
 assert(rhs.height == this->height);
 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)
 {
 assert(rhs.width == this->width);
 assert(rhs.height == this->height);
 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+(const Image<ElementT>& input1, const Image<ElementT>& input2);
 friend Image<ElementT> operator-(const Image<ElementT>& input1, const Image<ElementT>& input2);
 Image<ElementT>& operator=(Image<ElementT> const& input) = default; // Copy Assign
 Image<ElementT>& operator=(Image<ElementT>&& other) = default; // Move Assign
 Image(const Image<ElementT> &input) = default; // Copy Constructor
 Image(Image<ElementT> &&input) = default; // Move Constructor
 
 private:
 std::size_t width;
 std::size_t height;
 std::vector<ElementT> image_data;
 void checkBoundary(const size_t x, const size_t y) const
 {
 assert(x < width);
 assert(y < height);
 }
 };
 template<class ElementT>
 Image<ElementT> operator+(const Image<ElementT>& input1, const Image<ElementT>& input2)
 {
 return plus(input1, input2);
 }
 template<class ElementT>
 Image<ElementT> operator-(const Image<ElementT>& input1, const Image<ElementT>& input2)
 {
 return subtract(input1, input2);
 }
 template<typename T>
 T normalDistribution1D(const T x, const T standard_deviation)
 {
 return std::exp(-x * x / (2 * standard_deviation * standard_deviation));
 }
 template<std::size_t unwrap_level = 1, class... Args>
 constexpr static auto pixelwiseOperation(auto op, const Args&... inputs)
 {
 auto output = Image(
 recursive_transform<unwrap_level>(
 [&](auto&& element1, auto&&... elements) 
 {
 return op(element1, elements...);
 },
 inputs.getImageData()...),
 first_of(inputs...).getWidth(),
 first_of(inputs...).getHeight());
 return output;
 }
 template<std::size_t unwrap_level = 1, class ExPo, class InputT>
 requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
 constexpr static auto pixelwiseOperation(ExPo execution_policy, auto op, const Image<InputT>& input1)
 {
 auto output = Image(
 recursive_transform<unwrap_level>(
 execution_policy,
 [&](auto&& element1) 
 {
 return op(element1);
 },
 (input1.getImageData())),
 input1.getWidth(),
 input1.getHeight());
 return output;
 }
 template<class InputT>
 constexpr static Image<InputT> plus(const Image<InputT>& input1)
 {
 return input1;
 }
 template<class InputT, class... Args>
 constexpr static Image<InputT> plus(const Image<InputT>& input1, const Args&... inputs)
 {
 return pixelwiseOperation(std::plus<>{}, input1, plus(inputs...));
 }
 template<class InputT, class... Args>
 constexpr static auto plus(const std::vector<Image<InputT>>& input1, const Args&... inputs)
 {
 return recursive_transform<1>(
 [](auto&& input1_element, auto&&... inputs_element)
 {
 return plus(input1_element, inputs_element...);
 }, input1, inputs...);
 }
 template<class InputT>
 constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
 {
 is_size_same(input1, input2);
 return pixelwiseOperation(std::minus<>{}, input1, input2);
 }
 template<class InputT>
 constexpr static Image<InputT> subtract(const std::vector<Image<InputT>>& input1, const std::vector<Image<InputT>>& input2)
 {
 return recursive_transform<1>(
 [](auto&& input1_element, auto&& input2_element)
 {
 return subtract(input1_element, input2_element);
 }, input1, input2);
 }
 template<class InputT = RGB>
 requires (std::same_as<InputT, RGB>)
 constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
 {
 is_size_same(input1, input2);
 Image<InputT> output(input1.getWidth(), input1.getHeight());
 for (std::size_t y = 0; y < input1.getHeight(); ++y)
 {
 for (std::size_t x = 0; x < input1.getWidth(); ++x)
 {
 for(std::size_t channel_index = 0; channel_index < 3; ++channel_index)
 {
 output.at(x, y).channels[channel_index] = 
 std::clamp(
 input1.at(x, y).channels[channel_index] - 
 input2.at(x, y).channels[channel_index],
 0,
 255);
 }
 }
 }
 return output;
 }
 template<class InputT>
 constexpr static Image<InputT> multiplies(const Image<InputT>& input1, const Image<InputT>& input2)
 {
 return pixelwiseOperation(std::multiplies<>{}, input1, input2);
 }
 template<class ExPo, class InputT>
 requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
 constexpr static Image<InputT> multiplies(ExPo execution_policy, const Image<InputT>& input1, const Image<InputT>& input2)
 {
 return pixelwiseOperation(execution_policy, std::multiplies<>{}, input1, input2);
 }
 template<class InputT, class... Args>
 constexpr static Image<InputT> multiplies(const Image<InputT>& input1, const Args&... inputs)
 {
 return pixelwiseOperation(std::multiplies<>{}, input1, multiplies(inputs...));
 }
 template<class InputT, class... Args>
 constexpr static auto multiplies(const std::vector<Image<InputT>>& input1, const Args&... inputs)
 {
 return recursive_transform<1>(
 [](auto&& input1_element, auto&&... inputs_element)
 {
 return multiplies(input1_element, inputs_element...);
 }, input1, inputs...);
 }
 template<class InputT>
 constexpr static Image<InputT> divides(const Image<InputT>& input1, const Image<InputT>& input2)
 {
 return pixelwiseOperation(std::divides<>{}, input1, input2);
 }
 template<class InputT>
 constexpr static auto divides(const std::vector<Image<InputT>>& input1, const std::vector<Image<InputT>>& input2)
 {
 return recursive_transform<1>(
 [](auto&& input1_element, auto&& input2_element)
 {
 return divides(input1_element, input2_element);
 }, input1, input2);
 }
 template<class ExPo, class InputT>
 requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
 constexpr static Image<InputT> divides(ExPo execution_policy, const Image<InputT>& input1, const Image<InputT>& input2)
 {
 return pixelwiseOperation(execution_policy, std::divides<>{}, input1, input2);
 }
 template<arithmetic ElementT = double>
 constexpr static auto abs(const Image<ElementT>& input)
 {
 return pixelwiseOperation([](auto&& element) { return std::abs(element); }, input);
 }
 template<class ExPo, arithmetic ElementT = double>
 requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
 constexpr static auto abs(ExPo execution_policy, const Image<ElementT>& input)
 {
 return pixelwiseOperation(execution_policy, [](auto&& element) { return std::abs(element); }, input);
 }
 template<arithmetic ElementT = double>
 constexpr static auto difference(const Image<ElementT>& input1, const Image<ElementT>& input2)
 {
 return pixelwiseOperation([](auto&& element1, auto&& element2) { return std::abs(element1 - element2); }, input1, input2);
 }
 template<arithmetic ElementT = double>
 constexpr static ElementT manhattan_distance(const Image<ElementT>& input1, const Image<ElementT>& input2)
 {
 is_size_same(input1, input2);
 return recursive_reduce(difference(input1, input2).getImageData(), ElementT{});
 }
}
/* Matlab version:
function Dictionary = CreateDictionary(ND, xsize, ysize, zsize)
 Dictionary.X = zeros(xsize, ysize, zsize, ND);
 Dictionary.Y = zeros(xsize, ysize, zsize, ND);
 for i = 1:ND
 Dictionary.X(:, :, :, i) = ones(xsize, ysize, zsize) .* (i / ND);
 Dictionary.Y(:, :, :, i) = ones(xsize, ysize, zsize) .* (1 + i / ND);
 end
end
*/
template<class ElementT = double>
constexpr static auto create_dictionary(const std::size_t ND, const std::size_t xsize, const std::size_t ysize, const std::size_t zsize)
{
 auto code_words_x = TinyDIP::n_dim_vector_generator<1>(std::vector<TinyDIP::Image<ElementT>>(), ND);
 auto code_words_y = TinyDIP::n_dim_vector_generator<1>(std::vector<TinyDIP::Image<ElementT>>(), ND);
 for (std::size_t i = 0; i < ND; ++i)
 {
 code_words_x[i] = TinyDIP::n_dim_vector_generator<1>(
 TinyDIP::Image<ElementT>(xsize, ysize, static_cast<ElementT>(i) / static_cast<ElementT>(ND)), zsize);
 code_words_y[i] = TinyDIP::n_dim_vector_generator<1>(
 TinyDIP::Image<ElementT>(xsize, ysize, 1.0 + static_cast<ElementT>(i) / static_cast<ElementT>(ND)), zsize);
 }
 return std::make_tuple(code_words_x, code_words_y);
}
/* Matlab version:
function [output] = dictionaryBasedNonlocalMean(Dictionary, input)
 gaussian_sigma = 0.1;
 gaussian_mean = 0;
 if size(Dictionary.X) ~= size(Dictionary.Y)
 disp("Size of data in dictionary incorrect.");
 output = [];
 return
 end
 [X, Y, Z, DataCount] = size(Dictionary.X);
 weightOfEach = zeros(1, DataCount);
 for i = 1:DataCount
 % Gaussian of distance between X and input
 weightOfEach(i) = gaussmf(ManhattanDistance(input, Dictionary.X(:, :, :, i)), [gaussian_sigma gaussian_mean]); 
 end
 sumOfDist = sum(weightOfEach, 'all');
 output = zeros(X, Y, Z);
 %%% if sumOfDist too small
 if (sumOfDist < 1e-160)
 fprintf("sumOfDist = %d\n", sumOfDist);
 return;
 end
 for i = 1:DataCount
 output = output + Dictionary.Y(:, :, :, i) .* weightOfEach(i);
 end
 output = output ./ sumOfDist;
end
*/
template<class ExPo, class ElementT = double>
requires (std::is_execution_policy_v<std::remove_cvref_t<ExPo>>)
constexpr static auto dictionaryBasedNonlocalMean( ExPo execution_policy,
 const std::tuple<
 std::vector<std::vector<TinyDIP::Image<ElementT>>>,
 std::vector<std::vector<TinyDIP::Image<ElementT>>>
 >& dictionary,
 const std::vector<TinyDIP::Image<ElementT>>& input,
 const double gaussian_sigma = 3.0,
 const double gaussian_mean = 0,
 const double threshold = 1e-160)
{
 std::vector<TinyDIP::Image<ElementT>> output = 
 TinyDIP::n_dim_vector_generator<1>(
 TinyDIP::Image(input[0].getWidth(), input[0].getHeight(), 0.0), input.size());
 auto code_words_x = std::get<0>(dictionary);
 auto code_words_y = std::get<1>(dictionary);
 if (code_words_x.size() != code_words_y.size())
 {
 throw std::runtime_error("Size of data in dictionary incorrect.");
 }
 auto weights = TinyDIP::recursive_transform<1>(
 execution_policy,
 [&](auto&& element)
 {
 return TinyDIP::normalDistribution1D(
 TinyDIP::recursive_reduce(
 TinyDIP::recursive_transform<1>(
 [&](auto&& each_plane_x, auto&& each_plane_input) { return TinyDIP::manhattan_distance(each_plane_x, each_plane_input); },
 element, input),
 ElementT{}) + gaussian_mean,
 gaussian_sigma);
 }, code_words_x);
 auto sum_of_weights = TinyDIP::recursive_reduce(weights, ElementT{});
 std::cout << "sum_of_weights: " << std::to_string(sum_of_weights) << '\n';
 if (sum_of_weights < threshold)
 {
 return output;
 }
 auto outputs = TinyDIP::recursive_transform<1>(
 [&](auto&& input1, auto&& input2)
 {
 return TinyDIP::multiplies(input1, TinyDIP::n_dim_vector_generator<1>(TinyDIP::Image(input1[0].getWidth(), input1[0].getHeight(), input2), input1.size()));
 }, code_words_y, weights);
 
 for (std::size_t i = 0; i < outputs.size(); ++i)
 {
 output = TinyDIP::plus(output, outputs[i]);
 }
 output = TinyDIP::divides(output, TinyDIP::n_dim_vector_generator<1>(TinyDIP::Image(output[0].getWidth(), output[0].getHeight(), sum_of_weights), output.size()));
 return output;
}
void dictionaryBasedNonLocalMeanTest()
{
 std::size_t ND = 10;
 std::size_t xsize = 8;
 std::size_t ysize = 8;
 std::size_t zsize = 1;
 std::vector<TinyDIP::Image<double>> input;
 for (std::size_t z = 0; z < zsize; ++z)
 {
 input.push_back(TinyDIP::Image(xsize, ysize, 0.66));
 }
 dictionaryBasedNonlocalMean(
 std::execution::par,
 create_dictionary(ND, xsize, ysize, zsize),
 input
 ).at(0).print();
}
int main()
{
 auto start = std::chrono::system_clock::now();
 dictionaryBasedNonLocalMeanTest();
 auto end = std::chrono::system_clock::now();
 std::chrono::duration<double> elapsed_seconds = end - start;
 std::time_t end_time = std::chrono::system_clock::to_time_t(end);
 std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
 return 0;
}

The output of the testing code above:

sum_of_weights: 1.150129
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 1.66217 
Computation finished at Sat Dec 18 02:22:04 2021
elapsed time: 0.000108249

A Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Manhattan distance calculation between two images in C++ and

  Dictionary based non-local mean implementation in Matlab

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

  As the suggestions from G. Sliepen's answer and JDługosz's answer, these changes have been made:

  • Unnecessary usages of this-> and TinyDIP:: are removed

  • The updated version of Image class has been proposed.

• Why a new review is being asked for?

  If there is any possible improvement, please let me know.

Answer 1

I don't have the time at the moment for a full review, but here are some things I noticed that may help you improve your code.

Free your functions

The operator+ and operator- don't need to be friend functions because they don't rely on any internal knowledge of the Image class. See also Klaus Iglberger's 2017 CppCon talk. In fact, they can be made even simpler per the next suggestion.

Use the canonical definition for operators

The typical way to define a function such as operator+ is to define it in terms of operator+= which makes things simple for copyable types:

template<class ElementT>
Image<ElementT> operator+(Image<ElementT> input1, const Image<ElementT>& input2)
{
 return input1 += input2;
}

Note that input1 is passed by value so a copy is made. Then we use operator+= with input2.

• JDługosz: OTOH, "hidden friends" is a recommended style for compile-time performance and avoidance of ambiguities. The idea is that the operator can only be found via ADL.
• Edward: @JDługosz If the intent is to create a templated friend, isn't the syntax is wrong in the original? Wouldn't it need a <template...> declaration above it? I was unfamiliar with the "hidden friends" idea, but I'm reading about it now. Thanks!
• JDługosz: if the non-member operator+ is defined inside the class template (via friend), then there's no need to make that function a template as in your example. That's what you need if you pull it out of the class template. (What we need is a friend public that lets you define a non-member inside a class or class template for ADL-only, but does not grant private access to the class.)
• JDługosz: @JimmyHu that operator+ should not compile! You can't write Left += Right when Left is const. So, something is really strange about the operator+= it refers to. ... yes, it is not marked const, so now I wonder why trying to call += didn't give an error... You must not have actually used that function, so it was never instantiated from the template.
• JimmyHu: @JDługosz OK, I see. I'll try to fix the issue and thank you for the explanation.