散点拟合圆——最小二乘法

本文转载自这篇博客，修改了部分代码，并增加了算法展示效果
https://www.cnblogs.com/xiaxuexiaoab/p/16276402.html

最小二乘拟合原理

拟合效果展示
这里通过对原图中的圆形进行拟合，计算每个圆形的半径平均值。
原图

拟合结果图（原图放大2倍）

代码

• Fitter.h

class Fitter
{
public:
	static double FitCircle(const std::vector<Point2>& pointArray, Point2 & center, double& radius);
};

• Fitter.cpp

#include <limits>
#include <cmath>

#include "Fitter.h"

double Fitter::FitCircleByLeastSquare(const std::vector<Point2>& pointArray, Point2& center, double& radius)
{
	int N = pointArray.size();
	if (N < 3) {
		return std::numeric_limits<double>::max();
	}

	double sumX = 0.0; 
	double sumY = 0.0;
	double sumX2 = 0.0;
	double sumY2 = 0.0;
	double sumX3 = 0.0;
	double sumY3 = 0.0;
	double sumXY = 0.0;
	double sumXY2 = 0.0;
	double sumX2Y = 0.0;

	for (int pId = 0; pId < N; ++pId) {
		sumX += pointArray[pId].x;
		sumY += pointArray[pId].y;

		double x2 = pointArray[pId].x * pointArray[pId].x;
		double y2 = pointArray[pId].y * pointArray[pId].y;
		sumX2 += x2;
		sumY2 += y2;

		sumX3 += x2 * pointArray[pId].x;
		sumY3 += y2 * pointArray[pId].y;
		sumXY += pointArray[pId].x * pointArray[pId].y;
		sumXY2 += pointArray[pId].x * y2;
		sumX2Y += x2 * pointArray[pId].y;
 	}

	double C, D, E, G, H;
	double a, b, c;

	C = N * sumX2 - sumX * sumX;
	D = N * sumXY - sumX * sumY;
	E = N * sumX3 + N * sumXY2 - (sumX2 + sumY2) * sumX;
	G = N * sumY2 - sumY * sumY;
	H = N * sumX2Y + N * sumY3 - (sumX2 + sumY2) * sumY;

	a = (H * D - E * G) / (C * G - D * D);
	b = (H * C - E * D) / (D * D - G * C);
	c = -(a * sumX + b * sumY + sumX2 + sumY2) / N;

	center.x = -a / 2.0;
	center.y = -b / 2.0;
	radius = sqrt(a * a + b * b - 4 * c) / 2.0;

	double err = 0.0;
	double e;
	double r2 = radius * radius;
	for (int pId = 0; pId < N; ++pId){
		e = pow((pointArray[pId].x - center.x), 2) +  pow((pointArray[pId].y - center.y), 2) - r2;
		if (e > err) {
			err = e;
		}
	}
	return err;
}

• main.cpp

#include <iostream>
#include <vector>
#include "Fitter.h"

using namespace cv;
using namespace std;

int main() {

    Mat srcImg = imread("./9.bmp", 1);
    
    //将图像放大，减小边缘计算的误差
    Mat enlargeImg;
    int mulNum = 2;    //放大倍数
    resize(srcImg, enlargeImg, srcImg.size() * mulNum);

    Mat grayImg;
    cvtColor(enlargeImg, grayImg, COLOR_BGR2GRAY);
    Mat binaryImg;
    threshold(grayImg, binaryImg, 0, 255, THRESH_BINARY | THRESH_OTSU);

    vector<vector<cv::Point>> contours;
    vector<Vec4i> hierarchy;
    Mat edgeImg;
    Canny(binaryImg, edgeImg, 30, 100);
    findContours(edgeImg, contours, hierarchy, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE, cv::Point());
    
    Mat matDraw = Mat::zeros(binaryImg.size(), CV_8UC1);
    cv::Point center1;
    double r1, err1;
    double meanVal = 0.0;
    int countNum = 0;
    for (int i = 0; i < contours.size(); i++) {
        //进行最小二乘拟合
        err1 = Fitter().FitCircleByLeastSquare(contours[i], center1, r1);
        //若拟合的误差和半径在设定阈值范围内，认为这个圆拟合成功
        if (err1 < 3000 && r1 > 70 && r1 < 90) {
            circle(enlargeImg, center, r1, Scalar(255), 1);
            meanVal += r1;
            countNum += 1;
        }
        else {
            cout << i << " : " << err1 << endl;
        }
    }
    meanVal /= (countNum * mulNum);     //计算所有圆形的平均半径

    return 0;
}