namespace CGAL {
/*!

\mainpage User Manual 
\anchor Chapter_Poisson_Surface_Reconstruction
\cgalAutoToc

\authors Pierre Alliez, Laurent Saboret, Ga&euml;l Guennebaud

\section Poisson_surface_reconstruction_3Introduction Introduction

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
\cgalCite{cgal:ry-gsddrm-06} or potentially with any other surface
contouring algorithm.

\cgalFigureBegin{Poisson_surface_reconstruction_3figintroduction,introduction.jpg}
Poisson surface reconstruction.\n
Left: 17K points sampled on the statue of an elephant with a Minolta laser scanner. Right: reconstructed surface mesh.
\cgalFigureEnd

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.

\section Poisson_surface_reconstruction_3Common Common Reconstruction Pipeline

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 \ref chappoint_set_processing_3 "Point Set Processing"
describes algorithms to pre-process the point set before
reconstruction with functions devoted to the simplification, outlier
removal, smoothing, normal estimation and normal orientation.

\cgalFigureBegin{Poisson_surface_reconstruction_3figpipeline,pipeline.jpg}
Common surface reconstruction pipeline.
\cgalFigureEnd

\section Poisson_surface_reconstruction_3Poisson Poisson

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 \cgalCite{Kazhdan06} 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$ f\f$ represented as a piecewise linear function
over the refined triangulation. More specifically, it solves for the
Poisson equation \f$ \Delta f = div(\mathbf{n})\f$ 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$ f\f$ at all input points.

\section Poisson_surface_reconstruction_3Function Reconstruction Function

A global function `poisson_surface_reconstruction_delaunay()` is
provided. It takes points with normals as input and handles the whole
reconstruction pipeline :

- it computes the implicit function
- it reconstructs the surface with a given precision using the \cgal
surface mesh generator based on Delaunay refinement
\cgalCite{cgal:ry-gsddrm-06} \cgalCite{cgal:bo-pgsms-05}
- it outputs the result in a polygon mesh.

This function aims at providing a quick and user-friendly API for
Poisson reconstruction. Advanced users may be interested in using the
class (see \ref Poisson_surface_reconstruction_3Class) which allows
them, for example, to use another surface mesher or a different output
structure.

\subsection Poisson_surface_reconstruction_3Example_function Example

The following example reads a point set and reconstructs a surface using Poisson reconstruction.

\cgalExample{Poisson_surface_reconstruction_3/poisson_reconstruction_function.cpp}

\section Poisson_surface_reconstruction_3Class Reconstruction Class

The class template declaration is `template<class Gt> class Poisson_reconstruction_function` where
`Gt` is a geometric traits class.

For details see: `Poisson_reconstruction_function<GeomTraits>`

\subsection Poisson_surface_reconstruction_3Example_class Example

The following example reads a point set, creates a Poisson implicit function and reconstructs a surface.

\cgalExample{Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp}

\subsection Poisson_surface_reconstruction_3Contouring Contouring


The computed implicit functions can be iso-contoured to reconstruct a
surface by using the \cgal surface mesh generator 
\cgalCite{cgal:ry-gsddrm-06} \cgalCite{cgal:bo-pgsms-05} :

`make_surface_mesh()` 

The parameter `Tag` affects the behavior of `make_surface_mesh()`: 
- `Manifold_tag`: the output mesh is guaranteed to be a manifold surface without boundary.
- `Manifold_with_boundary_tag`: the output mesh is guaranteed to be manifold and may have boundaries.
- `Non_manifold_tag`: the output mesh has no guarantee and hence is outputted as a polygon soup.

\subsection Poisson_surface_reconstruction_3Output Output

The surface reconstructed by `make_surface_mesh()` is required to be a
model of the concept `SurfaceMeshComplex_2InTriangulation_3`, a data
structure devised to represent a two dimensional complex embedded into
a three dimensional triangulation.

`SurfaceMeshComplex_2InTriangulation_3` defines the methods to traverse the reconstructed surface, and e.g. convert it to a triangle soup.

Other \cgal components provide functions to write the reconstructed
surface mesh to the %Object File Format (OFF) \cgalCite{cgal:p-gmgv16-96}
and to convert it to a polyhedron (when it is manifold):
- `output_surface_facets_to_off()` 
- `output_surface_facets_to_polyhedron()` 

See \ref Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp "poisson_reconstruction_example.cpp" example above.

\section surface_reconstruction_section_case_studies Case Studies

The surface reconstruction problem being inherently ill-posed, the
proposed algorithm does not pretend to reconstruct all kinds of
surfaces with arbitrary sampling conditions. This section provides the
user with some hints about the ideal sampling and contouring
conditions, and depicts some failure cases when these conditions are
not matched.

\subsection Poisson_surface_reconstruction_3IdealConditions Ideal Conditions

The user must keep in mind that the Poisson surface reconstruction
algorithm comprises two phases (computing the implicit function from
the input point set and contouring an iso-surface of this
function). Both require some care in terms of sampling conditions and
parameter tuning.

\subsection Poisson_surface_reconstruction_3PointSet Point Set

Ideally the current implementation of the Poisson surface
reconstruction method expects a dense 3D oriented point set (typically
matching the epsilon-sampling condition \cgalCite{cgal:bo-pgsms-05}) and
sampled over a closed, smooth surface. Oriented herein means that all
3D points must come with consistently oriented normals to the inferred
surface.  \cgalFigureRef{Poisson_surface_reconstruction_3figbimba} and \cgalFigureRef{Poisson_surface_reconstruction_3figeros} 
illustrate cases where these ideal conditions are met.

\cgalFigureBegin{Poisson_surface_reconstruction_3figbimba,bimba.jpg}
Poisson reconstruction. Left: 120K points sampled on a statue (Minolta laser scanner). Right: reconstructed surface mesh.
\cgalFigureEnd

\cgalFigureBegin{Poisson_surface_reconstruction_3figeros,eros.jpg}
Left: 120K points sampled on a statue (Minolta laser scanner). Right: reconstructed surface mesh.
\cgalFigureEnd

The algorithm is fairly robust to anisotropic sampling and to
noise. It is also robust to missing data through filling the
corresponding holes as the algorithm is designed to reconstruct the
indicator function of an inferred solid (see
\cgalFigureRef{Poisson_surface_reconstruction_3figholes_good}).

\cgalFigureBegin{Poisson_surface_reconstruction_3figholes_good,holes_good.jpg}
Top left: 65K points sampled on a hand (Kreon laser scanner). Bottom left: the point set is highly anisotropic due to the scanning technology. Right: reconstructed surface mesh and closeup. The holes are properly closed.
\cgalFigureEnd

The algorithm is in general not robust to outliers, although a few
outliers do not always create a failure, see
\cgalFigureRef{Poisson_surface_reconstruction_3figoutliers}.

\cgalFigureBegin{Poisson_surface_reconstruction_3figoutliers,outliers.jpg}
Left: 70K points sampled on an elephant with few outliers emphasized with disks. Right: reconstructed surface mesh.
\cgalFigureEnd

The algorithm works well even when the inferred surface is composed of
several connected components, provided that both all normals are
properly estimated and oriented (the current \cgal normal orienter
algorithm may fail in some cases, see `mst_orient_normals()`),
and that the final contouring algorithm is properly seeded for each
component. When the inferred surface is composed of several nested
connected components care should be taken to orient the normals of
each component in alternation (inward/outward) so that the final
contouring stage picks a proper contouring value.

\subsection Poisson_surface_reconstruction_3ContouringP Contouring Parameters

Our implementation of the Poisson surface reconstruction algorithm
computes an implicit function represented as a piecewise linear
function over the tetrahedra of a 3D Delaunay triangulation
constructed from the input points then refined through Delaunay
refinement. For this reason, any iso-surface is also piecewise linear
and hence may contain sharp creases. As the contouring algorithm
`make_surface_mesh()` expects a smooth implicit function these
sharp creases may create spurious clusters of vertices in the final
reconstructed surface mesh when setting a small mesh sizing or surface
approximation error parameter (see
\cgalFigureRef{Poisson_surface_reconstruction_3figcontouring_bad}).

One way to avoid these spurious clusters consists of adjusting the
mesh sizing and surface approximation parameters large enough compared
to the average sampling density (obtained through
`compute_average_spacing()`) so that the contouring algorithm
<i>perceives</i> a smooth iso-surface. We recommend to use the following
contouring parameters:

- Max triangle radius: at least 100 times the average spacing.
- Approximation distance: at least 0.25 times the average spacing.

\cgalFigureBegin{Poisson_surface_reconstruction_3figcontouring_bad,contouring_bad.jpg}
Left: surface reconstructed with approximation distance = 0.25 * average spacing. Right: surface reconstructed with approximation distance = 0.15 * average spacing. Notice the spurious cluster on the cheek.
\cgalFigureEnd

\subsection Poisson_surface_reconstruction_3DegradedConditions Degraded Conditions

The conditions listed above are rather restrictive and in practice not
all of them are met in the applications. We now illustrates the
behavior of the algorithm when the conditions are not met in terms of
sampling, wrongly oriented normals, noise and sharp creases.

\subsection Poisson_surface_reconstruction_3SparseSampling Sparse Sampling

The reconstruction algorithm expects a sufficiently dense point
set. Although there is no formal proof of correctness of the algorithm
under certain density conditions due to its variational nature, our
experiments show that the algorithm reconstructs well all thin
features when the local spacing is at most one tenth of the local
feature size (the distance to the medial axis, which captures
altogether curvature, thickness and separation). When this condition
is not met the reconstruction does not reconstruct the thin
undersampled features (see
\cgalFigureRef{Poisson_surface_reconstruction_3figsampling}).

\cgalFigureBegin{Poisson_surface_reconstruction_3figsampling,sampling.jpg}
Left: 50K points sampled on the Neptune trident. The reconstruction (not shown) is successful in this case. Right: point set simplified to 1K points then reconstructed (all input points are depicted with normals). The thin feature is not reconstructed.
\cgalFigureEnd

\subsection Poisson_surface_reconstruction_3LargeHoles Large Holes

The reconstruction is devised to solve for an implicit function which
is an approximate indicator function of an inferred solid. For this
reason the contouring algorithm always extracts a closed surface mesh
and hence is able to fill the small holes where data are missing due,
e.g., to occlusions during acquisition (see
\cgalFigureRef{Poisson_surface_reconstruction_3figholes_bad}).

\cgalFigureBegin{Poisson_surface_reconstruction_3figholes_bad,holes_bad.jpg}
Left: 65K points sampled on a hand with no data captured at the wrist base. Right: reconstructed surface mesh. The surface is properly closed on the fingers and also closed at the wrist but in a less plausible manner.
\cgalFigureEnd



In case of large holes the algorithm still closes them all but the resulting piecewise
linear implicit function may exhibit large triangle patches and sharp creases as the 3D
Delaunay triangulation used for solving is very coarse where the holes are filled.
This can be avoided by a two pass approach. The first pass for a subset of the points
serves to get an approximation of the surface at the holes. This surface then serves to
compute a smoother 3D Delaunay triangulation for the second pass with the full set of points.

\cgalFigureBegin{Poisson_surface_reconstruction_3-fig-two_passes,two-passes.png}
Left: The wrist. Middle: one pass. Right: two passes.
\cgalFigureEnd

\subsection Poisson_surface_reconstruction_3WronglyOriented Wrongly Oriented Normals

The Poisson surface reconstruction approaches solves for an implicit
function whose gradient best matches a set of input normals. Because
it solves this problem in the least squares sense, it is robust to few
isolated wrongly oriented (flipped) normals. Nevertheless a cluster of
wrongly oriented normals leads to an incorrect implicit function and
hence to spurious geometric or even topological distortion (see
\cgalFigureRef{Poisson_surface_reconstruction_3figflipped_normals}).

\cgalFigureBegin{Poisson_surface_reconstruction_3figflipped_normals,flipped_normals.jpg}
Left: points sampled on a sphere with a cluster of wrongly oriented normals. Right: reconstructed surface mesh with a spurious bump.
\cgalFigureEnd

\subsection Poisson_surface_reconstruction_3NoiseandOutliers Noise and Outliers

A large amount of noise inevitably impacts on the reconstruction (see
\cgalFigureRef{Poisson_surface_reconstruction_3fignoise}, top) and the
current implementation does not provide any mean to trade data fitting
for smoothness. Nevertheless if the signal-to-noise ratio is
sufficiently high and/or the surface approximation and sizing
parameters set for contouring the iso-surface is large with respect to
the noise level the output surface mesh will appear smooth (not
shown). If the user wants to produce a smooth and detailed output
surface mesh, we recommend to apply smoothing through
`jet_smooth_point_set()` (see
\cgalFigureRef{Poisson_surface_reconstruction_3fignoise}, bottom).

\cgalFigureBegin{Poisson_surface_reconstruction_3fignoise,noise.jpg}
Top-left: points sampled on a sphere and corrupted with a lot of noise. Top-right: reconstructed surface mesh. Bottom-left: smoothed point set. Bottom-right: reconstructed surface mesh.
\cgalFigureEnd

For a large number of outliers the failure cases (not shown) translate
into spurious small connected components and massive distortion near
the inferred surface. In this case the outliers must be removed
through `remove_outliers()`.

\subsection Poisson_surface_reconstruction_3SharpCreases Sharp Creases

The current reconstruction algorithm is not able to recover the sharp
creases and corners present in the inferred surface. This translates
into smoothed sharp creases.

\cgalFigureBegin{Poisson_surface_reconstruction_3figsharp_features,sharp_features.jpg}
Left: 5K points sampled on a mechanical piece with sharp features (creases, darts and corners). Right: reconstructed surface mesh with smoothed creases.
\cgalFigureEnd

\section SurfReconstPerformances Performances

We provide some performance numbers for scanning data. We measure the Poisson implicit function computation time,
the contouring time for a range of approximation distances, the memory occupancy as well as the influence of
the point set simplification. The machine used is a PC running Windows 7 64 bits with an Intel CPU Core 2 Duo
processor clocked at 2.81 GHz and with 8 GB of RAM. The software is compiled with Visual \CC 2010 (VC9) compiler
with the 03 option which maximizes speed.  All measurements were done using the \ref thirdpartyEigen "Eigen" library.


\subsection SurfReconstPerfPIF Poisson Implicit Function

The point set chosen for benchmarking the Poisson implicit function is the Bimba con Nastrino point set
(1.6 million points) depicted by \cgalFigureRef{Poisson_surface_reconstruction_3-fig-contouring_bench}.
We measure the Poisson implicit function computation (i.e., the call to 
`Poisson_reconstruction_function::compute_implicit_function()` denoted by Poisson solve hereafter)
for this point set as well as for simplified versions obtained through random simplification.
The following table provides Poisson solve computation times in seconds for an increasing number of points.

<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>

Number of points (x1000) 
<TD class="math" ALIGN=CENTER NOWRAP>
Poisson solve duration (in s) 
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
15
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
100
<TD class="math" ALIGN=CENTER NOWRAP>
25
<TR>

<TD class="math" ALIGN=CENTER NOWRAP>
250
<TD class="math" ALIGN=CENTER NOWRAP>
96
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
<TD class="math" ALIGN=CENTER NOWRAP>
150
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
<TD class="math" ALIGN=CENTER NOWRAP>
249
<TR>

<TD class="math" ALIGN=CENTER NOWRAP>
1,800
<TD class="math" ALIGN=CENTER NOWRAP>
478
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
</TABLE>


\subsection SurfReconstPerfCont Contouring

The point set chosen for benchmarking the contouring stage is the Bimba con Nastrino point
set simplified to 100k points. We measure the contouring (i.e.\ the call to `make_surface_mesh()`)
duration and the reconstruction error for a range of approximation distances.
The reconstruction error is expressed as the average distance from input points to the reconstructed surface
in mm (the Bimba con Nastrino statue is 324 mm tall).

<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
Approx. distance (*average spacing)    
<TD class="math" ALIGN=CENTER NOWRAP>
Contouring duration (in s) 
<TD class="math" ALIGN=CENTER NOWRAP>
Reconstruction error (mm) 
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
<TR>

<TD class="math" ALIGN=CENTER NOWRAP>
0.1                                   
<TD class="math" ALIGN=CENTER NOWRAP>
19.2
<TD class="math" ALIGN=CENTER NOWRAP>
0.055
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.25
<TD class="math" ALIGN=CENTER NOWRAP>
6.9
<TD class="math" ALIGN=CENTER NOWRAP>
0.106
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
0.5
<TD class="math" ALIGN=CENTER NOWRAP>
3.2
<TD class="math" ALIGN=CENTER NOWRAP>
0.18
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1
<TD class="math" ALIGN=CENTER NOWRAP>
1.65
<TD class="math" ALIGN=CENTER NOWRAP>
0.36
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
2                                     
<TD class="math" ALIGN=CENTER NOWRAP>
0.8
<TD class="math" ALIGN=CENTER NOWRAP>
0.76
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=3><HR>
</TABLE>


\cgalFigureBegin{Poisson_surface_reconstruction_3-fig-contouring_bench,contouring_bench.jpg}
Contouring duration (in s) and reconstruction error (mm)
against several approximation distance parameters
for the Bimba con Nastrino point set simplified to 100k points.
\cgalFigureEnd

\subsection SurfReconstPerfMem Memory

We measure the memory occupancy for the reconstruction of the full Bimba con Nastrino point
set (1.8 millions points) as well as for simplified versions.\n
The Poisson implicit function computation has a memory peak when solving the Poisson linear 
system using the sparse linear solver.

<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
Number of points (x1000) 
<TD class="math" ALIGN=CENTER NOWRAP>
Memory occupancy (MBytes) 
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
180
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
100
<TD class="math" ALIGN=CENTER NOWRAP>
270
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
250
<TD class="math" ALIGN=CENTER NOWRAP>
790
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
<TD class="math" ALIGN=CENTER NOWRAP>
1300
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
<TD class="math" ALIGN=CENTER NOWRAP>
2200
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,800
<TD class="math" ALIGN=CENTER NOWRAP>
3800
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
</TABLE>


\subsection SurfReconstPerfPSS Point Set Simplification

Due to the memory limitations described above, we recommend to simplify the point sets captured by laser scanners.\n
We measure the reconstruction error for the Bimba con Nastrino point set (1.6M points) as well as for 
simplified versions. All reconstructions use the recommended contouring parameter 
`approximation distance = 0.25 * the input point` set's average spacing.
The reconstruction error is expressed as the average distance from input points to the reconstructed surface in mm 
(the Bimba con Nastrino statue is 324 mm tall).

<TABLE CELLSPACING=5 >
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
Number of points (x1000) 
<TD class="math" ALIGN=CENTER NOWRAP>
Reconstruction error (mm) 
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
60
<TD class="math" ALIGN=CENTER NOWRAP>
0.27
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
120                       
<TD class="math" ALIGN=CENTER NOWRAP>
0.15
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
250
<TD class="math" ALIGN=CENTER NOWRAP>
0.11
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
500
<TD class="math" ALIGN=CENTER NOWRAP>
0.079
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,000
<TD class="math" ALIGN=CENTER NOWRAP>
0.066
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,500
<TD class="math" ALIGN=CENTER NOWRAP>
0.061
<TR>
<TD class="math" ALIGN=CENTER NOWRAP>
1,600
<TD class="math" ALIGN=CENTER NOWRAP>
0.06
<TR><TD ALIGN=LEFT NOWRAP COLSPAN=2><HR>
</TABLE>


\cgalFigureBegin{Poisson_surface_reconstruction_3-fig-simplification_bench,simplification_bench.jpg}
Reconstruction error (mm) against number of points
for the Bimba con Nastrino point set with 1.6M points
as well as for simplified versions.
\cgalFigureEnd

\section SurfReconstDesignHistory Design and Implementation History

The initial implementation was essentially done by Laurent Saboret, guided by Pierre Alliez and Ga\"el Guennebaud. 
For later releases of the package Andreas Fabri worked on performance improvements, and Laurent Rineau added the
two passes for dealing with holes.


*/ 
} /* namespace CGAL */
