Return to Revisions

2 of 3

Update question description

edited Jul 19, 2021 at 23:47

5.5k
2
10
40

Image pixelwise operation function with multiple inputs in C++

This is a follow-up question for Tests for the operators of image template class in C++ and A recursive_transform template function for the multiple parameters cases in C++. I appreciated G. Sliepen's answer. I am attempting to extend the mentioned element-wise operations in TinyDIP::Image image class. In other words, not only +, -, * and / but also other customized calculations can be specified easily with the implemented pixelwiseOperation template function here. For example:

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

auto img1 = TinyDIP::Image<GrayScale>(10, 10, 4);
auto img2 = TinyDIP::Image<GrayScale>(10, 10, 3);
auto img3 = TinyDIP::Image<GrayScale>(10, 10, 2);
auto img4 = TinyDIP::Image<GrayScale>(10, 10, 1);

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:

auto output = TinyDIP::pixelwiseOperation(
    [](auto&& pixel_in_img1, auto&& pixel_in_img2, auto&& pixel_in_img3, auto&& pixel_in_img4)
    {
        return 2 * pixel_in_img1 + pixel_in_img2 - pixel_in_img3 * pixel_in_img4;
    },
    img1, img2, img3, img4
);

The result can be printed with output.print();:

9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9
9       9       9       9       9       9       9       9       9       9

The experimental implementation

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

The testing code

/* Developed by Jimmy Hu */

#include "image.h"
#include "basic_functions.h"

int main()
{
    auto img1 = TinyDIP::Image<GrayScale>(10, 10, 4);
    auto img2 = TinyDIP::Image<GrayScale>(10, 10, 3);
    auto img3 = TinyDIP::Image<GrayScale>(10, 10, 2);
    auto img4 = TinyDIP::Image<GrayScale>(10, 10, 1);
    auto output = TinyDIP::pixelwiseOperation(
        [](auto&& pixel_in_img1, auto&& pixel_in_img2, auto&& pixel_in_img3, auto&& pixel_in_img4)
        {
            return 2 * pixel_in_img1 + pixel_in_img2 - pixel_in_img3 * pixel_in_img4;
        },
        img1, img2, img3, img4
    );
    output.print();
    
    return 0;
}

Test platform

MacOS: g++-11 (Homebrew GCC 11.1.0_1) 11.1.0

All suggestions are welcome.

The summary information:

Which question it is a follow-up to?

Tests for the operators of image template class in C++ and

A recursive_transform template function for the multiple parameters cases in C++
What changes has been made in the code since last question?

I am attempting to extend the mentioned element-wise operations in this post.
Why a new review is being asked for?

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

c++image template lambda c++20

asked Jul 18, 2021 at 18:48

JimmyHu

5.5k
2
10
40