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 4d3b851fd82..c594461e9b1 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
@@ -279,6 +279,16 @@ namespace CGAL {
       return false;
     }
 
+    // U is longitude
+    virtual bool wraps_u() const {
+      return true;
+    }
+
+    // V is between caps
+    virtual bool wraps_v() const {
+      return false;
+    }
+
   private:
     FT m_angle;
     Point_3 m_apex;
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 a5273d95eab..c26d6b443e5 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
@@ -154,14 +154,12 @@ namespace CGAL {
         return;
 
       this->m_is_valid = true;
-      this->m_wrap_u = true;
     }
 
-    virtual void parameters(const std::vector<std::size_t> &indices,
-                            std::vector<std::pair<FT, FT> > &parameterSpace,
-                            FT &cluster_epsilon,
-                            FT min[2],
-                            FT max[2]) const {
+    void parameters(const std::vector<std::size_t> &indices,
+                    std::vector<std::pair<FT, FT> > &parameterSpace,                    
+                    FT min[2],
+                    FT max[2]) const {
       Vector_3 d1 = Vector_3((FT) 0, (FT) 0, (FT) 1);
       Vector_3 a = m_axis.to_vector();
       a = a * ((FT)1.0 / CGAL::sqrt(a.squared_length()));
@@ -189,120 +187,36 @@ namespace CGAL {
       length = CGAL::sqrt(vec.squared_length());
       vec = vec * (FT)1.0 / length;
 
-      FT a1 = vec * d1;
-      a1 = (a1 < (FT) -1.0) ? (FT) -1.0 : ((a1 > (FT) 1.0) ? (FT) 1.0 : a1);
-      a1 = acos(a1);
-      FT a2 = vec * d2;
-      a2 = (a2 < (FT) -1.0) ? (FT) -1.0 : ((a2 > (FT) 1.0) ? (FT) 1.0 : a2);
-      a2 = acos(a2);
+      FT a1 = acos(vec * d1);
+      FT a2 = acos(vec * d2);
 
       FT u = FT((a2 < M_PI_2) ? 2 * M_PI - a1 : a1) * m_radius;
 
       parameterSpace[0] = std::pair<FT, FT>(u, v);
 
-      min[0] = max[0] = u;
       min[1] = max[1] = v;
 
       for (std::size_t i = 0;i<indices.size();i++) {
-        vec = this->point(indices[i]) - m_point_on_axis;
-        v = vec * a;
+        Vector_3 vec = this->point(indices[i]) - m_point_on_axis;
+        FT v = vec * a;
         vec = vec - ((vec * a) * a);
         length = CGAL::sqrt(vec.squared_length());
         vec = vec * (FT)1.0 / length;
 
-        a1 = vec * d1;
-        a1 = (a1 < (FT) -1.0) ? (FT) -1.0 : ((a1 > (FT) 1.0) ? (FT) 1.0 : a1);
-        a1 = acos(a1);
-        a2 = vec * d2;
-        a2 = (a2 < (FT) -1.0) ? (FT) -1.0 : ((a2 > (FT) 1.0) ? (FT) 1.0 : a2);
-        a2 = acos(a2);
-
-        u = FT((a2 < M_PI_2) ? 2 * M_PI - a1 : a1) * m_radius;
-        min[0] = (std::min<FT>)(min[0], u);
-        max[0] = (std::max<FT>)(max[0], u);
+        FT a1 = acos(vec * d1);
+        FT a2 = acos(vec * d2);
 
+        FT u = FT((a2 < M_PI_2) ? 2 * M_PI - a1 : a1) * m_radius;
         min[1] = (std::min<FT>)(min[1], v);
         max[1] = (std::max<FT>)(max[1], v);
 
         parameterSpace[i] = std::pair<FT, FT>(u, v);
       }
 
-      // Is close to wraping around?
-      FT diff_to_full_range = min[0] + FT(M_PI * 2.0 * m_radius) - max[0];
-      if (diff_to_full_range < cluster_epsilon) {
-        m_wrap_u = true;
-        FT frac = (max[0] - min[0]) / cluster_epsilon;
-        FT trunc = floor(frac);
-        frac = frac - trunc;
-
-        if (frac < (FT) 0.5) {
-          cluster_epsilon = (max[0] - min[0]) / (trunc - (FT) 0.01);
-        }
-      }
-      else m_wrap_u = false;
-    }
-    
-    // The u coordinate corresponds to the rotation around the axis and
-    // therefore needs to be wrapped around.
-    virtual void post_wrap(const std::vector<unsigned int> &bitmap,
-      const std::size_t &u_extent,
-      const std::size_t &v_extent,
-      std::vector<unsigned int> &labels) const {
-        if (!m_wrap_u)
-          return;
-
-        // handle top index separately
-        unsigned int nw = bitmap[u_extent - 1];
-        unsigned int l = bitmap[0];
-
-        // Special case v_extent is just 1
-        if (v_extent == 1) {
-          if (nw && nw != l)
-            update_label(labels, (std::max<unsigned int>)(nw, l), l = (std::min<unsigned int>)(nw, l));
-
-          return;
-        }
-
-        unsigned int w = bitmap[2 * u_extent - 1];
-        unsigned int sw;
-
-        if (l) {
-          if (nw && nw != l)
-            update_label(labels, (std::max<unsigned int>)(nw, l), l = (std::min<unsigned int>)(nw, l));
-          else if (w && w != l)
-            update_label(labels, (std::max<unsigned int>)(w, l), l = (std::min<unsigned int>)(w, l));
-        }
-        
-        // handle mid indices
-        for (std::size_t y = 1;y<v_extent - 1;y++) {
-          l = bitmap[y * u_extent];
-          if (!l)
-            continue;
-
-          nw = bitmap[y * u_extent - 1];
-          w = bitmap[(y + 1) * u_extent - 1];
-          sw = bitmap[(y + 2) * u_extent - 1];
-
-          if (nw && nw != l)
-            update_label(labels, (std::max<unsigned int>)(nw, l), l = (std::min<unsigned int>)(nw, l));
-          if (w && w != l)
-            update_label(labels, (std::max<unsigned int>)(w, l), l = (std::min<unsigned int>)(w, l));
-          else if (sw && sw != l)
-            update_label(labels, (std::max<unsigned int>)(sw, l), l = (std::min<unsigned int>)(sw, l));
-        }
-
-        // handle last index
-        l = bitmap[(v_extent - 1) * u_extent];
-        if (!l)
-          return;
-
-        nw = bitmap[(v_extent - 1) * u_extent - 1];
-        w = bitmap[u_extent * v_extent - 1];
-
-        if (nw && nw != l)
-          update_label(labels, (std::max<unsigned int>)(nw, l), l = (std::min<unsigned int>)(nw, l));
-        else if (w && w != l)
-          update_label(labels, (std::max<unsigned int>)(w, l), l = (std::min<unsigned int>)(w, l));
+      // Due to wrapping, the u parameter 'rotation around axis' always needs
+      // to be the full extend.
+      min[0] = 0;
+      max[0] = FT(M_PI * 2.0 * m_radius);
     }
 
     virtual void squared_distance(const std::vector<std::size_t> &indices,
@@ -359,11 +273,20 @@ namespace CGAL {
       return true;
     }
 
+    // U is longitude
+    virtual bool wraps_u() const {
+      return true;
+    }
+
+    // V is between caps
+    virtual bool wraps_v() const {
+      return false;
+    }
+
   private:
     FT m_radius;
     Line_3 m_axis;
     Point_3 m_point_on_axis;
-    mutable bool m_wrap_u;
       
     /// \endcond
   };
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 26caf4fbce1..65c0f98a1b7 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
@@ -140,9 +140,8 @@ namespace CGAL {
       this->m_is_valid = true;
     }
 
-    virtual void parameters(const std::vector<std::size_t> &indices,
-                            std::vector<std::pair<FT, FT> > &parameterSpace,
-                            FT &cluster_epsilon,                    
+    virtual void parameters(const std::vector<std::size_t>& indices,
+                            std::vector<std::pair<FT, FT> >& parameterSpace,
                             FT min[2],
                             FT max[2]) const {
       // Transform first point before to initialize min/max
@@ -192,6 +191,14 @@ namespace CGAL {
       return true;
     }
 
+    virtual bool wraps_u() const {
+      return false;
+    }
+
+    virtual bool wraps_v() const {
+      return false;
+    }
+
   private:
     Point_3 m_point_on_primitive;
     Vector_3 m_base1, m_base2, m_normal;
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 4175d064ca1..510f64038a0 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
@@ -36,6 +36,11 @@
  \file Shape_base.h
  */
 
+// CODE REVIEW
+// make code more modular: connected_component()
+// use const where relevant, eg wrapU
+// initialize all variables including max
+
 
 namespace CGAL {
   namespace Shape_detection_3 {
@@ -131,113 +136,135 @@ namespace CGAL {
       m_has_connected_component = true;
       if (!this->supports_connected_component())
         return connected_component_kdTree(indices, cluster_epsilon);
-      
-      // Fetching parameters
+
       FT min[] = {0,0}, max[] = {0,0};
 
       std::vector<std::pair<FT, FT> > parameter_space;
       parameter_space.resize(indices.size());
 
-      parameters(m_indices, parameter_space, cluster_epsilon, min, max);
+      parameters(m_indices, parameter_space, min, max);
+      int i_min[2], i_max[2];
+      i_min[0] = (int) (min[0] / cluster_epsilon);
+      i_min[1] = (int) (min[1] / cluster_epsilon);
+      i_max[0] = (int) (max[0] / cluster_epsilon);
+      i_max[1] = (int) (max[1] / cluster_epsilon);
 
-      // Determine required size of bitmap
-      std::size_t u_extent = std::size_t(ceil((max[0] - min[0]) / cluster_epsilon));
-      std::size_t v_extent = std::size_t(ceil((max[1] - min[1]) / cluster_epsilon));
+      std::size_t u_extent = CGAL::abs(i_max[0] - i_min[0]) + 1;
+      std::size_t v_extent = CGAL::abs(i_max[1] - i_min[1]) + 1;
 
-      // Handle singular case
-      u_extent = (u_extent == 0) ? 1 : u_extent;
-      v_extent = (v_extent == 0) ? 1 : v_extent;
+      std::vector<std::vector<std::size_t> > bitmap;
+      std::vector<bool> visited;
+      bitmap.resize(u_extent * v_extent);
+      visited.resize(u_extent * v_extent, false);
 
-      std::vector<unsigned int> bitmap;
-      bitmap.resize(u_extent * v_extent, 0);
+      bool wrap_u = wraps_u();
+      bool wrap_v = wraps_v();
 
-      // Fill bitmap
       for (std::size_t i = 0;i<parameter_space.size();i++) {
         int u = (int)((parameter_space[i].first - min[0]) / cluster_epsilon);
         int v = (int)((parameter_space[i].second - min[1]) / cluster_epsilon);
-                
-        bitmap[v * int(u_extent) + u] = true;
-      }
 
-      // Iterate through the bitmap
-      std::vector<unsigned int> map;
-      map.reserve(64);
-      map.resize(2);
-      unsigned int label = 2;
-
-      for (std::size_t y = 0;y<v_extent;y++) {
-        for (std::size_t x = 0;x<u_extent;x++) {
-          if (!bitmap[y * u_extent + x])
-            continue;
-          unsigned int w = (x > 0) ? bitmap[y  * u_extent + x - 1] : 0;
-          unsigned int n = (y > 0) ? bitmap[(y - 1) * u_extent + x] : 0;
-          unsigned int nw = (x > 0 && y > 0) ? bitmap[(y - 1) * u_extent + x - 1] : 0;
-          unsigned int ne = ((x + 1 < u_extent) && y > 0) ? bitmap[(y - 1) * u_extent + x + 1] : 0;
-
-          // Find smallest set label;
-          unsigned int curLabel = map.size();
-          curLabel = (w != 0) ? (std::min<unsigned int>)(curLabel, w) : curLabel;
-          curLabel = (n != 0) ? (std::min<unsigned int>)(curLabel, n) : curLabel;
-          curLabel = (nw != 0) ? (std::min<unsigned int>)(curLabel, nw) : curLabel;
-          curLabel = (ne != 0) ? (std::min<unsigned int>)(curLabel, ne) : curLabel;
-
-          // Update merge map.
-          if (curLabel != map.size()) {
-            if (w > curLabel) update_label(map, w, curLabel);
-            if (nw > curLabel) update_label(map, nw, curLabel);
-            if (n > curLabel) update_label(map, n, curLabel);
-            if (ne > curLabel) update_label(map, ne, curLabel);
+        if (u < 0 || (std::size_t)u >= u_extent) {
+          if (wrap_u) {
+            while (u < 0) u += (int) u_extent;
+            while (u >= (int) u_extent) u-= (int)u_extent;
+          }
+          else {
+            u = (u < 0) ? 0 : (u >= (int) u_extent) ? (int)u_extent - 1 : u;
           }
-          else map.push_back(map.size());
-
-          bitmap[y * u_extent + x] = curLabel;
         }
+        if (v < 0 || v >= (int) v_extent) {
+          if (wrap_v) {
+            while (v < 0) v += (int) v_extent;
+            while (v >= (int) v_extent) v-= (int) v_extent;
+          }
+          else {
+            v = (v < 0) ? 0 : (v >= (int) v_extent) ? (int) v_extent - 1 : v;
+          }
+        }
+
+        bitmap[v * int(u_extent) + u].push_back(m_indices[i]);
       }
 
-      // post_wrap to handle boundaries in different shape types.
-      post_wrap(bitmap, u_extent, v_extent, map);
-      
-      // Update labels
-      for (std::size_t y = 0;y<v_extent;y++)
-        for (std::size_t x = 0;x<u_extent;x++) {
-          int label = bitmap[y * u_extent + x];
+      std::vector<std::vector<std::size_t> > cluster;
+      for (std::size_t i = 0;i<(u_extent * v_extent);i++) {
+        cluster.push_back(std::vector<std::size_t>());
+        if (bitmap[i].empty())
+          continue;
+        if (visited[i])
+          continue;
 
-          if (!label)
+        std::stack<std::size_t> fields;
+        fields.push(i);
+        while (!fields.empty()) {
+          std::size_t f = fields.top();
+          fields.pop();
+          if (visited[f])
+            continue;
+          visited[f] = true;
+          if (bitmap[f].empty())
             continue;
 
-          if (map[label] != label)
-            bitmap[y * u_extent + x] = map[label];
+          // copy indices
+          std::copy(bitmap[f].begin(), bitmap[f].end(),
+            std::back_inserter(cluster.back()));
+
+          // grow 8-neighborhood
+          int v_index = int(f / u_extent);
+          int u_index = int(f % u_extent);
+          bool upper_border = v_index == 0;
+          bool lower_border = v_index == ((int)v_extent - 1);
+          bool left_border = u_index == 0;
+          bool right_border = u_index == ((int)u_extent - 1);
+
+          int n;
+          if (!upper_border) {
+            n = int(f - u_extent);
+            if (!visited[n])
+              fields.push(n);
+          }
+          else if (wrap_v) {
+            n = int((f + v_extent - 1) * u_extent);
+            if (!visited[n]) fields.push(n);
+          }
+
+          if (!left_border) {
+            n = int(f - 1);
+            if (!visited[n]) fields.push(n);
+          }
+          else if (wrap_u) {
+            n = int(f + u_extent - 1);
+            if (!visited[n]) fields.push(n);
+          }
+
+          if (!lower_border) {
+            n = int(f + u_extent);
+            if (!visited[n]) fields.push(n);
+          }
+          else if (wrap_v) {
+            n = int((f - (v_extent - 1)) * u_extent);
+            if (!visited[n]) fields.push(n);
+          }
+
+          if (!right_border) {
+            n = int(f) + 1;
+            if (!visited[n]) fields.push(n);
+          }
+          else if (wrap_u) {
+            n = int(f - u_extent + 1);
+            if (!visited[n]) fields.push(n);
+          }
         }
-
-      // Count points per label.
-      std::vector<unsigned int> count(map.size(), 0);
-
-      for (std::size_t i = 0;i<parameter_space.size();i++) {
-        int u = (int)((parameter_space[i].first - min[0]) / cluster_epsilon);
-        int v = (int)((parameter_space[i].second - min[1]) / cluster_epsilon);
-
-        count[bitmap[v * int(u_extent) + u]]++;
       }
 
-      // Find largest component. Start at index 2 as 0/1 are reserved for
-      // basic free/occupied bitmap labels.
-      int largest = 2;
-      for (std::size_t i = 3;i<count.size();i++)
-        largest = (count[largest] < count[i]) ? i : largest;
-
-      // Extract sought-after indices.
-      std::vector<std::size_t> comp_indices;
-      comp_indices.reserve(count[largest]);
-
-      for (std::size_t i = 0;i<parameter_space.size();i++) {
-        int u = (int)((parameter_space[i].first - min[0]) / cluster_epsilon);
-        int v = (int)((parameter_space[i].second - min[1]) / cluster_epsilon);
-
-        if (bitmap[v * int(u_extent) + u] == largest)
-          comp_indices.push_back(indices[i]);
+      std::size_t max_cluster = 0;
+      for (std::size_t i = 1;i<cluster.size();i++) {
+        if (cluster[i].size() > cluster[max_cluster].size()) {
+          max_cluster = i;
+        }
       }
 
-      indices = comp_indices;
+      indices = cluster[max_cluster];
 
       return m_score = indices.size();
     }
@@ -396,17 +423,6 @@ namespace CGAL {
       return  m_upper_bound;
     }
 
-    virtual void post_wrap(const std::vector<unsigned int> &bitmap,
-                   const std::size_t &u_extent,
-                   const std::size_t &v_extent,
-                   std::vector<unsigned int> &labels) const {
-      // Avoid compiler warnings about unused parameters.
-      (void) bitmap;
-      (void) u_extent;
-      (void) v_extent;
-      (void) labels;
-    }
-
     // return last computed score, or -1 if no score yet
     FT inline score() const {
       return m_score;
@@ -421,16 +437,6 @@ namespace CGAL {
       return expected_value();
     }
 
-    void inline update_label(std::vector<unsigned int> &labels, unsigned int i, unsigned int &new_value) const {
-      if (labels[i] != i)
-        update_label(labels, labels[i], new_value);
-
-      if (new_value < labels[i])
-        labels[i] = new_value;
-      else
-        new_value = labels[i];
-    }
-
     void update_points(const std::vector<int> &shape_index) {
       if (!m_indices.size())
         return;
@@ -460,7 +466,6 @@ namespace CGAL {
 
     virtual void parameters(const std::vector<std::size_t>& indices,
                             std::vector<std::pair<FT, FT> >& parameter_space,
-                            FT &cluster_epsilon,
                             FT min[2],
                             FT max[2]) const {
       // Avoid compiler warnings about unused parameters.
@@ -552,6 +557,14 @@ namespace CGAL {
       return false;
     };
 
+    virtual bool wraps_u() const {
+      return false;
+    };
+
+    virtual bool wraps_v() const {
+      return false;
+    };
+
   protected:
     /// \endcond
     // 
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 46500a35321..ac0b3b9b806 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
@@ -223,7 +223,17 @@ namespace CGAL {
     virtual bool supports_connected_component() const {
       return false;
     }
-    
+
+    // U is longitude
+    virtual bool wraps_u() const {
+      return true;
+    }
+
+    // V is latitude
+    virtual bool wraps_v() const {
+      return false;
+    }
+
   private:
     Sphere_3 m_sphere;
 /// \endcond
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 77792308923..2306e68d356 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
@@ -340,6 +340,14 @@ namespace CGAL {
       return false;
     }
 
+    virtual bool wraps_u() const {
+      return false;
+    }
+
+    virtual bool wraps_v() const {
+      return false;
+    }
+
   private:
     FT getCircle(Point_3 &center, const Vector_3 &axis, std::vector<Point_3> p, FT &majorRad, FT &minorRad) const {
       // create spin image