OpenCV 最小二乘+距离最小拟合圆

我们经常需要由给定的点精确地拟合出一个圆, 下面讲解从 最小二乘算法 到 距离最小算法 的实现过程, 其中 距离最小算法 不使用 GSL 算法库

一. 最小二乘算法

最小二乘算法有两种方法可以实现, 一种是微积分算法, 一种是矩阵投影算法, 由于 OpenCV 中已经提供了矩阵运算功能, 所以这里使用矩阵方法来实现最小二乘拟合圆

假设已经有了一系列在圆上或者其附近的点, 想由这些点找出圆心和半径, 怎么实现呢? 我们先用代码生成需要拟合的点

所需要的头文件

#include <random>
#include <ctime>

using namespace std;

#include <opencv2\\opencv.hpp>
#include <opencv2\\features2d.hpp>

using namespace cv;

生成代码

vector<Point2f> pts;				// 代码生成的点
const Point2f center(200, 200);		// 圆心
const float std_r = 100;			// 标准半径

// 从 [0, 360) 生成圆周上的点, 都是利用 pt_start 旋转, 再加上一点随机偏移
for (int i = 0; i < 360; i += 8)
{
	default_random_engine e(time(nullptr) + rand());	// 随机 engine
	uniform_real_distribution<float> u(-8.0f, 8.0f);

	const float r = std_r + u(e);						// 半径 + 误差

	const Point2f pt_start(center.x + r, center.y);		// 旋转起点
	const Point2f ofst = pt_start - center;

	const float rad = (float)(i * CV_PI / 180);
	const float sin_val = (float)sin(rad);
	const float cos_val = (float)cos(rad);

	const Point2f delta = Point2f(ofst.x * cos_val - ofst.y * sin_val, ofst.x * sin_val + ofst.y * cos_val);

	pts.push_back(center + delta);
}

经过上面的步骤后, pts 中就有了 45 个随机点, 也就是需要拟合的数据, 可以显示出来看一下

Mat m = Mat::zeros(400, 400, CV_8UC3);

circle(m, Point2f(200, 200), 100, CV_RGB(0, 255, 0), 1);

const int pt_len = pts.size();

for (int i = 0; i < pt_len; i++)
{
	circle(m, pts[i], 2, CV_RGB(255, 0, 0), FILLED);
}

namedWindow("rand_circle", WINDOW_NORMAL);
imshow("rand_circle", m);

从图上可以看出, 生成的点还是比较随机的, 圆内和圆外都有, 接下来就是拟合了, 理论依据是矩阵投影, 可以看一下《Linear Algebra with Applications》第五章的最小二乘部分, 有一个例题就是关于拟合圆的, 这里就不抄公式了, 直接放码

Mat A(pt_len, 3, CV_32FC1);
Mat b(pt_len, 1, CV_32FC1);

// 下面的两个 for 循环初始化 A 和 b
for (int i = 0; i < pt_len; i++)
{
	float *pData = A.ptr<float>(i);

	pData[0] = pts[i].x * 2.0f;
	pData[1] = pts[i].y * 2.0f;

	pData[2] = 1.0f;
}

float *pb = (float *)b.data;

for (int i = 0; i < pt_len; i++)
{
	pb[i] = pts[i].x * pts[i].x + pts[i].y * pts[i].y;
}

// 下面的几行代码就是解超定方程的最小二乘解
Mat A_Trans;
transpose(A, A_Trans);

Mat Inv_A;
invert(A_Trans * A, Inv_A);

Mat res = Inv_A * A_Trans * b;

// 取出圆心和半径
float x = res.at<float>(0, 0);
float y = res.at<float>(1, 0);
float r = (float)sqrt(x * x + y * y + res.at<float>(2, 0));

上面很简单的代码就完成了, 因为是随机点, 所以你求出来的值可能和我的不一样, 不过相差应该不大. 我算出来的值是
x = 198.658, y = 199.511, r = 99.179
可以看出, 拟合的值和用来生成随机点的圆心和半径很接近, 但是还是不够理想. 主要是生成的点太少, 而且生成点时设置的误差范围(-8.0 ~ 8.0) 比较大, 设置这么大的误差是为了方便在图上可以看出点不在圆上, 你可以设置小一点试试看效果是很好的

最小二乘拟合有一个致命的缺陷就是拟合的结果会向杂点的位置偏移. 如下图, 我在刚才生成的点外面才加了 3 个杂点, 拟合结果就坏掉了. 其中绿色是生成随机点的圆, 黄色是加了杂点后拟合的结果, 直接被 3 个杂点带偏了

为什么最小二乘拟合会受杂点的干扰呢? 因为这三个杂点也是拟合数据的一部分, 最小化误差的时候也会把其考虑进去, 如果拟合出来是绿色的圆, 那这三个点的误差就会特别大. 所以就会向杂点偏移来减小拟合误差. 那这个问题有没有解呢, 让其不受或受杂点干扰很小? 答案是肯定的, 要不然我还写这些干什么呢

二. 距离最小算法

为了解决最小二乘拟合误差的问题, 距离最小拟合误差算法是让所有点到 圆上 的距离最小, 公式如下
$$Loss=|\sqrt{(x-x_0{)}^2+(y-y_0{)}^2}-R|$$
现在要想办法让这个 Loss 最小, 但是这个方程没有直接的方式可以解, 只有一次次迭代来求最优解. 开头已经说了不用 GSL 库, 所以我就用梯度下降法来解. 但是在开始之前迭代之前, 我们将要求解的 圆心和半径初始化为最小二乘拟合的结果, 可以方便计算和收敛. 如果不懂梯度下降算法的话, 可以看一下相关的文章或视频

const int lr = 1;				// 学习率, 一般拟合的点都比较多, 所以我设置的比较大, 可以根据你的情况来设置
const int iters = pt_len;		// 迭代次数, 我设置成了有多少个点就迭代多少次, 也可以根据实际情况设置

vector<float> losses(pt_len);	// 每次迭代后的 loss 值
vector<float> min_loss(pt_len);	// 每次迭代后的最小 loss
vector<float> root_val(pt_len);	// 每次迭代中的开平方值, 方便以后使用

for (int i = 0; i < iters; i++)
{
	float loop_loss = 0;

	for (int j = 0; j < pt_len; j++)
	{
		// 这里第一次迭代的 x, y, r 是最小二乘的结果, 第二次迭代开始就是修正后的结果
		root_val[j] = sqrt((pts[j].x - x) * (pts[j].x - x) + (pts[j].y - y) * (pts[j].y - y));

		const float loss = root_val[j] - r;

		losses[j] = loss;
		loop_loss += fabs(loss);
	}

	min_loss[i] = loop_loss;

	// 如果 loss 值不再减小, 就提前结束
	if (i > 0 && min_loss[i] > min_loss[i - 1])
	{
		break;
	}

	// 下面三个是梯度值
	float gx = 0;
	float gy = 0;
	float gr = 0;
	
	for (int j = 0; j < pt_len; j++)
	{
		// 在计算梯度时要先计算偏导数, 再将 x 代数公式得到
		float gxi = (x - pts[j].x) / root_val[j];

		if (losses[j] < 0)
		{
			gxi *= (-1);
		}

		float gyi = (y - pts[j].y) / root_val[j];

		if (losses[j] < 0)
		{
			gyi *= (-1);
		}

		float gri = -1;

		if (losses[j] < 0)
		{
			gri = 1;
		}

		gx += gxi;
		gy += gyi;
		gr += gri;
	}

	gx /= pt_len;
	gy /= pt_len;
	gr /= pt_len;

	x -= (lr * gx);
	y -= (lr * gy);
	r -= (lr * gr);
}

经过多次迭代之后, 拟合的结果比最小二乘法好了很多, 下图中白色的圆就是拟合结果
x = 200.954, y = 200.627, r = 100.747

三. 还可以优化吗

其实这个也算不上什么优化, 只是在上面迭代后, 把大于或者小于半径太多的点当杂点去掉, 再迭代一波就可以了. 至于大于或小于多少算杂点, 这个可以设置一个比例, 比如 1.05 和 0.95 之类的, 你自己决定. 这里代码我就不演示了, 你自己处理看看, 比较一下结果有没有更好

完成后, 将上面的代码组合成一个函数, 把迭代次数, 学习率, 去除杂点距离比例之类的设置成参数, 就完成一个拟合圆函数了, 是不是很简单!

还有一点要注意, 你几乎不可能拟合得到 x = 200, y = 200, r = 100, 除非你的点很有特殊性或者所有点就是在圆上. 因为要拟合的这些点不能完全代表原来的标准圆参数, 解本来就是有误差的

四. 代码下载

完整的代码可下载 VS2015 x64 代码示例