1 #include <boost/gil/image.hpp>
2 #include <boost/gil/image_view.hpp>
3 #include <boost/gil/image_processing/numeric.hpp>
4 #include <boost/gil/image_processing/hessian.hpp>
5 #include <boost/gil/extension/io/png.hpp>
12 namespace gil
= boost::gil
;
14 // some images might produce artifacts
15 // when converted to grayscale,
16 // which was previously observed on
17 // canny edge detector for test input
18 // used for this example.
19 // the algorithm here follows sRGB gamma definition
20 // taken from here (luminance calculation):
21 // https://en.wikipedia.org/wiki/Grayscale
22 gil::gray8_image_t
to_grayscale(gil::rgb8_view_t original
)
24 gil::gray8_image_t
output_image(original
.dimensions());
25 auto output
= gil::view(output_image
);
26 constexpr double max_channel_intensity
= (std::numeric_limits
<std::uint8_t>::max
)();
27 for (long int y
= 0; y
< original
.height(); ++y
)
29 for (long int x
= 0; x
< original
.width(); ++x
)
31 // scale the values into range [0, 1] and calculate linear intensity
32 auto& p
= original(x
, y
);
33 double red_intensity
= p
.at(std::integral_constant
<int, 0>{})
34 / max_channel_intensity
;
35 double green_intensity
= p
.at(std::integral_constant
<int, 1>{})
36 / max_channel_intensity
;
37 double blue_intensity
= p
.at(std::integral_constant
<int, 2>{})
38 / max_channel_intensity
;
39 auto linear_luminosity
= 0.2126 * red_intensity
40 + 0.7152 * green_intensity
41 + 0.0722 * blue_intensity
;
43 // perform gamma adjustment
44 double gamma_compressed_luminosity
= 0;
45 if (linear_luminosity
< 0.0031308)
47 gamma_compressed_luminosity
= linear_luminosity
* 12.92;
50 gamma_compressed_luminosity
= 1.055 * std::pow(linear_luminosity
, 1 / 2.4) - 0.055;
53 // since now it is scaled, descale it back
54 output(x
, y
) = gamma_compressed_luminosity
* max_channel_intensity
;
61 void apply_gaussian_blur(gil::gray8_view_t input_view
, gil::gray8_view_t output_view
)
63 constexpr static auto filter_height
= 5ull;
64 constexpr static auto filter_width
= 5ull;
65 constexpr static double filter
[filter_height
][filter_width
] =
73 constexpr double factor
= 1.0 / 159;
74 constexpr double bias
= 0.0;
76 const auto height
= input_view
.height();
77 const auto width
= input_view
.width();
78 for (std::ptrdiff_t x
= 0; x
< width
; ++x
)
80 for (std::ptrdiff_t y
= 0; y
< height
; ++y
)
82 double intensity
= 0.0;
83 for (std::ptrdiff_t filter_y
= 0; filter_y
< filter_height
; ++filter_y
)
85 for (std::ptrdiff_t filter_x
= 0; filter_x
< filter_width
; ++filter_x
)
87 int image_x
= x
- filter_width
/ 2 + filter_x
;
88 int image_y
= y
- filter_height
/ 2 + filter_y
;
89 if (image_x
>= input_view
.width() || image_x
< 0 ||
90 image_y
>= input_view
.height() || image_y
< 0)
94 const auto& pixel
= input_view(image_x
, image_y
);
95 intensity
+= pixel
.at(std::integral_constant
<int, 0>{})
96 * filter
[filter_y
][filter_x
];
99 auto& pixel
= output_view(gil::point_t(x
, y
));
100 pixel
= (std::min
)((std::max
)(int(factor
* intensity
+ bias
), 0), 255);
106 std::vector
<gil::point_t
> suppress(
107 gil::gray32f_view_t harris_response
,
108 double harris_response_threshold
)
110 std::vector
<gil::point_t
> corner_points
;
111 for (gil::gray32f_view_t::coord_t y
= 1; y
< harris_response
.height() - 1; ++y
)
113 for (gil::gray32f_view_t::coord_t x
= 1; x
< harris_response
.width() - 1; ++x
)
115 auto value
= [](gil::gray32f_pixel_t pixel
) {
116 return pixel
.at(std::integral_constant
<int, 0>{});
119 value(harris_response(x
- 1, y
- 1)),
120 value(harris_response(x
, y
- 1)),
121 value(harris_response(x
+ 1, y
- 1)),
122 value(harris_response(x
- 1, y
)),
123 value(harris_response(x
, y
)),
124 value(harris_response(x
+ 1, y
)),
125 value(harris_response(x
- 1, y
+ 1)),
126 value(harris_response(x
, y
+ 1)),
127 value(harris_response(x
+ 1, y
+ 1))
130 auto maxima
= *std::max_element(
133 [](double lhs
, double rhs
)
139 if (maxima
== value(harris_response(x
, y
))
140 && std::count(values
, values
+ 9, maxima
) == 1
141 && maxima
>= harris_response_threshold
)
143 corner_points
.emplace_back(x
, y
);
148 return corner_points
;
151 int main(int argc
, char* argv
[]) {
154 std::cout
<< "usage: " << argv
[0] << " <input.png> <odd-window-size>"
155 " <hessian-response-threshold> <output.png>\n";
159 std::size_t window_size
= std::stoul(argv
[2]);
160 long hessian_determinant_threshold
= std::stol(argv
[3]);
162 gil::rgb8_image_t input_image
;
164 gil::read_image(argv
[1], input_image
, gil::png_tag
{});
166 auto input_view
= gil::view(input_image
);
167 auto grayscaled
= to_grayscale(input_view
);
168 gil::gray8_image_t
smoothed_image(grayscaled
.dimensions());
169 auto smoothed
= gil::view(smoothed_image
);
170 apply_gaussian_blur(gil::view(grayscaled
), smoothed
);
171 gil::gray16s_image_t
x_gradient_image(grayscaled
.dimensions());
172 gil::gray16s_image_t
y_gradient_image(grayscaled
.dimensions());
174 auto x_gradient
= gil::view(x_gradient_image
);
175 auto y_gradient
= gil::view(y_gradient_image
);
176 auto scharr_x
= gil::generate_dx_scharr();
177 gil::detail::convolve_2d(smoothed
, scharr_x
, x_gradient
);
178 auto scharr_y
= gil::generate_dy_scharr();
179 gil::detail::convolve_2d(smoothed
, scharr_y
, y_gradient
);
181 gil::gray32f_image_t
m11(x_gradient
.dimensions());
182 gil::gray32f_image_t
m12_21(x_gradient
.dimensions());
183 gil::gray32f_image_t
m22(x_gradient
.dimensions());
184 gil::compute_hessian_entries(
192 gil::gray32f_image_t
hessian_response(x_gradient
.dimensions());
193 auto gaussian_kernel
= gil::generate_gaussian_kernel(window_size
, 0.84089642);
194 gil::compute_hessian_responses(
199 gil::view(hessian_response
)
202 auto corner_points
= suppress(gil::view(hessian_response
), hessian_determinant_threshold
);
203 for (auto point
: corner_points
) {
204 input_view(point
) = gil::rgb8_pixel_t(0, 0, 0);
205 input_view(point
).at(std::integral_constant
<int, 1>{}) = 255;
207 gil::write_view(argv
[4], input_view
, gil::png_tag
{});