图像配准经典算法_多项式图像配准原理

一、图像配准的定义

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

二、图像配准的算法流程

几何变换

线性变换

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

以投影变换为例处理图像存在倾斜的场景，两张影像之间没有平行线。变换矩阵中存在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);
	
		}
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

高斯滤波

 /*该函数根据窗口半径，标准差，生成圆形高斯滤波模板，不是正方向*/
/*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)
	{
		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;
		}
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

计算Harris角点响应值

/*************该类成员函数构建SAR_SIFT尺度空间*****************/
/*image表示输入的原始图像
 sar_harris_fun表示尺度空间的Sar_harris函数
 amplit表示尺度空间像素的梯度幅度
 orient表示尺度空间像素的梯度方向
 */
void shift::build_sar_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), Range::all()).copyTo(W12(Range(0, radius), Range::all()));
		// 处理部分列与全部行
		kernal(Range::all(), Range(radius+1,2*radius+1)).copyTo(W14(Range::all(), Range(radius + 1, 2 * radius + 1)));
		kernal(Range::all(), Range(0,radius)).copyTo(W23(Range::all(), Range(0, radius)));

		// 执行滤波函数
		Mat M34, M12, M14, M23;
		double eps = 0.00001;
		filter2D(float_img, M34, CV_32FC1, W34, Point(-1, -1), eps);
		filter2D(float_img, M12, CV_32FC1, W12, Point(-1, -1), eps);
		filter2D(float_img, M14, CV_32FC1, W14, Point(-1, -1), eps);
		filter2D(float_img, M23, CV_32FC1, W23, Point(-1, -1), eps);
		// 计算水平梯度与竖直梯度
		Mat Gx, Gy;
		log((M14) / (M23), Gx);
		log((M34) / (M12), Gy);

		// 计算幅度和相位方向
		magnitude(Gx, Gy, amplit[i]);
		phase(Gx, Gy, orient[i], true);

		//计算harris角点响应值
		Mat Csh_11 = scale * scale * Gx.mul(Gx);
		Mat Csh_12 = scale * scale * Gx.mul(Gy);
		Mat Csh_22 = scale * scale * Gy.mul(Gy);

		// 高斯权重
		float gauss_sigma = sqrt(2.f) * scale;
		int size = cvRound(3 * gauss_sigma);

		Size Kern_size(2 * size + 1, 2 * size + 1);
		// 使用自定义核的高斯滤波
		GaussianBlur(Csh_11, Csh_11, Kern_size, gauss_sigma);
		GaussianBlur(Csh_12, Csh_12, Kern_size, gauss_sigma);
		GaussianBlur(Csh_22, Csh_22, Kern_size, gauss_sigma);

		Mat Csh_21 = Csh_12;

		// 计算harris角点响应值
		Mat temp_add = Csh_11 + Csh_22;//迹
		// 行列式：Csh_11.mul(Csh_22) - Csh_21.mul(Csh_12)
		sar_harris_func[i] = Csh_11.mul(Csh_22) - Csh_21.mul(Csh_12) - (float)d * temp_add.mul(temp_add);
	}


}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

LoG特征检测

DoG特征检测

在某一 尺度层中采用两个相邻高斯尺度空间的图像作差得到响应图像。然后基于局部极值搜索策略在响应图像进行特征点检测。通过当前点的尺度空间和空间位置检测局部特征点。

1. 构建高斯金字塔
2. 特征点检测
3. 极值点选取（三维二次函数）

描述符构建

利用特征点周边像素特征构建特征描述向量对特征点进行表示所表示的信息需要具有显著性和代表性，既能够反映出不同特征点之间的差异，又能够反映出参考图像与待测图像之间相同特征点之间的一致性

SIFT描述符

又称为尺度不变特征变换描述符，是特征点周围梯度幅值与方向信息的统计量

计算步骤如下：

1. 首先统计特征点主方向并将特征点所在的邻域块旋转到主方向
   首先在高斯金字塔图像3σ邻域窗口中计算每个像素的梯度幅值和方向信息。然后将360°方向划分到36个方向块中，并在特征点所在的邻域中对每个像素的梯度幅值和方向信息进行统计，将幅值累加和最大的方向判定为特征点主方向，同时将主方向峰值80%的方向设定为辅方向
2. 将特征点所在的邻域块划分为4×4的子区域块，在每个子块中进行八个方向的梯度方向直方图统计， 最终生成128维的特征向量，构成描述符
   

SURF描述符

统计子区域块中所有像素点在水平和垂直方向上的Haar小波响应值之和以及响应值的绝对值之和

形状上下文

统计以任一点为中心，同心圆与径向区范围内的特征点数量
采用边缘检测算子提取形状特征，对特征进行采样，选取其中任意一点的坐标作为极坐标系中点，投建N个同心圆和N个径向区域通过统计不同位置点位数量来描述该点，得到一个M*N的特征向量。

特征匹配

通过特征描述符进行特征点之间的相似性度量，来寻找正确匹配的特征点位

NNDR匹配

通过描述符之间距离的计算，得到最近邻距离(有最小欧式距离的特征点)和次近邻距离之间的比值，将该比值的大小作为衡量标准进行描述符之间的相似性度量

基于k-d树的最近邻匹配

局部敏感哈希匹配

总结

2024年02月25日更新
参考文档如下：
1.可见光与SAR关键技术研究
2.图像配准代码