fix: Normalize Gaussian 2D kernel. (#725)

boostorg · Feb 7, 2023 · d6e67f3 · d6e67f3
1 parent 712b827
commit d6e67f3
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 12 deletions.
diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp
@@ -136,21 +136,26 @@ inline auto generate_gaussian_kernel(std::size_t side_length, double sigma)
         throw std::invalid_argument("kernel dimensions should be odd and equal");
 
     const double denominator = 2 * boost::gil::detail::pi * sigma * sigma;
-    auto middle = side_length / 2;
+    auto const middle = side_length / 2;
     std::vector<T, Allocator> values(side_length * side_length);
+    T sum{0};
     for (std::size_t y = 0; y < side_length; ++y)
     {
         for (std::size_t x = 0; x < side_length; ++x)
         {
-            const auto delta_x = middle > x ? middle - x : x - middle;
-            const auto delta_y = middle > y ? middle - y : y - middle;
-            const double power = (delta_x * delta_x +  delta_y * delta_y) / (2 * sigma * sigma);
+            const auto delta_x = x - middle;
+            const auto delta_y = y - middle;
+            const auto power = static_cast<double>(delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
             const double nominator = std::exp(-power);
-            const float value = static_cast<float>(nominator / denominator);
+            const auto value = static_cast<T>(nominator / denominator);
             values[y * side_length + x] = value;
+            sum += value;
         }
     }
 
+    // normalize so that Gaussian kernel sums up to 1.
+    std::transform(values.begin(), values.end(), values.begin(), [&sum](const auto & v) { return v/sum; });
+
     return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
 }
 

diff --git a/test/core/image_processing/simple_kernels.cpp b/test/core/image_processing/simple_kernels.cpp
@@ -38,15 +38,16 @@ void test_gaussian_kernel_generation()
     auto kernel = boost::gil::generate_gaussian_kernel(7, 0.84089642);
     const float expected_values[7][7] =
     {
-        {0.00000067f, 0.00002292f, 0.00019117f, 0.00038771f, 0.00019117f, 0.00002292f, 0.00000067f},
-        {0.00002292f, 0.00078633f, 0.00655965f, 0.01330373f, 0.00655965f, 0.00078633f, 0.00002292f},
-        {0.00019117f, 0.00655965f, 0.05472157f, 0.11098164f, 0.05472157f, 0.00655965f, 0.00019117f},
-        {0.00038771f, 0.01330373f, 0.11098164f, 0.25508352f, 0.11098164f, 0.01330373f, 0.00038711f},
-        {0.00019117f, 0.00655965f, 0.05472157f, 0.11098164f, 0.05472157f, 0.00655965f, 0.00019117f},
-        {0.00002292f, 0.00078633f, 0.00655965f, 0.01330373f, 0.00655965f, 0.00078633f, 0.00002292f},
-        {0.00000067f, 0.00002292f, 0.00019117f, 0.00038771f, 0.00019117f, 0.00002292f, 0.00000067f}
+        {6.67847706e-07f, 2.29160778e-05f, 1.91169232e-04f, 3.87713181e-04f, 1.91169232e-04f, 2.29160778e-05f, 6.67847706e-07f},
+        {2.29160778e-05f, 7.86326905e-04f, 6.55965268e-03f, 1.33037298e-02f, 6.55965268e-03f, 7.86326905e-04f, 2.29160778e-05f},
+        {1.91169232e-04f, 6.55965268e-03f, 5.47215706e-02f, 1.10981636e-01f, 5.47215706e-02f, 6.55965268e-03f, 1.91169232e-04f},
+        {3.87713181e-04f, 1.33037298e-02f, 1.10981636e-01f, 2.25083518e-01f, 1.10981636e-01f, 1.33037298e-02f, 3.87713181e-04f},
+        {1.91169232e-04f, 6.55965268e-03f, 5.47215706e-02f, 1.10981636e-01f, 5.47215706e-02f, 6.55965268e-03f, 1.91169232e-04f},
+        {2.29160778e-05f, 7.86326905e-04f, 6.55965268e-03f, 1.33037298e-02f, 6.55965268e-03f, 7.86326905e-04f, 2.29160778e-05f},
+        {6.67847706e-07f, 2.29160778e-05f, 1.91169232e-04f, 3.87713181e-04f, 1.91169232e-04f, 2.29160778e-05f, 6.67847706e-07f}
     };
 
+    double sum{0};
     for (gil::gray32f_view_t::coord_t y = 0; static_cast<std::size_t>(y) < kernel.size(); ++y)
     {
         for (gil::gray32f_view_t::coord_t x = 0; static_cast<std::size_t>(x) < kernel.size(); ++x)
@@ -55,8 +56,10 @@ void test_gaussian_kernel_generation()
             auto expected = expected_values[y][x];
             auto percent_difference = std::ceil(std::abs(expected - output) / expected);
             BOOST_TEST_LT(percent_difference, 5);
+            sum += output;
         }
     }
+    BOOST_TEST_LT(std::abs(sum-1), 1e-7);
 }
 
 int main()