图像配准经典算法

一、图像配准的定义

对包含同一场景的两幅或者多幅图像进行空间位置上的对齐， 使得多个图像中相同的特征尽量处在相同的位置上

二、图像配准的算法流程

几何变换

线性变换

通常包含平移、旋转、刚体运 动以及仿射变换等

以投影变换为例处理图像存在倾斜的场景，两张影像之间没有平行线。变换矩阵中存在8个未知参数，需要4对控制点才能进行求解，变换矩阵如下所示：

非线性变换

多项式变换多项式变换是一种典型的非线性变换模型，多项式的阶数越高，拟合效果越 好，然而需要的控制点对数量也随着增加

如下所示，二阶多项式变换模型公式存在12个未知数，需要6对控制点进行多项式参数计算。

特征点提取

Harris特征检测

通过在局部窗口中计算像素点在任意方向上的灰度变化量来寻找角点特征，若在任意方向移动窗口时灰度值均有明显变化，那么当前的像素点就被判别为角点。

ROEWA滤波模板

/*********该函数根据尺度和窗口半径生成ROEWA滤波模板************/
/*size表示核半径，因此核宽度是2*size+1
 scale表示指数权重参数
 kernel表示生成的滤波核
 */

static void roewa_kernal(int size, float scale, Mat& kernal)
{
	kernal.create(2 * size + 1, 2 * size + 1, CV_32FC1);
	for (int i = -size; i <= size; ++i)
	{
		// 里i + size的目的是将索引范围从[0, 2 * size]映射到[-size, size]。
		float* ptr_k = kernal.ptr<float>(i + size);
		for (int j = -size; j <= size; ++j)
		{
			ptr_k[j + size] = exp(-1.f * (abs(i) + abs(j)) / scale);
	
		}
	}
}

高斯滤波

 /*该函数根据窗口半径，标准差，生成圆形高斯滤波模板，不是正方向*/
/*size表示圆半径
 scale表示高斯标准差
 gauss_kernel表示生成的圆形高斯核
 二维高斯函数的形式：1/(2*pi*scale*scale)*exp(-(x*x+y*y)/(2*scale*scale))
 */

static void gauss_circle(int size, float scale, Mat& gauss_kernal)
{
	gauss_kernal.create(2 * size + 1, 2 * size + 1, CV_32FC1);
	float exp_temp = -1.f / (2 * scale * scale);
	float sum = 0;
	for (int i = -size; i <= size; ++i)
	{
	// 由于索引值没有负值，因此需要将索引+-size映射到0-2*size坐标系内，
		float* ptr = gauss_kernal.ptr<float>(i + size);
		for (int j = -size; j <= size; ++j)
		{
			if ((i * i + j * j) <= size * size)
			{
				float value = exp((i * i + j * j) * exp_temp);
				sum += value;
				ptr[j + size] = value;
			}
			else
				ptr[j + size] = 0.f;
		}
	}

	for (int i = -size; i <= size; ++i)
	{
		float* ptr = gauss_kernal.ptr<float>(i + size);
		for (int j = -size; j <= size; ++j)
		{
			ptr[j + size] = ptr[j + size] / sum;
		}
	}
}

计算Harris角点响应值
使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

/*************该类成员函数构建SIFT尺度空间*****************/
/*image表示输入的原始图像
 sar_harris_fun表示尺度空间的Sar_harris函数
 amplit表示尺度空间像素的梯度幅度
 orient表示尺度空间像素的梯度方向
 */
void shift::build_shift_space(const Mat& image, vector<Mat>& sar_harris_func, vector<Mat>& amplit, vector<Mat>& orient)
{
	// 转换输入图像格式
	Mat gray_img;
	if (image.channels() != 1)
		cvtColor(image, gray_img, CV_RGB2GRAY);
	else
		gray_img = image;
	// 更改图像类型
	Mat float_img;
	gray_img.convertTo(float_img, CV_32FC1, 1, 0.f);// 将整数数据转换位浮点数据数据

	//定义多尺度特征，相当于在第三维度repeat
	sar_harris_func.resize(Mmax);
	amplit.resize(Mmax);
	orient.resize(Mmax);

	for (int i = 0; i < Mmax; ++i)
	{
		//根据sigma和ratio计算
		float scale = (float)sigma * (float)pow(ratio, i);
		int radius = cvRound(2 * scale);
		Mat kernal;
		// 定义高斯滤波卷积核
		roewa_kernal(radius, scale, kernal);
		// 将高斯滤波分解为四个方向滤波
		// 四个滤波模版
		Mat W34 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);
		Mat W12 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);
		Mat W14 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);
		Mat W23 = Mat::zeros(2 * radius + 1, 2 * radius + 1, CV_32FC1);
		// 矩阵索引赋值实现滤波模版窗口的设定
		// 处理部分行与全部列
		kernal(Range(radius + 1, 2 * radius + 1), Range::all()).copyTo(W34(Range(radius + 1, 2 * radius + 1), Range::all()));
		kernal(Range(0, radius), Ra