# Gaussian Fisheye Image Generator Implementation in C++

This is a follow-up question for An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++ and Three dimensional gaussian image generator in C++. I am trying to make a Gaussian fisheye image generator in this post. The example output:

Image Input Output

The experimental implementation

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

namespace TinyDIP
{
//  gaussian_fisheye template function implementation
template<arithmetic ElementT, std::floating_point FloatingType = double>
constexpr static auto gaussian_fisheye(const Image<ElementT>& input, FloatingType D0)
{
if (input.getDimensionality()!=2)
{
throw std::runtime_error("Unsupported dimension!");
}

Image<ElementT> output(input.getWidth(), input.getHeight());
for (std::size_t y = 0; y < input.getHeight(); ++y)
{
for (std::size_t x = 0; x < input.getWidth(); ++x)
{
FloatingType distance_x = x - static_cast<FloatingType>(input.getWidth()) / 2.0;
FloatingType distance_y = y - static_cast<FloatingType>(input.getHeight()) / 2.0;
FloatingType distance = std::hypot(distance_x, distance_y);
FloatingType angle = std::atan2(distance_y, distance_x);
FloatingType weight = normalDistribution2D(std::fabs(distance_x), std::fabs(distance_y), D0) / normalDistribution2D(0.0, 0.0, D0);
FloatingType new_distance = distance * weight;
FloatingType new_distance_x = new_distance * std::cos(angle);
FloatingType new_distance_y = new_distance * std::sin(angle);
output.at(
static_cast<std::size_t>(new_distance_x + static_cast<FloatingType>(input.getWidth()) / 2.0),
static_cast<std::size_t>(new_distance_y + static_cast<FloatingType>(input.getHeight()) / 2.0)) =
input.at(x, y);
}
}
return output;
}

//  gaussian_fisheye template function implementation
template<typename ElementT, class FloatingType = double>
requires ((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, HSV>))
constexpr static auto gaussian_fisheye(const Image<ElementT>& input, FloatingType D0)
{
return apply_each(input, [&](auto&& planes) { return gaussian_fisheye(planes, D0); });
}
}

• normalDistribution2D function implementation (in file image_operations.h)

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 * std::numbers::pi * standard_deviation * standard_deviation);
}

• apply_each function implementation (in file image_operations.h)

template<class T = RGB, class F, class... Args>
requires (std::same_as<T, RGB>)
constexpr static auto apply_each(Image<T> input, F operation, Args&&... args)
{
return constructRGB(operation(getRplane(input), args...), operation(getGplane(input), args...), operation(getBplane(input), args...));
}

• 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!");
}
}

Image(std::vector<ElementT>&& 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";
}
}

//  Enable this function if ElementT = RGB
void print(std::string separator = "\t", std::ostream& os = std::cout) const
requires(std::same_as<ElementT, RGB>)
{
for (std::size_t y = 0; y < size[1]; ++y)
{
for (std::size_t x = 0; x < size[0]; ++x)
{
os << "( ";
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
os << +at(x, y).channels[channel_index] << separator;
}
os << ")" << separator;
}
os << "\n";
}
os << "\n";
return;
}

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

#endif


The usage of gaussian_fisheye function:

std::string file_path = "InputImages/1";
bmp1 = gaussian_fisheye(bmp1, 800.0);
TinyDIP::bmp_write("test", bmp1);


TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

Three dimensional gaussian image generator in C++

• What changes has been made in the code since last question?

I am trying to implement gaussian_fisheye template function in this post.

• Why a new review is being asked for?

Please review the implementation of gaussian_fisheye template function.

• Just like in your new question, here you should also iterate over the output image, and compute the corresponding location in the input. Otherwise either you write each output pixel multiple times, or you skip some output pixels. Commented Mar 27 at 1:41

# Unnecessary calculations

You calculate angle using std::atan2(), but later use std::cos() and std::sin() to convert it back into x and y coordinates. This is unnecessary, you already had those to begin with, you only need to scale them.

Furthermore, you don't need to use std::fabs() before passing the x and y position to normalDistribution2D(), since it will square those values anyway.

Finally, you don't need any of the static casts to FloatingType, only the final casts back to std::size_t. So:

FloatingType distance_x = x - input.getWidth() / 2.0;
FloatingType distance_y = y - input.getHeight() / 2.0;
FloatingType weight = normalDistribution2D(distance_x, distance_y, D0)
/ normalDistribution2D(0, 0, D0);
FloatingType new_distance_x = distance_x * weight;
FloatingType new_distance_y = distance_y * weight;
output.at(static_cast<std::size_t>(new_distance_x + input.getWidth() / 2.0),
static_cast<std::size_t>(new_distance_y + input.getHeight() / 2.0)) =
input.at(x, y);

• > normalDistribution2D(0, 0, …) is always equal to 1. Nope, its value still depends on the standard deviation Commented Mar 17 at 3:08
• Oh, I thought everything was inside the std::exp()... Commented Mar 17 at 8:38