From cea079532ae869afbcc987a3d70b7fe26f4b4fdc Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 30 Oct 2015 11:35:10 +0100 Subject: [PATCH 01/48] Plane regularization file and basic structure --- .../include/CGAL/Plane_regularization.h | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 Point_set_shape_detection_3/include/CGAL/Plane_regularization.h diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h new file mode 100644 index 00000000000..36da23911f5 --- /dev/null +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -0,0 +1,63 @@ +// 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) : +// + +/** +* \ingroup PkgPointSetShapeDetection3 +* \file CGAL/Plane_regularization.h +* +*/ + + +#ifndef CGAL_PLANE_REGULARIZATION_H +#define CGAL_PLANE_REGULARIZATION_H + +template +class Plane_regularization +{ +public: + + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point; + typedef typename GeomTraits::Vector_3 Vector; + typedef typename GeomTraits::Plane_3 Plane; + +private: + +public: + + Plane_regularization () + { + + } + + virtual ~Plane_regularization () + { + + } + + + + +}; + + + +#endif // CGAL_PLANE_REGULARIZATION_H From 4c601e00db2f4e9966c59663e9914cb102fa39a0 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 30 Oct 2015 13:00:21 +0100 Subject: [PATCH 02/48] Plug RANSAC in plane regularization class --- .../include/CGAL/Plane_regularization.h | 69 +++++++++++++++++-- 1 file changed, 63 insertions(+), 6 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 36da23911f5..202c0d4735f 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -29,35 +29,92 @@ #ifndef CGAL_PLANE_REGULARIZATION_H #define CGAL_PLANE_REGULARIZATION_H -template +#include + +#include + + +namespace CGAL { + +template class Plane_regularization { public: - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Point_3 Point; - typedef typename GeomTraits::Vector_3 Vector; - typedef typename GeomTraits::Plane_3 Plane; + typedef Plane_regularization Self; + + typedef typename Traits::FT FT; + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Vector; + typedef typename Traits::Plane_3 Plane; + + typedef typename Traits::Point_map Point_map; + typedef typename Traits::Normal_map Normal_map; + typedef typename Traits::Input_range Input_range; + typedef typename Input_range::iterator Input_iterator; + + typedef Shape_detection_3::Shape_base Shape; + typedef Shape_detection_3::Plane Plane_shape; + private: + Traits m_traits; + + Input_iterator m_input_begin; + Input_iterator m_input_end; + Point_map m_point_pmap; + Normal_map m_normal_pmap; + + std::vector > m_planes; + public: - Plane_regularization () + Plane_regularization (Traits t = Traits ()) + : m_traits (t) { } + Plane_regularization (Input_range& input_range, + Shape_detection_3::Efficient_RANSAC& shape_detection) + : m_traits (shape_detection.traits()) + { + m_input_begin = input_range.begin (); + m_input_end = input_range.end (); + + BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) + { + boost::shared_ptr pshape + = boost::dynamic_pointer_cast(shape); + + // Ignore all shapes other than plane + if (!pshape) + continue; + + m_planes.push_back (pshape); + + } + + } + virtual ~Plane_regularization () { + clear (); + } + void clear () + { + std::vector > ().swap (m_planes); } + }; +}; // namespace CGAL #endif // CGAL_PLANE_REGULARIZATION_H From 49017ae7b4a4a0490c9d7a9ea729771a55fbb453 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 30 Oct 2015 13:22:12 +0100 Subject: [PATCH 03/48] Commented structure of algorithm + two useful functions --- .../include/CGAL/Plane_regularization.h | 65 +++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 202c0d4735f..8ba247f5570 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -108,7 +108,72 @@ public: std::vector > ().swap (m_planes); } + + void run () + { + + } + + Vector regularize_normal (const Vector& n, FT cos_vertical) + { + FT A = 1 - cos_vertical * cos_vertical; + FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); + + FT vx = std::sqrt (A/B); + + if (n.x() < 0) + vx = -vx; + + FT vy = vx * (n.y() / n.x()); + + Vector res (vx, vy, cos_vertical); + + return res / std::sqrt (res * res); + } + + Vector regularize_normals_from_prior (const Vector& np, + const Vector& n, + FT cos_vertical) + { + FT vx, vy; + + if (np.x() != 0) + { + FT a = (np.y() * np.y()) / (np.x() * np.x()) + 1; + FT b = 2 * np.y() * np.z() * cos_vertical / np.x(); + FT c= cos_vertical * cos_vertical-1; + + if (4 * a * c > b * b) + return regularize_normal (n, cos_vertical); + else + { + FT delta = std::sqrt (b * b-4 * a * c); + FT vy1= (-b-delta) / (2 * a); + FT vy2= (-b+delta) / (2 * a); + + vy = (std::fabs(n.y()-vy1) < std::fabs(n.y()-vy2)) + ? vy1 : vy2; + + vx = (-np.y() * vy-np.z() * cos_vertical) / np.x(); + } + } + else if (np.y() != 0) + { + vy = -np.z() * cos_vertical / np.y(); + vx = std::sqrt (1 - cos_vertical * cos_vertical - vy * vy); + + if (n.x() < 0) + vx = -vx; + } + else + return regularize_normal (n, cos_vertical); + + Vector_3 res (vx, vy, cos_vertical); + FT norm = std::max(1e-5, 1. / sqrt(res.squared_length ())); + + return norm * res; + } From 971f445093a9f8b938a4592b065f4d4e7a177f06 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 30 Oct 2015 14:37:50 +0100 Subject: [PATCH 04/48] Dump raw research code (wrapper not complete yet) --- .../include/CGAL/Plane_regularization.h | 690 +++++++++++++++--- 1 file changed, 584 insertions(+), 106 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 8ba247f5570..4ecc3420121 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -20,10 +20,10 @@ // /** -* \ingroup PkgPointSetShapeDetection3 -* \file CGAL/Plane_regularization.h -* -*/ + * \ingroup PkgPointSetShapeDetection3 + * \file CGAL/Plane_regularization.h + * + */ #ifndef CGAL_PLANE_REGULARIZATION_H @@ -36,148 +36,626 @@ namespace CGAL { -template -class Plane_regularization -{ -public: + template + class Plane_regularization + { + public: - typedef Plane_regularization Self; + typedef Plane_regularization Self; - typedef typename Traits::FT FT; - typedef typename Traits::Point_3 Point; - typedef typename Traits::Vector_3 Vector; - typedef typename Traits::Plane_3 Plane; + typedef typename Traits::FT FT; + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Vector; + typedef typename Traits::Plane_3 Plane; - typedef typename Traits::Point_map Point_map; - typedef typename Traits::Normal_map Normal_map; - typedef typename Traits::Input_range Input_range; - typedef typename Input_range::iterator Input_iterator; + typedef typename Traits::Point_map Point_map; + typedef typename Traits::Normal_map Normal_map; + typedef typename Traits::Input_range Input_range; + typedef typename Input_range::iterator Input_iterator; - typedef Shape_detection_3::Shape_base Shape; - typedef Shape_detection_3::Plane Plane_shape; + typedef Shape_detection_3::Shape_base Shape; + typedef Shape_detection_3::Plane Plane_shape; -private: + private: - Traits m_traits; + struct HPoint { + Point position; + bool is_kept; + bool is_fitted_to_plane; + // bool is_fitted_to_plane_bordering_point; + bool is_edge; + bool is_corner; + int primitive_index; + std::vector < int > plane_index_list; + std::vector < int > list_camera_index; + float accuracy; + float weight; + int weight_visibility; - Input_iterator m_input_begin; - Input_iterator m_input_end; - Point_map m_point_pmap; - Normal_map m_normal_pmap; + }; - std::vector > m_planes; + Traits m_traits; -public: + Input_iterator m_input_begin; + Input_iterator m_input_end; + Point_map m_point_pmap; + Normal_map m_normal_pmap; - Plane_regularization (Traits t = Traits ()) - : m_traits (t) - { + std::vector > m_planes; - } + public: - Plane_regularization (Input_range& input_range, - Shape_detection_3::Efficient_RANSAC& shape_detection) - : m_traits (shape_detection.traits()) - { - m_input_begin = input_range.begin (); - m_input_end = input_range.end (); + Plane_regularization (Traits t = Traits ()) + : m_traits (t) + { - BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) - { - boost::shared_ptr pshape - = boost::dynamic_pointer_cast(shape); + } + + Plane_regularization (Input_range& input_range, + const Shape_detection_3::Efficient_RANSAC& shape_detection) + : m_traits (shape_detection.traits()) + { + m_input_begin = input_range.begin (); + m_input_end = input_range.end (); + + BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) + { + boost::shared_ptr pshape + = boost::dynamic_pointer_cast(shape); - // Ignore all shapes other than plane - if (!pshape) - continue; + // Ignore all shapes other than plane + if (pshape == boost::shared_ptr()) + continue; - m_planes.push_back (pshape); + m_planes.push_back (pshape); - } + } - } + } - virtual ~Plane_regularization () - { - clear (); - } + virtual ~Plane_regularization () + { + clear (); + } - void clear () - { - std::vector > ().swap (m_planes); - } + void clear () + { + std::vector > ().swap (m_planes); + } - void run () - { + void run (FT epsilon, FT tolerance_coplanarity) + { + // WRAPPER BEGIN + std::vector extracted_planes; + std::vector list_centroid; + std::vector list_areas; + std::vector > plane_point_index; + std::vector HPS; + std::vector label_plane; + + for (std::size_t i = 0; i < m_planes.size (); ++ i) + { + extracted_planes.push_back (static_cast (*(m_planes[i]))); + plane_point_index.push_back (std::vector()); + std::copy (m_planes[i]->indices_of_assigned_points().begin (), + m_planes[i]->indices_of_assigned_points().begin (), + std::back_inserter (plane_point_index.back ())); + } - } - Vector regularize_normal (const Vector& n, FT cos_vertical) - { - FT A = 1 - cos_vertical * cos_vertical; - FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); - - FT vx = std::sqrt (A/B); - - if (n.x() < 0) - vx = -vx; - - FT vy = vx * (n.y() / n.x()); + // WRAPPER END + + // find pairs of epsilon-parallel primitives and store them in table_parallel + std::vector < std::vector < bool > > table_parallel; + for( std::size_t i=0;i table_parallel_tmp; + for( std::size_t j=0;j1.-epsilon && i!=j) table_parallel_tmp.push_back(true); + else table_parallel_tmp.push_back(false); + } + table_parallel.push_back(table_parallel_tmp); + } + - Vector res (vx, vy, cos_vertical); - return res / std::sqrt (res * res); - } + // clustering the parallel primitives and store them in list_parallel_planes + // & compute the normal, size and cos angle to the vertical of each cluster, and store that in list_cluster_normales, list_cluster_angle and list_cluster_area + std::vector < std::vector < int > > list_parallel_planes; + std::vector < Vector > list_cluster_normales; + std::vector < double > list_cluster_cosangle_vertical; + std::vector < double > list_cluster_area; + std::vector < bool > is_available; for( std::size_t i=0;i b * b) - return regularize_normal (n, cos_vertical); - else - { - FT delta = std::sqrt (b * b-4 * a * c); - FT vy1= (-b-delta) / (2 * a); - FT vy2= (-b+delta) / (2 * a); + //initialization containers + std::vector < int > index_container_parallel; index_container_parallel.push_back(i); + std::vector < int > index_container_former_ring_parallel; index_container_former_ring_parallel.push_back(i); + std::list < int > index_container_current_ring_parallel; - vy = (std::fabs(n.y()-vy1) < std::fabs(n.y()-vy2)) - ? vy1 : vy2; + //propagation over the pairs of epsilon-parallel primitives + bool propagation=true; + Vector cluster_normal=extracted_planes[i].orthogonal_vector(); + double cumulated_area=list_areas[i]; + + do{ + propagation=false; - vx = (-np.y() * vy-np.z() * cos_vertical) / np.x(); + for (std::size_t k=0;k1.-epsilon ){ + + propagation=true; + index_container_current_ring_parallel.push_back(it); + is_available[it]=false; + + if(cluster_normal*normal_it <0) normal_it=-normal_it; + cluster_normal=(FT)cumulated_area*cluster_normal+(FT)list_areas[it]*normal_it; + FT norm=1./sqrt(cluster_normal.squared_length()); + cluster_normal=norm*cluster_normal; + cumulated_area+=list_areas[it]; + } + } + } + + //update containers + index_container_former_ring_parallel.clear(); + for(std::list < int >::iterator it = index_container_current_ring_parallel.begin(); it != index_container_current_ring_parallel.end(); ++it){ + index_container_former_ring_parallel.push_back(*it); + index_container_parallel.push_back(*it); + } + index_container_current_ring_parallel.clear(); + + }while(propagation); + + list_parallel_planes.push_back(index_container_parallel); + list_cluster_normales.push_back(cluster_normal); + list_cluster_area.push_back(cumulated_area); + Vector v_vertical(0.,0.,1.); + list_cluster_cosangle_vertical.push_back(abs(v_vertical*cluster_normal)); + } + } + is_available.clear(); + + //discovery orthogonal relationship between clusters + std::vector < std::vector < bool > > group_planes_orthogonal; + for( std::size_t i=0;i gp_tmp; for( std::size_t j=0;j cosangle_centroids; + std::vector < int > list_cluster_index; for( std::size_t i=0;i1.-epsilon) cosangle_centroids[i]=1; + } + for (std::size_t i=0; i "; + for (std::size_t j=0; j "; + for (std::size_t j=0;j "< "< > subgraph_clusters; + std::vector < int > subgraph_clusters_max_area_index; + std::vector < bool > is_free; for (std::size_t i=0; i index_container; index_container.push_back(i); + std::vector < int > index_container_former_ring; index_container_former_ring.push_back(i); + std::list < int > index_container_current_ring; + + //propagation + bool propagation=true; + do{ + propagation=false; + + //neighbors + for (std::size_t k=0;k::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); + } + } + is_free.clear(); + + + //create subgraphs of mutually orthogonal clusters in which the largest cluster is excluded and store in subgraph_clusters_prop + std::vector < std::vector < int > > subgraph_clusters_prop; + for (std::size_t i=0;i subgraph_clusters_prop_temp; + for (std::size_t j=0;j cluster_is_available; + for( std::size_t i=0;i index_container; index_container.push_back(index_current); + std::vector < int > index_container_former_ring; index_container_former_ring.push_back(index_current); + std::list < int > index_container_current_ring; + + //propagation + bool propagation=true; + do{ + propagation=false; + + //neighbors + for (std::size_t k=0;k::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); + } + + + + //recompute optimal plane for each primitive after normal regularization + for (std::size_t i=0; i 1. - epsilon) extracted_planes[index_prim]=plane_reg; + } + } + + + + + + //detecting co-planarity and store in list_coplanar_prim + std::vector< std::vector< std::vector < int > > > list_coplanar_prim; + for (std::size_t i=0; i > list_coplanar_prim_tmp; + Vector vec_reg=list_cluster_normales[i]; + std::vector < int > list_cop_index; for( std::size_t ip=0;ip list_coplanar_prim_tmp_tmp; + list_cop_index[j]=cop_index; + list_coplanar_prim_tmp_tmp.push_back(index_prim); + + Point pt_reg=extracted_planes[index_prim].projection(list_centroid[index_prim]); + Plane plan_reg(pt_reg,vec_reg); + + for (std::size_t k=j+1; k > list_primitive_reg_index_extracted_planes; + std::vector < Plane > list_primitive_reg; + + for (std::size_t i=0;i list_primitive_reg_index_extracted_planes_tmp1; + for (std::size_t k=0; k 1. - epsilon){ + if(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()<0) extracted_planes[index_prim]=plane_reg.opposite(); + else extracted_planes[index_prim]=plane_reg; + is_reg_used=true; + list_primitive_reg_index_extracted_planes_tmp1.push_back(index_prim); + } + else{ + list_primitive_reg.push_back(extracted_planes[index_prim]); + std::vector< int > list_primitive_reg_index_extracted_planes_tmp; + list_primitive_reg_index_extracted_planes_tmp.push_back(index_prim); + list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp); + } + } + if(is_reg_used) { + list_primitive_reg.push_back(plane_reg); + list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp1); + } + } + } + + std::cout< "; @@ -300,7 +301,7 @@ namespace CGAL { } std::cout<<" -> "< "< 1. - epsilon) extracted_planes[index_prim]=plane_reg; + if( std::fabs(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()) > 1. - epsilon) extracted_planes[index_prim]=plane_reg; } } @@ -521,7 +522,7 @@ namespace CGAL { std::vector< int > list_primitive_reg_index_extracted_planes_tmp1; for (std::size_t k=0; k 1. - epsilon){ + if( std::fabs(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()) > 1. - epsilon){ if(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()<0) extracted_planes[index_prim]=plane_reg.opposite(); else extracted_planes[index_prim]=plane_reg; is_reg_used=true; @@ -558,7 +559,7 @@ namespace CGAL { plane_point_index_temp.push_back(plane_point_index[i]); for (std::size_t k=0; k Date: Mon, 2 Nov 2015 11:17:34 +0100 Subject: [PATCH 06/48] Improve code readability --- .../include/CGAL/Plane_regularization.h | 658 ++++++++++-------- 1 file changed, 383 insertions(+), 275 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index d20abcacf9f..bdd34a4cd6d 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -151,18 +151,22 @@ namespace CGAL { // find pairs of epsilon-parallel primitives and store them in table_parallel std::vector < std::vector < bool > > table_parallel; - for( std::size_t i=0;i table_parallel_tmp; - for( std::size_t j=0;j table_parallel_tmp; + for( std::size_t j=0;j1.-epsilon && i!=j) table_parallel_tmp.push_back(true); - else table_parallel_tmp.push_back(false); + if (std::fabs(v1*v2)>1.-epsilon && i!=j) + table_parallel_tmp.push_back(true); + else + table_parallel_tmp.push_back(false); + } + table_parallel.push_back(table_parallel_tmp); } - table_parallel.push_back(table_parallel_tmp); - } @@ -172,81 +176,110 @@ namespace CGAL { std::vector < Vector > list_cluster_normales; std::vector < double > list_cluster_cosangle_vertical; std::vector < double > list_cluster_area; - std::vector < bool > is_available; for( std::size_t i=0;i is_available; + for( std::size_t i=0;i index_container_parallel; index_container_parallel.push_back(i); - std::vector < int > index_container_former_ring_parallel; index_container_former_ring_parallel.push_back(i); - std::list < int > index_container_current_ring_parallel; + //initialization containers + std::vector < int > index_container_parallel; + index_container_parallel.push_back(i); + + std::vector < int > index_container_former_ring_parallel; + index_container_former_ring_parallel.push_back(i); + + std::list < int > index_container_current_ring_parallel; - //propagation over the pairs of epsilon-parallel primitives - bool propagation=true; - Vector cluster_normal=extracted_planes[i].orthogonal_vector(); - double cumulated_area=list_areas[i]; + //propagation over the pairs of epsilon-parallel primitives + bool propagation=true; + Vector cluster_normal=extracted_planes[i].orthogonal_vector(); + double cumulated_area=list_areas[i]; - do{ - propagation=false; + do + { + propagation=false; - for (std::size_t k=0;k1.-epsilon ){ + if( table_parallel[plane_index][it] && is_available[it] + && std::fabs(normal_it*cluster_normal)>1.-epsilon ) + { - propagation=true; - index_container_current_ring_parallel.push_back(it); - is_available[it]=false; + propagation=true; + index_container_current_ring_parallel.push_back(it); + is_available[it]=false; - if(cluster_normal*normal_it <0) normal_it=-normal_it; - cluster_normal=(FT)cumulated_area*cluster_normal+(FT)list_areas[it]*normal_it; - FT norm=1./sqrt(cluster_normal.squared_length()); - cluster_normal=norm*cluster_normal; - cumulated_area+=list_areas[it]; - } - } - } + if(cluster_normal*normal_it <0) + normal_it=-normal_it; + + cluster_normal=(FT)cumulated_area*cluster_normal+(FT)list_areas[it]*normal_it; + FT norm=1./sqrt(cluster_normal.squared_length()); + cluster_normal=norm*cluster_normal; + cumulated_area+=list_areas[it]; + } + } + } - //update containers - index_container_former_ring_parallel.clear(); - for(std::list < int >::iterator it = index_container_current_ring_parallel.begin(); it != index_container_current_ring_parallel.end(); ++it){ - index_container_former_ring_parallel.push_back(*it); - index_container_parallel.push_back(*it); - } - index_container_current_ring_parallel.clear(); + //update containers + index_container_former_ring_parallel.clear(); + for (std::list < int >::iterator it = index_container_current_ring_parallel.begin(); + it != index_container_current_ring_parallel.end(); ++it) + { + index_container_former_ring_parallel.push_back(*it); + index_container_parallel.push_back(*it); + } + index_container_current_ring_parallel.clear(); - }while(propagation); + } + while(propagation); - list_parallel_planes.push_back(index_container_parallel); - list_cluster_normales.push_back(cluster_normal); - list_cluster_area.push_back(cumulated_area); - Vector v_vertical(0.,0.,1.); - list_cluster_cosangle_vertical.push_back(std::fabs(v_vertical*cluster_normal)); + list_parallel_planes.push_back(index_container_parallel); + list_cluster_normales.push_back(cluster_normal); + list_cluster_area.push_back(cumulated_area); + Vector v_vertical(0.,0.,1.); + list_cluster_cosangle_vertical.push_back(std::fabs(v_vertical*cluster_normal)); + } } - } is_available.clear(); //discovery orthogonal relationship between clusters std::vector < std::vector < bool > > group_planes_orthogonal; - for( std::size_t i=0;i gp_tmp; for( std::size_t j=0;j gp_tmp; + for( std::size_t j=0;j cosangle_centroids; - std::vector < int > list_cluster_index; for( std::size_t i=0;i list_cluster_index; + for( std::size_t i=0;i1.-epsilon) cosangle_centroids[i]=1; - } - for (std::size_t i=0; i1.-epsilon) + cosangle_centroids[i]=1; + } + for (std::size_t i=0; i "; - for (std::size_t j=0; j "; + for (std::size_t j=0; j "; - for (std::size_t j=0;j "; + for (std::size_t j=0;j "< "< "< "< > subgraph_clusters; std::vector < int > subgraph_clusters_max_area_index; - std::vector < bool > is_free; for (std::size_t i=0; i is_free; + for (std::size_t i=0; i index_container; + index_container.push_back(i); + std::vector < int > index_container_former_ring; + index_container_former_ring.push_back(i); + std::list < int > index_container_current_ring; - //initialization containers - std::vector < int > index_container; index_container.push_back(i); - std::vector < int > index_container_former_ring; index_container_former_ring.push_back(i); - std::list < int > index_container_current_ring; + //propagation + bool propagation=true; + do + { + propagation=false; - //propagation - bool propagation=true; - do{ - propagation=false; + //neighbors + for (std::size_t k=0;k::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); } - - //update containers - index_container_former_ring.clear(); - for(std::list < int >::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); } - } is_free.clear(); //create subgraphs of mutually orthogonal clusters in which the largest cluster is excluded and store in subgraph_clusters_prop std::vector < std::vector < int > > subgraph_clusters_prop; - for (std::size_t i=0;i subgraph_clusters_prop_temp; - for (std::size_t j=0;j subgraph_clusters_prop_temp; + for (std::size_t j=0;j cluster_is_available; - for( std::size_t i=0;i index_container; index_container.push_back(index_current); - std::vector < int > index_container_former_ring; index_container_former_ring.push_back(index_current); - std::list < int > index_container_current_ring; + //initialization containers + std::vector < int > index_container; + index_container.push_back(index_current); + std::vector < int > index_container_former_ring; + index_container_former_ring.push_back(index_current); + std::list < int > index_container_current_ring; - //propagation - bool propagation=true; - do{ - propagation=false; + //propagation + bool propagation=true; + do + { + propagation=false; - //neighbors - for (std::size_t k=0;k::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); - } + //update containers + index_container_former_ring.clear(); + for(std::list < int >::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); + } //recompute optimal plane for each primitive after normal regularization - for (std::size_t i=0; i 1. - epsilon) extracted_planes[index_prim]=plane_reg; + if( std::fabs(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()) > 1. - epsilon) + extracted_planes[index_prim]=plane_reg; + } } - } @@ -455,45 +537,54 @@ namespace CGAL { //detecting co-planarity and store in list_coplanar_prim std::vector< std::vector< std::vector < int > > > list_coplanar_prim; - for (std::size_t i=0; i > list_coplanar_prim_tmp; - Vector vec_reg=list_cluster_normales[i]; - std::vector < int > list_cop_index; for( std::size_t ip=0;ip > list_coplanar_prim_tmp; + Vector vec_reg=list_cluster_normales[i]; + std::vector < int > list_cop_index; + for( std::size_t ip=0;ip list_coplanar_prim_tmp_tmp; - list_cop_index[j]=cop_index; - list_coplanar_prim_tmp_tmp.push_back(index_prim); + std::vector < int > list_coplanar_prim_tmp_tmp; + list_cop_index[j]=cop_index; + list_coplanar_prim_tmp_tmp.push_back(index_prim); - Point pt_reg=extracted_planes[index_prim].projection(list_centroid[index_prim]); - Plane plan_reg(pt_reg,vec_reg); + Point pt_reg=extracted_planes[index_prim].projection(list_centroid[index_prim]); + Plane plan_reg(pt_reg,vec_reg); - for (std::size_t k=j+1; k > list_primitive_reg_index_extracted_planes; std::vector < Plane > list_primitive_reg; - for (std::size_t i=0;i list_primitive_reg_index_extracted_planes_tmp1; - for (std::size_t k=0; k 1. - epsilon){ - if(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()<0) extracted_planes[index_prim]=plane_reg.opposite(); - else extracted_planes[index_prim]=plane_reg; - is_reg_used=true; - list_primitive_reg_index_extracted_planes_tmp1.push_back(index_prim); + bool is_reg_used=false; + std::vector< int > list_primitive_reg_index_extracted_planes_tmp1; + + for (std::size_t k=0; k 1. - epsilon) + { + if(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()<0) + extracted_planes[index_prim]=plane_reg.opposite(); + else + extracted_planes[index_prim]=plane_reg; + is_reg_used=true; + list_primitive_reg_index_extracted_planes_tmp1.push_back(index_prim); + } + else{ + list_primitive_reg.push_back(extracted_planes[index_prim]); + std::vector< int > list_primitive_reg_index_extracted_planes_tmp; + list_primitive_reg_index_extracted_planes_tmp.push_back(index_prim); + list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp); + } + } + if(is_reg_used) { + list_primitive_reg.push_back(plane_reg); + list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp1); + } } - else{ - list_primitive_reg.push_back(extracted_planes[index_prim]); - std::vector< int > list_primitive_reg_index_extracted_planes_tmp; - list_primitive_reg_index_extracted_planes_tmp.push_back(index_prim); - list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp); - } - } - if(is_reg_used) { - list_primitive_reg.push_back(plane_reg); - list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp1); - } } - } std::cout<itemChanged(index); diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.ui b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.ui index e728862d1dd..6c236610d04 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.ui +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.ui @@ -7,7 +7,7 @@ 0 0 444 - 205 + 243 @@ -197,6 +197,16 @@ + + + + Regularize planes + + + true + + + From 455e341731edc4eebb295a1a9c1f60293415d89b Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 2 Nov 2015 11:36:10 +0100 Subject: [PATCH 08/48] Bugfix in wrapper for regularization --- .../include/CGAL/Plane_regularization.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index bdd34a4cd6d..ea5628a8e7f 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -126,7 +126,7 @@ namespace CGAL { extracted_planes.push_back (static_cast (*(m_planes[i]))); plane_point_index.push_back (std::vector()); std::copy (m_planes[i]->indices_of_assigned_points().begin (), - m_planes[i]->indices_of_assigned_points().begin (), + m_planes[i]->indices_of_assigned_points().end (), std::back_inserter (plane_point_index.back ())); std::vector < Point > listp; @@ -227,7 +227,7 @@ namespace CGAL { if(cluster_normal*normal_it <0) normal_it=-normal_it; - + cluster_normal=(FT)cumulated_area*cluster_normal+(FT)list_areas[it]*normal_it; FT norm=1./sqrt(cluster_normal.squared_length()); cluster_normal=norm*cluster_normal; @@ -343,7 +343,7 @@ namespace CGAL { for (std::size_t i=0; i "; - for (std::size_t j=0;j Date: Mon, 2 Nov 2015 15:15:27 +0100 Subject: [PATCH 09/48] Add some useful methods for RANSAC Planes --- .../include/CGAL/Shape_detection_3/Plane.h | 20 ++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h index 64af70b6aef..10f13fd75f6 100644 --- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h +++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h @@ -85,7 +85,25 @@ namespace CGAL { return d * d; } - + /*! + Computes the orthogonal projection of a query point on the shape. + */ + Point_3 projection (const Point_3 &p) const { + Vector_3 v (m_point_on_primitive, p); + return m_point_on_primitive + (v * m_base1) * m_base1 + (v * m_base2) * m_base2; + } + + /*! + Replaces the plane by p + */ + void update (const Plane_3& p) { + m_base1 = p.base1 (); m_base2 = p.base2 (); m_normal = 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. From a49b440ba03e42095ec2b3d28a05029d7cb99eec Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 2 Nov 2015 15:15:40 +0100 Subject: [PATCH 10/48] Use RANSAC planes directly (reduce wrapper) --- .../include/CGAL/Plane_regularization.h | 155 +++++++++--------- 1 file changed, 78 insertions(+), 77 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index ea5628a8e7f..5a82b2cbac7 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -111,10 +111,9 @@ namespace CGAL { - void run (FT epsilon, FT tolerance_coplanarity) + std::size_t run (FT epsilon, FT tolerance_coplanarity) { // WRAPPER BEGIN - std::vector extracted_planes; std::vector > plane_point_index; std::vector primitive_index (m_input_end - m_input_begin, -1); std::vector label_plane (m_input_end - m_input_begin, -1); @@ -123,7 +122,6 @@ namespace CGAL { for (std::size_t i = 0; i < m_planes.size (); ++ i) { - extracted_planes.push_back (static_cast (*(m_planes[i]))); plane_point_index.push_back (std::vector()); std::copy (m_planes[i]->indices_of_assigned_points().begin (), m_planes[i]->indices_of_assigned_points().end (), @@ -147,25 +145,25 @@ namespace CGAL { // WRAPPER END - + // std::cout< "< "<projection (list_centroid[index_prim]); + if( m_planes[index_prim]->plane_normal () * vec_reg < 0) vec_reg=-vec_reg; Plane plane_reg(pt_reg,vec_reg); - if( std::fabs(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()) > 1. - epsilon) - extracted_planes[index_prim]=plane_reg; + if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - epsilon) + { + m_planes[index_prim]->update (plane_reg); + } } } @@ -559,7 +559,7 @@ namespace CGAL { list_cop_index[j]=cop_index; list_coplanar_prim_tmp_tmp.push_back(index_prim); - Point pt_reg=extracted_planes[index_prim].projection(list_centroid[index_prim]); + Point pt_reg = m_planes[index_prim]->projection(list_centroid[index_prim]); Plane plan_reg(pt_reg,vec_reg); for (std::size_t k=j+1; kprojection(list_centroid[index_prim_next]); Point pt_proj=plan_reg.projection(pt_reg_next); double distance=distance_Point(pt_reg_next,pt_proj); @@ -603,12 +603,12 @@ namespace CGAL { for (std::size_t k=0; kprojection(list_centroid[index_prim]); pt_bary=barycenter(pt_bary, area,pt_reg,list_areas[index_prim]); area+=list_areas[index_prim]; } - Vector vec_reg=extracted_planes[list_coplanar_prim[i][j][0]].orthogonal_vector(); + Vector vec_reg = m_planes[list_coplanar_prim[i][j][0]]->plane_normal (); Plane plane_reg(pt_bary,vec_reg); @@ -618,17 +618,17 @@ namespace CGAL { for (std::size_t k=0; k 1. - epsilon) + if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector()) > 1. - epsilon) { - if(extracted_planes[index_prim].orthogonal_vector()*plane_reg.orthogonal_vector()<0) - extracted_planes[index_prim]=plane_reg.opposite(); + if(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector()<0) + m_planes[index_prim]->update (plane_reg.opposite()); else - extracted_planes[index_prim]=plane_reg; + m_planes[index_prim]->update (plane_reg); is_reg_used=true; list_primitive_reg_index_extracted_planes_tmp1.push_back(index_prim); } else{ - list_primitive_reg.push_back(extracted_planes[index_prim]); + list_primitive_reg.push_back(static_cast (*(m_planes[index_prim]))); std::vector< int > list_primitive_reg_index_extracted_planes_tmp; list_primitive_reg_index_extracted_planes_tmp.push_back(index_prim); list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp); @@ -641,59 +641,60 @@ namespace CGAL { } } - std::cout< > subgraph_clusters; std::vector < int > subgraph_clusters_max_area_index; std::vector < bool > is_free; - for (std::size_t i=0; i index_container; @@ -366,10 +365,10 @@ namespace CGAL { index_container_current_ring.push_back(j); is_free[j]=false; - if(max_area cluster_is_available; - for( std::size_t i=0;i index_container; @@ -462,10 +461,10 @@ namespace CGAL { index_container_current_ring.push_back(j); cluster_is_available[j]=false; - Vector new_vect=regularize_normals_from_prior(list_cluster_normales[cluster_index], - list_cluster_normales[j], - list_cluster_cosangle_vertical[j]); - list_cluster_normales[j]=new_vect; + Vector new_vect=regularize_normals_from_prior(clusters[cluster_index].normal, + clusters[j].normal, + clusters[j].cosangle_vertical); + clusters[j].normal = new_vect; } } } @@ -485,14 +484,14 @@ namespace CGAL { //recompute optimal plane for each primitive after normal regularization - for (std::size_t i=0; iprojection (m_centroids[index_prim]); if( m_planes[index_prim]->plane_normal () * vec_reg < 0) vec_reg=-vec_reg; @@ -511,24 +510,24 @@ namespace CGAL { //detecting co-planarity and store in list_coplanar_prim std::vector< std::vector< std::vector < int > > > list_coplanar_prim; - for (std::size_t i=0; i > list_coplanar_prim_tmp; - Vector vec_reg=list_cluster_normales[i]; + Vector vec_reg = clusters[i].normal; std::vector < int > list_cop_index; - for( std::size_t ip=0;ip list_coplanar_prim_tmp_tmp; list_cop_index[j]=cop_index; list_coplanar_prim_tmp_tmp.push_back(index_prim); @@ -536,12 +535,11 @@ namespace CGAL { Point pt_reg = m_planes[index_prim]->projection(m_centroids[index_prim]); Plane plan_reg(pt_reg,vec_reg); - for (std::size_t k=j+1; kprojection(m_centroids[index_prim_next]); Point pt_proj=plan_reg.projection(pt_reg_next); double distance=distance_Point(pt_reg_next,pt_proj); From 536abe1f3f8358f0f895f333fb3f86cba5792d7d Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 3 Nov 2015 18:02:15 +0100 Subject: [PATCH 14/48] Reorganize and simplify code --- .../include/CGAL/Plane_regularization.h | 420 +++++++++--------- 1 file changed, 209 insertions(+), 211 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 82dde7185a4..b9d28a927c9 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -62,10 +62,13 @@ namespace CGAL { struct Plane_cluster { + bool is_free; std::vector planes; + std::vector orthogonal_clusters; Vector normal; FT cosangle_vertical; FT area; + FT cosangle_centroid; }; Traits m_traits; @@ -127,214 +130,44 @@ namespace CGAL { { compute_centroids_and_areas (); - // find pairs of epsilon-parallel primitives and store them in parallel_planes - - std::vector < std::vector < std::size_t > > parallel_planes (m_planes.size ()); - - for (std::size_t i = 0; i < m_planes.size (); ++ i) - { - Vector v1 = m_planes[i]->plane_normal (); - - for (std::size_t j = 0; j < m_planes.size(); ++ j) - { - if (i == j) - continue; - - Vector v2 = m_planes[i]->plane_normal (); - - if (std::fabs (v1 * v2) > 1. - epsilon) - parallel_planes[i].push_back (j); - } - } - - - // clustering the parallel primitives and store them in clusters // & compute the normal, size and cos angle to the vertical of each cluster - std::vector clusters; - std::vector is_available (m_planes.size (), true); - - for (std::size_t i = 0; i < m_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 index_container_former_ring_parallel; - index_container_former_ring_parallel.push_back(i); - - std::list index_container_current_ring_parallel; - - //propagation over the pairs of epsilon-parallel primitives - bool propagation=true; - clu.normal = m_planes[i]->plane_normal (); - clu.area = m_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 = m_planes[it]->plane_normal (); - - if(is_available[it] - && std::fabs (normal_it*clu.normal) > 1. - epsilon ) - { - 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)m_areas[it] * normal_it; - FT norm = 1. / std::sqrt (clu.normal.squared_length()); - clu.normal = norm * clu.normal; - clu.area += m_areas[it]; - } - } - } - - //update containers - index_container_former_ring_parallel.clear(); - for (std::list < int >::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); - - Vector v_vertical(0.,0.,1.); - clu.cosangle_vertical = std::fabs(v_vertical*clu.normal); - } - } - is_available.clear(); + compute_parallel_clusters (clusters, epsilon); //discovery orthogonal relationship between clusters - std::vector < std::vector < bool > > group_planes_orthogonal; for (std::size_t i = 0; i < clusters.size(); ++ i) { - std::vector gp_tmp; - for (std::size_t j = 0; j < clusters.size(); ++ j) - gp_tmp.push_back(false); - group_planes_orthogonal.push_back (gp_tmp); - } - - for (std::size_t i = 0; i < group_planes_orthogonal.size(); ++ i) - { - for (std::size_t j = 0; j < group_planes_orthogonal.size(); ++ j) + for (std::size_t j = i + 1; j < clusters.size(); ++ j) { - - if (i != j && std::fabs (clusters[i].normal * clusters[j].normal) < epsilon) - { - group_planes_orthogonal[i][j]=true; - group_planes_orthogonal[j][i]=true; - } - } - } - - - - - - //clustering the vertical cosangle and store their centroids in cosangle_centroids and the centroid index of each cluster in list_cluster_index - std::vector < double > cosangle_centroids; - std::vector < int > list_cluster_index; - for( std::size_t i = 0; i < clusters.size(); ++ i) - list_cluster_index.push_back(-1); - - int mean_index = 0; - for (std::size_t i = 0; i < clusters.size(); ++ i) - { - if(list_cluster_index[i]<0) - { - list_cluster_index[i] = mean_index; - double mean = clusters[i].area * clusters[i].cosangle_vertical; - double mean_area = clusters[i].area; - for (std::size_t j = i+1; j < clusters.size(); ++ j) + if (std::fabs (clusters[i].normal * clusters[j].normal) < epsilon) { - if (list_cluster_index[j] < 0 && std::fabs (clusters[j].cosangle_vertical - - mean / mean_area) < epsilon) - { - list_cluster_index[j] = mean_index; - mean_area += clusters[j].area; - mean += clusters[j].area * clusters[j].cosangle_vertical; - } + clusters[i].orthogonal_clusters.push_back (j); + clusters[j].orthogonal_clusters.push_back (i); } - ++ mean_index; - mean /= mean_area; - cosangle_centroids.push_back (mean); } } - - //desactive Z-verticalité - for (std::size_t i = 0; i < cosangle_centroids.size(); ++ i) - { - if (cosangle_centroids[i] < epsilon) - cosangle_centroids[i] = 0; - else if (cosangle_centroids[i] > 1. - epsilon) - cosangle_centroids[i] = 1; - } - for (std::size_t i = 0; i < group_planes_orthogonal.size(); ++ i) - clusters[i].cosangle_vertical = cosangle_centroids[list_cluster_index[i]]; - - //display console - /* - std::cout< "; - for (std::size_t j=0; j "; - for (std::size_t j=0;j "< "< > subgraph_clusters; std::vector < int > subgraph_clusters_max_area_index; - std::vector < bool > is_free; + for (std::size_t i = 0; i < clusters.size(); ++ i) - is_free.push_back(true); + clusters[i].is_free = true; + for (std::size_t i = 0; i < clusters.size(); ++ i) { - if(is_free[i]) + if(clusters[i].is_free) { - is_free[i]=false; + clusters[i].is_free = false; double max_area = clusters[i].area; int index_max_area = i; @@ -357,13 +190,13 @@ namespace CGAL { int cluster_index=index_container_former_ring[k]; - for (std::size_t j=0;j > subgraph_clusters_prop; for (std::size_t i=0;i cluster_is_available; for (std::size_t i = 0; i < clusters.size(); ++ i) - cluster_is_available.push_back(true); + clusters[i].is_free = true; for (std::size_t i = 0; i < subgraph_clusters_prop.size(); ++ i) { @@ -430,7 +257,7 @@ namespace CGAL { Vector vec_current=regularize_normal(clusters[index_current].normal, clusters[index_current].cosangle_vertical); clusters[index_current].normal = vec_current; - cluster_is_available[index_current] = false; + clusters[index_current].is_free = false; //initialization containers std::vector < int > index_container; @@ -451,15 +278,15 @@ namespace CGAL { int cluster_index=index_container_former_ring[k]; - for (std::size_t j=0;jindices_of_assigned_points().size()) / 100.); } } - + + void compute_parallel_clusters (std::vector& clusters, double epsilon) + { + // find pairs of epsilon-parallel primitives and store them in parallel_planes + std::vector > parallel_planes (m_planes.size ()); + for (std::size_t i = 0; i < m_planes.size (); ++ i) + { + Vector v1 = m_planes[i]->plane_normal (); + + for (std::size_t j = 0; j < m_planes.size(); ++ j) + { + if (i == j) + continue; + + Vector v2 = m_planes[i]->plane_normal (); + + if (std::fabs (v1 * v2) > 1. - epsilon) + parallel_planes[i].push_back (j); + } + } + + + std::vector is_available (m_planes.size (), true); + + for (std::size_t i = 0; i < m_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 index_container_former_ring_parallel; + index_container_former_ring_parallel.push_back(i); + + std::list index_container_current_ring_parallel; + + //propagation over the pairs of epsilon-parallel primitives + bool propagation=true; + clu.normal = m_planes[i]->plane_normal (); + clu.area = m_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 = m_planes[it]->plane_normal (); + + if(is_available[it] + && std::fabs (normal_it*clu.normal) > 1. - epsilon ) + { + 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)m_areas[it] * normal_it; + FT norm = 1. / std::sqrt (clu.normal.squared_length()); + clu.normal = norm * clu.normal; + clu.area += m_areas[it]; + } + } + } + + //update containers + index_container_former_ring_parallel.clear(); + for (std::list::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); + + Vector v_vertical(0.,0.,1.); + clu.cosangle_vertical = std::fabs(v_vertical*clu.normal); + } + } + is_available.clear(); + } + + void cluster_vertical_cosangles (std::vector& clusters, double epsilon) + { + std::vector < double > cosangle_centroids; + std::vector < int > list_cluster_index; + for( std::size_t i = 0; i < clusters.size(); ++ i) + list_cluster_index.push_back(-1); + + int mean_index = 0; + for (std::size_t i = 0; i < clusters.size(); ++ i) + { + if(list_cluster_index[i]<0) + { + list_cluster_index[i] = mean_index; + double mean = clusters[i].area * clusters[i].cosangle_vertical; + double mean_area = clusters[i].area; + + for (std::size_t j = i+1; j < clusters.size(); ++ j) + { + if (list_cluster_index[j] < 0 && std::fabs (clusters[j].cosangle_vertical - + mean / mean_area) < epsilon) + { + list_cluster_index[j] = mean_index; + mean_area += clusters[j].area; + mean += clusters[j].area * clusters[j].cosangle_vertical; + } + } + ++ mean_index; + mean /= mean_area; + cosangle_centroids.push_back (mean); + } + } + + for (std::size_t i = 0; i < cosangle_centroids.size (); ++ i) + std::cerr << cosangle_centroids[i] << ", "; + std::cerr < 1. - epsilon) + cosangle_centroids[i] = 1; + } + for (std::size_t i = 0; i < clusters.size(); ++ i) + clusters[i].cosangle_vertical = cosangle_centroids[list_cluster_index[i]]; + } + + FT distance_Point (const Point& a, const Point& b) { return std::sqrt (CGAL::squared_distance (a, b)); From 752d9f2d16915e0f82ebdd585da183f47b89f65c Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 4 Nov 2015 11:23:39 +0100 Subject: [PATCH 15/48] More reorganization --- .../include/CGAL/Plane_regularization.h | 340 +++++++++--------- 1 file changed, 170 insertions(+), 170 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index b9d28a927c9..868a515c1da 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -157,159 +157,8 @@ namespace CGAL { //find subgraphs of mutually orthogonal clusters (store index of //clusters in subgraph_clusters), and select the cluster of //largest area - std::vector < std::vector < int > > subgraph_clusters; - std::vector < int > 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; - double max_area = clusters[i].area; - int index_max_area = i; - - //initialization containers - std::vector < int > index_container; - index_container.push_back(i); - std::vector < int > index_container_former_ring; - index_container_former_ring.push_back(i); - std::list < int > index_container_current_ring; - - //propagation - bool propagation=true; - do - { - propagation=false; - - //neighbors - for (std::size_t k=0;k::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 < int > > subgraph_clusters_prop; - for (std::size_t i=0;i subgraph_clusters_prop_temp; - for (std::size_t j=0;j index_container; - index_container.push_back(index_current); - std::vector < int > index_container_former_ring; - index_container_former_ring.push_back(index_current); - std::list < int > index_container_current_ring; - - //propagation - bool propagation=true; - do - { - propagation=false; - - //neighbors - for (std::size_t k=0;k::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_mutually_orthogonal_clusters (clusters); + //recompute optimal plane for each primitive after normal regularization for (std::size_t i=0; i < clusters.size(); ++ i) { @@ -325,18 +174,13 @@ namespace CGAL { Plane plane_reg(pt_reg,vec_reg); if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - epsilon) - { - m_planes[index_prim]->update (plane_reg); - } + m_planes[index_prim]->update (plane_reg); } } - - - //detecting co-planarity and store in list_coplanar_prim - std::vector< std::vector< std::vector < int > > > list_coplanar_prim; + std::vector > > list_coplanar_prim; for (std::size_t i = 0; i < clusters.size(); ++ i) { @@ -387,7 +231,7 @@ namespace CGAL { - //regularize primitive position by computing barycenter of coplanar planes + //regularize primitive position by computing barycenter of cplanar planes std::vector < std::vector < int > > list_primitive_reg_index_extracted_planes; std::vector < Plane > list_primitive_reg; @@ -442,15 +286,8 @@ namespace CGAL { // std::cout<& clusters) + { + std::vector < std::vector < int > > subgraph_clusters; + std::vector < int > 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; + double max_area = clusters[i].area; + int index_max_area = i; + + //initialization containers + std::vector < int > index_container; + index_container.push_back(i); + std::vector < int > index_container_former_ring; + index_container_former_ring.push_back(i); + std::list < int > index_container_current_ring; + + //propagation + bool propagation=true; + do + { + propagation=false; + + //neighbors + for (std::size_t k=0;k::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 < int > > subgraph_clusters_prop; + for (std::size_t i=0;i subgraph_clusters_prop_temp; + for (std::size_t j=0;j index_container; + index_container.push_back(index_current); + std::vector < int > index_container_former_ring; + index_container_former_ring.push_back(index_current); + std::list < int > index_container_current_ring; + + //propagation + bool propagation=true; + do + { + propagation=false; + + //neighbors + for (std::size_t k=0;k::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); + } + std::cerr << subgraph_clusters.size () << " subgraph clusters" << std::endl; + for (std::size_t i = 0; i < subgraph_clusters.size (); ++ i) + std::cerr << subgraph_clusters[i].size () << " "; + std::cerr << std::endl; + std::cerr << subgraph_clusters_max_area_index.size () << " subgraph clusters max area index" << std::endl + << subgraph_clusters_prop.size () << " subgraph clusters prop" << std::endl; + for (std::size_t i = 0; i < subgraph_clusters_prop.size (); ++ i) + std::cerr << subgraph_clusters_prop[i].size () << " "; + std::cerr << std::endl; + } FT distance_Point (const Point& a, const Point& b) From 7559892be39d4abee3e53b814de8e863c38d0843 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 4 Nov 2015 12:04:21 +0100 Subject: [PATCH 16/48] Remove useless vectors and apply directly transformations to planes --- .../include/CGAL/Plane_regularization.h | 187 ++++-------------- 1 file changed, 35 insertions(+), 152 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 868a515c1da..4b30d918518 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -64,6 +64,7 @@ namespace CGAL { { bool is_free; std::vector planes; + std::vector coplanar_group; std::vector orthogonal_clusters; Vector normal; FT cosangle_vertical; @@ -180,15 +181,12 @@ namespace CGAL { //detecting co-planarity and store in list_coplanar_prim - std::vector > > list_coplanar_prim; for (std::size_t i = 0; i < clusters.size(); ++ i) { - - std::vector< std::vector < int > > list_coplanar_prim_tmp; Vector vec_reg = clusters[i].normal; - std::vector < int > list_cop_index; + for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip) - list_cop_index.push_back(-1); + clusters[i].coplanar_group.push_back (-1); int cop_index=0; @@ -196,161 +194,61 @@ namespace CGAL { { int index_prim = clusters[i].planes[j]; - if (list_cop_index[j] < 0) + if (clusters[i].coplanar_group[j] == static_cast(-1)) { - - std::vector < int > list_coplanar_prim_tmp_tmp; - list_cop_index[j]=cop_index; - list_coplanar_prim_tmp_tmp.push_back(index_prim); + clusters[i].coplanar_group[j] = cop_index; Point pt_reg = m_planes[index_prim]->projection(m_centroids[index_prim]); Plane plan_reg(pt_reg,vec_reg); - for (std::size_t k = j+1; k < clusters[i].planes.size(); ++ k) + for (std::size_t k = j + 1; k < clusters[i].planes.size(); ++ k) { - if (list_cop_index[k] < 0) + if (clusters[i].coplanar_group[k] == static_cast(-1)) { int index_prim_next = clusters[i].planes[k]; Point pt_reg_next = m_planes[index_prim_next]->projection(m_centroids[index_prim_next]); Point pt_proj=plan_reg.projection(pt_reg_next); double distance=distance_Point(pt_reg_next,pt_proj); - if(distance > list_primitive_reg_index_extracted_planes; - std::vector < Plane > list_primitive_reg; - - for (std::size_t i=0;i pt_bary (cop_index); + std::vector 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 = m_planes[index_prim]->projection(m_centroids[index_prim]); - Point pt_bary(0.,0.,0.); - double area=0; - - for (std::size_t k=0; kprojection(m_centroids[index_prim]); - - pt_bary=barycenter(pt_bary, area,pt_reg,m_areas[index_prim]); - area+=m_areas[index_prim]; - } - Vector vec_reg = m_planes[list_coplanar_prim[i][j][0]]->plane_normal (); - - Plane plane_reg(pt_bary,vec_reg); - - bool is_reg_used=false; - std::vector< int > list_primitive_reg_index_extracted_planes_tmp1; - - for (std::size_t k=0; kplane_normal () * plane_reg.orthogonal_vector()) > 1. - epsilon) - { - if(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector()<0) - m_planes[index_prim]->update (plane_reg.opposite()); - else - m_planes[index_prim]->update (plane_reg); - is_reg_used=true; - list_primitive_reg_index_extracted_planes_tmp1.push_back(index_prim); - } - else{ - list_primitive_reg.push_back(static_cast (*(m_planes[index_prim]))); - std::vector< int > list_primitive_reg_index_extracted_planes_tmp; - list_primitive_reg_index_extracted_planes_tmp.push_back(index_prim); - list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp); - } - } - if(is_reg_used) { - list_primitive_reg.push_back(plane_reg); - list_primitive_reg_index_extracted_planes.push_back(list_primitive_reg_index_extracted_planes_tmp1); - } + pt_bary[group] = barycenter (pt_bary[group], area[group], pt_reg, m_areas[index_prim]); + area[group] += m_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 (m_planes[index_prim]->plane_normal () + * plane_reg.orthogonal_vector() < 0) + m_planes[index_prim]->update (plane_reg.opposite()); + else + m_planes[index_prim]->update (plane_reg); } } - - // std::cout< > plane_point_index_temp; - // std::vector < Plane > extracted_planes_temp; - // std::vector < bool > has_been_merged; - // for (std::size_t i=0; i< m_planes.size();i++) - // has_been_merged.push_back(false); - - // for (std::size_t i=0; i< m_planes.size();i++) - // { - - // if (!has_been_merged[i]) - // { - // extracted_planes_temp.push_back (m_planes[i]); - // int label_index=extracted_planes_temp.size()-1; - // plane_point_index_temp.push_back(m_planes[i]->indices_of_assigned_points[i]); - // for (std::size_t k=0; k< m_planes[i]->indices_of_assigned_points().size();k++) - // { - // int index_pt=m_planes[i]->indices_of_assigned_points[k]; - // primitive_index[index_pt]=label_index; - // label_plane[index_pt]=label_index; - // } - - // for (std::size_t j=i+1;j< m_planes.size();j++) - // { - - // if(m_planes[i]==m_planes[j]) - // { //if identical (do opposite plane too ?) then store the second in the first - - // has_been_merged[j]=true; - - // std::vector< int > plane_point_index_new - // = plane_point_index_temp[m_planes.size()-1]; - // for (std::size_t k=0; k< m_planes[j]->indices_of_assigned_points().size();k++) - // { - // int ind=m_planes[j]->indices_of_assigned_points[k]; - // plane_point_index_new.push_back(ind); - // primitive_index[ind]=label_index; - // label_plane[ind]=label_index; - // } - // plane_point_index_temp[plane_point_index_temp.size()-1]=plane_point_index_new; - - // } - // } - // } - // } - - // TODO - // m_planes=extracted_planes_temp; - // m_planes->indices_of_assigned_points=plane_point_index_temp; - - return list_primitive_reg.size (); + return clusters.size (); } void compute_centroids_and_areas () @@ -501,12 +399,6 @@ namespace CGAL { } } - for (std::size_t i = 0; i < cosangle_centroids.size (); ++ i) - std::cerr << cosangle_centroids[i] << ", "; - std::cerr < Date: Wed, 4 Nov 2015 14:11:59 +0100 Subject: [PATCH 17/48] Use same color for parallel planes --- ..._demo_point_set_shape_detection_plugin.cpp | 35 +++++++++++++++---- 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp index b67e1d5d52f..edf11eb4bb5 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp @@ -94,6 +94,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() { @@ -163,6 +164,18 @@ 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 (op.epsilon, 0.1); + + std::cerr << "done" << std::endl; + } + + std::map color_map; + //print_message(QString("%1 shapes found.").arg(shape_detection.number_of_shapes())); int index = 0; BOOST_FOREACH(boost::shared_ptr shape, shape_detection.shapes()) @@ -182,9 +195,11 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered point_item->point_set()->unselect_all (); unsigned char r, g, b; + r = static_cast(64 + rand.get_int(0, 192)); g = static_cast(64 + rand.get_int(0, 192)); b = static_cast(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 @@ -198,16 +213,29 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered { ss << item->name().toStdString() << "_plane_"; + + if (dialog.regularize ()) + { + Kernel::Point_3 ref + = CGAL::ORIGIN + (dynamic_cast*>(shape.get ()))->plane_normal (); + + if (color_map.find (ref) == color_map.end ()) + color_map[ref] = point_item->color (); + else + point_item->setColor (color_map[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*>(shape.get ()))); + std::cerr << plane << std::endl; build_alpha_shape (*(point_item->point_set()), plane, 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); @@ -235,11 +263,6 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered ++index; } - std::cerr << "Regularization of planes... " << std::endl; - Regularization regularization (*points, shape_detection); - regularization.run (op.epsilon, op.normal_threshold); - - std::cerr << "done" << std::endl; // Updates scene scene->itemChanged(index); From 58f4a7b74f32763d7357857ae01f79be3358ed7d Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 4 Nov 2015 15:08:59 +0100 Subject: [PATCH 18/48] Correction of variable name and meaning --- .../include/CGAL/Plane_regularization.h | 41 ++++++++++--------- ..._demo_point_set_shape_detection_plugin.cpp | 2 +- 2 files changed, 23 insertions(+), 20 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 4b30d918518..87bf0cc9327 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -106,9 +106,7 @@ namespace CGAL { // Ignore all shapes other than plane if (pshape == boost::shared_ptr()) continue; - m_planes.push_back (pshape); - } } @@ -127,14 +125,14 @@ namespace CGAL { - std::size_t run (FT epsilon, FT tolerance_coplanarity) + std::size_t run (FT tolerance_cosangle, FT tolerance_coplanarity) { compute_centroids_and_areas (); // clustering the parallel primitives and store them in clusters // & compute the normal, size and cos angle to the vertical of each cluster std::vector clusters; - compute_parallel_clusters (clusters, epsilon); + compute_parallel_clusters (clusters, tolerance_cosangle); //discovery orthogonal relationship between clusters for (std::size_t i = 0; i < clusters.size(); ++ i) @@ -142,7 +140,7 @@ namespace CGAL { for (std::size_t j = i + 1; j < clusters.size(); ++ j) { - if (std::fabs (clusters[i].normal * clusters[j].normal) < epsilon) + 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); @@ -153,7 +151,7 @@ namespace CGAL { //clustering the vertical cosangle and store their centroids in //cosangle_centroids and the centroid index of each cluster in //list_cluster_index - cluster_vertical_cosangles (clusters, epsilon); + cluster_vertical_cosangles (clusters, tolerance_cosangle); //find subgraphs of mutually orthogonal clusters (store index of //clusters in subgraph_clusters), and select the cluster of @@ -174,7 +172,7 @@ namespace CGAL { vec_reg=-vec_reg; Plane plane_reg(pt_reg,vec_reg); - if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - epsilon) + if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - tolerance_cosangle) m_planes[index_prim]->update (plane_reg); } } @@ -219,7 +217,7 @@ namespace CGAL { } //regularize primitive position by computing barycenter of cplanar planes - std::vector pt_bary (cop_index); + std::vector pt_bary (cop_index, Point (0., 0., 0.)); std::vector area (cop_index, 0.); for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) @@ -229,10 +227,14 @@ namespace CGAL { Point pt_reg = m_planes[index_prim]->projection(m_centroids[index_prim]); - pt_bary[group] = barycenter (pt_bary[group], area[group], pt_reg, m_areas[index_prim]); + pt_bary[group] = CGAL::barycenter (pt_bary[group], area[group], pt_reg, m_areas[index_prim]); area[group] += m_areas[index_prim]; } + + for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) + std::cerr << pt_bary[clusters[i].coplanar_group[j]] << " + " << vec_reg << std::endl; + for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) { std::size_t index_prim = clusters[i].planes[j]; @@ -267,7 +269,7 @@ namespace CGAL { } } - void compute_parallel_clusters (std::vector& clusters, double epsilon) + void compute_parallel_clusters (std::vector& clusters, double tolerance_cosangle) { // find pairs of epsilon-parallel primitives and store them in parallel_planes std::vector > parallel_planes (m_planes.size ()); @@ -282,7 +284,7 @@ namespace CGAL { Vector v2 = m_planes[i]->plane_normal (); - if (std::fabs (v1 * v2) > 1. - epsilon) + if (std::fabs (v1 * v2) > 1. - tolerance_cosangle) parallel_planes[i].push_back (j); } } @@ -329,7 +331,7 @@ namespace CGAL { Vector normal_it = m_planes[it]->plane_normal (); if(is_available[it] - && std::fabs (normal_it*clu.normal) > 1. - epsilon ) + && std::fabs (normal_it*clu.normal) > 1. - tolerance_cosangle ) { propagation = true; index_container_current_ring_parallel.push_back(it); @@ -367,7 +369,7 @@ namespace CGAL { is_available.clear(); } - void cluster_vertical_cosangles (std::vector& clusters, double epsilon) + void cluster_vertical_cosangles (std::vector& clusters, double tolerance_cosangle) { std::vector < double > cosangle_centroids; std::vector < int > list_cluster_index; @@ -386,7 +388,7 @@ namespace CGAL { for (std::size_t j = i+1; j < clusters.size(); ++ j) { if (list_cluster_index[j] < 0 && std::fabs (clusters[j].cosangle_vertical - - mean / mean_area) < epsilon) + mean / mean_area) < tolerance_cosangle) { list_cluster_index[j] = mean_index; mean_area += clusters[j].area; @@ -403,9 +405,9 @@ namespace CGAL { //desactive Z-verticalité for (std::size_t i = 0; i < cosangle_centroids.size(); ++ i) { - if (cosangle_centroids[i] < epsilon) + if (cosangle_centroids[i] < tolerance_cosangle) cosangle_centroids[i] = 0; - else if (cosangle_centroids[i] > 1. - epsilon) + else if (cosangle_centroids[i] > 1. - tolerance_cosangle) cosangle_centroids[i] = 1; } for (std::size_t i = 0; i < clusters.size(); ++ i) @@ -574,9 +576,10 @@ namespace CGAL { Vector regularize_normal (const Vector& n, FT cos_vertical) { + std::cerr << "N = " << n << std::endl; FT A = 1 - cos_vertical * cos_vertical; FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); - + FT vx = std::sqrt (A/B); if (n.x() < 0) @@ -585,7 +588,7 @@ namespace CGAL { FT vy = vx * (n.y() / n.x()); Vector res (vx, vy, cos_vertical); - + std::cerr << "Res = " << res << std::endl; return res / std::sqrt (res * res); } @@ -629,7 +632,7 @@ namespace CGAL { Vector res (vx, vy, cos_vertical); FT norm = std::max(1e-5, 1. / sqrt(res.squared_length ())); - + std::cerr << "NRes = " << norm * res << std::endl; return norm * res; } diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp index edf11eb4bb5..151ebd917e0 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp @@ -169,7 +169,7 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered { std::cerr << "Regularization of planes... " << std::endl; Regularization regularization (*points, shape_detection); - regularization.run (op.epsilon, 0.1); + regularization.run (1. - op.normal_threshold, op.epsilon); std::cerr << "done" << std::endl; } From 2f2910dd502b4c6f57fb91d5d36f9cd4ba072fad Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Thu, 5 Nov 2015 11:09:46 +0100 Subject: [PATCH 19/48] Bugfix (0/0 makes NaN) + adding parameters + start writing doc --- .../include/CGAL/Plane_regularization.h | 70 +++++++++++++++---- 1 file changed, 56 insertions(+), 14 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 87bf0cc9327..6e4855290a2 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -37,6 +37,35 @@ namespace CGAL { + +/*! +\ingroup PkgPointSetShapeDetection3 +\brief A plane regularization algorithm applied as a post-processing +to a shape detection algorithm. + +Given a set of detected planes with their respective inlier sets, this +class enables to regularize the planes: planes almost parallel are +made exactly parallel. In addition, some additional regularization can +be performed: + +- Plane clusters that are almost orthogonal can be made exactly + orthogonal. + +- Planes that are parallel and almost coplanar can be made exactly + coplanar. + +- Planes that are almost 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{verdie2015lod}. + +\tparam Traits a model of `EfficientRANSACTraits` + +*/ template class Plane_regularization { @@ -125,7 +154,10 @@ namespace CGAL { - std::size_t run (FT tolerance_cosangle, FT tolerance_coplanarity) + std::size_t run (FT tolerance_cosangle, FT tolerance_coplanarity, + bool regularize_coplanarity = true, + bool regularize_orthogonality = true, + bool z_symmetry = true) { compute_centroids_and_areas (); @@ -576,20 +608,30 @@ namespace CGAL { Vector regularize_normal (const Vector& n, FT cos_vertical) { - std::cerr << "N = " << n << std::endl; - FT A = 1 - cos_vertical * cos_vertical; - FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); - - FT vx = std::sqrt (A/B); - - if (n.x() < 0) - vx = -vx; - - FT vy = vx * (n.y() / n.x()); + FT vz = cos_vertical; + FT A = 1 - cos_vertical*cos_vertical; + + if (n.x() != 0) + { + FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); + FT vx = std::sqrt(A / B); + if (n.x() < 0) + vx = -vx; + FT vy = vx * (n.y() / n.x()); + + Vector res (vx, vy, vz); + return res / std::sqrt (res * res);; + } + else + { + FT vx = 0; + FT vy = sqrt(A); + if (n.y() < 0) + vy = -vy; + Vector res(vx, vy, vz); + return res / std::sqrt (res * res);; + } - Vector res (vx, vy, cos_vertical); - std::cerr << "Res = " << res << std::endl; - return res / std::sqrt (res * res); } From 1d44ead4bce773734a22c4ced38ad9ddf03e648d Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Thu, 5 Nov 2015 13:52:46 +0100 Subject: [PATCH 20/48] Documentation of plane regularization --- Documentation/biblio/cgal_manual.bib | 9 ++ .../PackageDescription.txt | 3 + .../include/CGAL/Plane_regularization.h | 86 +++++++++++++++++-- 3 files changed, 89 insertions(+), 9 deletions(-) diff --git a/Documentation/biblio/cgal_manual.bib b/Documentation/biblio/cgal_manual.bib index 2590e3c8f52..e5c6c3c687d 100644 --- a/Documentation/biblio/cgal_manual.bib +++ b/Documentation/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" diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt index a4ec80b21f0..7dd6c128def 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt @@ -47,5 +47,8 @@ - `CGAL::Shape_detection_3::Cylinder` - `CGAL::Shape_detection_3::Cone` - `CGAL::Shape_detection_3::Torus` + +## Regularization Classes ## +- `CGAL::Plane_regularization` */ diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 6e4855290a2..9f2562df74b 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -61,7 +61,7 @@ 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{verdie2015lod}. +The implementation follows \cgalCite{cgal:vla-lod-15}. \tparam Traits a model of `EfficientRANSACTraits` @@ -71,21 +71,33 @@ The implementation follows \cgalCite{verdie2015lod}. { public: + /// \cond SKIP_IN_MANUAL typedef Plane_regularization Self; + /// \endcond + /// \name Types + /// @{ + /// \cond SKIP_IN_MANUAL typedef typename Traits::FT FT; typedef typename Traits::Point_3 Point; typedef typename Traits::Vector_3 Vector; - typedef typename Traits::Plane_3 Plane; + + /// \endcond + typedef typename Traits::Plane_3 Plane; ///< Raw plane type typedef typename Traits::Point_map Point_map; + ///< property map to access the location of an input point. typedef typename Traits::Normal_map Normal_map; + ///< property map to access the unoriented normal of an input point typedef typename Traits::Input_range Input_range; - typedef typename Input_range::iterator Input_iterator; + ///< Model of the concept `Range` with random access iterators, providing input points and normals + /// through the following two property maps. - typedef Shape_detection_3::Shape_base Shape; - typedef Shape_detection_3::Plane Plane_shape; - + typedef typename Input_range::iterator Input_iterator; ///< Iterator on input data + + typedef Shape_detection_3::Shape_base Shape; ///< Shape type. + typedef Shape_detection_3::Plane Plane_shape; ///< Plane type. + /// @} private: @@ -114,12 +126,32 @@ The implementation follows \cgalCite{verdie2015lod}. public: + /// \name Initialization + /// @{ + /*! + Constructs an empty plane regularization engine. + */ Plane_regularization (Traits t = Traits ()) : m_traits (t) { } + /*! + + Constructs a plane regularization engine base on an input range + of points with its related shape detection engine. + + \param input_range Range of input data. + + \param shape_detection Shape detection engine 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 and must have been using `input_range` as input. + + */ Plane_regularization (Input_range& input_range, const Shape_detection_3::Efficient_RANSAC& shape_detection) : m_traits (shape_detection.traits()) @@ -140,22 +172,54 @@ The implementation follows \cgalCite{verdie2015lod}. } + /*! + Releases all memory allocated by this instance. + */ virtual ~Plane_regularization () { clear (); } + /// @} + /// \name Memory Management + /// @{ + /*! + Clear all internal structures. + */ void clear () { std::vector > ().swap (m_planes); std::vector ().swap (m_centroids); std::vector ().swap (m_areas); } + /// @} + /// \name Regularization + /// @{ + /*! - - std::size_t run (FT tolerance_cosangle, FT tolerance_coplanarity, - bool regularize_coplanarity = true, + Performs the plane regularization. Planes are directly modified. + + \param tolerance_cosangle Tolerance of deviation between normal + vectors of planes so that they are considered parallel, + expressed as the cosinus of their angle. + + \param tolerance_coplanarity Maximal distance between two + parallel planes such that they are considered coplanar. The + default value is 0, meaning that coplanarity is not taken into + account for regularization. + + \param regularize_orthogonality Make almost orthogonal clusters + of plane exactly orthogonal. + + \param z_symmetry Make almost Z-symmetrical clusters exactly + Z-symmetrical. + + \return The number of clusters of parallel planes found. + */ + + std::size_t run (FT tolerance_cosangle = 0.1, + FT tolerance_coplanarity = 0.0, bool regularize_orthogonality = true, bool z_symmetry = true) { @@ -284,7 +348,11 @@ The implementation follows \cgalCite{verdie2015lod}. return clusters.size (); } + /// @} + + private: + void compute_centroids_and_areas () { for (std::size_t i = 0; i < m_planes.size (); ++ i) From 06fb8a79164500c2c74ccb8444faaf975d5fafb1 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 6 Nov 2015 14:16:53 +0100 Subject: [PATCH 21/48] Generalize symmetry from Z-axis to any user-given vector --- .../include/CGAL/Plane_regularization.h | 126 +++++++++--------- 1 file changed, 65 insertions(+), 61 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 9f2562df74b..4b770bcbdc8 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -81,6 +81,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. typedef typename Traits::FT FT; typedef typename Traits::Point_3 Point; typedef typename Traits::Vector_3 Vector; + typedef typename Traits::Line_3 Line; /// \endcond typedef typename Traits::Plane_3 Plane; ///< Raw plane type @@ -221,14 +222,14 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. std::size_t run (FT tolerance_cosangle = 0.1, FT tolerance_coplanarity = 0.0, bool regularize_orthogonality = true, - bool z_symmetry = true) + Vector symmetry_direction = Vector (0., 0., 1.)) { compute_centroids_and_areas (); // clustering the parallel primitives and store them in clusters // & compute the normal, size and cos angle to the vertical of each cluster std::vector clusters; - compute_parallel_clusters (clusters, tolerance_cosangle); + compute_parallel_clusters (clusters, tolerance_cosangle, symmetry_direction); //discovery orthogonal relationship between clusters for (std::size_t i = 0; i < clusters.size(); ++ i) @@ -252,7 +253,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. //find subgraphs of mutually orthogonal clusters (store index of //clusters in subgraph_clusters), and select the cluster of //largest area - subgraph_mutually_orthogonal_clusters (clusters); + subgraph_mutually_orthogonal_clusters (clusters, symmetry_direction); //recompute optimal plane for each primitive after normal regularization for (std::size_t i=0; i < clusters.size(); ++ i) @@ -369,7 +370,8 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } } - void compute_parallel_clusters (std::vector& clusters, double tolerance_cosangle) + void compute_parallel_clusters (std::vector& clusters, double tolerance_cosangle, + const Vector& symmetry_direction) { // find pairs of epsilon-parallel primitives and store them in parallel_planes std::vector > parallel_planes (m_planes.size ()); @@ -462,8 +464,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } while(propagation); - Vector v_vertical(0.,0.,1.); - clu.cosangle_vertical = std::fabs(v_vertical*clu.normal); + clu.cosangle_vertical = std::fabs(symmetry_direction * clu.normal); } } is_available.clear(); @@ -514,7 +515,8 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. clusters[i].cosangle_vertical = cosangle_centroids[list_cluster_index[i]]; } - void subgraph_mutually_orthogonal_clusters (std::vector& clusters) + void subgraph_mutually_orthogonal_clusters (std::vector& clusters, + const Vector& symmetry_direction) { std::vector < std::vector < int > > subgraph_clusters; std::vector < int > subgraph_clusters_max_area_index; @@ -614,6 +616,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. int index_current=subgraph_clusters_max_area_index[i]; Vector vec_current=regularize_normal(clusters[index_current].normal, + symmetry_direction, clusters[index_current].cosangle_vertical); clusters[index_current].normal = vec_current; clusters[index_current].is_free = false; @@ -649,6 +652,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. Vector new_vect=regularize_normals_from_prior(clusters[cluster_index].normal, clusters[j].normal, + symmetry_direction, clusters[j].cosangle_vertical); clusters[j].normal = new_vect; } @@ -674,76 +678,76 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. return std::sqrt (CGAL::squared_distance (a, b)); } - Vector regularize_normal (const Vector& n, FT cos_vertical) + Vector regularize_normal (const Vector& n, const Vector& symmetry_direction, + FT cos_vertical) { - FT vz = cos_vertical; - FT A = 1 - cos_vertical*cos_vertical; - - if (n.x() != 0) - { - FT B = 1 + (n.y() * n.y()) / (n.x() * n.x()); - FT vx = std::sqrt(A / B); - if (n.x() < 0) - vx = -vx; - FT vy = vx * (n.y() / n.x()); + Point pt_symmetry = CGAL::ORIGIN + cos_vertical* symmetry_direction; - Vector res (vx, vy, vz); - return res / std::sqrt (res * res);; - } + 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; + + double delta = std::sqrt (1 - cos_vertical * cos_vertical); + + 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 - { - FT vx = 0; - FT vy = sqrt(A); - if (n.y() < 0) - vy = -vy; - Vector res(vx, vy, vz); - return res / std::sqrt (res * res);; - } - + return n; } Vector regularize_normals_from_prior (const Vector& np, const Vector& n, + const Vector& symmetry_direction, FT cos_vertical) { - FT vx, vy; + Plane plane_orthogonality (CGAL::ORIGIN, np); + Point pt_symmetry = CGAL::ORIGIN + cos_vertical* symmetry_direction; - if (np.x() != 0) - { - FT a = (np.y() * np.y()) / (np.x() * np.x()) + 1; - FT b = 2 * np.y() * np.z() * cos_vertical / np.x(); - FT c= cos_vertical * cos_vertical-1; + 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 (n, symmetry_direction, cos_vertical); - if (4 * a * c > b * b) - return regularize_normal (n, cos_vertical); - else - { - FT delta = std::sqrt (b * b-4 * a * c); - FT vy1= (-b-delta) / (2 * a); - FT vy2= (-b+delta) / (2 * a); + Point projected_origin = line.projection (CGAL::ORIGIN); + FT R = CGAL::squared_distance (Point (CGAL::ORIGIN), projected_origin); - vy = (std::fabs(n.y()-vy1) < std::fabs(n.y()-vy2)) - ? vy1 : vy2; - - vx = (-np.y() * vy-np.z() * cos_vertical) / np.x(); - } - } - else if (np.y() != 0) + if (R <= 1) // 2 (or 1) possible points intersecting the unit sphere and line { - vy = -np.z() * cos_vertical / np.y(); - vx = std::sqrt (1 - cos_vertical * cos_vertical - vy * vy); - - if (n.x() < 0) - vx = -vx; + double delta = std::sqrt (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 - return regularize_normal (n, cos_vertical); + else //no point intersecting the unit sphere and line + return regularize_normal (n,symmetry_direction, cos_vertical); - Vector res (vx, vy, cos_vertical); - FT norm = std::max(1e-5, 1. / sqrt(res.squared_length ())); - std::cerr << "NRes = " << norm * res << std::endl; - return norm * res; } From 82ee9c8ce07af5c83457cb19231b1f6780e1a110 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 6 Nov 2015 14:34:23 +0100 Subject: [PATCH 22/48] Default behavior: no symmetry required --- .../include/CGAL/Plane_regularization.h | 75 +++++++++++-------- 1 file changed, 43 insertions(+), 32 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 4b770bcbdc8..1f8b0a2bc95 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -109,7 +109,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. std::vector coplanar_group; std::vector orthogonal_clusters; Vector normal; - FT cosangle_vertical; + FT cosangle_symmetry; FT area; FT cosangle_centroid; }; @@ -222,38 +222,44 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. std::size_t run (FT tolerance_cosangle = 0.1, FT tolerance_coplanarity = 0.0, bool regularize_orthogonality = true, - Vector symmetry_direction = Vector (0., 0., 1.)) + Vector symmetry_direction = CGAL::NULL_VECTOR) { compute_centroids_and_areas (); // clustering the parallel primitives and store them in clusters - // & compute the normal, size and cos angle to the vertical of each cluster + // & compute the normal, size and cos angle to the symmetry + // direction of each cluster std::vector clusters; compute_parallel_clusters (clusters, tolerance_cosangle, symmetry_direction); - //discovery orthogonal relationship between clusters - for (std::size_t i = 0; i < clusters.size(); ++ i) + if (regularize_orthogonality) { - for (std::size_t j = i + 1; j < clusters.size(); ++ j) + //discovery orthogonal relationship between clusters + for (std::size_t i = 0; i < clusters.size(); ++ i) { - - if (std::fabs (clusters[i].normal * clusters[j].normal) < tolerance_cosangle) + for (std::size_t j = i + 1; j < clusters.size(); ++ j) { - clusters[i].orthogonal_clusters.push_back (j); - clusters[j].orthogonal_clusters.push_back (i); + + 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); + } } } } - - //clustering the vertical cosangle and store their centroids in + + //clustering the symmetry cosangle and store their centroids in //cosangle_centroids and the centroid index of each cluster in //list_cluster_index - cluster_vertical_cosangles (clusters, tolerance_cosangle); + if (symmetry_direction != CGAL::NULL_VECTOR) + cluster_symmetric_cosangles (clusters, tolerance_cosangle); //find subgraphs of mutually orthogonal clusters (store index of //clusters in subgraph_clusters), and select the cluster of //largest area - subgraph_mutually_orthogonal_clusters (clusters, symmetry_direction); + if (regularize_orthogonality) + subgraph_mutually_orthogonal_clusters (clusters, symmetry_direction); //recompute optimal plane for each primitive after normal regularization for (std::size_t i=0; i < clusters.size(); ++ i) @@ -464,13 +470,14 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } while(propagation); - clu.cosangle_vertical = std::fabs(symmetry_direction * clu.normal); + if (symmetry_direction != CGAL::NULL_VECTOR) + clu.cosangle_symmetry = std::fabs(symmetry_direction * clu.normal); } } is_available.clear(); } - void cluster_vertical_cosangles (std::vector& clusters, double tolerance_cosangle) + void cluster_symmetric_cosangles (std::vector& clusters, double tolerance_cosangle) { std::vector < double > cosangle_centroids; std::vector < int > list_cluster_index; @@ -483,17 +490,17 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if(list_cluster_index[i]<0) { list_cluster_index[i] = mean_index; - double mean = clusters[i].area * clusters[i].cosangle_vertical; + double mean = clusters[i].area * clusters[i].cosangle_symmetry; double mean_area = clusters[i].area; for (std::size_t j = i+1; j < clusters.size(); ++ j) { - if (list_cluster_index[j] < 0 && std::fabs (clusters[j].cosangle_vertical - + if (list_cluster_index[j] < 0 && 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_vertical; + mean += clusters[j].area * clusters[j].cosangle_symmetry; } } ++ mean_index; @@ -502,8 +509,6 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } } - - //desactive Z-verticalité for (std::size_t i = 0; i < cosangle_centroids.size(); ++ i) { if (cosangle_centroids[i] < tolerance_cosangle) @@ -512,7 +517,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. cosangle_centroids[i] = 1; } for (std::size_t i = 0; i < clusters.size(); ++ i) - clusters[i].cosangle_vertical = cosangle_centroids[list_cluster_index[i]]; + clusters[i].cosangle_symmetry = cosangle_centroids[list_cluster_index[i]]; } void subgraph_mutually_orthogonal_clusters (std::vector& clusters, @@ -606,7 +611,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. //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 vertical + //orthogonality and cosangle to symmetry direction for (std::size_t i = 0; i < clusters.size(); ++ i) clusters[i].is_free = true; @@ -617,7 +622,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. int index_current=subgraph_clusters_max_area_index[i]; Vector vec_current=regularize_normal(clusters[index_current].normal, symmetry_direction, - clusters[index_current].cosangle_vertical); + clusters[index_current].cosangle_symmetry); clusters[index_current].normal = vec_current; clusters[index_current].is_free = false; @@ -653,7 +658,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. Vector new_vect=regularize_normals_from_prior(clusters[cluster_index].normal, clusters[j].normal, symmetry_direction, - clusters[j].cosangle_vertical); + clusters[j].cosangle_symmetry); clusters[j].normal = new_vect; } } @@ -679,9 +684,12 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } Vector regularize_normal (const Vector& n, const Vector& symmetry_direction, - FT cos_vertical) + FT cos_symmetry) { - Point pt_symmetry = CGAL::ORIGIN + cos_vertical* symmetry_direction; + 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; @@ -694,7 +702,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (!assign(line, ob_1)) return n; - double delta = std::sqrt (1 - cos_vertical * cos_vertical); + double delta = std::sqrt (1 - cos_symmetry * cos_symmetry); Point projected_origin = line.projection (CGAL::ORIGIN); Vector line_vector (line); @@ -716,17 +724,20 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. Vector regularize_normals_from_prior (const Vector& np, const Vector& n, const Vector& symmetry_direction, - FT cos_vertical) + FT cos_symmetry) { + if (symmetry_direction == CGAL::NULL_VECTOR) + return n; + Plane plane_orthogonality (CGAL::ORIGIN, np); - Point pt_symmetry = CGAL::ORIGIN + cos_vertical* symmetry_direction; + 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 (n, symmetry_direction, cos_vertical); + return regularize_normal (n, symmetry_direction, cos_symmetry); Point projected_origin = line.projection (CGAL::ORIGIN); FT R = CGAL::squared_distance (Point (CGAL::ORIGIN), projected_origin); @@ -746,7 +757,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. return Vector (CGAL::ORIGIN, pt2); } else //no point intersecting the unit sphere and line - return regularize_normal (n,symmetry_direction, cos_vertical); + return regularize_normal (n,symmetry_direction, cos_symmetry); } From 939e541ba5b7152ae07f8a9e8e8cda1442bed95b Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 6 Nov 2015 18:02:19 +0100 Subject: [PATCH 23/48] More on manual and improve clarity of parameters --- .../Point_set_shape_detection_3.txt | 8 +++++++ .../include/CGAL/Plane_regularization.h | 21 ++++++++++--------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt index 35dc2696a89..27271847c4c 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt @@ -105,6 +105,14 @@ 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. In such input data, shapes often come with specific relationships between them: parallelism, coplanarity or orthogonality, for example. \cgal provides a tool to regularize planes detected on a point set. + +The class `CGAL::Plane_regularization` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). + + + \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}. diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 1f8b0a2bc95..2c7fe70bea9 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -201,9 +201,9 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. Performs the plane regularization. Planes are directly modified. - \param tolerance_cosangle Tolerance of deviation between normal - vectors of planes so that they are considered parallel, - expressed as the cosinus of their angle. + \param tolerance_angle Tolerance of deviation between normal + vectors of planes so that they are considered parallel (in + degrees). \param tolerance_coplanarity Maximal distance between two parallel planes such that they are considered coplanar. The @@ -213,18 +213,22 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. \param regularize_orthogonality Make almost orthogonal clusters of plane exactly orthogonal. - \param z_symmetry Make almost Z-symmetrical clusters exactly - Z-symmetrical. + \param symmetry_direction Make clusters that are almost + symmetrical in the symmetry direction exactly symmetrical. This + parameter is ignored if it is equal to `CGAL::NULL_VECTOR` + (default value). \return The number of clusters of parallel planes found. */ - std::size_t run (FT tolerance_cosangle = 0.1, + std::size_t run (FT tolerance_angle = 25.0, FT tolerance_coplanarity = 0.0, bool regularize_orthogonality = true, Vector symmetry_direction = CGAL::NULL_VECTOR) { compute_centroids_and_areas (); + + FT tolerance_cosangle = 1. - std::cos (tolerance_angle); // clustering the parallel primitives and store them in clusters // & compute the normal, size and cos angle to the symmetry @@ -274,7 +278,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if( m_planes[index_prim]->plane_normal () * vec_reg < 0) vec_reg=-vec_reg; Plane plane_reg(pt_reg,vec_reg); - + if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - tolerance_cosangle) m_planes[index_prim]->update (plane_reg); } @@ -334,9 +338,6 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. area[group] += m_areas[index_prim]; } - for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) - std::cerr << pt_bary[clusters[i].coplanar_group[j]] << " + " << vec_reg << std::endl; - for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) { From f93a21c0822f18f1236181cc0cadccbd45810cf4 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 9 Nov 2015 07:52:02 +0100 Subject: [PATCH 24/48] Better naming of RANSAC primitives in plugin --- ...n_demo_point_set_shape_detection_plugin.cpp | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp index 151ebd917e0..e64f228a872 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp @@ -169,7 +169,8 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered { std::cerr << "Regularization of planes... " << std::endl; Regularization regularization (*points, shape_detection); - regularization.run (1. - op.normal_threshold, op.epsilon); + + regularization.run (180 * std::acos (op.normal_threshold) / M_PI, op.epsilon); std::cerr << "done" << std::endl; } @@ -214,15 +215,24 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered ss << item->name().toStdString() << "_plane_"; - if (dialog.regularize ()) + // if (dialog.regularize ()) { Kernel::Point_3 ref = CGAL::ORIGIN + (dynamic_cast*>(shape.get ()))->plane_normal (); if (color_map.find (ref) == color_map.end ()) - color_map[ref] = point_item->color (); + { + ref + = CGAL::ORIGIN + (-1.) * (dynamic_cast*>(shape.get ()))->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 ()) @@ -231,7 +241,7 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered Scene_polyhedron_item* poly_item = new Scene_polyhedron_item; Plane_3 plane = (Plane_3)(*(dynamic_cast*>(shape.get ()))); - std::cerr << plane << std::endl; + build_alpha_shape (*(point_item->point_set()), plane, poly_item, dialog.cluster_epsilon()); From 192d8c00b4bba7dd4d21c1bb26c1453e6c2ca074 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 9 Nov 2015 12:26:40 +0100 Subject: [PATCH 25/48] Additional functions in plane shape to avoid constructing Plane_3 object --- .../include/CGAL/Shape_detection_3/Plane.h | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h index 10f13fd75f6..294c5ce4439 100644 --- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Plane.h +++ b/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 @@ -88,16 +89,27 @@ namespace CGAL { /*! Computes the orthogonal projection of a query point on the shape. */ - Point_3 projection (const Point_3 &p) const { + 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 m_point_on_primitive + (v * m_base1) * m_base1 + (v * m_base2) * m_base2; + 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 (); m_base2 = p.base2 (); m_normal = p.orthogonal_vector (); + 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)); From 9db8f461b142ed7614d974200e9ef33b656707b7 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 9 Nov 2015 12:31:18 +0100 Subject: [PATCH 26/48] Use functions of plane shape instead of Plane_3 --- ..._demo_point_set_shape_detection_plugin.cpp | 60 ++++++++++--------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp index e64f228a872..3cafd7a7594 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_point_set_shape_detection_plugin.cpp @@ -39,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 PointPMap; + typedef CGAL::Normal_of_point_with_normal_pmap NormalPMap; + + typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits Traits; + typedef CGAL::Shape_detection_3::Efficient_RANSAC Shape_detection; + typedef CGAL::Plane_regularization Regularization; + + public: void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface) { actionDetect = new QAction(tr("Point set shape detection"), mainWindow); @@ -66,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 > plane, Scene_polyhedron_item* item, double epsilon); }; // end Polyhedron_demo_point_set_shape_detection_plugin @@ -125,12 +133,6 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered QApplication::setOverrideCursor(Qt::WaitCursor); - typedef CGAL::Identity_property_map PointPMap; - typedef CGAL::Normal_of_point_with_normal_pmap NormalPMap; - - typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits Traits; - typedef CGAL::Shape_detection_3::Efficient_RANSAC Shape_detection; - typedef CGAL::Plane_regularization Regularization; Shape_detection shape_detection; shape_detection.set_input(*points); @@ -214,35 +216,30 @@ void Polyhedron_demo_point_set_shape_detection_plugin::on_actionDetect_triggered { ss << item->name().toStdString() << "_plane_"; + boost::shared_ptr > pshape + = boost::dynamic_pointer_cast > (shape); - // if (dialog.regularize ()) - { - Kernel::Point_3 ref - = CGAL::ORIGIN + (dynamic_cast*>(shape.get ()))->plane_normal (); + 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 ()) - { - ref - = CGAL::ORIGIN + (-1.) * (dynamic_cast*>(shape.get ()))->plane_normal (); - if (color_map.find (ref) == color_map.end ()) - color_map[ref] = point_item->color (); - else - point_item->setColor (color_map[ref]); - } + color_map[ref] = point_item->color (); else point_item->setColor (color_map[ref]); - - ss << "(" << 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*>(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->setColor(point_item->color ()); @@ -292,7 +289,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 > plane, + Scene_polyhedron_item* item, double epsilon) { typedef Kernel::Point_2 Point_2; typedef CGAL::Alpha_shape_vertex_base_2 Vb; @@ -306,7 +304,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); @@ -328,7 +326,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 ()); } } @@ -339,6 +337,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; } From d46674d0c1ece71a1a6d5032765d65c06d098e89 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 29 Dec 2015 15:27:54 +0100 Subject: [PATCH 27/48] Add regularization to test_scene --- .../test/Point_set_shape_detection_3/test_scene.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index 07fe4c9cceb..5282167816a 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -31,6 +32,8 @@ bool test_scene() { typedef CGAL::Shape_detection_3::Sphere Sphere; typedef CGAL::Shape_detection_3::Torus Torus; + typedef CGAL::Plane_regularization Regularization; + Pwn_vector points; // Loads point set from a file. @@ -126,6 +129,10 @@ bool test_scene() { return false; } + // Test regularization + Regularization regularization (points, ransac); + regularization.run (50., 0.01); + Point_index_range pts = ransac.indices_of_unassigned_points(); std::cout << " succeeded" << std::endl; From bb169a00d0e504ef1d9c1755a2a18bcc02831daf Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 08:31:02 +0100 Subject: [PATCH 28/48] Remove extra semicolon after namespace (pedantic error) --- Point_set_shape_detection_3/include/CGAL/Plane_regularization.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 2c7fe70bea9..b041410cbb1 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -767,6 +767,6 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. }; -}; // namespace CGAL +} // namespace CGAL #endif // CGAL_PLANE_REGULARIZATION_H From ad1feb89e1aecf197fba63c210e77c98930888c4 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 09:16:31 +0100 Subject: [PATCH 29/48] Warning corrections: incorrect types --- .../include/CGAL/Plane_regularization.h | 70 +++++++++---------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index b041410cbb1..8631777372a 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -273,7 +273,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) { - int index_prim = clusters[i].planes[j]; + std::size_t index_prim = clusters[i].planes[j]; Point pt_reg = m_planes[index_prim]->projection (m_centroids[index_prim]); if( m_planes[index_prim]->plane_normal () * vec_reg < 0) vec_reg=-vec_reg; @@ -293,11 +293,11 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip) clusters[i].coplanar_group.push_back (-1); - int cop_index=0; + std::size_t cop_index=0; for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) { - int index_prim = clusters[i].planes[j]; + std::size_t index_prim = clusters[i].planes[j]; if (clusters[i].coplanar_group[j] == static_cast(-1)) { @@ -310,10 +310,10 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. { if (clusters[i].coplanar_group[k] == static_cast(-1)) { - int index_prim_next = clusters[i].planes[k]; + std::size_t index_prim_next = clusters[i].planes[k]; Point pt_reg_next = m_planes[index_prim_next]->projection(m_centroids[index_prim_next]); Point pt_proj=plan_reg.projection(pt_reg_next); - double distance=distance_Point(pt_reg_next,pt_proj); + FT distance=distance_Point(pt_reg_next,pt_proj); if (distance < tolerance_coplanarity) clusters[i].coplanar_group[k] = cop_index; @@ -325,7 +325,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. //regularize primitive position by computing barycenter of cplanar planes std::vector pt_bary (cop_index, Point (0., 0., 0.)); - std::vector area (cop_index, 0.); + std::vector area (cop_index, 0.); for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) { @@ -368,16 +368,16 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. std::vector < Point > listp; for (std::size_t j = 0; j < m_planes[i]->indices_of_assigned_points ().size (); ++ j) { - int yy = m_planes[i]->indices_of_assigned_points()[j]; + std::size_t yy = m_planes[i]->indices_of_assigned_points()[j]; Point pt = get (m_point_pmap, *(m_input_begin + yy)); listp.push_back(pt); } m_centroids.push_back (CGAL::centroid (listp.begin (), listp.end ())); - m_areas.push_back ((double)(m_planes[i]->indices_of_assigned_points().size()) / 100.); + m_areas.push_back ((FT)(m_planes[i]->indices_of_assigned_points().size()) / 100.); } } - void compute_parallel_clusters (std::vector& clusters, double tolerance_cosangle, + void compute_parallel_clusters (std::vector& clusters, FT tolerance_cosangle, const Vector& symmetry_direction) { // find pairs of epsilon-parallel primitives and store them in parallel_planes @@ -478,21 +478,21 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. is_available.clear(); } - void cluster_symmetric_cosangles (std::vector& clusters, double tolerance_cosangle) + void cluster_symmetric_cosangles (std::vector& clusters, FT tolerance_cosangle) { - std::vector < double > cosangle_centroids; - std::vector < int > list_cluster_index; + 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(-1); - int mean_index = 0; + std::size_t mean_index = 0; for (std::size_t i = 0; i < clusters.size(); ++ i) { if(list_cluster_index[i]<0) { list_cluster_index[i] = mean_index; - double mean = clusters[i].area * clusters[i].cosangle_symmetry; - double mean_area = clusters[i].area; + 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) { @@ -524,8 +524,8 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. void subgraph_mutually_orthogonal_clusters (std::vector& clusters, const Vector& symmetry_direction) { - std::vector < std::vector < int > > subgraph_clusters; - std::vector < int > subgraph_clusters_max_area_index; + 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; @@ -535,15 +535,15 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if(clusters[i].is_free) { clusters[i].is_free = false; - double max_area = clusters[i].area; - int index_max_area = i; + FT max_area = clusters[i].area; + std::size_t index_max_area = i; //initialization containers - std::vector < int > index_container; + std::vector < std::size_t > index_container; index_container.push_back(i); - std::vector < int > index_container_former_ring; + std::vector < std::size_t > index_container_former_ring; index_container_former_ring.push_back(i); - std::list < int > index_container_current_ring; + std::list < std::size_t > index_container_current_ring; //propagation bool propagation=true; @@ -555,7 +555,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. for (std::size_t k=0;k::iterator it = index_container_current_ring.begin(); + 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); @@ -595,11 +595,11 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. //create subgraphs of mutually orthogonal clusters in which the //largest cluster is excluded and store in //subgraph_clusters_prop - std::vector < std::vector < int > > subgraph_clusters_prop; + std::vector < std::vector < std::size_t> > subgraph_clusters_prop; for (std::size_t i=0;i subgraph_clusters_prop_temp; + 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 index_container; + std::vector < std::size_t> index_container; index_container.push_back(index_current); - std::vector < int > index_container_former_ring; + std::vector < std::size_t> index_container_former_ring; index_container_former_ring.push_back(index_current); - std::list < int > index_container_current_ring; + std::list < std::size_t> index_container_current_ring; //propagation bool propagation=true; @@ -644,7 +644,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. for (std::size_t k=0;k::iterator it = index_container_current_ring.begin(); + 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); @@ -703,7 +703,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (!assign(line, ob_1)) return n; - double delta = std::sqrt (1 - cos_symmetry * cos_symmetry); + FT delta = std::sqrt (1 - cos_symmetry * cos_symmetry); Point projected_origin = line.projection (CGAL::ORIGIN); Vector line_vector (line); @@ -745,7 +745,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (R <= 1) // 2 (or 1) possible points intersecting the unit sphere and line { - double delta = std::sqrt (1 - R); + FT delta = std::sqrt (1 - R); Vector line_vector(line); line_vector = line_vector / std::sqrt (line_vector * line_vector); Point pt1 = projected_origin + delta * line_vector; From 4cd4d02b183ae12aaa93ee0929cfdd198098b6fd Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 13:42:24 +0100 Subject: [PATCH 30/48] Declare floats explicitly to avoid warnings on Visual Studio --- .../include/CGAL/Plane_regularization.h | 12 ++++++------ .../test/Point_set_shape_detection_3/test_scene.cpp | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 8631777372a..178fe336075 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -221,14 +221,14 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. \return The number of clusters of parallel planes found. */ - std::size_t run (FT tolerance_angle = 25.0, - FT tolerance_coplanarity = 0.0, + std::size_t run (FT tolerance_angle = (FT)25.0, + FT tolerance_coplanarity = (FT)0.0, bool regularize_orthogonality = true, Vector symmetry_direction = CGAL::NULL_VECTOR) { compute_centroids_and_areas (); - FT tolerance_cosangle = 1. - std::cos (tolerance_angle); + FT tolerance_cosangle = 1.f - std::cos (tolerance_angle); // clustering the parallel primitives and store them in clusters // & compute the normal, size and cos angle to the symmetry @@ -451,7 +451,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. clu.normal = (FT)clu.area * clu.normal + (FT)m_areas[it] * normal_it; - FT norm = 1. / std::sqrt (clu.normal.squared_length()); + FT norm = 1.f / std::sqrt (clu.normal.squared_length()); clu.normal = norm * clu.normal; clu.area += m_areas[it]; } @@ -703,7 +703,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (!assign(line, ob_1)) return n; - FT delta = std::sqrt (1 - cos_symmetry * cos_symmetry); + FT delta = std::sqrt (1.f - cos_symmetry * cos_symmetry); Point projected_origin = line.projection (CGAL::ORIGIN); Vector line_vector (line); @@ -745,7 +745,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (R <= 1) // 2 (or 1) possible points intersecting the unit sphere and line { - FT delta = std::sqrt (1 - R); + FT delta = std::sqrt (1.f - R); Vector line_vector(line); line_vector = line_vector / std::sqrt (line_vector * line_vector); Point pt1 = projected_origin + delta * line_vector; diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index 5282167816a..9330a007483 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -131,7 +131,7 @@ bool test_scene() { // Test regularization Regularization regularization (points, ransac); - regularization.run (50., 0.01); + regularization.run ((FT)50., (FT)0.01f); //, true, typename K::Vector_3(1., 0., 0.)); Point_index_range pts = ransac.indices_of_unassigned_points(); From 91aa735745eacdc43ab1f419222cb094aeef84a7 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 13:52:46 +0100 Subject: [PATCH 31/48] Bugfix: assigned -1 to unsigned and test if <0 never works --- .../include/CGAL/Plane_regularization.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 178fe336075..0f1754d89f2 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -291,7 +291,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. Vector vec_reg = clusters[i].normal; for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip) - clusters[i].coplanar_group.push_back (-1); + clusters[i].coplanar_group.push_back (static_cast(-1)); std::size_t cop_index=0; @@ -483,12 +483,12 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. 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(-1); + list_cluster_index.push_back(static_cast(-1)); std::size_t mean_index = 0; for (std::size_t i = 0; i < clusters.size(); ++ i) { - if(list_cluster_index[i]<0) + if(list_cluster_index[i] == static_cast(-1)) { list_cluster_index[i] = mean_index; FT mean = clusters[i].area * clusters[i].cosangle_symmetry; @@ -496,8 +496,9 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. for (std::size_t j = i+1; j < clusters.size(); ++ j) { - if (list_cluster_index[j] < 0 && std::fabs (clusters[j].cosangle_symmetry - - mean / mean_area) < tolerance_cosangle) + if (list_cluster_index[j] == static_cast(-1) + && std::fabs (clusters[j].cosangle_symmetry - + mean / mean_area) < tolerance_cosangle) { list_cluster_index[j] = mean_index; mean_area += clusters[j].area; From c51a9f5860b3c5fc17857c3ecfa3795816e16d74 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 13:53:35 +0100 Subject: [PATCH 32/48] Also test regularization of symmetry --- .../test/Point_set_shape_detection_3/test_scene.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index 9330a007483..ab494cc8320 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -131,7 +131,7 @@ bool test_scene() { // Test regularization Regularization regularization (points, ransac); - regularization.run ((FT)50., (FT)0.01f); //, true, typename K::Vector_3(1., 0., 0.)); + regularization.run ((FT)50., (FT)0.01f, true, typename K::Vector_3(1., 0., 0.)); Point_index_range pts = ransac.indices_of_unassigned_points(); From 1e079ed272f05bbb044090873c84d68c292fe59a Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2015 15:05:03 +0100 Subject: [PATCH 33/48] More on documentation --- .../Point_set_shape_detection_3.txt | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt index 27271847c4c..1306d3553e9 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt @@ -109,7 +109,17 @@ which is used by the example \ref Point_set_shape_detection_3/efficient_RANSAC_ Shape detection is very suited for man-made shapes such as urban scenes or scans of mechanical pieces. In such input data, shapes often come with specific relationships between them: parallelism, coplanarity or orthogonality, for example. \cgal provides a tool to regularize planes detected on a point set. -The class `CGAL::Plane_regularization` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). +The class `CGAL::Plane_regularization` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). More specifically: + +- Planes that are almost parallel are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal. + +- Planes that are almost coplanar are made coplanar. + +- Planes that are almost orthogonal are made exactly orthogonal. + +- Planes that are almost 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 depending on the chosen parameters (see reference manual). From 467e5992579f7e87e089b916577ba78dbb29ba24 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 4 Jan 2016 11:39:05 +0100 Subject: [PATCH 34/48] Explicit conversion of double/floats to FT --- .../include/CGAL/Plane_regularization.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 0f1754d89f2..09bc957c526 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -228,7 +228,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. { compute_centroids_and_areas (); - FT tolerance_cosangle = 1.f - std::cos (tolerance_angle); + 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 @@ -451,7 +451,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. clu.normal = (FT)clu.area * clu.normal + (FT)m_areas[it] * normal_it; - FT norm = 1.f / std::sqrt (clu.normal.squared_length()); + FT norm = (FT)1. / std::sqrt (clu.normal.squared_length()); clu.normal = norm * clu.normal; clu.area += m_areas[it]; } @@ -704,7 +704,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (!assign(line, ob_1)) return n; - FT delta = std::sqrt (1.f - cos_symmetry * cos_symmetry); + FT delta = std::sqrt ((FT)1. - cos_symmetry * cos_symmetry); Point projected_origin = line.projection (CGAL::ORIGIN); Vector line_vector (line); @@ -746,7 +746,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. if (R <= 1) // 2 (or 1) possible points intersecting the unit sphere and line { - FT delta = std::sqrt (1.f - R); + 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; From a0853b21a810c5dee512e147730d9fbabf182fda Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 4 Jan 2016 11:54:38 +0100 Subject: [PATCH 35/48] Authors --- Point_set_shape_detection_3/include/CGAL/Plane_regularization.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 09bc957c526..165393fdf82 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -16,7 +16,7 @@ // $Id$ // // -// Author(s) : +// Author(s) : Florent Lafarge, Simon Giraudot // /** From bd361eaa40c4513aca361e1728b07dbf2d8b72ba Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 4 Jan 2016 12:07:24 +0100 Subject: [PATCH 36/48] Update changes.html --- Installation/changes.html | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Installation/changes.html b/Installation/changes.html index 4c8ac5940ae..b7cf1e9474b 100644 --- a/Installation/changes.html +++ b/Installation/changes.html @@ -218,6 +218,13 @@ and src/ directories). CGAL::read_ply_points_and_normals(), CGAL::write_ply_points() and CGAL::write_ply_points_and_normals(). +

Point Set Shape Detection

+
    +
  • New post-processing + algorithm: Plane_regularization. This allows the user + to favor parallelism, orthogonality, coplanarity and/or axial + symmetry between detected planes.
    • +

Surface Mesh Parameterization

  • LSCM_parameterizer_3 now uses by default Eigen From 3ba780f49c8a64f0f0ff6ac4b158efc493f9c599 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 5 Jan 2016 10:33:45 +0100 Subject: [PATCH 37/48] Warning fix: explicit conversion of float to FT --- Point_set_shape_detection_3/include/CGAL/Plane_regularization.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 165393fdf82..33ceb102aec 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -373,7 +373,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. listp.push_back(pt); } m_centroids.push_back (CGAL::centroid (listp.begin (), listp.end ())); - m_areas.push_back ((FT)(m_planes[i]->indices_of_assigned_points().size()) / 100.); + m_areas.push_back ((FT)(m_planes[i]->indices_of_assigned_points().size()) / (FT)100.); } } From cf7d88c475c60835692767d7a0fbf6f50b64fb3f Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 6 Jan 2016 08:08:26 +0100 Subject: [PATCH 38/48] Explicit cast to FT for constructor point --- Point_set_shape_detection_3/include/CGAL/Plane_regularization.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 33ceb102aec..46dece94a9a 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -324,7 +324,7 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. } //regularize primitive position by computing barycenter of cplanar planes - std::vector pt_bary (cop_index, Point (0., 0., 0.)); + std::vector pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.)); std::vector area (cop_index, 0.); for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j) From 03decddb0bfbc811406f61a572b35e8238a5a9df Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 11 Mar 2016 15:29:15 +0100 Subject: [PATCH 39/48] Plane_regularization class is now a simple function regularize_planes --- .../include/CGAL/Plane_regularization.h | 1267 ++++++++--------- 1 file changed, 621 insertions(+), 646 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h index 46dece94a9a..b04636f3f4b 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h @@ -36,604 +36,340 @@ namespace CGAL { + +namespace internal { +namespace PlaneRegularization { + +template +struct Plane_cluster +{ + bool is_free; + std::vector planes; + std::vector coplanar_group; + std::vector orthogonal_clusters; + typename Traits::Vector_3 normal; + typename Traits::FT cosangle_symmetry; + typename Traits::FT area; + typename Traits::FT cosangle_centroid; +}; -/*! -\ingroup PkgPointSetShapeDetection3 -\brief A plane regularization algorithm applied as a post-processing -to a shape detection algorithm. +template +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; -Given a set of detected planes with their respective inlier sets, this -class enables to regularize the planes: planes almost parallel are -made exactly parallel. In addition, some additional regularization can -be performed: - -- Plane clusters that are almost orthogonal can be made exactly - orthogonal. - -- Planes that are parallel and almost coplanar can be made exactly - coplanar. - -- Planes that are almost 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` - -*/ - template - class Plane_regularization - { - public: - - /// \cond SKIP_IN_MANUAL - typedef Plane_regularization Self; - /// \endcond - - /// \name Types - /// @{ - /// \cond SKIP_IN_MANUAL - typedef typename Traits::FT FT; - typedef typename Traits::Point_3 Point; - typedef typename Traits::Vector_3 Vector; - typedef typename Traits::Line_3 Line; - - /// \endcond - typedef typename Traits::Plane_3 Plane; ///< Raw plane type - - typedef typename Traits::Point_map Point_map; - ///< property map to access the location of an input point. - typedef typename Traits::Normal_map Normal_map; - ///< property map to access the unoriented normal of an input point - typedef typename Traits::Input_range Input_range; - ///< Model of the concept `Range` with random access iterators, providing input points and normals - /// through the following two property maps. - - typedef typename Input_range::iterator Input_iterator; ///< Iterator on input data - - typedef Shape_detection_3::Shape_base Shape; ///< Shape type. - typedef Shape_detection_3::Plane Plane_shape; ///< Plane type. - /// @} - - private: - - struct Plane_cluster - { - bool is_free; - std::vector planes; - std::vector coplanar_group; - std::vector orthogonal_clusters; - Vector normal; - FT cosangle_symmetry; - FT area; - FT cosangle_centroid; - }; - - Traits m_traits; - - Input_iterator m_input_begin; - Input_iterator m_input_end; - Point_map m_point_pmap; - Normal_map m_normal_pmap; - - std::vector > m_planes; - std::vector m_centroids; - std::vector m_areas; - - public: - - /// \name Initialization - /// @{ - /*! - Constructs an empty plane regularization engine. - */ - Plane_regularization (Traits t = Traits ()) - : m_traits (t) - { - - } - - /*! - - Constructs a plane regularization engine base on an input range - of points with its related shape detection engine. - - \param input_range Range of input data. - - \param shape_detection Shape detection engine 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 and must have been using `input_range` as input. - - */ - Plane_regularization (Input_range& input_range, - const Shape_detection_3::Efficient_RANSAC& shape_detection) - : m_traits (shape_detection.traits()) - { - m_input_begin = input_range.begin (); - m_input_end = input_range.end (); - - BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) - { - boost::shared_ptr pshape - = boost::dynamic_pointer_cast(shape); - - // Ignore all shapes other than plane - if (pshape == boost::shared_ptr()) - continue; - m_planes.push_back (pshape); - } - - } - - /*! - Releases all memory allocated by this instance. - */ - virtual ~Plane_regularization () - { - clear (); - } - /// @} - - /// \name Memory Management - /// @{ - /*! - Clear all internal structures. - */ - void clear () - { - std::vector > ().swap (m_planes); - std::vector ().swap (m_centroids); - std::vector ().swap (m_areas); - } - /// @} - - /// \name Regularization - /// @{ - /*! - - Performs the plane regularization. Planes are directly modified. - - \param tolerance_angle Tolerance of deviation between normal - vectors of planes so that they are considered parallel (in - degrees). - - \param tolerance_coplanarity Maximal distance between two - parallel planes such that they are considered coplanar. The - default value is 0, meaning that coplanarity is not taken into - account for regularization. - - \param regularize_orthogonality Make almost orthogonal clusters - of plane exactly orthogonal. - - \param symmetry_direction Make clusters that are almost - symmetrical in the symmetry direction exactly symmetrical. This - parameter is ignored if it is equal to `CGAL::NULL_VECTOR` - (default value). - - \return The number of clusters of parallel planes found. - */ - - std::size_t run (FT tolerance_angle = (FT)25.0, - FT tolerance_coplanarity = (FT)0.0, - bool regularize_orthogonality = true, - Vector symmetry_direction = CGAL::NULL_VECTOR) - { - compute_centroids_and_areas (); - - FT tolerance_cosangle = (FT)1. - std::cos (tolerance_angle); + if (symmetry_direction == CGAL::NULL_VECTOR) + return n; - // 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 clusters; - compute_parallel_clusters (clusters, tolerance_cosangle, symmetry_direction); + Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction; - 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); - } - } - } - } - - //clustering the symmetry cosangle and store their centroids in - //cosangle_centroids and the centroid index of each cluster in - //list_cluster_index - if (symmetry_direction != CGAL::NULL_VECTOR) - cluster_symmetric_cosangles (clusters, tolerance_cosangle); + Plane plane_symmetry (pt_symmetry, symmetry_direction); + Point pt_normal = CGAL::ORIGIN + n; - //find subgraphs of mutually orthogonal clusters (store index of - //clusters in subgraph_clusters), and select the cluster of - //largest area - if (regularize_orthogonality) - subgraph_mutually_orthogonal_clusters (clusters, symmetry_direction); - - //recompute optimal plane for each primitive after normal regularization - for (std::size_t i=0; i < clusters.size(); ++ i) - { + 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; - Vector vec_reg = clusters[i].normal; + FT delta = std::sqrt ((FT)1. - cos_symmetry * cos_symmetry); - for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) - { - std::size_t index_prim = clusters[i].planes[j]; - Point pt_reg = m_planes[index_prim]->projection (m_centroids[index_prim]); - if( m_planes[index_prim]->plane_normal () * vec_reg < 0) - vec_reg=-vec_reg; - Plane plane_reg(pt_reg,vec_reg); - - if( std::fabs(m_planes[index_prim]->plane_normal () * plane_reg.orthogonal_vector ()) > 1. - tolerance_cosangle) - m_planes[index_prim]->update (plane_reg); - } - } + 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); - //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; + } + else + return n; +} - for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip) - clusters[i].coplanar_group.push_back (static_cast(-1)); +template +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; - std::size_t cop_index=0; + if (symmetry_direction == CGAL::NULL_VECTOR) + return n; - for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) - { - std::size_t index_prim = clusters[i].planes[j]; + Plane plane_orthogonality (CGAL::ORIGIN, np); + Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction; - if (clusters[i].coplanar_group[j] == static_cast(-1)) - { - clusters[i].coplanar_group[j] = cop_index; + 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 (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_reg = m_planes[index_prim]->projection(m_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(-1)) - { - std::size_t index_prim_next = clusters[i].planes[k]; - Point pt_reg_next = m_planes[index_prim_next]->projection(m_centroids[index_prim_next]); - Point pt_proj=plan_reg.projection(pt_reg_next); - FT distance=distance_Point(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 pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.)); - std::vector 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 = m_planes[index_prim]->projection(m_centroids[index_prim]); - - pt_bary[group] = CGAL::barycenter (pt_bary[group], area[group], pt_reg, m_areas[index_prim]); - area[group] += m_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 (m_planes[index_prim]->plane_normal () - * plane_reg.orthogonal_vector() < 0) - m_planes[index_prim]->update (plane_reg.opposite()); - else - m_planes[index_prim]->update (plane_reg); - } - } - - return clusters.size (); + 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 (n,symmetry_direction, cos_symmetry); +} - private: - - void compute_centroids_and_areas () +template +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) { - for (std::size_t i = 0; i < m_planes.size (); ++ i) + std::vector < Point > listp; + for (std::size_t j = 0; j < planes[i]->indices_of_assigned_points ().size (); ++ j) { - std::vector < Point > listp; - for (std::size_t j = 0; j < m_planes[i]->indices_of_assigned_points ().size (); ++ j) - { - std::size_t yy = m_planes[i]->indices_of_assigned_points()[j]; - Point pt = get (m_point_pmap, *(m_input_begin + yy)); - listp.push_back(pt); - } - m_centroids.push_back (CGAL::centroid (listp.begin (), listp.end ())); - m_areas.push_back ((FT)(m_planes[i]->indices_of_assigned_points().size()) / (FT)100.); + 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.); } +} - void compute_parallel_clusters (std::vector& clusters, FT tolerance_cosangle, - const Vector& symmetry_direction) + +template +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 > parallel_planes (planes.size ()); + for (std::size_t i = 0; i < planes.size (); ++ i) { - // find pairs of epsilon-parallel primitives and store them in parallel_planes - std::vector > parallel_planes (m_planes.size ()); - for (std::size_t i = 0; i < m_planes.size (); ++ i) - { - Vector v1 = m_planes[i]->plane_normal (); + Vector v1 = planes[i]->plane_normal (); - for (std::size_t j = 0; j < m_planes.size(); ++ j) - { - if (i == j) - continue; - - Vector v2 = m_planes[i]->plane_normal (); - - if (std::fabs (v1 * v2) > 1. - tolerance_cosangle) - parallel_planes[i].push_back (j); - } - } - - - std::vector is_available (m_planes.size (), true); - - for (std::size_t i = 0; i < m_planes.size(); ++ i) + for (std::size_t j = 0; j < planes.size(); ++ j) { - - 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); + if (i == j) + continue; - std::vector index_container_former_ring_parallel; - index_container_former_ring_parallel.push_back(i); + Vector v2 = planes[i]->plane_normal (); - std::list index_container_current_ring_parallel; - - //propagation over the pairs of epsilon-parallel primitives - bool propagation=true; - clu.normal = m_planes[i]->plane_normal (); - clu.area = m_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 = m_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)m_areas[it] * normal_it; - FT norm = (FT)1. / std::sqrt (clu.normal.squared_length()); - clu.normal = norm * clu.normal; - clu.area += m_areas[it]; - } - } - } - - //update containers - index_container_former_ring_parallel.clear(); - for (std::list::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); - } + if (std::fabs (v1 * v2) > 1. - tolerance_cosangle) + parallel_planes[i].push_back (j); } - is_available.clear(); } - void cluster_symmetric_cosangles (std::vector& clusters, FT tolerance_cosangle) - { - 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(-1)); + + std::vector is_available (planes.size (), true); - std::size_t mean_index = 0; - for (std::size_t i = 0; i < clusters.size(); ++ i) - { - if(list_cluster_index[i] == static_cast(-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(-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]]; - } - - void subgraph_mutually_orthogonal_clusters (std::vector& clusters, - const Vector& symmetry_direction) + for (std::size_t i = 0; i < planes.size(); ++ i) { - 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(is_available[i]) { - if(clusters[i].is_free) - { - clusters[i].is_free = false; - FT max_area = clusters[i].area; - std::size_t index_max_area = i; + is_available[i] = false; - //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::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_prop_temp; - for (std::size_t j=0;j 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; + clu.planes.push_back (i); + + std::vector index_container_former_ring_parallel; + index_container_former_ring_parallel.push_back(i); + + std::list 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::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 +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(-1)); + + std::size_t mean_index = 0; + for (std::size_t i = 0; i < clusters.size(); ++ i) + { + if(list_cluster_index[i] == static_cast(-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(-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 +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; @@ -649,23 +385,21 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. 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(clusters[cluster_index].normal, - clusters[j].normal, - symmetry_direction, - clusters[j].cosangle_symmetry); - clusters[j].normal = new_vect; - } - } + 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(); @@ -675,97 +409,338 @@ The implementation follows \cgalCite{cgal:vla-lod-15}. index_container.push_back(*it); } index_container_current_ring.clear(); - }while(propagation); + + } + 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_prop_temp; + for (std::size_t j=0;j + (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 + (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); + } +} - FT distance_Point (const Point& a, const Point& b) - { - return std::sqrt (CGAL::squared_distance (a, b)); - } + +} // namespace PlaneRegularization +} // namespace internal + - Vector regularize_normal (const Vector& n, const Vector& symmetry_direction, - FT cos_symmetry) + /*! + + Given a set of detected planes with their respective inlier sets, + this function enables to regularize the planes: planes almost + parallel are made exactly parallel. In addition, some additional + regularization can be performed: + + - Plane clusters that are almost orthogonal can be made exactly + orthogonal. + + - Planes that are parallel and almost coplanar can be made exactly + coplanar. + + - Planes that are almost 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 engine 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 and must have been using `input_range` as input. + + \param tolerance_angle Tolerance of deviation between normal + vectors of planes so that they are considered parallel (in + degrees). + + \param tolerance_coplanarity Maximal distance between two + parallel planes such that they are considered coplanar. The + default value is 0, meaning that coplanarity is not taken into + account for regularization. + + \param regularize_orthogonality Make almost orthogonal clusters + of plane exactly orthogonal. + + \param symmetry_direction Make clusters that are almost + symmetrical in the symmetry direction exactly symmetrical. This + parameter is ignored if it is equal to `CGAL::NULL_VECTOR` + (default value). + + \return The number of clusters of parallel planes found. +*/ + +template +void regularize_planes (RandomAccessIterator input_begin, + RandomAccessIterator /*input_end*/, + const Shape_detection_3::Efficient_RANSAC& shape_detection, + typename EfficientRANSACTraits::FT tolerance_angle + = (typename EfficientRANSACTraits::FT)25.0, + typename EfficientRANSACTraits::FT tolerance_coplanarity + = (typename EfficientRANSACTraits::FT)0.0, + bool regularize_orthogonality = true, + typename EfficientRANSACTraits::Vector_3 symmetry_direction + = CGAL::NULL_VECTOR) +{ + 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 Shape; + typedef Shape_detection_3::Plane Plane_shape; + + typedef typename internal::PlaneRegularization::Plane_cluster + Plane_cluster; + + std::vector > planes; + + BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) { - if (symmetry_direction == CGAL::NULL_VECTOR) - return n; + boost::shared_ptr pshape + = boost::dynamic_pointer_cast(shape); + + // Ignore all shapes other than plane + if (pshape == boost::shared_ptr()) + continue; + planes.push_back (pshape); + } + + + /* + * Compute centroids and areas + */ + std::vector centroids; + std::vector areas; + internal::PlaneRegularization::compute_centroids_and_areas + (input_begin, planes, Point_map(), centroids, areas); + + FT tolerance_cosangle = (FT)1. - std::cos (tolerance_angle); - Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction; + // 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 clusters; + internal::PlaneRegularization::compute_parallel_clusters + (planes, clusters, areas, tolerance_cosangle, 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; - } - - - Vector regularize_normals_from_prior (const Vector& np, - const Vector& n, - const Vector& symmetry_direction, - FT cos_symmetry) + if (regularize_orthogonality) { - 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 (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 + //discovery orthogonal relationship between clusters + for (std::size_t i = 0; i < clusters.size(); ++ i) { - 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); + 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); + } + } } - else //no point intersecting the unit sphere and line - return regularize_normal (n,symmetry_direction, cos_symmetry); + } + + //clustering the symmetry cosangle and store their centroids in + //cosangle_centroids and the centroid index of each cluster in + //list_cluster_index + if (symmetry_direction != CGAL::NULL_VECTOR) + internal::PlaneRegularization::cluster_symmetric_cosangles + (clusters, tolerance_cosangle); + //find subgraphs of mutually orthogonal clusters (store index of + //clusters in subgraph_clusters), and select the cluster of + //largest area + if (regularize_orthogonality) + internal::PlaneRegularization::subgraph_mutually_orthogonal_clusters + (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); + } } - - }; + //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(-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(-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(-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 pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.)); + std::vector 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 From cae23c7fe99f3b4385cfa286d137c9114a323d5d Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 07:12:55 +0100 Subject: [PATCH 40/48] Update API of plane regularization --- .../test/Point_set_shape_detection_3/test_scene.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index ab494cc8320..7b622f4ab30 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -32,8 +32,6 @@ bool test_scene() { typedef CGAL::Shape_detection_3::Sphere Sphere; typedef CGAL::Shape_detection_3::Torus Torus; - typedef CGAL::Plane_regularization Regularization; - Pwn_vector points; // Loads point set from a file. @@ -130,8 +128,8 @@ bool test_scene() { } // Test regularization - Regularization regularization (points, ransac); - regularization.run ((FT)50., (FT)0.01f, true, typename K::Vector_3(1., 0., 0.)); + CGAL::regularize_planes (points.begin(), points.end(), ransac, + (FT)50., (FT)0.01f, true, typename K::Vector_3(1., 0., 0.)); Point_index_range pts = ransac.indices_of_unassigned_points(); From e30e2e278ffa1becebbdda9004410fc61072a868 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 09:06:45 +0100 Subject: [PATCH 41/48] Propagate name change (Plane_regularization -> regularize_planes) --- .../CGAL/{Plane_regularization.h => regularize_planes.h} | 8 ++++---- .../test/Point_set_shape_detection_3/test_scene.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) rename Point_set_shape_detection_3/include/CGAL/{Plane_regularization.h => regularize_planes.h} (99%) diff --git a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h similarity index 99% rename from Point_set_shape_detection_3/include/CGAL/Plane_regularization.h rename to Point_set_shape_detection_3/include/CGAL/regularize_planes.h index b04636f3f4b..1e318783c27 100644 --- a/Point_set_shape_detection_3/include/CGAL/Plane_regularization.h +++ b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h @@ -21,13 +21,13 @@ /** * \ingroup PkgPointSetShapeDetection3 - * \file CGAL/Plane_regularization.h + * \file CGAL/regularize_planes.h * */ -#ifndef CGAL_PLANE_REGULARIZATION_H -#define CGAL_PLANE_REGULARIZATION_H +#ifndef CGAL_REGULARIZE_PLANES_H +#define CGAL_REGULARIZE_PLANES_H #include #include @@ -745,4 +745,4 @@ void regularize_planes (RandomAccessIterator input_begin, } // namespace CGAL -#endif // CGAL_PLANE_REGULARIZATION_H +#endif // CGAL_REGULARIZE_PLANES_H diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index 7b622f4ab30..e8f975a5ed4 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -5,7 +5,7 @@ #include #include -#include +#include #include #include From cab047817a968bff6a00eb361366d20fd760df8b Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 10:16:43 +0100 Subject: [PATCH 42/48] Clarify API of plane regularization --- .../CGAL/Shape_detection_3/Efficient_RANSAC.h | 9 + .../include/CGAL/regularize_planes.h | 159 +++++++++--------- .../test_scene.cpp | 4 +- 3 files changed, 94 insertions(+), 78 deletions(-) diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Efficient_RANSAC.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Efficient_RANSAC.h index a5821ac5984..d51fe46a3fe 100644 --- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Efficient_RANSAC.h +++ b/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 diff --git a/Point_set_shape_detection_3/include/CGAL/regularize_planes.h b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h index 1e318783c27..bf3245eb71f 100644 --- a/Point_set_shape_detection_3/include/CGAL/regularize_planes.h +++ b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h @@ -563,18 +563,18 @@ void subgraph_mutually_orthogonal_clusters (PlaneClusterContainer& clusters, \return The number of clusters of parallel planes found. */ -template -void regularize_planes (RandomAccessIterator input_begin, - RandomAccessIterator /*input_end*/, - const Shape_detection_3::Efficient_RANSAC& shape_detection, +template +void regularize_planes (const Shape_detection_3::Efficient_RANSAC& 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.0, - bool regularize_orthogonality = true, + = (typename EfficientRANSACTraits::FT)0.01, typename EfficientRANSACTraits::Vector_3 symmetry_direction - = CGAL::NULL_VECTOR) + = typename EfficientRANSACTraits::Vector_3 (0., 0., 1.)) { typedef typename EfficientRANSACTraits::FT FT; typedef typename EfficientRANSACTraits::Point_3 Point; @@ -588,6 +588,8 @@ void regularize_planes (RandomAccessIterator input_begin, typedef typename internal::PlaneRegularization::Plane_cluster Plane_cluster; + typename EfficientRANSACTraits::Input_range::iterator input_begin = shape_detection.input_iterator_first(); + std::vector > planes; BOOST_FOREACH (boost::shared_ptr shape, shape_detection.shapes()) @@ -617,7 +619,9 @@ void regularize_planes (RandomAccessIterator input_begin, // direction of each cluster std::vector clusters; internal::PlaneRegularization::compute_parallel_clusters - (planes, clusters, areas, tolerance_cosangle, symmetry_direction); + (planes, clusters, areas, + (regularize_parallelism ? tolerance_cosangle : (FT)0.0), + (regularize_axis_symmetry ? symmetry_direction : CGAL::NULL_VECTOR)); if (regularize_orthogonality) { @@ -636,19 +640,20 @@ void regularize_planes (RandomAccessIterator input_begin, } } - //clustering the symmetry cosangle and store their centroids in - //cosangle_centroids and the centroid index of each cluster in - //list_cluster_index - if (symmetry_direction != CGAL::NULL_VECTOR) - internal::PlaneRegularization::cluster_symmetric_cosangles - (clusters, tolerance_cosangle); - + 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 + (clusters, tolerance_cosangle); + } + //find subgraphs of mutually orthogonal clusters (store index of //clusters in subgraph_clusters), and select the cluster of //largest area - if (regularize_orthogonality) - internal::PlaneRegularization::subgraph_mutually_orthogonal_clusters - (clusters, symmetry_direction); + internal::PlaneRegularization::subgraph_mutually_orthogonal_clusters + (clusters, symmetry_direction); //recompute optimal plane for each primitive after normal regularization for (std::size_t i=0; i < clusters.size(); ++ i) @@ -670,75 +675,77 @@ void regularize_planes (RandomAccessIterator input_begin, } - //detecting co-planarity and store in list_coplanar_prim - for (std::size_t i = 0; i < clusters.size(); ++ i) + if (regularize_coplanarity) { - 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(-1)); - - std::size_t cop_index=0; - - for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) + //detecting co-planarity and store in list_coplanar_prim + for (std::size_t i = 0; i < clusters.size(); ++ i) { - std::size_t index_prim = clusters[i].planes[j]; + Vector vec_reg = clusters[i].normal; - if (clusters[i].coplanar_group[j] == static_cast(-1)) + for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip) + clusters[i].coplanar_group.push_back (static_cast(-1)); + + std::size_t cop_index=0; + + for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j) { - clusters[i].coplanar_group[j] = cop_index; - - Point pt_reg = planes[index_prim]->projection(centroids[index_prim]); - Plane plan_reg(pt_reg,vec_reg); + std::size_t index_prim = clusters[i].planes[j]; - for (std::size_t k = j + 1; k < clusters[i].planes.size(); ++ k) + if (clusters[i].coplanar_group[j] == static_cast(-1)) { - if (clusters[i].coplanar_group[k] == static_cast(-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) { - 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 (clusters[i].coplanar_group[k] == static_cast(-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; + if (distance < tolerance_coplanarity) + clusters[i].coplanar_group[k] = cop_index; + } } + cop_index++; } - cop_index++; + } + + //regularize primitive position by computing barycenter of cplanar planes + std::vector pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.)); + std::vector 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); } } - - //regularize primitive position by computing barycenter of cplanar planes - std::vector pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.)); - std::vector 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); - } - } - + } } /// @} diff --git a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp index e8f975a5ed4..14e7b3c6742 100644 --- a/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp +++ b/Point_set_shape_detection_3/test/Point_set_shape_detection_3/test_scene.cpp @@ -128,8 +128,8 @@ bool test_scene() { } // Test regularization - CGAL::regularize_planes (points.begin(), points.end(), ransac, - (FT)50., (FT)0.01f, true, typename K::Vector_3(1., 0., 0.)); + CGAL::regularize_planes (ransac, true, true, true, true, + (FT)50., (FT)0.01f); Point_index_range pts = ransac.indices_of_unassigned_points(); From 24d1c662ce2865f73681d379a84d6a5c7195e39d Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 11:53:41 +0100 Subject: [PATCH 43/48] Update doc with new API --- .../PackageDescription.txt | 4 +- .../Point_set_shape_detection_3.txt | 12 ++-- .../include/CGAL/regularize_planes.h | 69 +++++++++++-------- 3 files changed, 49 insertions(+), 36 deletions(-) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt index 7dd6c128def..d2ff9f2880f 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/PackageDescription.txt @@ -48,7 +48,7 @@ - `CGAL::Shape_detection_3::Cone` - `CGAL::Shape_detection_3::Torus` -## Regularization Classes ## -- `CGAL::Plane_regularization` +## Functions ## +- `CGAL::regularize_planes()` */ diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt index 1306d3553e9..3e749c5408d 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt @@ -109,17 +109,17 @@ which is used by the example \ref Point_set_shape_detection_3/efficient_RANSAC_ Shape detection is very suited for man-made shapes such as urban scenes or scans of mechanical pieces. In such input data, shapes often come with specific relationships between them: parallelism, coplanarity or orthogonality, for example. \cgal provides a tool to regularize planes detected on a point set. -The class `CGAL::Plane_regularization` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). More specifically: +The function `CGAL::regularize_planes()` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). More specifically: -- Planes that are almost parallel are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal. +- Planes that are near parallel are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal. -- Planes that are almost coplanar are made coplanar. +- Planes that are near coplanar are made coplanar. -- Planes that are almost orthogonal are made exactly orthogonal. +- Planes that are near orthogonal are made exactly orthogonal. -- Planes that are almost symmetrical with respect to a user-defined axis are made symmetrical. +- 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 depending on the chosen parameters (see reference manual). +The user can choose to only regularize one or several of these 4 properties (see reference manual). diff --git a/Point_set_shape_detection_3/include/CGAL/regularize_planes.h b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h index bf3245eb71f..0697cc9262f 100644 --- a/Point_set_shape_detection_3/include/CGAL/regularize_planes.h +++ b/Point_set_shape_detection_3/include/CGAL/regularize_planes.h @@ -37,6 +37,10 @@ namespace CGAL { +// ---------------------------------------------------------------------------- +// Private section +// ---------------------------------------------------------------------------- +/// \cond SKIP_IN_MANUAL namespace internal { namespace PlaneRegularization { @@ -510,23 +514,28 @@ void subgraph_mutually_orthogonal_clusters (PlaneClusterContainer& clusters, } // 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 almost - parallel are made exactly parallel. In addition, some additional - regularization can be performed: + this function enables to regularize the planes: - - Plane clusters that are almost orthogonal can be made exactly - orthogonal. + - Planes near parallel can be made exactly parallel. - - Planes that are parallel and almost coplanar can be made exactly - coplanar. + - Planes near orthogonal can be made exactly orthogonal. - - Planes that are almost symmetrical with a user-defined axis can be - made exactly symmetrical. + - 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 @@ -536,31 +545,36 @@ void subgraph_mutually_orthogonal_clusters (PlaneClusterContainer& clusters, \tparam Traits a model of `EfficientRANSACTraits` - \param shape_detection Shape detection engine used to detect + \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. + primitive shapes but only planes will be regularized. \warning The `shape_detection` parameter must have already - detected shapes and must have been using `input_range` as input. + 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 so that they are considered parallel (in - degrees). + 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. The - default value is 0, meaning that coplanarity is not taken into - account for regularization. + \param tolerance_coplanarity Maximal distance between two parallel + planes such that they are considered coplanar. Default value is + 0.01. - \param regularize_orthogonality Make almost orthogonal clusters - of plane exactly orthogonal. - - \param symmetry_direction Make clusters that are almost - symmetrical in the symmetry direction exactly symmetrical. This - parameter is ignored if it is equal to `CGAL::NULL_VECTOR` - (default value). - - \return The number of clusters of parallel planes found. + \param symmetry_direction Chosen axis for symmetry + regularization. Default value is the Z axis. */ template @@ -747,7 +761,6 @@ void regularize_planes (const Shape_detection_3::Efficient_RANSAC Date: Mon, 14 Mar 2016 14:23:57 +0100 Subject: [PATCH 44/48] Update manual --- .../Point_set_shape_detection_3.txt | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt index 3e749c5408d..bda0b638e55 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt @@ -107,22 +107,21 @@ which is used by the example \ref Point_set_shape_detection_3/efficient_RANSAC_ \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. In such input data, shapes often come with specific relationships between them: parallelism, coplanarity or orthogonality, for example. \cgal provides a tool to regularize planes detected on a point set. +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. -The function `CGAL::regularize_planes()` can be used to postprocess the planes detected by `CGAL::Shape_detection_3` (other primitives are left unchanged). More specifically: +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` (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. +- Planes that are near __parallel__ are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal. -- Planes that are near coplanar are made coplanar. +- Parallel planes that are near __coplanar__ are made coplanar. -- Planes that are near orthogonal are made exactly orthogonal. +- Planes that are near __orthogonal__ are made exactly orthogonal. -- Planes that are near symmetrical with respect to a user-defined axis are made symmetrical. +- 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). - \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}. From e3f9bac2c78e5d3865c4199beb972f437ee31c30 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 14:42:21 +0100 Subject: [PATCH 45/48] New example for plane regularization --- .../CMakeLists.txt | 1 + .../plane_regularization.cpp | 53 +++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 Point_set_shape_detection_3/examples/Point_set_shape_detection_3/plane_regularization.cpp diff --git a/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/CMakeLists.txt b/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/CMakeLists.txt index 6e1d154f8a0..27ed2648705 100644 --- a/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/CMakeLists.txt +++ b/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() diff --git a/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/plane_regularization.cpp b/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/plane_regularization.cpp new file mode 100644 index 00000000000..239b6e8670a --- /dev/null +++ b/Point_set_shape_detection_3/examples/Point_set_shape_detection_3/plane_regularization.cpp @@ -0,0 +1,53 @@ +#include +#include +#include +#include + +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef std::pair Point_with_normal; +typedef std::vector Pwn_vector; +typedef CGAL::First_of_pair_property_map Point_map; +typedef CGAL::Second_of_pair_property_map Normal_map; + +typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits + Traits; +typedef CGAL::Shape_detection_3::Efficient_RANSAC Efficient_ransac; +typedef CGAL::Shape_detection_3::Plane 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(); + 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; +} From 27773111137d9fcafd02cd4a78dffcbf2fe387fe Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 15:25:26 +0100 Subject: [PATCH 46/48] Add example to documentation --- .../doc/Point_set_shape_detection_3/examples.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/examples.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/examples.txt index b045616177b..1a558b75778 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/examples.txt +++ b/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 */ From 64c69ebf97df1a89e1a5272c747289992b8001d0 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 15:25:38 +0100 Subject: [PATCH 47/48] Update manual --- .../Point_set_shape_detection_3.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt index bda0b638e55..f59878aa9ab 100644 --- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt +++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Point_set_shape_detection_3.txt @@ -119,8 +119,9 @@ Among all these possibilities, \cgal provides a means to regularize four propert - 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 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 From dd235c6b928a3c0f120be9566901b326106d460a Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 14 Mar 2016 15:31:30 +0100 Subject: [PATCH 48/48] Update changes.html --- Installation/changes.html | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Installation/changes.html b/Installation/changes.html index b7cf1e9474b..549942a494b 100644 --- a/Installation/changes.html +++ b/Installation/changes.html @@ -221,7 +221,7 @@ and src/ directories).

    Point Set Shape Detection

    • New post-processing - algorithm: Plane_regularization. This allows the user + algorithm: CGAL::regularize_planes(). This allows the user to favor parallelism, orthogonality, coplanarity and/or axial symmetry between detected planes.