基於OpenCV資料結構最小二乘法擬合圓-程式碼部分
阿新 • • 發佈：2018-12-12
對於網上常用的擬合圓程式碼(經過修改， 因為除數可能為0)
/* * 參考: http://blog.csdn.net/liyuanbhu/article/details/50889951 * 通過最小二乘法來擬合圓的資訊 * pts: 所有點座標 * center: 得到的圓心座標 * radius: 圓的半徑 */ static bool CircleInfo2(std::vector<cv::Point2d>& pts, cv::Point2d& center, double& radius) { center = cv::Point2d(0, 0); radius = 0.0; if (pts.size() < 3) return false;; 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 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; double y = pts.at(i).y; double x2 = x * x; double y2 = y * y; double x3 = x2 *x; double y3 = y2 *y; double xy = x * y; double x1y2 = x * y2; double x2y1 = x2 * y; sumX += x; sumY += y; sumX2 += x2; sumY2 += y2; sumX3 += x3; sumY3 += y3; sumXY += xy; sumX1Y2 += x1y2; sumX2Y1 += x2y1; } double C = N * sumX2 - sumX * sumX; double D = N * sumXY - sumX * sumY; double E = N * sumX3 + N * sumX1Y2 - (sumX2 + sumY2) * sumX; double G = N * sumY2 - sumY * sumY; double H = N * sumX2Y1 + N * sumY3 - (sumX2 + sumY2) * sumY; double denominator = C * G - D * D; if (std::abs(denominator) < DBL_EPSILON) return false; 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 * sumX + b * sumY + sumX2 + sumY2) / N; center.x = a / (-2); center.y = b / (-2); radius = std::sqrt(a * a + b * b - 4 * c) / 2; return true; } --------------------- 本文來自 wfh2015 的CSDN 部落格 ，全文地址請點選：https://blog.csdn.net/wfh2015/article/details/79626581?utm_source=copy
然而，很多人對這部分程式碼表示有疑問，所以我根據網上的資料另外寫了一個部分，但是速度差很多，僅僅用於學習程式碼公式是按照這個部落格來的。個人認為，這份部落格的大部分是正確的，比如求圓的圓心是正確的。但是最後面的計算圓的半徑部落格公式可能有問題，我根據下面網友評論修改了一下，得到正確的結果。下面是我寫的程式碼
#include <iostream> #include <string> #include "opencv2/opencv.hpp" /////////////////////////////////////根據公式編寫程式碼////////////////////////////////////////// /** 計算所有點的X座標均值 pts: 點座標 return: 均值 */ static double meanX(std::vector<cv::Point2d>& pts) { if (pts.size() <= 0) return 0.0; const double N = (double)pts.size(); double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { cv::Point2d pt = pts.at(i); sum += (pt.x); } double avg = sum / N; return avg; } /** 計算所有點的Y座標均值 pts: 點座標 return: 均值 */ static double meanY(std::vector<cv::Point2d>& pts) { if (pts.size() <= 0) return 0.0; const double N = (double)pts.size(); double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { cv::Point2d pt = pts.at(i); sum += (pt.y); } double avg = sum / N; return avg; } static double Ui(std::vector<cv::Point2d>& pts, int index, double meanXValue) { double xi = pts.at(index).x; return xi - meanXValue; } static double Vi(std::vector<cv::Point2d>& pts, int index, double meanYValue) { double yi = pts.at(index).y; return yi - meanYValue; } static double Suvv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double u = Ui(pts, i, mean_x); double v = Vi(pts, i, mean_y); double cur = u * v *v; sum += cur; } return sum; } static double Suuv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double u = Ui(pts, i, mean_x); double v = Vi(pts, i, mean_y); double cur = u * u * v; sum += cur; } return sum; } static double Suv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double u = Ui(pts, i, mean_x); double v = Vi(pts, i, mean_y); double cur = u * v; sum += cur; } return sum; } static double Svv(std::vector<cv::Point2d>& pts, double mean_y) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double value = Vi(pts, i, mean_y); double cur = value * value; sum += cur; } return sum; } static double Suu(std::vector<cv::Point2d>& pts, double mean_x) { double sum = 0; for (int i = 0; i < pts.size(); ++i) { double value = Ui(pts, i, mean_x); double cur = value * value; sum += cur; } return sum; } static double Svvv(std::vector<cv::Point2d>& pts, double mean_y) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double v = Vi(pts, i, mean_y); double cur = v * v * v; sum += cur; } return sum; } static double Suuu(std::vector<cv::Point2d>& pts, double mean_x) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double u = Ui(pts, i, mean_x); double cur = u * u * u; sum += cur; } return sum; } static double Uc(double A, double B, double C, double D, double E, double F, double G) { double numerator = A * B - C * D - E * D + B * F; double denominator = 2 * (B * B - G * D); double result = numerator / denominator; return result; } static double Vc(double A, double B, double C, double D, double E, double F, double G) { double numerator = -1 * G * A + C * B + B * E - G * F; double denominator = 2 * (B * B - G * D); double result = numerator / denominator; return result; } static double GetRadius(std::vector<cv::Point2d>& pts, cv::Point2d center) { double sum = 0.0; for (int i = 0; i < pts.size(); ++i) { double p1 = pts.at(i).x - center.x; double p2 = pts.at(i).y - center.y; double cur = p1*p1 + p2*p2; sum += cur; } double radius = std::sqrt(sum / (pts.size())); return radius; } static void CircleInfo(std::vector<cv::Point2d>& pts, cv::Point2d& center, double& radius) { double mean_x = meanX(pts); double mean_y = meanY(pts); double A = Suuv(pts, mean_x, mean_y); double B = Suv(pts, mean_x, mean_y); double C = Suuu(pts, mean_x); double D = Svv(pts, mean_y); double E = Suvv(pts, mean_x, mean_y); double F = Svvv(pts, mean_y); double G = Suu(pts, mean_x); double uc = Uc(A, B, C, D, E, F, G); double vc = Vc(A, B, C, D, E, F, G); center.x = uc + mean_x; center.y = vc + mean_y; radius = GetRadius(pts, center); } --------------------- 本文來自 wfh2015 的CSDN 部落格 ，全文地址請點選：https://blog.csdn.net/wfh2015/article/details/79626581?utm_source=copy
另外，附上所有的測試程式碼下載。 最後，感謝我參考的部落格作者的公式推導。謝謝！