8
\$\begingroup\$

This is a follow-up question for Gaussian Fisheye Image Generator Implementation in C++ and An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. I am trying to implement image rotation function in this post. The example output:

Image Input Output
Input rotatedImage

The experimental implementation

  • rotate, rotate_detail and rotate_degree template functions implementation (in file image_operations.h)

    rotate_detail template function performs image rotation operation between 0° to 90°.

    namespace TinyDIP
    {
        //  rotate_detail template function implementation
        //  rotate_detail template function performs image rotation between 0° to 90°
        template<arithmetic ElementT, std::floating_point FloatingType = double>
        constexpr static auto rotate_detail(const Image<ElementT>& input, FloatingType radians)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            //  if 0° rotation case
            if (radians == 0)
            {
                return input;
            }
            //  if 90° rotation case
            if(radians == std::numbers::pi_v<long double> / 2.0)
            {
                Image<ElementT> output(input.getHeight(), input.getWidth());
                for (std::size_t y = 0; y < input.getHeight(); ++y)
                {
                    for (std::size_t x = 0; x < input.getWidth(); ++x)
                    {
                        output.at(input.getHeight() - y - 1, x) = 
                            input.at(x, y);
                    }
                }
                return output;
            }
    
            FloatingType half_width = static_cast<FloatingType>(input.getWidth()) / 2.0;
            FloatingType half_height = static_cast<FloatingType>(input.getHeight()) / 2.0;
            FloatingType new_width = 2 * 
                std::hypot(half_width, half_height) *
                std::abs(std::sin(std::atan2(half_width, half_height) + radians));
            FloatingType new_height = 2 *
                std::hypot(half_width, half_height) *
                std::abs(std::sin(std::atan2(half_height, half_width) + radians));
    
            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 - half_width;
                    FloatingType distance_y = y - half_height;
                    FloatingType distance = std::hypot(distance_x, distance_y);
                    FloatingType angle = std::atan2(distance_y, distance_x) + radians;
    
                    FloatingType width_ratio = new_width / static_cast<FloatingType>(input.getWidth());
                    FloatingType height_ratio = new_height / static_cast<FloatingType>(input.getHeight());
                    FloatingType distance_weight = (input.getWidth() > input.getHeight())?(1 / height_ratio):(1 / width_ratio);
                    FloatingType new_distance = 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 + half_width),
                        static_cast<std::size_t>(new_distance_y + half_height)) = 
                        input.at(x, y);
                }
            }
            return output;
        }
    
        //  rotate template function implementation
        template<arithmetic ElementT, std::floating_point FloatingType = double>
        constexpr static auto rotate(const Image<ElementT>& input, FloatingType radians)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            auto output = input;
            while(radians >= 2 * std::numbers::pi_v<long double>)
            {
                radians = radians - 2 * std::numbers::pi_v<long double>;
            }
            while(radians > std::numbers::pi_v<long double> / 2.0)
            {
                output = rotate_detail(output, std::numbers::pi_v<long double> / 2.0);
                radians-=(std::numbers::pi_v<long double> / 2.0);
            }
            output = rotate_detail(output, radians);
            return output;
        }
    
        //  rotate template function implementation
        template<typename ElementT, class FloatingType = double>
        requires ((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, HSV>))
        constexpr static auto rotate(const Image<ElementT>& input, FloatingType radians)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            return apply_each(input, [&](auto&& planes) { return rotate(planes, radians); });
        }
    
        //  rotate template function implementation
        template<arithmetic ElementT, std::integral T = int>
        constexpr static auto rotate(const Image<ElementT>& input, T radians)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            return rotate(input, static_cast<double>(radians));
        }
    
        //  rotate_degree template function implementation
        template<typename ElementT, class T = double>
        constexpr static auto rotate_degree(const Image<ElementT>& input, T degrees)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            return rotate(input, static_cast<double>(degrees) * std::numbers::pi_v<long double> / 180.0);
        }
    }
    
  • transpose template function implementation (in file image_operations.h)

    namespace TinyDIP
    {
        //  transpose template function implementation
        template<typename ElementT>
        constexpr static auto transpose(const Image<ElementT>& input)
        {
            if (input.getDimensionality()!=2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            Image<ElementT> output(input.getHeight(), input.getWidth());
            for (std::size_t y = 0; y < input.getHeight(); ++y)
            {
                for (std::size_t x = 0; x < input.getWidth(); ++x)
                {
                    output.at(y, x) = 
                        input.at(x, y);
                }
            }
            return output;
        }
    }
    
  • 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 rotate function:

std::string file_path = "InputImages/1";
auto bmp1 = TinyDIP::bmp_read(file_path.c_str(), false);
bmp1 = rotate(bmp1, 1.0);
TinyDIP::bmp_write("test", bmp1);

The usage transpose function:

template<typename ElementT>
void transposeTest(const size_t size = 5)
{
    std::cout << "Test with 2D image:\n";
    auto image2d = TinyDIP::Image<ElementT>(size, size);
    image2d.setAllValue(1);
    image2d.at(0, 1) = 3;
    image2d.print();
    TinyDIP::transpose(image2d).print();
    
    return;
}

Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

\$\endgroup\$

2 Answers 2

5
\$\begingroup\$

Avoid trigonometric functions where possible

Trigonometric functions like atan2(), sin() and cos() are very slow. You do three of those operations per pixel, which is going to be the major performance issue with your code.

There is a way to avoid doing any per-pixel trigonometry though, and that is by decomposing the rotation into three shear operations, like for example shown in this StackOverflow question.

Loop over the output pixels

You loop over the input pixels, and remap their coordinates to the output. However, you will then have a risk that you miss some output pixels, as some might be skipped due to floating point rounding. Maybe you get lucky most of the time, but to avoid the problem, loop over the output pixels, and remap their coordinates to the input to find out which color to copy.

Even better, if you use the shearing technique then you can avoid this issue completely.

Handle negative rotation angles

rotate() seems to only check for radians being too large in the positive direction, but doesn't check for potentially large negative values.

I also recommend you use std::fmod() or something similar to truncate the angle here, instead of using a while-loop.

Avoid rotating up to four times

If the angle is between \$1\frac{1}{2} \pi\$ and \$2 \pi\$, you are performing four rotations. You should be able to rotate everything in one go (or if you use shear operations, using not more than 3 shears).

If your algorithm really can't handle rotations larger than 90 degrees, then at least implement a fast algorithm to do 90, 180 and 270 degree rotations, as these are trivial to implement.

\$\endgroup\$
3
  • 2
    \$\begingroup\$ "decomposing the rotation into three shear operations, like for example shown in this StackOverflow question." - Matt Parker did a video on this very recently: youtube.com/watch?v=1LCEiVDHJmc. No affiliation, just a demonstration I thought was neat. \$\endgroup\$ Mar 26 at 19:11
  • 1
    \$\begingroup\$ Another advantage of the three sheer technique is that you interpolate in 1D only, which is much simpler than interpolating in 2D when sampling a pixel at an arbitrary location in the image. \$\endgroup\$ Mar 27 at 1:23
  • \$\begingroup\$ @G. Sliepen Thank you for letting me know the shearing technique. It looks good to deal with image rotation with it. :) \$\endgroup\$
    – JimmyHu
    Mar 27 at 5:48
2
\$\begingroup\$

G. Sliepen already pointed out the three sheer technique, the need to iterate over the output pixels rather than the input, and the expense of doing large rotations 90 degrees at the time. There are a few things to add here:

  1. The three sheer technique is not only cheaper, it also makes adding interpolation very easy, as all you do is shift lines of pixels over by a certain distance.

  2. You do a lot of computation repeatedly inside the loop, you should move these outside the loop. For example, rotated coordinates are computed by multiplying the input coordinates by a rotation matrix, which is fixed for the whole operation. So inside the loop all you have is four multiplications and two additions to rotate each coordinate. Though of course you should iterate over the output image, and so you use the inverse rotation matrix (or negate the angle) to find the input coordinate corresponding to each output pixel.

  3. Indexing using .at(x,y) is quite expensive. To access an arbitrary element it is great, but in a loop that iterates over each pixel, it is much cheaper to add one to the pointer for one pixel to get to the next pixel. You also don’t want to continuously test for indices being inside the bounds, you already know they are.

  4. requires ((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, HSV>)) means that you need to rewrite this function (and all the other ones in the library that do a similar thing) if you ever decide to add a color space. I would suggest making this more generic for all color images. Maybe have a base class for both RGB and HSV, that you can use in the requirement?


PS: looking up the definition of apply_each in your repo (it’s not included in the question), I see:

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)
    {…}

That requirement makes T not a template parameter. Prefer:

template<class F, class... Args>
    constexpr static auto apply_each(Image<RGB> input, F operation, Args&&... args)
    {…}

Just because it’s so much simpler. Also:

template<class T = RGB>
    requires (std::same_as<T, RGB>)
    constexpr static auto getRplane(Image<T> input)
    {…}

This is not a template at all, just do:

    constexpr static auto getRplane(Image<RGB> input)
    {…}
\$\endgroup\$
8
  • \$\begingroup\$ > (Image<RGB> input...) makes error: 'input' has incomplete type error. \$\endgroup\$
    – JimmyHu
    Mar 27 at 4:46
  • \$\begingroup\$ I think it's called "partial template specialization" and it is not allowed. \$\endgroup\$
    – JimmyHu
    Mar 27 at 4:53
  • \$\begingroup\$ @JimmyHu Do you mean in apply_each? You'd have to update each of the templates in the same way \$\endgroup\$ Mar 27 at 6:34
  • \$\begingroup\$ Yes, I already updated each of apply_each templates in the same way but the error is still occurred. \$\endgroup\$
    – JimmyHu
    Mar 27 at 6:59
  • 1
    \$\begingroup\$ @JimmyHu The templated image processing libraries I know do something like Image<vec3<float>> for color images, without specifying the color space. In your case I was thinking of struct Color {}; struct RGB : Color { uint8 value[3]; }; and so on. Then you can still do Image<RGB>, but in your templated functions you can restrict to T being derived from Color. There might be other solutions. But you do want to avoid listing all the color spaces your library knows for every function that applies generically to color images. \$\endgroup\$ Mar 28 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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