libs/gil/example/hessian.cpp

   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>
   6 #include <vector>
   7 #include <functional>
   8 #include <set>
   9 #include <iostream>
  10 #include <fstream>
  11
  12 namespace gil = boost::gil;
  13
  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 auto filter_height = 5ull;
  64     constexpr static auto 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 std::vector<gil::point_t> suppress(
 107     gil::gray32f_view_t harris_response,
 108     double harris_response_threshold)
 109 {
 110     std::vector<gil::point_t> corner_points;
 111     for (gil::gray32f_view_t::coord_t y = 1; y < harris_response.height() - 1; ++y)
 112     {
 113         for (gil::gray32f_view_t::coord_t x = 1; x < harris_response.width() - 1; ++x)
 114         {
 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))
 128             };
 129
 130             auto maxima = *std::max_element(
 131                 values,
 132                 values + 9,
 133                 [](double lhs, double rhs)
 134                 {
 135                     return lhs < rhs;
 136                 }
 137             );
 138
 139             if (maxima == value(harris_response(x, y))
 140                 && std::count(values, values + 9, maxima) == 1
 141                 && maxima >= harris_response_threshold)
 142             {
 143                 corner_points.emplace_back(x, y);
 144             }
 145         }
 146     }
 147
 148     return corner_points;
 149 }
 150
 151 int main(int argc, char* argv[]) {
 152     if (argc != 5)
 153     {
 154         std::cout << "usage: " << argv[0] << " <input.png> <odd-window-size>"
 155             " <hessian-response-threshold> <output.png>\n";
 156         return -1;
 157     }
 158
 159     std::size_t window_size = std::stoul(argv[2]);
 160     long hessian_determinant_threshold = std::stol(argv[3]);
 161
 162     gil::rgb8_image_t input_image;
 163
 164     gil::read_image(argv[1], input_image, gil::png_tag{});
 165
 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());
 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());
 184     gil::compute_hessian_entries(
 185         x_gradient,
 186         y_gradient,
 187         gil::view(m11),
 188         gil::view(m12_21),
 189         gil::view(m22)
 190     );
 191
 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(
 195         gil::view(m11),
 196         gil::view(m12_21),
 197         gil::view(m22),
 198         gaussian_kernel,
 199         gil::view(hessian_response)
 200     );
 201
 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;
 206     }
 207     gil::write_view(argv[4], input_view, gil::png_tag{});
 208 }