diff --git a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Concepts/EfficientRANSACTraits.h b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Concepts/EfficientRANSACTraits.h
index 4a2e9d25314..b032b6f13db 100644
--- a/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Concepts/EfficientRANSACTraits.h
+++ b/Point_set_shape_detection_3/doc/Point_set_shape_detection_3/Concepts/EfficientRANSACTraits.h
@@ -60,7 +60,9 @@ public:
   
   /*!
    * Function object type that provides
-   * `Point_3 operator()(FT x, FT y, FT z)`
+   * `Point_3 operator()(Origin p)`
+   * returning the point with 0, 0, 0 as Cartesian coordinates
+   * and `Point_3 operator()(FT x, FT y, FT z)`
    * returning the point with `x`, `y` and `z` as Cartesian coordinates.
    */
   typedef unspecified_type Construct_point_3;
diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cone.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cone.h
index 18fdc0181bd..4fac3649848 100644
--- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cone.h
+++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cone.h
@@ -361,11 +361,12 @@ namespace CGAL {
       FT max[2]) const {
 
       // Create basis d1, d2
-      Vector_3 d1 = Vector_3((FT) 0, (FT) 0, (FT) 1);
+      Vector_3 d1 = this->constr_vec(
+        ORIGIN, this->constr_pt(FT(0), FT(0), FT(1)));
       Vector_3 d2 = this->cross_pdct(m_axis, d1);
       FT l = this->sqlen(d2);
       if (l < (FT)0.0001) {
-        d1 = Vector_3((FT) 1, (FT) 0, (FT) 0);
+        d1 = this->constr_vec(ORIGIN, this->constr_pt(FT(1), FT(0), FT(0)));
         d2 = this->cross_pdct(m_axis, d1);
         l = this->sqlen(d2);
       }
diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cylinder.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cylinder.h
index 285c0c31158..c0e7c31c789 100644
--- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cylinder.h
+++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Cylinder.h
@@ -168,14 +168,15 @@ namespace CGAL {
                             FT &cluster_epsilon,
                             FT min[2],
                             FT max[2]) const {
-      Vector_3 d1 = Vector_3((FT) 0, (FT) 0, (FT) 1);
+      Vector_3 d1 = this->constr_vec(
+        ORIGIN, this->constr_pt(FT(0), FT(0), FT(1)));
       Vector_3 a = this->constr_vec(m_axis);
       a = this->scale(a, (FT)1.0 / CGAL::sqrt(this->sqlen(a)));
 
       Vector_3 d2 = this->cross_pdct(a, d1);
       FT l = this->sqlen(d2);
       if (l < (FT)0.0001) {
-        d1 = Vector_3((FT) 1, (FT) 0, (FT) 0);
+        d1 = this->constr_vec(ORIGIN, this->constr_pt(FT(1), FT(0), FT(0)));
         d2 = this->cross_pdct(this->constr_vec(m_axis), d1);
         l = this->sqlen(d2);
       }
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 c43f860cc91..64af70b6aef 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
@@ -123,7 +123,7 @@ namespace CGAL {
         + this->get_z(p1) * this->get_z(m_normal));
 
       //check deviation of the 3 normal
-      Vector_3 l_v;
+      Vector_3 l_v = this->constr_vec();
       for (std::size_t i = 0;i<3;i++) {
         l_v = this->normal(indices[i]);
 
diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h
index 73c1c82bd3d..a1c381abb96 100644
--- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h
+++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Shape_base.h
@@ -127,7 +127,8 @@ namespace CGAL {
       Determines the largest cluster of inlier points. A point belongs to a cluster
       if there is a point in the cluster closer than `cluster_epsilon` distance.
      */
-    std::size_t connected_component(std::vector<std::size_t> &indices, FT cluster_epsilon) {
+    virtual std::size_t connected_component(
+      std::vector<std::size_t> &indices, FT cluster_epsilon) {
       if (indices.size() == 0)
         return 0;
 
@@ -638,7 +639,13 @@ namespace CGAL {
     FT get_x(const Point_3& p) const { return m_traits.compute_x_3_object()(p); }
     FT get_y(const Point_3& p) const { return m_traits.compute_y_3_object()(p); }
     FT get_z(const Point_3& p) const { return m_traits.compute_z_3_object()(p); }
-      
+    
+    Point_3 constr_pt() const
+    { return m_traits.construct_point_3_object()(ORIGIN); }
+    Point_3 constr_pt(FT x, FT y, FT z) const
+    { return m_traits.construct_point_3_object()(x, y, z); }
+    Vector_3 constr_vec() const
+    { return m_traits.construct_vector_3_object()(NULL_VECTOR); }
     Vector_3 constr_vec(const Point_3& p, const Point_3& q) const
     { return m_traits.construct_vector_3_object()(p, q); }
 
diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h
index 142b9d2c6a3..93c9b6659f2 100644
--- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h
+++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Sphere.h
@@ -286,7 +286,7 @@ namespace CGAL {
                             FT &cluster_epsilon,
                             FT min[2],
                             FT max[2]) const {
-      Vector_3 axis;
+      Vector_3 axis = this->constr_vec();
       FT rad = radius();
       // Take average normal as axis
       for (std::size_t i = 0;i<indices.size();i++)
@@ -294,11 +294,12 @@ namespace CGAL {
       axis = this->scale(axis, FT(1) / CGAL::sqrt(this->sqlen(axis)));
 
       // create basis d1, d2
-      Vector_3 d1 = Vector_3((FT) 0, (FT) 0, (FT) 1);
+      Vector_3 d1 = this->constr_vec(
+        ORIGIN, this->constr_pt(FT(0), FT(0), FT(1)));
       Vector_3 d2 = this->cross_pdct(axis, d1);
       FT l = this->sqlen(d2);
       if (l < (FT)0.0001) {
-        d1 = Vector_3((FT) 1, (FT) 0, (FT) 0);
+        d1 = this->constr_vec(ORIGIN, this->constr_pt(FT(1), FT(0), FT(0)));
         d2 = this->cross_pdct(axis, d1);
         l = this->sqlen(d2);
       }
@@ -320,6 +321,7 @@ namespace CGAL {
       FT phi = atan2(this->scalar_pdct(vec, d2), this->scalar_pdct(vec, d1));
       FT x = FT(0), y = FT(0);
       concentric_mapping(phi, proj, rad, x, y);
+      CGAL_assertion( x==x && y==y); // check not nan's
 
       min[0] = max[0] = x;
       min[1] = max[1] = y;
@@ -334,6 +336,7 @@ namespace CGAL {
         phi = atan2(this->scalar_pdct(vec, d2), this->scalar_pdct(vec, d1));
 
         concentric_mapping(phi, proj, rad, x, y);
+        CGAL_assertion( x==x && y==y); // check not nan's
 
         min[0] = (std::min<FT>)(min[0], x);
         max[0] = (std::max<FT>)(max[0], x);
diff --git a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Torus.h b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Torus.h
index a55a6002852..15a04e9db4a 100644
--- a/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Torus.h
+++ b/Point_set_shape_detection_3/include/CGAL/Shape_detection_3/Torus.h
@@ -115,7 +115,7 @@ namespace CGAL {
     FT squared_distance(const Point_3 &p) const {
       const Vector_3 d = this->constr_vec(m_center, p);
       
-	    // height over symmetry plane
+      // height over symmetry plane
       const FT height = this->scalar_pdct(d, m_axis);
 
       // distance from axis in plane
@@ -204,8 +204,8 @@ namespace CGAL {
 
       // 1. center + axis
       FT majorRad1 = FLT_MAX, minorRad1 = FLT_MAX, dist1 = FLT_MAX;
-      Point_3 c1;
-      Vector_3 axis1;
+      Point_3 c1 = this->constr_pt();
+      Vector_3 axis1 = this->constr_vec();
       if (is_finite(x1) && is_finite(y1)) {
         c1 = this->transl(p[0], this->scale(n[0], x1));
         axis1 = this->constr_vec(this->transl(p[1], this->scale(n[1], y1)), c1);
@@ -219,8 +219,8 @@ namespace CGAL {
 
       // 2. center + axis
       FT majorRad2 = 0, minorRad2 = 0, dist2 = FLT_MAX;
-      Point_3 c2;
-      Vector_3 axis2;
+      Point_3 c2 = this->constr_pt();
+      Vector_3 axis2 = this->constr_vec();
       if (is_finite(x2) && is_finite(y2)) {
         c2 = this->transl(p[0], this->scale(n[0], x2));
         axis2 = this->constr_vec(this->transl(p[1], this->scale(n[1], y2)), c2);