B样条拟合

徐开元

开发环境

操作系统：windows 11

编程语言：C++ & CUDA

编译器：MSVC v142（Visual Studio 2019）& nvcc

CUDA Runtime版本：11.0

CPU：AMD Ryzen 3900X

内存： DDR4@3000Mhz 32GB

GPU：NVIDIA RTX 3080

使用说明

BSplineFitting是用B样条拟合封闭平面曲线的程序，输入平面上构成封闭曲线的数据点集和初始控制点集，输出拟合数据点的3次B样条控制点，并将曲线绘制出来。

可以使用数据点到曲线的代数距离或几何距离（一阶/二阶）作为优化的目标函数。代数距离为PDM，一阶几何距离为TDM，二阶几何距离为SDM。

本程序使用GPU做曲线投影点匹配和数值求解，需要CUDA运行时版本不低于11.0才能运行。

工程源文件在./src目录，可执行程序./bin目录，资源文件在./res目录。

运行程序：执行以下命令

./bin/BSplineFitting.exe config.json

其中config.json是json格式的配置文件，示例如下：

{
    "pts": "../res/tian.pts",   // 需要拟合的数据点
    "control_pts": "../res/tian.cpts",  // 初始控制点
    "sample_per_range": 1024,   // 投影点采样数，要求为2的n次方，不超过1024
    "max_iter_projection": 30,  // 计算数据点在曲线上投影的最大迭代次数
    "solver": "PDM",  // 优化目标函数，可选PDM，TDM或SDM
    "max_iter_solver":20,  // 曲线拟合的最大迭代次数
    "errot_target": 0.1   // RMSE标准误差的阈值，低于该值停止迭代
}

执行完成得到每轮迭代的运行报告，包括算出数据点在曲线上投影的迭代次数、RMSE误差。最后输出最终的控制点坐标并绘制曲线。以上述配置文件为例，运行得到如下曲线：

算法

B样条曲线

本项目选择使用unit knot vector的3次B样条曲线。代码位于./src/Geometry目录下以cubic_bspline_curve开头的文件中。

数据点与初始控制点

这一部分并不是算法的重点，本人是从网上找图片，使用python的图像处理库计算图案的轮廓点，并采样轮廓点，加上一些噪音，生成b样条控制点，使控制点数量大致为数据点的几十分之一。资源文件目录res中，包含多个数据点与控制点文件。数据点后缀为.pts，控制点后缀为.cpts。

投影计算

无论以代数距离还是几何距离作为目标函数，都需要计算每个数据点d_k在当前B样条曲线上的对应投影点p(t_k)（t_k为参数），使得向量p(t_k)d_k垂直于曲线在t_k处的切向量。

本人采用一种迭代式算法计算投影点。

为每个数据点d_k（0<=k<=nd-1）记录一个在b样条上的采样参数区间[tl[k],tr[k])，采样点个数为配置文件中指定的sample_per_range个，简记为spr。

1.采样阶段：根据d_k对应的采样区间，在b样条曲线上均匀采样spr个点，记录这些点的坐标、参数、切向量和法向量；

2.比较与挑选阶段：将每个数据点d_k和与它对应的spr个曲线上的采样点p(t_ki)(0<=i<=spr-1) 连成向量，判断向量p(t_ki)d_k和曲线在t_ki处切向量的夹角。若夹角<epsilon，则视为向量与曲线垂直，找到投影点，直接返回；否则，按夹角是锐角还是钝角将t_ki分为两类，每一类中挑选出参数tm，使曲线在tm处法向量与d_k的距离最小。从锐角类和钝角类中挑选出的两个tm构成了d_k的新的采样参数区间[tl[k]_new, tr[k]_new)。如下图所示。

3.若仍存在未找到曲线投影点的数据点，则重复1、2。本程序只考虑拟合平面曲线为封闭曲线的情形，因此随着每个数据点的采样参数区间越来越小，投影点最终一定会收敛到某个B样条曲线的点处。

代码（./src/Algorithm/CurveFitting/pt2curve_projecion.cu的computeProjectionDirect函数）：

while (num_not_converged > 0 && _iter < max_iter) {
	sampleCubicBSplineCurveLinesLocal << <num_pts, spr >> > (
		curve,
		spr, getRawPtr(found), getRawPtr(tl), getRawPtr(tr),
		tsamples, p_tsamples,
		tangents, normals
		);
	cudaDeviceSynchronize();

	pickClosestTksDirect << <num_pts, spr, sizeof(Float_t)* spr * 4 >> > (
		curve,
		getRawPtr(pts),
		tsamples, p_tsamples,
		tangents, normals,
		getRawPtr(tk), getRawPtr(pts_projected),
		getRawPtr(found), getRawPtr(tl), getRawPtr(tr)
		);
	cudaDeviceSynchronize();
	num_not_converged = thrust::count(found.begin(), found.end(), false);
	_iter++;
}

目标函数

这部分算法参考论文Fitting B-spline curves to point clouds by curvature-based squared distance minimization[1].

目标函数可以统一写成如下形式：

其中X_k为数据点，P(t)为X_k在参数曲线上的投影点。

1.代数距离

表示为Point Distance Metrics，PDM.目标函数为

这是课上所讲的逼近算法。该问题可以转换成求解

其中A为系数矩阵，P为控制点，D为数据点。记控制点个数为nc，数据点个数为nd，则A的维度是nd x nc，形如下图：

（直接用了课件第十课的图，N替换为A即可）。

本程序没有使用课件上的法方程法（A^T*AP=A^T*D)，而是使用QR分解（householder变换），经测试数值稳定性更好。

由于本程序求解封闭曲线，对于3次B样条，需要让首尾的3个控制点相同，因此将A后三列的bspline系数值加到前三列，通过最小二乘解出nc-3个控制点，再将前3个控制点复制到索引nc,nc+1,nc+2即可。

实现位于./src/Algorithm/CurveFitting/pdm.cu

2.几何距离：TDM

表示为Tangent Distance Metrics，TDM

目标函数为

这个error term代表数据点和投影点连成向量与曲线在投影点处切向量的距离。这属于一阶逼近方法。

实现上与PDM类似，error term对控制点C_i求导得到系数矩阵B，维度是nc x nd。其中

求解线性方程组BAP=BD.由于B的元素为矩阵，实际的实现是将矩阵写成分块矩阵的形式。

3.几何距离：SDM

这是论文[1]提出的新方法，即Squared Distance Metrics，SDM。

error term为

其中ρ为曲线在t_k处的曲率半径。论文作者证明该方法是一种近似二阶逼近。

算法与TDM形式相近，区别是矩阵B的元素变为

TDM和SDM两种方法的代码在./src/Algorithm/CurveFitting/sdm.cu。

测试结果

投影点最大迭代次数统一设置为30. RMSE取曲线拟合n轮迭代中的最小值。

图样：“天”字（2241个数据点，71个控制点）

tian.pts/tian.cpts

初始：

SDM 10轮：

SDM 50轮：

SDM 100轮：

图样：点赞（149个数据点，16个控制点）

thumb.pts/thumb.cpts

可以看出：

1.三种算法随着迭代次数的增加，标准误差都呈下降趋势；误差下降速度上，SDM（近似二阶逼近）>TDM（一阶逼近）>PDM（0阶逼近）。

2.基于几何距离的SDM在较少迭代次数的情况下已经可以得到相对满意的结果。

带*的是异常的结果，如点赞图案的TDM运行20轮得到的参数曲线：

这符合论文[1]中对TDM的结论：不够稳定，往往会收敛到局部最优，造成图样变形。

相比之下，SDM就比较稳定。以下是SDM运行20轮的结果

参考文献

[1]Wang W, Pottmann H, Liu Y. Fitting B-spline curves to point clouds by curvature-based squared distance minimization[J]. ACM Transactions on Graphics (ToG), 2006, 25(2): 214-238.