Compiler Explorer

Source code

#include <algorithm>
#include <array>
#include <cassert>
#include <chrono>
#include <cmath>
#include <concepts>
#include <exception>
#include <execution>
#include <fstream>
#include <functional>
#include <iostream>
#include <iterator>
#include <numbers>
#include <numeric>
#include <ranges>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>

using BYTE = unsigned char;

struct RGB
{
    BYTE channels[3];
};

using GrayScale = BYTE;

namespace TinyDIP
{
    #define is_size_same(x, y) {assert(x.getWidth() == y.getWidth()); assert(x.getHeight() == y.getHeight());}

//  Reference: https://stackoverflow.com/a/58067611/6667035
    template <typename T>
    concept arithmetic = std::is_arithmetic_v<T>;

//  recursive_depth function implementation
    template<typename T>
    constexpr std::size_t recursive_depth()
    {
        return 0;
    }

template<std::ranges::input_range Range>
    constexpr std::size_t recursive_depth()
    {
        return recursive_depth<std::ranges::range_value_t<Range>>() + 1;
    }
    
    //  recursive_invoke_result_t implementation
    template<std::size_t, typename, typename>
    struct recursive_invoke_result { };

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

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

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

//  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)
        {
            static_assert(unwrap_level <= recursive_depth<Arg1>(),
                "unwrap level higher than recursion depth of input");
            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...);
        }
    }

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(std::vector<ElementT> input, std::size_t newWidth, std::size_t newHeight):
            width(newWidth),
            height(newHeight)
        {
            if (input.size() != newWidth * newHeight)
            {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
            image_data = std::move(input);
        }

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

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

constexpr std::size_t getWidth() const
        {
            return width;
        }

constexpr std::size_t getHeight() const noexcept
        {
            return height;
        }

constexpr auto getSize() noexcept
        {
            return std::make_tuple(width, height);
        }

std::vector<ElementT> const& getImageData() const { return image_data; }      //  expose the internal data

void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            for (std::size_t y = 0; y < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++x)
                {
                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                    os << +at(x, y) << separator;
                }
                os << "\n";
            }
            os << "\n";
            return;
        }

//  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 < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++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;
        }

friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
        {
            const std::string separator = "\t";
            for (std::size_t y = 0; y < rhs.height; ++y)
            {
                for (std::size_t x = 0; x < rhs.width; ++x)
                {
                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                    os << +rhs.at(x, y) << separator;
                }
                os << "\n";
            }
            os << "\n";
            return os;
        }

Image<ElementT>& operator+=(const Image<ElementT>& rhs)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            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)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            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)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            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)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            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;
        }

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:
        std::size_t width;
        std::size_t height;
        std::vector<ElementT> image_data;

void checkBoundary(const size_t x, const size_t y) const
        {
            assert(x < width);
            assert(y < height);
        }

};

template<std::floating_point ElementT = double, std::floating_point OutputT = ElementT>
    Image<OutputT> dct3_detail(const std::vector<Image<ElementT>>& input, const std::size_t plane_index)
    {
        auto N1 = static_cast<OutputT>(input[0].getWidth());
        auto N2 = static_cast<OutputT>(input[0].getHeight());
        auto N3 = input.size();
        auto alpha1 = (plane_index == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
        auto output = Image<OutputT>(input[plane_index].getWidth(), input[plane_index].getHeight());
        for (std::size_t y = 0; y < output.getHeight(); ++y)
        {
            OutputT alpha2 = (y == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
            for (std::size_t x = 0; x < output.getWidth(); ++x)
            {
                OutputT sum{};
                OutputT alpha3 = (x == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
                for (std::size_t inner_z = 0; inner_z < N3; ++inner_z)
                {
                    auto plane = input[inner_z];
                    for (std::size_t inner_y = 0; inner_y < plane.getHeight(); ++inner_y)
                    {
                        for (std::size_t inner_x = 0; inner_x < plane.getWidth(); ++inner_x)
                        {
                            auto l1 = (std::numbers::pi_v<OutputT> / (2 * N1) * (2 * static_cast<OutputT>(inner_x) + 1) * x);
                            auto l2 = (std::numbers::pi_v<OutputT> / (2 * N2) * (2 * static_cast<OutputT>(inner_y) + 1) * y);
                            auto l3 = (std::numbers::pi_v<OutputT> / (2 * static_cast<OutputT>(N3)) * (2 * static_cast<OutputT>(inner_z) + 1) * static_cast<OutputT>(plane_index));
                            sum += static_cast<OutputT>(plane.at(inner_x, inner_y)) *
                                std::cos(l1) * std::cos(l2) * std::cos(l3);
                        }
                    }
                }
                output.at(x, y) = 8 * alpha1 * alpha2 * alpha3 * sum / (N1 * N2 * N3);
            }
        }
        return output;
    }

template<std::floating_point ElementT = double, std::floating_point OutputT = ElementT>
    std::vector<Image<OutputT>> dct3(const std::vector<Image<ElementT>>& input)
    {
        std::vector<Image<OutputT>> output;
        output.resize(input.size());
        for (std::size_t i = 0; i < input.size(); ++i)
        {
            output[i] = dct3_detail<ElementT, OutputT>(input, i);
        }
        return output;
    }

template<std::floating_point ElementT = double, std::floating_point OutputT = ElementT>
    Image<OutputT> idct3_detail(const std::vector<Image<ElementT>>& input, const std::size_t plane_index)
    {
        auto N1 = static_cast<OutputT>(input[0].getWidth());
        auto N2 = static_cast<OutputT>(input[0].getHeight());
        auto N3 = input.size();
        auto output = Image<OutputT>(input[plane_index].getWidth(), input[plane_index].getHeight());
        for (std::size_t y = 0; y < output.getHeight(); ++y)
        {
            for (std::size_t x = 0; x < output.getWidth(); ++x)
            {
                OutputT sum{};
                for (std::size_t inner_z = 0; inner_z < N3; ++inner_z)
                {
                    auto plane = input[inner_z];
                    for (std::size_t inner_y = 0; inner_y < plane.getHeight(); ++inner_y)
                    {
                        for (std::size_t inner_x = 0; inner_x < plane.getWidth(); ++inner_x)
                        {
                            auto l1 = (std::numbers::pi_v<OutputT> / (2 * N1) * (2 * x + 1) * static_cast<OutputT>(inner_x));
                            auto l2 = (std::numbers::pi_v<OutputT> / (2 * N2) * (2 * y + 1) * static_cast<OutputT>(inner_y));
                            auto l3 = (std::numbers::pi_v<OutputT> / (2 * static_cast<OutputT>(N3)) * (2 * static_cast<OutputT>(plane_index) + 1) * static_cast<OutputT>(inner_z));
                            OutputT alpha1 = (inner_x == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
                            OutputT alpha2 = (inner_y == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
                            OutputT alpha3 = (inner_z == 0) ? (std::numbers::sqrt2_v<OutputT> / 2) : (OutputT{1.0});
                            sum += alpha1 * alpha2 * alpha3 * static_cast<OutputT>(plane.at(inner_x, inner_y)) *
                                std::cos(l1) * std::cos(l2) * std::cos(l3);
                        }
                    }
                }
                output.at(x, y) = sum;
            }
        }
        return output;
    }

template<std::floating_point ElementT = double, std::floating_point OutputT = ElementT>
    std::vector<Image<OutputT>> idct3(const std::vector<Image<ElementT>>& input)
    {
        std::vector<Image<OutputT>> output;
        output.resize(input.size());
        for (std::size_t i = 0; i < input.size(); ++i)
        {
            output[i] = idct3_detail<ElementT, OutputT>(input, i);
        }
        return output;
    }
}

template<typename ElementT>
void print3(std::vector<TinyDIP::Image<ElementT>> input)
{
    for (std::size_t i = 0; i < input.size(); i++)
    {
        input[i].print();
        std::cout << "*******************\n";
    }
}

void idct3Test()
{
    std::size_t N1 = 10, N2 = 10, N3 = 10;
    std::vector<TinyDIP::Image<double>> test_input;
    for (std::size_t z = 0; z < N3; z++)
    {
        test_input.push_back(TinyDIP::Image<double>(N1, N2));
    }
    for (std::size_t z = 1; z <= N3; z++)
    {
        for (std::size_t y = 1; y <= N2; y++)
        {
            for (std::size_t x = 1; x <= N1; x++)
            {
                test_input[z - 1].at(y - 1, x - 1) = x * 100 + y * 10 + z;
            }
        }
    }
    print3(test_input);

auto dct3_output = TinyDIP::dct3(test_input);
    print3(dct3_output);

auto idct3_output = TinyDIP::idct3(dct3_output);
    print3(TinyDIP::recursive_transform(
        [](auto&& input1, auto&& input2)
        {
            return input1 - input2;
        }, idct3_output, test_input));
}

int main()
{
    auto start = std::chrono::system_clock::now();
    idct3Test();
    auto end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end - start;
    std::time_t end_time = std::chrono::system_clock::to_time_t(end);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}