Two dimensional bicubic interpolation implementation in C++

Question

This is a follow-up question for Two dimensional bicubic interpolation implementation in C and A recursive_transform Template Function with Unwrap Level for Various Type Arbitrary Nested Iterable Implementation in C++. Besides the C version code, I am attempting to make a C++ version bicubic interpolation function bicubicInterpolation which can be applied on two dimensional nested vectors std::vector<std::vector<>> structure.

Example input matrix:

1       1       1
1       100     1
1       1       1

Output matrix (bicubic interpolation result from the input matrix above):

1       1       1       1       1       1       1       1       1       1       1       1
1       9       19      27      30      27      19      9       1       0       0       0
1       19      39      56      62      56      39      19      1       0       0       0
1       27      56      79      89      79      56      27      1       0       0       0
1       30      62      89      100     89      62      30      1       0       0       0
1       27      56      79      89      79      56      27      1       0       0       0
1       19      39      56      62      56      39      19      1       0       0       0
1       9       19      27      30      27      19      9       1       0       0       0
1       1       1       1       1       1       1       1       1       1       1       1
1       0       0       0       0       0       0       0       1       2       2       1
1       0       0       0       0       0       0       0       1       2       2       1
1       0       0       0       0       0       0       0       1       1       1       1

The experimental implementation

namespace: TinyDIP

bicubicInterpolation function implementation:

constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY)
{
    auto output = Image<ElementT>(newSizeX, newSizeY);
    auto ratiox = (float)this->getSizeX() / (float)newSizeX;
    auto ratioy = (float)this->getSizeY() / (float)newSizeY;

    for (size_t y = 0; y < newSizeY; y++)
    {
        for (size_t x = 0; x < newSizeX; x++)
        {
            float xMappingToOrigin = (float)x * ratiox;
            float yMappingToOrigin = (float)y * ratioy;
            float xMappingToOriginFloor = floor(xMappingToOrigin);
            float yMappingToOriginFloor = floor(yMappingToOrigin);
            float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
            float 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)] = this->get(
                        clip(xMappingToOriginFloor + ndatax, 0, this->getSizeX() - 1), 
                        clip(yMappingToOriginFloor + ndatay, 0, this->getSizeY() - 1));
                }

            }
            output.set(x, y, bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac));
        }
    }
    return output;
}

Helper functions for bicubicInterpolation function:

template<class InputT>
constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
}

template<class InputT1, class InputT2>
constexpr 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 InputT1, class InputT2, class InputT3>
constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
{
    if (input < lowerbound)
    {
        return static_cast<InputT1>(lowerbound);
    }
    if (input > upperbound)
    {
        return static_cast<InputT1>(upperbound);
    }
    return input;
}

Image template class implementation (image.h):

/* Develop by Jimmy Hu */

#ifndef Image_H
#define Image_H

#include <algorithm>
#include <array>
#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 "basic_functions.h"

namespace TinyDIP
{
    template <typename ElementT>
    class Image
    {
    public:
        Image()
        {
        }

        Image(const int newWidth, const int newHeight)
        {
            this->image_data.resize(newHeight);
            for (size_t i = 0; i < newHeight; ++i) {
                this->image_data[i].resize(newWidth);
            }
            this->image_data = recursive_transform<2>(this->image_data, [](ElementT element) { return ElementT{}; });
            return;
        }

        Image(const int newWidth, const int newHeight, const ElementT initVal)
        {
            this->image_data.resize(newHeight);
            for (size_t i = 0; i < newHeight; ++i) {
                this->image_data[i].resize(newWidth);
            }
            this->image_data = recursive_transform<2>(this->image_data, [initVal](ElementT element) { return initVal; });
            return;
        }

        Image(const std::vector<std::vector<ElementT>>& input)
        {
            this->image_data = recursive_transform<2>(input, [](ElementT element) {return element; } ); //  Deep copy
            return;
        }

        template<class OutputT>
        constexpr auto cast()
        {
            return this->transform([](ElementT element) { return static_cast<OutputT>(element); });
        }

        constexpr auto get(const unsigned int locationx, const unsigned int locationy)
        {
            return this->image_data[locationy][locationx];
        }

        constexpr auto set(const unsigned int locationx, const unsigned int locationy, const ElementT& element)
        {
            this->image_data[locationy][locationx] = element;
            return *this;
        }

        template<class InputT>
        constexpr auto set(const unsigned int locationx, const unsigned int locationy, const InputT& element)
        {
            this->image_data[locationy][locationx] = static_cast<ElementT>(element);
            return *this;
        }

        constexpr auto getSizeX()
        {
            return this->image_data[0].size();
        }

        constexpr auto getSizeY()
        {
            return this->image_data.size();
        }

        constexpr auto getData()
        {
            return this->transform([](ElementT element) { return element; });   //  Deep copy
        }

        void print()
        {
            for (auto& row_element : this->toString())
            {
                for (auto& element : row_element)
                {
                    std::cout << element << "\t";
                }
                std::cout << "\n";
            }
            std::cout << "\n";
            return;
        }

        constexpr auto toString()
        {
            return this->transform([](ElementT element) { return std::to_string(element); });
        }

        constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY)
        {
            auto output = Image<ElementT>(newSizeX, newSizeY);
            auto ratiox = (float)this->getSizeX() / (float)newSizeX;
            auto ratioy = (float)this->getSizeY() / (float)newSizeY;

            for (size_t y = 0; y < newSizeY; y++)
            {
                for (size_t x = 0; x < newSizeX; x++)
                {
                    float xMappingToOrigin = (float)x * ratiox;
                    float yMappingToOrigin = (float)y * ratioy;
                    float xMappingToOriginFloor = floor(xMappingToOrigin);
                    float yMappingToOriginFloor = floor(yMappingToOrigin);
                    float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                    float 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)] = this->get(
                                clip(xMappingToOriginFloor + ndatax, 0, this->getSizeX() - 1), 
                                clip(yMappingToOriginFloor + ndatay, 0, this->getSizeY() - 1));
                        }

                    }
                    output.set(x, y, bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac));
                }
            }
            return output;
        }

        Image<ElementT>& operator=(Image<ElementT> const& input)  //  Copy Assign
        {
            this->image_data = input.getData();
            return *this;
        }

        Image<ElementT>& operator=(Image<ElementT>&& other)       //  Move Assign
        {
            this->image_data = std::move(other.image_data);
            std::cout << "move assigned\n";
            return *this;
        }

        Image(const Image<ElementT> &input)                      //  Copy Constructor
        {
            this->image_data = input.getData();
        }

        /*    Move Constructor
         */
        Image(Image<ElementT> &&input) : image_data(std::move(input.image_data))
        {
        }

    private:
        std::vector<std::vector<ElementT>> image_data;

        template<class F>
        constexpr auto transform(const F& f)
        {
            return recursive_transform<2>(this->image_data, f);
        }

        template<class InputT>
        constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
        }

        template<class InputT1, class InputT2>
        constexpr 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 InputT1, class InputT2, class InputT3>
        constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
        {
            if (input < lowerbound)
            {
                return static_cast<InputT1>(lowerbound);
            }
            if (input > upperbound)
            {
                return static_cast<InputT1>(upperbound);
            }
            return input;
        }
    };
}

#endif

base_types.h: The base types

/* Develop by Jimmy Hu */

#ifndef BASE_H
#define BASE_H

#include <cmath>
#include <cstdbool>
#include <cstdio>
#include <cstdlib>
#include <string>

#define MAX_PATH 256
#define FILE_ROOT_PATH "./"

#define True true
#define False false

typedef unsigned char BYTE;

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

typedef BYTE GrayScale;

typedef struct HSV
{
    long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
}HSV;

#endif

basic_functions.h: The basic functions

/* Develop 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
{
    template<typename T>
    concept is_back_inserterable = requires(T x)
    {
        std::back_inserter(x);
    };

    template<typename T>
    concept is_inserterable = requires(T x)
    {
        std::inserter(x, std::ranges::end(x));
    };

    //  recursive_invoke_result_t implementation
    template<typename, typename>
    struct recursive_invoke_result { };

    template<typename T, std::invocable<T> F>
    struct recursive_invoke_result<F, T> { using type = std::invoke_result_t<F, T>; };

    template<typename F, template<typename...> typename Container, typename... Ts>
    requires (
        !std::invocable<F, Container<Ts...>>&&
        std::ranges::input_range<Container<Ts...>>&&
        requires { typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type; })
        struct recursive_invoke_result<F, Container<Ts...>>
    {
        using type = Container<typename recursive_invoke_result<F, std::ranges::range_value_t<Container<Ts...>>>::type>;
    };

    template<typename F, typename T>
    using recursive_invoke_result_t = typename recursive_invoke_result<F, T>::type;

    //  recursive_transform implementation (the version with unwrap_level)
    template<std::size_t unwrap_level = 1, class T, class F>
    constexpr auto recursive_transform(const T& input, const F& f)
    {
        if constexpr (unwrap_level > 0)
        {
            recursive_invoke_result_t<F, T> output{};
            std::ranges::transform(
                std::ranges::cbegin(input),
                std::ranges::cend(input),
                std::inserter(output, std::ranges::end(output)),
                [&f](auto&& element) { return recursive_transform<unwrap_level - 1>(element, f); }
            );
            return output;
        }
        else
        {
            return f(input);
        }
    }

    template<std::size_t dim, class T>
    constexpr auto n_dim_vector_generator(T input, std::size_t times)
    {
        if constexpr (dim == 0)
        {
            return input;
        }
        else
        {
            auto element = n_dim_vector_generator<dim - 1>(input, times);
            std::vector<decltype(element)> output(times, element);
            return output;
        }
    }

    template<std::size_t dim, std::size_t times, class T>
    constexpr auto n_dim_array_generator(T input)
    {
        if constexpr (dim == 0)
        {
            return input;
        }
        else
        {
            auto element = n_dim_array_generator<dim - 1, times>(input);
            std::array<decltype(element), times> output;
            std::fill(std::ranges::begin(output), std::ranges::end(output), element);
            return output;
        }
    }

    template<std::size_t dim, class T>
    constexpr auto n_dim_deque_generator(T input, std::size_t times)
    {
        if constexpr (dim == 0)
        {
            return input;
        }
        else
        {
            auto element = n_dim_deque_generator<dim - 1>(input, times);
            std::deque<decltype(element)> output(times, element);
            return output;
        }
    }

    template<std::size_t dim, class T>
    constexpr auto n_dim_list_generator(T input, std::size_t times)
    {
        if constexpr (dim == 0)
        {
            return input;
        }
        else
        {
            auto element = n_dim_list_generator<dim - 1>(input, times);
            std::list<decltype(element)> output(times, element);
            return output;
        }
    }

    template<std::size_t dim, template<class...> class Container = std::vector, class T>
    constexpr auto n_dim_container_generator(T input, std::size_t times)
    {
        if constexpr (dim == 0)
        {
            return input;
        }
        else
        {
            return Container(times, n_dim_container_generator<dim - 1, Container, T>(input, times));
        }
    }
}
#endif

The full testing code

The grayscale type data has been tested here.

/* Develop by Jimmy Hu */

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

void bicubicInterpolationTest();

int main()
{
    bicubicInterpolationTest();
    return 0;
}

void bicubicInterpolationTest()
{
    TinyDIP::Image<GrayScale> image1(3, 3, 1);
    std::cout << "Width: " + std::to_string(image1.getSizeX()) + "\n";
    std::cout << "Height: " + std::to_string(image1.getSizeY()) + "\n";
    image1 = image1.set(1, 1, 100);
    image1.print();

    auto image2 = image1.bicubicInterpolation(12, 12);
    image2.print();
}

All suggestions are welcome.

The summary information:

Which question it is a follow-up to?

Two dimensional bicubic interpolation implementation in C and

A recursive_transform Template Function with Unwrap Level for Various Type Arbitrary Nested Iterable Implementation in C++.
What changes has been made in the code since last question?

I am attempting to make a C++ version bicubic interpolation function bicubicInterpolation which can be applied on two dimensional nested vectors std::vector<std::vector<>> structure.
Why a new review is being asked for?

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

Can't help but notice that the bug from last time (filling 9 out of 16 entries in ndata) is still in here — user555045, Commented Jun 21, 2021 at 18:31
@harold I apologize for the missing point you mentioned. The question has been updated. — JimmyHu, Commented Jun 22, 2021 at 0:08

user673679 · Accepted Answer · 2021-06-22 10:27:27Z

I like auto for lots of things, but not as a return type unless it's actually necessary.

For example, to understand what cast returns it's necessary to dig through several template functions, layers of meta-templates, and then understand what std::invoke_result_t does in this particular case.

Please just write out the return type.

    template<class InputT1, class InputT2, class InputT3>
    constexpr auto clip(const InputT1& input, const InputT2& lowerbound, const InputT3& upperbound)
    {
        if (input < lowerbound)
        {
            return static_cast<InputT1>(lowerbound);
        }
        if (input > upperbound)
        {
            return static_cast<InputT1>(upperbound);
        }
        return input;
    }

The standard library already provides std::clamp so this isn't necessary.

Note that providing a single template argument is safer than separate arguments for all the inputs. Conversions (and thus comparisons) between different types may not be well-defined since types have different ranges.

#define True true
#define False false

Why? C++ uses true and false.

    std::vector<std::vector<ElementT>> image_data;

It's faster (and usually easier) to use a one-dimensional vector: std::vector<ElementT> image_data; containing width * height elements, and calculate indices when necessary as y * width + x.

So this:

    Image(const int newWidth, const int newHeight)
    {
        this->image_data.resize(newHeight);
        for (size_t i = 0; i < newHeight; ++i) {
            this->image_data[i].resize(newWidth);
        }
        this->image_data = recursive_transform<2>(this->image_data, [](ElementT element) { return ElementT{}; });
        return;
    }

Should simply be:

    Image(const unsigned int width, const unsigned int height):
        m_width(width),
        m_height(height),
        m_image_data(width * height) { }

Note that the std::vector constructor (and resize function) will "value-initialize" elements, so we don't need to do that separately.

We should also be using an unsigned type for the size, or at least checking that the integer isn't less than or equal to zero.

    template<class OutputT>
    constexpr auto cast()
    {
        return this->transform([](ElementT element) { return static_cast<OutputT>(element); });
    }

I don't think this belongs inside the class. If we want to do something with the internal data, it would be better to expose the internal data in a function like:

    std::vector<ElementT> const& getImageData() const { return m_imageData; }

The user can then use std::transform or whatever as they need to.

    constexpr auto get(const unsigned int locationx, const unsigned int locationy)
    {
        return this->image_data[locationy][locationx];
    }

    constexpr auto set(const unsigned int locationx, const unsigned int locationy, const ElementT& element)
    {
        this->image_data[locationy][locationx] = element;
        return *this;
    }

    template<class InputT>
    constexpr auto set(const unsigned int locationx, const unsigned int locationy, const InputT& element)
    {
        this->image_data[locationy][locationx] = static_cast<ElementT>(element);
        return *this;
    }

It's more idiomatic in C++ to provide two versions of a get function, one returning a mutable reference, and the other a const-reference.

constexpr ElementT& at(const unsigned int x, const unsigned int y) { return m_image_data[y * width + x]; }
constexpr ElementT const& at(const unsigned int x, const unsigned int y) const { return m_image_data[y * width + x]; }

    constexpr auto getData()
    {
        return this->transform([](ElementT element) { return element; });   //  Deep copy
    }

As with cast above we don't need this function. We can just expose a const& of the image data.

    void print()
    {
        for (auto& row_element : this->toString())
        {
            for (auto& element : row_element)
            {
                std::cout << element << "\t";
            }
            std::cout << "\n";
        }
        std::cout << "\n";
        return;
    }

It's more flexible to provide an output stream operator, since we can specify the stream to print to.

We definitely don't want to be creating a temporary copy filled with strings to do this!!!

We can simply use ElementT's operator<< and send the elements directly to the output stream.

    constexpr auto toString()
    {
        return this->transform([](ElementT element) { return std::to_string(element); });
    }

This seems an unlikely thing for anyone to want to do. Again, we can expose the image data by const& and let the user use std::transform if they really want to.

    constexpr auto bicubicInterpolation(const int& newSizeX, const int& newSizeY) ...

There's no point in using references for the size parameters. Again, they should be unsigned, or we have to check that the values are above 0.

This function would probably be better as a free function, e.g.:

   template<ElementT>
   Image<ElementT> copyResizeBicubic(Image<ElementT> const& image, int width, int height);

    Image<ElementT>& operator=(Image<ElementT> const& input)  //  Copy Assign
    {
        this->image_data = input.getData();
        return *this;
    }

    Image<ElementT>& operator=(Image<ElementT>&& other)       //  Move Assign
    {
        this->image_data = std::move(other.image_data);
        std::cout << "move assigned\n";
        return *this;
    }

    Image(const Image<ElementT> &input)                      //  Copy Constructor
    {
        this->image_data = input.getData();
    }

    /*    Move Constructor
     */
    Image(Image<ElementT> &&input) : image_data(std::move(input.image_data))
    {
    }

These can all be = default no?

Note that we can write Image<ElementT> as just Image inside the class.

    template<class InputT>
    constexpr 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 clip(cubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
    }

    template<class InputT1, class InputT2>
    constexpr 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));
    }

These functions don't need to access any class data. They be static, or even non-member helper functions in another file.

As small POD types, the parameters should be passed by value, not const&.

Thank you for the answer. We can simply use ElementT's operator<< and send the elements directly to the output stream. How about TinyDIP::Image<GrayScale> image1 here? The GrayScale type is defined from unsigned char (8-bits per pixel), but operator<< in unsigned char can't print number correctly. The std::to_string is still needed here. — JimmyHu, Commented Jun 24, 2021 at 7:34
True. I'd forgotten about that. It is possible to use the + trick to force numeric output instead, see: isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number — user673679, Commented Jun 24, 2021 at 7:51

Stack Exchange Network

Two dimensional bicubic interpolation implementation in C++

1 Answer 1

Your Answer

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

Linked

Hot Network Questions

Two dimensional bicubic interpolation implementation 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++reinventing-the-wheelimagetemplatec++20 or ask your own question.

Linked

Related

Hot Network Questions

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