OpenCV 最小二乘+距离最小拟合圆
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()); // 随机 engineuniform_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);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);circle(m, Point2f(200, 200), 100, CV_RGB(0, 255, 0), 1);
}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=∣(x−x0)2+(y−y0)2−R∣Loss = |\sqrt{(x - x_0 )^2 + (y - y_0 )^2} - R| Loss=∣(x−x0)2+(y−y0)2−R∣
现在要想办法让这个 Loss 最小, 但是这个方程没有直接的方式可以解, 只有一次次迭代来求最优解. 开头已经说了不用 GSL 库, 所以我就用梯度下降法来解. 但是在开始之前迭代之前, 我们将要求解的 圆心和半径初始化为最小二乘拟合的结果, 可以方便计算和收敛. 如果不懂梯度下降算法的话, 可以看一下相关的文章或视频
const int lr = 1; // 学习率, 一般拟合的点都比较多, 所以我设置的比较大, 可以根据你的情况来设置
const int nIters = 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 < nIters; 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, 除非你的点很有特殊性或者所有点就是在圆上. 因为要拟合的这些点不能完全代表原来的标准圆参数, 解本来就是有误差的
OpenCV 最小二乘+距离最小拟合圆相关推荐
- OpenCV入门系列 —— cv::minEnclosingCircle 随机生成点坐标并计算最小包围圆
OpenCV入门系列 -- cv::minEnclosingCircle 随机生成点坐标并计算最小包围圆 前言 程序说明 输出结果 代码示例 前言 随着工业自动化.智能化的不断推进,机器视觉(2D/3 ...
- python最小二乘法拟合圆_最小二乘法拟合圆
有一系列的数据点 {xi,yi}.我们知道这些数据点近似的落在一个圆上.依据这些数据预计这个圆的參数就是一个非常有意义的问题.今天就来讲讲怎样来做圆的拟合.圆拟合的方法有非常多种,最小二乘法属于比較简 ...
- 最小二乘法拟合圆公式推导及vc实现
最小二乘法(least squares analysis)是一种 数学 优化 技术,它通过 最小化 误差 的平方和找到一组数据的最佳 函数 匹配. 最小二乘法是用最简的方法求得一些绝对不可知的真值,而 ...
- python 拟合圆_最小二乘法拟合圆 转
有一系列的数据点 {xi,yi}{xi,yi},我们知道这些数据点近似的落在一个圆上,根据这些数据估计这个圆的参数就是一个很有意义的问题.今天就来讲讲如何来做圆的拟合.圆拟合的方法有很多种,最小二乘法 ...
- python最小二乘法拟合圆_最小二乘法拟合圆(示例代码)
有一系列的数据点 {xi,yi}.我们知道这些数据点近似的落在一个圆上.依据这些数据预计这个圆的參数就是一个非常有意义的问题.今天就来讲讲怎样来做圆的拟合.圆拟合的方法有非常多种,最小二乘法属于比較简 ...
- OpenCV 学习(直线拟合)
Hough 变换可以提取图像中的直线.但是提取的直线的精度不高.而很多场合下,我们需要精确的估计直线的参数,这时就需要进行直线拟合. 直线拟合的方法很多,比如一元线性回归就是一种最简单的直线拟合方法. ...
- opencv4 c++ 提取图片中的白色区域_【从零学习OpenCV 4】点集拟合
本文首发于"小白学视觉"微信公众号,欢迎关注公众号 本文作者为小白,版权归人民邮电出版社发行所有,禁止转载,侵权必究! 经过几个月的努力,小白终于完成了市面上第一本OpenCV 4 ...
- OpenCV基于Python霍夫圆检测—基于梯度的霍夫圆检测
基于梯度的霍夫圆检测 1. 回顾与目标 2. 基于梯度的霍夫圆检测 2.1 问题分析 2.2 基于梯度的霍夫圆检测步骤 3. 基于梯度的霍夫圆检测函数HoughCircles 3.1 函数HoughC ...
- 最小二乘法拟合圆公式推导及其实现
1.1最小二乘拟合圆介绍与推导 最小二乘法(least squares analysis)是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配.最小二乘法是用最简的方法求得一些绝对不 ...
最新文章
- 实战centos6安装zabbix-2.4版(终极版)
- 【图论专题】最小生成树的扩展应用
- 第13届景驰-埃森哲杯广东工业大学ACM程序设计大赛 L-回旋星空
- Android开发之使用Preferences设计软件设置界面(源代码分享)
- KeeperErrorCode = Unimplemented for /test
- [深入学习C#]LINQ查询表达式详解(2)——查询表达式的转换
- 黄聪:ThinkSAAS开发文档 常用函数 模版修改
- 什么是多核电脑?什么是64位电脑?
- 目标检测: Anchor-Free 时代
- Windows核心编程_让窗口跟随系统样式变化
- 【Robot 学习1】 机器人平台搭建
- html注册手机号验证,js正则表达式验证手机号码,用户名和邮箱
- 易语言版{大智慧/分析家/飞狐交易师}DLL插件接口开发模块(beta5),自定义股票软件公式扩展函数
- 图神经网络详解(四)
- 超声波模块SRF05
- 18.8.17 考试总结
- 报Keystore was tampered with, or password was incorret的原因
- Android 360度全景图片 源码
- 软件开发工程师招聘笔试题面试题223套和招聘考察内容
- 互联网下半场新征程启航,AI、大数据等前沿科技助力传统零售产业转型
热门文章
- JS高级程序设计(14)
- 实用 通用Adb无线调试开发Android应用程序
- 黑暗背景(所有暗主题cobalt,dracula...)Rstudio查看对象窗口viewer没有滚动条,白亮背景就有(所有白主题chrome,cloud)。R版本[64-bit] R-3.6.0
- 计算机专业的研究生专业方向
- 常用稳压二极管参数表
- MOOC(幕课)的到来!
- python 处理excel文件,按某一列值生成多个excel文件
- 微软计算机键盘上Tab,电脑键盘上的tab键有什么功能
- java中接口的优点和缺点
- Chrome浏览器双击无反应