Merge pull request #367 from sgiraudot/Point_set_processing-hierarchical_clustering-sgiraudot

New point set processing algorithm: Hierarchical clustering

Documentation/biblio/cgal_manual.bib

@@ -1550,6 +1550,15 @@ ABSTRACT = {We present the first complete, exact and efficient C++ implementatio
 }
 
 
+@inproceedings{cgal:pgk-esops-02,
+  title={Efficient simplification of point-sampled surfaces},
+  author={Pauly, Mark and Gross, Markus and Kobbelt, Leif P},
+  booktitle={Proceedings of the conference on Visualization'02},
+  pages={163--170},
+  year={2002},
+  organization={IEEE Computer Society}
+}
+
 @article{ cgal:pp-cdmsc-93,
    author  = "U. Pinkall and K. Polthier",
    title   = "Computing discrete minimal surfaces and their conjugates",

Installation/changes.html

@@ -154,6 +154,13 @@ and <code>src/</code> directories).
     or <code>CGAL::Parallel_tag</code> when calling one of these functions.</li>
     <li> <code>CGAL::Parallel_tag</code> can no longer be used in
     Point Set Processing algorithms if TBB is not available.</li>
+    <li>
+    Add a new simplification algorithm based on hierarchical
+    clustering: <code>CGAL::hierarchy_simplify_point_set()</code>. It
+    allows either to uniformly simplify the point set or to
+    automatically adapt the local density of points to the local
+    variation of the input computed by principal component analysis.
+    </li>
   </ul>
   <h3>Surface Mesh Parameterization</h3>
   <ul>

Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt

@@ -30,6 +30,7 @@
 - `CGAL::remove_outliers()`
 - `CGAL::grid_simplify_point_set()`
 - `CGAL::random_simplify_point_set()`
+- `CGAL::hierarchy_simplify_point_set()`	
 - `CGAL::wlop_simplify_and_regularize_point_set()`
 - `CGAL::jet_smooth_point_set()`
 - `CGAL::bilateral_smooth_point_set()`

Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt

@@ -138,7 +138,7 @@ functions in this component.)
 \section Point_set_processing_3Simplification Simplification
 
 
-Three simplification functions are devised to reduce an input point set.
+Four simplification functions are devised to reduce an input point set.
 
 Function `random_simplify_point_set()` randomly deletes a
 user-specified fraction of points from the input point set. This
@@ -150,6 +150,12 @@ points sharing the same cell of the grid by picking as representant
 one arbitrarily chosen point. This algorithm is slower than
 `random_simplify_point_set()`.
 
+Function `hierarchy_simplify_point_set()` provides an adaptive
+simplification of the point set through local clusters
+\cgalCite{cgal:pgk-esops-02}. The size of the clusters is either
+directly selected by the user or it automatically adapts to the local
+variation of the point set.
+
 Function `wlop_simplify_and_regularize_point_set()` not only simplifies, 
 but also regularizes downsampled points. This is an implementation of
 the Weighted Locally Optimal Projection (WLOP) algorithm \cgalCite{wlop-2009}.
@@ -164,6 +170,40 @@ The following example reads a point set and simplifies it by clustering.
 Point set simplification through grid-based clustering. Removed points are depicted in red. Notice how low-density areas (in green) are not simplified.
 \cgalFigureEnd
 
+\subsection Point_set_processing_3Example_9 Hierarchy Simplification Example
+The following example reads a point set and produces a set of clusters.
+
+\cgalExample{Point_set_processing_3/hierarchy_simplification_example.cpp}
+
+\subsubsection Point_set_processing_3Hierarchy_simplification_parameter_size Parameter: size
+The hierarchy simplification algorithm recursively split the point set
+in two until each cluster's size is less than the parameter `size`.
+
+\cgalFigureBegin{Point_set_processing_3figHierarchy_simplification_size, hierarchical_clustering_size.jpg}
+Input point set and hierarchy simplification with different `size`
+parameter: \f$10\f$, \f$100\f$ and \f$1000\f$. In the 3 cases,
+`var_max`\f$=1/3\f$.  \cgalFigureEnd
+
+
+\subsubsection Point_set_processing_3Hierarchy_simplification_parameter_var_max Parameter: var_max
+In addition to the size parameter, a variation parameter allows to
+increase simplification in monotoneous regions. For each cluster, a
+surface variation measure is computed using the sorted eigenvalues of
+the covariance matrix: \f[ \sigma(p) = \frac{\lambda_0}{\lambda_0 +
+    \lambda_1 + \lambda_2}. \f]
+
+This function goes from \f$0\f$ if the cluster is coplanar to
+\f$1/3\f$ if it is fully isotropic. If a cluster's variation is above
+`var_max`, it is splitted. If `var_max` is equal to \f$1/3\f$, this
+parameter has no effect and the clustering is regular on the whole
+point set.
+
+\cgalFigureBegin{Point_set_processing_3figHierarchical_clustering_var_max, hierarchical_clustering_var_max.jpg}
+Input point set and hierarchy simplification with different `var_max`
+parameter: \f$0.00001\f$, \f$0.001\f$ and \f$0.1\f$. In the 3 cases,
+`size`\f$=1000\f$.  \cgalFigureEnd
+
+
 \subsection Point_set_processing_3Example_4 WLOP Simplification Example
 The following example reads a point set, simplifies and regularizes it by WLOP.
 
@@ -199,6 +239,7 @@ for more details. We provide below a speed-up chart generated using the parallel
 Parallel WLOP speed-up, compared to the sequential version of the algorithm. 
 \cgalFigureEnd
 
+
 \section Point_set_processing_3Smoothing Smoothing
 
 Two smoothing functions are devised to smooth an input point set.

Point_set_processing_3/doc/Point_set_processing_3/examples.txt

@@ -4,6 +4,7 @@
 \example Point_set_processing_3/remove_outliers_example.cpp
 \example Point_set_processing_3/grid_simplification_example.cpp
 \example Point_set_processing_3/grid_simplify_indices.cpp
+\example Point_set_processing_3/hierarchy_simplification_example.cpp	
 \example Point_set_processing_3/jet_smoothing_example.cpp
 \example Point_set_processing_3/normals_example.cpp
 \example Point_set_processing_3/wlop_simplify_and_regularize_point_set_example.cpp

Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt

@@ -57,6 +57,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "bilateral_smooth_point_set_example.cpp" )
   create_single_source_cgal_program( "grid_simplification_example.cpp" )
   create_single_source_cgal_program( "grid_simplify_indices.cpp" )
+  create_single_source_cgal_program( "hierarchy_simplification_example.cpp" )
   create_single_source_cgal_program( "normals_example.cpp" )
   create_single_source_cgal_program( "property_map.cpp" )
   create_single_source_cgal_program( "random_simplification_example.cpp" )

Point_set_processing_3/examples/Point_set_processing_3/hierarchy_simplification_example.cpp

@@ -0,0 +1,48 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/hierarchy_simplify_point_set.h>
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/IO/write_xyz_points.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Memory_sizer.h>
+
+#include <vector>
+#include <fstream>
+
+// types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3 Point;
+
+int main(int argc, char*argv[])
+{
+  // Reads a .xyz point set file in points[].
+  std::vector<Point> points;
+  const char* fname = (argc>1)?argv[1]:"data/oni.xyz";
+  std::ifstream stream(fname);
+  if (!stream ||
+      !CGAL::read_xyz_points(stream, std::back_inserter(points)))
+  {
+    std::cerr << "Error: cannot read file " << fname << std::endl;
+    return EXIT_FAILURE;
+  }
+  std::cout << "Read " << points.size () << " point(s)" << std::endl;
+
+  CGAL::Timer task_timer; task_timer.start();
+
+  // simplification by clustering using erase-remove idiom
+  points.erase (CGAL::hierarchy_simplify_point_set (points.begin (), points.end (),
+						    100, // Max cluster size
+						    0.01), // Max surface variation
+		points.end ());
+
+  std::size_t memory = CGAL::Memory_sizer().virtual_size();
+  
+  std::cout << points.size () << " point(s) kept, computed in "
+	    << task_timer.time() << " seconds, "
+	    << (memory>>20) << " Mib allocated." << std::endl;
+
+  std::ofstream f ("out.xyz");
+  CGAL::write_xyz_points (f, points.begin (), points.end ());
+  
+  return EXIT_SUCCESS;
+}
+

Point_set_processing_3/include/CGAL/hierarchy_simplify_point_set.h

@@ -0,0 +1,380 @@
+// Copyright (c) 2012  INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:
+// $Id: 
+//
+// Author(s) : Simon Giraudot, Pierre Alliez
+
+
+#ifndef HIERARCHY_SIMPLIFY_POINT_SET_H
+#define HIERARCHY_SIMPLIFY_POINT_SET_H
+
+#include <cmath>
+#include <stack>
+
+#include <CGAL/property_map.h>
+#include <CGAL/basic.h>
+#include <CGAL/Dimension.h>
+#include <CGAL/Object.h>
+#include <CGAL/centroid.h>
+#include <CGAL/point_set_processing_assertions.h>
+#include <CGAL/Default_diagonalize_traits.h>
+#include <CGAL/PCA_util.h>
+
+namespace CGAL {
+
+
+  namespace internal {
+
+    template < typename InputIterator,
+	       typename PointPMap,
+	       typename K >
+    typename K::Point_3
+    hsps_centroid(InputIterator begin, 
+		 InputIterator end,
+		 PointPMap& point_pmap,	     
+		 const K&)
+    {
+      typedef typename K::Point_3 Point;
+      typedef typename K::FT FT;
+
+      CGAL_precondition(begin != end);
+
+      FT x = (FT)0., y = (FT)0., z = (FT)0.;
+      unsigned int nb_pts = 0;
+      while(begin != end) 
+	{
+#ifdef CGAL_USE_PROPERTY_MAPS_API_V1
+	  const Point& point = get(point_pmap, begin);
+#else
+	  const Point& point = get(point_pmap, *begin);
+#endif
+	  x += point.x ();  y += point.y ();  z += point.z ();
+	  ++ nb_pts;
+	  ++ begin;
+	}
+      return Point (x/nb_pts, y/nb_pts, z/nb_pts);
+    }
+
+    template < typename Input_type,
+	       typename PointPMap,
+	       typename K >
+    void
+    hsc_terminate_cluster (std::list<Input_type>& cluster,
+			   std::list<Input_type>& points_to_keep,
+			   std::list<Input_type>& points_to_remove,
+			   PointPMap& point_pmap,
+			   const typename K::Point_3& centroid,
+			   const K&)
+    {
+      typedef typename std::list<Input_type>::iterator Iterator;
+      typedef typename K::FT FT;
+      typedef typename K::Point_3 Point;
+
+      FT dist_min = (std::numeric_limits<FT>::max)();
+
+      typename std::list<Input_type>::iterator point_min;
+      for (Iterator it = cluster.begin (); it != cluster.end (); ++ it)
+	{
+#ifdef CGAL_USE_PROPERTY_MAPS_API_V1
+	  const Point& point = get(point_pmap, it);
+#else
+	  const Point& point = get(point_pmap, *it);
+#endif
+	  FT dist = CGAL::squared_distance (point, centroid);
+	  if (dist < dist_min)
+	    {
+	      dist_min = dist;
+	      point_min = it;
+	    }
+	}
+
+      points_to_keep.splice (points_to_keep.end (), cluster, point_min);
+      points_to_remove.splice (points_to_remove.end (), cluster, cluster.begin (), cluster.end ());
+    }			       
+			    
+
+
+
+  } // namespace internal
+    
+
+  /// \ingroup PkgPointSetProcessing
+  
+  /// Recursively split the point set in smaller clusters until the
+  /// clusters have less than `size` elements or until their variation
+  /// factor is below `var_max`. 
+  ///
+  /// This method modifies the order of input points so as to pack all remaining points first,
+  /// and returns an iterator over the first point to remove (see erase-remove idiom).
+  /// For this reason it should not be called on sorted containers.
+  ///
+  /// \pre `0 < var_max < 1/3`
+  /// \pre `size > 0`
+  ///
+  /// @tparam ForwardIterator iterator over input points.
+  /// @tparam PointPMap is a model of `ReadablePropertyMap` with value type `Point_3<Kernel>`.
+  ///        It can be omitted if the value type of `ForwardIterator` is convertible to `Point_3<Kernel>`.
+  /// @tparam DiagonalizeTraits is a model of `DiagonalizeTraits`. It
+  ///        can be omitted: if Eigen 3 (or greater) is available and
+  ///        `CGAL_EIGEN3_ENABLED` is defined then an overload using
+  ///        `Eigen_diagonalize_traits` is provided. Otherwise, the internal
+  ///        implementation `Internal_diagonalize_traits` is used.
+  /// @tparam Kernel Geometric traits class.
+  ///        It can be omitted and deduced automatically from the value type of `PointPMap`.
+  ///
+  /// @return iterator over the first point to remove.
+
+
+  // This variant requires all parameters.
+
+  template <typename ForwardIterator,
+	    typename PointPMap,
+	    typename DiagonalizeTraits,
+	    typename Kernel>
+  ForwardIterator hierarchy_simplify_point_set (ForwardIterator begin,
+						ForwardIterator end,
+						PointPMap point_pmap,
+						const unsigned int size,
+						const double var_max,
+						const DiagonalizeTraits&,
+						const Kernel&)
+  {
+    typedef typename std::iterator_traits<ForwardIterator>::value_type Input_type;
+    typedef typename Kernel::FT FT;
+    typedef typename Kernel::Point_3  Point;
+    typedef typename Kernel::Vector_3 Vector;
+
+    // We define a cluster as a point set + its centroid (useful for
+    // faster computations of centroids - to be implemented)
+    typedef std::pair< std::list<Input_type>, Point > cluster;
+
+    std::list<cluster> clusters_stack;
+    typedef typename std::list<cluster>::iterator cluster_iterator;
+
+    CGAL_precondition (begin != end);
+    CGAL_point_set_processing_precondition (size > 0);
+    CGAL_point_set_processing_precondition (var_max > 0.0);
+
+    // The first cluster is the whole input point set
+    clusters_stack.push_front (cluster (std::list<Input_type>(), Point (0., 0., 0.)));
+    std::copy (begin, end, std::back_inserter (clusters_stack.front ().first));
+    
+    clusters_stack.front ().second = internal::hsps_centroid (clusters_stack.front ().first.begin (),
+							     clusters_stack.front ().first.end (),
+							     point_pmap, Kernel());
+
+    std::list<Input_type> points_to_keep;
+    std::list<Input_type> points_to_remove;
+    
+    while (!(clusters_stack.empty ()))
+      {
+	cluster_iterator current_cluster = clusters_stack.begin ();
+
+	// If the cluster only has 1 element, we add it to the list of
+	// output points
+	if (current_cluster->first.size () == 1)
+	  {
+	    points_to_keep.splice (points_to_keep.end (), current_cluster->first,
+				   current_cluster->first.begin ());
+	    clusters_stack.pop_front ();
+	    continue;
+	  }
+
+	// Compute the covariance matrix of the set
+	cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
+
+	for (typename std::list<Input_type>::iterator it = current_cluster->first.begin ();
+	     it != current_cluster->first.end (); ++ it)
+	  {
+#ifdef CGAL_USE_PROPERTY_MAPS_API_V1
+	    const Point& point = get(point_pmap, it);
+#else
+	    const Point& point = get(point_pmap, *it);
+#endif
+	    Vector d = point - current_cluster->second;
+	    covariance[0] += d.x () * d.x ();
+	    covariance[1] += d.x () * d.y ();
+	    covariance[2] += d.x () * d.z ();
+	    covariance[3] += d.y () * d.y ();
+	    covariance[4] += d.y () * d.z ();
+	    covariance[5] += d.z () * d.z ();
+	  }
+
+	cpp11::array<FT, 3> eigenvalues = {{ 0., 0., 0. }};
+	cpp11::array<FT, 9> eigenvectors = {{ 0., 0., 0.,
+					      0., 0., 0.,
+					      0., 0., 0. }};
+	// Linear algebra = get eigenvalues and eigenvectors for
+	// PCA-like analysis
+	DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
+	  (covariance, eigenvalues, eigenvectors);
+	
+	// Variation of the set defined as lambda_min / (lambda_0 + lambda_1 + lambda_2)
+	double var = eigenvalues[0] / (eigenvalues[0] + eigenvalues[1] + eigenvalues[2]);
+
+	// Split the set if size OR variance of the cluster is too large
+	if (current_cluster->first.size () > size || var > var_max)
+	  {
+	    clusters_stack.push_front (cluster (std::list<Input_type>(), Point (0., 0., 0.)));
+	    cluster_iterator negative_side = clusters_stack.begin ();
+	    // positive_side is built directly from current_cluster
+	    
+	    // The plane which splits the point set into 2 point sets:
+	    //  * Normal to the eigenvector with highest eigenvalue
+	    //  * Passes through the centroid of the set
+	    Vector v (eigenvectors[6], eigenvectors[7], eigenvectors[8]);
+
+	    std::size_t current_cluster_size = 0;
+	    typename std::list<Input_type>::iterator it = current_cluster->first.begin ();
+	    while (it != current_cluster->first.end ())
+	      {
+		typename std::list<Input_type>::iterator current = it ++;
+
+		// Test if point is on negative side of plane and
+		// transfer it to the negative_side cluster if it is
+		if (Vector (current_cluster->second, *current) * v < 0)
+		  negative_side->first.splice (negative_side->first.end (),
+					       current_cluster->first, current);
+		++ current_cluster_size;
+	      }
+
+	    // If one of the clusters is empty, stop to avoid infinite
+	    // loop and keep the non-empty one
+	    if (current_cluster->first.empty () || negative_side->first.empty ())
+	      {
+		cluster_iterator nonempty = (current_cluster->first.empty ()
+					     ? negative_side : current_cluster);
+
+		// Compute the centroid
+		nonempty->second = internal::hsps_centroid (nonempty->first.begin (),
+							    nonempty->first.end (),
+							    point_pmap, Kernel());
+		
+		internal::hsc_terminate_cluster (nonempty->first,
+						 points_to_keep,
+						 points_to_remove,
+						 point_pmap,
+						 nonempty->second,
+						 Kernel ());
+		clusters_stack.pop_front ();
+		clusters_stack.pop_front ();
+	      }
+	    else
+	      {
+		// Save old centroid for faster computation
+		Point old_centroid = current_cluster->second;
+		
+		// Compute the first centroid
+		current_cluster->second = internal::hsps_centroid (current_cluster->first.begin (),
+								   current_cluster->first.end (),
+								   point_pmap, Kernel());
+
+		// The second centroid can be computed with the first and
+		// the old ones :
+		// centroid_neg = (n_total * old_centroid - n_pos * first_centroid)
+		//                 / n_neg;
+		negative_side->second = Point ((current_cluster_size * old_centroid.x ()
+						- current_cluster->first.size () * current_cluster->second.x ())
+					       / negative_side->first.size (),
+					       (current_cluster_size * old_centroid.y ()
+						- current_cluster->first.size () * current_cluster->second.y ())
+					       / negative_side->first.size (),
+					       (current_cluster_size * old_centroid.z ()
+						- current_cluster->first.size () * current_cluster->second.z ())
+					       / negative_side->first.size ());
+	      }
+	  }
+	// If the size/variance are small enough, add the centroid as
+	// and output point
+	else
+	  {
+	    internal::hsc_terminate_cluster (current_cluster->first,
+					     points_to_keep,
+					     points_to_remove,
+					     point_pmap,
+					     current_cluster->second,
+					     Kernel ());
+	    clusters_stack.pop_front ();
+	  }
+      }
+    ForwardIterator first_point_to_remove =
+      std::copy (points_to_keep.begin(), points_to_keep.end(), begin);
+    std::copy (points_to_remove.begin(), points_to_remove.end(), first_point_to_remove);
+
+    return first_point_to_remove;
+
+  }
+  /// @endcond
+
+
+  /// @cond SKIP_IN_MANUAL
+  // This variant deduces the kernel from the iterator type.
+  template <typename ForwardIterator,
+	    typename PointPMap,
+	    typename DiagonalizeTraits>
+  ForwardIterator hierarchy_simplify_point_set (ForwardIterator begin,
+						ForwardIterator end,
+						PointPMap point_pmap,
+						const unsigned int size,
+						const double var_max,
+						const DiagonalizeTraits& diagonalize_traits)
+  {
+    typedef typename boost::property_traits<PointPMap>::value_type Point;
+    typedef typename Kernel_traits<Point>::Kernel Kernel;
+    return hierarchy_simplify_point_set (begin, end, point_pmap, size, var_max,
+					 diagonalize_traits, Kernel());
+  }
+  /// @endcond
+
+  /// @cond SKIP_IN_MANUAL  
+  // This variant uses default diagonalize traits
+  template <typename ForwardIterator,
+	    typename PointPMap >
+  ForwardIterator hierarchy_simplify_point_set (ForwardIterator begin,
+						ForwardIterator end,
+						PointPMap point_pmap,
+						const unsigned int size,
+						const double var_max)
+  {
+    typedef typename boost::property_traits<PointPMap>::value_type Point;
+    typedef typename Kernel_traits<Point>::Kernel Kernel;
+    return hierarchy_simplify_point_set (begin, end, point_pmap, size, var_max,
+					 Default_diagonalize_traits<double, 3> (), Kernel());
+  }
+  /// @endcond  
+
+  /// @cond SKIP_IN_MANUAL
+  // This variant creates a default point property map = Identity_property_map.
+  template <typename ForwardIterator >
+  ForwardIterator hierarchy_simplify_point_set (ForwardIterator begin,
+						ForwardIterator end,
+						const unsigned int size = 10,
+						const double var_max = 0.333)
+  {
+    return hierarchy_simplify_point_set
+      (begin, end,
+#ifdef CGAL_USE_PROPERTY_MAPS_API_V1
+       make_dereference_property_map(first),
+#else
+       make_identity_property_map (typename std::iterator_traits<ForwardIterator>::value_type()),
+#endif
+       size, var_max);
+  }
+  /// @endcond  
+
+} // namespace CGAL
+
+#endif // HIERARCHY_SIMPLIFY_POINT_SET_H

Point_set_processing_3/test/Point_set_processing_3/CMakeLists.txt

@@ -71,6 +71,7 @@ if ( CGAL_FOUND )
   if(EIGEN3_FOUND OR LAPACK_FOUND)
     # Executables that require Eigen or BLAS and LAPACK
     create_single_source_cgal_program( "normal_estimation_test.cpp" )
+    create_single_source_cgal_program( "hierarchy_simplification_test.cpp" )
     create_single_source_cgal_program( "smoothing_test.cpp" )
     create_single_source_cgal_program( "vcm_plane_test.cpp" )
     create_single_source_cgal_program( "vcm_all_test.cpp" )

Point_set_processing_3/test/Point_set_processing_3/hierarchy_simplification_test.cpp

@@ -0,0 +1,90 @@
+#include <limits>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/hierarchy_simplify_point_set.h>
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/IO/write_xyz_points.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Memory_sizer.h>
+
+#include <vector>
+#include <fstream>
+
+// types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::FT FT;
+
+void test (std::vector<Point>& input,
+	   int result0 = 1, int result1 = 1, int result2 = 1, int result3 = 1, int result4 = 1)
+{
+  std::vector<Point>::iterator it = 
+    CGAL::hierarchy_simplify_point_set (input.begin (), input.end (), 1);
+  if (result0 > 0 && std::distance (input.begin (), it) != (result0))
+    exit (EXIT_FAILURE);
+
+  it = CGAL::hierarchy_simplify_point_set (input.begin (), input.end ());
+  if (result1 > 0 && std::distance (input.begin (), it) != (result1))
+    exit (EXIT_FAILURE);
+
+  it = CGAL::hierarchy_simplify_point_set (input.begin (), input.end (), 100);
+  if (result2 > 0 && std::distance (input.begin (), it) != (result2))
+    exit (EXIT_FAILURE);
+
+
+  it = CGAL::hierarchy_simplify_point_set (input.begin (), input.end (), 1000, 0.1);
+  if (result3 > 0 && std::distance (input.begin (), it) != (result3))
+    exit (EXIT_FAILURE);
+
+
+  it = CGAL::hierarchy_simplify_point_set (input.begin (), input.end (),
+					   CGAL::Identity_property_map<Point>(),
+					   (std::numeric_limits<unsigned int>::max)(),
+					   0.0001);
+  if (result4 > 0 && std::distance (input.begin (), it) != (result4))
+    exit (EXIT_FAILURE);
+
+  input.clear ();
+}
+
+
+int main(void)
+{
+
+  std::vector<Point> input;
+
+  // Test 1 point
+  input.push_back (Point (0., 0., 0.));
+  test (input);
+
+  // Test twice the same point
+  input.push_back (Point (0., 0., 0.));
+  input.push_back (Point (0., 0., 0.));
+  test (input);
+
+  // Test 2 points
+  input.push_back (Point (0., 0., 0.));
+  input.push_back (Point (1., 0., 0.));
+  test (input, 2);
+
+  // Test line
+  for (std::size_t i = 0; i < 1000; ++ i)
+    input.push_back (Point (0., 0., i));
+  test (input, input.size (), 128, 16, 1, 1);
+  
+  // Test plane
+  for (std::size_t i = 0; i < 128; ++ i)
+    for (std::size_t j = 0; j < 128; ++ j)
+      input.push_back (Point (0., j, i));
+  test (input, input.size (), 2048, 256, 32, 1);
+  
+  // Test random
+  for (std::size_t i = 0; i < 10000; ++ i)
+    input.push_back (Point (rand() / (FT)RAND_MAX,
+  			    rand() / (FT)RAND_MAX,
+  			    rand() / (FT)RAND_MAX));
+  test (input, input.size (), -1, -1, -1, -1);
+  
+  return EXIT_SUCCESS;
+}
+

Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_simplification_plugin.cpp

@@ -5,6 +5,7 @@
 
 #include <CGAL/grid_simplify_point_set.h>
 #include <CGAL/random_simplify_point_set.h>
+#include <CGAL/hierarchy_simplify_point_set.h>
 #include <CGAL/compute_average_spacing.h>
 #include <CGAL/Timer.h>
 #include <CGAL/Memory_sizer.h>
@@ -55,6 +56,7 @@ public:
 public Q_SLOTS:
   void on_actionSimplify_triggered();
 
+
 }; // end Polyhedron_demo_point_set_simplification_plugin
 
 class Point_set_demo_point_set_simplification_dialog : public QDialog, private Ui::PointSetSimplificationDialog
@@ -66,9 +68,44 @@ class Point_set_demo_point_set_simplification_dialog : public QDialog, private U
       setupUi(this);
     }
 
-    QString simplificationMethod() const { return m_simplificationMethod->currentText(); }
-    double randomSimplificationPercentage() const { return m_randomSimplificationPercentage->value(); }
-    double gridCellSize() const { return m_gridCellSize->value(); }
+  unsigned int simplificationMethod() const
+  {
+    if (Random->isChecked())
+      return 0;
+    else if (Grid->isChecked())
+      return 1;
+    else
+      return 2;
+  }
+  double randomSimplificationPercentage() const { return m_randomSimplificationPercentage->value(); }
+  double gridCellSize() const { return m_gridCellSize->value(); }
+  unsigned int maximumClusterSize() const { return m_maximumClusterSize->value(); }
+  double maximumSurfaceVariation() const { return m_maximumSurfaceVariation->value(); }
+
+public Q_SLOTS:
+  
+  void on_Random_toggled (bool toggled)
+  {
+    m_randomSimplificationPercentage->setEnabled (toggled);
+    m_gridCellSize->setEnabled (!toggled);
+    m_maximumClusterSize->setEnabled (!toggled);
+    m_maximumSurfaceVariation->setEnabled (!toggled);
+  }
+  void on_Grid_toggled (bool toggled)
+  {
+    m_randomSimplificationPercentage->setEnabled (!toggled);
+    m_gridCellSize->setEnabled (toggled);
+    m_maximumClusterSize->setEnabled (!toggled);
+    m_maximumSurfaceVariation->setEnabled (!toggled);
+  }
+  void on_Hierarchy_toggled (bool toggled)
+  {
+    m_randomSimplificationPercentage->setEnabled (!toggled);
+    m_gridCellSize->setEnabled (!toggled);
+    m_maximumClusterSize->setEnabled (toggled);
+    m_maximumSurfaceVariation->setEnabled (toggled);
+  }
+
 };
 
 void Polyhedron_demo_point_set_simplification_plugin::on_actionSimplify_triggered()
@@ -97,18 +134,19 @@ void Polyhedron_demo_point_set_simplification_plugin::on_actionSimplify_triggere
     // First point to delete
     Point_set::iterator first_point_to_remove = points->end();
 
-    if (dialog.simplificationMethod() == "Random")
+    unsigned int method = dialog.simplificationMethod ();
+    if (method == 0)
     {
-      std::cerr << "Random point cloud simplification (" << dialog.randomSimplificationPercentage() <<"%)...\n";
+      std::cerr << "Point set random simplification (" << dialog.randomSimplificationPercentage() <<"%)...\n";
 
       // Computes points to remove by random simplification
       first_point_to_remove =
         CGAL::random_simplify_point_set(points->begin(), points->end(),
                                         dialog.randomSimplificationPercentage());
     }
-    else if (dialog.simplificationMethod() == "Grid Clustering")
+    else if (method == 1)
     {
-      std::cerr << "Point cloud simplification by clustering (cell size = " << dialog.gridCellSize() <<" * average spacing)...\n";
+      std::cerr << "Point set grid simplification (cell size = " << dialog.gridCellSize() <<" * average spacing)...\n";
 
       // Computes average spacing
       double average_spacing = CGAL::compute_average_spacing<Concurrency_tag>(
@@ -120,6 +158,17 @@ void Polyhedron_demo_point_set_simplification_plugin::on_actionSimplify_triggere
         CGAL::grid_simplify_point_set(points->begin(), points->end(),
                                       dialog.gridCellSize()*average_spacing);
     }
+    else
+    {
+      std::cerr << "Point set hierarchy simplification (cluster size = " << dialog.maximumClusterSize()
+		<< ", maximum variation = " << dialog.maximumSurfaceVariation() << ")...\n";
+
+      // Computes points to remove by Grid Clustering
+      first_point_to_remove =
+        CGAL::hierarchy_simplify_point_set(points->begin(), points->end(),
+					   dialog.maximumClusterSize(),
+					   dialog.maximumSurfaceVariation());
+    }
 
     std::size_t nb_points_to_remove = std::distance(first_point_to_remove, points->end());
     std::size_t memory = CGAL::Memory_sizer().virtual_size();

Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_simplification_plugin.ui

@@ -1,104 +1,166 @@
-<ui version="4.0" >
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
  <class>PointSetSimplificationDialog</class>
- <widget class="QDialog" name="PointSetSimplificationDialog" >
-  <property name="geometry" >
+ <widget class="QDialog" name="PointSetSimplificationDialog">
+  <property name="geometry">
    <rect>
     <x>0</x>
     <y>0</y>
-    <width>403</width>
-    <height>153</height>
+    <width>450</width>
+    <height>251</height>
    </rect>
   </property>
-  <property name="windowTitle" >
+  <property name="windowTitle">
    <string>Simplification</string>
   </property>
-  <layout class="QGridLayout" >
-   <item row="0" column="0" >
-    <widget class="QLabel" name="label" >
-     <property name="text" >
-      <string>Method:</string>
+  <layout class="QGridLayout">
+   <item row="7" column="0">
+    <widget class="QLabel" name="label_4">
+     <property name="text">
+      <string>Maximum cluster size</string>
      </property>
     </widget>
    </item>
-   <item row="0" column="1" >
-    <widget class="QComboBox" name="m_simplificationMethod" >
-     <item>
-      <property name="text" >
-       <string>Random</string>
-      </property>
-     </item>
-     <item>
-      <property name="text" >
-       <string>Grid Clustering</string>
-      </property>
-     </item>
-    </widget>
-   </item>
-   <item row="1" column="0" >
-    <widget class="QLabel" name="label_2" >
-     <property name="text" >
-      <string>Points to Remove Randomly</string>
+   <item row="9" column="0" colspan="2">
+    <widget class="QDialogButtonBox" name="buttonBox">
+     <property name="orientation">
+      <enum>Qt::Horizontal</enum>
+     </property>
+     <property name="standardButtons">
+      <set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
      </property>
     </widget>
    </item>
-   <item row="1" column="1" >
-    <widget class="QDoubleSpinBox" name="m_randomSimplificationPercentage" >
-     <property name="suffix" >
+   <item row="2" column="0">
+    <widget class="QRadioButton" name="Random">
+     <property name="text">
+      <string>Random</string>
+     </property>
+     <property name="checked">
+      <bool>true</bool>
+     </property>
+    </widget>
+   </item>
+   <item row="4" column="0">
+    <widget class="QRadioButton" name="Grid">
+     <property name="text">
+      <string>Grid</string>
+     </property>
+    </widget>
+   </item>
+   <item row="3" column="1">
+    <widget class="QDoubleSpinBox" name="m_randomSimplificationPercentage">
+     <property name="suffix">
       <string> %</string>
      </property>
-     <property name="decimals" >
+     <property name="decimals">
       <number>2</number>
      </property>
-     <property name="minimum" >
+     <property name="minimum">
       <double>0.100000000000000</double>
      </property>
-     <property name="maximum" >
+     <property name="maximum">
       <double>100.000000000000000</double>
      </property>
-     <property name="singleStep" >
+     <property name="singleStep">
       <double>0.100000000000000</double>
      </property>
-     <property name="value" >
+     <property name="value">
       <double>50.000000000000000</double>
      </property>
     </widget>
    </item>
-   <item row="2" column="0" >
-    <widget class="QLabel" name="label_3" >
-     <property name="text" >
-      <string>Grid Cell Size</string>
+   <item row="6" column="0">
+    <widget class="QRadioButton" name="Hierarchy">
+     <property name="text">
+      <string>Hierarchy</string>
      </property>
     </widget>
    </item>
-   <item row="2" column="1" >
-    <widget class="QDoubleSpinBox" name="m_gridCellSize" >
-     <property name="suffix" >
+   <item row="5" column="1">
+    <widget class="QDoubleSpinBox" name="m_gridCellSize">
+     <property name="enabled">
+      <bool>false</bool>
+     </property>
+     <property name="suffix">
       <string> * average spacing</string>
      </property>
-     <property name="decimals" >
+     <property name="decimals">
       <number>2</number>
      </property>
-     <property name="minimum" >
+     <property name="minimum">
       <double>0.100000000000000</double>
      </property>
-     <property name="maximum" >
+     <property name="maximum">
       <double>10.000000000000000</double>
      </property>
-     <property name="singleStep" >
+     <property name="singleStep">
       <double>0.100000000000000</double>
      </property>
-     <property name="value" >
+     <property name="value">
       <double>1.000000000000000</double>
      </property>
     </widget>
    </item>
-   <item row="3" column="0" colspan="2" >
-    <widget class="QDialogButtonBox" name="buttonBox" >
-     <property name="orientation" >
-      <enum>Qt::Horizontal</enum>
+   <item row="3" column="0">
+    <widget class="QLabel" name="label_2">
+     <property name="text">
+      <string>Points to Remove Randomly</string>
      </property>
-     <property name="standardButtons" >
-      <set>QDialogButtonBox::Cancel|QDialogButtonBox::NoButton|QDialogButtonBox::Ok</set>
+    </widget>
+   </item>
+   <item row="5" column="0">
+    <widget class="QLabel" name="label_3">
+     <property name="text">
+      <string>Grid Cell Size</string>
+     </property>
+    </widget>
+   </item>
+   <item row="7" column="1">
+    <widget class="QSpinBox" name="m_maximumClusterSize">
+     <property name="enabled">
+      <bool>false</bool>
+     </property>
+     <property name="minimum">
+      <number>1</number>
+     </property>
+     <property name="maximum">
+      <number>2147483647</number>
+     </property>
+     <property name="value">
+      <number>10</number>
+     </property>
+    </widget>
+   </item>
+   <item row="8" column="0">
+    <widget class="QLabel" name="label_5">
+     <property name="text">
+      <string>Maximum surface variation</string>
+     </property>
+    </widget>
+   </item>
+   <item row="8" column="1">
+    <widget class="QDoubleSpinBox" name="m_maximumSurfaceVariation">
+     <property name="enabled">
+      <bool>false</bool>
+     </property>
+     <property name="suffix">
+      <string/>
+     </property>
+     <property name="decimals">
+      <number>5</number>
+     </property>
+     <property name="minimum">
+      <double>0.000010000000000</double>
+     </property>
+     <property name="maximum">
+      <double>0.333330000000000</double>
+     </property>
+     <property name="singleStep">
+      <double>0.012340000000000</double>
+     </property>
+     <property name="value">
+      <double>0.333330000000000</double>
      </property>
     </widget>
    </item>
@@ -112,11 +174,11 @@
    <receiver>PointSetSimplificationDialog</receiver>
    <slot>accept()</slot>
    <hints>
-    <hint type="sourcelabel" >
+    <hint type="sourcelabel">
      <x>248</x>
      <y>254</y>
     </hint>
-    <hint type="destinationlabel" >
+    <hint type="destinationlabel">
      <x>157</x>
      <y>274</y>
     </hint>
@@ -128,11 +190,11 @@
    <receiver>PointSetSimplificationDialog</receiver>
    <slot>reject()</slot>
    <hints>
-    <hint type="sourcelabel" >
+    <hint type="sourcelabel">
      <x>316</x>
      <y>260</y>
     </hint>
-    <hint type="destinationlabel" >
+    <hint type="destinationlabel">
      <x>286</x>
      <y>274</y>
     </hint>