基于OpenCV相位相关算法解决图像的刚性平移问题

参考文章：https://blog.csdn.net/u010440456/article/details/95343110

在做光学测试的实验，计算散斑位移，遂用到此方法。

参考文章中有阐述原理和matlab代码等。

下面是可实现的代码，用了库函数和 Phase Correlation自实现 两种对比

#include <stdio.h>
#include <iostream>
#include "fftw3.h"
#include "vector"

#include <opencv2/opencv.hpp>
#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"


using namespace cv;
using namespace std;

void PhaseCorrelation2D(const Mat &signal,//原图像信号
    const  Mat &pattern,//带配准图像信号
    int &height_offset,//用来获取估计的高度偏移量
    int &width_offset);//获取宽的偏移量


//本程序做了一个很简单的测试：计算“配准样图”中那些样图的偏移量
//（这些图都是通过matlab理想化的平移图）    
//int main()
//{
//    Mat srcImage11 = imread("2-al-60-15.bmp");
//    Mat srcImage21 = imread("2-al-60-16.bmp");
//
//    Mat srcImage1, srcImage2;
//    cvtColor(srcImage11, srcImage1, CV_BGR2GRAY);
//    cvtColor(srcImage21, srcImage2, CV_BGR2GRAY);
//
//    int height_offset = 0;
//    int width_offset = 0;
//    PhaseCorrelation2D(srcImage1, srcImage2, height_offset, width_offset);//宽的偏移量
//    cout << "Phase Correlation法 :  高度偏移量" << -height_offset << "   宽度偏移量" << -width_offset << endl;
//    getchar();
//    return 0;
//}

int main(int argc, char** argv)
{
    Mat srcImage11 = imread("2-al-60-19.bmp");
    Mat srcImage21 = imread("2-al-60-20.bmp");
    Mat dst1, dst2;
    Mat srcImage1, srcImage2;
    cvtColor(srcImage11, srcImage1, CV_BGR2GRAY);
    cvtColor(srcImage21, srcImage2, CV_BGR2GRAY);


    srcImage1.convertTo(dst1, CV_32FC1);
    srcImage2.convertTo(dst2, CV_32FC1);


    Point2d phase_shift;
    phase_shift = phaseCorrelate(dst1, dst2);
    cout << "OpneCV库函数实现 :" << endl << "\tX方向的偏移 : " << phase_shift.x << ",\tY方向的偏移 : " << phase_shift.y << endl;


    int height_offset = 0;
    int width_offset = 0;
    PhaseCorrelation2D(srcImage1, srcImage2, height_offset, width_offset);//宽的偏移量  
    cout << "Phase Correlation自实现 :" << endl << "高度偏移量" << -height_offset << ",   宽度偏移量" << -width_offset << endl;

    waitKey(0);
    getchar();
    return 0;

}

//该函数返回图像的刚性平移量
void PhaseCorrelation2D(const Mat &signal,//原信号
    const  Mat &pattern,//带配准信号
    int &height_offset,//高的偏移量
    int &width_offset)//宽的偏移量
{
    int nRow = signal.rows;
    int nCol = signal.cols;
    fftw_complex *signal_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);
    fftw_complex *pattern_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);

    for (int i = 0; i < nRow; i++)
    {
        for (int j = 0; j < nCol; j++)
        {
            signal_img[i*nCol + j][0] = signal.at<uchar>(i, j);
            signal_img[i*nCol + j][1] = 0;
            pattern_img[i*nCol + j][0] = pattern.at<uchar>(i, j);
            pattern_img[i*nCol + j][1] = 0;
        }
    }

    // 对两幅图像进行傅里叶变换
    fftw_plan signal_forward_plan = fftw_plan_dft_2d(nRow, nCol, signal_img, signal_img,
        FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan pattern_forward_plan = fftw_plan_dft_2d(nRow, nCol, pattern_img, pattern_img,
        FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(signal_forward_plan);
    fftw_execute(pattern_forward_plan);

    // 求互功率谱
    fftw_complex *cross_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);
    double temp;
    for (int i = 0; i < nRow*nCol; i++)
    {
        cross_img[i][0] = (signal_img[i][0] * pattern_img[i][0]) -
            (signal_img[i][1] * (-1.0*pattern_img[i][1]));
        cross_img[i][1] = (signal_img[i][0] * (-1.0*pattern_img[i][1])) +
            (signal_img[i][1] * pattern_img[i][0]);
        temp = sqrt(cross_img[i][0] * cross_img[i][0] + cross_img[i][1] * cross_img[i][1]) + 0.001;
        cross_img[i][0] /= temp;
        cross_img[i][1] /= temp;
    }

    // backward fft，对互功率谱求反变换
    // FFTW computes an unnormalized transform,
    // in that there is no coefficient in front of
    // the summation in the DFT.
    // In other words, applying the forward and then
    // the backward transform will multiply the input by n.

    // BUT, we only care about the maximum of the inverse DFT,
    // so we don't need to normalize the inverse result.

    // the storage order in FFTW is row-order
    fftw_plan cross_backward_plan = fftw_plan_dft_2d(nRow, nCol, cross_img, cross_img,
        FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(cross_backward_plan);

    // free memory
    fftw_destroy_plan(signal_forward_plan);
    fftw_destroy_plan(pattern_forward_plan);
    fftw_destroy_plan(cross_backward_plan);
    fftw_free(signal_img);
    fftw_free(pattern_img);

    double *cross_real = new double[nRow*nCol];
    for (int i = 0; i < nRow*nCol; i++)
        cross_real[i] = cross_img[i][0];
    //根据反变换求出平移量
    int max_loc = 0;//准备存放最大值的位置坐标（注意，只有一个值）
    double max_vlaue = 0.0;
    for (int i = 0; i < nRow*nCol; i++)
    {
        if (max_vlaue < cross_real[i])
        {
            max_vlaue = cross_real[i];
            max_loc = i;
        }
    }

    height_offset = (int)floor(((int)max_loc) / nCol);
    width_offset = (int)max_loc - nCol * height_offset;

    if (height_offset > 0.5*nRow)
        height_offset = height_offset - nRow;
    if (width_offset > 0.5*nCol)
        width_offset = width_offset - nCol;

    delete[] cross_real;
    cross_real = NULL;
}

需要自行配置FFTW库，前面的文章有配置教程