最小二乘法曲面拟合及C++代码，用于求取图像定位的亚像素级精度_曲面拟合求亚像素点位置

曲面拟合是指在三维坐标系中，将多个点拟合为一个曲面的过程。曲面拟合的主要作用是可以找到极值点。在图像定位领域，用于寻找亚像素级精度。

在本里中，曲面的公式采用：

（也可以采用其他的曲面公式，原理一样）。
并使用最小二乘法进行求解：是下面式子的误差最小：
则应该满足：







写为矩阵的表达式：

下一步就是要根据拟合点的坐标信息，解出[a,b,c,d,e,f]。在本例子中，共选取9个点，对曲面进行拟合（具体可以根据实际情况进行选择）。9点的坐标可以表示为。

则[a,b,c,d,e,f]的解可以直接根据9个点的z坐标进行解出。C++代码如下所示：（以下代码是采用opencv库进行的矩阵运算，进攻参考）。

#include <iostream>

#include <vector>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
using std::cout;
using std::endl;
using namespace cv;

//曲面拟合函数，最终得到亚像素级精度
Point SurfaceFitting(Mat &CorrResult);

Point SurfaceFitting(Mat &CorrResult)
{
	/*
	原始坐标矩阵
	float matrixdata[] = {  51,25,27,27,15,15,
						    25,51,27,15,27,15,
							27,27,25,15,15, 9,
							27,15,15,15, 9, 9,
							15,27,15, 9,15, 9,
							15,15, 9, 9, 9, 9 };
	*/
	double Time = (double)cvGetTickCount();
	//此值为matrixdata矩阵计算出来的,可以作为初始参数使用
	float Edata[36] = { 0.50000107, 1.7764522e-07, -2.3374371e-08, -1.0000021, -3.3659094e-07, 0.16666715,
						-2.7998499e-07, 0.50000054, 8.6149221e-08, 4.4305315e-07, -1.0000013, 0.16666681,
						-1.5373146e-08, 2.3572157e-07, 0.25000003, -0.25000003, -0.25000057, 0.25000012,
						-1.0000023, -6.5050045e-07, -0.25, 2.4166713, 0.25000131, -0.75000113,
						5.8114591e-07, -1.0000014, -0.25000021, 0.24999917, 2.4166701, -0.75000048,
						0.166667, 0.16666727, 0.25000006, -0.75000083, -0.75000137, 0.805556 };

	Mat E(6, 6, CV_32F, Edata);


	float n1, n2, n3, n4, n5, n6;

	n1 = CorrResult.at<float>(0, 1) + 4 * CorrResult.at<float>(0, 2)
		+ CorrResult.at<float>(1, 1) + 4 * CorrResult.at<float>(1, 2)
		+ CorrResult.at<float>(2, 1) + 4 * CorrResult.at<float>(2, 2);

	n2 = CorrResult.at<float>(1, 0) + CorrResult.at<float>(1, 1) + CorrResult.at<float>(1, 2)
		+ 4 * CorrResult.at<float>(2, 0)+4 * CorrResult.at<float>(2, 1) + 4 * CorrResult.at<float>(2, 2);

	n3 = CorrResult.at<float>(1, 1) + 2 * CorrResult.at<float>(1, 2)
		+ 2 * CorrResult.at<float>(2, 1) + 4 * CorrResult.at<float>(2, 2);

	n4= CorrResult.at<float>(0, 1) + 2 * CorrResult.at<float>(0, 2)
		+ CorrResult.at<float>(1, 1) + 2 * CorrResult.at<float>(1, 2)
		+ CorrResult.at<float>(2, 1) + 2 * CorrResult.at<float>(2, 2);

	n5 = CorrResult.at<float>(1, 0) + CorrResult.at<float>(1, 1) + CorrResult.at<float>(1, 2)
		+ 2 * CorrResult.at<float>(2, 0) + 2 * CorrResult.at<float>(2, 1) + 2 * CorrResult.at<float>(2, 2);

	n6 = CorrResult.at<float>(0, 0) + CorrResult.at<float>(0, 1) + CorrResult.at<float>(0, 2)
		+ CorrResult.at<float>(1, 0) + CorrResult.at<float>(1, 1) + CorrResult.at<float>(1, 2)
		+ CorrResult.at<float>(2, 0) + CorrResult.at<float>(2, 1) + CorrResult.at<float>(2, 2);



	float Ndata[6] = { n1, n2, n3, n4, n5, n6 };

	Mat N(6, 1, CV_32F, Ndata);

	Mat A(6, 1, CV_32F);
	Time = (double)cvGetTickCount() - Time;
	printf("run time = %gms\n", Time / (cvGetTickFrequency() * 1000));//毫秒
	A = E*N;    //主要是这一步非常耗时

	
	float a_ = A.ptr<float>(0)[0];
	float b_ = A.ptr<float>(1)[0];
	float c_ = A.ptr<float>(2)[0];
	float d_ = A.ptr<float>(3)[0];
	float e_ = A.ptr<float>(4)[0];

	float max_xpod = (2 * b_*d_ - c_*e_) / (c_*c_ - 4 * a_*b_);
	float max_ypod = (2 * a_*e_ - d_*c_) / (c_*c_ - 4 * a_*b_);

	Point result;
	result.x = max_xpod;
	result.y = max_ypod;

	return result;
}

int main()
{
	//9个点的z坐标
	float corrresult[9] = { 10,20,10,
							40,100,50,
							10,80,40 };
	//相关值计算构成的矩阵
	Mat CorrResult(3, 3, CV_32F, corrresult);

	Point result = SurfaceFitting(CorrResult);


	cout << result << endl;

	return 0;

}
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

总结：在使用最小二乘法进行曲线拟合的原理和曲面拟合的原理是一样的。本人使用曲面拟合主要是在图像匹配中像得到亚像素级的精度，事实证明效果不错。