Fix pixel value evaluation overflow in bit-exact GaussianBlur implementation

modules/imgproc/src/smooth.simd.hpp

@@ -334,7 +334,7 @@ void hlineSmooth3Naba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const
         {
             int src_idx = borderInterpolate(-1, len, borderType);
             for (int k = 0; k < cn; k++)
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[1] * src[k] + ((uint16_t*)m)[0] * ((uint16_t)(src[cn + k]) + (uint16_t)(src[src_idx*cn + k]));
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * (uint32_t)(src[k]) + ((uint16_t*)m)[0] * ((uint32_t)(src[cn + k]) + (uint32_t)(src[src_idx*cn + k])));
         }
         else
         {
@@ -354,14 +354,14 @@ void hlineSmooth3Naba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const
                                     v_mul_wrap(vx_load_expand(src), v_mul1));
 #endif
         for (; i < lencn; i++, src++, dst++)
-            *((uint16_t*)dst) = ((uint16_t*)m)[1] * src[0] + ((uint16_t*)m)[0] * ((uint16_t)(src[-cn]) + (uint16_t)(src[cn]));
+            *((uint16_t*)dst) = saturate_cast<uint16_t>(((uint16_t*)m)[1] * (uint32_t)(src[0]) + ((uint16_t*)m)[0] * ((uint32_t)(src[-cn]) + (uint32_t)(src[cn])));
 
         // Point that fall right from border
         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
         {
             int src_idx = (borderInterpolate(len, len, borderType) - (len - 1))*cn;
             for (int k = 0; k < cn; k++)
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[1] * src[k] + ((uint16_t*)m)[0] * ((uint16_t)(src[k - cn]) + (uint16_t)(src[src_idx + k]));
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * (uint32_t)(src[k]) + ((uint16_t*)m)[0] * ((uint32_t)(src[k - cn]) + (uint32_t)(src[src_idx + k])));
         }
         else
         {
@@ -896,8 +896,8 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             int idxp2 = borderInterpolate(3, len, borderType)*cn;
             for (int k = 0; k < cn; k++)
             {
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[1] * ((uint16_t)(src[k + idxm1]) + (uint16_t)(src[k + cn])) + ((uint16_t*)m)[2] * src[k] + ((uint16_t*)m)[0] * ((uint16_t)(src[k + idxp1]) + (uint16_t)(src[k + idxm2]));
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[0] * ((uint16_t)(src[k + idxm1]) + (uint16_t)(src[k + idxp2])) + ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[k + idxp1])) + ((uint16_t*)m)[2] * src[k + cn];
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * ((uint32_t)(src[k + idxm1]) + (uint32_t)(src[k + cn])) + ((uint16_t*)m)[2] * (uint32_t)(src[k]) + ((uint16_t*)m)[0] * ((uint32_t)(src[k + idxp1]) + (uint32_t)(src[k + idxm2])));
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[0] * ((uint32_t)(src[k + idxm1]) + (uint32_t)(src[k + idxp2])) + ((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[k + idxp1])) + ((uint16_t*)m)[2] * (uint32_t)(src[k + cn]));
             }
         }
     }
@@ -907,7 +907,7 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             for (int k = 0; k < cn; k++)
             {
                 dst[k] = m[2] * src[k] + m[1] * src[k + cn] + m[0] * src[k + 2 * cn];
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[k + 2 * cn])) + ((uint16_t*)m)[2] * src[k + cn];
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[k + 2 * cn])) + ((uint16_t*)m)[2] * (uint32_t)(src[k + cn]));
                 dst[k + 2 * cn] = m[0] * src[k] + m[1] * src[k + cn] + m[2] * src[k + 2 * cn];
             }
         else
@@ -918,9 +918,9 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             int idxp2 = borderInterpolate(4, len, borderType)*cn;
             for (int k = 0; k < cn; k++)
             {
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[2] * src[k] + ((uint16_t*)m)[1] * ((uint16_t)(src[k + cn]) + (uint16_t)(src[k + idxm1])) + ((uint16_t*)m)[0] * ((uint16_t)(src[k + 2 * cn]) + (uint16_t)(src[k + idxm2]));
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[2] * src[k + cn] + ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[k + 2 * cn])) + ((uint16_t*)m)[0] * ((uint16_t)(src[k + idxm1]) + (uint16_t)(src[k + idxp1]));
-                ((uint16_t*)dst)[k + 2 * cn] = ((uint16_t*)m)[0] * ((uint16_t)(src[k]) + (uint16_t)(src[k + idxp2])) + ((uint16_t*)m)[1] * ((uint16_t)(src[k + cn]) + (uint16_t)(src[k + idxp1])) + ((uint16_t*)m)[2] * src[k + 2 * cn];
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[2] * (uint32_t)(src[k]) + ((uint16_t*)m)[1] * ((uint32_t)(src[k + cn]) + (uint32_t)(src[k + idxm1])) + ((uint16_t*)m)[0] * ((uint32_t)(src[k + 2 * cn]) + (uint32_t)(src[k + idxm2])));
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[2] * (uint32_t)(src[k + cn]) + ((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[k + 2 * cn])) + ((uint16_t*)m)[0] * ((uint32_t)(src[k + idxm1]) + (uint32_t)(src[k + idxp1])));
+                ((uint16_t*)dst)[k + 2 * cn] = saturate_cast<uint16_t>(((uint16_t*)m)[0] * ((uint32_t)(src[k]) + (uint32_t)(src[k + idxp2])) + ((uint16_t*)m)[1] * ((uint32_t)(src[k + cn]) + (uint32_t)(src[k + idxp1])) + ((uint16_t*)m)[2] * (uint32_t)(src[k + 2 * cn]));
             }
         }
     }
@@ -933,8 +933,8 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
             for (int k = 0; k < cn; k++)
             {
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[2] * src[k] + ((uint16_t*)m)[1] * ((uint16_t)(src[cn + k]) + (uint16_t)(src[idxm1 + k])) + ((uint16_t*)m)[0] * ((uint16_t)(src[2 * cn + k]) + (uint16_t)(src[idxm2 + k]));
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[2 * cn + k])) + ((uint16_t*)m)[2] * src[cn + k] + ((uint16_t*)m)[0] * ((uint16_t)(src[3 * cn + k]) + (uint16_t)(src[idxm1 + k]));
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[2] * (uint32_t)(src[k]) + ((uint16_t*)m)[1] * ((uint32_t)(src[cn + k]) + (uint32_t)(src[idxm1 + k])) + ((uint16_t*)m)[0] * ((uint32_t)(src[2 * cn + k]) + (uint32_t)(src[idxm2 + k])));
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[2 * cn + k])) + ((uint16_t*)m)[2] * (uint32_t)(src[cn + k]) + ((uint16_t*)m)[0] * ((uint32_t)(src[3 * cn + k]) + (uint32_t)(src[idxm1 + k])));
             }
         }
         else
@@ -942,7 +942,7 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             for (int k = 0; k < cn; k++)
             {
                 dst[k] = m[2] * src[k] + m[1] * src[cn + k] + m[0] * src[2 * cn + k];
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[2 * cn + k])) + ((uint16_t*)m)[2] * src[cn + k] + ((uint16_t*)m)[0] * src[3 * cn + k];
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[2 * cn + k])) + ((uint16_t*)m)[2] * (uint32_t)(src[cn + k]) + ((uint16_t*)m)[0] * (uint32_t)(src[3 * cn + k]));
             }
         }
 
@@ -960,7 +960,7 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
                                     v_mul_wrap(vx_load_expand(src), v_mul2));
 #endif
         for (; i < lencn; i++, src++, dst++)
-            *((uint16_t*)dst) = ((uint16_t*)m)[0] * ((uint16_t)(src[-2 * cn]) + (uint16_t)(src[2 * cn])) + ((uint16_t*)m)[1] * ((uint16_t)(src[-cn]) + (uint16_t)(src[cn])) + ((uint16_t*)m)[2] * src[0];
+            *((uint16_t*)dst) = saturate_cast<uint16_t>(((uint16_t*)m)[0] * ((uint32_t)(src[-2 * cn]) + (uint32_t)(src[2 * cn])) + ((uint16_t*)m)[1] * ((uint32_t)(src[-cn]) + (uint32_t)(src[cn])) + ((uint16_t*)m)[2] * (uint32_t)(src[0]));
 
         // Points that fall right from border
         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
@@ -969,15 +969,15 @@ void hlineSmooth5Nabcba<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, cons
             int idxp2 = (borderInterpolate(len + 1, len, borderType) - (len - 2))*cn;
             for (int k = 0; k < cn; k++)
             {
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[0] * ((uint16_t)(src[k - 2 * cn]) + (uint16_t)(src[idxp1 + k])) + ((uint16_t*)m)[1] * ((uint16_t)(src[k - cn]) + (uint16_t)(src[k + cn])) + ((uint16_t*)m)[2] * src[k];
-                ((uint16_t*)dst)[k + cn] = ((uint16_t*)m)[0] * ((uint16_t)(src[k - cn]) + (uint16_t)(src[idxp2 + k])) + ((uint16_t*)m)[1] * ((uint16_t)(src[k]) + (uint16_t)(src[idxp1 + k])) + ((uint16_t*)m)[2] * src[k + cn];
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[0] * ((uint32_t)(src[k - 2 * cn]) + (uint32_t)(src[idxp1 + k])) + ((uint16_t*)m)[1] * ((uint32_t)(src[k - cn]) + (uint32_t)(src[k + cn])) + ((uint16_t*)m)[2] * (uint32_t)(src[k]));
+                ((uint16_t*)dst)[k + cn] = saturate_cast<uint16_t>(((uint16_t*)m)[0] * ((uint32_t)(src[k - cn]) + (uint32_t)(src[idxp2 + k])) + ((uint16_t*)m)[1] * ((uint32_t)(src[k]) + (uint32_t)(src[idxp1 + k])) + ((uint16_t*)m)[2] * (uint32_t)(src[k + cn]));
             }
         }
         else
         {
             for (int k = 0; k < cn; k++)
             {
-                ((uint16_t*)dst)[k] = ((uint16_t*)m)[0] * src[k - 2 * cn] + ((uint16_t*)m)[1] * ((uint16_t)(src[k - cn]) + (uint16_t)(src[k + cn])) + ((uint16_t*)m)[2] * src[k];
+                ((uint16_t*)dst)[k] = saturate_cast<uint16_t>(((uint16_t*)m)[0] * (uint32_t)(src[k - 2 * cn]) + ((uint16_t*)m)[1] * ((uint32_t)(src[k - cn]) + (uint32_t)(src[k + cn])) + ((uint16_t*)m)[2] * (uint32_t)(src[k]));
                 dst[k + cn] = m[0] * src[k - cn] + m[1] * src[k] + m[2] * src[k + cn];
             }
         }

modules/imgproc/test/test_smooth_bitexact.cpp

@@ -158,4 +158,12 @@ TEST(GaussianBlur_Bitexact, Linear8U)
     }
 }
 
+TEST(GaussianBlur_Bitexact, regression_15015)
+{
+    Mat src(100,100,CV_8UC3,Scalar(255,255,255));
+    Mat dst;
+    GaussianBlur(src, dst, Size(5, 5), 9);
+    ASSERT_EQ(0.0, cvtest::norm(dst, src, NORM_INF));
+}
+
 }} // namespace