精华内容
下载资源
问答
  • 2019-08-13 21:42:25
    1. NuGet->安装MathNet.Numerics ;
    2. 引用“System.Drawing”;
    3. Using指令:
    using System.Drawing;
    using MathNet.Numerics.LinearAlgebra;
    
    1. 代码块(转)
    /// <summary>
    /// 最小二乘法拟合圆获取圆心坐标、圆半径
    /// </summary>
    /// <param name="pointFs">输入的List<PointF>类型的变量</param>
    /// <param name="CenterX">圆心x坐标</param>
    /// <param name="CenterY">圆心y坐标</param>
    /// <param name="CenterR">圆半径</param>
    /// <returns></returns>
    public static int FitCircle(List<PointF> pointFs, out float CenterX, out float CenterY, out float CenterR)
    {
        Matrix<float> YMat;
        Matrix<float> RMat;
        Matrix<float> AMat;
        List<float> YLit = new List<float>();
        List<float[]> RLit = new List<float[]>();
    
        foreach (var pointF in pointFs)
            YLit.Add(pointF.X * pointF.X + pointF.Y * pointF.Y);
    
        float[,] Yarray = new float[YLit.Count, 1];
        for (int i = 0; i < YLit.Count; i++)
            Yarray[i, 0] = YLit[i];
    
        YMat = CreateMatrix.DenseOfArray<float>(Yarray);
    
        foreach (var pointF in pointFs)
            RLit.Add(new float[] { -pointF.X, -pointF.Y, -1 });
    
        float[,] Rarray = new float[RLit.Count, 3];
        for (int i = 0; i < RLit.Count; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                Rarray[i, j] = RLit[i][j];
            }
        }
    
        RMat = CreateMatrix.DenseOfArray<float>(Rarray);
        Matrix<float> RTMat = RMat.Transpose();
        Matrix<float> RRTInvMat = (RTMat.Multiply(RMat)).Inverse();
        AMat = RRTInvMat.Multiply(RTMat.Multiply(YMat));
    
        float[,] Array = AMat.ToArray();
        float A = Array[0, 0];
        float B = Array[1, 0];
        float C = Array[2, 0];
        CenterX = A / -2.0f;
        CenterY = B / -2.0f;
        CenterR = (float)(Math.Sqrt(A * A + B * B - 4 * C) / 2.0f);
        return 0x0000;
    }
    
    更多相关内容
  • 最小二乘法拟合圆心,基于Hough变换的圆心检测,基于harris亚像素棋盘格检测,对三种方法角点检测进度进行对比分析。
  • 最小二乘法拟合圆心和半径 ,转化为所有点到圆心(xc, yc)距离与半径Rc的差最小问题,表达式如下,再通过scipy中的最小二乘算法进行求解即可。 因为圆半径等于所有点到圆心的平均距离,所以一旦圆心得到了,就可以...

    代码实现

    scipy.optimize.leastsq

    最小二乘法拟合圆心和半径 ,转化为所有点到圆心(xc, yc)距离与半径Rc的差最小问题,表达式如下,再通过scipy中的最小二乘算法进行求解即可。
    因为圆半径等于所有点到圆心的平均距离,所以一旦圆心得到了,就可以直接求得半径

    Ri = sqrt( (x - xc)**2 + (y - yc)**2)
    residu = sum( (Ri - Rc)**2)
    

    展示一个例子

    from numpy import *
    
    # Coordinates of the 2D points
    
    x = r_[  9, 35, -13,  10,  23,   0]
    y = r_[ 34, 10,   6, -14,  27, -10]
    basename = 'circle'
    x_m = mean(x)
    y_m = mean(y)
    
    # Decorator to count functions calls
    import functools
    def countcalls(fn):
        "decorator function count function calls "
    
        @functools.wraps(fn)
        def wrapped(*args):
            wrapped.ncalls +=1
            return fn(*args)
    
        wrapped.ncalls = 0
        return wrapped
    
    #  == METHOD 2 ==
    # Basic usage of optimize.leastsq
    from scipy      import optimize
    
    
    method_2  = "leastsq"
    
    def calc_R(xc, yc):
        """ calculate the distance of each 2D points from the center (xc, yc) """
        return sqrt((x-xc)**2 + (y-yc)**2)
    
    @countcalls
    def f_2(c):
        """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
        Ri = calc_R(*c)
        return Ri - Ri.mean()
    
    center_estimate = x_m, y_m
    center_2, ier = optimize.leastsq(f_2, center_estimate)
    
    xc_2, yc_2 = center_2
    Ri_2       = calc_R(xc_2, yc_2)
    R_2        = Ri_2.mean()
    residu_2   = sum((Ri_2 - R_2)**2)
    residu2_2  = sum((Ri_2**2-R_2**2)**2)
    ncalls_2   = f_2.ncalls
    
    # == METHOD 2b ==
    # Advanced usage, with jacobian
    method_2b  = "leastsq with jacobian"
    
    def calc_R(xc, yc):
        """ calculate the distance of each 2D points from the center c=(xc, yc) """
        return sqrt((x-xc)**2 + (y-yc)**2)
    
    @countcalls
    def f_2b(c):
        """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
        Ri = calc_R(*c)
        return Ri - Ri.mean()
    
    @countcalls
    def Df_2b(c):
        """ Jacobian of f_2b
        The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
        xc, yc     = c
        df2b_dc    = empty((len(c), x.size))
    
        Ri = calc_R(xc, yc)
        df2b_dc[ 0] = (xc - x)/Ri                   # dR/dxc
        df2b_dc[ 1] = (yc - y)/Ri                   # dR/dyc
        df2b_dc       = df2b_dc - df2b_dc.mean(axis=1)[:, newaxis]
    
        return df2b_dc
    
    center_estimate = x_m, y_m
    center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True)
    
    xc_2b, yc_2b = center_2b
    Ri_2b        = calc_R(xc_2b, yc_2b)
    R_2b         = Ri_2b.mean()
    residu_2b    = sum((Ri_2b - R_2b)**2)
    residu2_2b   = sum((Ri_2b**2-R_2b**2)**2)
    ncalls_2b    = f_2b.ncalls
    
    print ("""
    Method 2b :
    print "Functions calls : f_2b=%d Df_2b=%d""" % ( f_2b.ncalls, Df_2b.ncalls))
    
    print  (method_2 , xc_2 , yc_2 , R_2 , ncalls_2 , Ri_2.std() , residu_2 , residu2_2 )
    

    参考
    https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html

    展开全文
  • 最小二乘法拟合圆心 文章为个人学习过程中笔记,原理部分参考其他作者内容,侵权必删 最小二乘法(least squares analysis)是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配。最小二乘法是用...

    最小二乘法拟合圆心

    文章为个人学习过程中笔记,原理部分参考其他作者内容,侵权必删
    最小二乘法(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^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,使用此行判断须将函数改为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);
    //****函数声明end
    
    int main(int argc, char**argv) {
    
    	//【0】输入图像并进行高斯滤波和canny边缘检测
    	cv::Mat srcImg = cv::imread("Led_Point.jpg", 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::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());//创建与原图相同尺寸和类型的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;
    	}
    }
    
    

    测试程序用图:在这里插入图片描述
    程序改进:将以上功能封装进一个函数中

    #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());//创建与原图相同尺寸和类型的dstImg
    	dstImg = 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^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;
    
    		Circle_Center[k].x = a / (-2);
    		Circle_Center[k].y = b / (-2);
    		Circle_R[k] = std::sqrt(a * a + b * b - 4 * c) / 2;
    	}
    }
    
    展开全文
  • 使用最小二乘法对一个点集拟合一个圆,返回该拟合圆的圆心和半径。c++的opencv版本。 原理代码:https://blog.csdn.net/qq_42914355/article/details/104927023
  • 该方法采用提取图像轮廓,并进行最小二乘法计算拟合出圆。最后得出圆心和半径。
  • 最小二乘法拟合圆.zip

    2020-12-15 23:16:55
    matlab版图像读取, 颜色转换 ,灰度化, 二值化, 边缘检测, 进行最小二乘法拟合。将不标准圆进行拟合。获得最后的圆心坐标以及半径大小。
  • 最小二乘法拟合

    2015-09-06 18:22:08
    for(int i=0;i;i++){ int x=samples[i].x; int y=samples[i].y; X1=X1+x; Y1=Y1+y; X2=X2+x*x; Y2 = Y2 +y*y; X3 = X3 + x*x*x; Y3 = Y3 + y*y*y; X1Y1 = X1Y1 + x*y; X1Y2 = X1Y2 + x*y*y;...
  • 拟合的方法有很多种,最小二乘法属于比较简单的一种。今天就先将这种。我们知道圆方程可以写为:(x?xc)2+(y?yc)2=R2(x?xc)2+(y?yc)2=R2通常的最小二乘拟合要求距离的平方和最小。也就是f=∑((xi?xc)...

    有一系列的数据点 {xi,yi}{xi,yi},我们知道这些数据点近似的落在一个圆上,根据这些数据估计这个圆的参数就是一个很有意义的问题。今天就来讲讲如何来做圆的拟合。圆拟合的方法有很多种,最小二乘法属于比较简单的一种。今天就先将这种。

    我们知道圆方程可以写为:

    (x?xc)2+(y?yc)2=R2(x?xc)2+(y?yc)2=R2

    通常的最小二乘拟合要求距离的平方和最小。也就是

    f=∑((xi?xc)2+(yi?yc)2??????????????????√?R)2f=∑((xi?xc)2+(yi?yc)2?R)2

    最小。这个算起来会很麻烦。也得不到解析解。所以我们退而求其次。

    f=∑((xi?xc)2+(yi?yc)2?R2)2f=∑((xi?xc)2+(yi?yc)2?R2)2

    这个式子要简单的多。我们定义一个辅助函数:

    g(x,y)=(x?xc)2+(y?yc)2?R2g(x,y)=(x?xc)2+(y?yc)2?R2

    那么上面的式子可以表示为:

    f=∑g(xi,yi)2f=∑g(xi,yi)2

    按照最小二乘法的通常的步骤,可知ff 取极值时对应下面的条件。

    ?f?xc=0?f?yc=0?f?R=0?f?xc=0?f?yc=0?f?R=0

    先来化简 ?f?R=0?f?R=0

    ?f?R=?2R×∑((xi?xc)2+(yi?yc)2?R2)=?2R×∑g(xi,yi)=0?f?R=?2R×∑((xi?xc)2+(yi?yc)2?R2)=?2R×∑g(xi,yi)=0

    我们知道半径 RR 是不能为 0 的。所以必然有:

    ∑g(xi,yi)=0∑g(xi,yi)=0

    这是个非常有用的结论。剩下的两个式子:

    ?f?xc=?4∑((xi?xc)2+(yi?yc)2?R2)(xi?xc)=?4∑g(xi,yi)(xi?xc)=?4∑xig(xi,yi)=0?f?yc=?4∑((xi?xc)2+(yi?yc)2?R2)(yi?yc)=?4∑g(xi,yi)(yi?yc)=?4∑yig(xi,yi)=0?f?xc=?4∑((xi?xc)2+(yi?yc)2?R2)(xi?xc)=?4∑g(xi,yi)(xi?xc)=?4∑xig(xi,yi)=0?f?yc=?4∑((xi?xc)2+(yi?yc)2?R2)(yi?yc)=?4∑g(xi,yi)(yi?yc)=?4∑yig(xi,yi)=0

    也就是下面两个等式:

    ∑xig(xi,yi)=0∑yig(xi,yi)=0∑xig(xi,yi)=0∑yig(xi,yi)=0

    这里设:

    ui=xi?xˉuc=xc?xˉvi=yi?yˉvc=yc?yˉui=xi?xˉuc=xc?xˉvi=yi?yˉvc=yc?yˉ

    其中:

    xˉ=∑xi/Nyˉ=∑yi/Nxˉ=∑xi/Nyˉ=∑yi/N

    那么简单的推导一下,就会发现:

    ∑uig(ui,vi)=0∑vig(ui,vi)=0∑uig(ui,vi)=0∑vig(ui,vi)=0

    其实都不需要推导,这个变量替换只不过是修改了坐标原点的位置。而我们的等式根本就与坐标原点的具体位置无关。所以必然成立。

    这两个式子展开写是这样的:

    ∑((ui?uc)2+(vi?vc)2?R2)ui=0∑((ui?uc)2+(vi?vc)2?R2)vi=0∑((ui?uc)2+(vi?vc)2?R2)ui=0∑((ui?uc)2+(vi?vc)2?R2)vi=0

    进一步展开:

    ∑(u3i?2u2iuc+uiu2c+uiv2i?2uivivc+uiv2c?uiR2)=0∑(u2ivi?2uiviuc+viu2c+v3i?2v2ivc+viv2c?viR2)=0∑(ui3?2ui2uc+uiuc2+uivi2?2uivivc+uivc2?uiR2)=0∑(ui2vi?2uiviuc+viuc2+vi3?2vi2vc+vivc2?viR2)=0

    我们知道 ∑ui=0∑ui=0, ∑vi=0∑vi=0。所以上面两个式子是可以化简的。

    ∑(u3i?2u2iuc+uiv2i?2uivivc)=0∑(u2ivi?2uiviuc+v3i?2v2ivc)=0∑(ui3?2ui2uc+uivi2?2uivivc)=0∑(ui2vi?2uiviuc+vi3?2vi2vc)=0

    为了简化式子,我们定义几个参数:

    Suuu=∑u3iSvvv=∑v3iSuu=∑u2iSvv=∑v2iSuv=∑uiviSuuv=∑u2iviSuvv=∑uiv2iSuuu=∑ui3Svvv=∑vi3Suu=∑ui2Svv=∑vi2Suv=∑uiviSuuv=∑ui2viSuvv=∑uivi2

    那么上面的式子可以写为:

    Suuuc+Suvvc=Suuu+Suvv2Suvuc+Svvvc=Suuv+Svvv2Suuuc+Suvvc=Suuu+Suvv2Suvuc+Svvvc=Suuv+Svvv2

    至此,就可以解出ucuc 和vcvc 了。

    uc=SuuvSuv?SuuuSvv?SuvvSvv+SuvSvvv2(S2uv?SuuSvv)vc=?SuuSuuv+SuuuSuv+SuvSuvv?SuuSvvv2(S2uv?SuuSvv)uc=SuuvSuv?SuuuSvv?SuvvSvv+SuvSvvv2(Suv2?SuuSvv)vc=?SuuSuuv+SuuuSuv+SuvSuvv?SuuSvvv2(Suv2?SuuSvv)

    那么:

    xc=uc+xˉyc=vc+yˉxc=uc+xˉyc=vc+yˉ

    还剩下个 RR 没求出来, 也很简单。

    ∑g(xi,yi)=0∑((xi?xc)2+(yi?yc)2?R2)=0∑g(xi,yi)=0∑((xi?xc)2+(yi?yc)2?R2)=0

    所以:

    R2=∑((xi?xc)2+(yi?yc)2)R2=∑((xi?xc)2+(yi?yc)2)

    好了。下面给出个代码,这个代码的具体公式和我这里给出的有一点小差异,但是原理是相同的。

    /**

    * 最小二乘法拟合圆

    * 拟合出的圆以圆心坐标和半径的形式表示

    * 此代码改编自 newsmth.net 的 jingxing 在 Graphics 版贴出的代码。

    * 版权归 jingxing, 我只是搬运工外加一些简单的修改工作。

    */

    typedef complex POINT;

    bool circleLeastFit(const std::vector &points, double &center_x, double &center_y, double &radius)

    {

    center_x = 0.0f;

    center_y = 0.0f;

    radius = 0.0f;

    if (points.size() < 3)

    {

    return false;

    }

    double sum_x = 0.0f, sum_y = 0.0f;

    double sum_x2 = 0.0f, sum_y2 = 0.0f;

    double sum_x3 = 0.0f, sum_y3 = 0.0f;

    double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;

    int N = points.size();

    for (int i = 0; i < N; i++)

    {

    double x = points[i].real();

    double y = points[i].imag();

    double x2 = x * x;

    double y2 = y * y;

    sum_x += x;

    sum_y += y;

    sum_x2 += x2;

    sum_y2 += y2;

    sum_x3 += x2 * x;

    sum_y3 += y2 * y;

    sum_xy += x * y;

    sum_x1y2 += x * y2;

    sum_x2y1 += x2 * y;

    }

    double C, D, E, G, H;

    double a, b, c;

    C = N * sum_x2 - sum_x * sum_x;

    D = N * sum_xy - sum_x * sum_y;

    E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;

    G = N * sum_y2 - sum_y * sum_y;

    H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;

    a = (H * D - E * G) / (C * G - D * D);

    b = (H * C - E * D) / (D * D - G * C);

    c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;

    center_x = a / (-2);

    center_y = b / (-2);

    radius = sqrt(a * a + b * b - 4 * c) / 2;

    return true;

    }

    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

    下图是个实际测试的结果。对这种均匀分布的噪声数据计算的结果还是很准确的。

    但是当数据中有部分偏向某一个方向的干扰数据时。结果就不那么乐观了。下图就很说明问题。

    数据点中有 20% 是干扰数据。拟合出的圆就明显被拽偏了。

    之所以出现这个问题就是因为最小二乘拟合的平方项对离群点非常敏感。解决这个问题就要用其他的拟合算法,比如用距离之和作为拟合判据。等我有空了再写一篇博客介绍其他几种方法。

    展开全文
  • 我试图拟合一些数据点以找到圆的中心。以下所有点都是围绕圆周长的噪声数据点:data = [(2.2176383052987667, 4.218574252410221),(3.3041214516913033, 5.223500807396272),(4.280815855023374, 6.461487709813785)...
  • 最小二乘法拟合圆——MATLAB和Qt-C++实现

    万次阅读 多人点赞 2018-07-01 21:24:46
    本节Jungle尝试用最小二乘法拟合圆,并用MATLAB和C++实现。 1.理论知识 根据圆心(A,B)和半径R可确定平面上一个圆。平面上圆方程的通式为 其中 第一个圆的通式是关于a、b和c的线性方程。利用最小二乘法...
  • 最小二乘法拟合圆弧

    千次阅读 2022-04-05 17:55:31
    (4)最小二乘法拟合圆弧。 本文要拟合的边缘如实验图中绿色区域,根据检测出的坐标特点提取出所要拟合的边缘坐标。 实验图和结果图如下所示: MATLAB实现: clear; close all; clc; Origin=im2gray(imread('circle...
  • 最小二乘法拟合(附完整代码)

    千次阅读 2021-12-04 21:14:41
    一种简洁、高效的2D/3D圆弧拟合方法,圆弧拟合可以选择是否精确经过起点和终点。附完整代码,方便学习。
  • https://blog.csdn.net/Jacky_Ponder/article/details/703149191.1最小二乘拟合圆介绍与推导最小二乘法(least squares analysis)是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配。最小...
  • 前言 椭圆在高中数学里就开始提到,都是从标准方程开始如: x2a2+y2b2=1(a>b>0) \frac{x^2}{a^2}+\frac{y^2}{b^2}=1(a>b>0) a2x2​+b2y2​=1(a>b>0) 其中 aaa 为椭圆的长轴, bbb 为椭圆短轴,通常默认椭圆圆心 ...
  • 拟合的方法有非常多种,最小二乘法属于比較简单的一种。今天就先将这样的。我们知道圆方程能够写为:(x?xc)2+(y?yc)2=R2通常的最小二乘拟合要求距离的平方和最小。也就是f=∑((xi?xc)2+(yi?yc)2????????????...
  • C++:最小二乘法 拟合

    千次阅读 2019-12-30 15:39:51
    一、推导 ... //圆心XY 半径 double xCenter = -0.5*a; double yCenter = -0.5*b; radius = 0.5 * sqrt(a * a + b * b - 4 * c); centerP[0] = xCenter; centerP[1] = yCenter; }  
  • 最小二乘法拟合圆(Python)

    万次阅读 多人点赞 2019-03-26 17:16:01
    上文已经对比了三种数据点拟合圆的方法,本文分享最小二乘法拟合过程。旨在了解如何用Python编程拟合圆。 #! /usr/bin/env python # -*- coding: utf-8 -*- """ This program is debugged by Harden Qiu """ ...
  • 最小二乘法圆拟合(含matlab程序及说明).ppt* 最小二乘法 圆拟合 (含matlab程序及说明) 最小二乘法拟合圆曲线:令a=-2A,b=-2B,则:圆的另一形式为: * 只需求出参数a,b,c即可以求的圆半径参数: * 样本集点 到圆边缘...
  • 之前介绍了用矩阵形式最小二乘法拟合平平面的方式,这次来做用矩阵形式来拟合空间球(圆)(由于测量误差的存在,认为空间圆的测量值分布在一个球面上,且球心与空间圆心重合)。二维空间中求圆心和求球的方法相同,只是...
  • 学过统计学的朋友们应该知道最小二乘法至关重要的地位,目前我的研究方面要对数据进行拟合,具体任务是在平面坐标系中拟合圆心位置和半径来,首当其冲的想到了这一种方法,因此便把所学内容整理出来和大家分享。...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 721
精华内容 288
关键字:

最小二乘法拟合圆心

友情链接: ex04_adc.zip