最小二乘圆心拟合--11-27ECYXNH（CSDN）

 能跑

无图，只会算点

#include "iostream"
#include "opencv2\opencv.hpp"


using namespace std;
using namespace cv;

//****函数声明圆心拟合函数
void CircleFitting(
	std::vector<std::vector<cv::Point>>& pts,
	std::vector<cv::Point2d>& center, std::vector<double>& radius);
//****函数声明end

int main(int argc, char** argv) {

	//【0】输入图像并进行高斯滤波和canny边缘检测
	cv::Mat srcImg = cv::imread("E:\\img\\img1129-01.bmp", 1);//以原图的形式输入RGB
	if (!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::GaussianBlur(gray, gray, cv::Size(3, 3), 0, 0, cv::BORDER_DEFAULT);

	cv::Canny(gray, gray, 3, 9, 3);//边缘检测后的图存放于gray中,单通道
	imshow("canny", gray);
	//std::cout << "gray.channels=" << gray.channels() << std::endl;
	dstImg.create(srcImg.size(), srcImg.type());//创建与原图相同尺寸和类型的dstImg
	dstImg = 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^2
			double y2 = y * y;				//计算y^2
			double x3 = x2 * x;				//计算x^3
			double y3 = y2 * y;				//计算y^3
			double xy = x * y;				//计算xy
			double x1y2 = x * y2;			//计算x*y^2
			double x2y1 = x2 * y;			//计算x^2*y

			sumX1 += 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
		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[k].x = a / (-2);
		center[k].y = b / (-2);
		radius[k] = std::sqrt(a * a + b * b - 4 * c) / 2;
		//	return true;
	}
}