1
1
#include < boost/gil/image.hpp>
2
2
#include < boost/gil/image_view.hpp>
3
+ #include < boost/gil/extension/numeric/matrix.hpp>
4
+ #include < boost/gil/extension/numeric/arithmetics.hpp>
3
5
#include < boost/gil/image_processing/numeric.hpp>
4
6
#include < boost/gil/image_processing/hessian.hpp>
5
7
#include < boost/gil/extension/io/png.hpp>
11
13
12
14
namespace gil = boost::gil;
13
15
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)
23
- {
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)
28
- {
29
- for (long int x = 0 ; x < original.width (); ++x)
30
- {
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;
42
-
43
- // perform gamma adjustment
44
- double gamma_compressed_luminosity = 0 ;
45
- if (linear_luminosity < 0.0031308 )
46
- {
47
- gamma_compressed_luminosity = linear_luminosity * 12.92 ;
48
- } else
49
- {
50
- gamma_compressed_luminosity = 1.055 * std::pow (linear_luminosity, 1 / 2.4 ) - 0.055 ;
51
- }
52
-
53
- // since now it is scaled, descale it back
54
- output (x, y) = gamma_compressed_luminosity * max_channel_intensity;
55
- }
56
- }
57
-
58
- return output_image;
59
- }
60
-
61
- void apply_gaussian_blur (gil::gray8_view_t input_view, gil::gray8_view_t output_view)
62
- {
63
- constexpr static std::ptrdiff_t filter_height = 5ull ;
64
- constexpr static std::ptrdiff_t filter_width = 5ull ;
65
- constexpr static double filter[filter_height][filter_width] =
66
- {
67
- { 2 , 4 , 6 , 4 , 2 },
68
- { 4 , 9 , 12 , 9 , 4 },
69
- { 5 , 12 , 15 , 12 , 5 },
70
- { 4 , 9 , 12 , 9 , 4 },
71
- { 2 , 4 , 5 , 4 , 2 }
72
- };
73
- constexpr double factor = 1.0 / 159 ;
74
- constexpr double bias = 0.0 ;
75
-
76
- const auto height = input_view.height ();
77
- const auto width = input_view.width ();
78
- for (std::ptrdiff_t x = 0 ; x < width; ++x)
79
- {
80
- for (std::ptrdiff_t y = 0 ; y < height; ++y)
81
- {
82
- double intensity = 0.0 ;
83
- for (std::ptrdiff_t filter_y = 0 ; filter_y < filter_height; ++filter_y)
84
- {
85
- for (std::ptrdiff_t filter_x = 0 ; filter_x < filter_width; ++filter_x)
86
- {
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 )
91
- {
92
- continue ;
93
- }
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];
97
- }
98
- }
99
- auto & pixel = output_view (gil::point_t (x, y));
100
- pixel = (std::min)((std::max)(int (factor * intensity + bias), 0 ), 255 );
101
- }
102
-
103
- }
104
- }
105
-
106
16
std::vector<gil::point_t > suppress (
107
- gil::gray32f_view_t harris_response,
17
+ gil::matrix harris_response,
108
18
double harris_response_threshold)
109
19
{
110
20
std::vector<gil::point_t > corner_points;
111
21
for (gil::gray32f_view_t ::coord_t y = 1 ; y < harris_response.height () - 1 ; ++y)
112
22
{
113
23
for (gil::gray32f_view_t ::coord_t x = 1 ; x < harris_response.width () - 1 ; ++x)
114
24
{
115
- auto value = [](gil::gray32f_pixel_t pixel) {
116
- return pixel.at (std::integral_constant<int , 0 >{});
117
- };
118
- double values[9 ] = {
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 ))
25
+ gil::elem_t values[9 ] = {
26
+ harris_response (x - 1 , y - 1 ),
27
+ harris_response (x, y - 1 ),
28
+ harris_response (x + 1 , y - 1 ),
29
+ harris_response (x - 1 , y),
30
+ harris_response (x, y),
31
+ harris_response (x + 1 , y),
32
+ harris_response (x - 1 , y + 1 ),
33
+ harris_response (x, y + 1 ),
34
+ harris_response (x + 1 , y + 1 ),
128
35
};
129
36
130
37
auto maxima = *std::max_element (
131
38
values,
132
39
values + 9 ,
133
- [](double lhs, double rhs)
40
+ [](gil:: elemc_ref_t lhs, gil:: elemc_ref_t rhs)
134
41
{
135
- return lhs < rhs;
42
+ return std::abs ( lhs[ 0 ]) < std::abs ( rhs[ 0 ]) ;
136
43
}
137
44
);
138
45
139
- if (maxima == value (harris_response (x, y))
140
- && std::count (values, values + 9 , maxima) == 1
141
- && maxima >= harris_response_threshold)
46
+ if (maxima[0 ] == harris_response (x, y)[0 ]
47
+ && std::abs (maxima[0 ]) >= harris_response_threshold)
142
48
{
143
49
corner_points.emplace_back (x, y);
144
50
}
@@ -148,6 +54,13 @@ std::vector<gil::point_t> suppress(
148
54
return corner_points;
149
55
}
150
56
57
+ std::pair<float , float > get_min_max_elements (const gil::matrix& m) {
58
+ auto min_element = *std::min_element (m.begin (), m.end ());
59
+ auto max_element = *std::max_element (m.begin (), m.end ());
60
+
61
+ return {min_element, max_element};
62
+ }
63
+
151
64
int main (int argc, char * argv[]) {
152
65
if (argc != 5 )
153
66
{
@@ -157,39 +70,48 @@ int main(int argc, char* argv[]) {
157
70
}
158
71
159
72
std::size_t window_size = std::stoul (argv[2 ]);
160
- long hessian_determinant_threshold = std::stol (argv[3 ]);
73
+ double hessian_determinant_threshold = std::stod (argv[3 ]);
161
74
162
75
gil::rgb8_image_t input_image;
163
76
164
77
gil::read_image (argv[1 ], input_image, gil::png_tag{});
165
78
166
79
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 ());
173
-
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);
180
-
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 ());
80
+ gil::gray8_image_t grayscaled_image (input_image.dimensions ());
81
+ auto grayscaled = gil::view (grayscaled_image);
82
+ gil::copy_and_convert_pixels (input_view, grayscaled);
83
+ gil::write_view (" grayscale.png" , grayscaled, gil::png_tag{});
84
+
85
+ gil::matrix_storage input_matrix_storage (grayscaled_image.dimensions ());
86
+ auto input_matrix = gil::view (input_matrix_storage);
87
+ gil::value_cast (grayscaled, input_matrix);
88
+
89
+ gil::matrix_storage grayscaled_matrix_storage (grayscaled.dimensions ());
90
+ auto grayscaled_matrix = gil::view (grayscaled_matrix_storage);
91
+ gil::value_cast (grayscaled, grayscaled_matrix);
92
+
93
+ gil::matrix_storage smoothed_matrix (grayscaled.dimensions ());
94
+ auto smoothed = gil::view (smoothed_matrix);
95
+ gil::convolve_2d<gil::elem_t >(grayscaled_matrix, gil::generate_gaussian_kernel (3 , 1.0 ), smoothed);
96
+ gil::matrix_storage dx_matrix (grayscaled.dimensions ());
97
+ auto dx = gil::view (dx_matrix);
98
+ gil::matrix_storage dy_matrix (grayscaled.dimensions ());
99
+ auto dy = gil::view (dy_matrix);
100
+ gil::convolve_2d<gil::elem_t >(smoothed, gil::generate_dx_scharr (), dx);
101
+ gil::convolve_2d<gil::elem_t >(smoothed, gil::generate_dy_sobel (), dy);
102
+
103
+ gil::matrix_storage m11 (dx.dimensions ());
104
+ gil::matrix_storage m12_21 (dx.dimensions ());
105
+ gil::matrix_storage m22 (dy.dimensions ());
184
106
gil::compute_hessian_entries (
185
- x_gradient ,
186
- y_gradient ,
107
+ dx ,
108
+ dy ,
187
109
gil::view (m11),
188
110
gil::view (m12_21),
189
111
gil::view (m22)
190
112
);
191
113
192
- gil::gray32f_image_t hessian_response (x_gradient .dimensions ());
114
+ gil::matrix_storage hessian_response (dx .dimensions ());
193
115
auto gaussian_kernel = gil::generate_gaussian_kernel (window_size, 0.84089642 );
194
116
gil::compute_hessian_responses (
195
117
gil::view (m11),
@@ -199,7 +121,12 @@ int main(int argc, char* argv[]) {
199
121
gil::view (hessian_response)
200
122
);
201
123
124
+ auto minmax_elements = get_min_max_elements (gil::view (hessian_response));
125
+ std::cout << " min Hessian response: " << minmax_elements.first
126
+ << " max Hessian response: " << minmax_elements.second << ' \n ' ;
127
+
202
128
auto corner_points = suppress (gil::view (hessian_response), hessian_determinant_threshold);
129
+ // std::vector<gil::point_t> corner_points = {};
203
130
for (auto point: corner_points) {
204
131
input_view (point) = gil::rgb8_pixel_t (0 , 0 , 0 );
205
132
input_view (point).at (std::integral_constant<int , 1 >{}) = 255 ;
0 commit comments