From c5d70a7f229ed3d9eae75fe2258d3b5ece507fb5 Mon Sep 17 00:00:00 2001 From: zdenyhraz Date: Sat, 13 Dec 2025 10:45:05 +0100 Subject: [PATCH] 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 | 10 + modules/imgproc/include/opencv2/imgproc.hpp | 16 ++ modules/imgproc/src/phasecorr_iterative.cpp | 191 ++++++++++++++++++++ modules/imgproc/test/test_ipc.cpp | 117 ++++++++++++ 4 files changed, 334 insertions(+) create mode 100644 modules/imgproc/src/phasecorr_iterative.cpp create mode 100644 modules/imgproc/test/test_ipc.cpp diff --git a/doc/opencv.bib b/doc/opencv.bib index 97ad726dbd..5c8898f120 100644 --- a/doc/opencv.bib +++ b/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} +} diff --git a/modules/imgproc/include/opencv2/imgproc.hpp b/modules/imgproc/include/opencv2/imgproc.hpp index 35420b6f53..e02d57875e 100644 --- a/modules/imgproc/include/opencv2/imgproc.hpp +++ b/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) diff --git a/modules/imgproc/src/phasecorr_iterative.cpp b/modules/imgproc/src/phasecorr_iterative.cpp new file mode 100644 index 0000000000..ced9339c87 --- /dev/null +++ b/modules/imgproc/src/phasecorr_iterative.cpp @@ -0,0 +1,191 @@ +#include "precomp.hpp" +#include + +namespace { + +template +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>(row); + const auto* dft1p = dft1.ptr>(row); + const auto* dft2p = dft2.ptr>(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(dft1, dft2, cps); + else if (dft1.type() == CV_64FC2) + calculateCrossPowerSpectrum(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(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(L3peak.x - L2size / 2), + static_cast(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(L2Upeak.x - L1size / 2), + static_cast(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(L2Usize) / L2size); + } + } + + // iterative refinement failed to converge, return non-iterative subpixel shift + return getSubpixelShift(L3, Point2d(L3peak), L3mid, L2size); +} diff --git a/modules/imgproc/test/test_ipc.cpp b/modules/imgproc/test/test_ipc.cpp new file mode 100644 index 0000000000..321b1f2fe3 --- /dev/null +++ b/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 + +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(size.width * 0.1), static_cast(size.height * 0.1)), + Point(static_cast(size.width * 0.9), static_cast(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 pcErrors; + std::vector 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_(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