diff --git a/modules/imgproc/src/templmatch.cpp b/modules/imgproc/src/templmatch.cpp
index 827af3eb47..ac01c85d74 100644
--- a/modules/imgproc/src/templmatch.cpp
+++ b/modules/imgproc/src/templmatch.cpp
@@ -975,13 +975,14 @@ static void common_matchTemplate( Mat& img, Mat& templ, Mat& result, int method,
 
     for( i = 0; i < result.rows; i++ )
     {
-        float* rrow = result.ptr<float>(i);
+        float*  rrow  = result.depth() == CV_32F ? result.ptr<float>(i)  : nullptr;
+        double* drow  = result.depth() == CV_64F ? result.ptr<double>(i) : nullptr;
         int idx = i * sumstep;
         int idx2 = i * sqstep;
 
         for( j = 0; j < result.cols; j++, idx += cn, idx2 += cn )
         {
-            double num = rrow[j], t;
+            double num = rrow ? (double)rrow[j] : drow[j], t;
             double wndMean2 = 0, wndSum2 = 0;
 
             if( numType == 1 )
@@ -1027,7 +1028,8 @@ static void common_matchTemplate( Mat& img, Mat& templ, Mat& result, int method,
                     num = method != cv::TM_SQDIFF_NORMED ? 0 : 1;
             }
 
-            rrow[j] = (float)num;
+            if (rrow) rrow[j] = (float)num;
+            else       drow[j] = num;
         }
     }
 }
@@ -1120,6 +1122,11 @@ static bool ipp_matchTemplate( Mat& img, Mat& templ, Mat& result, int method)
     if(templ.size().area()*4 > img.size().area())
         return false;
 
+    // CV_8U SQDIFF/SQDIFF_NORMED suffer from float32 catastrophic cancellation
+    // in IPP's internal accumulators; fall through to the double-precision path instead.
+    if(img.depth() == CV_8U && (method == cv::TM_SQDIFF || method == cv::TM_SQDIFF_NORMED))
+        return false;
+
     if(method == cv::TM_SQDIFF)
     {
         if(ipp_sqrDistance(img, templ, result))
@@ -1192,9 +1199,18 @@ void cv::matchTemplate( InputArray _img, InputArray _templ, OutputArray _result,
 
     CV_IPP_RUN_FAST(ipp_matchTemplate(img, templ, result, method))
 
-    crossCorr( img, templ, result, Point(0,0), 0, 0);
+    bool use64f = (depth == CV_8U) && (method == cv::TM_SQDIFF || method == cv::TM_SQDIFF_NORMED);
+    Mat result64f;
+    Mat& workResult = use64f ? result64f : result;
+    if (use64f)
+        result64f.create(corrSize, CV_64F);
 
-    common_matchTemplate(img, templ, result, method, cn);
+    crossCorr( img, templ, workResult, Point(0,0), 0, 0);
+
+    common_matchTemplate(img, templ, workResult, method, cn);
+
+    if (use64f)
+        result64f.convertTo(result, CV_32F);
 }
 
 CV_IMPL void
diff --git a/modules/imgproc/test/test_templmatch.cpp b/modules/imgproc/test/test_templmatch.cpp
index a10f44ce0f..76229100fa 100644
--- a/modules/imgproc/test/test_templmatch.cpp
+++ b/modules/imgproc/test/test_templmatch.cpp
@@ -342,4 +342,19 @@ INSTANTIATE_TEST_CASE_P(/**/,
             testing::Values(TM_SQDIFF, TM_SQDIFF_NORMED, TM_CCORR, TM_CCORR_NORMED, TM_CCOEFF, TM_CCOEFF_NORMED)));
 
 
+TEST(Imgproc_MatchTemplate, bug_21786)
+{
+    // CV_8U identical image/template with large patch sums triggers float32
+    // catastrophic cancellation in TM_SQDIFF. Result must be exactly zero.
+    Mat img(100, 100, CV_8U, Scalar(255));
+    Mat templ(25, 25, CV_8U, Scalar(255));
+    Mat result;
+
+    matchTemplate(img, templ, result, TM_SQDIFF);
+    EXPECT_NEAR(0.0, cv::norm(result, cv::NORM_INF), 1e-6);
+
+    matchTemplate(img, templ, result, TM_SQDIFF_NORMED);
+    EXPECT_NEAR(0.0, cv::norm(result, cv::NORM_INF), 1e-6);
+}
+
 }} // namespace