conv3 Template Function Implementation for Image in C++

Question

This is a follow-up question for conv2 Template Function Implementation for Image in C++. Besides two-dimensional convolution, I am trying to implement three-dimensional convolution for 3D image in C++. Ideally, n-dimensional convolution convn like Matlab should be implemented; however, I have no idea how to implement n-dimensional case gracefully with existed multi-dimensional Image class now.

The experimental implementation

• conv3 template function implementation (in file image_operations.h)

  namespace TinyDIP
  {
      //  conv3 template function implementation
      template<typename ElementT>
      requires(std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value)
      constexpr static auto conv3(const Image<ElementT>& x, const Image<ElementT>& y)
      {
          auto output = Image<ElementT>(
                                  x.getSize(0) + y.getSize(0) - 1,
                                  x.getSize(1) + y.getSize(1) - 1,
                                  x.getSize(2) + y.getSize(2) - 1);
  
          for (std::size_t z2 = 0; z2 < y.getSize(2); ++z2)
          {
              for (std::size_t z1 = 0; z1 < x.getSize(2); ++z1)
              {
                  for (std::size_t y2 = 0; y2 < y.getSize(1); ++y2)
                  {
                      for (std::size_t y1 = 0; y1 < x.getSize(1); ++y1)
                      {
                          for (std::size_t x2 = 0; x2 < y.getSize(0); ++x2)
                          {
                              for (std::size_t x1 = 0; x1 < x.getSize(0); ++x1)
                              {
                                  output.at(x1 + x2, y1 + y2, z1 + z2) = output.at(x1 + x2, y1 + y2, z1 + z2) + x.at(x1, y1, z1) * y.at(x2, y2, z2);
                              }
                          }
                      }
                  }
              }
          }
          return output;
      }
  }
  

Full Testing Code

The full testing code:

//  conv3 Template Function Implementation for Image in C++

#include <algorithm>
#include <cassert>
#include <cctype>
#include <chrono>
#include <cmath>
#include <complex>
#include <concepts>
#include <execution>
#include <format>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <optional>
#include <queue>
#include <ranges>
#include <stack>
#include <string>

//  From https://stackoverflow.com/a/37264642/6667035
#ifndef NDEBUG
#   define M_Assert(Expr, Msg) \
    M_Assert_Helper(#Expr, Expr, __FILE__, __LINE__, Msg)
#else
#   define M_Assert(Expr, Msg) ;
#endif

void M_Assert_Helper(const char* expr_str, bool expr, const char* file, int line, std::string msg)
{
    if (!expr)
    {
        std::cerr << "Assert failed:\t" << msg << "\n"
            << "Expected:\t" << expr_str << "\n"
            << "Source:\t\t" << file << ", line " << line << "\n";
        abort();
    }
}

struct recursive_print_fn
{
    template<std::ranges::input_range T>
    constexpr auto operator()(const T& input, const int level = 0) const
    {
        T output = input;
        std::cout << std::string(level, ' ') << "Level " << level << ":" << std::endl;
        std::ranges::transform(std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
            [&](auto&& x)
            {
                std::cout << std::string(level, ' ') << x << std::endl;
                return x;
            }
        );
        return output;
    }

    template<std::ranges::input_range T>
    requires (std::ranges::input_range<std::ranges::range_value_t<T>>)
    constexpr auto operator()(const T& input, const int level = 0) const
    {
        T output = input;
        std::cout << std::string(level, ' ') << "Level " << level << ":" << std::endl;
        std::ranges::transform(std::ranges::cbegin(input), std::ranges::cend(input), std::ranges::begin(output),
            [&](auto&& element)
            {
                return operator()(element, level + 1);
            }
        );
        return output;
    }
};

inline constexpr recursive_print_fn recursive_print;

struct RGB
{
    std::uint8_t channels[3];
};

struct RGB_DOUBLE
{
    double channels[3];
};

using GrayScale = std::uint8_t;

struct HSV
{
    double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
};

namespace TinyDIP
{
    //  Reference: https://stackoverflow.com/a/64287611/6667035
    template <typename T>
    struct is_complex : std::false_type {};

    template <typename T>
    struct is_complex<std::complex<T>> : std::true_type {};

    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.push_back(sizes), ...);
            image_data.resize(
                std::reduce(
                    std::ranges::cbegin(size),
                    std::ranges::cend(size),
                    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 (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)
        {
            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)
        {
            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)
        {
            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];
        }

        //  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];
        }

        std::vector<ElementT> const& getImageData() const noexcept { return image_data; }      //  expose the internal data

        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";
            }
        }

        //  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), ...);
        }

#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

    };

    template<typename T, typename ElementT>
    concept is_Image = std::is_same_v<T, Image<ElementT>>;

    //  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;
    }

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

    //  conv3 template function implementation
    template<typename ElementT>
    requires(std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value)
    constexpr static auto conv3(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        auto output = Image<ElementT>(
                                x.getSize(0) + y.getSize(0) - 1,
                                x.getSize(1) + y.getSize(1) - 1,
                                x.getSize(2) + y.getSize(2) - 1);

        for (std::size_t z2 = 0; z2 < y.getSize(2); ++z2)
        {
            for (std::size_t z1 = 0; z1 < x.getSize(2); ++z1)
            {
                for (std::size_t y2 = 0; y2 < y.getSize(1); ++y2)
                {
                    for (std::size_t y1 = 0; y1 < x.getSize(1); ++y1)
                    {
                        for (std::size_t x2 = 0; x2 < y.getSize(0); ++x2)
                        {
                            for (std::size_t x1 = 0; x1 < x.getSize(0); ++x1)
                            {
                                output.at(x1 + x2, y1 + y2, z1 + z2) = output.at(x1 + x2, y1 + y2, z1 + z2) + x.at(x1, y1, z1) * y.at(x2, y2, z2);
                            }
                        }
                    }
                }
            }
        }
        return output;
    }
}

int main()
{
    auto start = std::chrono::system_clock::now();
    auto image1 = TinyDIP::ones<double>(std::size_t{2}, std::size_t{2}, std::size_t{2});
    TinyDIP::conv3(image1, image1).print();
    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 EXIT_SUCCESS;
}

The output of the test code above:

1   2   1   
2   4   2   
1   2   1   

2   4   2   
4   8   4   
2   4   2   

1   2   1   
2   4   2   
1   2   1   


Computation finished at Sun Jul 21 08:39:45 2024
elapsed time: 0.000163883

Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  conv2 Template Function Implementation for Image in C++

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

  I am trying to implement conv3 template function in this post.

• Why a new review is being asked for?

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

Answer 1

In this line:

auto output = Image<ElementT>(

auto is not needed, because you only need to write the type once for a variable assignment that calls the constructor:

Image<ElementT> output(

In this line:

constexpr static auto conv3

auto needlessly obfuscates the signature. auto has its uses, especially for complex iterator types, etc., but it's often abused, and the cases in this code are both abuses.