最小二乘法拟合圆心

文章为个人学习过程中笔记，原理部分参考其他作者内容，侵权必删
最小二乘法(least squares analysis)是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小来寻找一组数据的最佳匹配函数的计算方法，最小二乘法通常用于曲线拟合 (least squares fitting) 。最小二乘圆拟合方法是一种基于统计的检测方法，即便是图像中圆形目标受光照强度不均等因素的影响而产生边缘缺失，也不会影响圆心的定位和半径的检测，若边缘定位精确轮廓清晰，最小二乘法可实现亚像素级别的精确拟合定位。

1.1公式推导过程

原理部分来自博客：Jacky_Ponder

2.基于opencv数据结构的c++代码实现

下面是一个最小二乘拟合圆心的函数：
参数输入1.存储2d点坐标的向量输出2.圆心坐标点输出double类型的圆的半径

void CircleFitting(std::vector<cv::Point> &pts,        // 2d点坐标存储向量cv::Point2d& center,                //圆心坐标double& radius)                       //圆的半径
{center = cv::Point2d(0, 0);           //初始化圆心为(0,0)radius = 0.0;                     //半径为0
//  if (pts.size() < 3) return false;    //判断输入点的个数若小于三个则不能拟合//定义计算中间变量double sumX1 = 0.0; //代表Xi的和(从1~n) ，X1代表X的1次方double sumY1 = 0.0;double sumX2 = 0.0; //代表(Xi)^2的和(i从1~n)，X2代表X的二次方double sumY2 = 0.0;double sumX3 = 0.0;double sumY3 = 0.0;double sumX1Y1 = 0.0;double sumX1Y2 = 0.0;double sumX2Y1 = 0.0;const double N = (double)pts.size();//获得输入点的个数for (int i = 0; i < pts.size(); ++i){double x = pts.at(i).x;            //获得第i个点的x坐标double y = pts.at(i).y;            //获得第i个点的y坐标double x2 = x * x;             //计算x^2double y2 = y * y;              //计算y^2double x3 = x2 * x;             //计算x^3double y3 = y2 * y;             //计算y^3double xy = x * y;              //计算xydouble x1y2 = x * y2;            //计算x*y^2double x2y1 = x2 * y;         //计算x^2*ysumX1 += x;                      //sumX=sumX+x;计算x坐标的和sumY1 += y;                        //sumY=sumY+y;计算y坐标的和sumX2 += x2;                   //计算x^2的和sumY2 += y2;                 //计算各个点的y^2的和sumX3 += x3;                 //计算x^3的和sumY3 += y3;sumX1Y1 += xy;sumX1Y2 += x1y2;sumX2Y1 += x2y1;}double C = N * sumX2 - sumX1 * sumX1;double D = N * sumX1Y1 - sumX1 * sumY1;double E = N * sumX3 + N * sumX1Y2 - (sumX2 + sumY2) * sumX1;double G = N * sumY2 - sumY1 * sumY1;double H = N * sumX2Y1 + N * sumY3 - (sumX2 + sumY2) * sumY1;double denominator = C * G - D * D;
//  if (std::abs(denominator) < DBL_EPSILON) return false;//判断分母的绝对值是否接近于（等于）0，使用此行判断须将函数改为bool类型，才能有返回值double a = (H * D - E * G) / (denominator);//denominator = D * D - G * C;//if (std::abs(denominator) < DBL_EPSILON) return false;double b = (H * C - E * D) / (-denominator);double c = -(a * sumX1 + b * sumY1 + sumX2 + sumY2) / N;center.x = a / (-2);center.y = b / (-2);radius = std::sqrt(a * a + b * b - 4 * c) / 2;
//  return true;将函数改为bool类型需要加返回值
}

以下为自己修改的程序，主要功能如下：
测试图片在下方
1.输入一张图片进行高斯滤波和边缘检测检测出圆形led的边缘
2.使用opencv提供的findContours()函数将提取的边缘存储于Circle_Data中
3.调用CircleFitting()函数进行最小二乘圆拟合，计算的结果存放于向量Circle_Center 和向量Circle_R中

#include<iostream>
#include<opencv2/opencv.hpp>//****函数声明圆心拟合函数
void CircleFitting(std::vector<std::vector<cv::Point>>& pts,std::vector<cv::Point2d>& center, std::vector<double>& radius);
//****函数声明endint main(int argc, char**argv) {//【0】输入图像并进行高斯滤波和canny边缘检测cv::Mat srcImg = cv::imread("Led_Point.jpg", 1);//以原图的形式输入RGBif (!srcImg.data) { std::cout << "enter srcImg wrong" << std::endl;    return false; }cv::Mat srcImg_1 = srcImg.clone();  //复制到srcImg_1中cv::Mat gray, dstImg;     //原图的灰度图，canny检测后输出cv::cvtColor(srcImg_1, gray, cv::COLOR_BGR2GRAY);//gray为滤波之前的灰度图cv::GssianBlur(gray, gray, cv::Size(3, 3), 0, 0, cv::BORDER_DEFAULT);cv::Canny(gray, gray, 3, 9, 3);//边缘检测后的图存放于gray中,单通道//std::cout << "gray.channels=" << gray.channels() << std::endl;dstImg.create(srcImg.size(), srcImg.type());//创建与原图相同尺寸和类型的dstImgdstImg = cv::Scalar::all(0);srcImg_1.copyTo(dstImg, gray);//测试通道数用//std::cout << "dstImg.channels=" << dstImg.channels() << std::endl;//三通道cv::cvtColor(dstImg, dstImg, cv::COLOR_BGR2GRAY);//转化为单通道灰度图cv::threshold(dstImg, dstImg, 100, 255, cv::THRESH_BINARY);//转化为二值图（单通道）//std::cout << "dstImg=" << dstImg.channels() << std::endl;//不使用cvtColor则为三通道//[1]使用opencv提供的findContours()函数将检测到的边缘按每个边缘一组存放于Circle_Data中std::vector<std::vector<cv::Point>> Circle_Data;//Circle_Data用于存放多个led边缘的坐标std::vector<cv::Vec4i> hierarchy;//函数findContours参数:dstImg必须为单通道图像cv::findContours(dstImg, Circle_Data, hierarchy, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE, cv::Point());//测试函数findContours输出是否正确//for (int i = 0; i < Circle_Data.size(); i++) {//  //Circle_Data[i]代表的是第i个轮廓，Circle_Data[i].size()代表的是第i个轮廓上所有的像素点数  // for (int j = 0; j < Circle_Data[i].size(); j++) {//       //每个点坐标为： cv::point(contours[i][j].x, contours[i][j].y);//       std::cout << "【" << i << "】" << "(" << Circle_Data[i][j].x << "," << Circle_Data[i][j].y << ")" << std::endl;// }//}//[3]创建存储每个led轮廓的圆心坐标向量和半径存储向量进行最小二乘曲线拟合拟合圆心坐标和圆的半径std::vector<cv::Point2d> Circle_Center;    //声明用于存储圆心坐标的向量std::vector<double> Circle_R;          //声明用于存储半径的向量//函数可处理多个圆轮廓的圆心和半径的拟合CircleFitting(Circle_Data, Circle_Center, Circle_R);//测试用：输出圆的坐标和半径for (int i = 0; i < Circle_Data.size(); i++) {std::cout << "Circle_Center[" << i << "]= (" << Circle_Center[i].x << "," << Circle_Center[i].y << ")" << std::endl;std::cout << "Circle_R[" << i << "]=" << Circle_R[i] << std::endl;}system("PAUSE");//用于保持输出终端return 0;
}//最小二乘拟合圆心原理https://blog.csdn.net/Jacky_Ponder/article/details/70314919
void CircleFitting(std::vector<std::vector<cv::Point>> &pts,        // 按组存储2d点坐标std::vector<cv::Point2d>& center,             //圆心坐标std::vector<double>& radius)                        //圆的半径
{for (int k = 0; k < pts.size(); k++) {// std::cout << "k=" << pts.size() << std::endl;center.push_back(cv::Point2d(0, 0));          //初始化圆心为(0,0)radius.push_back(0.0);                     //半径为0//    if (pts.size() < 3) return false;    //判断输入点的个数若小于三个则不能拟合//定义计算中间变量double sumX1 = 0.0; //代表Xi的和(从1~n) ，X1代表X的1次方double sumY1 = 0.0;double sumX2 = 0.0; //代表(Xi)^2的和(i从1~n)，X2代表X的二次方double sumY2 = 0.0;double sumX3 = 0.0;double sumY3 = 0.0;double sumX1Y1 = 0.0;double sumX1Y2 = 0.0;double sumX2Y1 = 0.0;const double N = (double)pts[k].size();//获得第k组输入点的个数for (int i = 0; i < pts[k].size(); ++i)//遍历第k组中所有数据{double x = 0;double y = 0;x = pts[k][i].x;          //获得第k组中第i个点的x坐标y = pts[k][i].y;           //获得第k组中第i个点的y坐标double x2 = x * x;             //计算x^2double y2 = y * y;              //计算y^2double x3 = x2 * x;             //计算x^3double y3 = y2 * y;             //计算y^3double xy = x * y;              //计算xydouble x1y2 = x * y2;            //计算x*y^2double x2y1 = x2 * y;         //计算x^2*ysumX1 += x;                      //sumX=sumX+x;计算x坐标的和sumY1 += y;                        //sumY=sumY+y;计算y坐标的和sumX2 += x2;                   //计算x^2的和sumY2 += y2;                 //计算各个点的y^2的和sumX3 += x3;                 //计算x^3的和sumY3 += y3;sumX1Y1 += xy;sumX1Y2 += x1y2;sumX2Y1 += x2y1;}double C = N * sumX2 - sumX1 * sumX1;double D = N * sumX1Y1 - sumX1 * sumY1;double E = N * sumX3 + N * sumX1Y2 - (sumX2 + sumY2) * sumX1;double G = N * sumY2 - sumY1 * sumY1;double H = N * sumX2Y1 + N * sumY3 - (sumX2 + sumY2) * sumY1;double denominator = C * G - D * D;//  if (std::abs(denominator) < DBL_EPSILON) return false;//判断分母的绝对值是否接近于（等于）0double a = (H * D - E * G) / (denominator);//denominator = D * D - G * C;//if (std::abs(denominator) < DBL_EPSILON) return false;double b = (H * C - E * D) / (-denominator);double c = -(a * sumX1 + b * sumY1 + sumX2 + sumY2) / N;center[k].x = a / (-2);center[k].y = b / (-2);radius[k] = std::sqrt(a * a + b * b - 4 * c) / 2;// return true;}
}

测试程序用图：
程序改进：将以上功能封装进一个函数中

#include<iostream>
#include<opencv2/opencv.hpp>void LedCircleCenter(cv::Mat &srcImg, std::vector<std::vector<cv::Point>> &Circle_Data,std::vector<cv::Vec4i> &hierarchy, std::vector<cv::Point2d> &Circle_Center,std::vector<double> &Circle_R);int main() {cv::Mat srcImg = cv::imread("Led_Point.jpg", 1);std::vector<std::vector<cv::Point>> CircleData;std::vector<cv::Vec4i> hierarchy;std::vector<cv::Point2d> CircleCenter;std::vector<double> CircleR;LedCircleCenter(srcImg,CircleData,hierarchy,CircleCenter,CircleR);for (int i = 0; i < CircleCenter.size(); i++) {std::cout<<"CircleCenter["<<i<<"]="<<CircleCenter[i]<<"  "<<"CircleR["<<i<<"]="<<CircleR[i]<<std::endl;}system("Pause");return 0;}//void LedCircleCenter(cv::Mat &srcImg, std::vector<std::vector<cv::Point>> &Circle_Data,
//  std::vector<cv::Vec4i> &hierarchy, std::vector<cv::Point2d> &Circle_Center,
//  std::vector<double> &Circle_R)
//function:函数用于将输入图片中的led灯的边缘提取后进行最小二乘拟合圆心和半径
//step:1.高斯滤波 2.canny边缘提取 3.存储每个边缘的坐标点 4.采用最小二乘法拟合圆心
//prgm 1: cv::Mat 类型的输入图像srcImg（格式不限）
/*prgm 2: std::vector<std::vector<cv::Point>> 类型的Circle_Data，存储了n组圆的边缘坐标，第一维代表圆的组号i,
第二维代表第i组中的j个坐标点,第三维代表第i个圆中的第j个点的x坐标和y坐标*/
//prgm 3:std::vector<cv::Vec4i> 检测用hierarchy(必须声明)
//prgm 4:std::vector<cv::Point2d> 类型的Circle_Center，用于存放i组圆圆心坐标的向量
//prgm 5:std::vector<double> 类型的Circle_R，用于存放i组圆对应的半径
void LedCircleCenter(cv::Mat &srcImg, std::vector<std::vector<cv::Point>> &Circle_Data,std::vector<cv::Vec4i> &hierarchy, std::vector<cv::Point2d> &Circle_Center,std::vector<double> &Circle_R) {cv::Mat srcImg_1 = srcImg.clone(); //复制到srcImg_1中cv::Mat gray, dstImg;     //原图的灰度图，canny检测后输出cv::cvtColor(srcImg_1, gray, cv::COLOR_BGR2GRAY);//gray为滤波之前的灰度图cv::GaussianBlur(gray, gray, cv::Size(5, 5), 0, 0, cv::BORDER_DEFAULT);//高斯滤波cv::Canny(gray, gray, 40, 120, 3);//边缘检测后的图存放于gray中,单通道，数据来源于dstImg.create(srcImg.size(), srcImg.type());//创建与原图相同尺寸和类型的dstImgdstImg = cv::Scalar::all(0);srcImg_1.copyTo(dstImg, gray);cv::cvtColor(dstImg, dstImg, cv::COLOR_BGR2GRAY);//转化为单通道灰度图cv::threshold(dstImg, dstImg, 100, 255, cv::THRESH_BINARY);//转化为二值图（单通道）//函数findContours,将检测到的边缘存放于Circle_Data中 参数:dstImg必须为单通道图像cv::findContours(dstImg, Circle_Data, hierarchy, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE, cv::Point());for (int k = 0; k < Circle_Data.size(); k++) {//    std::cout << "k=" << pts.size() << std::endl;Circle_Center.push_back(cv::Point2d(0, 0));           //初始化圆心为(0,0)Circle_R.push_back(0.0);                       //半径为0//    if (pts.size() < 3) return false;    //判断输入点的个数若小于三个则不能拟合//定义计算中间变量double sumX1 = 0.0; //代表Xi的和(从1~n) ，X1代表X的1次方double sumY1 = 0.0;double sumX2 = 0.0; //代表(Xi)^2的和(i从1~n)，X2代表X的二次方double sumY2 = 0.0;double sumX3 = 0.0;double sumY3 = 0.0;double sumX1Y1 = 0.0;double sumX1Y2 = 0.0;double sumX2Y1 = 0.0;const double N = (double)Circle_Data[k].size();//获得第k组输入点的个数for (int i = 0; i < Circle_Data[k].size(); ++i)//遍历第k组中所有数据{double x = 0;double y = 0;x = Circle_Data[k][i].x;          //获得第k组中第i个点的x坐标y = Circle_Data[k][i].y;           //获得第k组中第i个点的y坐标double x2 = x * x;             //计算x^2double y2 = y * y;              //计算y^2double x3 = x2 * x;             //计算x^3double y3 = y2 * y;             //计算y^3double xy = x * y;              //计算xydouble x1y2 = x * y2;            //计算x*y^2double x2y1 = x2 * y;         //计算x^2*ysumX1 += x;                      //sumX=sumX+x;计算x坐标的和sumY1 += y;                        //sumY=sumY+y;计算y坐标的和sumX2 += x2;                   //计算x^2的和sumY2 += y2;                 //计算各个点的y^2的和sumX3 += x3;                 //计算x^3的和sumY3 += y3;sumX1Y1 += xy;sumX1Y2 += x1y2;sumX2Y1 += x2y1;}double C = N * sumX2 - sumX1 * sumX1;double D = N * sumX1Y1 - sumX1 * sumY1;double E = N * sumX3 + N * sumX1Y2 - (sumX2 + sumY2) * sumX1;double G = N * sumY2 - sumY1 * sumY1;double H = N * sumX2Y1 + N * sumY3 - (sumX2 + sumY2) * sumY1;double denominator = C * G - D * D;//if (std::abs(denominator) < DBL_EPSILON) return false;//判断分母的绝对值是否接近于（等于）0double a = (H * D - E * G) / (denominator);//denominator = D * D - G * C;//if (std::abs(denominator) < DBL_EPSILON) return false;double b = (H * C - E * D) / (-denominator);double c = -(a * sumX1 + b * sumY1 + sumX2 + sumY2) / N;Circle_Center[k].x = a / (-2);Circle_Center[k].y = b / (-2);Circle_R[k] = std::sqrt(a * a + b * b - 4 * c) / 2;}
}