From cb18c331b233eb16e07fb79a46d08cc8b2484cea Mon Sep 17 00:00:00 2001
From: Jocelyn MEYRON <Jocelyn.Meyron@ensimag.grenoble-inp.fr>
Date: Thu, 17 Jul 2014 16:11:52 +0200
Subject: [PATCH] Fix crash in halfspaces_intersection_3

---
 .../CGAL/dual/halfspaces_intersection_3.h     | 48 ++++++++++++++++---
 .../include/CGAL/vcm_estimate_normals.h       |  8 ++--
 2 files changed, 45 insertions(+), 11 deletions(-)

diff --git a/Convex_hull_3/include/CGAL/dual/halfspaces_intersection_3.h b/Convex_hull_3/include/CGAL/dual/halfspaces_intersection_3.h
index 591b0dc7755..250a62dd3da 100644
--- a/Convex_hull_3/include/CGAL/dual/halfspaces_intersection_3.h
+++ b/Convex_hull_3/include/CGAL/dual/halfspaces_intersection_3.h
@@ -8,7 +8,6 @@
 #include <CGAL/convex_hull_3.h>
 #include <CGAL/intersections.h>
 #include <CGAL/assertions.h>
-#include <CGAL/Point_inside_polyhedron_3.h>
 
 namespace CGAL
 {
@@ -107,9 +106,10 @@ namespace CGAL
                             // To do this, for each dual vertex, we circulate around this vertex
                             // and we add an edge between each facet we encounter
                             for (Vertex_const_iterator it = _dual.vertices_begin();
-                                 it != _dual.vertices_end(); ++it, ++n) {
+                                 it != _dual.vertices_end(); ++it) {
                                 typename Polyhedron_dual::Halfedge_around_vertex_const_circulator
                                     h0 = it->vertex_begin(), hf = h0;
+
                                 B.begin_facet();
                                 do {
                                     B.add_vertex_to_facet(primal_vertices[hf->facet()]);
@@ -120,6 +120,43 @@ namespace CGAL
                             B.end_surface();
                         }
                     };
+
+            // Functor used during the computation of the equations
+            // of the facets of a convex polyhedron
+            template <class Facet>
+            struct Plane_equation_convex_polyhedron {
+                typename Facet::Plane_3 operator()(Facet& f) {
+                    typename Facet::Halfedge_handle h = f.halfedge();
+                    typedef typename Facet::Plane_3 Plane;
+                    return Plane(h->vertex()->point(),
+                                 h->next()->vertex()->point(),
+                                 h->next()->next()->vertex()->point());
+                }
+            };
+
+            // Test if a point is inside a convex polyhedron
+            template <class Polyhedron>
+            bool point_inside_convex_polyhedron (Polyhedron &P,
+                                                 typename Polyhedron::Traits::Point_3 const& p) {
+                // Compute the equations of the facets of the polyhedron
+                typedef typename Polyhedron::Traits::Kernel K;
+                typedef typename Polyhedron::Plane_iterator Plane_iterator;
+                typedef typename Polyhedron::Facet Facet;
+
+                std::transform(P.facets_begin(), P.facets_end(), P.planes_begin(),
+                               Plane_equation_convex_polyhedron<Facet>());
+
+                // Check if the point is inside the polyhdreon
+                for (Plane_iterator pit = P.planes_begin();
+                     pit != P.planes_end();
+                     ++pit) {
+                    if (! pit->has_on_negative_side(p)) {
+                        return false;
+                    }
+                }
+
+                return true;
+            }
         } // namespace internal
     } // namespace Convex_hull_3
 
@@ -140,10 +177,9 @@ namespace CGAL
             Builder build_primal(dual_convex_hull, origin);
             P.delegate(build_primal);
 
-            // Posterior check for the origin inside the cmputed polyhedron
-            Point_inside_polyhedron_3<Polyhedron, K> is_inside(P);
-            CGAL_assertion_msg(is_inside(origin) == CGAL::ON_BOUNDED_SIDE,
-                               "halfspaces_intersection_3: origin not in the polyhedron");
+            // Posterior check if the origin is inside the computed polyhedron
+            Polyhedron Q(P);
+            CGAL_assertion_msg(!Convex_hull_3::internal::point_inside_convex_polyhedron(Q, origin), "halfspaces_intersection_3: origin not in the polyhedron");
         }
 } // namespace CGAL
 
diff --git a/Point_set_processing_3/include/CGAL/vcm_estimate_normals.h b/Point_set_processing_3/include/CGAL/vcm_estimate_normals.h
index eb1e8d1425f..f8caded7d5d 100644
--- a/Point_set_processing_3/include/CGAL/vcm_estimate_normals.h
+++ b/Point_set_processing_3/include/CGAL/vcm_estimate_normals.h
@@ -2,12 +2,14 @@
 #define CGAL_VCM_ESTIMATE_NORMALS_H
 
 #include <CGAL/property_map.h>
+#include <CGAL/assertions.h>
 
 #include <CGAL/Delaunay_triangulation_3.h>
 #include <CGAL/voronoi_covariance_3.h>
 #include <CGAL/vcm_utilities.h>
 #include <CGAL/Kd_tree.h>
 #include <CGAL/Search_traits_3.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
 #include <CGAL/Fuzzy_sphere.h>
 
 #include <Eigen/Dense>
@@ -284,11 +286,7 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
     for (ForwardIterator it = first; it != beyond; ++it) {
         Eigen::Vector3f enormal;
         Eigen::Vector3f dir;
-        if (! internal::extract_greater_eigenvector(cov[i], enormal)) {
-            std::cerr << "Error during extraction of normal: " <<
-                "the covariance matrix is not diagonalizable!\n";
-            exit(1);
-        }
+        internal::extract_greater_eigenvector(cov[i], enormal);
 
         CGAL::Vector_3<Kernel> normal(enormal[0],
                                       enormal[1],