6
\$\begingroup\$

This is a follow-up question for SIFT Keypoint Detection for Image in C++, difference_of_gaussian Template Function Implementation for Image in C++, conv2 Template Function Implementation for Image in C++ and imgaussfilt Template Function Implementation for Image in C++. I am trying to implement circle drawing feature in TinyDIP library in C++. The usage of draw_circle template function here is to plot the maximum gradient magnitude of the neighborhood of each SIFT keypoint in the test image.

The experimental implementation

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

    namespace TinyDIP
{
    //  draw_circle template function implementation
    template<typename ElementT>
    constexpr static auto draw_circle(
        const Image<ElementT>& input,
        std::tuple<std::size_t, std::size_t> central_point,
        std::size_t radius = 2,
        ElementT draw_value = ElementT{}
    )
    {
        if (input.getDimensionality() != 2)
        {
            throw std::runtime_error("Unsupported dimension!");
        }
        auto point_x = std::get<0>(central_point);
        auto point_y = std::get<1>(central_point);
        auto output = input;
        auto height = input.getHeight();
        auto width = input.getWidth();
        #pragma omp parallel for collapse(2)
        for (std::size_t y = point_y - radius; y <= point_y + radius; ++y)
        {
            for (std::size_t x = point_x - radius; x <= point_x + radius; ++x)
            {
                if (x >= width || y >= height)
                {
                    continue;
                }
                if(std::abs(std::pow(static_cast<double>(x) - static_cast<double>(point_x), 2.0) +
                std::pow(static_cast<double>(y) - static_cast<double>(point_y), 2.0) - std::pow(radius, 2)) < radius * 2)
                {
                    output.at_without_boundary_check(x, y) = draw_value;
                }
            }
        }
        return output;
    }
}

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

    namespace TinyDIP
{
    namespace SIFT_impl {
        //  get_orientation_histogram template function implementation
        template<typename ElementT>
        requires((std::floating_point<ElementT> || std::integral<ElementT>))
        constexpr static auto get_orientation_histogram(
            const Image<ElementT>& input,
            std::tuple<std::size_t, std::size_t> point,
            std::size_t block_size = 3
        )
        {
            if (input.getDimensionality() != 2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            std::vector<double> raw_histogram;
            raw_histogram.resize(37);
            for (std::size_t y = std::get<1>(point) - block_size; y <= std::get<1>(point) + block_size; ++y)
            {
                for (std::size_t x = std::get<0>(point) - block_size; x <= std::get<0>(point) + block_size; ++x)
                {
                    if (x >= input.getWidth() || y >= input.getHeight())
                    {
                        continue;
                    }
                    auto each_pixel_orientation = compute_each_pixel_orientation(subimage(input, 3, 3, x, y));
                    std::size_t bin_index = static_cast<std::size_t>(std::get<1>(each_pixel_orientation) / 10.0);
                    raw_histogram[bin_index] += std::get<0>(each_pixel_orientation);
                }
            }
            return raw_histogram;
        }
    }
}

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

    namespace TinyDIP
{
    namespace SIFT_impl {
        //  compute_each_pixel_orientation template function implementation
        /*  input is 3 * 3 image, calculate the gradient magnitude
        *    M(1, 1) = ((input(2, 1) - input(0, 1))^(2) + (input(1, 2) - input(1, 0))^(2))^(1/2)
        *   orientation
        *    θ(1, 1) = tan^(-1)((input(1, 2) - input(1, 0)) / (input(2, 1) - input(0, 1)))
        *   the value range of orientation is 0° ~ 360°
        */
        template<typename ElementT>
        constexpr static auto compute_each_pixel_orientation(const Image<ElementT>& input)
        {
            if (input.getDimensionality() != 2)
            {
                throw std::runtime_error("Unsupported dimension!");
            }
            if (input.getWidth() != 3 || input.getHeight() != 3)
                throw std::runtime_error("Input size error!");
            double gradient_magnitude =
                std::sqrt(
                    std::pow((static_cast<double>(input.at_without_boundary_check(2, 1)) - static_cast<double>(input.at_without_boundary_check(0, 1))), 2.0) +
                    std::pow((static_cast<double>(input.at_without_boundary_check(1, 2)) - static_cast<double>(input.at_without_boundary_check(1, 0))), 2.0)
                );
            double orientation = std::atan2(
                    static_cast<double>(input.at_without_boundary_check(1, 2)) - static_cast<double>(input.at_without_boundary_check(1, 0)),
                    static_cast<double>(input.at_without_boundary_check(2, 1)) - static_cast<double>(input.at_without_boundary_check(0, 1))
                );
            orientation *= (180.0 / std::numbers::pi_v<double>);
            orientation += 180;
            return std::make_tuple(gradient_magnitude, orientation);
        }
    }
}

The usage of draw_circle template function:

int main()
{
    auto start = std::chrono::system_clock::now();
    std::string file_path = "InputImages/1";
    auto bmp1 = TinyDIP::bmp_read(file_path.c_str(), false);
    bmp1 = copyResizeBicubic(bmp1, bmp1.getWidth() * 2, bmp1.getHeight() * 2);
    auto v_plane = TinyDIP::getVplane(TinyDIP::rgb2hsv(bmp1));
    auto SIFT_keypoints = TinyDIP::SIFT_impl::get_potential_keypoint(v_plane);
    std::cout << "SIFT_keypoints = " << SIFT_keypoints.size() << "\n";
    bmp1 = TinyDIP::draw_points(bmp1, SIFT_keypoints);
    for (auto&& each_SIFT_keypoint : SIFT_keypoints)
    {
        auto orientation_histogram = TinyDIP::SIFT_impl::get_orientation_histogram(v_plane, each_SIFT_keypoint);
        RGB rgb{ 255, 255, 255 };
        bmp1 = TinyDIP::draw_circle(bmp1, each_SIFT_keypoint, TinyDIP::recursive_max(orientation_histogram), rgb);
    }
    TinyDIP::bmp_write("test20240816", bmp1);
    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() << " seconds\n";
    
    return EXIT_SUCCESS;
}

The output image:

OutputImage

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

\$\endgroup\$

1 Answer 1

Reset to default
7
\$\begingroup\$

  1. You’re ensuring, inside the loop, that you’re not outside the image on the right and bottom, but not on the left or top (x or y can be negative in that loop. Also, it would be a lot more efficient if you changed the loop limits instead of adding that test to the loop: for (std::size_t y = std::max(point_y - radius, 0); y <= std::min(point_y + radius, height - 1); ++y) etc.

  2. You’re making the loop unconditionally parallel. But the default case of radius == 2 leads to a tiny loop that will only slow down if parallelized (because creating threads is expensive).

  3. Your condition “| distance^2 - radius^2 | < 2 radius” is unclear to me. The thickness of the drawn circle is related to the square root of the radius? Where does this come from? Why does that not match what I see in your example output?

  4. There’s a much simpler and efficient algorithm for drawing a circle. Your algorithm iterates over all pixels in a square around the circle, and does a test for each. Instead, you can just loop over the pixels that form the circle: https://en.m.wikipedia.org/wiki/Midpoint_circle_algorithm

\$\endgroup\$

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.