Repository: opencv/opencv
File: modules/imgproc/src/phasecorr_iterative.cpp (C++)

Code before:
    // 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 the L2 is out of bounds of L3
    int L2size = 7;
    while (isOutOfBounds(L3peak, L3, L2size))
        if (!reduceL2size(L2size))
            return Point2d(L3peak) - L3mid;

    // extract the L2 maximum correlation neighborhood
    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 = 357;
    resize(L2, L2U, {L2Usize, L2Usize}, 0, 0, INTER_LINEAR);
    const Point2d L2Umid(L2U.cols / 2, L2U.rows / 2);

    const int MaxIter = 10;  // maximum number of iterative refinement iterations
    const double L1ratioBase = 0.45;  // base ratio of L1 size (% of L2Usize)
    const double L1ratioStep = 0.025;  // reduction step if failed to converge

    // run the iterative refinement algorithm using the specified L1 ratio,
    // gradually decrease L1 ratio if convergence is not achieved
    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 < MaxIter; ++iter)
        {
            // verify that the L1 region is withing the upsampled L2U region
            if (isOutOfBounds(L2Upeak, L2U, L1size))
                break;

            // extract the L1 region
            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

Review comment (refactor):
It makes sense to extract it as function interface. OpenCV has TermCriteria for it: https://docs.opencv.org/4.x/d9/d5d/classcv_1_1TermCriteria.html

Code after:
    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