1 //
2 // Copyright 2019 Olzhas Zhumabek <anonymous.from.applecity@gmail.com>
3 //
4 // Use, modification and distribution are subject to the Boost Software License,
5 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8 #include <boost/gil/image.hpp>
9 #include <boost/gil/image_view.hpp>
10 #include <boost/gil/extension/io/png.hpp>
11 #include <boost/gil/image_processing/numeric.hpp>
12 #include <boost/gil/image_processing/harris.hpp>
13 #include <boost/gil/extension/numeric/convolve.hpp>
14 #include <vector>
15 #include <functional>
16 #include <set>
17 #include <iostream>
18 #include <fstream>
19 
20 namespace gil = boost::gil;
21 
22 // some images might produce artifacts
23 // when converted to grayscale,
24 // which was previously observed on
25 // canny edge detector for test input
26 // used for this example
to_grayscale(gil::rgb8_view_t original)27 gil::gray8_image_t to_grayscale(gil::rgb8_view_t original)
28 {
29     gil::gray8_image_t output_image(original.dimensions());
30     auto output = gil::view(output_image);
31     constexpr double max_channel_intensity = (std::numeric_limits<std::uint8_t>::max)();
32     for (long int y = 0; y < original.height(); ++y) {
33         for (long int x = 0; x < original.width(); ++x) {
34             // scale the values into range [0, 1] and calculate linear intensity
35             double red_intensity = original(x, y).at(std::integral_constant<int, 0>{})
36                 / max_channel_intensity;
37             double green_intensity = original(x, y).at(std::integral_constant<int, 1>{})
38                 / max_channel_intensity;
39             double blue_intensity = original(x, y).at(std::integral_constant<int, 2>{})
40                 / max_channel_intensity;
41             auto linear_luminosity = 0.2126 * red_intensity
42                                     + 0.7152 * green_intensity
43                                     + 0.0722 * blue_intensity;
44 
45             // perform gamma adjustment
46             double gamma_compressed_luminosity = 0;
47             if (linear_luminosity < 0.0031308) {
48                 gamma_compressed_luminosity = linear_luminosity * 12.92;
49             } else {
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 
apply_gaussian_blur(gil::gray8_view_t input_view,gil::gray8_view_t output_view)61 void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_view)
62 {
63     constexpr static auto filterHeight = 5ull;
64     constexpr static auto filterWidth = 5ull;
65     constexpr static double filter[filterHeight][filterWidth] =
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 (long x = 0; x < width; ++x) {
79         for (long y = 0; y < height; ++y) {
80             double intensity = 0.0;
81             for (size_t filter_y = 0; filter_y < filterHeight; ++filter_y) {
82                 for (size_t filter_x = 0; filter_x < filterWidth; ++filter_x) {
83                     int image_x = x - filterWidth / 2 + filter_x;
84                     int image_y = y - filterHeight / 2 + filter_y;
85                     if (image_x >= input_view.width() || image_x < 0
86                         || image_y >= input_view.height() || image_y < 0) {
87                         continue;
88                     }
89                     auto& pixel = input_view(image_x, image_y);
90                     intensity += pixel.at(std::integral_constant<int, 0>{})
91                         * filter[filter_y][filter_x];
92                 }
93             }
94             auto& pixel = output_view(gil::point_t(x, y));
95             pixel = (std::min)((std::max)(int(factor * intensity + bias), 0), 255);
96         }
97 
98     }
99 }
100 
suppress(gil::gray32f_view_t harris_response,double harris_response_threshold)101 std::vector<gil::point_t> suppress(
102     gil::gray32f_view_t harris_response,
103     double harris_response_threshold)
104 {
105     std::vector<gil::point_t> corner_points;
106     for (gil::gray32f_view_t::coord_t y = 1; y < harris_response.height() - 1; ++y)
107     {
108         for (gil::gray32f_view_t::coord_t x = 1; x < harris_response.width() - 1; ++x)
109         {
110             auto value = [](gil::gray32f_pixel_t pixel) {
111                 return pixel.at(std::integral_constant<int, 0>{});
112             };
113             double values[9] = {
114                 value(harris_response(x - 1, y - 1)),
115                 value(harris_response(x, y - 1)),
116                 value(harris_response(x + 1, y - 1)),
117                 value(harris_response(x - 1, y)),
118                 value(harris_response(x, y)),
119                 value(harris_response(x + 1, y)),
120                 value(harris_response(x - 1, y + 1)),
121                 value(harris_response(x, y + 1)),
122                 value(harris_response(x + 1, y + 1))
123             };
124 
125             auto maxima = *std::max_element(
126                 values,
127                 values + 9,
128                 [](double lhs, double rhs)
129                 {
130                     return lhs < rhs;
131                 }
132             );
133 
134             if (maxima == value(harris_response(x, y))
135                 && std::count(values, values + 9, maxima) == 1
136                 && maxima >= harris_response_threshold)
137             {
138                 corner_points.emplace_back(x, y);
139             }
140         }
141     }
142 
143     return corner_points;
144 }
145 
main(int argc,char * argv[])146 int main(int argc, char* argv[])
147 {
148     if (argc != 6)
149     {
150         std::cout << "usage: " << argv[0] << " <input.png> <odd-window-size>"
151             " <discrimination-constant> <harris-response-threshold> <output.png>\n";
152         return -1;
153     }
154 
155     std::size_t window_size = std::stoul(argv[2]);
156     double discrimnation_constant = std::stof(argv[3]);
157     long harris_response_threshold = std::stol(argv[4]);
158 
159     gil::rgb8_image_t input_image;
160 
161     gil::read_image(argv[1], input_image, gil::png_tag{});
162 
163     auto input_view = gil::view(input_image);
164     auto grayscaled = to_grayscale(input_view);
165     gil::gray8_image_t smoothed_image(grayscaled.dimensions());
166     auto smoothed = gil::view(smoothed_image);
167     apply_gaussian_blur(gil::view(grayscaled), smoothed);
168     gil::gray16s_image_t x_gradient_image(grayscaled.dimensions());
169     gil::gray16s_image_t y_gradient_image(grayscaled.dimensions());
170 
171     auto x_gradient = gil::view(x_gradient_image);
172     auto y_gradient = gil::view(y_gradient_image);
173     auto scharr_x = gil::generate_dx_scharr();
174     gil::detail::convolve_2d(smoothed, scharr_x, x_gradient);
175     auto scharr_y = gil::generate_dy_scharr();
176     gil::detail::convolve_2d(smoothed, scharr_y, y_gradient);
177 
178     gil::gray32f_image_t m11(x_gradient.dimensions());
179     gil::gray32f_image_t m12_21(x_gradient.dimensions());
180     gil::gray32f_image_t m22(x_gradient.dimensions());
181     gil::compute_tensor_entries(
182         x_gradient,
183         y_gradient,
184         gil::view(m11),
185         gil::view(m12_21),
186         gil::view(m22)
187     );
188 
189     gil::gray32f_image_t harris_response(x_gradient.dimensions());
190     auto gaussian_kernel = gil::generate_gaussian_kernel(window_size, 0.84089642);
191     gil::compute_harris_responses(
192         gil::view(m11),
193         gil::view(m12_21),
194         gil::view(m22),
195         gaussian_kernel,
196         discrimnation_constant,
197         gil::view(harris_response)
198     );
199 
200     auto corner_points = suppress(gil::view(harris_response), harris_response_threshold);
201     for (auto point: corner_points)
202     {
203         input_view(point) = gil::rgb8_pixel_t(0, 0, 0);
204         input_view(point).at(std::integral_constant<int, 1>{}) = 255;
205     }
206     gil::write_view(argv[5], input_view, gil::png_tag{});
207 }
208