namespace CGAL {
/*!

\mainpage User Manual 
\anchor Chapter_Optimal_Transportation_Curve_Reconstruction
\cgalAutoToc

\authors Pierre Alliez, David Cohen-Steiner, Fernando de Goes, Clément Jamin and Ivo Vigan

\section Optimal_transportation_reconstruction_2Introduction Introduction

This package implements a method to reconstruct and simplify 2D point sets \cgalCite{degoes:hal-00758019}. The input is a set of 2D points with mass attributes, possibly hampered by noise and outliers. The output is a set of <em>line segments</em> and <em>isolated points</em> which approximate the input points, as illustrated in \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_summary}. The mass attribute relates to the importance given to each point for approximation.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_summary,summary.png}
Left: input point set hampered by noise. Right: The corresponding reconstructed shape consisting of line segments.
\cgalFigureEnd

Internally, the algorithm constructs an initial 2D Delaunay triangulation from all the input points, then simplifies the triangulation so that a subset of the edges and vertices of the triangulation approximate well the input points. Approximate herein refers to a robust distance based on <em>optimal transportation</em> (see section \ref Optimal_transportation_reconstruction_2HowDoesItWork for more details). The triangulation is simplified using a combination of <em>half edge collapse</em>, <em>edge flips</em> and <em>vertex relocation</em> operators. The triangulation remains valid during simplification, i.e., with neither overlaps nor fold-overs. 

The output of the reconstruction algorithm is a subset of edges and vertices of the triangulation. \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_algorithm} depicts an example where the output is composed of green edges and one isolated vertex. The green edges are considered relevant as they approximate well many of the input points. The edges depicted in grey, referred to as <em>ghost edges</em> and discarded, approximate none of the input points. The edges depicted in red, referred to as irrelevant and also discarded, approximate some of the input points but not enough to be considered relevant.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_algorithm,algorithm.png}
(a) Input points. (b) Delaunay triangulation of the input points. (c) After simplification, with ghost edges in grey, relevant edges in green, and irrelevant edges in red. (d) Final reconstruction made up of several edges and one isolated vertex.
\cgalFigureEnd


\subsection Optimal_transportation_reconstruction_2Simplest_example Simplest Examples

The following example first generates a set of input points on a
square. The points, with no mass attribute, are then passed to the
Optimal_transportation_reconstruction_2 object. After initialization,
100 iterations of the reconstruction process are performed.
\cgalExample{Optimal_transportation_reconstruction_2/otr2_simplest_example.cpp}

There is no direct relationship between the target number of edges for
the output reconstruction and a notion of approximation error. The
simplification algorithm can thus be stopped via a maximum tolerance
error homogeneous to a distance. More specifically, the tolerance
error is specified as the maximum square root of transport cost per
unit of mass, which is homogeneous to a distance.

Note that the tolerance is given in the sense of the Wasserstein
distance (see \ref Optimal_transportation_reconstruction_2Wasserstein). It is _not_ a
Hausdorff tolerance: it does not mean that the distance between the
input samples and the output polyline is guaranteed to be less than
`tolerance`. It means that the square root of transport cost per mass
(homogeneous to a distance) is at most `tolerance`.
\cgalExample{Optimal_transportation_reconstruction_2/otr2_simplest_example_with_tolerance.cpp}

\subsection Optimal_transportation_reconstruction_2Output_example Output Examples
The output of the reconstruction can be obtained in two ways: either as a sequence of 2D points and segments, or as an indexed format where the connectivity of the segments is encoded, hence the terms vertices and edges. The indexed format records a list of points, then pairs of point indices in the said list for the edges, and point indices for the isolated vertices. 
\cgalExample{Optimal_transportation_reconstruction_2/otr2_list_output_example.cpp}
\cgalExample{Optimal_transportation_reconstruction_2/otr2_indexed_output_example.cpp}

\subsection Optimal_transportation_reconstruction_2Mass_example Example with Mass Attributes
The following example first reads a set of input points and masses from an ASCII file. Using two property maps, the points and their initial mass are passed to the Optimal_transportation_reconstruction_2 object. After initialization 100 iterations of the reconstruction process are performed, then the segments and isolated points of the reconstructed shape are extracted and printed to the console. 
\cgalExample{Optimal_transportation_reconstruction_2/otr2_mass_example.cpp}


\section Optimal_transportation_reconstruction_2API API

The only class exposed to the user is the <I> Optimal_transportation_reconstruction_2 </I> class.

\subsection Optimal_transportation_reconstruction_2Sample Sample Call

\code{.cpp}
/*
K : a geometric kernel.
*/

Optimal_transportation_reconstruction_2<K>
  otr2(points.begin(), points.end());

otr2.run(100); // perform 100 edge collapse operators
\endcode

In case the input is not just points without masses, one can provide 
a property map that matches this input.
\code{.cpp}
/*
K                  : a geometric kernel.
Point_property_map : a PropertyMap for accessing the input points.
Mass_property_map  : a PropertyMap for accessing the mass attributes of the 
                     input points.
*/

Optimal_transportation_reconstruction_2<K, Point_property_map, Mass_property_map>
  otr2(points.begin(), points.end(), point_pmap, mass_pmap);

otr2.run(100); // perform 100 edge collapse operators
\endcode

Alternatively to calling \link Optimal_transportation_reconstruction_2::run `run()`\endlink, 
one can call \link Optimal_transportation_reconstruction_2::run_until `run_until()` \endlink
and specify the number of output vertices one wants to keep as illustrated in
 \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_twenty_vertices}.
\code{.cpp}
  otr2.run_until(20); // perform edge collapse operators until 20 vertices remain.
\endcode

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_twenty_vertices,twenty_vertices.png}
Examples of 20-vertex reconstructions from datasets consisting of 2000, 400 
and 200 input points respectively. These examples illustrate the behavior of
the algorithm when the input point density decreases. 
\cgalFigureEnd

Finally, calling \link Optimal_transportation_reconstruction_2::run_under_wasserstein_tolerance `run_under_wasserstein_tolerance()` \endlink
allows the user to run the algorithm based on a distance criterion
when the number of output edges is not known as illustrated in
 \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_tolerance}.
\code{.cpp}
otr2.run_under_wasserstein_tolerance(0.1); // perform collapses until no more can be done within tolerance
\endcode

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_tolerance,tolerances.png}
Examples of reconstructions with different Wasserstein tolerance threshold. Top: input point set and reconstruction with a tolerance of 0.005. Bottom: reconstructions with a tolerance of 0.05 and with a tolerance of 0.1.
\cgalFigureEnd


\subsection Optimal_transportation_reconstruction_2Global_relocation Global Point Relocation

As noise and missing data may prevent the reconstructed shape to have sharp
corners at the correct places, the algorithm offers a function to relocate 
all points of the reconstruction:
\code{.cpp}
  otr2.relocate_all_points();
\endcode
Note that these points coincide with the vertices of the underlying triangulation. This function can be called either after one run of simplification, or interleaved with several runs of simplification.

The new point locations are chosen such that the approximation of the output segments and isolated points to the input points is improved. More specifically, the relocation process iterates between computing the best transport plan given the current reconstruction, and relocating the triangulation vertices while keeping the current transport plan unchanged. The vertices are relocated so as to minimize the transport cost induced by the current transport plan \cgalCite{degoes:hal-00758019}.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_relocation,relocation.png}
Left: before point relocation. Right: after point relocation.
\cgalFigureEnd


\section Optimal_transportation_reconstruction_2Parameters User Parameters

The behavior of the algorithm is controlled via the following parameters. 

\subsection Optimal_transportation_reconstruction_2Flip Edge Flipping

During simplification of the internal triangulation some recursive edge flip operators are required to guarantee that the triangulation remain valid when applying a half edge collapse operator (see \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_edgeflip}). Calling set_use_flip(false) prevents the algorithm from using edge flips, yielding shorter computational times at the price of suboptimal results as not all edges can be considered for being collapsible.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_edgeflip,edgeflip.jpg}
Edge flipping. Left: the blue edge creates fold-overs because of blocking edges shown in black. Middle: after running the recursive edge flipping procedure the blue edge is collapsible. Right: triangulation after edge collapse.
\cgalFigureEnd

\subsection Optimal_transportation_reconstruction_2Relevance Edge Relevance

An edge is relevant from the approximation point of view if (1) it is long, (2) approximates a large number of points (or a large amount of mass when points have mass attributes), and (3) has a small approximation error. More specifically, the notion of relevance is defined as \f$ m(e) * |e|^2 / cost(e) \f$, where \f$ m(e) \f$ denotes the mass of the points approximated by the edge, \f$ |e| \f$ denotes the edge length and \f$ cost(e) \f$ its approximation error. As the error is defined by mass time squared distance the relevance is unitless.
The default value is 1, so that all edges approximating some input points are considered relevant. A larger relevance value provides us with a means to increase resilience to outliers.


\subsection Optimal_transportation_reconstruction_2Random Random Sample Size

By default the simplification relies upon an exhaustive priority queue of half edge collapse operators during decimation. For improved efficiency, a parameter <I> sample size </I> strictly greater than 0 switches to a multiple choice approach, i.e., a best-choice selection in a random sample of edge collapse operators, of size <I>sample size</I>. A typical value for the sample size is 15, but this value must be enlarged when targeting a very coarse simplification. 

\subsection Optimal_transportation_reconstruction_2Local_relocation Local Point Relocation

In addition to the global relocation function described above, an optional parameter of the constructor of the <I> Optimal_transportation_reconstruction_2 </I> class provides a means to relocate the points <em>locally</em> after each edge collapse operator (possibly combined with edge flips). Locally herein means that only the vertices of a local stencil in the triangulation around the each edge collapse operator are relocated, with a process similar to the one described above in the global relocation function. The local stencil is chosen as the one-ring neighborhood of the vertex remaining after collapsing an edge. The relocation process being iterative, one parameter controls the number of relocation steps. 


\subsection Optimal_transportation_reconstruction_2Verbose Verbose Output

The verbose parameter, between 0 and 2, determines how much console output the algorithm generates. A 0 value generates no output to the standard output. A value greater than 0 generates output to the standard output `std::cerr`.


\section Optimal_transportation_reconstruction_2Robustness Robustness
A virtue of the algorithm is its robustness to noise and outliers. 
\cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_robustness} shows that
the output of the algorithm is robust in the sense that it is hardly affected by noise and/or outliers, as long as the density of outliers is small compared to the density of the input points.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_robustness,robustness.png}
Robustness to noise and outliers. Left: noise-free point set. Middle: noisy point set. Right: point set hampered by noise and outliers.
\cgalFigureEnd

\section Optimal_transportation_reconstruction_2Density Variable Density
\cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_density} illustrates the behavior of the algorithm on a point set with uniform mass attributes, versus variable density. As the algorithm gives more importance to densely sampled areas, this translates into smaller edges on densely sampled areas. On sparsely sampled areas the algorithm initially approximates each point by one isolated vertex, then progressively approximates the points with edges. 
\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_density,density.png}
Variable density. Left: input point set. The three other pictures show
how the approximation evolves when pushing the simplification forward.
\cgalFigureEnd

\section Optimal_transportation_reconstruction_2Mixed Mixed Dimension
\cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_mixed} depicts an input point set sampling a set of line segments and a solid area. Depending on the targeted number of points in the output, the solid area is approximated by a set of evenly sampled isolated vertices.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_mixed,mixed.png}
Mixed dimension. Left: input point set. Middle: Isolated vertices in blue, relevant edges in green and irrelevant edges in red. Right: final output.
\cgalFigureEnd

\section Optimal_transportation_reconstruction_2Variable_mass Variable Mass

The mass attributes provides a means to adjust the importance given to each point for approximation. \cgalFigureRef{Optimal_Transportation_Curve_Reconstruction_variable} depicts a reconstruction from a gray level image after thresholding, where the gray level of the pixels are used as mass attribute.

\cgalFigureBegin{Optimal_Transportation_Curve_Reconstruction_variable,variable.png}
Variable mass. Left: input gray level image. Middle: image after thresholding to reduce the number of pixels used as points with non-zero mass. Right: final reconstruction.
\cgalFigureEnd

\section Optimal_transportation_reconstruction_2HowDoesItWork How Does It Work?

The task addressed here is to reconstruct a shape from a noisy point set \f$ S \f$ in \f$ \mathbb{R}^2 \f$, i.e., given a set of points in the plane, find a set of points and segments (more formally, a <I> 0-1 simplicial complex </I>) which best approximates \f$ S \f$. 

The approximation error is derived from the theory of optimal transportation between geometric measures \cgalCite{degoes:hal-00758019}. More specifically, the input point set is seen as a discrete measure, i.e., a set of pointwise masses. The goal is to find a 0-1 simplicial complex where the edges are the support of a <em>piecewise uniform</em> continuous measure (i.e., a line density of masses) and the vertices are the support of a discrete measure. Approximating the input point set in our context translates into approximating the input discrete measure by another measure composed of line segments and points. 

\subsection Optimal_transportation_reconstruction_2Wasserstein Wasserstein Distance 

Intuitively, the optimal transportation distance (Wasserstein-2 distance in our context) measures the amount of work that it takes to transport the input measure onto the vertices and edges of the triangulation, where the measure is constrained to be uniform (and greater or equal to zero) on each edge, and just greater or equal to zero on each vertex. Note that the Wasserstein distance is symmetric.

When all vertices of the triangulation coincide with the input points (after full initialization) the total transport cost is zero as each input point relocates to a vertex at no cost, and the reconstruction is made of isolated vertices only.

Assume for now that the input point set is composed of 10 points (each with mass 1) uniformly sampled on a line segment, and that the triangulation contains a single edge coinciding with the line segment. Although the (one-sided Euclidean) distance from the points to the edge is zero (the converse being not zero), the Wasserstein distance from the points to the edge is non zero, as we constrain the mass density of the edge to be uniform, and the total mass (integral of density) of the edge to equal 10, i.e., the total mass of the input points. The input points should hence be transported <em> tangentially </em> on the edge in order to match the uniform density, the optimal <em> transport plan </em> for the input points being described as smaller line segments with equal lengths, covering the edge.

If one now samples uniformly 20 points on the same edge (each with mass 0.5) the Wasserstein distance is smaller (albeit the total mass is 10 as before), as the transport plan is described by smaller line segments. 
In a slightly different configuration where the 20 input points have different masses, the optimal transport plan is described by small line segments whose length is proportional to the mass of the associated input point.
When the input points are not strictly on the edge, the transport plan has both tangential and normal components. 

In other words, one approximates well the input points by a single edge when the input points are densely and uniformly sampling that edge. Beyond symmetry, one virtue of the Wasserstein distance is thus to quantize both the deviation from the points to the edge, and the non-uniformity of the points on that edge. This distance is also resilient to outliers (points far from the edge) when the mass of these outliers is small compared to the total mass of the input points.

\subsection Optimal_transportation_reconstruction_2Reconstruction Reconstruction

The algorithm performs a fine to coarse simplification of a triangulation. It starts by constructing a box around the input points \f$S\f$ and computes the Delaunay Triangulation \f$ T_0 \f$ on all or a subset of \f$ S \f$. \f$ T_0 \f$ is the first output simplex, which is simplified in subsequent iterations by repeated edge collapses. To chose the next edge, an edge collapse operator is simulated for all feasible edges, i.e., for edges which introduce neither overlaps nor fold-overs in the triangulation. The next edge \f$ e \f$ for collapse is chosen according to the overall cost of the transportation plan for \f$ T \setminus e \f$, where the cheapest overall cost is preferred. Since disregarding edges which do not preserve the embedding of the triangulation can severely affect the performance of the greedy approach to optimal transport, the collapse operator is modified by adding a local flip procedure which makes every edge collapsible.

The transportation plan is approximated by assigning each input point temporarily to the closest simplex edge. After this partitioning of the input points w.r.t. the edges, all the points temporarily assigned to a given edge are being assigned to it permanently if and only if the corresponding transportation costs are less than the transportation cost for each of the two end vertices of the edge. Otherwise each of the points is assigned to the cheapest of the two vertices. This process of edge collapse and transportation plan update is repeated until the desired number of vertices, specified by the users, has been reached. 

At the end of the process, the edges which carry little mass can be filtered out and the remaining relevant edges and isolated vertices are reported as reconstructing the input shape. 

\section Optimal_transportation_reconstruction_2History Design and Implementation History
This package is the result of the work of Ivo Vigan during the 2014 season
of the Google Summer of Code. He has been mentored by Clément Jamin, 
Pierre Alliez and Fernando de Goes. The code is based on an initial research 
code written by Pierre Alliez and Fernando de Goes in 2010, for a paper published 
at the EUROGRAPHICS Symposium on Geometry Processing 2011 \cgalCite{degoes:hal-00758019}.

*/ 
} /* namespace CGAL */

