1 介绍
Fitting trimmed B-splines（修剪B样条曲线拟合）是一种用于对给定的点云数据进行曲线拟合的算法。该算法使用B样条曲线模型来逼近给定的点云数据，并通过对模型进行修剪来提高拟合的精度和准确性。
B样条曲线是一种常用的曲线表示方法，它通过一组控制点和节点向量来定义曲线的形状。B样条曲线具有局部控制性和平滑性，因此在曲线拟合问题中被广泛应用。
修剪B样条曲线拟合算法的基本步骤如下：
初始化：定义一个B样条曲线模型，并设置模型的阶数、节点向量和控制点。
迭代优化：通过迭代优化的方法，不断调整模型的控制点，使得模型更好地逼近给定的点云数据。在每次迭代中，根据给定的拟合参数和优化目标函数，计算出新的控制点位置。
修剪：根据给定的修剪参数，对拟合的B样条曲线进行修剪。修剪可以通过删除不需要的曲线段或调整曲线段的权重来实现。修剪的目的是提高拟合的准确性和精度。
终止条件：根据设定的终止条件，判断是否终止迭代优化过程。终止条件可以是达到最大迭代次数、拟合精度满足要求或其他自定义条件。
输出结果：输出最终的修剪B样条曲线拟合结果。结果可以是修剪后的B样条曲线模型，也可以是曲线的控制点和节点向量等表示。
2 效果
3 说明
在进行B样条曲面拟合时，通常需要先对曲面进行拟合，然后再对曲面上的曲线进行拟合。这是因为曲面上的曲线通常是曲面的边界或者用于修剪曲面的曲线，它们与曲面的形状紧密相关，因此需要单独进行拟合和调整。
首先，进行B样条曲面拟合时，目标是通过一组控制点和节点向量来逼近给定的点云数据，以生成一个平滑的曲面模型。这个曲面模型可以用来表示曲面的整体形状。
然后，在曲面拟合的基础上，进行B样条曲线拟合。曲线拟合通常是对曲面的边界曲线或者修剪曲线进行拟合。这些曲线与曲面的形状紧密相关，因此需要单独进行拟合和调整。通过对曲线进行拟合，可以更好地捕捉曲面的边界形状或者修剪曲线的形状，从而提高拟合的精度和准确性。
总结起来，进行B样条曲面拟合后，还需要对曲面上的曲线进行拟合的原因是曲线与曲面的形状紧密相关，需要单独进行拟合和调整以提高拟合的精度和准确性。通过分别拟合曲面和曲线，可以更好地捕捉曲面的整体形状和边界形状，从而得到更好的拟合结果。
曲面拟合时的参数
order：曲面的阶数。该参数定义了B样条曲面的阶数，即每个控制点的影响范围。较高的阶数可以提供更大的自由度，但也会增加计算时间和内存消耗。
refinement：细化次数。该参数定义了在拟合过程中进行的细化次数。通过细化，可以在拟合过程中增加新的控制点，以提高拟合的精度和准确性。
iterations：迭代次数。该参数定义了拟合过程中的迭代次数。增加迭代次数可以提高拟合的精度，但也会增加计算时间和内存消耗。
mesh_resolution：网格分辨率。该参数定义了生成曲面网格的分辨率。较高的分辨率会导致更细致的曲面表示，但也会增加计算时间和内存消耗。
params.interior_smoothness：内部平滑度。该参数用于调整曲面内部的平滑度。较大的值会使曲面更平滑。
params.interior_weight：内部权重。该参数用于调整曲面内部点对拟合结果的权重。较大的权重会使内部点对拟合结果的影响更大。
params.boundary_smoothness：边界平滑度。该参数用于调整曲面边界的平滑度。较大的值会使曲面边界更平滑。
params.boundary_weight：边界权重。该参数用于调整曲面边界点对拟合结果的权重。较大的权重会使边界点对拟合结果的影响更大。
曲线拟合的参数
curve_params.addCPsAccuracy：添加控制点的精度。控制点是用于定义B样条曲线形状的关键点。该参数指定了在拟合过程中添加新的控制点的精度。较小的值会导致更精确的拟合结果，但可能会增加计算时间和内存消耗。
curve_params.addCPsIteration：添加控制点的迭代次数。在拟合过程中，可以通过迭代的方式不断添加新的控制点来改进拟合结果。该参数指定了添加控制点的迭代次数。增加迭代次数可以提高拟合的精度，但也会增加计算时间和内存消耗。
curve_params.maxCPs：最大控制点数。该参数限制了拟合过程中允许的最大控制点数。超过这个数目的控制点将被丢弃。通过调整这个参数，可以控制拟合结果的复杂度和平滑度。
curve_params.accuracy：拟合的精度。该参数指定了拟合结果与原始数据之间的误差阈值。较小的值会导致更精确的拟合结果，但可能会增加计算时间和内存消耗。
curve_params.iterations：迭代次数。该参数指定了拟合过程中的迭代次数。增加迭代次数可以提高拟合的精度，但也会增加计算时间和内存消耗。
curve_params.param.closest_point_resolution：最近点的分辨率。该参数用于计算曲线上最近点的分辨率。较小的值会导致更精确的最近点计算，但可能会增加计算时间和内存消耗。
curve_params.param.closest_point_weight：最近点的权重。该参数用于计算曲线上最近点的权重。较大的权重会使最近点对拟合结果的影响更大。
curve_params.param.closest_point_sigma2：最近点的方差。该参数用于计算曲线上最近点的方差。较小的方差会使最近点对拟合结果的影响更大。
curve_params.param.interior_sigma2：内部点的方差。该参数用于计算曲线上内部点的方差。较小的方差会使内部点对拟合结果的影响更大。
curve_params.param.smooth_concavity：平滑凹度。该参数用于调整曲线的平滑凹度。较大的值会使曲线更平滑。
curve_params.param.smoothness：平滑度。该参数用于调整曲线的平滑度。较大的值会使曲线更平滑。
4 代码
#include <pcl/point_cloud.h>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/surface/on_nurbs/fitting_surface_tdm.h>
#include <pcl/surface/on_nurbs/fitting_curve_2d_asdm.h>
#include <pcl/surface/on_nurbs/triangulation.h>
typedef pcl::PointXYZ Point;
std::ostringstream os;
// 将点云数据中的有效点（即非NaN值）提取出来，并以Eigen::Vector3d类型的形式保存起来
void PointCloud2Vector3d(pcl::PointCloud<Point>::Ptr cloud, pcl::on_nurbs::vector_vec3d & data)
{
for (unsigned i = 0; i< cloud->size(); i++)
{
Point &p = cloud->at(i);
if(!std::isnan(p.x) && !std::isnan(p.z) && !std::isnan(p.z))
data.push_back(Eigen::Vector3d(p.x, p.y, p.z));
}
}
// 将ON_NurbsCurve和ON_NurbsSurface的曲线和控制点可视化显示在PCLVisualizer中，方便用户观察和分析曲线的形状和控制点的位置。
void visualizeCurve (ON_NurbsCurve &curve,
ON_NurbsSurface &surface,
pcl::visualization::PCLVisualizer &viewer)
{
// 将曲线转换为点云数据
pcl::PointCloud<pcl::PointXYZRGB>::Ptr curve_cloud(new pcl::PointCloud<pcl::PointXYZRGB>);
pcl::on_nurbs::Triangulation::convertCurve2PointCloud (curve, surface, curve_cloud, 4);//该函数会将曲线上的点均匀地采样，并将采样点作为点云数据的点。
// 可视化
for(std::size_t i=0; i< curve_cloud->size() - 1; i++)
{
pcl::PointXYZRGB &p1 = curve_cloud->at(i);
pcl::PointXYZRGB &p2 = curve_cloud->at(i+1);
os << "line" << i;
viewer.removeShape(os.str());
viewer.addLine<pcl::PointXYZRGB>(p1, p2, 1.0, 0.0, 0.0, os.str());
}
//
pcl::PointCloud<pcl::PointXYZRGB>::Ptr curve_cps (new pcl::PointCloud<pcl::PointXYZRGB>);
for(int i=0; i< curve.CVCount(); i++)
{
ON_3dPoint p1;
curve.GetCV(i, p1); // 曲线的一个控制点
double pnt[3];
surface.Evaluate(p1.x ,p1.y, 0, 3, pnt); // 加usn曲面上对应的点的坐标
pcl::PointXYZRGB p2;
p2.x = float (pnt[0]);
p2.y = float (pnt[1]);
p2.z = float (pnt[2]);
p2.r = 255;
p2.g = 0;
p2.b = 0;
curve_cps->push_back (p2);
}
viewer.removePointCloud ("cloud_cps");
viewer.addPointCloud (curve_cps, "cloud_cps");
}
int main()
{
pcl::visualization::PCLVisualizer viewer("B-spline surface fitting");
viewer.setSize(800, 600);
pcl::PCLPointCloud2 cloud2;
pcl::PointCloud<Point>::Ptr cloud(new pcl::PointCloud<Point>);
if(pcl::io::loadPCDFile("/home/lrj/work/pointCloudData/bun0.pcd", cloud2) == -1)
throw std::runtime_error(" PCD file not found.");
pcl::fromPCLPointCloud2(cloud2, *cloud);
pcl::on_nurbs::NurbsDataSurface data;
PointCloud2Vector3d(cloud, data.interior);
pcl::visualization::PointCloudColorHandlerCustom<Point> handler(cloud, 0, 255, 0);
viewer.addPointCloud<Point>(cloud, handler, "cloud_cylinder");
std::printf(" %lu points in data set\n", cloud->size());
// ############################################################################
// fit B-spline surface
/*
* 整个过程的目的是通过ON-Nurbs算法对给定的点云数据进行曲面拟合，并将拟合结果以三角网格的形式可视化显示出来。
* 通过多次细化和迭代，逐步优化曲面拟合结果，使其更加接近原始点云数据。
*/
// parameters
unsigned order(3); // 曲面的阶数
unsigned refinement(5); // 细化次数
unsigned iterations(10); // 迭代次数
unsigned mesh_resolution(256); // 网格分辨率
pcl::on_nurbs::FittingSurface::Parameter params;
params.interior_smoothness = 0.2; // 内部平滑度
params.interior_weight = 1.0; // 内部权重
params.boundary_smoothness = 0.2; // 边界平滑度
params.boundary_weight = 0.0; // 边界权重
// 生成初始的曲面拟合结果
printf(" surface fitting ...\n");
ON_NurbsSurface nurbs = pcl::on_nurbs::FittingSurface::initNurbsPCABoundingBox(order, &data);
pcl::on_nurbs::FittingSurface fit(&data, nurbs);
// fit.setQuiet (false); // enable/disable debug output
// 将拟合结果转为三角网格，并将其添加到可视化窗口进行现实
pcl::PolygonMesh mesh;
pcl::PointCloud<pcl::PointXYZ>::Ptr mesh_cloud(new pcl::PointCloud<pcl::PointXYZ>);
std::vector<pcl::Vertices> mesh_vertices;
std::string mesh_id = "mesh_nurbs";
pcl::on_nurbs::Triangulation::convertSurface2PolygonMesh(fit.m_nurbs, mesh, mesh_resolution);
viewer.addPolygonMesh(mesh, mesh_id);
// 进行多次曲面细化和求解
for (unsigned i=0; i< refinement; i++)
{
fit.refine(0);
fit.refine(1);
fit.assemble(params);
fit.solve();
pcl::on_nurbs::Triangulation::convertSurface2Vertices (fit.m_nurbs, mesh_cloud, mesh_vertices, mesh_resolution); // 将曲面转为顶点数据
viewer.updatePolygonMesh<pcl::PointXYZ> (mesh_cloud, mesh_vertices, mesh_id); // 视窗刷新，以便观察到曲面拟合的过程
viewer.spinOnce ();
}
// 进行一定次数的曲面拟合迭代
for (unsigned i = 0; i < iterations; i++)
{
fit.assemble (params);
fit.solve ();
pcl::on_nurbs::Triangulation::convertSurface2Vertices (fit.m_nurbs, mesh_cloud, mesh_vertices, mesh_resolution);
viewer.updatePolygonMesh<pcl::PointXYZ> (mesh_cloud, mesh_vertices, mesh_id);
viewer.spinOnce();
}
// ############################################################################
// fit B-spline curve
/*
* 整个曲线拟合的过程的目的是使用ON-Nurbs算法对给定的点云数据进行曲线拟合，并将拟合结果可视化显示出来。
* 通过调整拟合参数，可以控制拟合的精度和平滑度，以得到最优的拟合结果。
*/
// parameters
pcl::on_nurbs::FittingCurve2dAPDM::FitParameter curve_params;
curve_params.addCPsAccuracy = 5e-2; // 添加控制点的精度
curve_params.addCPsIteration = 3; // 添加控制点的迭代次数
curve_params.maxCPs = 200; // 最大控制点数
curve_params.accuracy = 1e-3; // 拟合的精度
curve_params.iterations = 100; // 迭代次数
curve_params.param.closest_point_resolution = 0; // 最近点的分辨率
curve_params.param.closest_point_weight = 1.0; // 最近点的权重
curve_params.param.closest_point_sigma2 = 0.1; // 最近点的方差
curve_params.param.interior_sigma2 = 0.00001; // 内部点的方差
curve_params.param.smooth_concavity = 1.0; // 平滑凹度
curve_params.param.smoothness = 1.0; // 平滑度
// initialisation (circular)
printf (" curve fitting ...\n");
pcl::on_nurbs::NurbsDataCurve2d curve_data;
curve_data.interior = data.interior_param;
curve_data.interior_weight_function.push_back (true);
ON_NurbsCurve curve_nurbs = pcl::on_nurbs::FittingCurve2dAPDM::initNurbsCurve2D (order, curve_data.interior);
// curve fitting
pcl::on_nurbs::FittingCurve2dASDM curve_fit (&curve_data, curve_nurbs);
// curve_fit.setQuiet (false); // enable/disable debug output
curve_fit.fitting (curve_params);
visualizeCurve (curve_fit.m_nurbs, fit.m_nurbs, viewer);
// ############################################################################
// triangulation of trimmed surface
printf (" triangulate trimmed surface ...\n");
viewer.removePolygonMesh (mesh_id);
pcl::on_nurbs::Triangulation::convertTrimmedSurface2PolygonMesh (fit.m_nurbs, curve_fit.m_nurbs, mesh,
mesh_resolution); // 将拟合的曲面、曲线转为三角网格
viewer.addPolygonMesh (mesh, mesh_id);
// save trimmed B-spline surface
if ( fit.m_nurbs.IsValid() )
{
ONX_Model model;
ONX_Model_Object& surf = model.m_object_table.AppendNew();
surf.m_object = new ON_NurbsSurface(fit.m_nurbs);
surf.m_bDeleteObject = true;
surf.m_attributes.m_layer_index = 1;
surf.m_attributes.m_name = "surface";
ONX_Model_Object& curv = model.m_object_table.AppendNew();
curv.m_object = new ON_NurbsCurve(curve_fit.m_nurbs);
curv.m_bDeleteObject = true;
curv.m_attributes.m_layer_index = 2;
curv.m_attributes.m_name = "trimming curve";
// model.Write(file_3dm.c_str());
// printf(" model saved: %s\n", file_3dm.c_str());
}
printf (" ... done.\n");
viewer.spin();
return 0;
}
