diff --git a/doc/tutorials/content/bspline_fitting.rst b/doc/tutorials/content/bspline_fitting.rst
new file mode 100644
index 00000000000..6d7d76e2bd5
--- /dev/null
+++ b/doc/tutorials/content/bspline_fitting.rst
@@ -0,0 +1,276 @@
+.. _bspline_fitting:
+
+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):
+
+.. raw:: html
+
+  <iframe title="Trimmed B-spline surface fitting" width="480" height="390" src="http://www.youtube.com/embed/trH2kWELvyw?rel=0" frameborder="0" allowfullscreen></iframe>
+
+
+Theoretical background
+----------------------
+
+Theoretical information on the algorithm can be found in this `report
+<http://pointclouds.org/blog/trcs/moerwald/index.php>`_ and in my `PhD thesis
+<http://users.acin.tuwien.ac.at/tmoerwald/?site=3>`_.
+
+
+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*.
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :linenos:
+   :lines: 1-169
+
+
+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:
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :linenos:
+   :lines: 56-66
+
+* *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
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 68-72
+
+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.
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 74-80
+
+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.
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 82-92
+
+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.
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 94-102
+
+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.
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 107-120
+
+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.
+
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 122-127
+
+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
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 129-133
+
+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
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 136-142
+
+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:
+
+.. literalinclude:: sources/bspline_fitting/CMakeLists.txt
+   :language: cmake
+   :linenos:   
+
+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:
+
+.. literalinclude:: sources/bspline_fitting/bspline_fitting.cpp
+   :language: cpp
+   :lines: 145-163
+
+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).
+
diff --git a/doc/tutorials/content/images/bspline_bunny.png b/doc/tutorials/content/images/bspline_bunny.png
new file mode 100644
index 00000000000..90fbd1b9cc6
Binary files /dev/null and b/doc/tutorials/content/images/bspline_bunny.png differ
diff --git a/doc/tutorials/content/index.rst b/doc/tutorials/content/index.rst
index 9ddf91a544e..d894e224d69 100644
--- a/doc/tutorials/content/index.rst
+++ b/doc/tutorials/content/index.rst
@@ -1025,6 +1025,22 @@ Surface
      .. |su_3| image:: images/greedy_triangulation.png
                :height: 100px
 
+  * :ref:`bspline_fitting`
+
+     ======  ======
+     |su_4|  Title: **Fitting trimmed B-splines to unordered point clouds**
+
+             Author: *Thomas Mörwald*
+
+             Compatibility: > PCL 1.7
+
+             In this tutorial we will learn how to reconstruct a smooth surface from an unordered point-cloud by fitting trimmed B-splines.
+     ======  ======
+     
+     .. |su_4| image:: images/bspline_bunny.png
+               :height: 100px
+
+
 .. _visualization_tutorial:
 
 Visualization
diff --git a/doc/tutorials/content/sources/bspline_fitting/CMakeLists.txt b/doc/tutorials/content/sources/bspline_fitting/CMakeLists.txt
new file mode 100644
index 00000000000..332515885ea
--- /dev/null
+++ b/doc/tutorials/content/sources/bspline_fitting/CMakeLists.txt
@@ -0,0 +1,13 @@
+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})
+
diff --git a/doc/tutorials/content/sources/bspline_fitting/bspline_fitting.cpp b/doc/tutorials/content/sources/bspline_fitting/bspline_fitting.cpp
new file mode 100644
index 00000000000..b5628ee64fc
--- /dev/null
+++ b/doc/tutorials/content/sources/bspline_fitting/bspline_fitting.cpp
@@ -0,0 +1,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 ("  %zu 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");
+}
diff --git a/examples/surface/CMakeLists.txt b/examples/surface/CMakeLists.txt
index 1507f7871a8..9d105618d00 100644
--- a/examples/surface/CMakeLists.txt
+++ b/examples/surface/CMakeLists.txt
@@ -19,6 +19,10 @@ if(BUILD_surface_on_nurbs)
     PCL_ADD_EXAMPLE(pcl_example_nurbs_fitting_surface
             FILES example_nurbs_fitting_surface.cpp
                     LINK_WITH pcl_common pcl_io pcl_surface pcl_visualization)
+
+    PCL_ADD_EXAMPLE(pcl_example_nurbs_viewer_surface
+            FILES example_nurbs_viewer_surface.cpp
+                    LINK_WITH pcl_common pcl_io pcl_surface pcl_visualization)
                     
     PCL_ADD_EXAMPLE(pcl_example_nurbs_fitting_closed_curve
         	FILES example_nurbs_fitting_closed_curve.cpp
diff --git a/examples/surface/example_nurbs_fitting_surface.cpp b/examples/surface/example_nurbs_fitting_surface.cpp
index 91e61112003..b5628ee64fc 100644
--- a/examples/surface/example_nurbs_fitting_surface.cpp
+++ b/examples/surface/example_nurbs_fitting_surface.cpp
@@ -10,78 +10,32 @@
 typedef pcl::PointXYZ Point;
 
 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));
-  }
-}
+PointCloud2Vector3d (pcl::PointCloud<Point>::Ptr cloud, pcl::on_nurbs::vector_vec3d &data);
 
 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");
-}
+visualizeCurve (ON_NurbsCurve &curve,
+                ON_NurbsSurface &surface,
+                pcl::visualization::PCLVisualizer &viewer);
 
 int
 main (int argc, char *argv[])
 {
-  std::string pcd_file;
+  std::string pcd_file, file_3dm;
 
-  if (argc < 2)
+  if (argc < 3)
   {
-    printf ("\nUsage: pcl_example_nurbs_fitting_surface pcd<PointXYZ>-file\n\n");
+    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];
 
-  unsigned order (3);
-  unsigned refinement (6);
-  unsigned iterations (10);
-  unsigned mesh_resolution (256);
-
-  pcl::visualization::PCLVisualizer viewer ("Test: NURBS surface fitting");
+  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;
@@ -97,27 +51,35 @@ main (int argc, char *argv[])
   printf ("  %zu points in data set\n", cloud->size ());
 
   // ############################################################################
-  // fit NURBS surface
+  // 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);
+  //  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);
 
-  pcl::on_nurbs::FittingSurface::Parameter params;
-  params.interior_smoothness = 0.15;
-  params.interior_weight = 1.0;
-  params.boundary_smoothness = 0.15;
-  params.boundary_weight = 0.0;
-
-  // NURBS refinement
+  // surface refinement
   for (unsigned i = 0; i < refinement; i++)
   {
     fit.refine (0);
@@ -129,7 +91,7 @@ main (int argc, char *argv[])
     viewer.spinOnce ();
   }
 
-  // fitting iterations
+  // surface fitting with final refinement level
   for (unsigned i = 0; i < iterations; i++)
   {
     fit.assemble (params);
@@ -140,44 +102,118 @@ main (int argc, char *argv[])
   }
 
   // ############################################################################
-  // fit NURBS curve
+  // fit B-spline curve
+
+  // parameters
   pcl::on_nurbs::FittingCurve2dAPDM::FitParameter curve_params;
-  curve_params.addCPsAccuracy = 3e-2;
+  curve_params.addCPsAccuracy = 5e-2;
   curve_params.addCPsIteration = 3;
   curve_params.maxCPs = 200;
   curve_params.accuracy = 1e-3;
-  curve_params.iterations = 10000;
+  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.0001;
+  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);
-
-  // ############################### FITTING ###############################
+  // 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");
+}
diff --git a/examples/surface/example_nurbs_viewer_surface.cpp b/examples/surface/example_nurbs_viewer_surface.cpp
new file mode 100644
index 00000000000..710cfac9ddc
--- /dev/null
+++ b/examples/surface/example_nurbs_viewer_surface.cpp
@@ -0,0 +1,88 @@
+#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;
+
+int
+main (int argc, char *argv[])
+{
+  std::string file_3dm;
+
+  if (argc < 2)
+  {
+    printf ("\nUsage: pcl_example_nurbs_viewer_surface 3dm-out-file\n\n");
+    exit (0);
+  }
+  file_3dm = argv[1];
+
+  pcl::visualization::PCLVisualizer viewer ("B-spline surface viewer");
+  viewer.setSize (800, 600);
+
+  int mesh_resolution = 128;
+
+  ON::Begin();
+
+  // load surface
+  ONX_Model on_model;
+  bool rc = on_model.Read(file_3dm.c_str());
+
+  // print diagnostic
+  if ( rc )
+    std::cout << "Successfully read: " << file_3dm << std::endl;
+  else
+    std::cout << "Errors during reading: " << file_3dm << std::endl;
+
+//  ON_TextLog out;
+//  on_model.Dump(out);
+
+  if(on_model.m_object_table.Count()==0)
+  {
+    std::cout << "3dm file does not contain any objects: " << file_3dm << std::endl;
+    return -1;
+  }
+
+  const ON_Object* on_object = on_model.m_object_table[0].m_object;
+  if(on_object==NULL)
+  {
+    std::cout << "object[0] not valid." << std::endl;
+    return -1;
+  }
+
+  const ON_NurbsSurface& on_surf = *(ON_NurbsSurface*)on_object;
+
+  pcl::PolygonMesh mesh;
+  std::string mesh_id = "mesh_nurbs";
+  if(on_model.m_object_table.Count()==1)
+  {
+    std::cout << "3dm file does not contain a trimming curve: " << file_3dm << std::endl;
+
+
+    pcl::on_nurbs::Triangulation::convertSurface2PolygonMesh (on_surf, mesh, mesh_resolution);
+  }
+  else
+  {
+    on_object = on_model.m_object_table[1].m_object;
+    if(on_object==NULL)
+    {
+      std::cout << "object[1] not valid." << std::endl;
+      return -1;
+    }
+
+    const ON_NurbsCurve& on_curv = *(ON_NurbsCurve*)on_object;
+
+    pcl::on_nurbs::Triangulation::convertTrimmedSurface2PolygonMesh (on_surf, on_curv, mesh,
+                                                                     mesh_resolution);
+  }
+
+  viewer.addPolygonMesh (mesh, mesh_id);
+
+  viewer.spin ();
+  return 0;
+}
+