Image pixelwise operation function with multiple inputs in C++

Question

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)
    {
        assert(input1.getWidth() == input2.getWidth());
        assert(input1.getHeight() == input2.getHeight());
        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.

I’ve commented this in an answer to you before, include only the include files you actually need in that source file. Your basic_functions.h includes everything that you never need. There are more lines of include statements than actual code. There is a reason there’s not just one <std> header. It is very inefficient to include everything. — Cris Luengo, Commented Jul 22, 2021 at 3:05

JDługosz · Accepted Answer · 2021-07-21 00:28:23Z

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

Don't use #define, and since you marked this as C++20 use the constants supplied in the standard library.

constexpr int MAX_PATH = 256;
#define FILE_ROOT_PATH "./"

Use the path type, not fixed-size arrays of characters, for filename manipulation. Don't use #define.

using BYTE = unsigned char;

There is a std::byte type now. If that is not what you want, you ought to use uint_8 to avoid confusion and be consistent with standard code.

struct RGB
{
    unsigned char channels[3];
};

Come on, you just defined BYTE on the previous line! Why don't you use it?

struct BMPIMAGE
{
    char FILENAME[MAX_PATH];

    unsigned int XSIZE;
    unsigned int YSIZE;
    unsigned char FILLINGBYTE;
    unsigned char *IMAGE_DATA;
};

Use std::filesystem::path rather than a fixed array of characters. Use int32_t or whatever, rather than the implementation-dependant unsigned int. Prefer signed types unless you really need that one more bit for the range, or are doing bitwise operations.

You defined BYTE and probably mean that here, but don't use it. Use uint8_t or std::byte for these.

You might consider using a 2-D point class rather than separate height and width fields. You may find it is common to have x and y things always being used together, for positions and sizes.

Stack Exchange Network

Image pixelwise operation function with multiple inputs in C++

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
c++
image
template
lambda
c++20
or ask your own question.

Linked

Hot Network Questions

Image pixelwise operation function with multiple inputs in C++

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

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

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
c++
image
template
lambda
c++20
or ask your own question.