1,对绝对理想的图像进行配准模拟

（代码见后面）

绝对理想的图像：图中整体的图像都没有因为平移而带来像素的损失，即再往下移动越界后像素偏移到了上面，往右移动越界后像素偏移到左边，所以图片中模拟的图像没有任何像素损失。

%读取
Image1 = (imread('image_gray.jpg')); %
Image2 = (imread('img_result_8_15.jpg')); % 带有偏移量的图像，高度向下偏移8，宽度向右偏移15%显示
subplot(1,2,1);imshow(Image1);  title('原参考图像');
subplot(1,2,2);imshow(Image2);  title('带有偏移量的图像');FFT1 = fft2(Image1); % 二维 FFT
FFT2 = conj(fft2(Image2)); %共轭复数  %求复功率谱
FFTR = FFT1.*FFT2; % 复共轭（复功率谱的分子）
magFFTR = abs(FFTR); %sqrt(real^2 + imag^2)取模 （复功率谱的分母）
FFTRN = (FFTR./magFFTR);  %复功率谱的反变换
result = ifft2(double(FFTRN));   M = max(max(result));
[i,j] = find(result == M)figure;  colormap('gray');imagesc(result);
figure;  colormap(jet);  mesh(result);  

int main(int argc, char** argv)
{Mat srcImage11 = imread("Brain13x17y.bmp");Mat srcImage21 = imread("BrainP.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);QueryPerformanceFrequency(&freq);QueryPerformanceCounter(&startTime);//开始计时Point2d phase_shift;phase_shift = phaseCorrelate(dst1, dst2);cout << "OpneCV库函数实现 :" << endl << "\tX方向的偏移 : " << phase_shift.x << ",\tY方向的偏移 : " << phase_shift.y << endl;QueryPerformanceCounter(&stopTime);time = 1e3*(stopTime.QuadPart - startTime.QuadPart) / freq.QuadPart;cout << "你电脑的频率为：" << freq.QuadPart << endl;cout << "opencv库函数消耗时间为：" << time << "毫秒" << endl;QueryPerformanceCounter(&startTime);//开始计时int height_offset = 0;int width_offset = 0;PhaseCorrelation2D(srcImage1, srcImage2, height_offset, width_offset);//宽的偏移量  cout << "Phase Correlation自实现 :" << endl << "高度偏移量" << -height_offset << ",   宽度偏移量" << -width_offset << endl;QueryPerformanceCounter(&stopTime);//结束计时time = 1e3*(stopTime.QuadPart - startTime.QuadPart) / freq.QuadPart;cout << "自实现相位相关函数消耗时间为：" << time << "毫秒" << endl;waitKey(0);getchar();return 0;
}

参考代码3（需要配置opencv库和MIT的fftw库）：

#include <stdio.h>
#include <iostream>
#include "fftw3.h"
#include "vector"#include <opencv2/legacy/legacy.hpp>
#include <opencv2/nonfree/nonfree.hpp>//opencv_nonfree模块：包含一些拥有专利的算法，如SIFT、SURF函数源码。
#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <opencv2/nonfree/features2d.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("image_gray.jpg");Mat srcImage21 = imread("img_result_-8_-25.jpg");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;
}//该函数返回图像的刚性平移量
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-orderfftw_plan cross_backward_plan = fftw_plan_dft_2d(nRow, nCol, cross_img, cross_img,FFTW_BACKWARD, FFTW_ESTIMATE);fftw_execute(cross_backward_plan);// free memoryfftw_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;
}

文章收录：

这俩天一直在做关于物体匹配的方面的工作，前几天朋友推荐我看西安电子科技大学张瑞娟的一篇硕士论文“图像配准理论及算法研究”，我收获很大，所以我也总结一些对我有用的算法，将来便于查找应用。
我做的目标追踪这一块，虽然和图像配准不是一个方向，但是前期工作都是一样的，首先都需要物体检测，特征检测和匹配。这里我总结一些对我有用的，也希望对和我一样研究方向的人有帮助。这里图像配准可以换成物体匹配的。
1， 图像配准要素结合：特征空间，搜索空间，搜索策略，近似性度量
2， 图像配准方法：
2.1基于灰度信息的方法，
交叉相关（互相关）方法，相关系数度量，序贯相似检测算法，信息理论的交换信息相似性准则
2.2基于变换域的方法
相位相关法，Walsh Transform变换
2.3基于特征的方法
常用的图像特征有：特征点（包括角点、高曲率点等）、直线段、边缘（Robert、高斯-拉普拉斯LoG、Canny、Gabor滤波等边缘检测算子）或轮廓、闭合区域、特征结构以及统计特征如矩不变量等
注：像素灰度信息的互相关算法相比，特征提取包含了高层信号信息，所以该类算法对光照、噪声等的抗干扰能力强。
3，常用的空间变换模型
刚体变换（平移、旋转与缩放的组合）、仿射变换、透射变换、投影变换、非线性变换
4， 常用的相似性测度
4.1距离测度
均方根误差，差绝对值和误差，兰氏距离，Mahalanobis距离，绝对差，Hausdorff距离等
4.2角度度量法（概率测度）。
4.3 相关度量法
5，配准算法的评价标准
配准时间、配准率、算法复杂度、算法的可移植性、算法的适用性、图像数据对算法的影响等（这里虽然不是目标追踪的评价标准，但是我们可以借鉴这些评价算法的标准）

收录文章2：

cv::Point2d cv::phaseCorrelateRes(InputArray _src1, InputArray _src2, InputArray _window, double* response)
{  //分别得到两幅输入图像和窗函数的矩阵形式  Mat src1 = _src1.getMat();  Mat src2 = _src2.getMat();  Mat window = _window.getMat();  //输入图像的类型和大小的判断，必须一致，而且类型必须是32位或64位浮点灰度图像  CV_Assert( src1.type() == src2.type());  CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );  CV_Assert( src1.size == src2.size);  //如果使用窗函数，则窗函数的大小和类型必须与输入图像的一致  if(!window.empty())  {  CV_Assert( src1.type() == window.type());  CV_Assert( src1.size == window.size);  }  //因为要进行离散傅立叶变换，所以为了提高效率，就要得到最佳的图像尺寸  int M = getOptimalDFTSize(src1.rows);  int N = getOptimalDFTSize(src1.cols);  Mat padded1, padded2, paddedWin;  //生成尺寸修改以后的矩阵  if(M != src1.rows || N != src1.cols)   //最佳尺寸不是原图像的尺寸  {  //通过补零的方式填充多出来的像素  copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));  copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));  if(!window.empty())  {  copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));  }  }  else    //最佳尺寸与原图像的尺寸一致  {  padded1 = src1;  padded2 = src2;  paddedWin = window;  }  Mat FFT1, FFT2, P, Pm, C;  // perform window multiplication if available  //执行步骤1，两幅输入图像分别与窗函数逐点相乘  if(!paddedWin.empty())  {  // apply window to both images before proceeding...  multiply(paddedWin, padded1, padded1);  multiply(paddedWin, padded2, padded2);  }  // execute phase correlation equation  // Reference: http://en.wikipedia.org/wiki/Phase_correlation  //执行步骤2，分别对两幅图像取傅立叶变换  dft(padded1, FFT1, DFT_REAL_OUTPUT);  dft(padded2, FFT2, DFT_REAL_OUTPUT);  //执行步骤3  //计算互功率谱的分子部分，即公式3中的分子，其中P为输出结果，true表示的是对FF2取共轭，所以得到的结果为：P=FFT1×FFT2*，mulSpectrums函数为通用函数  mulSpectrums(FFT1, FFT2, P, 0, true);  //计算互功率谱的分母部分，即公式3中的分母，结果为：Pm=|P|，magSpectrums函数就是在phasecorr.cpp文件内给出的，它的作用是对复数取模。  magSpectrums(P, Pm);  //计算互功率谱，即公式3，结果为：C=P / Pm，divSpectrums函数也是在phasecorr.cpp文件内给出的，它仿照mulSpectrums函数的写法，其中参数false表示不取共轭  divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)  //执行步骤4，傅立叶逆变换  idft(C, C); // gives us the nice peak shift location...  /*平移处理，fftShift函数也是在phasecorr.cpp文件内给出的，它的作用是把图像平均分割成——左上、左下、右上、右下，把左上和右下对调，把右上和左下对调。它的目的是把能量调整到图像的中心，也就是图像的中心对应于两幅图像相频差为零的地方，即没有发生位移的地方。*/  fftShift(C); // shift the energy to the center of the frame.  //执行步骤5  // locate the highest peak  //找到最大点处的像素位置，minMaxLoc为通用函数  Point peakLoc;  minMaxLoc(C, NULL, NULL, NULL, &peakLoc);  // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...  //在5×5的窗体内确定亚像素精度的坐标位置  Point2d t;  // weightedCentroid也是在phasecorr.cpp文件内给出的，它是利用公式4来计算更精确的坐标位置  t = weightedCentroid(C, peakLoc, Size(5, 5), response);  // max response is M*N (not exactly, might be slightly larger due to rounding errors)  //求最大响应值  if(response)  *response /= M*N;  // adjust shift relative to image center...  //最终确定位移量  //先找到图像中点，然后用中点减去由步骤5得到的坐标位置  Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);  return (center - t);
}