diff --git a/.gitignore b/.gitignore
index 4c4769bf1ab..1bb95887866 100644
--- a/.gitignore
+++ b/.gitignore
@@ -373,6 +373,7 @@ Mesh_3/examples/Mesh_3/*.mesh.*
 Mesh_3/examples/Mesh_3/*.off
 Mesh_3/examples/Mesh_3/*.png
 Mesh_3/examples/Mesh_3/.*.deps
+Mesh_3/examples/Mesh_3/random-image.inr
 Mesh_3/examples/Mesh_3/Makefile
 Mesh_3/examples/Mesh_3/applications
 Mesh_3/examples/Mesh_3/cgal_test_with_cmake
diff --git a/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp b/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp
index 6893f5d9594..15400fa0d77 100644
--- a/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp
+++ b/Mesh_3/examples/Mesh_3/mesh_3D_image.cpp
@@ -1,4 +1,3 @@
-
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 #include <CGAL/Mesh_triangulation_3.h>
diff --git a/Mesh_3/examples/Mesh_3/random_labeled_image.h b/Mesh_3/examples/Mesh_3/random_labeled_image.h
index 3452b1ff004..b01a87f1edd 100644
--- a/Mesh_3/examples/Mesh_3/random_labeled_image.h
+++ b/Mesh_3/examples/Mesh_3/random_labeled_image.h
@@ -7,36 +7,47 @@ CGAL::Image_3 random_labeled_image()
   const int dim = 400;
   const unsigned char number_of_spheres = 50;
   const int max_radius_of_spheres = 10;
+  const int radius_of_big_sphere = 80;
   _image* image = _createImage(dim, dim, dim, 1,
                                1.f, 1.f, 1.f, 1,
                                WK_FIXED, SGN_UNSIGNED);
   unsigned char* ptr = (unsigned char*)(image->data);
   std::fill(ptr, ptr+dim*dim*dim, '\0');
 
+  std::ptrdiff_t center = dim / 2;
   CGAL::Random rand(0);
   for(unsigned char n = number_of_spheres; n > 0 ; --n) {
-    std::size_t i = rand.uniform_smallint(    1 + max_radius_of_spheres,
-                                          dim-2 - max_radius_of_spheres);
-    std::size_t j = rand.uniform_smallint(    1 + max_radius_of_spheres,
-                                          dim-2 - max_radius_of_spheres);
-    std::size_t k = rand.uniform_smallint(    1 + max_radius_of_spheres,
-                                          dim-2 - max_radius_of_spheres);
-    for(std::ptrdiff_t ii = - max_radius_of_spheres;
-        ii <= max_radius_of_spheres; ++ii)
+    std::size_t i, j, k;
+    do {
+      i = rand.uniform_smallint(1 + max_radius_of_spheres,
+                                dim-2 - max_radius_of_spheres);
+      j = rand.uniform_smallint(1 + max_radius_of_spheres,
+                                dim-2 - max_radius_of_spheres);
+      k = rand.uniform_smallint(1 + max_radius_of_spheres,
+                                dim-2 - max_radius_of_spheres);
+    } while ( ( CGAL::square(double(center) - i) +
+                CGAL::square(double(center) - j) +
+                CGAL::square(double(center) - k) )
+              <
+              CGAL::square(double(radius_of_big_sphere) + 4 * max_radius_of_spheres) );
+    std::ptrdiff_t radius = max_radius_of_spheres;
+    if(n==1) {
+      i = j = k = center;
+      radius = radius_of_big_sphere;
+    }
+    for(std::ptrdiff_t ii = - radius; ii <= radius; ++ii)
     {
-      for(std::ptrdiff_t jj = - max_radius_of_spheres;
-          jj <= max_radius_of_spheres; ++jj)
+      for(std::ptrdiff_t jj = - radius; jj <= radius; ++jj)
       {
-        for(std::ptrdiff_t kk = - max_radius_of_spheres;
-            kk <= max_radius_of_spheres; ++kk)
+        for(std::ptrdiff_t kk = - radius; kk <= radius; ++kk)
         {
-          if(ii*ii + jj*jj + kk*kk >
-             max_radius_of_spheres * max_radius_of_spheres) continue;
+          if(ii*ii + jj*jj + kk*kk > radius * radius) continue;
           using CGAL::IMAGEIO::static_evaluate;
           static_evaluate<unsigned char>(image, i+ii, j+jj, k+kk) = n;
         }
       }
     }
   }
+  _writeImage(image, "random-image.inr");
   return CGAL::Image_3(image);
 }