This is a follow-up question for [Tests for the operators of image template class in C++](https://codereview.stackexchange.com/q/263729/231235) and [A recursive_transform template function for the multiple parameters cases in C++](https://codereview.stackexchange.com/q/263956/231235). I appreciated [G. Sliepen's answer](https://codereview.stackexchange.com/a/263988/231235). 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.
```C++
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:
```C++
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`
```C++
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
```C++
/* 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.
```C++
/* 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
```C++
/* 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`
```C++
/* 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**
```C++
/* 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++](https://codereview.stackexchange.com/q/263729/231235) and
[A recursive_transform template function for the multiple parameters cases in C++](https://codereview.stackexchange.com/q/263956/231235)
- 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.