diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
index 31dab91f2b68..0b67ab004114 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
@@ -584,7 +584,7 @@ insert_corners()
     Index p_index = domain_.index_from_corner_index(cit->first);
 
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
-      std::cerr << "** treat corner #" << CGAL::IO::oformat(p_index) << std::endl;
+      std::cerr << "\n** treat corner #" << CGAL::IO::oformat(p_index) << std::endl;
 #endif
 
     // Get weight (the ball radius is given by the 'query_size' function)
@@ -639,6 +639,13 @@ insert_point(const Bare_point& p, const Weight& w, int dim, const Index& index,
 {
   using CGAL::Mesh_3::internal::weight_modifier;
 
+#if CGAL_MESH_3_PROTECTION_DEBUG & 1
+  std::cerr << "insert_point( (" << p
+            << "), w=" << w
+            << ", dim=" << dim
+            << ", index=" << CGAL::IO::oformat(index) << ")\n";
+#endif
+
   // Convert the dimension if it was set to a negative value (marker for special balls).
   if(dim < 0)
     dim = -1 - dim;
@@ -715,6 +722,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
             << "), w=" << w
             << ", dim=" << dim
             << ", index=" << CGAL::IO::oformat(index) << ")\n";
+  std::cerr << "triangulation dimension is " << c3t3_.triangulation().dimension() << std::endl;
 #endif
   const Tr& tr = c3t3_.triangulation();
 
@@ -749,8 +757,8 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
 #endif
 
     // if sq_d < nearest_vh's weight
-    while ( cwsr(c3t3_.triangulation().point(nearest_vh), - sq_d) == CGAL::SMALLER &&
-            ! is_special(nearest_vh) )
+    while ( ! is_special(nearest_vh) &&
+            cwsr(c3t3_.triangulation().point(nearest_vh), - sq_d) == CGAL::SMALLER )
     {
       CGAL_assertion( minimal_size_ > 0 || sq_d > 0 );
 
@@ -839,7 +847,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
       std::cerr << "smart_insert_point: weight " << w
                 << " reduced to " << min_sq_d
-                << "\n (near existing point: " << nearest_point << " )\n";
+                << " (near existing point: " << nearest_point << " )\n";
 #endif
       w = min_sq_d;
       add_handle_to_unchecked = true;
@@ -860,38 +868,35 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
   else // tr.dimension() <= 2
   {
     // change size of existing balls which include p
-    bool restart = true;
-    while ( restart )
+    for ( typename Tr::Finite_vertices_iterator it = tr.finite_vertices_begin(),
+          end = tr.finite_vertices_end() ; it != end ; ++it )
     {
-      restart = false;
-      for ( typename Tr::Finite_vertices_iterator it = tr.finite_vertices_begin(),
-           end = tr.finite_vertices_end() ; it != end ; ++it )
+      const Weighted_point& it_wp = tr.point(it);
+      FT sq_d = tr.min_squared_distance(p, cp(it_wp));
+      if ( cwsr(it_wp, - sq_d) == CGAL::SMALLER )
       {
-        const Weighted_point& it_wp = tr.point(it);
-        FT sq_d = tr.min_squared_distance(p, cp(it_wp));
-        if ( cwsr(it_wp, - sq_d) == CGAL::SMALLER )
+        bool special_ball = false;
+        if(minimal_weight_ != Weight() && sq_d < minimal_weight_)
         {
-          bool special_ball = false;
-          if(minimal_weight_ != Weight() && sq_d > minimal_weight_) {
-            sq_d = minimal_weight_;
-            w = minimal_weight_;
-            special_ball = true;
-            insert_a_special_ball = true;
-          }
-          if( ! is_special(it) ) {
-            *out++ = it;
-            change_ball_size(it, sq_d, special_ball);
-            restart = true;
-          }
-          break;
+          sq_d = minimal_weight_;
+          w = minimal_weight_;
+          special_ball = true;
+          insert_a_special_ball = true;
+        }
+
+        if( ! is_special(it) ) {
+          *out++ = it;
+          change_ball_size(it, sq_d, special_ball);
         }
       }
     }
 
+    // Change w in order to be sure that no existing point will be included in (p,w)
     FT min_sq_d = w;
+#if CGAL_MESH_3_PROTECTION_DEBUG & 1
     typename Tr::Point nearest_point;
-    // Change w in order to be sure that no existing point will be included
-    // in (p,w)
+#endif
+
     for ( typename Tr::Finite_vertices_iterator it = tr.finite_vertices_begin(),
          end = tr.finite_vertices_end() ; it != end ; ++it )
     {
@@ -899,7 +904,9 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
       FT sq_d = tr.min_squared_distance(p, cp(it_wp));
       if(sq_d < min_sq_d) {
         min_sq_d = sq_d;
+#if CGAL_MESH_3_PROTECTION_DEBUG & 1
         nearest_point = c3t3_.triangulation().point(it);
+#endif
       }
     }
 
@@ -908,7 +915,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
       std::cerr << "smart_insert_point: weight " << w
                 << " reduced to " << min_sq_d
-                << "\n (near existing point: " << nearest_point << " )\n";
+                << " (near existing point: " << nearest_point << " )\n";
 #endif
       w = min_sq_d;
       add_handle_to_unchecked = true;
@@ -927,6 +934,11 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
   }
 
   if( w < minimal_weight_) {
+#if CGAL_MESH_3_PROTECTION_DEBUG & 1
+    std::cerr << "smart_insert_point: weight " << w
+              << " was smaller than minimal weight (" << minimal_weight_ << ")\n";
+#endif
+
     w = minimal_weight_;
     insert_a_special_ball = true;
   }
@@ -962,7 +974,7 @@ insert_balls_on_edges()
     if ( ! is_treated(curve_index) )
     {
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
-      std::cerr << "** treat curve #" << curve_index << std::endl;
+      std::cerr << "\n** treat curve #" << curve_index << std::endl;
 #endif
       const Bare_point& p = std::get<1>(*fit).first;
       const Bare_point& q = std::get<2>(*fit).first;
diff --git a/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h b/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h
index 584e148f7265..4bb670cbaa98 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h
@@ -418,14 +418,14 @@ tr_intersection(const typename K::Triangle_3  &t,
 
   typedef typename K::Point_3 Point_3;
 
-  typename K::Construct_vertex_3 vertex_on =
-    k.construct_vertex_3_object();
-
+  typename K::Do_intersect_3 do_intersect =
+    k.do_intersect_3_object();
   typename K::Orientation_3 orientation =
     k.orientation_3_object();
-
   typename K::Construct_point_on_3 point_on =
     k.construct_point_on_3_object();
+  typename K::Construct_vertex_3 vertex_on =
+    k.construct_vertex_3_object();
 
   typename Mesh_3::Vector_plane_orientation_3_static_filter<K> vector_plane_orient;
 
@@ -439,17 +439,24 @@ tr_intersection(const typename K::Triangle_3  &t,
   const Point_3& q = point_on(r,1);
 
   const Orientation ray_direction = vector_plane_orient(p, q, a, b, c);
-
-  if(ray_direction == COPLANAR) return result_type();
+  if(ray_direction == COPLANAR)
+    return result_type();
 
   const Orientation abcp = orientation(a,b,c,p);
+  if(abcp == COPLANAR) // p belongs to the triangle's supporting plane
+  {
+    if(do_intersect(t, p))
+      return result_type(p);
+    else
+      return result_type();
+  }
 
-  if(abcp == COPLANAR) return result_type(); // p belongs to the triangle's
-                                        // supporting plane
-
-  if(ray_direction == abcp) return result_type();
-  // The ray lies entirely in one of the two open halfspaces defined by the
-  // triangle's supporting plane.
+  if(ray_direction == abcp)
+  {
+    // The ray lies entirely in one of the two open halfspaces defined by the
+    // triangle's supporting plane.
+    return result_type();
+  }
 
   // Here we know that the ray crosses the plane (abc)
 
diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h
index 2d9bf625f57d..8c7e57e1dd93 100644
--- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h
+++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h
@@ -691,13 +691,14 @@ Construct_initial_points::operator()(OutputIterator pts,
   typename IGT::Construct_vector_3 vector = IGT().construct_vector_3_object();
 
   const Bounding_box bbox = r_domain_.tree_.bbox();
-  const Point_3 center( FT( (bbox.xmin() + bbox.xmax()) / 2),
-                        FT( (bbox.ymin() + bbox.ymax()) / 2),
-                        FT( (bbox.zmin() + bbox.zmax()) / 2) );
+  Point_3 center( FT( (bbox.xmin() + bbox.xmax()) / 2),
+                  FT( (bbox.ymin() + bbox.ymax()) / 2),
+                  FT( (bbox.zmin() + bbox.zmax()) / 2) );
 
   CGAL::Random& rng = *(r_domain_.p_rng_ != 0 ?
                         r_domain_.p_rng_ :
                         new Random(0));
+
   Random_points_on_sphere_3<Point_3> random_point(1., rng);
 
   int i = n;
@@ -727,6 +728,15 @@ Construct_initial_points::operator()(OutputIterator pts,
         % (n - i)
         % n;
 # endif
+
+      // If the source of the ray is on the surface, every ray will return its source
+      // so change the source to a random point in the bounding box
+      if(std::get<0>(intersection) == ray_shot.source())
+      {
+        center = Point_3(rng.get_double(bbox.xmin(), bbox.xmax()),
+                         rng.get_double(bbox.ymin(), bbox.ymax()),
+                         rng.get_double(bbox.zmin(), bbox.zmax()));
+      }
     }
     ++random_point;
   }
diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
index c851268390e4..a74c0ed2fe7d 100644
--- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
+++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
@@ -1690,8 +1690,8 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
 #endif
 
   // This will never happen for a dummy point
-  while(cwsr(c3t3_.triangulation().point(nearest_vh), - sq_d) == CGAL::SMALLER &&
-        ! is_special(nearest_vh))
+  while(! is_special(nearest_vh) &&
+        cwsr(c3t3_.triangulation().point(nearest_vh), - sq_d) == CGAL::SMALLER)
   {
     CGAL_assertion(minimal_size_ > 0 || sq_d > 0);