Fitting trimmed B-splines to unordered point clouds

This tutorial explains how to run a B-spline fitting algorithm on a point-cloud, to obtain a smooth, parametric surface representation. The algorithm consists of the following steps:

  • Initialization of the B-spline surface by using the Principal Component Analysis (PCA). This assumes that the point-cloud has two main orientations, i.e. that it is roughly planar.
  • Refinement and fitting of the B-spline surface.
  • Circular initialization of the B-spline curve. Here we assume that the point-cloud is compact, i.e. no separated clusters.
  • Fitting of the B-spline curve.
  • Triangulation of the trimmed B-spline surface.

In this video, the algorithm is applied to the frontal scan of the stanford bunny (204800 points):

Theoretical background

Theoretical information on the algorithm can be found in this report and in my PhD thesis.

PCL installation settings

Please note that the modules for NURBS and B-splines are not enabled by default. Make sure you enable “BUILD_surface_on_nurbs” in your ccmake configuration, by setting it to ON.

If your license permits, also enable “USE_UMFPACK” for sparse linear solving. This requires SuiteSparse (libsuitesparse-dev in Ubuntu) which is faster, allows more degrees of freedom (i.e. control points) and more data points.

The program created during this tutorial is available in pcl/examples/surface/example_nurbs_fitting_surface.cpp and is built when “BUILD_examples” is set to ON. This will create the binary called pcl_example_nurbs_fitting_surface in your bin folder.

The code

The cpp file used in this tutorial can be found in pcl/doc/tutorials/content/sources/bspline_fitting/bspline_fitting.cpp. You can find the input file at pcl/test/bunny.pcd.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
#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;

void
PointCloud2Vector3d (pcl::PointCloud<Point>::Ptr cloud, pcl::on_nurbs::vector_vec3d &data);

void
visualizeCurve (ON_NurbsCurve &curve,
                ON_NurbsSurface &surface,
                pcl::visualization::PCLVisualizer &viewer);

int
main (int argc, char *argv[])
{
  std::string pcd_file, file_3dm;

  if (argc < 3)
  {
    printf ("\nUsage: pcl_example_nurbs_fitting_surface pcd<PointXYZ>-in-file 3dm-out-file\n\n");
    exit (0);
  }
  pcd_file = argv[1];
  file_3dm = argv[2];

  pcl::visualization::PCLVisualizer viewer ("B-spline surface fitting");
  viewer.setSize (800, 600);

  // ############################################################################
  // load point cloud

  printf ("  loading %s\n", pcd_file.c_str ());
  pcl::PointCloud<Point>::Ptr cloud (new pcl::PointCloud<Point>);
  pcl::PCLPointCloud2 cloud2;
  pcl::on_nurbs::NurbsDataSurface data;

  if (pcl::io::loadPCDFile (pcd_file, cloud2) == -1)
    throw std::runtime_error ("  PCD file not found.");

  fromPCLPointCloud2 (cloud2, *cloud);
  PointCloud2Vector3d (cloud, data.interior);
  pcl::visualization::PointCloudColorHandlerCustom<Point> handler (cloud, 0, 255, 0);
  viewer.addPointCloud<Point> (cloud, handler, "cloud_cylinder");
  printf ("  %lu points in data set\n", cloud->size ());

  // ############################################################################
  // fit B-spline surface

  // 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;

  // initialize
  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

  // mesh for visualization
  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);

  // surface refinement
  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 ();
  }

  // surface fitting with final refinement level
  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

  // 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;
}

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 (!pcl_isnan (p.x) && !pcl_isnan (p.y) && !pcl_isnan (p.z))
      data.push_back (Eigen::Vector3d (p.x, p.y, p.z));
  }
}

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);
    std::ostringstream os;
    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);
    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");
}

The explanation

Now, let’s break down the code piece by piece. Lets start with the choice of the parameters for B-spline surface fitting:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
  // 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;
  • order is the polynomial order of the B-spline surface.
  • refinement is the number of refinement iterations, where for each iteration control-points are inserted, approximately doubling the control points in each parametric direction of the B-spline surface.
  • iterations is the number of iterations that are performed after refinement is completed.
  • mesh_resolution the number of vertices in each parametric direction, used for triangulation of the B-spline surface.

Fitting:

  • interior_smoothness is the smoothness of the surface interior.
  • interior_weight is the weight for optimization for the surface interior.
  • boundary_smoothness is the smoothness of the surface boundary.
  • boundary_weight is the weight for optimization for the surface boundary.

Note, that the boundary in this case is not the trimming curve used later on. The boundary can be used when a point-set exists that defines the boundary. Those points can be declared in pcl::on_nurbs::NurbsDataSurface::boundary. In that case, when the boundary_weight is greater than 0.0, the algorithm tries to align the domain boundaries to these points. In our example we are trimming the surface anyway, so there is no need for aligning the boundary.

Initialization of the B-spline surface

  // initialize
  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

The command initNurbsPCABoundingBox uses PCA to create a coordinate systems, where the principal eigenvectors point into the direction of the maximum, middle and minimum extension of the point-cloud. The center of the coordinate system is located at the mean of the points. To estimate the extension of the B-spline surface domain, a bounding box is computed in the plane formed by the maximum and middle eigenvectors. That bounding box is used to initialize the B-spline surface with its minimum number of control points, according to the polynomial degree chosen.

The surface fitting class pcl::on_nurbs::FittingSurface is initialized with the point data and the initial B-spline.

  // mesh for visualization
  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);

The on_nurbs::Triangulation class allows easy conversion between the ON_NurbsSurface and the PolygonMesh class, for visualization of the B-spline surfaces. Note that NURBS are a generalization of B-splines, and are therefore a valid container for B-splines, with all control-point weights = 1.

  // surface refinement
  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 ();
  }

Refinement and fitting of the B-spline surface

At this point of the code we have a B-spline surface with minimal number of control points. Typically they are not enough to represent finer details of the underlying geometry of the point-cloud. However, if we increase the control-points to our desired level of detail and subsequently fit the refined B-spline, we run into problems. For robust fitting B-spline surfaces the rule is: “The higher the degree of freedom of the B-spline surface, the closer we have to be to the points to be approximated”.

This is the reason why we iteratively increase the degree of freedom by refinement in both directions (line 85-86), and fit the B-spline surface to the point-cloud, getting closer to the final solution.

  // surface fitting with final refinement level
  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 ();
  }

After we reached the final level of refinement, the surface is further fitted to the point-cloud for a pleasing end result.

Initialization of the B-spline curve

Now that we have the surface fitted to the point-cloud, we want to cut off the overlapping regions of the surface. To achieve this we project the point-cloud into the parametric domain using the closest points to the B-spline surface. In this domain of R^2 we perform the weighted B-spline curve fitting, that creates a closed trimming curve that approximately contains all the points.

  // 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;

The topic of curve fitting goes a bit deeper into the thematics of B-splines. Here we assume that you are familiar with the concept of B-splines, knot vectors, control-points, and so forth. Please consider the curve being split into supporting regions which is bound by consecutive knots. Also note that points that are inside and outside the curve are distinguished.

  • addCPsAccuracy the distance of the supporting region of the curve to the closest data points has to be below this value, otherwise a control point is inserted.
  • addCPsIteration inner iterations without inserting control points.
  • maxCPs the maximum total number of control-points.
  • accuracy the average fitting accuracy of the curve, w.r.t. the supporting regions.
  • iterations maximum number of iterations performed.
  • closest_point_resolution number of control points that must lie within each supporting region. (0 turns this constraint off)
  • closest_point_weight weight for fitting the curve to its closest points.
  • closest_point_sigma2 threshold for closest points (disregard points that are further away from the curve).
  • interior_sigma2 threshold for interior points (disregard points that are further away from and lie within the curve).
  • smooth_concavity value that leads to inward bending of the curve (0 = no bending; <0 inward bending; >0 outward bending).
  • smoothness weight of smoothness term.
  // 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);

The curve is initialized using a minimum number of control points to represent a circle, with the center located at the mean of the point-cloud and the radius of the maximum distance of a point to the center. Please note that interior weighting is enabled for all points with the command curve_data.interior_weight_function.push_back (true).

Fitting of the B-spline curve

  // 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);

Similar to the surface fitting approach, the curve is iteratively fitted and refined, as shown in the video. Note how the curve tends to bend inwards at regions where it is not supported by any points.

Triangulation of the trimmed B-spline surface

  // 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);

After the curve fitting terminated, our geometric representation consists of a B-spline surface and a closed B-spline curved, defined within the parametric domain of the B-spline surface. This is called trimmed B-spline surface. In line 140 we can use the trimmed B-spline to create a triangular mesh. The triangulation algorithm first triangulates the whole domain and afterwards removes triangles that lie outside of the trimming curve. Vertices of triangles that intersect the trimming curve are clamped to the curve.

When running this example and switch to wire-frame mode (w), you will notice that the triangles are ordered in a rectangular way, which is a result of the rectangular domain of the surface.

Some hints

Please bear in mind that the robustness of this algorithm heavily depends on the underlying data. The parameters for B-spline fitting are designed to model the characteristics of this data.

  • If you have holes or steps in your data, you might want to work with lower refinement levels and lower accuracy to prevent the B-spline from folding and twisting. Moderately increasing of the smoothness might also work.
  • Try to introduce as much pre-conditioning and constraints to the parameters. E.g. if you know, that the trimming curve is rather simple, then limit the number of maximum control points.
  • Start simple! Before giving up on gaining control over twisting and bending B-splines, I highly recommend to start your fitting trials with a small number of control points (low refinement), low accuracy but also low smoothness (B-splines have implicit smoothing property).

Compiling and running the program

Add the following lines to your CMakeLists.txt file:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)

project(bspline_fitting)

find_package(PCL 1.7 REQUIRED)

include_directories(${PCL_INCLUDE_DIRS})
link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})

add_executable (bspline_fitting bspline_fitting.cpp)
target_link_libraries (bspline_fitting ${PCL_LIBRARIES})

After you have made the executable, you can run it. Simply do:

$ ./bspline_fitting ${PCL_ROOT}/test/bunny.pcd

Saving and viewing the result

  • Saving as OpenNURBS (3dm) file

You can save the B-spline surface by using the commands provided by OpenNurbs:

  // 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());
  }

The files generated can be viewed with the pcl/examples/surface/example_nurbs_viewer_surface.cpp.

  • Saving as triangle mesh into a vtk file

You can save the triangle mesh for example by saving into a VTK file by:

#include <pcl/io/vtk_io.h> … pcl::io::saveVTKFile (“mesh.vtk”, mesh);

PCL also provides vtk conversion into other formats (PLY, OBJ).