Merge pull request #28146 from zdenyhraz:iterative-phase-correlation

Iterative Phase Correlation #28146

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [x] I agree to contribute to the project under Apache 2 License.
- [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [x] The PR is proposed to the proper branch
- [x] There is a reference to the original bug report and related work
- [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [x] The feature is well documented and sample code can be built with the project CMake

doc/opencv.bib

@@ -1561,3 +1561,13 @@
   volume = {7},
   url = {https://doi.org/10.1007/BF01898354}
 }
+@article{hrazdira2020iterative,
+  title={Iterative phase correlation algorithm for high-precision subpixel image registration},
+  author={Hrazd{\'\i}ra, Zdenek and Druckm{\"u}ller, Miloslav and Habbal, Shadia},
+  journal={The Astrophysical Journal Supplement Series},
+  volume={247},
+  number={1},
+  pages={8},
+  year={2020},
+  publisher={IOP Publishing}
+}

modules/imgproc/include/opencv2/imgproc.hpp

@@ -3051,6 +3051,22 @@ peak) and will be smaller when there are multiple peaks.
 CV_EXPORTS_W Point2d phaseCorrelate(InputArray src1, InputArray src2,
                                     InputArray window = noArray(), CV_OUT double* response = 0);
 
+/** @brief Detects translational shifts between two images.
+
+This function extends the standard @ref phaseCorrelate method by improving sub-pixel accuracy
+through iterative shift refinement in the phase-correlation space, as described in
+@cite hrazdira2020iterative.
+
+@param src1 Source floating point array (CV_32FC1 or CV_64FC1)
+@param src2 Source floating point array (CV_32FC1 or CV_64FC1)
+@param L2size The size of the correlation neighborhood used by the iterative shift refinement algorithm.
+@param maxIters The maximum number of iterations the iterative refinement algorithm will run.
+@returns detected sub-pixel shift between the two arrays.
+
+@sa phaseCorrelate, dft, idft, createHanningWindow
+ */
+CV_EXPORTS_W Point2d phaseCorrelateIterative(InputArray src1, InputArray src2, int L2size = 7, int maxIters = 10);
+
 /** @brief This function computes a Hanning window coefficients in two dimensions.
 
 See (https://en.wikipedia.org/wiki/Hann_function) and (https://en.wikipedia.org/wiki/Window_function)

modules/imgproc/src/phasecorr_iterative.cpp

@@ -0,0 +1,191 @@
+#include "precomp.hpp"
+#include <cmath>
+
+namespace {
+
+template <typename T>
+void calculateCrossPowerSpectrum(const cv::Mat& dft1, const cv::Mat& dft2, cv::Mat& cps)
+{
+    for (int row = 0; row < dft1.rows; ++row)
+    {
+        auto* cpsp = cps.ptr<cv::Vec<T, 2>>(row);
+        const auto* dft1p = dft1.ptr<cv::Vec<T, 2>>(row);
+        const auto* dft2p = dft2.ptr<cv::Vec<T, 2>>(row);
+        for (int col = 0; col < dft1.cols; ++col)
+        {
+            const T re = dft1p[col][0] * dft2p[col][0] + dft1p[col][1] * dft2p[col][1];
+            const T im = dft1p[col][0] * dft2p[col][1] - dft1p[col][1] * dft2p[col][0];
+            const T mag = std::sqrt(re * re + im * im);
+            cpsp[col][0] = re / mag;
+            cpsp[col][1] = im / mag;
+        }
+    }
+}
+
+cv::Mat calculateCrossPowerSpectrum(const cv::Mat& dft1, const cv::Mat& dft2)
+{
+    cv::Mat cps(dft1.rows, dft1.cols, dft1.type());
+    if (dft1.type() == CV_32FC2)
+        calculateCrossPowerSpectrum<float>(dft1, dft2, cps);
+    else if (dft1.type() == CV_64FC2)
+        calculateCrossPowerSpectrum<double>(dft1, dft2, cps);
+    else
+        CV_Error(cv::Error::StsNotImplemented, "Only CV_32FC2 and CV_64FC2 types are supported");
+    return cps;
+}
+
+void fftshift(cv::Mat& out)
+{
+    int cx = out.cols / 2;
+    int cy = out.rows / 2;
+    cv::Mat q0(out, cv::Rect(0, 0, cx, cy));
+    cv::Mat q1(out, cv::Rect(cx, 0, cx, cy));
+    cv::Mat q2(out, cv::Rect(0, cy, cx, cy));
+    cv::Mat q3(out, cv::Rect(cx, cy, cx, cy));
+
+    cv::Mat tmp;
+    q0.copyTo(tmp);
+    q3.copyTo(q0);
+    tmp.copyTo(q3);
+    q1.copyTo(tmp);
+    q2.copyTo(q1);
+    tmp.copyTo(q2);
+}
+
+bool isOutOfBounds(const cv::Point2i& peak, const cv::Mat& mat, int size)
+{
+    return peak.x - size / 2 < 0 || peak.y - size / 2 < 0 || peak.x + size / 2 >= mat.cols ||
+           peak.y + size / 2 >= mat.rows;
+}
+
+bool reduceL2size(int& L2size)
+{
+    L2size -= 2;
+    return L2size >= 3;
+}
+
+int getL1size(int L2Usize, double L1ratio)
+{
+    int L1size = static_cast<int>(std::floor(L1ratio * L2Usize));
+    return (L1size % 2) ? L1size : L1size + 1;
+}
+
+cv::Point2d getPeakSubpixel(const cv::Mat& mat)
+{
+    const auto m = moments(mat);
+    return cv::Point2d(m.m10 / m.m00, m.m01 / m.m00);
+}
+
+bool accuracyReached(const cv::Point2d& L1peak, const cv::Point2d& L1mid)
+{
+    return std::abs(L1peak.x - L1mid.x) < 0.5 && std::abs(L1peak.y - L1mid.y) < 0.5;
+}
+
+cv::Point2d getSubpixelShift(const cv::Mat& L3,
+                             const cv::Point2d& L3peak,
+                             const cv::Point2d& L3mid,
+                             int L2size)
+{
+    while (isOutOfBounds(L3peak, L3, L2size))
+        if (!reduceL2size(L2size))
+            return L3peak - L3mid;
+
+    cv::Mat L2 = L3(cv::Rect(static_cast<int>(L3peak.x - L2size / 2),
+                             static_cast<int>(L3peak.y - L2size / 2),
+                             L2size,
+                             L2size));
+    cv::Point2d L2peak = getPeakSubpixel(L2);
+    cv::Point2d L2mid(L2.cols / 2, L2.rows / 2);
+    return L3peak - L3mid + L2peak - L2mid;
+}
+
+}  // namespace
+
+cv::Point2d
+    cv::phaseCorrelateIterative(InputArray _src1, InputArray _src2, int L2size, int maxIters)
+{
+    CV_INSTRUMENT_REGION();
+
+    Mat src1 = _src1.getMat();
+    Mat src2 = _src2.getMat();
+
+    CV_Assert(src1.type() == src2.type());
+    CV_Assert(src1.size() == src2.size());
+    CV_Assert(src1.type() == CV_32FC1 || src1.type() == CV_64FC1);
+
+    // apply DFT window to input images
+    Mat window;
+    createHanningWindow(window, src1.size(), _src1.type());
+    Mat image1, image2;
+    multiply(_src1, window, image1);
+    multiply(_src2, window, image2);
+
+    // compute the DFTs of input images
+    dft(image1, image1, DFT_COMPLEX_OUTPUT);
+    dft(image2, image2, DFT_COMPLEX_OUTPUT);
+
+    // compute the phase correlation landscape L3
+    Mat L3 = calculateCrossPowerSpectrum(image1, image2);
+    dft(L3, L3, DFT_INVERSE | DFT_SCALE | DFT_REAL_OUTPUT);
+    fftshift(L3);
+    Point2d L3mid(L3.cols / 2, L3.rows / 2);
+
+    // calculate the maximum correlation location
+    Point2i L3peak;
+    minMaxLoc(L3, nullptr, nullptr, nullptr, &L3peak);
+
+    // reduce the L2size as long as L2 is out of bounds of L3
+    while (isOutOfBounds(L3peak, L3, L2size))
+        if (!reduceL2size(L2size))
+            return Point2d(L3peak) - L3mid;
+
+    // extract the L2 maximum correlation neighborhood from L3
+    Mat L2 = L3(Rect(L3peak.x - L2size / 2, L3peak.y - L2size / 2, L2size, L2size));
+
+    // upsample L2 maximum correlation neighborhood to get L2U
+    Mat L2U;
+    const int L2Usize = 223;  // empirically determined optimal constant
+    resize(L2, L2U, {L2Usize, L2Usize}, 0, 0, INTER_LINEAR);
+    const Point2d L2Umid(L2U.cols / 2, L2U.rows / 2);
+
+    // run the iterative refinement algorithm using the specified L1 ratio
+    // gradually decrease L1 ratio if convergence is not achieved
+    const double L1ratioBase = 0.45;  // empirically determined optimal constant
+    const double L1ratioStep = 0.025;
+    for (double L1ratio = L1ratioBase; getL1size(L2U.cols, L1ratio) > 0; L1ratio -= L1ratioStep)
+    {
+        Point2d L2Upeak = L2Umid;  // reset the accumulated L2U peak position
+        const int L1size = getL1size(L2U.cols, L1ratio);  // calculate the current L1 size
+        const Point2d L1mid(L1size / 2, L1size / 2);  // update the L1 mid position
+
+        // perform the iterative refinement algorithm
+        for (int iter = 0; iter < maxIters; ++iter)
+        {
+            // verify that the L1 region is within the L2U region
+            if (isOutOfBounds(L2Upeak, L2U, L1size))
+                break;
+
+            // extract the L1 region from L2U
+            const Mat L1 = L2U(Rect(static_cast<int>(L2Upeak.x - L1size / 2),
+                                    static_cast<int>(L2Upeak.y - L1size / 2),
+                                    L1size,
+                                    L1size));
+
+            // calculate the centroid location
+            const Point2d L1peak = getPeakSubpixel(L1);
+
+            // add the contribution of the current iteration to the accumulated L2U peak location
+            L2Upeak += Point2d(std::round(L1peak.x - L1mid.x), std::round(L1peak.y - L1mid.y));
+
+            // check for convergence
+            if (accuracyReached(L1peak, L1mid))
+                // return the refined subpixel image shift
+                return Point2d(L3peak) - L3mid +
+                       (L2Upeak - L2Umid + L1peak - L1mid) /
+                           (static_cast<double>(L2Usize) / L2size);
+        }
+    }
+
+    // iterative refinement failed to converge, return non-iterative subpixel shift
+    return getSubpixelShift(L3, Point2d(L3peak), L3mid, L2size);
+}

modules/imgproc/test/test_ipc.cpp

@@ -0,0 +1,117 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "test_precomp.hpp"
+#include <vector>
+
+namespace opencv_test { namespace {
+
+Mat CropMid(InputArray src, int w, int h)
+{
+    Mat mat = src.getMat();
+    return mat(Rect(mat.cols / 2 - w / 2, mat.rows / 2 - h / 2, w, h));
+}
+
+Mat GenerateTestImage(Size size)
+{
+    Mat image = Mat::zeros(size.height * 2, size.width * 2, CV_32F);
+    rectangle(image,
+              Point(static_cast<int>(size.width * 0.1), static_cast<int>(size.height * 0.1)),
+              Point(static_cast<int>(size.width * 0.9), static_cast<int>(size.height * 0.9)),
+              Scalar(1),
+              -1);
+    return image;
+}
+
+void TestPhaseCorrelationIterative(const Size& size, const double maxShift)
+{
+    const auto iters = std::max(201., maxShift * 10 + 1);
+    const Point2d shiftOffset(-maxShift * 0.5, -maxShift * 0.5);
+    Mat image1 = GenerateTestImage(size);
+    Mat crop1 = CropMid(image1, size.width, size.height);
+    Mat image2 = image1.clone();
+
+    std::vector<double> pcErrors;
+    std::vector<double> ipcErrors;
+
+    for (int i = 0; i < iters; ++i)
+    {
+        const auto shift =
+            Point2d(maxShift * i / (iters - 1), maxShift * i / (iters - 1)) + shiftOffset;
+        const Mat Tmat = (Mat_<double>(2, 3) << 1., 0., shift.x, 0., 1., shift.y);
+        warpAffine(image1, image2, Tmat, image2.size());
+        Mat crop2 = CropMid(image2, size.width, size.height);
+        const auto ipcshift = phaseCorrelateIterative(crop1, crop2);
+        const auto pcshift = phaseCorrelate(crop1, crop2);
+
+        pcErrors.push_back(
+            0.5 * std::abs(pcshift.x - shift.y) + 0.5 * std::abs(pcshift.y - shift.x));
+        ipcErrors.push_back(
+            0.5 * std::abs(ipcshift.x - shift.y) + 0.5 * std::abs(ipcshift.y - shift.x));
+
+        // error should be low
+        EXPECT_NEAR(ipcshift.x - shift.x, 0.0, 0.1);
+        EXPECT_NEAR(ipcshift.y - shift.y, 0.0, 0.1);
+    }
+
+    cv::Scalar pcMean, pcStddev, ipcMean, ipcStddev;
+    meanStdDev(ipcErrors, ipcMean, ipcStddev);
+    meanStdDev(pcErrors, pcMean, pcStddev);
+
+    // average error should be low
+    ASSERT_LT(ipcMean[0], 0.03);
+    // average error should be less than non-iterative average error
+    ASSERT_LT(ipcMean[0], pcMean[0]);
+    // error stddev should be less than non-iterative error stddev
+    ASSERT_LT(ipcStddev[0], pcStddev[0]);
+}
+
+
+TEST(Imgproc_PhaseCorrelationIterative, 256x128_accuracy)
+{
+    TestPhaseCorrelationIterative(Size(256, 128), 1);
+}
+
+TEST(Imgproc_PhaseCorrelationIterative, 64x64_accuracy_shift_1)
+{
+    TestPhaseCorrelationIterative(Size(64, 64), 1);
+}
+
+TEST(Imgproc_PhaseCorrelationIterative, 64x64_accuracy_shift_16)
+{
+    TestPhaseCorrelationIterative(Size(64, 64), 16);
+}
+
+TEST(Imgproc_PhaseCorrelationIterative, 0x0_image)
+{
+    ASSERT_ANY_THROW(TestPhaseCorrelationIterative(Size(0, 0), 1));
+}
+
+TEST(Imgproc_PhaseCorrelationIterative, 1x1_image)
+{
+    ASSERT_ANY_THROW(TestPhaseCorrelationIterative(Size(1, 1), 1));
+}
+
+TEST(Imgproc_PhaseCorrelationIterative, accuracy_real_img)
+{
+    Mat img = imread(cvtest::TS::ptr()->get_data_path() + "shared/airplane.png", IMREAD_GRAYSCALE);
+    if (img.empty())
+        return;
+    img.convertTo(img, CV_64FC1);
+
+    const int xLen = 256;
+    const int yLen = 256;
+    const int xShift = 40;
+    const int yShift = 14;
+
+    Mat roi1 = img(Rect(xShift, yShift, xLen, yLen));
+    Mat roi2 = img(Rect(0, 0, xLen, yLen));
+
+    const Point2d ipcShift = phaseCorrelateIterative(roi1, roi2);
+
+    ASSERT_NEAR(ipcShift.x, (double)xShift, 1.);
+    ASSERT_NEAR(ipcShift.y, (double)yShift, 1.);
+}
+
+}}  // namespace opencv_test