Merge branch 'Point_set_shape_detection_3-Plane_regularization-GF-old' into Point_set_shape_detection_3-Plane_regularization-GF

Documentation/doc/biblio/cgal_manual.bib

@@ -1941,6 +1941,15 @@ location    = {Salt Lake City, Utah, USA}
   ,update       = "98.01 kettner"
 }
 
+
+@techreport{cgal:vla-lod-15,
+  title={LOD Generation for Urban Scenes},
+  author={Verdie, Yannick and Lafarge, Florent and Alliez, Pierre},
+  year={2015},
+  institution={Association for Computing Machinery}
+}
+
+
 @book{ cgal:vj-ctcg-03
   ,author       = "David Vandevoorde and Nicolai M. Josuttis"
   ,title        = "{C}++ Templates: The Complete Guide"

Installation/changes.html

@@ -269,6 +269,13 @@ and <code>src/</code> directories).
       <code>CGAL::write_ply_points()</code>
       and <code>CGAL::write_ply_points_and_normals()</code>.</li>
   </ul>
+  <h3>Point Set Shape Detection</h3>
+  <ul>
+    <li>New post-processing
+      algorithm: <code>CGAL::regularize_planes()</code>. This allows the user
+      to favor parallelism, orthogonality, coplanarity and/or axial
+      symmetry between detected planes.</li>
+  </ul>
   <h3>Surface Mesh Parameterization</h3>
   <ul>
     <li><code>LSCM_parameterizer_3</code> now uses by default Eigen

Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt

@@ -47,5 +47,8 @@
 - `CGAL::Shape_detection_3::Cylinder<Traits>`
 - `CGAL::Shape_detection_3::Cone<Traits>`
 - `CGAL::Shape_detection_3::Torus<Traits>`
+
+## Functions ##
+- `CGAL::regularize_planes()`
 */
 

Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt

@@ -105,6 +105,24 @@ which is used by the example  \ref Point_set_shape_detection_3/efficient_RANSAC_
 
 \cgalExample{Point_set_shape_detection_3/efficient_RANSAC_custom_shape.h}
 
+\section Point_set_shape_detection_3Plane_regularization Plane Regularization
+
+Shape detection is very suited for man-made shapes such as urban scenes or scans of mechanical pieces. Such scenes may contain a wide variety of regularities that user may want to favor: parallelism, orthogonality, concentricity, coaxiality, etc., not to mention orbits and composed relationships.
+
+Among all these possibilities, \cgal provides a means to regularize four properties of regularity through a function `CGAL::regularize_planes()`. It only postprocesses planes detected by `CGAL::Shape_detection_3<Traits>` (other primitives are left unchanged):
+
+- Planes that are near __parallel__ are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal.
+
+- Parallel planes that are near __coplanar__ are made coplanar.
+
+- Planes that are near __orthogonal__ are made exactly orthogonal.
+
+- Planes that are near __symmetrical__ with respect to a user-defined axis are made symmetrical.
+
+The user can choose to only regularize one or several of these 4 properties (see reference manual). The process is greedy and based on a hierarchical decomposition (coplanar clusters are subgroups of parallel clusters which are subgroups of axis-symmetric and orthogonal clusters) as described by Verdie et al. \cgalCite{cgal:vla-lod-15}
+
+\cgalExample{Point_set_shape_detection_3/plane_regularization.cpp}
+
 \section Point_set_shape_detection_3Performance Performance
  
 The running time and detection performance depend on the chosen parameters. A selective error tolerance parameter leads to higher running times and smaller shapes, as many shape candidates are generated to find the largest shape. We plot the detection performance against the epsilon error tolerance parameter for detecting planes in a complex scene with 5M points, see \cgalFigureRef{Point_set_shape_detection_3_performace_epsilon}. The probability parameter controls the endurance when searching for the largest candidate at each iteration. It barely impacts the number of detected shapes, has a moderate impact on the size of the detected shapes and increases the running times. We plot the performance against the probability parameter, see \cgalFigureRef{Point_set_shape_detection_3_performace_probability}.

Point_set_shape_detection_3/doc/Point_set_shape_detection_3/examples.txt

@@ -4,4 +4,5 @@
 \example Point_set_shape_detection_3/efficient_RANSAC_point_access.cpp
 \example Point_set_shape_detection_3/efficient_RANSAC_custom_shape.cpp
 \example Point_set_shape_detection_3/efficient_RANSAC_custom_shape.h
+\example Point_set_shape_detection_3/plane_regularization.cpp
 */

Point_set_shape_detection_3/examples/Point_set_shape_detection_3/CMakeLists.txt

@@ -27,6 +27,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "efficient_RANSAC_custom_shape.cpp" )
   create_single_source_cgal_program( "efficient_RANSAC_parameters.cpp" )
   create_single_source_cgal_program( "efficient_RANSAC_point_access.cpp" )
+  create_single_source_cgal_program( "plane_regularization.cpp" )
 
 else()
   

Point_set_shape_detection_3/examples/Point_set_shape_detection_3/plane_regularization.cpp

@@ -0,0 +1,53 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/Point_with_normal_3.h>
+#include <CGAL/property_map.h>
+
+#include <CGAL/Shape_detection_3.h>
+#include <CGAL/regularize_planes.h>
+
+#include <iostream>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
+typedef std::pair<Kernel::Point_3, Kernel::Vector_3>         Point_with_normal;
+typedef std::vector<Point_with_normal>                       Pwn_vector;
+typedef CGAL::First_of_pair_property_map<Point_with_normal>  Point_map;
+typedef CGAL::Second_of_pair_property_map<Point_with_normal> Normal_map;
+
+typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits
+  <Kernel, Pwn_vector, Point_map, Normal_map>                Traits;
+typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits>    Efficient_ransac;
+typedef CGAL::Shape_detection_3::Plane<Traits>               Plane;
+
+int main() 
+{
+  Pwn_vector points;
+  std::ifstream stream("data/cube.pwn");
+
+  if (!stream || 
+    !CGAL::read_xyz_points_and_normals(stream,
+      std::back_inserter(points),
+      Point_map(),
+      Normal_map()))
+  {
+      std::cerr << "Error: cannot read file cube.pwn" << std::endl;
+      return EXIT_FAILURE;
+  }
+
+  // Call RANSAC shape detection with planes
+  Efficient_ransac ransac;
+  ransac.set_input(points);
+  ransac.add_shape_factory<Plane>();
+  ransac.detect();
+
+  // Regularize detected planes
+  CGAL::regularize_planes (ransac,
+                           true, // Regularize parallelism
+                           true, // Regularize orthogonality
+                           false, // Do not regularize coplanarity
+                           true, // Regularize Z-symmetry (default)
+                           10); // 10 degrees of tolerance for parallelism/orthogonality
+  
+  return EXIT_SUCCESS;
+}

Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Efficient_RANSAC.h

@@ -209,6 +209,15 @@ shape. The implementation follows \cgalCite{schnabel2007efficient}.
       return m_traits;
     }
 
+    Input_iterator input_iterator_first() const
+    {
+      return m_input_iterator_first;
+    }
+    Input_iterator input_iterator_beyond() const
+    {
+      return m_input_iterator_beyond;
+    }
+
     /*!
       Sets the input data. The range must stay valid
       until the detection has been performed and the access to the

Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h

@@ -46,6 +46,7 @@ namespace CGAL {
      ///< property map to access the unoriented normal of an input point.
     typedef typename Traits::FT FT; ///< number type.
     typedef typename Traits::Point_3 Point_3; ///< point type.
+    typedef typename Traits::Point_2 Point_2; ///< point 2D type.
     typedef typename Traits::Vector_3 Vector_3;
     /// \endcond
 
@@ -85,7 +86,36 @@ namespace CGAL {
       return d * d;
     }
 
-    
+    /*!
+      Computes the orthogonal projection of a query point on the shape.
+     */
+    Point_3 projection (const Point_3& p) const {
+      return to_3d (to_2d (p));
+    }
+
+    Point_2 to_2d (const Point_3& p) const {
+      Vector_3 v (m_point_on_primitive, p);
+      return Point_2 (v * m_base1, v * m_base2);
+    }
+
+    Point_3 to_3d (const Point_2& p) const {
+      return m_point_on_primitive + p.x () * m_base1 + p.y () * m_base2;
+    }
+
+    /*!
+      Replaces the plane by p
+    */
+    void update (const Plane_3& p) {
+      m_base1 = p.base1 () / std::sqrt (p.base1() * p.base1 ());
+      m_base2 = p.base2 () / std::sqrt (p.base2() * p.base2 ());
+      m_normal = p.orthogonal_vector () / std::sqrt (p.orthogonal_vector () * p.orthogonal_vector ());
+      
+      m_d = -(this->get_x(m_point_on_primitive) * this->get_x(m_normal) 
+        + this->get_y(m_point_on_primitive) * this->get_y(m_normal) 
+        + this->get_z(m_point_on_primitive) * this->get_z(m_normal));
+
+    }
+
     /*!
       Helper function to write the plane equation and
       number of assigned points into a string.

Point_set_shape_detection_3/include/CGAL/regularize_planes.h

@@ -0,0 +1,768 @@
+// Copyright (c) 2015 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)     : Florent Lafarge, Simon Giraudot
+//
+
+/**
+ * \ingroup PkgPointSetShapeDetection3
+ * \file CGAL/regularize_planes.h
+ *
+ */
+
+
+#ifndef CGAL_REGULARIZE_PLANES_H
+#define CGAL_REGULARIZE_PLANES_H
+
+#include <CGAL/Shape_detection_3.h>
+#include <CGAL/centroid.h>
+
+#include <boost/foreach.hpp>
+
+
+namespace CGAL {
+  
+// ----------------------------------------------------------------------------
+// Private section
+// ----------------------------------------------------------------------------
+/// \cond SKIP_IN_MANUAL
+namespace internal {
+namespace PlaneRegularization {
+
+template <typename Traits>
+struct Plane_cluster
+{
+  bool is_free;
+  std::vector<std::size_t> planes;
+  std::vector<std::size_t> coplanar_group;
+  std::vector<std::size_t> orthogonal_clusters;
+  typename Traits::Vector_3 normal;
+  typename Traits::FT cosangle_symmetry;
+  typename Traits::FT area;
+  typename Traits::FT cosangle_centroid;
+};
+
+  
+template <typename Traits>
+typename Traits::Vector_3 regularize_normal
+  (const typename Traits::Vector_3& n,
+   const typename Traits::Vector_3& symmetry_direction,
+   typename Traits::FT cos_symmetry)
+{
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Point_3 Point;
+  typedef typename Traits::Vector_3 Vector;
+  typedef typename Traits::Line_3 Line;
+  typedef typename Traits::Plane_3 Plane;
+
+  if (symmetry_direction == CGAL::NULL_VECTOR)
+    return n;
+      
+  Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction;
+
+  Plane plane_symmetry (pt_symmetry, symmetry_direction);
+  Point pt_normal = CGAL::ORIGIN + n;
+
+  if (n != symmetry_direction || n != -symmetry_direction)
+    {
+      Plane plane_cut (CGAL::ORIGIN, pt_normal, CGAL::ORIGIN + symmetry_direction);
+      Line line;
+      CGAL::Object ob_1 = CGAL::intersection(plane_cut, plane_symmetry);
+      if (!assign(line, ob_1))
+        return n;
+
+      FT delta = std::sqrt ((FT)1. - cos_symmetry * cos_symmetry);
+
+      Point projected_origin = line.projection (CGAL::ORIGIN);
+      Vector line_vector (line);
+      line_vector = line_vector / std::sqrt (line_vector * line_vector);
+      Point pt1 = projected_origin + delta * line_vector;
+      Point pt2 = projected_origin - delta * line_vector;
+
+      if (CGAL::squared_distance (pt_normal, pt1) <= CGAL::squared_distance (pt_normal, pt2))
+        return Vector (CGAL::ORIGIN, pt1);
+      else
+        return Vector (CGAL::ORIGIN, pt2);
+
+    }
+  else
+    return n;
+}
+
+template <typename Traits>  
+typename Traits::Vector_3 regularize_normals_from_prior
+  (const typename Traits::Vector_3& np,
+   const typename Traits::Vector_3& n,
+   const typename Traits::Vector_3& symmetry_direction,
+   typename Traits::FT cos_symmetry)
+{
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Point_3 Point;
+  typedef typename Traits::Vector_3 Vector;
+  typedef typename Traits::Line_3 Line;
+  typedef typename Traits::Plane_3 Plane;
+
+  if (symmetry_direction == CGAL::NULL_VECTOR)
+    return n;
+
+  Plane plane_orthogonality (CGAL::ORIGIN, np);
+  Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction;
+
+  Plane plane_symmetry (pt_symmetry, symmetry_direction);
+		
+  Line line;
+  CGAL::Object ob_1 = CGAL::intersection (plane_orthogonality, plane_symmetry);
+  if (!assign(line, ob_1))
+    return regularize_normal<Traits> (n, symmetry_direction, cos_symmetry);
+
+  Point projected_origin = line.projection (CGAL::ORIGIN);
+  FT R = CGAL::squared_distance (Point (CGAL::ORIGIN), projected_origin);
+
+  if (R <= 1)  // 2 (or 1) possible points intersecting the unit sphere and line
+    {
+      FT delta = std::sqrt ((FT)1. - R);
+      Vector line_vector(line); 
+      line_vector = line_vector / std::sqrt (line_vector * line_vector);
+      Point pt1 = projected_origin + delta * line_vector;
+      Point pt2 = projected_origin - delta * line_vector;
+			
+      Point pt_n = CGAL::ORIGIN + n;
+      if (CGAL::squared_distance (pt_n, pt1) <= CGAL::squared_distance (pt_n, pt2))
+        return Vector (CGAL::ORIGIN, pt1);
+      else
+        return Vector (CGAL::ORIGIN, pt2);
+    }
+  else //no point intersecting the unit sphere and line
+    return regularize_normal<Traits> (n,symmetry_direction, cos_symmetry);
+
+}
+
+template <typename Traits,
+          typename RandomAccessIterator,
+          typename PlaneContainer,
+          typename PointPMap,
+          typename CentroidContainer,
+          typename AreaContainer>
+void compute_centroids_and_areas (RandomAccessIterator input_begin,
+                                  PlaneContainer& planes,
+                                  PointPMap point_pmap,
+                                  CentroidContainer& centroids,
+                                  AreaContainer& areas)
+{
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Point_3 Point;
+  
+  for (std::size_t i = 0; i < planes.size (); ++ i)
+    {
+      std::vector < Point > listp;
+      for (std::size_t j = 0; j < planes[i]->indices_of_assigned_points ().size (); ++ j)
+        {
+          std::size_t yy = planes[i]->indices_of_assigned_points()[j];
+          Point pt = get (point_pmap, *(input_begin + yy));
+          listp.push_back(pt);
+        }
+      centroids.push_back (CGAL::centroid (listp.begin (), listp.end ()));
+      areas.push_back ((FT)(planes[i]->indices_of_assigned_points().size()) / (FT)100.);
+    }
+}
+
+
+template <typename Traits,
+          typename PlaneContainer,
+          typename PlaneClusterContainer,
+          typename AreaContainer>
+void compute_parallel_clusters (PlaneContainer& planes,
+                                PlaneClusterContainer& clusters,
+                                AreaContainer& areas,
+                                typename Traits::FT tolerance_cosangle,
+                                const typename Traits::Vector_3& symmetry_direction)
+{
+
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Vector_3 Vector;
+  
+  typedef typename PlaneClusterContainer::value_type Plane_cluster;
+  
+  // find pairs of epsilon-parallel primitives and store them in parallel_planes
+  std::vector<std::vector<std::size_t> > parallel_planes (planes.size ());
+  for (std::size_t i = 0; i < planes.size (); ++ i)
+    {
+      Vector v1 = planes[i]->plane_normal ();
+          
+      for (std::size_t j = 0; j < planes.size(); ++ j)
+        {
+          if (i == j)
+            continue;
+              
+          Vector v2 = planes[i]->plane_normal ();
+              
+          if (std::fabs (v1 * v2) > 1. - tolerance_cosangle)
+            parallel_planes[i].push_back (j);
+        }
+    }
+
+
+  std::vector<bool> is_available (planes.size (), true);
+      
+  for (std::size_t i = 0; i < planes.size(); ++ i)
+    {
+
+      if(is_available[i])
+        {
+          is_available[i] = false;
+
+          clusters.push_back (Plane_cluster());
+          Plane_cluster& clu = clusters.back ();
+
+          //initialization containers
+          clu.planes.push_back (i);
+              
+          std::vector<std::size_t> index_container_former_ring_parallel;
+          index_container_former_ring_parallel.push_back(i);
+              
+          std::list<std::size_t> index_container_current_ring_parallel;
+
+          //propagation over the pairs of epsilon-parallel primitives
+          bool propagation=true;
+          clu.normal = planes[i]->plane_normal ();
+          clu.area = areas[i];
+			
+          do
+            {
+              propagation = false;
+
+              for (std::size_t k = 0; k < index_container_former_ring_parallel.size(); ++ k)
+                {
+
+                  std::size_t plane_index = index_container_former_ring_parallel[k];
+
+                  for (std::size_t l = 0; l < parallel_planes[plane_index].size(); ++ l)
+                    {
+                      std::size_t it = parallel_planes[plane_index][l];
+                          
+                      Vector normal_it =  planes[it]->plane_normal ();
+
+                      if(is_available[it]
+                         && std::fabs (normal_it*clu.normal) > 1. - tolerance_cosangle )
+                        {	
+                          propagation = true;
+                          index_container_current_ring_parallel.push_back(it);
+                          is_available[it]=false;
+                              
+                          if(clu.normal * normal_it <0)
+                            normal_it = -normal_it;
+
+                          clu.normal = (FT)clu.area * clu.normal
+                            + (FT)areas[it] * normal_it;
+                          FT norm = (FT)1. / std::sqrt (clu.normal.squared_length()); 
+                          clu.normal = norm * clu.normal;
+                          clu.area += areas[it];
+                        }	
+                    }
+                }
+
+              //update containers
+              index_container_former_ring_parallel.clear();
+              for (std::list<std::size_t>::iterator it = index_container_current_ring_parallel.begin();
+                   it != index_container_current_ring_parallel.end(); ++it)
+                {
+                  index_container_former_ring_parallel.push_back(*it);
+                  clu.planes.push_back(*it);
+                }
+              index_container_current_ring_parallel.clear();
+
+            }
+          while(propagation);
+
+          if (symmetry_direction != CGAL::NULL_VECTOR)
+            clu.cosangle_symmetry = std::fabs(symmetry_direction * clu.normal);
+        }
+    }
+  is_available.clear();
+}
+
+template <typename Traits,
+          typename PlaneClusterContainer>
+void cluster_symmetric_cosangles (PlaneClusterContainer& clusters,
+                                  typename Traits::FT tolerance_cosangle)
+{
+  typedef typename Traits::FT FT;
+  
+  std::vector < FT > cosangle_centroids;
+  std::vector < std::size_t> list_cluster_index;
+  for( std::size_t i = 0; i < clusters.size(); ++ i)
+    list_cluster_index.push_back(static_cast<std::size_t>(-1));
+      
+  std::size_t mean_index = 0;
+  for (std::size_t i = 0; i < clusters.size(); ++ i)
+    {
+      if(list_cluster_index[i] == static_cast<std::size_t>(-1))
+        {
+          list_cluster_index[i] = mean_index;
+          FT mean = clusters[i].area * clusters[i].cosangle_symmetry;
+          FT mean_area = clusters[i].area;
+              
+          for (std::size_t j = i+1; j < clusters.size(); ++ j)
+            {
+              if (list_cluster_index[j] == static_cast<std::size_t>(-1)
+                  && std::fabs (clusters[j].cosangle_symmetry -
+                                mean / mean_area) < tolerance_cosangle)
+                {
+                  list_cluster_index[j] = mean_index;
+                  mean_area += clusters[j].area;
+                  mean += clusters[j].area * clusters[j].cosangle_symmetry;
+                }
+            }
+          ++ mean_index;
+          mean /= mean_area;
+          cosangle_centroids.push_back (mean);
+        }
+    }
+
+  for (std::size_t i = 0; i < cosangle_centroids.size(); ++ i)
+    {
+      if (cosangle_centroids[i] < tolerance_cosangle)
+        cosangle_centroids[i] = 0;
+      else if (cosangle_centroids[i] > 1. - tolerance_cosangle)
+        cosangle_centroids[i] = 1;
+    }
+  for (std::size_t i = 0; i < clusters.size(); ++ i)
+    clusters[i].cosangle_symmetry = cosangle_centroids[list_cluster_index[i]];
+}
+
+
+template <typename Traits,
+          typename PlaneClusterContainer>
+void subgraph_mutually_orthogonal_clusters (PlaneClusterContainer& clusters,
+                                            const typename Traits::Vector_3& symmetry_direction)
+{
+  typedef typename Traits::FT FT;
+  typedef typename Traits::Vector_3 Vector;
+  
+  std::vector < std::vector < std::size_t> > subgraph_clusters;
+  std::vector < std::size_t> subgraph_clusters_max_area_index;
+
+  for (std::size_t i = 0; i < clusters.size(); ++ i)
+    clusters[i].is_free = true;
+
+  for (std::size_t i = 0; i < clusters.size(); ++ i)
+    {
+      if(clusters[i].is_free)
+        {
+          clusters[i].is_free = false;
+          FT max_area = clusters[i].area;
+          std::size_t index_max_area = i;
+
+          //initialization containers
+          std::vector < std::size_t > index_container;
+          index_container.push_back(i);
+          std::vector < std::size_t > index_container_former_ring;
+          index_container_former_ring.push_back(i);
+          std::list < std::size_t > index_container_current_ring;
+
+          //propagation
+          bool propagation=true;
+          do
+            {
+              propagation=false;
+
+              //neighbors
+              for (std::size_t k=0;k<index_container_former_ring.size();k++)
+                {
+
+                  std::size_t cluster_index=index_container_former_ring[k];
+
+                  for (std::size_t j = 0; j < clusters[cluster_index].orthogonal_clusters.size(); ++ j)
+                    {
+                      if(clusters[j].is_free)
+                        { 	
+                          propagation = true;
+                          index_container_current_ring.push_back(j);
+                          clusters[j].is_free = false;
+
+                          if(max_area < clusters[j].area)
+                            {
+                              max_area = clusters[j].area;
+                              index_max_area = j;
+                            }
+                        }	
+                    }
+                }
+
+              //update containers
+              index_container_former_ring.clear();
+              for(std::list < std::size_t>::iterator it = index_container_current_ring.begin();
+                  it != index_container_current_ring.end(); ++it)
+                {
+                  index_container_former_ring.push_back(*it);
+                  index_container.push_back(*it);
+                }
+              index_container_current_ring.clear();
+
+            }
+          while(propagation);
+          subgraph_clusters.push_back(index_container);
+          subgraph_clusters_max_area_index.push_back(index_max_area);
+        }
+    }
+
+
+  //create subgraphs of mutually orthogonal clusters in which the
+  //largest cluster is excluded and store in
+  //subgraph_clusters_prop
+  std::vector < std::vector < std::size_t> > subgraph_clusters_prop;
+  for (std::size_t i=0;i<subgraph_clusters.size(); i++)
+    {
+      std::size_t index=subgraph_clusters_max_area_index[i];
+      std::vector < std::size_t> subgraph_clusters_prop_temp;
+      for (std::size_t j=0;j<subgraph_clusters[i].size(); j++)
+        if(subgraph_clusters[i][j]!=index)
+          subgraph_clusters_prop_temp.push_back(subgraph_clusters[i][j]);
+
+      subgraph_clusters_prop.push_back(subgraph_clusters_prop_temp);
+    }
+
+
+
+  //regularization of cluster normals : in eachsubgraph, we start
+  //from the largest area cluster and we propage over the subgraph
+  //by regularizing the normals of the clusters accorting to
+  //orthogonality and cosangle to symmetry direction
+
+  for (std::size_t i = 0; i < clusters.size(); ++ i)
+    clusters[i].is_free = true;
+
+  for (std::size_t i = 0; i < subgraph_clusters_prop.size(); ++ i)
+    {
+	
+      std::size_t index_current=subgraph_clusters_max_area_index[i];
+      Vector vec_current=regularize_normal<Traits>
+        (clusters[index_current].normal,
+         symmetry_direction,
+         clusters[index_current].cosangle_symmetry);
+      clusters[index_current].normal = vec_current;
+      clusters[index_current].is_free = false;
+
+      //initialization containers
+      std::vector < std::size_t> index_container;
+      index_container.push_back(index_current);
+      std::vector < std::size_t> index_container_former_ring;
+      index_container_former_ring.push_back(index_current);
+      std::list < std::size_t> index_container_current_ring;
+
+      //propagation
+      bool propagation=true;
+      do
+        {
+          propagation=false;
+
+          //neighbors
+          for (std::size_t k=0;k<index_container_former_ring.size();k++)
+            {
+
+              std::size_t cluster_index=index_container_former_ring[k];
+
+              for (std::size_t j = 0; j < clusters[cluster_index].orthogonal_clusters.size(); ++ j)
+                {
+						
+                  if(clusters[j].is_free)
+                    { 	
+							
+                      propagation = true;
+                      index_container_current_ring.push_back(j);
+                      clusters[j].is_free = false;
+
+                      Vector new_vect=regularize_normals_from_prior<Traits>
+                        (clusters[cluster_index].normal,
+                         clusters[j].normal,
+                         symmetry_direction,
+                         clusters[j].cosangle_symmetry);
+                      clusters[j].normal = new_vect;
+                    }
+                }	
+            }
+			
+          //update containers
+          index_container_former_ring.clear();
+          for(std::list < std::size_t>::iterator it = index_container_current_ring.begin();
+              it != index_container_current_ring.end(); ++it)
+            {
+              index_container_former_ring.push_back(*it);
+              index_container.push_back(*it);
+            }
+          index_container_current_ring.clear();
+        }while(propagation);
+    }
+}
+                                    
+
+
+} // namespace PlaneRegularization
+} // namespace internal
+/// \endcond
+
+
+// ----------------------------------------------------------------------------
+// Public section
+// ----------------------------------------------------------------------------
+
+/// \ingroup PkgPointSetShapeDetection3
+  
+  /*! 
+
+    Given a set of detected planes with their respective inlier sets,
+    this function enables to regularize the planes: 
+
+    - Planes near parallel can be made exactly parallel.
+
+    - Planes near orthogonal can be made exactly orthogonal.
+
+    - Planes parallel and near coplanar can be made exactly coplanar.
+
+    - Planes near symmetrical with a user-defined axis can be made
+    exactly symmetrical.
+
+    Planes are directly modified. Points are left unaltered, as well as
+    their relationships to planes (no transfer of point from a primitive
+    plane to another).
+
+    The implementation follows \cgalCite{cgal:vla-lod-15}.
+
+    \tparam Traits a model of `EfficientRANSACTraits`
+
+    \param shape_detection Shape detection object used to detect
+    shapes from the input data. This engine may handle any types of
+    primitive shapes but only planes will be regularized. 
+
+    \warning The `shape_detection` parameter must have already
+    detected shapes. If no plane exists in it, the regularization
+    function doesn't do anything.
+
+    \param regularize_parallelism Select whether parallelism is
+    regularized or not.
+
+    \param regularize_orthogonality Select whether orthogonality is
+    regularized or not.
+
+    \param regularize_coplanarity Select whether coplanarity is
+    regularized or not.
+
+    \param regularize_axis_symmetry Select whether axis symmetry is
+    regularized or not.
+
+    \param tolerance_angle Tolerance of deviation between normal
+    vectors of planes (in degrees) used for parallelism, orthogonality
+    and axis symmetry. Default value is 25 degrees.
+
+    \param tolerance_coplanarity Maximal distance between two parallel
+    planes such that they are considered coplanar. Default value is
+    0.01.
+
+    \param symmetry_direction Chosen axis for symmetry
+    regularization. Default value is the Z axis.
+*/ 
+
+template <typename EfficientRANSACTraits>
+void regularize_planes (const Shape_detection_3::Efficient_RANSAC<EfficientRANSACTraits>& shape_detection,
+                        bool regularize_parallelism,
+                        bool regularize_orthogonality,
+                        bool regularize_coplanarity,
+                        bool regularize_axis_symmetry,
+                        typename EfficientRANSACTraits::FT tolerance_angle
+                        = (typename EfficientRANSACTraits::FT)25.0,
+                        typename EfficientRANSACTraits::FT tolerance_coplanarity
+                        = (typename EfficientRANSACTraits::FT)0.01,
+                        typename EfficientRANSACTraits::Vector_3 symmetry_direction
+                        = typename EfficientRANSACTraits::Vector_3 (0., 0., 1.))
+{
+  typedef typename EfficientRANSACTraits::FT FT;
+  typedef typename EfficientRANSACTraits::Point_3 Point;
+  typedef typename EfficientRANSACTraits::Vector_3 Vector;
+  typedef typename EfficientRANSACTraits::Plane_3 Plane;
+  typedef typename EfficientRANSACTraits::Point_map Point_map;
+
+  typedef Shape_detection_3::Shape_base<EfficientRANSACTraits> Shape;
+  typedef Shape_detection_3::Plane<EfficientRANSACTraits> Plane_shape;
+
+  typedef typename internal::PlaneRegularization::Plane_cluster<EfficientRANSACTraits>
+    Plane_cluster;
+
+  typename EfficientRANSACTraits::Input_range::iterator input_begin = shape_detection.input_iterator_first();
+
+  std::vector<boost::shared_ptr<Plane_shape> > planes;
+    
+  BOOST_FOREACH (boost::shared_ptr<Shape> shape, shape_detection.shapes())
+    {
+      boost::shared_ptr<Plane_shape> pshape
+        = boost::dynamic_pointer_cast<Plane_shape>(shape);
+        
+      // Ignore all shapes other than plane
+      if (pshape == boost::shared_ptr<Plane_shape>())
+        continue;
+      planes.push_back (pshape);
+    }
+
+
+  /*
+   * Compute centroids and areas
+   */
+  std::vector<Point> centroids;
+  std::vector<FT> areas;
+  internal::PlaneRegularization::compute_centroids_and_areas<EfficientRANSACTraits>
+    (input_begin, planes, Point_map(), centroids, areas);
+
+  FT tolerance_cosangle = (FT)1. - std::cos (tolerance_angle);
+      
+  // clustering the parallel primitives and store them in clusters
+  // & compute the normal, size and cos angle to the symmetry
+  // direction of each cluster
+  std::vector<Plane_cluster> clusters;
+  internal::PlaneRegularization::compute_parallel_clusters<EfficientRANSACTraits>
+    (planes, clusters, areas,
+     (regularize_parallelism ? tolerance_cosangle : (FT)0.0),
+     (regularize_axis_symmetry ? symmetry_direction : CGAL::NULL_VECTOR));
+
+  if (regularize_orthogonality)
+    {
+      //discovery orthogonal relationship between clusters 
+      for (std::size_t i = 0; i < clusters.size(); ++ i)
+        {
+          for (std::size_t j = i + 1; j < clusters.size(); ++ j)
+            {
+              
+              if (std::fabs (clusters[i].normal * clusters[j].normal) < tolerance_cosangle)
+                {
+                  clusters[i].orthogonal_clusters.push_back (j);
+                  clusters[j].orthogonal_clusters.push_back (i);
+                }
+            }
+        }
+    }
+      
+  if (regularize_axis_symmetry)
+    {
+      //clustering the symmetry cosangle and store their centroids in
+      //cosangle_centroids and the centroid index of each cluster in
+      //list_cluster_index
+      internal::PlaneRegularization::cluster_symmetric_cosangles<EfficientRANSACTraits>
+        (clusters, tolerance_cosangle);
+    }
+  
+  //find subgraphs of mutually orthogonal clusters (store index of
+  //clusters in subgraph_clusters), and select the cluster of
+  //largest area
+  internal::PlaneRegularization::subgraph_mutually_orthogonal_clusters<EfficientRANSACTraits>
+    (clusters, symmetry_direction);
+      
+  //recompute optimal plane for each primitive after normal regularization
+  for (std::size_t i=0; i < clusters.size(); ++ i)
+    {
+
+      Vector vec_reg = clusters[i].normal;
+
+      for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j)
+        {
+          std::size_t index_prim = clusters[i].planes[j];
+          Point pt_reg = planes[index_prim]->projection (centroids[index_prim]);
+          if( planes[index_prim]->plane_normal () * vec_reg < 0)
+            vec_reg=-vec_reg;
+          Plane plane_reg(pt_reg,vec_reg);
+              
+          if( std::fabs(planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - tolerance_cosangle)
+            planes[index_prim]->update (plane_reg);
+        }
+    }
+
+
+  if (regularize_coplanarity)
+    {
+      //detecting co-planarity and store in list_coplanar_prim
+      for (std::size_t i = 0; i < clusters.size(); ++ i)
+        {
+          Vector vec_reg = clusters[i].normal;
+
+          for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip)
+            clusters[i].coplanar_group.push_back (static_cast<std::size_t>(-1));
+
+          std::size_t cop_index=0;
+
+          for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j)
+            {
+              std::size_t index_prim = clusters[i].planes[j];
+
+              if (clusters[i].coplanar_group[j] == static_cast<std::size_t>(-1))
+                {
+                  clusters[i].coplanar_group[j] = cop_index;
+			
+                  Point pt_reg = planes[index_prim]->projection(centroids[index_prim]);
+                  Plane plan_reg(pt_reg,vec_reg);
+
+                  for (std::size_t k = j + 1; k < clusters[i].planes.size(); ++ k)
+                    {
+                      if (clusters[i].coplanar_group[k] == static_cast<std::size_t>(-1))
+                        {
+                          std::size_t index_prim_next = clusters[i].planes[k];
+                          Point pt_reg_next = planes[index_prim_next]->projection(centroids[index_prim_next]);
+                          Point pt_proj=plan_reg.projection(pt_reg_next);
+                          FT distance = std::sqrt (CGAL::squared_distance(pt_reg_next,pt_proj));
+					
+                          if (distance < tolerance_coplanarity)
+                            clusters[i].coplanar_group[k] = cop_index;
+                        }
+                    }
+                  cop_index++; 
+                }
+            }
+
+          //regularize primitive position by computing barycenter of cplanar planes
+          std::vector<Point> pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.));
+          std::vector<FT> area (cop_index, 0.);
+      
+          for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j)
+            {
+              std::size_t index_prim = clusters[i].planes[j];
+              std::size_t group = clusters[i].coplanar_group[j];
+              
+              Point pt_reg = planes[index_prim]->projection(centroids[index_prim]);
+
+              pt_bary[group] = CGAL::barycenter (pt_bary[group], area[group], pt_reg, areas[index_prim]); 
+              area[group] += areas[index_prim];
+            }
+
+
+          for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j)
+            {
+              std::size_t index_prim = clusters[i].planes[j];
+              std::size_t group = clusters[i].coplanar_group[j];
+              
+              Plane plane_reg (pt_bary[group], vec_reg);
+              
+              if (planes[index_prim]->plane_normal ()
+                  * plane_reg.orthogonal_vector() < 0)
+                planes[index_prim]->update (plane_reg.opposite());
+              else
+                planes[index_prim]->update (plane_reg);
+            }
+        }
+    } 
+}
+
+
+} // namespace CGAL
+
+#endif // CGAL_REGULARIZE_PLANES_H

Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp

@@ -5,6 +5,7 @@
 #include <CGAL/Simple_cartesian.h>
 
 #include <CGAL/Shape_detection_3.h>
+#include <CGAL/regularize_planes.h>
 #include <CGAL/Point_with_normal_3.h>
 #include <CGAL/property_map.h>
 
@@ -126,6 +127,10 @@ bool test_scene() {
     return false;
   }
 
+  // Test regularization
+  CGAL::regularize_planes (ransac, true, true, true, true,
+                           (FT)50., (FT)0.01f);
+  
   Point_index_range pts = ransac.indices_of_unassigned_points();
 
   std::cout << " succeeded" << std::endl;

Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.cpp

@@ -8,6 +8,7 @@
 #include <CGAL/Random.h>
 
 #include <CGAL/Shape_detection_3.h>
+#include <CGAL/Plane_regularization.h>
 #include <CGAL/Delaunay_triangulation_2.h>
 #include <CGAL/Alpha_shape_2.h>
 
@@ -38,6 +39,14 @@ class Polyhedron_demo_point_set_shape_detection_plugin :
     Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
     QAction* actionDetect;
 
+  typedef CGAL::Identity_property_map<Point_set::Point_with_normal> PointPMap;
+  typedef CGAL::Normal_of_point_with_normal_pmap<Point_set::Geom_traits> NormalPMap;
+
+  typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits<Epic_kernel, Point_set, PointPMap, NormalPMap> Traits;
+  typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits> Shape_detection;
+  typedef CGAL::Plane_regularization<Traits>                   Regularization;
+
+  
 public:
   void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface) {
     actionDetect = new QAction(tr("Point Set Shape Detection"), mainWindow);
@@ -65,7 +74,7 @@ private:
 
   typedef Kernel::Plane_3 Plane_3;
   
-  void build_alpha_shape (Point_set& points, const Plane_3& plane,
+  void build_alpha_shape (Point_set& points, boost::shared_ptr<CGAL::Shape_detection_3::Plane<Traits> > plane,
                           Scene_polyhedron_item* item, double epsilon);
 
 }; // end Polyhedron_demo_point_set_shape_detection_plugin
@@ -93,6 +102,7 @@ public:
   bool detect_cone() const { return coneCB->isChecked(); }
   bool generate_alpha() const { return m_generate_alpha->isChecked(); }
   bool generate_subset() const { return !(m_do_not_generate_subset->isChecked()); }
+  bool regularize() const { return m_regularize->isChecked(); }
 };
 
 void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered() {
@@ -123,12 +133,7 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered
 
     QApplication::setOverrideCursor(Qt::WaitCursor);
 
-    typedef CGAL::Identity_property_map<Point_set::Point_with_normal> PointPMap;
-    typedef CGAL::Normal_of_point_with_normal_pmap<Point_set::Geom_traits> NormalPMap;
-
-    typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits<Epic_kernel, Point_set, PointPMap, NormalPMap> Traits;
-    typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits> Shape_detection;
-
+    
     Shape_detection shape_detection;
     shape_detection.set_input(*points);
 
@@ -161,6 +166,19 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered
     shape_detection.detect(op);
 
     std::cout << shape_detection.shapes().size() << " shapes found" << std::endl;
+
+    if (dialog.regularize ())
+      {
+        std::cerr << "Regularization of planes... " << std::endl;
+        Regularization regularization (*points, shape_detection);
+
+        regularization.run (180 * std::acos (op.normal_threshold) / M_PI, op.epsilon);
+    
+        std::cerr << "done" << std::endl;
+      }
+
+    std::map<Kernel::Point_3, QColor> color_map;
+    
     //print_message(QString("%1 shapes found.").arg(shape_detection.number_of_shapes()));
     int index = 0;
     BOOST_FOREACH(boost::shared_ptr<Shape_detection::Shape> shape, shape_detection.shapes())
@@ -178,9 +196,11 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered
         point_item->point_set()->push_back((*points)[i]);
       
       unsigned char r, g, b;
+
       r = static_cast<unsigned char>(64 + rand.get_int(0, 192));
       g = static_cast<unsigned char>(64 + rand.get_int(0, 192));
       b = static_cast<unsigned char>(64 + rand.get_int(0, 192));
+
       point_item->setRbgColor(r, g, b);
 
       // Providing a useful name consisting of the order of detection, name of type and number of inliers
@@ -194,16 +214,33 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered
         {
           ss << item->name().toStdString() << "_plane_";
 
+          boost::shared_ptr<CGAL::Shape_detection_3::Plane<Traits> > pshape
+            = boost::dynamic_pointer_cast<CGAL::Shape_detection_3::Plane<Traits> > (shape);
+          
+          Kernel::Point_3 ref = CGAL::ORIGIN + pshape->plane_normal ();
+
+          if (color_map.find (ref) == color_map.end ())
+            {
+              ref = CGAL::ORIGIN + (-1.) * pshape->plane_normal ();
+              if (color_map.find (ref) == color_map.end ())
+                color_map[ref] = point_item->color ();
+              else
+                point_item->setColor (color_map[ref]);
+            }
+          else
+            point_item->setColor (color_map[ref]);
+
+          ss << "(" << ref << ")_";
+      
           if (dialog.generate_alpha ())
             {
               // If plane, build alpha shape
               Scene_polyhedron_item* poly_item = new Scene_polyhedron_item;
 
-              Plane_3 plane = (Plane_3)(*(dynamic_cast<CGAL::Shape_detection_3::Plane<Traits>*>(shape.get ())));
-              build_alpha_shape (*(point_item->point_set()), plane,
+              build_alpha_shape (*(point_item->point_set()), pshape,
                                  poly_item, dialog.cluster_epsilon());
           
-              poly_item->setRbgColor(r-32, g-32, b-32);
+              poly_item->setColor(point_item->color ());
               poly_item->setName(QString("%1%2_alpha_shape").arg(QString::fromStdString(ss.str()))
                                  .arg (QString::number (shape->indices_of_assigned_points().size())));
               poly_item->setRenderingMode (Flat);
@@ -250,7 +287,8 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered
 }
 
 void Polyhedron_demo_point_set_shape_detection_plugin::build_alpha_shape
-(Point_set& points, const Plane_3& plane, Scene_polyhedron_item* item, double epsilon)
+(Point_set& points,  boost::shared_ptr<CGAL::Shape_detection_3::Plane<Traits> > plane,
+ Scene_polyhedron_item* item, double epsilon)
 {
   typedef Kernel::Point_2  Point_2;
   typedef CGAL::Alpha_shape_vertex_base_2<Kernel> Vb;
@@ -264,7 +302,7 @@ void Polyhedron_demo_point_set_shape_detection_plugin::build_alpha_shape
   projections.reserve (points.size ());
 
   for (std::size_t i = 0; i < points.size (); ++ i)
-    projections.push_back (plane.to_2d (points[i]));
+    projections.push_back (plane->to_2d (points[i]));
 
   Alpha_shape_2 ashape (projections.begin (), projections.end (), epsilon);
   
@@ -286,7 +324,7 @@ void Polyhedron_demo_point_set_shape_detection_plugin::build_alpha_shape
           if (map_v2i.find (it->vertex (i)) == map_v2i.end ())
             {
               map_v2i.insert (std::make_pair (it->vertex (i), current_index ++));
-              Point p = plane.to_3d (it->vertex (i)->point ());
+              Point p = plane->to_3d (it->vertex (i)->point ());
               soup_item->new_vertex (p.x (), p.y (), p.z ());
             }
         }
@@ -297,6 +335,14 @@ void Polyhedron_demo_point_set_shape_detection_plugin::build_alpha_shape
 
   soup_item->orient();
   soup_item->exportAsPolyhedron (item->polyhedron());
+
+  if (soup_item->isEmpty ())
+    {
+      std::cerr << "POLYGON SOUP EMPTY" << std::endl;
+      for (std::size_t i = 0; i < projections.size (); ++ i)
+        std::cerr << projections[i] << std::endl;
+      
+    }
   
   delete soup_item;
 }

Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_shape_detection_plugin.ui

@@ -211,6 +211,16 @@
      </property>
     </widget>
    </item>
+   <item>
+    <widget class="QCheckBox" name="m_regularize">
+     <property name="text">
+      <string>Regularize planes</string>
+     </property>
+     <property name="checked">
+      <bool>true</bool>
+     </property>
+    </widget>
+   </item>
    <item>
     <widget class="QDialogButtonBox" name="buttonBox">
      <property name="orientation">