利用多个二维点拟合圆、多个三维点拟合球面

C++代码，用最小二乘法拟合圆、拟合球面方程。之前在网上找拟合球面的代码，没找到C++代码，就自己写了个拟合球面的

#include
#include
using namespace std;

struct Point3f
{
	double x;
	double y;
	double z;
};

struct Point2f
{
	double x;
	double y;
};

bool ballFit（const vector &points double ¢er_x double ¢er_y double ¢er_z double &radius）
{
	center_x = 0.0f;
	center_y = 0.0f;
	center_z = 0.0f;
	radius = 0.0f;
	if （points.size（） < 3）
	{
		return false;
	}

	double sum_x = 0.0f sum_y = 0.0f sum_z = 0.0f;
	double sum_x2 = 0.0f sum_y2 = 0.0f sum_z2 = 0.0f;
	double sum_x3 = 0.0f sum_y3 = 0.0f sum_z3 = 0.0f;
	double sum_xy = 0.0f sum_x1y2 = 0.0f sum_x2y1 = 0.0f;
	double sum_xz = 0.0f sum_x1z2 = 0.0f sum_x2z1 = 0.0f;
	double sum_yz = 0.0f sum_y1z2 = 0.0f sum_y2z1 = 0.0f;

	int N = points.size（）;
	for （int i = 0; i < N; i++）
	{
		double x = points[i].x;
		double y = points[i].y;
		double z = points[i].z;
		double x2 = x * x;
		double y2 = y * y;
		double z2 = z * z;
		sum_x += x;
		sum_y += y;
		sum_z += z;

		sum_x2 += x2;
		sum_y2 += y2;
		sum_z2 += z2;

		sum_x3 += x2 * x;
		sum_y3 += y2 * y;
		sum_z3 += z2 * z;

		sum_xy += x * y;
		sum_x1y2 += x * y2;
		sum_x2y1 += x2 * y;

		sum_xz += x * z;
		sum_x1z2 += x * z2;
		sum_x2z1 += x2 * z;

		sum_yz += y * z;
		sum_y1z2 += y * z2;
		sum_y2z1 += y2 * z;

	}

	double C D E P G F H J K;
	double a b c d;

	C = N * sum_x2 - sum_x * sum_x;
	D = N * sum_xy - sum_x * sum_y;
	P = N * sum_xz - sum_x * sum_z;
	E = N * sum_x3 + N * sum_x1y2 + N * sum_x1z2 - （sum_x2 + sum_y2 + sum_z2） * sum_x;

	G = N * sum_y2 - sum_y * sum_y;
	F = N * sum_yz - sum_y * sum_z;
	H = N * sum_x2y1 + N * sum_y3 + N * sum_y1z2 - （sum_x2 + sum_y2 + sum_z2） * sum_y;

	J = N * sum_z2 - sum_z * sum_z;
	K = N * sum_x2z1 + N * sum_y2z1 + N * sum_z3 - （sum_x2 + sum_y2 + sum_z2） * sum_z;

	double CC DD EE GG HH;
	CC = C*J - P*P;
	DD = D*J - F*P;
	EE = E*J - P*K;
	GG = G*J - F*F;
	HH = H*J - F*K;

	a = （HH * DD - EE * GG） / （CC * GG - DD * DD）;
	b = （HH * CC - EE * DD） / （DD * DD - GG * CC）;
	c = 0.0f;
	if （P != 0）
	{
		c = （-C*a - D*b - E） / P;
	}
	else if （F != 0）
	{
		c = （-D*a - G*b - H） / F;
	}
	else if （J != 0）
	{
		c = （-P*a - F*b - K） / J;
	}
	d = -（a * sum_x + b * sum_y + c * sum_z + sum_x2 + sum_y2 + sum_z2） / N;

	center_x = a / （-2）;
	center_y = b / （-2）;
	center_z = c / （-2）;

	radius = sqrt（a * a + b * b + c * c - 4 * d） / 2;
	return true;
}

int main（）
{
	vector point_3d;
	Point3f xyz;

	for （float i = 0.5;i < 5.5;i += 0.1）
	{
		float j_nub = sqrt（25 - （i - 0.5）*（i - 0.5））;
		for （float j = 1.5;j < j_nub;j += 0.1）
		{
			xyz.x = i;
			xyz.y = j;
			float z = sqrt（25 - （i - 0.5）*（i - 0.5） - （j - 1.5）*（j - 1.5）） + 2;
			xyz.z = z;
			point_3d.push_back（xyz）;
		}
	}

	double x_3d y_3d z_3d r_3d;
	ballFit（point_3d x_3d y_3d z_3d r_3d）;
	cout << x_3d << “  x  “ << y_3d << “   y    “ << z_3d << “   z   “ << r_3d << endl;
	system（“p

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3022  2017-12-21 11:31  test.cpp

----------- ---------  ---------- -----  ----

                 3022                    1