Pierre Alliez, Laurent Saboret,
Gaël Guennebaud
45.1 | Introduction | ||||
45.2 | Common Reconstruction Pipeline | ||||
45.3 | Poisson | ||||
45.3.1 Interface | |||||
45.3.2 Example | |||||
45.4 | Contouring | ||||
45.5 | Output |
This Cgal component implements a surface reconstruction method which takes as input point sets with oriented normals and computes an implicit function. We assume that the input points contain no outliers and little noise. The output surface mesh is generated by extracting an isosurface of this function with the Cgal Surface Mesh Generator [RY07] or potentially with any other surface contouring algorithm.
Figure 45.1: Poisson surface reconstruction. Left: 17K points sampled on the statue of an elephant with a Minolta laser scanner. Right: reconstructed surface mesh.
More specifically, the core surface reconstruction algorithm consists of computing an implicit function which is an approximate indicator function of the inferred solid (Poisson Surface Reconstruction - referred to as Poisson). Poisson is a two steps process: it requires solving for the implicit function before function evaluation.
Surface reconstruction from point sets is often a sequential process with the following steps: 1) Scanning and scan alignment produce a set of points or points with normals; 2) Outlier removal; 3) Simplification to reduce the number of input points; 4) Smoothing to reduce noise in the input data; 5) Normal estimation and orientation when the normals are not already provided by the acquisition device; and 6) Surface reconstruction.
Cgal provides algorithms for all steps listed above except alignment.
Chapter Point_set_processing_3 54 describes algorithms to pre-process the point set before reconstruction with functions devoted to the simplification, outlier removal, smoothing, normal estimation and normal orientation.
Figure 45.2: Common surface reconstruction pipeline.
Given a set of 3D points with oriented normals (denoted oriented points in the sequel) sampled on the boundary of a 3D solid, the Poisson Surface Reconstruction method [KBH05] solves for an approximate indicator function of the inferred solid, whose gradient best matches the input normals. The output scalar function, represented in an adaptive octree, is then iso-contoured using an adaptive marching cubes.
Cgal implements a variant of this algorithm which solves for a piecewise linear function on a 3D Delaunay triangulation instead of an adaptive octree. The algorithm takes as input a set of 3D oriented points. It builds a 3D Delaunay triangulation from these points and refines it by Delaunay refinement so as to remove all badly shaped (non isotropic) tetrahedra and to tessellate a loose bounding box of the input oriented points. The normal of each Steiner point added during refinement is set to zero. It then solves for a scalar indicator function f represented as a piecewise linear function over the refined triangulation. More specifically, it solves for the Poisson equation Δf = div(n) at each vertex of the triangulation using a sparse linear solver. Eventually, the Cgal surface mesh generator extracts an isosurface with function value set by default to be the median value of f at all input points.
The class template declaration is:
template<
class Gt>
class Poisson_reconstruction_function;
with
Gt: Geometric traits class.
Creation:
| ||||||
|
| |||||
Creates a Poisson implicit function from the [first, beyond) range of points.
|
The main operations are:
See details in
CGAL::Poisson_reconstruction_function<GeomTraits>
poisson_reconstruction_example.cpp reads a point set, creates a Poisson implicit function and reconstructs a surface.
File: examples/Surface_reconstruction_points_3/poisson_reconstruction_example.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Surface_mesh_default_triangulation_3.h> #include <CGAL/make_surface_mesh.h> #include <CGAL/Implicit_surface_3.h> #include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h> #include <CGAL/Poisson_reconstruction_function.h> #include <CGAL/Point_with_normal_3.h> #include <CGAL/property_map.h> #include <CGAL/IO/read_xyz_points.h> #include <vector> #include <fstream> // Types typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::FT FT; typedef Kernel::Point_3 Point; typedef CGAL::Point_with_normal_3<Kernel> Point_with_normal; typedef Kernel::Sphere_3 Sphere; typedef std::vector<Point_with_normal> PointList; typedef CGAL::Poisson_reconstruction_function<Kernel> Poisson_reconstruction_function; typedef CGAL::Surface_mesh_default_triangulation_3 STr; typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3; typedef CGAL::Implicit_surface_3<Kernel, Poisson_reconstruction_function> Surface_3; int main(void) { // Poisson options FT sm_angle = 20.0; // Min triangle angle in degrees. FT sm_radius = 0.03; // Max triangle size w.r.t. point set radius. FT sm_distance = 0.003; // Surface approximation error w.r.t. p.s.r. // Reads the point set file in points[]. // Note: read_xyz_points_and_normals() requires an iterator over points // + property maps to access each point's position and normal. // The position property map can be omitted here as we use iterators over Point_3 elements. PointList points; std::ifstream stream("data/kitten.xyz"); if (!stream || !CGAL::read_xyz_points_and_normals( stream, std::back_inserter(points), CGAL::make_normal_of_point_with_normal_pmap(std::back_inserter(points)))) { std::cerr << "Error: cannot read file data/kitten.xyz" << std::endl; return EXIT_FAILURE; } // Creates implicit function from the read points. // Note: this method requires an iterator over points // + property maps to access each point's position and normal. // The position property map can be omitted here as we use iterators over Point_3 elements. Poisson_reconstruction_function function( points.begin(), points.end(), CGAL::make_normal_of_point_with_normal_pmap(points.begin())); // Computes the Poisson indicator function f() // at each vertex of the triangulation. if ( ! function.compute_implicit_function() ) return EXIT_FAILURE; // Gets one point inside the implicit surface // and computes implicit function bounding sphere radius. Point inner_point = function.get_inner_point(); Sphere bsphere = function.bounding_sphere(); FT radius = std::sqrt(bsphere.squared_radius()); // Defines the implicit surface: requires defining a // conservative bounding sphere centered at inner point. FT sm_sphere_radius = 2.01 * radius; FT sm_dichotomy_error = sm_distance/10.0; // Dichotomy error must be << sm_distance Surface_3 surface(function, Sphere(inner_point,sm_sphere_radius*sm_sphere_radius), sm_dichotomy_error); // Defines surface mesh generation criteria CGAL::Surface_mesh_default_criteria_3<STr> criteria(sm_angle, // Min triangle angle (degrees) sm_radius*radius, // Max triangle size sm_distance*radius); // Approximation error // Generates surface mesh with manifold option STr tr; // 3D Delaunay triangulation for surface mesh generation C2t3 c2t3(tr); // 2D complex in 3D Delaunay triangulation CGAL::make_surface_mesh(c2t3, // reconstructed mesh surface, // implicit surface criteria, // meshing criteria CGAL::Manifold_tag()); // require manifold mesh with no boundary if(tr.number_of_vertices() == 0) return EXIT_FAILURE; // saves reconstructed surface mesh std::ofstream out("kitten_poisson_0.003.off"); CGAL::output_surface_facets_to_off(out, c2t3); return EXIT_SUCCESS; }