SIFT Keypoint Detection for Image in C++

Question 1

This is a follow-up question for 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 perform SIFT (Scale Invariant Feature Transform) keypoint detection (with get_potential_keypoint template function). There are several steps included:

Build several difference-of-Gaussian (DoG) images in several octaves (with generate_octave template function).
Find local maxima and minima pixels as keypoints in difference-of-Gaussian (DoG) images (with is_it_extremum template function and find_local_extrema template function).
Refine these keypoints' location using Taylor expansion (with keypoint_refinement template function).
Filter keypoints based on contrast and edge response (with keypoint_filtering template function).
Plot keypoints on the image (with draw_point, draw_points template functions).

The output example image (The white points are the detected keypoints):

OutputImage

The experimental implementation

is_it_extremum template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 /* is_it_extremum template function implementation
 input1, input2 and input3 are 3 * 3 images. If the center pixel (at location (1,1) ) of input2 is the largest / smallest one, 
 return true; otherwise, return false.
 Test: https://godbolt.org/z/Kb34EW5Yj
 */
 template<typename ElementT>
 constexpr static bool is_it_extremum(
 const Image<ElementT>& input1,
 const Image<ElementT>& input2,
 const Image<ElementT>& input3,
 const double threshold = 0.03)
 {
 if (input1.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input2.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input3.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input1.getWidth() != 3 || input1.getHeight() != 3)
 {
 throw std::runtime_error("Size error!");
 }
 if (input2.getWidth() != 3 || input2.getHeight() != 3)
 {
 throw std::runtime_error("Size error!");
 }
 if (input3.getWidth() != 3 || input3.getHeight() != 3)
 {
 throw std::runtime_error("Size error!");
 }
 auto center_pixel = input2.at(static_cast<std::size_t>(1), static_cast<std::size_t>(1));
 auto input2_img_data = input2.getImageData();
 input2_img_data.erase(input2_img_data.begin() + 4); // https://stackoverflow.com/a/875117/6667035
 if (std::abs(center_pixel) > threshold)
 {
 if (std::ranges::all_of(input1.getImageData(), [&](ElementT i) { return center_pixel > i; }) &&
 std::ranges::all_of(input3.getImageData(), [&](ElementT i) { return center_pixel > i; }) &&
 std::ranges::all_of(input2_img_data, [&](ElementT i) { return center_pixel > i; }))
 {
 return true;
 }
 if (std::ranges::all_of(input1.getImageData(), [&](ElementT i) { return center_pixel < i; }) &&
 std::ranges::all_of(input3.getImageData(), [&](ElementT i) { return center_pixel < i; }) &&
 std::ranges::all_of(input2_img_data, [&](ElementT i) { return center_pixel < i; }))
 {
 return true;
 }
 }
 return false;
 }
 }
}

keypoint_refinement template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // keypoint_refinement template function implementation
 // refine the given keypoint's location using Taylor expansion
 template<typename ElementT = double>
 constexpr static auto keypoint_refinement(
 const Image<ElementT>& input,
 const std::tuple<std::size_t, std::size_t>& point
 )
 {
 // Calculate the gradient at the keypoint (x, y)
 const ElementT first_derivative_x = (input.at(std::size_t{ 2 }, std::size_t{ 1 }) - input.at(std::size_t{ 0 }, std::size_t{ 1 })) / 2.0;
 const ElementT first_derivative_y = (input.at(std::size_t{ 1 }, std::size_t{ 2 }) - input.at(std::size_t{ 1 }, std::size_t{ 0 })) / 2.0;
 const ElementT second_derivative_x = (input.at(std::size_t{ 2 }, std::size_t{ 1 }) + input.at(std::size_t{ 0 }, std::size_t{ 1 }) - 2.0 * input.at(std::size_t{ 1 }, std::size_t{ 1 }));
 const ElementT second_derivative_y = (input.at(std::size_t{ 1 }, std::size_t{ 2 }) + input.at(std::size_t{ 1 }, std::size_t{ 0 }) - 2.0 * input.at(std::size_t{ 1 }, std::size_t{ 1 }));
 const ElementT second_derivative_xy = (input.at(std::size_t{ 2 }, std::size_t{ 2 }) - input.at(std::size_t{ 2 }, std::size_t{ 0 }) - input.at(std::size_t{ 0 }, std::size_t{ 2 }) + input.at(std::size_t{ 0 }, std::size_t{ 0 })) / 4.0;
 const ElementT A = -second_derivative_x / (second_derivative_x * second_derivative_y - second_derivative_xy * second_derivative_xy);
 const ElementT B = second_derivative_xy / (second_derivative_x * second_derivative_y - second_derivative_xy * second_derivative_xy);
 const ElementT C = B;
 const ElementT D = -second_derivative_y / (second_derivative_x * second_derivative_y - second_derivative_xy * second_derivative_xy);
 const ElementT offset_x = A * first_derivative_x + B * first_derivative_y;
 const ElementT offset_y = C * first_derivative_x + D * first_derivative_y;
 return std::make_tuple(
 static_cast<ElementT>(std::get<0>(point)) + offset_x,
 static_cast<ElementT>(std::get<1>(point)) + offset_y);
 }
 }
}

keypoint_filtering template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // keypoint_filtering template function implementation
 template<typename ElementT = double>
 constexpr static bool keypoint_filtering(
 const Image<ElementT>& input,
 const ElementT contrast_check_threshold = 8,
 const ElementT edge_response_threshold = 12.1)
 {
 // Calculate Hessian matrix at the keypoint (x, y)
 const ElementT second_derivative_x = (input.at(std::size_t{ 2 }, std::size_t{ 1 }) + input.at(std::size_t{ 0 }, std::size_t{ 1 }) - 2.0 * input.at(std::size_t{ 1 }, std::size_t{ 1 })); // D_{xx}
 const ElementT second_derivative_y = (input.at(std::size_t{ 1 }, std::size_t{ 2 }) + input.at(std::size_t{ 1 }, std::size_t{ 0 }) - 2.0 * input.at(std::size_t{ 1 }, std::size_t{ 1 })); // D_{yy}
 const ElementT second_derivative_xy = (input.at(std::size_t{ 2 }, std::size_t{ 2 }) - input.at(std::size_t{ 2 }, std::size_t{ 0 }) - input.at(std::size_t{ 0 }, std::size_t{ 2 }) + input.at(std::size_t{ 0 }, std::size_t{ 0 })) / 4.0; // D_{xy}
 const ElementT trace_Hessian_matrix = second_derivative_x + second_derivative_y;
 const ElementT determinant_Hessian_matrix = second_derivative_x * second_derivative_y - second_derivative_xy * second_derivative_xy;
 const ElementT principal_curvature_ratio = trace_Hessian_matrix * trace_Hessian_matrix / determinant_Hessian_matrix;
 if (input.at(std::size_t{ 1 }, std::size_t{ 1 }) <= contrast_check_threshold || principal_curvature_ratio >= edge_response_threshold)
 {
 return false;
 }
 return true;
 }
 }
}

mapping_point template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // mapping_point template function implementation
 template<class FloatingType = double>
 constexpr static std::tuple<std::size_t, std::size_t> mapping_point(
 const std::tuple<FloatingType, FloatingType>& input_location,
 const std::size_t input_width,
 const std::size_t input_height,
 const std::size_t target_width,
 const std::size_t target_height)
 {
 FloatingType width_percentage = static_cast<FloatingType>(std::get<0>(input_location)) / static_cast<FloatingType>(input_width);
 FloatingType height_percentage = static_cast<FloatingType>(std::get<1>(input_location)) / static_cast<FloatingType>(input_height);
 return std::make_tuple(
 static_cast<std::size_t>(width_percentage * static_cast<FloatingType>(target_width)),
 static_cast<std::size_t>(height_percentage * static_cast<FloatingType>(target_height))
 );
 }
 }
}

find_local_extrema template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // find_local_extrema template function implementation
 template<typename ElementT>
 constexpr static auto find_local_extrema(
 const Image<ElementT>& input1,
 const Image<ElementT>& input2,
 const Image<ElementT>& input3,
 const std::size_t octave_index,
 const std::size_t scale_index,
 const ElementT contrast_check_threshold = 8,
 const ElementT edge_response_threshold = 12.1)
 {
 if (input1.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input2.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input3.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 if (input1.getSize() != input2.getSize())
 {
 throw std::runtime_error("Size mismatched!");
 }
 if (input2.getSize() != input3.getSize())
 {
 throw std::runtime_error("Size mismatched!");
 }
 const int block_size = 3;
 std::vector<std::tuple<std::size_t, std::size_t, ElementT, ElementT>> output;
 auto width = input1.getWidth() - 1;
 auto height = input1.getHeight() - 1;
 #pragma omp parallel for collapse(2)
 for (std::size_t y = 1; y < height; ++y)
 {
 for (std::size_t x = 1; x < width; ++x)
 {
 auto subimage1 = subimage(input1, block_size, block_size, x, y);
 auto subimage2 = subimage(input2, block_size, block_size, x, y);
 auto subimage3 = subimage(input3, block_size, block_size, x, y);
 if (is_it_extremum(subimage1, subimage2, subimage3, contrast_check_threshold) && keypoint_filtering(subimage2, contrast_check_threshold, edge_response_threshold))
 {
 auto new_location = keypoint_refinement(subimage2, std::make_tuple(x, y));
 output.emplace_back(
 std::make_tuple(
 octave_index,
 scale_index,
 std::get<0>(new_location),
 std::get<1>(new_location)));
 }
 }
 }
 return output;
 }
 }
}

generate_octave template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // generate_octave template function implementation
 template<typename ElementT, typename SigmaT = double>
 requires(std::floating_point<SigmaT> || std::integral<SigmaT>)
 constexpr static auto generate_octave(
 const Image<ElementT>& input,
 std::size_t number_of_scale_levels = 5,
 SigmaT initial_sigma = 1.6,
 double k = 1.4142135623730950488016887242097)
 {
 std::vector<Image<ElementT>> difference_of_gaussian_images;
 difference_of_gaussian_images.resize(number_of_scale_levels - 1);
 #pragma omp parallel for
 for (int index = 0; index < number_of_scale_levels - 1; ++index)
 {
 difference_of_gaussian_images[index] = (difference_of_gaussian(input, initial_sigma * std::pow(k, index), initial_sigma * std::pow(k, index + 1)));
 }
 return difference_of_gaussian_images;
 }
 }
}

get_potential_keypoint template function implementation (in file image_operations.h)

namespace TinyDIP
{
 namespace SIFT_impl
 {
 // get_potential_keypoint template function implementation
 template<typename ElementT = double, typename SigmaT = double>
 requires( (std::floating_point<ElementT> || std::integral<ElementT>) &&
 (std::floating_point<SigmaT> || std::integral<SigmaT>))
 constexpr static auto get_potential_keypoint(
 const Image<ElementT>& input,
 std::size_t octaves_count = 4,
 std::size_t number_of_scale_levels = 5,
 SigmaT initial_sigma = 1.6,
 double k = 1.4142135623730950488016887242097,
 ElementT contrast_check_threshold = 8,
 ElementT edge_response_threshold = 12.1)
 {
 if (input.getDimensionality() != 2)
 {
 throw std::runtime_error("Unsupported dimension!");
 }
 // Generate octaves
 std::vector<std::vector<Image<ElementT>>> octaves;
 octaves.reserve(octaves_count);
 Image<ElementT> base_image = input;
 for (std::size_t octave_index = 0; octave_index < octaves_count; ++octave_index)
 {
 octaves.emplace_back(generate_octave(base_image, number_of_scale_levels, initial_sigma, k));
 base_image = copyResizeBicubic(
 base_image,
 static_cast<std::size_t>(static_cast<double>(base_image.getWidth()) / 2.0),
 static_cast<std::size_t>(static_cast<double>(base_image.getHeight()) / 2.0));
 }
 // Find potential KeyPoints
 /* KeyPoint structure: octave_index, scale_index + 1, location_x, location_y
 */
 std::vector<std::tuple<std::size_t, std::size_t, ElementT, ElementT>> keypoints;
 #pragma omp parallel for
 for (int octave_index = 0; octave_index < octaves_count; ++octave_index)
 {
 auto each_octave = octaves[octave_index];
 #pragma omp parallel for
 for (int scale_index = 0; scale_index < each_octave.size() - 2; ++scale_index)
 {
 /* if `append_range` function is supported
 keypoints.append_range(
 find_local_extrema(
 each_octave[scale_index],
 each_octave[scale_index + 1],
 each_octave[scale_index + 2],
 octave_index,
 scale_index,
 contrast_check_threshold,
 edge_response_threshold)
 );
 */
 for (auto&& element : find_local_extrema(
 each_octave[scale_index],
 each_octave[scale_index + 1],
 each_octave[scale_index + 2],
 octave_index,
 scale_index,
 contrast_check_threshold,
 edge_response_threshold))
 {
 keypoints.emplace_back(element);
 }
 }
 }
 // mapping keypoints in different scale to the original image
 std::vector<std::tuple<std::size_t, std::size_t>> mapped_keypoints;
 for (auto&& each_keypoint : keypoints)
 {
 auto input_width = octaves[std::get<0>(each_keypoint)][0].getWidth();
 auto input_height = octaves[std::get<0>(each_keypoint)][0].getHeight();
 mapped_keypoints.emplace_back(
 mapping_point(
 std::make_tuple(std::get<2>(each_keypoint), std::get<3>(each_keypoint)),
 input_width,
 input_height,
 input.getWidth(),
 input.getHeight()
 ));
 }
 return mapped_keypoints;
 }
 }
}

copyResizeBicubic template function implementation (in file image_operations.h)

namespace TinyDIP
{
 // copyResizeBicubic template function implementation
 template<class FloatingType = float, arithmetic ElementT>
 constexpr static Image<ElementT> copyResizeBicubic(const Image<ElementT>& image, const std::size_t width, const std::size_t height)
 {
 Image<ElementT> output(width, height);
 // get used to the C++ way of casting
 const auto ratiox = static_cast<FloatingType>(image.getWidth()) / static_cast<FloatingType>(width);
 const auto ratioy = static_cast<FloatingType>(image.getHeight()) / static_cast<FloatingType>(height);
 #pragma omp parallel for collapse(2)
 for (std::size_t y = 0; y < height; ++y)
 {
 for (std::size_t x = 0; x < width; ++x)
 {
 const FloatingType xMappingToOrigin = static_cast<FloatingType>(x) * ratiox;
 const FloatingType yMappingToOrigin = static_cast<FloatingType>(y) * ratioy;
 const FloatingType xMappingToOriginFloor = std::floor(xMappingToOrigin);
 const FloatingType yMappingToOriginFloor = std::floor(yMappingToOrigin);
 const FloatingType xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
 const FloatingType yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
 ElementT ndata[4 * 4];
 for (int ndatay = -1; ndatay <= 2; ++ndatay)
 {
 for (int ndatax = -1; ndatax <= 2; ++ndatax)
 {
 ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image.at(
 static_cast<std::size_t>(std::clamp(xMappingToOriginFloor + ndatax, static_cast<FloatingType>(0), static_cast<FloatingType>(image.getWidth()) - static_cast<FloatingType>(1))),
 static_cast<std::size_t>(std::clamp(yMappingToOriginFloor + ndatay, static_cast<FloatingType>(0), static_cast<FloatingType>(image.getHeight()) - static_cast<FloatingType>(1))));
 }
 }
 output.at(x, y) = bicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
 }
 }
 return output;
 }
 // copyResizeBicubic template function implementation for color image
 template<class FloatingType = double, class ElementT>
 requires (std::same_as<ElementT, RGB> || (std::same_as<ElementT, RGB_DOUBLE>) || (std::same_as<ElementT, HSV>) || is_MultiChannel<ElementT>::value)
 Image<ElementT> copyResizeBicubic(
 const Image<ElementT>& image,
 const std::size_t width,
 const std::size_t height)
 {
 return apply_each(image, [&](auto&& each_plane)
 {
 return copyResizeBicubic<FloatingType>(each_plane, width, height);
 });
 }
}

bicubicPolate template function implementation (in file image_operations.h)

namespace TinyDIP
{
 template<class InputT = float, class ElementT>
 constexpr static auto bicubicPolate(const ElementT* const ndata, const InputT fracx, const InputT fracy)
 {
 auto x1 = cubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
 auto x2 = cubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
 auto x3 = cubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
 auto x4 = cubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
 return std::clamp(
 static_cast<InputT>(cubicPolate(x1, x2, x3, x4, fracy)),
 static_cast<InputT>(std::numeric_limits<ElementT>::min()),
 static_cast<InputT>(std::numeric_limits<ElementT>::max()));
 }
}

cubicPolate template function implementation (in file image_operations.h)

namespace TinyDIP
{
 template<class InputT1, class InputT2>
 constexpr static auto cubicPolate(const InputT1 v0, const InputT1 v1, const InputT1 v2, const InputT1 v3, const InputT2 frac)
 {
 auto A = (v3-v2)-(v0-v1);
 auto B = (v0-v1)-A;
 auto C = v2-v0;
 auto D = v1;
 return D + frac * (C + frac * (B + frac * A));
 }
}

draw_point template function implementation (in file image_operations.h)

namespace TinyDIP
{
 // draw_point template function implementation
 template<typename ElementT>
 constexpr static auto draw_point(
 const Image<ElementT>& input,
 std::tuple<std::size_t, std::size_t> point,
 ElementT draw_value = ElementT{},
 std::size_t radius = 3
 )
 {
 auto point_x = std::get<0>(point);
 auto point_y = std::get<1>(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::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))
 {
 output.at_without_boundary_check(x, y) = draw_value;
 }
 }
 }
 return output;
 }
}

draw_points template function implementation (in file image_operations.h)

namespace TinyDIP
{
 // draw_points template function implementation
 template<typename ElementT>
 constexpr static auto draw_points(
 const Image<ElementT>& input,
 std::vector<std::tuple<std::size_t, std::size_t>> points,
 ElementT draw_value = ElementT{},
 std::size_t radius = 3
 )
 {
 auto output = input;
 for (auto&& each_point : points)
 {
 output = draw_point(output, each_point, draw_value, radius);
 }
 return output;
 }
}

apply_each template function implementation (in file image_operations.h)

namespace TinyDIP
{
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<RGB>& input, F operation, Args&&... args)
 {
 auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input), args...); });
 auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input), args...); });
 auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input), args...); });
 if constexpr ((std::is_same<Image<GrayScale>, decltype(Rplane.get())>::value) and
 (std::is_same<Image<GrayScale>, decltype(Gplane.get())>::value) and
 (std::is_same<Image<GrayScale>, decltype(Bplane.get())>::value))
 {
 return constructRGB(Rplane.get(), Gplane.get(), Bplane.get());
 }
 else if constexpr ((std::is_same<Image<double>, decltype(Rplane.get())>::value) and
 (std::is_same<Image<double>, decltype(Gplane.get())>::value) and
 (std::is_same<Image<double>, decltype(Bplane.get())>::value))
 {
 return constructRGBDOUBLE(Rplane.get(), Gplane.get(), Bplane.get());
 }
 else
 {
 return constructMultiChannel(Rplane.get(), Gplane.get(), Bplane.get());
 }
 }
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<RGB_DOUBLE>& input, F operation, Args&&... args)
 {
 auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input), args...); });
 auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input), args...); });
 auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input), args...); });
 if constexpr ((std::is_same<Image<double>, decltype(Rplane.get())>::value) and
 (std::is_same<Image<double>, decltype(Gplane.get())>::value) and
 (std::is_same<Image<double>, decltype(Bplane.get())>::value))
 {
 return constructRGBDOUBLE(Rplane.get(), Gplane.get(), Bplane.get());
 }
 else
 {
 return constructMultiChannel(Rplane.get(), Gplane.get(), Bplane.get());
 }
 }
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<HSV>& input, F operation, Args&&... args)
 {
 auto Hplane = std::async(std::launch::async, [&] { return std::invoke(operation, getHplane(input), args...); });
 auto Splane = std::async(std::launch::async, [&] { return std::invoke(operation, getSplane(input), args...); });
 auto Vplane = std::async(std::launch::async, [&] { return std::invoke(operation, getVplane(input), args...); });
 if constexpr ((std::is_same<Image<double>, decltype(Hplane.get())>::value) and
 (std::is_same<Image<double>, decltype(Splane.get())>::value) and
 (std::is_same<Image<double>, decltype(Vplane.get())>::value))
 {
 return constructHSV(Hplane.get(), Splane.get(), Vplane.get());
 }
 else
 {
 return constructMultiChannel(Hplane.get(), Splane.get(), Vplane.get());
 }
 }
 // apply_each_impl template function implementation
 // Helper function implementation using index_sequence
 template<class ElementT, std::size_t Channels, class F, class... Args, std::size_t... Is>
 constexpr static auto apply_each_impl(const Image<MultiChannel<ElementT, Channels>>& input, F operation, Args&&... args, std::index_sequence<Is...>) {
 return constructMultiChannel(
 std::async(std::launch::async, [&] {
 return std::invoke(operation, getPlane(input, Is), std::forward<Args>(args)...);
 }).get()...
 );
 }
 // apply_each template function implementation
 template<class ElementT, std::size_t Channels, class F, class... Args>
 constexpr static auto apply_each(const Image<MultiChannel<ElementT, Channels>>& input, F operation, Args&&... args)
 {
 return apply_each_impl(input, operation, std::forward<Args>(args)..., std::make_index_sequence<Channels>{});
 }
 // apply_each template function implementation
 template<class ElementT, class F, class... Args>
 constexpr static auto apply_each(const Image<MultiChannel<ElementT, 3>>& input, F operation, Args&&... args)
 {
 return constructMultiChannel(
 std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input, 0), args...); }).get(),
 std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input, 1), args...); }).get(),
 std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input, 2), args...); }).get()
 );
 }
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<RGB>& input1, const Image<RGB>& input2, F operation, Args&&... args)
 {
 auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input1), getRplane(input2), args...); });
 auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input1), getGplane(input2), args...); });
 auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input1), getBplane(input2), args...); });
 return constructRGB(Rplane.get(), Gplane.get(), Bplane.get());
 }
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<RGB_DOUBLE>& input1, const Image<RGB_DOUBLE>& input2, F operation, Args&&... args)
 {
 auto Rplane = std::async(std::launch::async, [&] { return std::invoke(operation, getRplane(input1), getRplane(input2), args...); });
 auto Gplane = std::async(std::launch::async, [&] { return std::invoke(operation, getGplane(input1), getGplane(input2), args...); });
 auto Bplane = std::async(std::launch::async, [&] { return std::invoke(operation, getBplane(input1), getBplane(input2), args...); });
 return constructRGBDOUBLE(Rplane.get(), Gplane.get(), Bplane.get());
 }
 // apply_each template function implementation
 template<class F, class... Args>
 constexpr static auto apply_each(const Image<HSV> input1, const Image<HSV> input2, F operation, Args&&... args)
 {
 auto Hplane = std::async(std::launch::async, [&] { return std::invoke(operation, getHplane(input1), getHplane(input2), args...); });
 auto Splane = std::async(std::launch::async, [&] { return std::invoke(operation, getSplane(input1), getSplane(input2), args...); });
 auto Vplane = std::async(std::launch::async, [&] { return std::invoke(operation, getVplane(input1), getVplane(input2), args...); });
 return constructHSV(Hplane.get(), Splane.get(), Vplane.get());
 }
 // apply_each template function implementation
 template<class ElementT, class F, class... Args>
 constexpr static auto apply_each(const Image<MultiChannel<ElementT, 3>> input1, const Image<MultiChannel<ElementT, 3>> input2, F operation, Args&&... args)
 {
 auto plane1 = std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input1, 0), getPlane(input2, 0), args...); });
 auto plane2 = std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input1, 1), getPlane(input2, 1), args...); });
 auto plane3 = std::async(std::launch::async, [&] { return std::invoke(operation, getPlane(input1, 2), getPlane(input2, 2), args...); });
 return constructMultiChannel(plane1.get(), plane2.get(), plane3.get());
 }
}

getRplane function implementation (in file image_operations.h)

namespace TinyDIP
{
 // getRplane function implementation
 constexpr static auto getRplane(const Image<RGB>& input)
 {
 return getPlane(input, 0);
 }
 // getRplane function implementation
 constexpr static auto getRplane(const Image<RGB_DOUBLE>& input)
 {
 return getPlane(input, 0);
 }
}

getGplane function implementation (in file image_operations.h)

namespace TinyDIP
{
 // getGplane function implementation
 constexpr static auto getGplane(const Image<RGB>& input)
 {
 return getPlane(input, 1);
 }
 // getGplane function implementation
 constexpr static auto getGplane(const Image<RGB_DOUBLE>& input)
 {
 return getPlane(input, 1);
 }
}

getBplane function implementation (in file image_operations.h)

namespace TinyDIP
{
 // getBplane function implementation
 constexpr static auto getBplane(const Image<RGB>& input)
 {
 return getPlane(input, 2);
 }
 // getBplane function implementation
 constexpr static auto getBplane(const Image<RGB_DOUBLE>& input)
 {
 return getPlane(input, 2);
 }
}

constructRGB function implementation (in file image_operations.h)

namespace TinyDIP
{
 // constructRGB template function implementation
 template<typename OutputT = RGB>
 constexpr static auto constructRGB(const Image<GrayScale>& r, const Image<GrayScale>& g, const Image<GrayScale>& b)
 {
 check_size_same(r, g);
 check_size_same(g, b);
 auto image_data_r = r.getImageData();
 auto image_data_g = g.getImageData();
 auto image_data_b = b.getImageData();
 std::vector<OutputT> new_data;
 new_data.resize(r.count());
 #pragma omp parallel for
 for (std::size_t index = 0; index < r.count(); ++index)
 {
 OutputT rgb { image_data_r[index],
 image_data_g[index],
 image_data_b[index]};
 new_data[index] = rgb;
 }
 Image<OutputT> output(new_data, r.getSize());
 return output;
 }
}

The usage of get_potential_keypoint 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 SIFT_keypoints = TinyDIP::SIFT_impl::get_potential_keypoint(
 TinyDIP::getVplane(TinyDIP::rgb2hsv(bmp1)));
 std::cout << "SIFT_keypoints count = " << SIFT_keypoints.size() << "\n";
 RGB rgb{ 255, 255, 255 };
 bmp1 = TinyDIP::draw_points(bmp1, SIFT_keypoints, rgb);
 TinyDIP::bmp_write("test20240809", 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() << '\n';
 return EXIT_SUCCESS;
}

The output of the test code above:

Width of the input image: 512
Height of the input image: 512
Size of the input image(Byte): 786432
SIFT_keypoints count = 430
Computation finished at Fri Aug 9 15:04:24 2024
elapsed time: 60.1763

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

Which question it is a follow-up to?

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 perform SIFT (Scale Invariant Feature Transform) keypoint detection in this post.
Why a new review is being asked for?

Please review the implementation of get_potential_keypoint template function and the related functions.

Question 2

Avoid repeating yourself

Code duplication increases the maintenance burden and the chance of introducing bugs. Try to avoid it as much as possible. For exampe, in is_it_extremum(), you can write this instead:

template<typename ElementT>
constexpr static bool is_it_extremum(
 const Image<ElementT>& input1,
 const Image<ElementT>& input2,
 const Image<ElementT>& input3,
 const double threshold = 0.03)
{
 for (auto* input: {&input1, &input2, &input3}) {
 if (input->getDimensionality() != 2) {
 throw std::runtime_error("Unsupported dimensionality");
 }
 if (input->getWidth() != 3 || input->getHeight() != 3) {
 throw std::runtime_error("Unsupported size");
 }
 }
 ...
 if (std::abs(center_pixel) <= threshold) {
 return false;
 }
 
 auto image_datas = {&input1.getImageData(), &input2.img_data, &input3.getImageData()};
 if (std::ranges::all_of(imageDatas, [&](auto* image_data) {
 return std::ranges::all_of(*image_data, [&](auto i) {
 return center_pixel > i;
 });
 })
 {
 return true;
 }
 ...
 return false;
}

The nested std::ranges::all_of() are a bit uncommon, and might be a bit hard to understand at first. Although perhaps you could create a recursive_all_of()? That way you can write:

auto image_datas = {&input1.getImageData(), &input2.img_data, &input3.getImageData()};
if (recursive_all_of(image_datas, [&](auto i) { return center_pixel > i; })) {
 return true;
}
if (recursive_all_of(image_datas, [&](auto i) { return center_pixel < i; })) {
 return true;
}

I used initializer lists here to quickly generate something that can be iterated over. You can't store references in them (unless you use std::reference_wrapper of course), and the other drawback is that it only works for objects of the same type. If you want to "loop" over multiple heterogenous objects, you can std::invoke(), pack expansion and immediately invoked lambdas:

std::invoke([](auto&... inputs) {
 ([](auto& inputs) {
 if (inputs->getDimensionality != 2) {
 throw std::runtime_error("Unsupported dimensionality");
 }
 ...
 }, ...);
}, input1, input2, input3);

It is a bit unfortunate that C++ doesn't have a universal syntax to loop over both homogenous and heterogenous collections of objects.

There are more opportunities for removing code duplication in your code, not just repeated lines, but also repeated patterns:

Avoid needing to `static_cast<std::size>()` everywhere

Instead of:

auto center_pixel = input2.at(static_cast<std::size_t>(1), static_cast<std::size_t>(1));

Write the following:

auto center_pixel = input2.at(1, 1);

If that doesn't compile, fix at() so it does. Having to write static_cast<std::size_t>() everywhere for constants is very annoying.

Create a `struct Point`

Instead of passing coordinates around in a std::tuple<>, create a struct Point. Not only does that save a lot of typing (no more need for lengthy declarations, for std::make_tuple or for std::get<>()), it is also more type-safe, as it cannot be confused with possible other uses of a tuple of two std::size_ts. Bonus points if you allow math operations on Points. And probably even better is just using a math library that provides a 2D vector type.

Try to use std::tuple only for generic code where you don't know the size of the tuple up front.

Don't write out mathematical constants

I see k = 1.4142135623730950488016887242097 a few times. Is that enough digits? Did you make a mistake in one of the digits? It's hard to tell. Why not write k = std::sqrt(2.0) or k = std::numbers::sqrt2_v<double>? Even that last one is less typing.

Explain magic numbers

There are a lot of magic numbers in your code. For example, in get_potential_keypoint() there is 4, 5, 1.6, \$\sqrt{2}\$, 8 and 12.1 for various parameters. Why those values? At the very least, add some comments explaining where you got these numbers from, preferably with a link to some paper or article.

Naming

Most of your code uses good naming practices. However, in some functions you seem to have given up. The most egregious example is cubicPolate():

cubicPolate(): it's doing some cubic interpolation, so why not write cubicInterpolate()? But wait, most of your functions use snake_case instead of camelCase, so it should be cubic_interpolate().
v0...v3, A...D: these are some very generic names. For implementations of mathematical formulas, it makes sense though, but if you don't already know the formula for cubic interpolation by heart, the code doesn't make any sense. The solution here is to link to a paper or article that explains the formula, for example: https://en.wikipedia.org/wiki/Monotone_cubic_interpolation. But then also use the same names as used in that paper or article.
InputT1, InputT2: it's not very clear what those types are for. Maybe VT and FracT, so it's clear they apply to v0...v3 and frac respectively?

G. Sliepen G. Sliepen 68.7k3 gold badges74 silver badges179 bronze badges · Accepted Answer · 2025-04-11 19:39:44Z

Avoid repeating yourself

Code duplication increases the maintenance burden and the chance of introducing bugs. Try to avoid it as much as possible. For exampe, in is_it_extremum(), you can write this instead:

template<typename ElementT>
constexpr static bool is_it_extremum(
 const Image<ElementT>& input1,
 const Image<ElementT>& input2,
 const Image<ElementT>& input3,
 const double threshold = 0.03)
{
 for (auto* input: {&input1, &input2, &input3}) {
 if (input->getDimensionality() != 2) {
 throw std::runtime_error("Unsupported dimensionality");
 }
 if (input->getWidth() != 3 || input->getHeight() != 3) {
 throw std::runtime_error("Unsupported size");
 }
 }
 ...
 if (std::abs(center_pixel) <= threshold) {
 return false;
 }
 
 auto image_datas = {&input1.getImageData(), &input2.img_data, &input3.getImageData()};
 if (std::ranges::all_of(imageDatas, [&](auto* image_data) {
 return std::ranges::all_of(*image_data, [&](auto i) {
 return center_pixel > i;
 });
 })
 {
 return true;
 }
 ...
 return false;
}

The nested std::ranges::all_of() are a bit uncommon, and might be a bit hard to understand at first. Although perhaps you could create a recursive_all_of()? That way you can write:

auto image_datas = {&input1.getImageData(), &input2.img_data, &input3.getImageData()};
if (recursive_all_of(image_datas, [&](auto i) { return center_pixel > i; })) {
 return true;
}
if (recursive_all_of(image_datas, [&](auto i) { return center_pixel < i; })) {
 return true;
}

I used initializer lists here to quickly generate something that can be iterated over. You can't store references in them (unless you use std::reference_wrapper of course), and the other drawback is that it only works for objects of the same type. If you want to "loop" over multiple heterogenous objects, you can std::invoke(), pack expansion and immediately invoked lambdas:

std::invoke([](auto&... inputs) {
 ([](auto& inputs) {
 if (inputs->getDimensionality != 2) {
 throw std::runtime_error("Unsupported dimensionality");
 }
 ...
 }, ...);
}, input1, input2, input3);

It is a bit unfortunate that C++ doesn't have a universal syntax to loop over both homogenous and heterogenous collections of objects.

There are more opportunities for removing code duplication in your code, not just repeated lines, but also repeated patterns:

Avoid needing to `static_cast<std::size>()` everywhere

Instead of:

auto center_pixel = input2.at(static_cast<std::size_t>(1), static_cast<std::size_t>(1));

Write the following:

auto center_pixel = input2.at(1, 1);

If that doesn't compile, fix at() so it does. Having to write static_cast<std::size_t>() everywhere for constants is very annoying.

Create a `struct Point`

Instead of passing coordinates around in a std::tuple<>, create a struct Point. Not only does that save a lot of typing (no more need for lengthy declarations, for std::make_tuple or for std::get<>()), it is also more type-safe, as it cannot be confused with possible other uses of a tuple of two std::size_ts. Bonus points if you allow math operations on Points. And probably even better is just using a math library that provides a 2D vector type.

Try to use std::tuple only for generic code where you don't know the size of the tuple up front.

Don't write out mathematical constants

I see k = 1.4142135623730950488016887242097 a few times. Is that enough digits? Did you make a mistake in one of the digits? It's hard to tell. Why not write k = std::sqrt(2.0) or k = std::numbers::sqrt2_v<double>? Even that last one is less typing.

Explain magic numbers

There are a lot of magic numbers in your code. For example, in get_potential_keypoint() there is 4, 5, 1.6, \$\sqrt{2}\$, 8 and 12.1 for various parameters. Why those values? At the very least, add some comments explaining where you got these numbers from, preferably with a link to some paper or article.

Naming

Most of your code uses good naming practices. However, in some functions you seem to have given up. The most egregious example is cubicPolate():

cubicPolate(): it's doing some cubic interpolation, so why not write cubicInterpolate()? But wait, most of your functions use snake_case instead of camelCase, so it should be cubic_interpolate().
v0...v3, A...D: these are some very generic names. For implementations of mathematical formulas, it makes sense though, but if you don't already know the formula for cubic interpolation by heart, the code doesn't make any sense. The solution here is to link to a paper or article that explains the formula, for example: https://en.wikipedia.org/wiki/Monotone_cubic_interpolation. But then also use the same names as used in that paper or article.
InputT1, InputT2: it's not very clear what those types are for. Maybe VT and FracT, so it's clear they apply to v0...v3 and frac respectively?

Stack Exchange Network

SIFT Keypoint Detection for Image in C++

1 Answer 1

Avoid repeating yourself

Avoid needing to `static_cast<std::size>()` everywhere

Create a `struct Point`

Don't write out mathematical constants

Explain magic numbers

Naming

Your Answer

Sign up or log in

Post as a guest

Post as a guest

Linked

Hot Network Questions

SIFT Keypoint Detection for Image in C++

1 Answer 1

Avoid repeating yourself

Avoid needing to static_cast<std::size>() everywhere

Create a struct Point

Don't write out mathematical constants

Explain magic numbers

Naming

Your Answer

Sign up or log in

Post as a guest

Post as a guest

Linked

Related

Hot Network Questions

Avoid needing to `static_cast<std::size>()` everywhere

Create a `struct Point`