Why the corners detected by Fast algorithm in OPENCV and my own implementation are not same?

Question

I have implemented a naive FAST corner detector using C++ boost algorithm. Currently it does not have non max suppression. I ran my implementation and the OPENCV implementation on this photo below

The output of OPENCV is

The output from my implementation is

Here is my code

#include <bits/stdc++.h>
#include <boost/gil.hpp>
#include <boost/gil/extension/io/jpeg.hpp>
#include <boost/gil/extension/io/png.hpp>
#include <boost/gil/image_processing/circle.hpp>
namespace gil = boost::gil;
template <typename T, typename U>
bool FastCornerDetector(const T& buffer, int r, int c, std::vector<U> points,
                        int t = 20) {
    int valid_points_count = 16;
    for (auto& u :
         points) // checking whether all the 16 points lie within the image
    {
        if (u[1] + r >= buffer.height() || u[0] + c >= buffer.width()) {
            return false;
        } else if (u[1] + r < 0 || u[0] + c < 0) {
            return false;
        } else {
            u += gil::point_t(c, r);
        }
    }

    int threshold_indicator[16], index = -1;
    memset(threshold_indicator, -1, sizeof(threshold_indicator));
    // marking the pixel value of every pixel as >I_p+t or <I_p-t or between
    // these two
    for (const auto& point : points) {
        if (buffer(point) < buffer(gil::point_t(c, r)) - t) {
            threshold_indicator[++index] = -1;
        } else if (buffer(point) > buffer(gil::point_t(c, r)) + t) {
            threshold_indicator[++index] = 1;
        } else {
            threshold_indicator[++index] = 0;
        }
    }

    bool is_feature_point = false;
    // threshold check for n=9 consecutive points. I am not doing the speed test
    std::ptrdiff_t count = 0;
    for (int i = 0; i < 25; i++) {
        if (threshold_indicator[i] == -1) {
            if (++count > 8) {
                is_feature_point = true;
                break;
            }
        } else
            count = 0;
    }

    if (!is_feature_point) {
        count = 0;
        for (int i = 0; i < 25; i++) {
            if (threshold_indicator[i] == 1) {
                if (++count > 8) {
                    is_feature_point = true;
                    break;
                }
            } else
                count = 0;
        }
    }
    return is_feature_point;
}

int main() {
    const std::ptrdiff_t radius = 3;
    const auto rasterizer       = gil::midpoint_circle_rasterizer{};
    int number_of_points        = rasterizer.point_count(radius);
    std::vector<gil::point_t> circle_points(number_of_points);
    rasterizer(radius, {0, 0}, circle_points.begin());

    auto cmp = [](const auto& a, const auto& b) {
        return a[0] < b[0] || a[1] < b[1];
    };
    std::set<gil::point_t, decltype(cmp)> s(cmp);
    for (const auto& point : circle_points) {
        s.insert(point);
    }

    // std::cout<<s.size()<<std::endl;
    std::vector<gil::point_t> points_clockwise[4];
    for (const auto& point : s) {
        if (point[1] > 0 && point[0] >= 0) {
            points_clockwise[0].push_back(point);
        } else if (point[1] <= 0 && point[0] > 0) {
            points_clockwise[1].push_back(point);
        } else if (point[1] < 0 && point[0] <= 0) {
            points_clockwise[2].push_back(point);
        } else if (point[1] >= 0 && point[0] < 0) {
            points_clockwise[3].push_back(point);
        }
    }

    std::sort(points_clockwise[0].begin(), points_clockwise[0].end(),
              [](const auto& a, const auto& b) -> bool {
                  if (a[1] == b[1])
                      return a[0] < b[0];
                  return a[1] > b[1];
              });

    std::sort(points_clockwise[1].begin(), points_clockwise[1].end(),
              [](const auto& a, const auto& b) -> bool {
                  if (a[1] == b[1])
                      return a[0] > b[0];
                  return a[1] > b[1];
              });

    std::sort(points_clockwise[2].begin(), points_clockwise[2].end(),
              [](const auto& a, const auto& b) -> bool {
                  if (a[1] == b[1])
                      return a[0] > b[0];
                  return a[1] < b[1];
              });

    std::sort(points_clockwise[3].begin(), points_clockwise[3].end(),
              [](const auto& a, const auto& b) -> bool {
                  if (a[1] == b[1])
                      return a[0] < b[0];
                  return a[1] < b[1];
              });

    std::vector<gil::point_t> final_points_clockwise;

    // all the 16 points on the circumference of the circle are present in
    // clockwise format in final_points_clockwise
    final_points_clockwise.insert(final_points_clockwise.end(),
                                  points_clockwise[0].begin(),
                                  points_clockwise[0].end());
    final_points_clockwise.insert(final_points_clockwise.end(),
                                  points_clockwise[1].begin(),
                                  points_clockwise[1].end());
    final_points_clockwise.insert(final_points_clockwise.end(),
                                  points_clockwise[2].begin(),
                                  points_clockwise[2].end());
    final_points_clockwise.insert(final_points_clockwise.end(),
                                  points_clockwise[3].begin(),
                                  points_clockwise[3].end());

    gil::gray8_image_t input_image;
    gil::rgb8_image_t input_color_image;
    gil::rgb8_view_t input_color_image_view;
    gil::read_image("box_image.png", input_image, gil::png_tag{});
    gil::read_image("box.jpg", input_color_image, gil::jpeg_tag{});
    input_color_image_view = gil::view(input_color_image);
    auto input_image_view =
        gil::color_converted_view<gil::gray8_pixel_t>(gil::view(input_image));
    s.clear();
    for (int i = 0; i < input_image_view.height(); i++) {
        for (int j = 0; j < input_image_view.width(); j++) {
            if (FastCornerDetector(input_image_view, i, j,
                                   final_points_clockwise)) {
                // if it is a corner draw a circle against it
                for (auto u : final_points_clockwise) {
                    input_color_image_view(u + gil::point_t(j, i)) =
                        gil::rgb8_pixel_t(0, 255, 0);
                }
            }
        }
    }

    gil::write_view("Fast_Output.jpeg", input_color_image_view,
                    gil::jpeg_tag{});
}

I have not used speed test and non max suppression but I do not think that is the reason behind the corners not getting detected. The coordinates for the bresenham circle of radius 3 is correct calculated in my code as I verified. I am unable to find the reason then why is the output differing. Why less corners are detected in my output? The threshold in both the cases have been set to 20. Please , if anyone can kindly help me...

Answer 1

I can't compile your code. I can only spot a few random things, that might help you fix things:

1. Your threshold_indicator array is sixteen elements, but your loops go to index 24 This will lead to Undefined Behaviour (unless the loop incidentally breaks before that point)

2. valid_points_count is unused, yet it is being initialized to a magic number

3. You can reduce a lot of code by not repeating yourself quite as much. E.g your code for the following bit took 45 lines of code. That's >4x as much.

   std::vector<gil::point_t> final_points_clockwise;
   for (auto& pc : points_clockwise) {
       std::sort(pc.begin(), pc.end(),
                 [](const auto& a, const auto& b) -> bool {
                     if (a[1] == b[1])
                         return a[0] < b[0];
                     return a[1] > b[1];
                 });
       final_points_clockwise.insert(final_points_clockwise.end(), pc.begin(),
                                     pc.end());
   }
   // all the 16 points on the circumference of the circle are present in
   // clockwise format in final_points_clockwise
   

   In fact, further down it becomes just

    std::vector<gil::point_t> final_cw;
    for (auto& pc : points_clockwise) {
          std::sort(pc.begin(), pc.end(), ClockwiseCmp{});
          final_cw.insert(final_cw.end(), pc.begin(), pc.end());
    }
   

4. Similarly, all of

   std::ptrdiff_t count = 0;
   for (int i = 0; i < 25; i++) {
       if (threshold_indicator[i] == -1) {
           if (++count > 8) {
               is_feature_point = true;
               break;
           }
       } else
           count = 0;
   }
   if (!is_feature_point) {
       count = 0;
       for (int i = 0; i < 25; i++) {
           if (threshold_indicator[i] == 1) {
               if (++count > 8) {
                   is_feature_point = true;
                   break;
               }
           } else
               count = 0;
       }
   }
   return is_feature_point;
   

   Can be replaced by

   bool is_feature_point =
       indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, -1) ||
       indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, 1);
   

   This at once removes the need for initializing the indicator array and the problem that the loop index goes beyond the indicator size.

5. In fact I'd probably write the entire FastCornerDetector as

   template <typename T, typename U>
   bool FastCornerDetector(T const& buffer, int r, int c,
                           std::vector<U> points, int t = 20)
   {
       assert(points.size() == 16);
   
       if (std::any_of(points.begin(), points.end(), [&](auto& u) {
               u += gil::point_t(c, r);
               return (u[1] >= buffer.height() || u[0] >= buffer.width()) ||
                   (u[1] < 0 || u[0] < 0);
           })) {
           return false;
       }
   
       // marking every pixel as >I_p+t or <I_p-t
       auto const I_p = buffer(gil::point_t(c, r));
   
       std::vector<int> indicator;
       std::transform(
           points.begin(), points.end(), back_inserter(indicator),
           [&buffer, low = I_p - t, hi = I_p + t](auto const& point) {
               if (buffer(point) < low)     return -1;
               else if (buffer(point) > hi) return 1;
               else                         return 0;
           });
   
       bool is_feature_point =
           indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, -1) ||
           indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, 1);
       return is_feature_point;
   }
   

   That's

   • less than half the number of lines of code
   • more than double the readability¹
   • much less error prone because
     • no manual (loop) indexing
     • no magic numbers
   • easier to maintain; e.g. it's now easy to rephrase the feature point criterium as something else

6. This comparison looks broken:

   auto cmp = [](const auto& a, const auto& b) {
       return a[0] < b[0] || a[1] < b[1];
   };
   

   I'd expect

   auto cmp = [](const auto& a, const auto& b) {
       if (a[0] < b[0])  return true;
       if (a[0] == b[0]) return a[1] < b[1];
       return false;
   };
   

   Or, you know, less error prone:

   auto cmp = [](const auto& a, const auto& b) {
       return std::tie(a[0], a[1]) < std::tie(b[0], b[1]);
   };
   

   In fact I'd probably define it out-of-line:

7. And since you don't need to know the number of rasterized points beforehand, and don't ever use the circle_points vector except to fill the set, why not just

   constexpr std::ptrdiff_t radius = 3;
   auto const rasterizer = gil::midpoint_circle_rasterizer{};
   
   std::set<gil::point_t, PointCmp> s;
   rasterizer(radius, {0, 0}, inserter(s, s.end()));
   

8. Trying to simplify the point separation:

   std::vector<gil::point_t> points_clockwise[4];
   for (const auto& point : s) {
       auto index = 0;
       if (false);
       else if (point[1] >  0 && point[0] >= 0) { index = 0; }
       else if (point[1] <= 0 && point[0] >  0) { index = 1; }
       else if (point[1] <  0 && point[0] <= 0) { index = 2; }
       else if (point[1] >= 0 && point[0] <  0) { index = 3; }
       points_clockwise[index].push_back(point);
   }
   

   Makes me question whether the inconsistent(?) boundary handling is on purpose. I'm assuming it might be. My sense of geometry is not very good.

9. The (counter)clockwise sort looks interesting. How does it work? I've not found any code on the internet that doesn't look at the angles to achieve this. Regardless, you can at least refactor the code into an equivalent comparator again using the trick from above (proof of quivalence):

   struct ClockwiseCmp {
       bool operator()(gil::point_t const& a, gil::point_t const& b) const {
           return std::make_tuple(-a[1], a[0]) < std::make_tuple(-b[1], b[0]);
       }
   };
   
   std::vector<gil::point_t> final_cw;
   for (auto& pc : points_clockwise) {
       std::sort(pc.begin(), pc.end(), ClockwiseCmp{});
       final_cw.insert(final_cw.end(), pc.begin(), pc.end());
   }
   

LISTING

Fully reviewed/refactored listing. Note that I wasn't able to run any of the code because of the missing GIL headers.

"Live" On Compiler Explorer

#include <boost/gil.hpp>
#include <boost/gil/extension/io/jpeg.hpp>
#include <boost/gil/extension/io/png.hpp>
#include <set>

#if 0
    #include <boost/gil/image_processing/circle.hpp>
#else
    namespace boost::gil { // mockups to satisfy the compiler
        struct midpoint_circle_rasterizer{
            void operator()(ptrdiff_t, gil::point_t, auto iterator) const {}
        };
    }
#endif

namespace gil = boost::gil;

namespace /*file static*/ {
    template <typename T, typename U>
    bool FastCornerDetector(T const& buffer, int r, int c,
                            std::vector<U> points, int t = 20)
    {
        assert(points.size() == 16);

        if (std::any_of(points.begin(), points.end(), [&](auto& u) {
                u += gil::point_t(c, r);
                return (u[1] >= buffer.height() || u[0] >= buffer.width()) ||
                    (u[1] < 0 || u[0] < 0);
            })) {
            return false;
        }

        // marking every pixel as >I_p+t or <I_p-t
        auto const I_p = buffer(gil::point_t(c, r));

        std::vector<int> indicator;
        std::transform(
            points.begin(), points.end(), back_inserter(indicator),
            [&buffer, low = I_p - t, hi = I_p + t](auto const& point) {
                if (buffer(point) < low)     return -1;
                else if (buffer(point) > hi) return 1;
                else                         return 0;
            });

        bool is_feature_point =
            indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, -1) ||
            indicator.end() != std::search_n(indicator.begin(), indicator.end(), 8, 1);
        return is_feature_point;
    }

    struct PointCmp {
        bool operator()(const gil::point_t& a, const gil::point_t& b) const {
            return std::tie(a[0], a[1]) < std::tie(b[0], b[1]);
        }
    };

    struct ClockwiseCmp {
        bool operator()(gil::point_t const& a, gil::point_t const& b) const {
            return std::make_tuple(-a[1], a[0]) < std::make_tuple(-b[1], b[0]);
        }
    };
} // namespace

int main() {
    constexpr std::ptrdiff_t radius = 3;
    auto const rasterizer = gil::midpoint_circle_rasterizer{};

    std::set<gil::point_t, PointCmp> s;
    rasterizer(radius, {0, 0}, inserter(s, s.end()));

    // std::cout<<s.size()<<std::endl;
    std::vector<gil::point_t> points_clockwise[4];
#if 0 // old boundary spellings:
    for (const auto& point : s) {
        auto index = 0;
        if (false);
        else if (point[1] >  0 && point[0] >= 0) { index = 0; }
        else if (point[1] <= 0 && point[0] >  0) { index = 1; }
        else if (point[1] <  0 && point[0] <= 0) { index = 2; }
        else if (point[1] >= 0 && point[0] <  0) { index = 3; }
        points_clockwise[index].push_back(point);
    }
#else // unified version (UNTESTED):
    for (const auto& point : s) {
        bool const a = point[0] < 0;
        bool const b = point[1] > 0;

        auto index = 0;
            if (!a &&  b) { index = 0; }
        else if (!a && !b) { index = 1; }
        else if ( a && !b) { index = 2; }
        else if ( a &&  b) { index = 3; }
        points_clockwise[index].push_back(point);
    }
#endif

    std::vector<gil::point_t> final_cw;
    for (auto& pc : points_clockwise) {
        std::sort(pc.begin(), pc.end(), ClockwiseCmp{});
        final_cw.insert(final_cw.end(), pc.begin(), pc.end());
    }
    // all the 16 points on the circumference of the circle are present in
    // clockwise format in final_cw

    gil::gray8_image_t input_image;
    gil::rgb8_image_t  input_color_image;
    gil::rgb8_view_t   input_color_image_view;

    gil::read_image("box_image.png", input_image, gil::png_tag{});
    gil::read_image("box.jpg", input_color_image, gil::jpeg_tag{});

    input_color_image_view = gil::view(input_color_image);
    auto input_image_view =
        gil::color_converted_view<gil::gray8_pixel_t>(gil::view(input_image));

    for (int i = 0; i < input_image_view.height(); i++) {
        for (int j = 0; j < input_image_view.width(); j++) {
            if (FastCornerDetector(input_image_view, i, j, final_cw)) {
                // if it is a corner draw a circle against it
                for (auto& u : final_cw) {
                    input_color_image_view(u + gil::point_t(j, i)) =
                        gil::rgb8_pixel_t(0, 255, 0);
                }
            }
        }
    }
    gil::write_view("Fast_Output.jpeg", input_color_image_view,
                    gil::jpeg_tag{});
}

---

¹ (citation needed)