rand Template Function Implementation for 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 mimic Matlab's rand function into TinyDIP library.

The experimental implementation

• rand template function implementation

  namespace TinyDIP { template<typename ElementT = double, std::same_as<std::size_t>... Sizes> constexpr static auto rand(Sizes... sizes) { auto output = zeros<ElementT>(sizes...); auto image_data = output.getImageData(); // Reference: https://stackoverflow.com/a/23143753/6667035 // First create an instance of an engine. std::random_device rnd_device; // Specify the engine and distribution. std::mt19937 mersenne_engine {rnd_device()}; // Generates random integers std::uniform_real_distribution<double> dist {0, 1}; auto gen = [&](){ return dist(mersenne_engine); }; std::generate(std::ranges::begin(image_data), std::ranges::end(image_data), gen); output = Image<ElementT>(image_data, sizes...); return output; } } 

• zeros template function implementation

  namespace TinyDIP { // zeros template function implementation template<typename ElementT, std::same_as<std::size_t>... Sizes> constexpr static auto zeros(Sizes... sizes) { auto output = Image<ElementT>(sizes...); return output; } } 

• Image class implementation (in file image.h)

  namespace TinyDIP { 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(const std::vector<std::size_t>& sizes) { if (sizes.empty()) { throw std::runtime_error("Image size vector is empty!"); } size = std::move(sizes); image_data.resize( std::reduce( std::ranges::cbegin(sizes), std::ranges::cend(sizes), std::size_t{1}, std::multiplies<>() )); } template<std::ranges::input_range Range, std::same_as<std::size_t>... Sizes> Image(const Range& input, Sizes... sizes): size{sizes...}, image_data(begin(input), end(input)) { if (image_data.size() != (1 * ... * sizes)) { throw std::runtime_error("Image data input and the given size are mismatched!"); } } template<std::same_as<std::size_t>... Sizes> Image(std::vector<ElementT>&& input, Sizes... sizes): size{sizes...}, image_data(begin(input), end(input)) { if (input.empty()) { throw std::runtime_error("Input vector is empty!"); } if (image_data.size() != (1 * ... * sizes)) { throw std::runtime_error("Image data input and the given size are mismatched!"); } } Image(const std::vector<ElementT>& input, const std::vector<std::size_t>& sizes) { if (input.empty()) { throw std::runtime_error("Input vector is empty!"); } size = std::move(sizes); image_data = std::move(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<ElementT>& input, std::size_t newWidth, std::size_t newHeight) { if (input.empty()) { throw std::runtime_error("Input vector is empty!"); } size.reserve(2); size.emplace_back(newWidth); size.emplace_back(newHeight); if (input.size() != newWidth * newHeight) { throw std::runtime_error("Image data input and the given size are mismatched!"); } image_data = std::move(input); // Reference: https://stackoverflow.com/a/51706522/6667035 } 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(input), std::ranges::end(input)); // flatten } return; } // at template function implementation template<typename... 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<typename... Args> constexpr ElementT const& at(const Args... indexInput) const { checkBoundary(indexInput...); constexpr std::size_t n = sizeof...(Args); if(n != size.size()) { throw std::runtime_error("Dimensionality mismatched!"); } 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++]; }; (update_position(indexInput), ...); return image_data[position]; } // at_without_boundary_check template function implementation template<typename... Args> constexpr ElementT& at_without_boundary_check(const Args... indexInput) { return const_cast<ElementT&>(static_cast<const Image &>(*this).at_without_boundary_check(indexInput...)); } template<typename... Args> constexpr ElementT const& at_without_boundary_check(const Args... indexInput) const { 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++]; }; (update_position(indexInput), ...); return image_data[position]; } // 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) noexcept { return image_data[index]; } // set template function implementation template<class TupleT> requires(is_tuple<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), [](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), 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 size[1]; } // 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 void print(std::string separator = "\t", std::ostream& os = std::cout) const { if(size.size() == 1) { for(std::size_t x = 0; x < size[0]; ++x) { // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number os << +at(x) << separator; } os << "\n"; } else if(size.size() == 2) { for (std::size_t y = 0; y < size[1]; ++y) { for (std::size_t x = 0; x < size[0]; ++x) { // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number os << +at(x, y) << separator; } os << "\n"; } os << "\n"; } else if (size.size() == 3) { for(std::size_t z = 0; z < size[2]; ++z) { for (std::size_t y = 0; y < size[1]; ++y) { for (std::size_t x = 0; x < size[0]; ++x) { // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number os << +at(x, y, z) << separator; } os << "\n"; } os << "\n"; } os << "\n"; } else if (size.size() == 4) { for(std::size_t a = 0; a < size[3]; ++a) { os << "group = " << a << "\n"; for(std::size_t z = 0; z < size[2]; ++z) { for (std::size_t y = 0; y < size[1]; ++y) { for (std::size_t x = 0; x < size[0]; ++x) { // Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number os << +at(x, y, z, a) << separator; } os << "\n"; } os << "\n"; } os << "\n"; } os << "\n"; } } // 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 < size[1]; ++y) { for (std::size_t x = 0; x < size[0]; ++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; } Image<ElementT>& setAllValue(const ElementT input) { std::fill(std::ranges::begin(image_data), std::ranges::end(image_data), input); 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); } #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; 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(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(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 }; } 

The usage of rand template function:

void randFunctionTest( const std::size_t sizex = 10, const std::size_t sizey = 10) { auto image1 = TinyDIP::rand(sizex, sizey); image1.print(); return; } int main() { auto start = std::chrono::system_clock::now(); randFunctionTest(); 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); if (elapsed_seconds.count() != 1) { std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << " seconds.\n"; } else { std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << " second.\n"; } return EXIT_SUCCESS; } 

The output of the test code above:

0.581644 0.684941 0.515366 0.304047 0.521844 0.8678 0.183192 0.073638 0.0496283 0.312203 0.39244 0.415687 0.465747 0.859637 0.407833 0.25225 0.582787 0.0317016 0.181481 0.505746 0.660929 0.487896 0.767568 0.0589383 0.636377 0.499133 0.00633735 0.602423 0.520545 0.713378 0.372946 0.642914 0.560535 0.250965 0.0199227 0.876541 0.460131 0.114159 0.110129 0.698122 0.9975 0.723339 0.31677 0.45846 0.89344 0.811311 0.405475 0.944949 0.068755 0.542503 0.682108 0.389317 0.105687 0.792593 0.867703 0.167464 0.614758 0.618699 0.241497 0.320009 0.618365 0.236864 0.430887 0.463439 0.340202 0.547175 0.575041 0.650054 0.530186 0.0566746 0.0861456 0.27647 0.469531 0.175167 0.639957 0.478262 0.967421 0.835256 0.538194 0.171175 0.176912 0.747396 0.485279 0.632251 0.581173 0.477933 0.472728 0.15672 0.520322 0.199587 0.88586 0.480489 0.207542 0.502268 0.99326 0.384412 0.484337 0.775651 0.219304 0.666175 Computation finished at Tue Dec 17 12:12:36 2024 elapsed time: 0.000118999 seconds. 

Godbolt link is here.

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 am trying to implement rand template function in this post.

• Why a new review is being asked for?

  Please review the implementation of rand template function and its tests.

Answer 1

I would suggest splitting rand() into two functions, with one taking a random number generator and the other one using that:

template<typename T = double, typename Urbg, std::same_as<std::size_t>... Size> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen, Size... size) { // ... } template<typename T = double, std::same_as<std::size_t>... Size> // can't be constexpr, because of std::random_device inline auto rand(Size... size) { return rand<T>(std::mt19937{std::random_device{}()}, size...); } 

This adds another function, but simplifies both functions, and gives you the freedom to use other random number generators. It also gives you a function interface similar to the seventh one listed in the MATLAB reference.

Now, right now, rand() (and zeroes()) takes any number of size arguments… from zero to infinity. If you really want to mimic the MATLAB interface, their rand() comes in 3 different flavours (to start):

1. takes no arguments, returns a random number
2. takes one argument n, returns a \$n \times n\$ 2D matrix of random numbers; and
3. takes two or more arguments, uses them as extents.

Personally, I don’t think that’s a great interface, but if you want to mimic it, you can do so with three functions (plus a fourth without the random generator that delegates to the three):

template<typename T = double, typename Urbg> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen) -> T; template<typename T = double, typename Urbg> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen, std::size_t n) -> Image<T>; template<typename T = double, typename Urbg, std::same_as<std::size_t>... Size> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> and (sizeof(Size) >= 2) constexpr auto rand(Urbg&& gen, Size... size) -> Image<T>; // OR: template<typename T = double, typename Urbg, std::same_as<std::size_t>... Size> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen, std::size_t size_1, std::size_t size_2, Size... size) -> Image<T>; 

This plus the function with out the random number generator (from above) gives you effectively the entire MATLAB API of rand(), except for the prototyping overload. (There is no need for the typename overload, because C++ is intrinsically strongly-typed, but if you really want, you could make an analogue of it.) You could also add the prototyping overload, if you wanted.

There is only one more little wrinkle, and that is how to constrain the element type. If you are using std::uniform_real_distribution, then the element type must be float, double, or long double. Not int, not std::complex<double>, not any pixel descriptor class… justfloat, double, or long double. Whether this is a problem is for you to decide.

What I would suggest is making a concept to check for valid image element types, and have a sub-concept for just the floating point types:

template<typename T> concept image_element_standard_floating_point_type = std::same_as<double, T> or std::same_as<float, T> or std::same_as<long double, T> ; template<typename T> concept image_element_type = image_element_standard_floating_point_type<T> // or image_element_any_other_category_of_type<T> ; 

Now you can constrain with the sub-concept:

template<image_element_standard_floating_point_type T = double, typename Urbg> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen) -> T; template<image_element_standard_floating_point_type T = double, typename Urbg> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& gen, std::size_t n) -> Image<T>; template<image_element_standard_floating_point_type T = double, typename Urbg, std::same_as<std::size_t>... Size> requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> and (sizeof(Size) >= 2) constexpr auto rand(Urbg&& gen, Size... size) -> Image<T>; 

If you want to be able to generate other image types, you can use other sub-concepts, and other random number distributions, or maybe generator functions. For example, you could have a image_element_standard_integer_type concept that checks for int, unsigned, etc., and an overload of rand() constrained with that that uses std::uniform_int_distribution.

On to reviewing the actual code:

 template<typename ElementT = double, std::same_as<std::size_t>... Sizes> constexpr static auto rand(Sizes... sizes) { auto output = zeros<ElementT>(sizes...); auto image_data = output.getImageData(); // Reference: https://stackoverflow.com/a/23143753/6667035 // First create an instance of an engine. std::random_device rnd_device; // Specify the engine and distribution. std::mt19937 mersenne_engine {rnd_device()}; // Generates random integers std::uniform_real_distribution<double> dist {0, 1}; auto gen = [&](){ return dist(mersenne_engine); }; std::generate(std::ranges::begin(image_data), std::ranges::end(image_data), gen); output = Image<ElementT>(image_data, sizes...); return output; } 

The static there serves absolutely no purpose except to potentially unnecessarily duplicate the function in multiple translation units.

So, the way the function is implemented is it uses zeros() generate an all-zero image of the right size… then it copies the image data out of that… overwrites it with random values… then copies the data yet again into a new image… and finally overwrites the zero image. This seems incredibly wasteful.

A better plan seems to be to create a vector of the right size (which, unfortunately, also gets unnecessarily zeroed, but that is a shortcoming of std::vector that we can’t easily avoid here), write the random values into that, and them move that vector directly into the image:

template< image_element_standard_floating_point_type ElementT = double, typename Urbg, std::same_as<std::size_t>... Sizes > requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> and (sizeof(Size) >= 2) constexpr auto rand(Urbg&& urbg, Sizes... sizes) { auto image_data = std::vector<ElementT>((... * sizes)); // if you want to allow empty sizes lists, then: `(1uz * ... * sizes)` auto dist = std::uniform_real_distribution<ElementT>{}; std::ranges::generate(image_data, [&dist, &urbg]() { return dist(urbg); }); return Image<ElementT>{std::move(image_data), sizes...}; } 

Assuming RVO, no image data is ever copied. The only “waste” is that the image data is set to all zeros before being overwritten with random data; this is a shortcoming of vector that can be worked around… but not here. (You’d probably need to change to std::pmr::vector everywhere for that.)

As for testingrand()… it is generally impossible to write comprehensive unit tests for anything that generates randomness. The best you can do is test that the generated image size(s) is correct, and that every element is in the range \$[0,1)\$. Unless you want to expose and codify the internals of the function, that’s about the best you can do in a proper unit test. You can do statistical checks on the data to see if it is probably random, but you don’t want to do that in a unit test, because it is not deterministic.

Just for clarity because I suggested the other functions, I’ll show what their implementation might look like, because I find code generally more clear than words:

template<typename T> concept image_element_standard_floating_point_type = std::same_as<double, T> or std::same_as<float, T> or std::same_as<long double, T> ; template<typename T> concept image_element_type = image_element_standard_floating_point_type<T> // or image_element_any_other_category_of_type<T> ; template< image_element_standard_floating_point_type ElementT = double, typename Urbg, std::same_as<std::size_t>... Sizes > requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> and (sizeof(Size) >= 2) constexpr auto rand(Urbg&& urbg, Sizes... sizes) { auto image_data = std::vector<ElementT>((... * sizes)); // if you want to allow empty sizes lists, then: `(1uz * ... * sizes)` auto dist = std::uniform_real_distribution<ElementT>{}; std::ranges::generate(image_data, [&dist, &urbg]() { return dist(urbg); }); return Image<ElementT>{std::move(image_data), sizes...}; } template< image_element_standard_floating_point_type ElementT = double, typename Urbg, > requires std::uniform_random_bit_generator<std::remove_reference_t<Urbg>> constexpr auto rand(Urbg&& urbg, std::size_t size) { return rand<ElementT>(std::forward<Urbg>(urbg), size, size); } template< image_element_standard_floating_point_type ElementT = double, typename Urbg > constexpr auto rand(Urbg&& urbg) { auto dist = std::uniform_real_distribution<ElementT>{}; return dist(urbg); } template<image_element_type T = double, std::same_as<std::size_t>... Size> inline auto rand(Size... size) { return rand<T>(std::mt19937{std::random_device{}()}, size...); } 

If you also want to support extended floating point types, integer types, complex types, or anything else, it’s very possible, but you would need to specify what you would expect. (For example, would rand<int>(...) produce an image where each element is in the range [0,INT_MAX], [0,INT_MAX), [INT_MIN,INT_MAX], (INT_MIN,INT_MAX], or something else? Would complex numbers have only the real part randomized, or both the real and imaginary parts?) What you could do is have another overload set that, instead of taking a uniform random bit generator, takes a generator function. So rand(gen, sizes...) would generate an image with the sizes given, where each element is generated by gen().