Midpoint Circle Algorithm Implementation for Image in C++

Question

This is a follow-up question for draw_circle Template Function Implementation for Image in C++, 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 follow the concept of Midpoint Circle Algorithm to implement draw_circle template function in this post.

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();
          std::size_t l = static_cast<std::size_t>(radius) * std::cos(std::numbers::pi_v<double> / 4.0);
          for (std::size_t x = 0; x <= l; ++x)
          {
              std::size_t y = static_cast<std::size_t>(std::sqrt(radius * radius - x * x));
              std::size_t location_x1 = point_x + x;
              std::size_t location_y1 = point_y + y;
              std::size_t location_x2 = point_x + x;
              std::size_t location_y2 = point_y - y;
              std::size_t location_x3 = point_x - x;
              std::size_t location_y3 = point_y + y;
              std::size_t location_x4 = point_x - x;
              std::size_t location_y4 = point_y - y;
              std::size_t location_x5 = point_x + y;
              std::size_t location_y5 = point_y + x;
              std::size_t location_x6 = point_x + y;
              std::size_t location_y6 = point_y - x;
              std::size_t location_x7 = point_x - y;
              std::size_t location_y7 = point_y + x;
              std::size_t location_x8 = point_x - y;
              std::size_t location_y8 = point_y - x;
  
  
              if (location_x1 >= output.getWidth() ||
                  location_y1 >= output.getHeight() ||
                  location_x2 >= output.getWidth() ||
                  location_y2 >= output.getHeight() ||
                  location_x3 >= output.getWidth() ||
                  location_y3 >= output.getHeight() ||
                  location_x4 >= output.getWidth() ||
                  location_y4 >= output.getHeight() ||
                  location_x5 >= output.getWidth() ||
                  location_y5 >= output.getHeight() ||
                  location_x6 >= output.getWidth() ||
                  location_y6 >= output.getHeight() ||
                  location_x7 >= output.getWidth() ||
                  location_y7 >= output.getHeight() ||
                  location_x8 >= output.getWidth() ||
                  location_y8 >= output.getHeight())
              {
                  continue;
              }
              output.at_without_boundary_check(location_x1, location_y1) = draw_value;
              output.at_without_boundary_check(location_x2, location_y2) = draw_value;
              output.at_without_boundary_check(location_x3, location_y3) = draw_value;
              output.at_without_boundary_check(location_x4, location_y4) = draw_value;
              output.at_without_boundary_check(location_x5, location_y5) = draw_value;
              output.at_without_boundary_check(location_x6, location_y6) = draw_value;
              output.at_without_boundary_check(location_x7, location_y7) = draw_value;
              output.at_without_boundary_check(location_x8, location_y8) = draw_value;
          }
          return output;
      }
  }
  

The example output:

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

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  draw_circle Template Function Implementation for Image in C++,

  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++

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

  I am trying to follow the concept of Midpoint Circle Algorithm to implement draw_circle template function in this post.

• Why a new review is being asked for?

  Please review the implementation of draw_circle template function and its tests.

Answer 1

You have a serious bug. Each loop iteration draws 8 pixels. If one of those pixels fails a bounds check, you skip the entire loop iteration. Correct behavior is to draw the other 7 pixels as normal. The bug is clearly visible in the image you posted. Entire segments are missing in some of the circles.

Use std::ptrdiff_t for quantities that can be negative, as std::size_t is unsigned. If you subtract two unsigned quantities, you might get a negative result which then wraps into a large positive number. You deliberately allowed this to happen to avoid a < 0 check, but that's bad style in my opinion and makes the code harder to read.

Don't repeat yourself. Write a helper function that performs the bounds check and then draws the pixel, or does nothing if it would be out of bounds. I called it draw_if_possible here. I didn't implement it, I'll leave that as an exercise to the reader.

Finally, don't mix floating point and integer calculations. The midpoint circle algorithm uses integers. If you implement the conditions correctly you don't need floating point.

Here is the improved code:

if (r <= 0)
{
    // early out avoids y going negative in loop
    return output;
}
for (std::ptrdiff_t x = 0, y = r; x <= y; x++)
{
    // try to decrement y, then accept or revert
    y--;
    if (x * x + y * y < r * r)
    {
        y++;
    }
    // do nothing if out of bounds, otherwise draw
    draw_if_possible(output, value, mx + x, my + y);
    draw_if_possible(output, value, mx - x, my + y);
    draw_if_possible(output, value, mx + x, my - y);
    draw_if_possible(output, value, mx - x, my - y);
    draw_if_possible(output, value, mx + y, my + x);
    draw_if_possible(output, value, mx - y, my + x);
    draw_if_possible(output, value, mx + y, my - x);
    draw_if_possible(output, value, mx - y, my - x);
}
return output;

Answer 2

The very first thing is that chunks of (almost) repeated lines is a big code smell. There is definitely a way to rewrite it to avoid code duplication, which also prevents bugs. In really old computational C/FORTRAN code sometimes people will avoid using a struct and pass around ten arguments, but I would really suggest using a struct or class for Point.x and Point.y. Then you can put the points in an array and loop over the points.

Also in your main function I would suggest judicious use of newlines and comments to make the code easier to read. Right now it appears as a "wall of text" which isn't good for code or prose. You can use the time utility from the command line to avoid manually writing a timing computation (neither are very accurate for precise benchmarking).