From e1dba2c6d28c21cf0661005cd961a670e46375c2 Mon Sep 17 00:00:00 2001
From: Robert Lexmann <rlexmann@gmail.com>
Date: Mon, 10 Jun 2024 15:47:29 +0200
Subject: [PATCH] Perform arithmetics in CV_32F branch of divSpectrums() with
 doubles to prevent infs/NaNs (+ corresponding test).

---
 modules/imgproc/src/phasecorr.cpp | 12 ++++++------
 modules/imgproc/test/test_pc.cpp  | 29 +++++++++++++++++++++++++++++
 2 files changed, 35 insertions(+), 6 deletions(-)

diff --git a/modules/imgproc/src/phasecorr.cpp b/modules/imgproc/src/phasecorr.cpp
index dadc3a3da7..d2f88420be 100644
--- a/modules/imgproc/src/phasecorr.cpp
+++ b/modules/imgproc/src/phasecorr.cpp
@@ -252,18 +252,18 @@ void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int fla
             if( !conjB )
                 for( j = j0; j < j1; j += 2 )
                 {
-                    double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps);
-                    double re = (double)(dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]);
-                    double im = (double)(dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]);
+                    double denom = (double)dataB[j]*dataB[j] + (double)dataB[j+1]*dataB[j+1] + (double)eps;
+                    double re = (double)dataA[j]*dataB[j] + (double)dataA[j+1]*dataB[j+1];
+                    double im = (double)dataA[j+1]*dataB[j] - (double)dataA[j]*dataB[j+1];
                     dataC[j] = (float)(re / denom);
                     dataC[j+1] = (float)(im / denom);
                 }
             else
                 for( j = j0; j < j1; j += 2 )
                 {
-                    double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps);
-                    double re = (double)(dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]);
-                    double im = (double)(dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]);
+                    double denom = (double)dataB[j]*dataB[j] + (double)dataB[j+1]*dataB[j+1] + (double)eps;
+                    double re = (double)dataA[j]*dataB[j] - (double)dataA[j+1]*dataB[j+1];
+                    double im = (double)dataA[j+1]*dataB[j] + (double)dataA[j]*dataB[j+1];
                     dataC[j] = (float)(re / denom);
                     dataC[j+1] = (float)(im / denom);
                 }
diff --git a/modules/imgproc/test/test_pc.cpp b/modules/imgproc/test/test_pc.cpp
index 87f0d804b2..969f5bcfa1 100644
--- a/modules/imgproc/test/test_pc.cpp
+++ b/modules/imgproc/test/test_pc.cpp
@@ -120,6 +120,35 @@ TEST(Imgproc_PhaseCorrelatorTest, accuracy_1d_odd_fft) {
     ASSERT_NEAR(phaseShift.x, (double)xShift, 1.);
 }
 
+TEST(Imgproc_PhaseCorrelatorTest, float32_overflow) {
+    // load
+    Mat im = imread(cvtest::TS::ptr()->get_data_path() + "shared/baboon.png", IMREAD_GRAYSCALE);
+    ASSERT_EQ(im.type(), CV_8UC1);
+
+    // convert to 32F, scale values as if original image was 16U
+    constexpr auto u8Max = std::numeric_limits<std::uint8_t>::max();
+    constexpr auto u16Max = std::numeric_limits<std::uint16_t>::max();
+    im.convertTo(im, CV_32FC1, double(u16Max) / double(u8Max));
+
+    // enlarge and create ROIs
+    const auto w = im.cols * 5;
+    const auto h = im.rows * 5;
+    const auto roiW = (w * 2) / 3; // 50% overlap
+    Mat imLarge;
+    resize(im, imLarge, { w, h });
+    const auto roiLeft = imLarge(Rect(0, 0, roiW, h));
+    const auto roiRight = imLarge(Rect(w - roiW, 0, roiW, h));
+
+    // correlate
+    double response = 0.0;
+    Point2d phaseShift = phaseCorrelate(roiLeft, roiRight, cv::noArray(), &response);
+    ASSERT_TRUE(std::isnormal(phaseShift.x) || 0.0 == phaseShift.x);
+    ASSERT_TRUE(std::isnormal(phaseShift.y) || 0.0 == phaseShift.y);
+    ASSERT_TRUE(std::isnormal(response) || 0.0 == response);
+    EXPECT_NEAR(std::abs(phaseShift.x), w / 3.0, 1.0);
+    EXPECT_NEAR(std::abs(phaseShift.y), 0.0, 1.0);
+}
+
 ////////////////////// DivSpectrums ////////////////////////
 class CV_DivSpectrumsTest : public cvtest::ArrayTest
 {