# BSplineFitting

**Repository Path**: SouthernHermit/bspline-fitting

## Basic Information

- **Project Name**: BSplineFitting
- **Description**: 使用B样条拟合任意闭合平面曲线
- **Primary Language**: C++
- **License**: 0BSD
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 8
- **Forks**: 0
- **Created**: 2022-06-23
- **Last Updated**: 2024-09-10

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# 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格式的配置文件，示例如下：

```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误差。最后输出最终的控制点坐标并绘制曲线。以上述配置文件为例，运行得到如下曲线：

![1642862389153](readme.asset/1642862389153.png)

## 算法

### 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)`。如下图所示。

![1642865003756](readme.assets/1642865003756.png)

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

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

```c++
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].

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

![1642865450338](readme.asset/1642865450338.png)

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

#### 1.代数距离

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

![1642865481195](readme.assets/1642865481195.png)

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

![1642865915869](readme.assets/1642865915869.png)

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

![1642866040912](readme.assets/1642866040912.png)

（直接用了课件第十课的图，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

目标函数为

![1642866359153](readme.assets/1642866359153.png)

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

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

![1642866881916](readme.assets/1642866881916.png)

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

#### 3.几何距离：SDM

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

error term为

![1642867020608](readme.assets/1642867020608.png)

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

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

![1642867144252](readme.assets/1642867144252.png)

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

## 测试结果

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

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

`tian.pts/tian.cpts`

| 迭代次数/ RMSE误差 | PDM      | TDM      | SDM      |
| ------------------ | -------- | -------- | -------- |
| 10                 | 4.144932 | 4.001418 | 1.678290 |
| 50                 | 1.393888 | 1.461655 | 1.424733 |
| 200                | 1.340575 | 1.137016 | 0.404641 |

初始：

![1642869778667](readme.assets/1642869778667.png)

SDM 10轮：

![1642870054889](readme.assets/1642870054889.png)

SDM 50轮：

![1642870108092](readme.assets/1642870108092.png)

SDM 100轮：

![1642870011585](readme.assets/1642870011585.png)

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

`thumb.pts/thumb.cpts`

| 迭代次数/ RMSE误差 | PDM      | TDM       | SDM      |
| ------------------ | -------- | --------- | -------- |
| 5                  | 0.287540 | 0.287544  | 0.285964 |
| 20                 | 0.224001 | 0.224009* | 0.163134 |
| 50                 | 0.123805 | 0.199373* | 0.129117 |

可以看出：

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

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

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

![1642869284062](readme.assets/1642869284062.png)

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

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

![1642869735840](readme.assets/1642869735840.png)



## 参考文献

[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.