球拟合算法

一、 介绍

1. 本章节中的球拟合算法基于GSL、OCC等第三方库完成，代码已通过PTB认证测试：

2. 使用了GSL库中提供的单纯形法进行最小二乘距离拟合

二、 主要函数介绍

代价函数：

代价函数很简单，计算点集到球面的平方距离之和就行。

static double  CalcSphereCost(const vector<gp_Pnt>& points, const gp_Pnt& center, double radius)
{
	double cost = 0.0;
	double currentCost = 0.0;
	for (const auto& p : points) {
		currentCost = p.Distance(center) - radius;
		cost += currentCost * currentCost;
	}

	return cost;
}

迭代函数：

参数介绍：

gsl_vector 中包含了4个变量：dx dy dz dr

dx dy dz: 本此迭代 球的中心点X Y Z坐标值相对于初始值的行进距离

dr : 本此迭代 球的半径相对于初始值的行进距离

迭代过程：

1. 获取此次迭代的 dx dy dz dr值
2. 计算当前的球中心点坐标和半径
3. 调用代价函数计算当前球的代价

迭代目的

使得代价函数计算出来的代价最小

函数代码

void IterUpdate(const gsl_vector* v) {
	// 获取实参；
	double dx = gsl_vector_get(v, 0);
	double dy = gsl_vector_get(v, 1);
	double dz = gsl_vector_get(v, 2);
	double dr = gsl_vector_get(v, 3);

	// 当前球心；
	this->currentCenter = this->initCenter.XYZ() + gp_XYZ(dx,dy,dz);

	// 当前半径；
	this->currentRadius = this->initRadius + dr;

	// 当前代价；
	this->currentCost = CalcSphereCost(points, this->currentCenter, this->currentRadius);
}

单纯形法的调用函数

GSL单纯形法的使用参考手册：GNU Scientific Library — GSL 2.8 documentation

函数代码

	static const size_t ARGUE_NUM = 4;


	// 构建GSL params；
	SphereFitterParams theParams(measuredPoints);

	// 代价函数；
	gsl_multimin_function f;
	f.n = ARGUE_NUM;
	f.f = CostFunction;
	f.params = &theParams;

	// 初值；
	gsl_vector* x = gsl_vector_alloc(ARGUE_NUM);
	gsl_vector_set_all(x, 0.0);

	// 设置初始步长；
	gsl_vector* ss = gsl_vector_alloc(ARGUE_NUM);
	gsl_vector_set_all(ss, this->iterationStepSize);

	// 算法；
	const gsl_multimin_fminimizer_type* T = gsl_multimin_fminimizer_nmsimplex;

	// 创建&设置最小化器；
	gsl_multimin_fminimizer* s = gsl_multimin_fminimizer_alloc(T, ARGUE_NUM);
	gsl_multimin_fminimizer_set(s, &f, x, ss);

	// 开始迭代；
	int status = GSL_SUCCESS;
	for (size_t iterationNum = 0; iterationNum < this->iterationMax; iterationNum++)
	{
		// 迭代；
		status = gsl_multimin_fminimizer_iterate(s);

		// 检查迭代状态；
		if (status) {
			break;
		}

		// 评估迭代精度；
		double size = gsl_multimin_fminimizer_size(s);
		status = gsl_multimin_test_size(size, this->iterationPrecision);
		if (status == GSL_SUCCESS) {
			break;
		}
	}

	// 释放GSL解构；
	gsl_multimin_fminimizer_free(s);
	gsl_vector_free(ss);
	gsl_vector_free(x);

	// 迭代球；
	auto sphere = theParams.CurrentSphere();
	if (sphere == nullptr) {
		return make_tuple(gp_Sphere(), 0.0, false);
	}
	else {
		return make_tuple(*sphere, theParams.currentCost, true);
	}

其他

计算球的初始值：

// 计算点集的中心点；
this->initCenter = FitFacility::CalcCenter(points);

// 计算点集到中心点的平均距离.做初始半径
this->initRadius = FitFacility::CalcAverageDistance(points, this->initCenter);