7
\$\begingroup\$

This is a follow-up question for An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. I am trying to implement conv2 template function like Matlab's conv2.

The experimental implementation

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

    namespace TinyDIP
    {
        //  conv2 template function implementation
        template<typename ElementT>
        requires(std::floating_point<ElementT> || std::integral<ElementT>)
        constexpr auto conv2(const Image<ElementT>& x, const Image<ElementT>& y)
        {
            auto output = Image<ElementT>(x.getWidth() + y.getWidth() - 1, x.getHeight() + y.getHeight() - 1);
            for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1)
            {
                for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1)
                {
                    for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2)
                    {
                        for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2)
                        {
                            output.at(x1 + x2, y1 + y2) = output.at(x1 + x2, y1 + y2) + x.at(x1, y1) * y.at(x2, y2);
                        }
                    }
                }
            }
            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.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::ranges::input_range Range>
            requires(std::same_as<std::ranges::range_value_t<Range>, ElementT>)
            Image(Range&& 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];
            }
    
            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];
            }
    
            constexpr auto getSize() noexcept
            {
                return size;
            }
    
            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";
                }
            }
    
            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>>;
    }
    

Full Testing Code

The full testing code:

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

#include <algorithm>
#include <cassert>
#include <cctype>
#include <chrono>
#include <cmath>
#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;

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.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::ranges::input_range Range>
        requires(std::same_as<std::ranges::range_value_t<Range>, ElementT>)
        Image(Range&& 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];
        }

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

        constexpr auto getSize() noexcept
        {
            return size;
        }

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

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

    //  conv2 template function implementation
    template<typename ElementT>
    requires(std::floating_point<ElementT> || std::integral<ElementT>)
    constexpr auto conv2(const Image<ElementT>& x, const Image<ElementT>& y)
    {
        auto output = Image<ElementT>(x.getWidth() + y.getWidth() - 1, x.getHeight() + y.getHeight() - 1);
        for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1)
        {
            for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1)
            {
                for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2)
                {
                    for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2)
                    {
                        output.at(x1 + x2, y1 + y2) = output.at(x1 + x2, y1 + y2) + x.at(x1, y1) * y.at(x2, y2);
                    }
                }
            }
        }
        return output;
    }
}

int main()
{
    auto start = std::chrono::system_clock::now();

    TinyDIP::Image<double> I1(3, 3);
    for(std::size_t y = 0; y < I1.getHeight(); ++y)
    {
        for(std::size_t x = 0; x < I1.getWidth(); ++x)
        {
            I1.at(x, y) = x + 1;
        }
    }
    I1.print();
    std::cout << "convolution with \n";
    I1.print();
    std::cout << "is \n";
    TinyDIP::conv2(I1, I1).print();
    std::cout << "========================\n";
    std::vector<double> data = {1, 2};
    TinyDIP::Image<double> I2(data, 1, 2);
    I1.print();
    std::cout << "convolution with \n";
    I2.print();
    std::cout << "is \n";
    conv2(I1, I2).print();
    std::cout << "========================\n";
    data = {1, 2, 3, 4};
    TinyDIP::Image<double> I3(data, 4, 1);
    I1.print();
    std::cout << "convolution with \n";
    I3.print();
    std::cout << "is \n";
    conv2(I1, I3).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 0;
}

The output of the test code above:

1   2   3   
1   2   3   
1   2   3   

convolution with 
1   2   3   
1   2   3   
1   2   3   

is 
1   4   10  12  9   
2   8   20  24  18  
3   12  30  36  27  
2   8   20  24  18  
1   4   10  12  9   

========================
1   2   3   
1   2   3   
1   2   3   

convolution with 
1   
2   

is 
1   2   3   
3   6   9   
3   6   9   
2   4   6   

========================
1   2   3   
1   2   3   
1   2   3   

convolution with 
1   2   3   4   

is 
1   4   10  16  17  12  
1   4   10  16  17  12  
1   4   10  16  17  12  

Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

\$\endgroup\$

3 Answers 3

6
\$\begingroup\$

Algorithm

I like the algorithm you use. Usually people will implement convolution by doing all computations for one output pixel in the inner two loops. I honestly don't know which is more efficient, but the code for your algorithm is very clean and simple.

However, in image processing one usually wants an output image that is the same size as the input image, this implementation always returns something larger. And this implementation has a strong boundary effect because the output near the image boundary mixes in the assumed zeros outside the image domain (they are there implicitly). Whether these conditions are OK or not depends on your application.

API

Your function takes two images with the same type as input:

template<typename ElementT>
requires(std::floating_point<ElementT> || std::integral<ElementT>)
constexpr auto conv2(const Image<ElementT>& x, const Image<ElementT>& y)

This keeps it simple, but it also makes it impossible to implement a smoothing filter on a uint8 image. Because if the image is uint8, then the kernel must be uint8 as well. The most trivial smoothing filter then would be a 3x3 array of 1s. A convolution with such an array would produce values that exceed the range of uint8. A smoothing filter is normalized, and so needs floating-point values. This function therefore forces the user to convert their image to floating-point before computing a simple smoothing filter.

Ideally, x and y can have different types. This will make your code more complex (there'll be explicit casting, rounding for integer output types, and whatnot). A simpler solution would rename y to h or kernel and force it to be a floating-point type.

Your function I guess is a basic building block for filters. Higher-level functions would allow various forms of convolution, add padding with something other than zeros for boundary extension, cropping of the output to the input size, etc. A separable convolution would be implemented by calling conv2 twice.

Efficiency

output.at(x1 + x2, y1 + y2) = output.at(x1 + x2, y1 + y2) + x.at(x1, y1) * y.at(x2, y2);

Here you compute the indexing into the output twice. Prefer

output.at(x1 + x2, y1 + y2) += x.at(x1, y1) * y.at(x2, y2);

Furthermore, your .at() does bounds checking. This is a huge performance hit for this algorithm. By construction, you know that you will never index out of bounds in this algorithm. So implement an unchecked version of .at(), a simple computation y * width + x. The faster this computation can be, the faster your convolution will be. The computation of the indexing is the most expensive part here. I would even suggest moving the y * width out of the loop over x (see the next point).

Unless you store the images column-wise (which no image processing package does except MATLAB, don't do that!), you should order your loops so that pixels are accessed in memory order. Your current order,

for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1)
    for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1)
        for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2)
            for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2)

loops over y1 inside the loop over x1. This means that in.at(x1, y1) will be jumping all over the memory. This means most reads will be cache misses, and you'll read data very slowly. Make the loop over x1 the inner one.

And because none of these four loops depend on each other, we can reorder all four any way we please. I would suggest putting both y loops first:

for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1)
    for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2)
        for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1)
            for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2)

Now the inner two loops both run over rows in the two images, which are contiguous in memory (the stride is 1). So:

for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1) {
    auto* x_row = &(x.at(0, y1));
    for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2) {
        auto* y_row = &(y.at(0, y2));
        auto* out_row = &(out.at(0, y1 + y2));
        for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1) {
            for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2) {
                out_row[x1 + x2] += x_row[x1] * y_row[x2];
            }
        }
    }
}

This can probably be made even more efficient, but this is a good compromise. Without having tested this, I'd be willing to bet it is significantly faster than the original code. It is also much less clean than the original code... :)

\$\endgroup\$
0
2
\$\begingroup\$

Think about the names you use. The function's input images are called x and y, which makes the subsequent x1, y1, x2 and y2 variables confusing.

If we called the inputs a and b, we could call their corresponding coordinates ax, ay, bx and by for code that's more comprehensible and less prone to mistakes of cognition.


This convolution seems only to work with images of two dimensions, yet Image supports an arbitrary number of dimensions. If we're not going to support arbitrary-dimension convolutions, perhaps we need more constraints so that only 2-dimensional images can be passed (or perhaps a runtime check to throw an exception, as it appears the dimensionality is variable).

\$\endgroup\$
2
\$\begingroup\$

I see several problems with the overall design. I'll refer to issues that @CrisLuengo didn't cover.

Interface: I am not entirely sure what exactly the aim of this conv2 function is, but the most typical in computer vision is convolving an image with a kernel. The image's elements are typically uint8_t and sometimes come with 3 channels. The kernel is a small image of floating-point type elements or integral elements, with the assumption that the output is divided by the total sum of the kernel's elements. It is performed this way because uint8_t easily overflows.

The problem is that this conv2 method doesn't allow one to perform the convolution unless one first converts the data into a float/double element type, which is less than ideal in terms of performance.

If you want to perform convolutions with both images being large, you can do so via FFT2 transformation rather than a direct computation. This method takes O(n^2 ln n) operations instead of O(n^4).

Efficiency: The code's efficiency is somewhat problematic. In modern computers, many performance issues occur due to relatively slow memory loads from RAM. It's crucial to have code more cache-friendly, which you can achieve if you use data that was recently used rather than one used a long while ago.

Presuming that you have the scenario where an image is convolved with a kernel, the kernel is small, so it all will likely be permanently in the cache, while the image is large, so just going over a single row might result in the cache losing the previous line. For this reason, it is usually more desirable to compute the convolutions in subimages of some fixed size, say, 64x64, as then one more frequently refers to the same elements, leading to improved efficiency.

The last problem is that you don't use any SIMD operations. Incorporating them is very troublesome, but it is one of the significant performance improvements you can add. Unless you utilize it, your code's performance won't reach anywhere near what's available in image processing libraries available on the internet.

Edit: @TobySpeight asked why SIMD won't happen automatically. The compiler is sometimes capable of implementing it, but frequently, it is unable to, as from its perspective, it might potentially change the program's behavior.

Quite often, the compiler assumes that the input and output array can intersect, and thus, adding such modifications as SIMD or other reordering operations would lead to behavior change. In this particular case, it should be feasible, but assuming that the compiler will succeed is still unreliable. Furthermore, in the case of operations over floating points, automatic SIMD and other improvements will likely require commutativity, which isn't the case. The compiler can produce such behavior-changing alterations with the fast-math flag, but this may cause damage to the program elsewhere in the code where commutativity cannot be assumed.

Besides, when writing SIMD operations, one tends to write differently for the best performance, which would be hard for the compiler to replicate in such a quadruple for scenario.

\$\endgroup\$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.