Image pixelwise operation function with multiple inputs in C++

There are four images and each pixel value in these images are set to 4, 3, 2 and 1, respectively.

If we want to perform the element-wise calculation “Two times of img1 plus img2 then minus the result of img3 multiply img4, this task could be done with the following code:

  • pixelwiseOperation template function implementation: based on recursive_transform

    template<typename Op, class InputT, class... Args>
    constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
    {
        Image<InputT> output(
            recursive_transform<1>(
                (&)(auto&& element1, auto&&... elements) 
                    {
                        auto result = op(element1, elements...);
                        return static_cast<InputT>(std::clamp(
                            result,
                            static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                            static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                    },
                (input1.getImageData()),
                (inputs.getImageData())...),
            input1.getWidth(),
            input1.getHeight());
        return output;
    }
    
  • image_operations.h: The file contains pixelwiseOperation template function and other image processing functions

    /* Developed by Jimmy Hu */
    
    #ifndef ImageOperations_H
    #define ImageOperations_H
    
    #include <string>
    #include "base_types.h"
    #include "image.h"
    
    //  Reference: https://stackoverflow.com/a/26065433/6667035
    #ifndef M_PI
        #define M_PI 3.14159265358979323846
    #endif
    
    
    namespace TinyDIP
    {
        // Forward Declaration class Image
        template <typename ElementT>
        class Image;
    
        template<typename T>
        T normalDistribution1D(const T x, const T standard_deviation)
        {
            return std::exp(-x * x / (2 * standard_deviation * standard_deviation));
        }
    
        template<typename T>
        T normalDistribution2D(const T xlocation, const T ylocation, const T standard_deviation)
        {
            return std::exp(-(xlocation * xlocation + ylocation * ylocation) / (2 * standard_deviation * standard_deviation)) / (2 * M_PI * standard_deviation * standard_deviation);
        }
    
        template<class InputT1, class InputT2>
        constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
        {
            auto A = (v3-v2)-(v0-v1);
            auto B = (v0-v1)-A;
            auto C = v2-v0;
            auto D = v1;
            return D + frac * (C + frac * (B + frac * A));
        }
    
        template<class InputT = float, class ElementT>
        constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
        {
            auto x1 = cubicPolate( ndata(0), ndata(1), ndata(2), ndata(3), fracx );
            auto x2 = cubicPolate( ndata(4), ndata(5), ndata(6), ndata(7), fracx );
            auto x3 = cubicPolate( ndata(8), ndata(9), ndata(10), ndata(11), fracx );
            auto x4 = cubicPolate( ndata(12), ndata(13), ndata(14), ndata(15), fracx );
    
            return std::clamp(
                        cubicPolate( x1, x2, x3, x4, fracy ), 
                        static_cast<InputT>(std::numeric_limits<ElementT>::min()),
                        static_cast<InputT>(std::numeric_limits<ElementT>::max()));
        }
    
        template<class FloatingType = float, class ElementT>
        Image<ElementT> copyResizeBicubic(Image<ElementT>& image, size_t width, size_t height)
        {
            auto output = Image<ElementT>(width, height);
            //  get used to the C++ way of casting
            auto ratiox = static_cast<FloatingType>(image.getWidth()) / static_cast<FloatingType>(width);
            auto ratioy = static_cast<FloatingType>(image.getHeight()) / static_cast<FloatingType>(height);
    
            for (size_t y = 0; y < height; ++y)
            {
                for (size_t x = 0; x < width; ++x)
                {
                    FloatingType xMappingToOrigin = static_cast<FloatingType>(x) * ratiox;
                    FloatingType yMappingToOrigin = static_cast<FloatingType>(y) * ratioy;
                    FloatingType xMappingToOriginFloor = std::floor(xMappingToOrigin);
                    FloatingType yMappingToOriginFloor = std::floor(yMappingToOrigin);
                    FloatingType xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                    FloatingType yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
    
                    ElementT ndata(4 * 4);
                    for (int ndatay = -1; ndatay <= 2; ++ndatay)
                    {
                        for (int ndatax = -1; ndatax <= 2; ++ndatax)
                        {
                            ndata((ndatay + 1) * 4 + (ndatax + 1)) = image.at(
                                std::clamp(xMappingToOriginFloor + ndatax, static_cast<FloatingType>(0), image.getWidth() - static_cast<FloatingType>(1)), 
                                std::clamp(yMappingToOriginFloor + ndatay, static_cast<FloatingType>(0), image.getHeight() - static_cast<FloatingType>(1)));
                        }
    
                    }
                    output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
                }
            }
            return output;
        }
    
        //  multiple standard deviations
        template<class InputT>
        constexpr static Image<InputT> gaussianFigure2D(
            const size_t xsize, const size_t ysize, 
            const size_t centerx, const size_t centery,
            const InputT standard_deviation_x, const InputT standard_deviation_y)
        {
            auto output = TinyDIP::Image<InputT>(xsize, ysize);
            auto row_vector_x = TinyDIP::Image<InputT>(xsize, 1);
            for (size_t x = 0; x < xsize; ++x)
            {
                row_vector_x.at(x, 0) = normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x);
            }
    
            auto row_vector_y = TinyDIP::Image<InputT>(ysize, 1);
            for (size_t y = 0; y < ysize; ++y)
            {
                row_vector_y.at(y, 0) = normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
            }
    
            for (size_t y = 0; y < ysize; ++y)
            {
                for (size_t x = 0; x < xsize; ++x)
                {
                    output.at(x, y) = row_vector_x.at(x, 0) * row_vector_y.at(y, 0);
                }
            }
            return output;
        }
    
        //  single standard deviation
        template<class InputT>
        constexpr static Image<InputT> gaussianFigure2D(
            const size_t xsize, const size_t ysize,
            const size_t centerx, const size_t centery,
            const InputT standard_deviation)
        {
            return gaussianFigure2D(xsize, ysize, centerx, centery, standard_deviation, standard_deviation);
        }
    
        template<typename Op, class InputT, class... Args>
        constexpr static Image<InputT> pixelwiseOperation(Op op, const Image<InputT>& input1, const Args&... inputs)
        {
            Image<InputT> output(
                recursive_transform<1>(
                    (&)(auto&& element1, auto&&... elements) 
                        {
                            auto result = op(element1, elements...);
                            return static_cast<InputT>(std::clamp(
                                result,
                                static_cast<decltype(result)>(std::numeric_limits<InputT>::min()),
                                static_cast<decltype(result)>(std::numeric_limits<InputT>::max())));
                        },
                    (input1.getImageData()),
                    (inputs.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 TinyDIP::pixelwiseOperation(std::plus<>{}, input1, plus(inputs...));
        }
    
        template<class InputT>
        constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2)
        {
            return TinyDIP::pixelwiseOperation(std::minus<>{}, input1, input2);
        }
    
        template<class InputT = RGB>
        requires (std::same_as<InputT, RGB>)
        constexpr static Image<InputT> subtract(Image<InputT>& input1, Image<InputT>& input2)
        {
            assert(input1.getWidth() == input2.getWidth());
            assert(input1.getHeight() == input2.getHeight());
            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;
        }
    }
    
    #endif
    
  • image.h: The file contains the definition of Image class.

    /* Developed by Jimmy Hu */
    
    #ifndef Image_H
    #define Image_H
    
    #include <algorithm>
    #include <array>
    #include <cassert>
    #include <chrono>
    #include <complex>
    #include <concepts>
    #include <functional>
    #include <iostream>
    #include <iterator>
    #include <list>
    #include <numeric>
    #include <string>
    #include <type_traits>
    #include <variant>
    #include <vector>
    #include "image_operations.h"
    
    namespace TinyDIP
    {
        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)
            {
                assert(input.size() == newWidth * newHeight);
                this->image_data = input;   //  Deep copy
            }
    
            Image(const std::vector<std::vector<ElementT>>& input)
            {
                this->height = input.size();
                this->width = input(0).size();
    
                for (auto& rows : input)
                {
                    this->image_data.insert(this->image_data.end(), std::begin(input), std::end(input));    //  flatten
                }
                return;
            }
    
            constexpr ElementT& at(const unsigned int x, const unsigned int y)
            { 
                checkBoundary(x, y);
                return this->image_data(y * width + x);
            }
    
            constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
            {
                checkBoundary(x, y);
                return this->image_data(y * width + x);
            }
    
            constexpr std::size_t getWidth() const
            {
                return this->width;
            }
    
            constexpr std::size_t getHeight() const
            {
                return this->height;
            }
    
            std::vector<ElementT> const& getImageData() const { return this->image_data; }      //  expose the internal data
    
            void print()
            {
                for (std::size_t y = 0; y < this->height; ++y)
                {
                    for (std::size_t x = 0; x < this->width; ++x)
                    {
                        //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                        std::cout << +this->at(x, y) << "t";
                    }
                    std::cout << "n";
                }
                std::cout << "n";
                return;
            }
    
            //  Enable this function if ElementT = RGB
            void print() requires(std::same_as<ElementT, RGB>)
            {
                for (std::size_t y = 0; y < this->height; ++y)
                {
                    for (std::size_t x = 0; x < this->width; ++x)
                    {
                        std::cout << "( ";
                        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
                            std::cout << +this->at(x, y).channels(channel_index) << "t";
                        }
                        std::cout << ")t";
                    }
                    std::cout << "n";
                }
                std::cout << "n";
                return;
            }
    
            Image<ElementT>& operator+=(const Image<ElementT>& rhs)
            {
                assert(rhs.width == this->width);
                assert(rhs.height == this->height);
                std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                       image_data.begin(), std::plus<>{});
                return *this;
            }
    
            Image<ElementT>& operator-=(const Image<ElementT>& rhs)
            {
                assert(rhs.width == this->width);
                assert(rhs.height == this->height);
                std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                       image_data.begin(), std::minus<>{});
                return *this;
            }
    
            Image<ElementT>& operator*=(const Image<ElementT>& rhs)
            {
                assert(rhs.width == this->width);
                assert(rhs.height == this->height);
                std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                       image_data.begin(), std::multiplies<>{});
                return *this;
            }
    
            Image<ElementT>& operator/=(const Image<ElementT>& rhs)
            {
                assert(rhs.width == this->width);
                assert(rhs.height == this->height);
                std::transform(image_data.cbegin(), image_data.cend(), rhs.image_data.cbegin(),
                       image_data.begin(), std::divides<>{});
                return *this;
            }
    
            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:
            size_t width;
            size_t height;
            std::vector<ElementT> image_data;
    
            void checkBoundary(const size_t x, const size_t y)
            {
                assert(x < width);
                assert(y < height);
            }
        };
    }
    #endif
    
  • base_types.h: The base types declaration

    /* Developed by Jimmy Hu */
    
    #ifndef BASE_H
    #define BASE_H
    
    #define _USE_MATH_DEFINES
    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include <string>
    #include <utility>
    
    constexpr int MAX_PATH = 256;
    #define FILE_ROOT_PATH "./"
    
    using BYTE = unsigned char;
    
    struct RGB
    {
        unsigned char channels(3);
    };
    
    using GrayScale = BYTE;
    
    struct HSV
    {
        double channels(3);    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
    };
    
    struct BMPIMAGE
    {
        char FILENAME(MAX_PATH);
    
        unsigned int XSIZE;
        unsigned int YSIZE;
        unsigned char FILLINGBYTE;
        unsigned char *IMAGE_DATA;
    };
    #endif
    
  • basic_functions.h: The file contains the definition of recursive_transform

    /* Developed by Jimmy Hu */
    
    #ifndef BasicFunctions_H
    #define BasicFunctions_H
    
    #include <algorithm>
    #include <array>
    #include <cassert>
    #include <chrono>
    #include <complex>
    #include <concepts>
    #include <deque>
    #include <execution>
    #include <exception>
    #include <functional>
    #include <iostream>
    #include <iterator>
    #include <list>
    #include <map>
    #include <mutex>
    #include <numeric>
    #include <optional>
    #include <ranges>
    #include <stdexcept>
    #include <string>
    #include <tuple>
    #include <type_traits>
    #include <utility>
    #include <variant>
    #include <vector>
    
    namespace TinyDIP
    {
        //  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)
            {
                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...);
            }
        }
    }
    
    #endif
    
  • All suggestions are welcome.