Documentation updates.

Filters are moved to another file (Filters.h)
Surface_mesh_segmentation.h moved to internal folder (and namespace).

.gitattributes

@@ -4089,13 +4089,14 @@ Surface_mesh_segmentation/examples/Surface_mesh_segmentation/camel.off -text
 Surface_mesh_segmentation/examples/Surface_mesh_segmentation/sdf_values_computation_example.cpp -text
 Surface_mesh_segmentation/examples/Surface_mesh_segmentation/surface_mesh_segmentation_example.cpp -text
 Surface_mesh_segmentation/examples/Surface_mesh_segmentation/surface_mesh_segmentation_from_sdf_values_example.cpp -text
-Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h -text
 Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_traversal_traits.h -text
 Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h -text
-Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_sampling.h -text
+Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h -text
 Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h -text
+Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Filters.h -text
 Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h -text
 Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h -text
+Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h -text
 Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h -text
 Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/copyright -text
 Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/license.txt -text

Surface_mesh_segmentation/doc/Surface_Mesh_Segmentation.txt

@@ -16,7 +16,7 @@ The algorithm presented here is based on [reference to paper] which introduces S
 From the API perspective, the user can access two main outputs of the algorithm: SDF values and segmentation of the mesh. SDF values might be useful for some other purposes than segmentation such as skeletonization, and part-retrieval.
 
 # Overview of the Segmentation Process #
-The segmentation algorithm consists of three major parts. 
+The segmentation algorithm consists of three major parts: Shape Diameter Function (SDF), soft clustering, and graph-cut for hard clustering.
 
 ## Shape Diameter Function ##
 First part is computation of Shape Diameter Function (SDF), which represents a connection between the surface and its volume. More precisely, SDF is a scalar-valued function defined on the surface which measures the corresponding local volume.
@@ -34,7 +34,7 @@ Note that there is no direct relation between number of clusters (parameter for
 At the end of this part, we have soft clustering result which is a matrix that contains probability values for each facet to each cluster. These probability values are used in graph-cut part.
 
 ## Graph-Cut ##
-For hard clustering, which gives us the final partitioning of the mesh, an energy function that combines probability matrix and boundary properties is minimized. The aim of energy function is while trying to assign facets to their most probable clusters, it is also trying to minimize cost of boundary edges by considering their concavity and dihedral-angle (concave edges with large dihedral angles are considered to be less costly, which means that they are good candidates for segment boundaries).
+For hard clustering, which gives us the final partitioning of the mesh, an energy function that combines probability matrix and surface features is minimized. The aim of energy function is while trying to assign facets to their most probable clusters, it is also trying to minimize cost of boundary edges by considering their concavity and dihedral-angle (concave edges with large dihedral angles are considered to be less costly, which means that they are good candidates for segment boundaries).
 
 Energy function that is minimized using alpha-expansion graph cut algorithm is:
 
@@ -67,7 +67,7 @@ where:
 Having facets assigned to corresponding clusters, segments are determined by giving a unique id to each set of connected facets which are placed under same cluster.
 
 # API #
-The SDF computation and segmentation algorithms are implemented as free template functions. The API provides three functions listed below:
+The API provides three free template functions listed below:
 	- CGAL::sdf_values_computation : computes SDF values.
 	- CGAL::surface_mesh_segmentation_from_sdf_values : computes mesh segmentation from SDF values.
 	- CGAL::surface_mesh_segmentation : computes SDF values and mesh segmentation in one go (i.e. combines two functions listed above).	
@@ -99,5 +99,6 @@ Talk about effect of "number of levels" on segmentation results ? giving a few s
 
 Talk about pose-invariant feature of segmentation by giving a few screen shots from "benchmark" models (same models with different poses) ?
 
+Adding results of "benchmark" at the end ? Screenshots ? 
 */ 
 } /* namespace CGAL */

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h

@@ -17,13 +17,16 @@ namespace CGAL
 namespace internal
 {
 
+/**
+ * @brief Implements alpha-expansion graph cut algorithm.
+ */
 class Alpha_expansion_graph_cut
 {
-public:
+private:
-  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS, boost::directedS>
+  typedef boost::adjacency_list_traits<boost::listS, boost::vecS, boost::directedS>
   Adjacency_list_traits;
 
-  typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
+  typedef boost::adjacency_list<boost::listS, boost::vecS, boost::directedS,
           // 4 vertex properties (nested)
           boost::property<boost::vertex_index_t, int,
           boost::property<boost::vertex_color_t, boost::default_color_type,
@@ -44,18 +47,20 @@ public:
   typedef Traits::vertex_iterator Vertex_iterator;
   typedef Traits::edge_iterator Edge_iterator;
 
+public:
   Alpha_expansion_graph_cut(
     std::vector<std::pair<int, int> >& edges,
     const std::vector<double>& edge_weights,
     const std::vector<std::vector<double> >& probability_matrix,
     std::vector<int>& labels, double* result = NULL) {
-    double min_cut = apply_alpha_expansion(edges, edge_weights, probability_matrix,
+    double min_cut = apply_alpha_expansion_2(edges, edge_weights,
-                                           labels);
+                     probability_matrix, labels);
     if(result != NULL) {
       *result = min_cut;
     }
   }
 
+private:
   boost::tuple<Edge_descriptor, Edge_descriptor>
   add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1,
                        double w2, Graph& graph) {
@@ -81,38 +86,37 @@ public:
                                const std::vector<double>& edge_weights,
                                const std::vector<std::vector<double> >& probability_matrix,
                                std::vector<int>& labels) {
-    std::ofstream log_file("log_file.txt");
+    //std::ofstream log_file("log_file_2.txt");
-
+    //// logging input
-    // logging input
+    //log_file << "edges: " << std::endl;
-    log_file << "edges: " << std::endl;
+    //for(std::vector<std::pair<int, int> >::const_iterator edge_it = edges.begin(); edge_it != edges.end();
-    for(std::vector<std::pair<int, int> >::const_iterator edge_it = edges.begin();
+    //    ++edge_it)
-        edge_it != edges.end();
+    //{
-        ++edge_it) {
+    //    log_file << edge_it->first << " " << edge_it->second << std::endl;
-      log_file << edge_it->first << " " << edge_it->second << std::endl;
+    //}
-    }
+    //log_file << "edge weights: " << std::endl;
-    log_file << "edge weights: " << std::endl;
+    //for(std::vector<double>::const_iterator w_it = edge_weights.begin(); w_it != edge_weights.end(); ++w_it)
-    for(std::vector<double>::const_iterator w_it = edge_weights.begin();
+    //{
-        w_it != edge_weights.end(); ++w_it) {
+    //    log_file << (*w_it) << std::endl;
-      log_file << (*w_it) << std::endl;
+    //}
-    }
+    //log_file << "prob matrix: " << std::endl;
-    log_file << "prob matrix: " << std::endl;
+    //for(std::vector< std::vector<double> >::const_iterator v_it = probability_matrix.begin();
-    for(std::vector< std::vector<double> >::const_iterator v_it =
+    //    v_it != probability_matrix.end(); ++v_it)
-          probability_matrix.begin();
+    //{
-        v_it != probability_matrix.end(); ++v_it) {
+    //    for(std::vector<double>::const_iterator p_it = v_it->begin(); p_it != v_it->end(); ++p_it)
-      for(std::vector<double>::const_iterator p_it = v_it->begin();
+    //    {
-          p_it != v_it->end(); ++p_it) {
+    //        log_file << (*p_it) << " ";
-        log_file << (*p_it) << " ";
+    //    }
-      }
+    //    log_file << std::endl;
-      log_file << std::endl;
+    //}
-    }
+    //log_file << "labels-input:" << std::endl;
-    log_file << "labels-input:" << std::endl;
+    //for(std::vector<int>::const_iterator l_it = labels.begin(); l_it != labels.end(); ++l_it)
-    for(std::vector<int>::const_iterator l_it = labels.begin();
+    //{
-        l_it != labels.end(); ++l_it) {
+    //    log_file << (*l_it) << std::endl;
-      log_file << (*l_it) << std::endl;
+    //}
-    }
 
     ////////////////////////////////////////////////////////////
-    int number_of_clusters = probability_matrix.size();
+    const int number_of_clusters = probability_matrix.size();
     double min_cut = (std::numeric_limits<double>::max)();
     bool success;
     Timer gt;
@@ -159,10 +163,9 @@ public:
                             v2 = inserted_vertices[edge_it->second];
           int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
           if(label_1 == label_2) {
-            // actually no need for this, since two alpha labeled vertices will not be seperated
+            if(label_1 != alpha) {
-            // (their edges between sink is infitinity)
+              add_edge_and_reverse(v1, v2, *weight_it, *weight_it, graph);
-            double w1 = (label_1 == alpha) ? 0 : *weight_it;
+            }
-            add_edge_and_reverse(v1, v2, w1, w1, graph);
           } else {
             Vertex_descriptor inbetween = boost::add_vertex(graph);
 
@@ -182,7 +185,7 @@ public:
           continue;
         }
 
-        log_file << "prev flow: " << min_cut << " new flow: " << flow << std::endl;
+        // log_file << "prev flow: " << min_cut << " new flow: " << flow << std::endl;
         min_cut = flow;
         success = true;
         //update labeling
@@ -197,11 +200,11 @@ public:
 
       }
     } while(success);
-    log_file << "labels-output:" << std::endl;
+    // log_file << "labels-output:" << std::endl;
-    for(std::vector<int>::const_iterator l_it = labels.begin();
+    //for(std::vector<int>::const_iterator l_it = labels.begin(); l_it != labels.end(); ++l_it)
-        l_it != labels.end(); ++l_it) {
+    //{
-      log_file << (*l_it) << std::endl;
+    //    log_file << (*l_it) << std::endl;
-    }
+    //}
 
     std::cout << "Graph-cut time: " << gt.time() <<  std::endl;
     return min_cut;
@@ -216,7 +219,7 @@ public:
     bool success;
     Timer gt;
     gt.start();
-    double total_time = 0.0;
+
     do {
       success = false;
       for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
@@ -244,8 +247,8 @@ public:
           add_edge_and_reverse(cluster_source, new_vertex, source_weight, 0.0, graph);
           add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, 0.0, graph);
         }
-        total_time += t.time();
+
-        //std::cout << "vertex time: " << t.time() << std::endl;
+        std::cout << "vertex time: " << t.time() << std::endl;
         t.reset();
         // For E-Smooth
         // add edge between every vertex,
@@ -257,10 +260,9 @@ public:
           Vertex_descriptor v1 = edge_it->first + 2, v2 = edge_it->second + 2;
           int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
           if(label_1 == label_2) {
-            // actually no need for this, since two alpha labeled vertices will not be seperated
+            if(label_1 != alpha) {
-            // (their edges between sink is infitinity)
+              add_edge_and_reverse(v1, v2, *weight_it, *weight_it, graph);
-            double w1 = (label_1 == alpha) ? 0 : *weight_it;
+            }
-            add_edge_and_reverse(v1, v2, w1, w1, graph);
           } else {
             Vertex_descriptor inbetween = b_index++;
 
@@ -271,15 +273,14 @@ public:
             add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
           }
         }
-        total_time += t.time();
+
-        //std::cout << "edge time: " << t.time() << std::endl;
+        std::cout << "edge time: " << t.time() << std::endl;
         t.reset();
 
         double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source,
                       cluster_sink);
 
-        total_time += t.time();
+        std::cout << "flow time: " << t.time() << std::endl;
-        //std::cout << "flow time: " << t.time() << std::endl;
         t.reset();
         if(min_cut - flow < flow * 1e-10) {
           continue;

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h

@@ -1,8 +1,8 @@
-#ifndef CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLING_H
+#ifndef CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLERS_H
-#define CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLING_H
+#define CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLERS_H
 /**
- * @file Disk_sampling.h
+ * @file Disk_samplers.h
- * This file contains 3 sampling methods, which can be used as a template parameter for `SDF_calculation<Polyhedron, DiskSampling>`.
+ * This file contains 3 sampling methods, which can be used as a template parameter for SDF_calculation.
  */
 #include <CGAL/number_type_basic.h>
 
@@ -43,11 +43,11 @@ namespace internal
 // Uniform //                                  // Custom power (biased to center) //
 /**
  * @brief Uses Vogel's method to sample points from unit-disk.
- * @see Disk_sampling.h
+ * @see Disk_samplers.h, SDF_calculation
  */
 class Vogel_disk_sampling
 {
-protected:
+private:
   typedef boost::tuple<double, double, double> Disk_sample;
   typedef std::vector<Disk_sample>             Disk_samples_list;
 public:
@@ -116,7 +116,7 @@ public:
  */
 class Polar_disk_sampling
 {
-protected:
+private:
   typedef boost::tuple<double, double, double> Disk_sample;
   typedef std::vector<Disk_sample>             Disk_samples_list;
 public:
@@ -181,7 +181,7 @@ public:
  */
 class Concentric_disk_sampling
 {
-protected:
+private:
   typedef boost::tuple<double, double, double> Disk_sample;
   typedef std::vector<Disk_sample>             Disk_samples_list;
 public:
@@ -242,4 +242,4 @@ public:
 }//namespace internal
 }//namespace CGAL
 #undef CGAL_ANGLE_ST_DEV_DIVIDER
-#endif //CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLING_H
+#endif //CGAL_SURFACE_MESH_SEGMENTATION_DISK_SAMPLERS_H

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h

@@ -42,7 +42,7 @@ namespace internal
  */
 class Expectation_maximization
 {
-protected:
+private:
   /**
    * @brief Represents centers in Expectation Maximization algorithm.
    * @see Expectation_maximization
@@ -91,15 +91,14 @@ public:
   };
 
   double final_likelihood;
-protected:
+private:
-  std::vector<Gaussian_center> centers;
+  std::vector<Gaussian_center>      centers;
-  std::vector<double>  points;
+  std::vector<double>               points;
+  std::vector<std::vector<double> > responsibility_matrix;
 
   double threshold;
   int    maximum_iteration;
 
-  std::vector<std::vector<double> > responsibility_matrix;
-  unsigned int seed; /**< Seed for random initializations */
   Initialization_types init_type;
 
 public:
@@ -126,8 +125,7 @@ public:
     init_type(init_type),
     responsibility_matrix(std::vector<std::vector<double> >(number_of_centers,
                           std::vector<double>(points.size()))),
-    final_likelihood(-(std::numeric_limits<double>::max)()),
+    final_likelihood(-(std::numeric_limits<double>::max)()) {
-    seed(CGAL_DEFAULT_SEED) {
     // For initialization with k-means, with one run
     if(init_type == K_MEANS_INITIALIZATION) {
       K_means_clustering k_means(number_of_centers, data);
@@ -139,7 +137,7 @@ public:
     }
     // For initialization with random center selection, with multiple run
     else {
-      srand(seed);
+      srand(CGAL_DEFAULT_SEED);
       calculate_clustering_with_multiple_run(number_of_centers, number_of_runs);
     }
     sort(centers.begin(), centers.end());
@@ -187,29 +185,7 @@ public:
     }
   }
 
-  /** Going to be removed */
+private:
-  double calculate_distortion() {
-    double distortion = 0.0;
-
-    for(std::vector<double>::iterator it = points.begin(); it!= points.end();
-        ++it) {
-      int closest_center = 0;
-      double min_distance = std::abs(centers[0].mean - *it);
-      for(std::size_t i = 1; i < centers.size(); ++i) {
-        double distance = std::abs(centers[i].mean - *it);
-        if(distance < min_distance) {
-          min_distance = distance;
-          closest_center = i;
-        }
-      }
-      double distance = (centers[closest_center].mean - *it) /
-                        centers[closest_center].deviation;
-      distortion += distance * distance;
-    }
-    return std::sqrt(distortion / points.size());
-  }
-
-protected:
   /**
    * Calculates deviation for each center.
    * Initial deviation for a center is equal to deviation of the points whose closest center is the current center.
@@ -286,9 +262,9 @@ protected:
         distance_square_cumulative[j] = cumulative_distance_square;
       }
 
-      double random_ds = (rand() / static_cast<double>(RAND_MAX)) *
+      double zero_one = rand() / (RAND_MAX + 1.0); // [0,1) random number
-                         (distance_square_cumulative.back());
+      double random_ds =  zero_one * (distance_square_cumulative.back());
-      int selection_index = std::lower_bound(distance_square_cumulative.begin(),
+      int selection_index = std::upper_bound(distance_square_cumulative.begin(),
                                              distance_square_cumulative.end(), random_ds)
                             - distance_square_cumulative.begin();
       double initial_mean = points[selection_index];

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Filters.h

@@ -0,0 +1,213 @@
+#ifndef CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H
+#define CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H
+
+/**
+ * @file Filters.h
+ * This file contains 2 filtering methods, which can be used as a template parameter for Surface_mesh_segmentation.
+ */
+#include <vector>
+#include <map>
+#include <algorithm>
+
+namespace CGAL
+{
+namespace internal
+{
+
+/** Applies bilateral filtering on values which are associated with polyhedron facets. */
+template <class Polyhedron, class NeighborSelector = Neighbor_selector_by_edge<Polyhedron> >
+class Bilateral_filtering
+{
+public:
+  /**
+   * Bilateral filtering for values associated with facets.
+   * Takes neighbors in @a window_size, and assigns weighted average of neighbors as filtered result.
+   * For weighting two weights are multiplied:
+   *   - spatial: over geodesic distances (number of edges)
+   *   - domain : over value distances
+   * @param mesh `CGAL Polyhedron` on which @a values are defined
+   * @param window_size range of effective neighbors
+   * @param values `ReadablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
+   * @param smoothed_values `WritablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
+   */
+  template<class InputPropertyMap, class OutputPropertyMap>
+  void operator()(const Polyhedron& mesh,
+                  int window_size,
+                  InputPropertyMap values,
+                  OutputPropertyMap smoothed_values) const {
+    typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
+    typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator;
+
+    for(Facet_const_iterator facet_it = mesh.facets_begin();
+        facet_it != mesh.facets_end(); ++facet_it) {
+      std::map<Facet_const_handle, int> neighbors;
+      NeighborSelector()(facet_it, window_size,
+                         neighbors); // gather neighbors in the window
+
+      double total_sdf_value = 0.0, total_weight = 0.0;
+      double current_sdf_value = values[facet_it];
+      // calculate deviation for range weighting.
+      double deviation = 0.0;
+      for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
+          it != neighbors.end(); ++it) {
+        deviation += std::pow(values[it->first] - current_sdf_value, 2);
+      }
+      deviation = std::sqrt(deviation / neighbors.size());
+      if(deviation == 0.0) {
+        deviation = std::numeric_limits<double>::epsilon();  //this might happen
+      }
+      for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
+          it != neighbors.end(); ++it) {
+        double spatial_weight = gaussian_function(it->second,
+                                window_size / 2.0); // window_size => 2*sigma
+        double domain_weight = gaussian_function(values[it->first] -  current_sdf_value,
+                               1.5 * deviation);
+        // we can use just spatial_weight for Gauissian filtering
+        double weight = spatial_weight * domain_weight;
+
+        total_sdf_value += values[it->first] * weight;
+        total_weight += weight;
+      }
+      smoothed_values[facet_it] = total_sdf_value / total_weight;
+    }
+  }
+private:
+  /** Gauissian function for weighting. */
+  double gaussian_function(double value, double deviation) const {
+    return exp(-0.5 * (std::pow(value / deviation, 2)));
+  }
+};
+
+/** Applies median filtering on values which are associated with polyhedron facets. */
+template <class Polyhedron, class NeighborSelector = Neighbor_selector_by_vertex<Polyhedron> >
+class Median_filtering
+{
+public:
+  /**
+   * Median filtering for values associated with facets.
+   * Takes neighbors in @a window_size, and assigns median of values of neighbors as filtered result.
+   *
+   * @param mesh `CGAL Polyhedron` on which @a values are defined
+   * @param window_size range of effective neighbors
+   * @param values `ReadablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
+   * @param smoothed_values `WritablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
+   */
+  template<class InputPropertyMap, class OutputPropertyMap>
+  void operator()(const Polyhedron& mesh,
+                  int window_size,
+                  InputPropertyMap values,
+                  OutputPropertyMap smoothed_values) const {
+    typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
+    typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator;
+
+    for(Facet_const_iterator facet_it = mesh.facets_begin();
+        facet_it != mesh.facets_end(); ++facet_it) {
+      //Find neighbors and put their sdf values into a list
+      std::map<Facet_const_handle, int> neighbors;
+      NeighborSelector()(facet_it, window_size, neighbors);
+
+      std::vector<double> neighbor_values;
+      for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
+          it != neighbors.end(); ++it) {
+        neighbor_values.push_back(values[it->first]);
+      }
+      // Find median.
+      int half_neighbor_count = neighbor_values.size() / 2;
+      std::nth_element(neighbor_values.begin(),
+                       neighbor_values.begin() + half_neighbor_count, neighbor_values.end());
+      double median_sdf = neighbor_values[half_neighbor_count];
+      if( half_neighbor_count % 2 == 0) {
+        median_sdf += *std::max_element(neighbor_values.begin(),
+                                        neighbor_values.begin() + half_neighbor_count);
+        median_sdf /= 2;
+      }
+      smoothed_values[facet_it] = median_sdf;
+    }
+  }
+};
+
+
+/** Gathers neighbors of a facet for a given window range. @see Bilateral_filtering, Median_filtering */
+template<class Polyhedron>
+class Neighbor_selector_by_edge
+{
+private:
+  typedef typename Polyhedron::Facet::Halfedge_around_facet_const_circulator
+  Halfedge_around_facet_const_circulator;
+public:
+  typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
+  /**
+   * Breadth-first traversal on facets by treating facets, which share a common edge, are 1-level neighbors.
+   * @param facet root facet
+   * @param max_level maximum distance (number of levels) between root facet and visited facet
+   * @param[out] neighbors visited facets and their distances to root facet
+   */
+  void operator()(Facet_const_handle& facet, int max_level,
+                  std::map<Facet_const_handle, int>& neighbors) const {
+    typedef std::pair<Facet_const_handle, int> Facet_level_pair;
+    std::queue<Facet_level_pair> facet_queue;
+    facet_queue.push(Facet_level_pair(facet, 0));
+    while(!facet_queue.empty()) {
+      const Facet_level_pair& pair = facet_queue.front();
+      bool inserted = neighbors.insert(pair).second;
+      if(inserted && pair.second < max_level) {
+        Halfedge_around_facet_const_circulator facet_circulator =
+          pair.first->facet_begin();
+        do {
+          if(!facet_circulator->opposite()->is_border()) {
+            facet_queue.push(Facet_level_pair(facet_circulator->opposite()->facet(),
+                                              pair.second + 1));
+          }
+        } while(++facet_circulator != pair.first->facet_begin());
+      }
+      facet_queue.pop();
+    }
+  }
+};
+
+/** Gathers neighbors of a facet for a given window range. @see Bilateral_filtering, Median_filtering */
+template<class Polyhedron>
+class Neighbor_selector_by_vertex
+{
+private:
+  typedef typename Polyhedron::Facet::Halfedge_around_vertex_const_circulator
+  Halfedge_around_vertex_const_circulator;
+  typedef typename Polyhedron::Halfedge_const_iterator Halfedge_const_iterator;
+  typedef typename Polyhedron::Vertex_const_iterator   Vertex_const_iterator;
+public:
+  typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
+  /**
+   * Breadth-first traversal on facets by treating facets, which share a common vertex, are 1-level neighbors.
+   * @param facet root facet
+   * @param max_level maximum distance (number of levels) between root facet and visited facet
+   * @param[out] neighbors visited facets and their distances to root facet
+   */
+  void operator()(Facet_const_handle& facet, int max_level,
+                  std::map<Facet_const_handle, int>& neighbors) const {
+    typedef std::pair<Facet_const_handle, int> Facet_level_pair;
+    std::queue<Facet_level_pair> facet_queue;
+    facet_queue.push(Facet_level_pair(facet, 0));
+    while(!facet_queue.empty()) {
+      const Facet_level_pair& pair = facet_queue.front();
+      bool inserted = neighbors.insert(pair).second;
+      if(inserted && pair.second < max_level) {
+        Facet_const_handle facet = pair.first;
+        Halfedge_const_iterator edge = facet->halfedge();
+        do { // loop on three vertices of the facet
+          Vertex_const_iterator vertex = edge->vertex();
+          Halfedge_around_vertex_const_circulator vertex_circulator =
+            vertex->vertex_begin();
+          do { // for each vertex loop on incoming edges (through those edges loop on neighbor facets which includes the vertex)
+            if(!vertex_circulator->is_border()) {
+              facet_queue.push(Facet_level_pair(vertex_circulator->facet(), pair.second + 1));
+            }
+          } while(++vertex_circulator != vertex->vertex_begin());
+        } while((edge = edge->next()) != facet->halfedge());
+      }
+      facet_queue.pop();
+    }
+  }
+};
+}//namespace internal
+}//namespace CGAL
+#endif //CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h

@@ -24,7 +24,7 @@ namespace internal
 class K_means_clustering
 {
 // Nested classes
-protected:
+private:
   /**
    * @brief Represents centers in k-means algorithm.
    * @see K_means_point, K_means_clustering
@@ -33,7 +33,7 @@ protected:
   {
   public:
     double mean; /**< Mean of the center */
-  protected:
+  private:
     double new_mean;
     int    new_number_of_points;
 
@@ -106,12 +106,11 @@ public:
     PLUS_INITIALIZATION    /**< place initial centers using k-means++ algorithm */
   };
 
-protected:
+private:
   std::vector<K_means_center> centers;
   std::vector<K_means_point>  points;
   int  maximum_iteration;
 
-  unsigned int seed; /**< Seed for random initializations */
   Initialization_types init_type;
 
 public:
@@ -130,9 +129,8 @@ public:
                      int number_of_run = CGAL_DEFAULT_NUMBER_OF_RUN,
                      int maximum_iteration = CGAL_DEFAULT_MAXIMUM_ITERATION)
     :
-    points(data.begin(), data.end()), maximum_iteration(maximum_iteration),
+    points(data.begin(), data.end()), maximum_iteration(maximum_iteration) {
-    seed(CGAL_DEFAULT_SEED) { //seed(static_cast<unsigned int>(time(NULL)))
+    srand(CGAL_DEFAULT_SEED); //(static_cast<unsigned int>(time(NULL)))
-    srand(seed);
     calculate_clustering_with_multiple_run(number_of_centers, number_of_run);
     sort(centers.begin(), centers.end());
   }
@@ -149,7 +147,7 @@ public:
     }
   }
 
-protected:
+private:
   /**
    * Initializes centers by choosing random points from data.
    * @param number_of_centers
@@ -190,10 +188,9 @@ protected:
         cumulative_distance_square += distance_square[j];
         distance_square_cumulative[j] = cumulative_distance_square;
       }
-
+      double zero_one = rand() / (RAND_MAX + 1.0); // [0,1) random number
-      double random_ds = (rand() / static_cast<double>(RAND_MAX)) *
+      double random_ds =  zero_one * (distance_square_cumulative.back());
-                         (distance_square_cumulative.back());
+      int selection_index = std::upper_bound(distance_square_cumulative.begin(),
-      int selection_index = std::lower_bound(distance_square_cumulative.begin(),
                                              distance_square_cumulative.end(), random_ds)
                             - distance_square_cumulative.begin();
       double initial_mean = points[selection_index].data;

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h

@@ -5,7 +5,7 @@
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_polyhedron_triangle_primitive.h>
 #include <CGAL/internal/Surface_mesh_segmentation/AABB_traversal_traits.h>
-#include <CGAL/internal/Surface_mesh_segmentation/Disk_sampling.h>
+#include <CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h>
 #include <map>
 #include <vector>
 #include <algorithm>
@@ -28,24 +28,24 @@ namespace internal
  * @code
  *  // template parameters
  *  // Polyhedron          : CGAL Polyhedron (type of the mesh)
- *  // Vogel_disk_sampling : Functor type which samples rays from cone (for other functors, or for writing your own method, check Disk_sampling.h)
+ *  // Vogel_disk_sampling : Functor type which samples rays from cone (for other functors, or for writing your own method, check Disk_samplers.h)
  *
  *  // function parameters
  *  // mesh       : CGAL Polyhedron to calculate SDF
- *  // sdf_values : Writable Property-Map where results are
+ *  // sdf_values : Writable Property-Map where results are stored
  *
  *  SDF_calculation<Polyhedron, Vogel_disk_sampling>().
- *      calculate_sdf_values(120.0, 25, mesh, sdf_values);
+ *      calculate_sdf_values(mesh, sdf_values, 120.0, 25);
  * @endcode
  *
- * @see Disk_sampling.h
+ * @see Disk_samplers.h
  */
 template <class Polyhedron, class DiskSampling = Vogel_disk_sampling>
 class SDF_calculation
 {
 
 //type definitions
-protected:
+private:
   typedef typename Polyhedron::Traits Kernel;
   typedef typename Polyhedron::Facet  Facet;
   typedef typename Polyhedron::Facet  Vertex;
@@ -71,7 +71,7 @@ protected:
   typedef std::vector<Disk_sample>             Disk_samples_list;
 
 // member variables
-protected:
+private:
   double cone_angle;
 
   Disk_samples_list disk_samples_sparse;
@@ -85,21 +85,20 @@ public:
    * Assign default values to member variables.
    */
   SDF_calculation()
-    : use_minimum_segment(false),
+    : use_minimum_segment(false), multiplier_for_segment(1) {
-      multiplier_for_segment(1) {
   }
 
   /**
    * Calculates SDF values for each facet, and stores them in @a sdf_values. Note that sdf values are neither smoothed nor normalized.
    * @pre parameter @a mesh should consist of triangles.
-   * @param cone_angle opening angle for cone
+   * @param mesh `CGAL Polyhedron` on which SDF values are computed
-   * @param number_of_rays number of rays picked from cone for each facet
-   * @param mesh polyhedron that SDF values are calculated on
    * @param[out] sdf_values `WritablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
+   * @param cone_angle opening angle for cone, expressed in radians
+   * @param number_of_rays number of rays picked from cone for each facet
    */
   template <class FacetValueMap>
-  void calculate_sdf_values(double cone_angle, int number_of_rays,
+  void calculate_sdf_values(const Polyhedron& mesh, FacetValueMap sdf_values,
-                            const Polyhedron& mesh, FacetValueMap sdf_values) {
+                            double cone_angle, int number_of_rays) {
     this->cone_angle = cone_angle;
 
     const int sparse_ray_count = number_of_rays;
@@ -117,7 +116,7 @@ public:
     }
   }
 
-protected:
+private:
   /**
    * Calculates SDF value for parameter @a facet.
    * @param facet

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h

@@ -16,8 +16,9 @@
 
 #include <CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h>
 #include <CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h>
-#include <CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h>
+#include <CGAL/internal/Surface_mesh_segmentation/Filters.h>
-//#include "Alpha_expansion_graph_cut.h"
+//#include <CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h>
+#include "Alpha_expansion_graph_cut.h"
 #include <CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h>
 
 #include <CGAL/utility.h>
@@ -58,15 +59,16 @@
 
 namespace CGAL
 {
-
+namespace internal
+{
 /**
  * @brief Main entry point for mesh segmentation algorithm.
- * It is a connector class which uses
+ * It is a connector class which uses:
- *   - `SDF_calculation` for calculating sdf values
+ *   -  SDF_calculation for calculating sdf values
- *   - `Expectation_maximization` for soft clustering
+ *   -  Expectation_maximization for soft clustering
- *   - `Alpha_expansion_graph_cut` for hard clustering
+ *   -  Alpha_expansion_graph_cut for hard clustering
  *
- * Other than being a connector it is responsable for preprocess and postprocess on intermadiate data, which are:
+ * Other than being a connector, it is also responsable for preprocess and postprocess on intermadiate data, which are:
  *   - log-normalizing probabilities received from soft clustering
  *   - log-normalizing and calculating dihedral-angle based weights for edges
  *   - smoothing and log-normalizing sdf values received from sdf calculation
@@ -75,14 +77,15 @@ namespace CGAL
 template <
 class Polyhedron,
       class FacetIndexMap =
-      boost::associative_property_map<std::map<typename Polyhedron::Facet_const_handle, int> >
+      boost::associative_property_map<std::map<typename Polyhedron::Facet_const_handle, int> >,
+      class Filter = Bilateral_filtering<Polyhedron>
       >
 class Surface_mesh_segmentation
 {
-protected:
+private:
   /**
    * An adaptor for Lvalue property-map. It stores a pointer to vector for underlying data-structure,
-   * and also stores another property-map which maps `key` to vector index.
+   * and also stores another property-map which maps each `key` to vector index.
    */
   template<class Polyhedron, class ValueType, class FacetIdPropertyMap>
   struct Polyhedron_property_map_for_facet
@@ -122,13 +125,13 @@ public:
   typedef typename Polyhedron::Facet_const_handle      Facet_const_handle;
   typedef typename Polyhedron::Vertex_const_iterator   Vertex_const_iterator;
 
-  typedef internal::SDF_calculation<Polyhedron> SDF_calculation;
+  typedef SDF_calculation<Polyhedron> SDF_calculation;
 
-protected:
+private:
   typedef typename Kernel::Plane_3   Plane;
 
 // member variables
-public: // going to be protected !
+public: // going to be private !
   const Polyhedron& mesh;
 
   FacetIndexMap
@@ -138,9 +141,10 @@ public: // going to be protected !
   std::vector<double> sdf_values;
   std::vector<int>    centers;
   std::vector<int>    segments;
-// member functions
+
   bool is_pmmap_custom;
-  //std::ofstream log_file;
+
+// member functions
 public:
   Surface_mesh_segmentation(const Polyhedron& mesh)
     : mesh(mesh), facet_index_map(facet_index_map_internal),
@@ -205,14 +209,14 @@ public:
     Polyhedron_property_map_for_facet<Polyhedron, double, FacetIndexMap>
     sdf_pmap(&sdf_values, facet_index_map);
 
-    SDF_calculation().calculate_sdf_values(cone_angle, number_of_ray, mesh,
+    SDF_calculation().calculate_sdf_values(mesh, sdf_pmap, cone_angle,
-                                           sdf_pmap);
+                                           number_of_ray);
 
     SEG_DEBUG(std::cerr << "SDF computation time: " << t.time() << std::endl)
     SEG_DEBUG(t.reset())
 
     check_zero_sdf_values();
-    smooth_sdf_values_with_bilateral();
+    smooth_sdf_values();
     normalize_sdf_values();
 
     SEG_DEBUG(std::cerr <<"Normalization and smoothing time: " << t.time() <<
@@ -229,8 +233,8 @@ public:
     SEG_DEBUG(Timer t)
     SEG_DEBUG(t.start())
     // soft clustering using GMM-fitting initialized with k-means
-    internal::Expectation_maximization fitter(number_of_centers, sdf_values,
+    Expectation_maximization fitter(number_of_centers, sdf_values,
-        internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1);
+                                    Expectation_maximization::K_MEANS_INITIALIZATION, 1);
 
     std::vector<int> labels;
     fitter.fill_with_center_ids(labels);
@@ -249,8 +253,7 @@ public:
         edge_weights);
 
     // apply graph cut
-    internal::Alpha_expansion_graph_cut(edges, edge_weights, probability_matrix,
+    Alpha_expansion_graph_cut(edges, edge_weights, probability_matrix, labels);
-                                        labels);
     centers = labels;
     // assign a segment id for each facet
     assign_segments();
@@ -268,13 +271,13 @@ public:
     Polyhedron_property_map_for_facet<Polyhedron, double, FacetIndexMap>
     sdf_pmap_internal(&sdf_values, facet_index_map);
 
-    SDF_calculation().calculate_sdf_values(cone_angle, number_of_rays, mesh,
+    SDF_calculation().calculate_sdf_values(mesh, sdf_pmap_internal, cone_angle,
-                                           sdf_pmap_internal);
+                                           number_of_rays);
 
     SEG_DEBUG(std::cout << t.time() << std::endl)
 
     check_zero_sdf_values();
-    smooth_sdf_values_with_bilateral();
+    smooth_sdf_values();
     linear_normalize_sdf_values();
 
     for(Facet_const_iterator facet_it = mesh.facets_begin();
@@ -303,8 +306,8 @@ public:
     // log normalize sdf values
     normalize_sdf_values();
     // soft clustering using GMM-fitting initialized with k-means
-    internal::Expectation_maximization fitter(number_of_centers, sdf_values,
+    Expectation_maximization fitter(number_of_centers, sdf_values,
-        internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1);
+                                    Expectation_maximization::K_MEANS_INITIALIZATION, 1);
 
     std::vector<int> labels;
     fitter.fill_with_center_ids(labels);
@@ -320,8 +323,7 @@ public:
         edge_weights);
 
     // apply graph cut
-    internal::Alpha_expansion_graph_cut(edges, edge_weights, probability_matrix,
+    Alpha_expansion_graph_cut(edges, edge_weights, probability_matrix, labels);
-                                        labels);
     centers = labels;
     // assign a segment id for each facet
     assign_segments();
@@ -339,7 +341,7 @@ public:
     return segments[facet_index_map[facet]];
   }
 
-protected:
+private:
   /**
    * Going to be removed
    */
@@ -478,139 +480,23 @@ protected:
   }
 
   /**
-   * Going to be removed
+   * Smoothing sdf values.
-   */
-  void smooth_sdf_values() {
-    std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
-    for(Facet_const_iterator facet_it = mesh.facets_begin();
-        facet_it != mesh.facets_end(); ++facet_it) {
-      typename Facet::Halfedge_around_facet_const_circulator facet_circulator =
-        facet_it->facet_begin();
-      double total_neighbor_sdf = 0.0;
-      do {
-        if(!facet_circulator->opposite()->is_border()) {
-          total_neighbor_sdf += get(sdf_values, facet_circulator->opposite()->facet());
-        }
-      } while( ++facet_circulator !=  facet_it->facet_begin());
-
-      total_neighbor_sdf /= 3.0;
-      get(smoothed_sdf_values, facet_it) = (get(sdf_values,
-                                            facet_it) + total_neighbor_sdf) / 2.0;
-    }
-    sdf_values = smoothed_sdf_values;
-  }
-
-  /**
-   * Going to be removed
-   */
-  void smooth_sdf_values_with_gaussian() {
-    // take neighbors, use weighted average of neighbors as filtered result. (for weights use gaussian kernel with sigma = window_size/2)
-    const int window_size = get_window_size();
-    const int iteration = 1;
-
-    for(int i = 0; i < iteration; ++i) {
-      std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
-      for(Facet_const_iterator facet_it = mesh.facets_begin();
-          facet_it != mesh.facets_end(); ++facet_it) {
-        std::map<Facet_const_handle, int> neighbors;
-        get_neighbors_by_vertex(facet_it, neighbors, window_size);
-
-        double total_sdf_value = 0.0;
-        double total_weight = 0.0;
-        for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
-            it != neighbors.end(); ++it) {
-          double weight = gaussian_function(it->second,
-                                            window_size/2.0); // window_size => 2*sigma
-          total_sdf_value += get(sdf_values, it->first) * weight;
-          total_weight += weight;
-        }
-        get(smoothed_sdf_values, facet_it) = total_sdf_value / total_weight;
-      }
-      sdf_values = smoothed_sdf_values;
-    }
-  }
-
-  /**
-   * Going to be removed
-   */
-  void smooth_sdf_values_with_median() {
-    // take neighbors, use median sdf_value as filtered one.
-    const int window_size = get_window_size();
-    const int iteration = 1;
-
-    for(int i = 0; i < iteration; ++i) {
-      std::vector<double> smoothed_sdf_values(mesh.size_of_facets());;
-      for(Facet_const_iterator facet_it = mesh.facets_begin();
-          facet_it != mesh.facets_end(); ++facet_it) {
-        //Find neighbors and put their sdf values into a list
-        std::map<Facet_const_handle, int> neighbors;
-        get_neighbors_by_edge(facet_it, window_size, neighbors);
-        std::vector<double> sdf_of_neighbors;
-        for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
-            it != neighbors.end(); ++it) {
-          sdf_of_neighbors.push_back(get(sdf_values, it->first));
-        }
-        // Find median.
-        int half_neighbor_count = sdf_of_neighbors.size() / 2;
-        std::nth_element(sdf_of_neighbors.begin(),
-                         sdf_of_neighbors.begin() + half_neighbor_count, sdf_of_neighbors.end());
-        double median_sdf = sdf_of_neighbors[half_neighbor_count];
-        if( half_neighbor_count % 2 == 0) {
-          median_sdf += *std::max_element(sdf_of_neighbors.begin(),
-                                          sdf_of_neighbors.begin() + half_neighbor_count);
-          median_sdf /= 2;
-        }
-        get(smoothed_sdf_values, facet_it) = median_sdf;
-      }
-      sdf_values = smoothed_sdf_values;
-    }
-  }
-
-  /**
-   * Bilateral smoothing for sdf values.
    * Takes neighbors using `get_window_size()`, and assigns weighted average of neighbors as filtered result.
    * For weighting two weights are multiplied:
    *   - spatial: over geodesic distances (number of edges)
    *   - domain : over sdf value distances
    */
-  void smooth_sdf_values_with_bilateral() {
+  void smooth_sdf_values() {
-    const int window_size = get_window_size();
+    typedef Polyhedron_property_map_for_facet<Polyhedron, double, FacetIndexMap>
-    const int iteration = 1;
+    Facet_vector_map;
 
-    for(int i = 0; i < iteration; ++i) {
+    Facet_vector_map sdf_input(&sdf_values, facet_index_map);
-      std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
-      for(Facet_const_iterator facet_it = mesh.facets_begin();
-          facet_it != mesh.facets_end(); ++facet_it) {
-        std::map<Facet_const_handle, int> neighbors;
-        get_neighbors_by_edge(facet_it, window_size, neighbors);
 
-        double total_sdf_value = 0.0, total_weight = 0.0;
+    std::vector<double> smoothed_sdf(mesh.size_of_facets());
-        double current_sdf_value = get(sdf_values, facet_it);
+    Facet_vector_map sdf_output(&smoothed_sdf, facet_index_map);
-        // calculate deviation for range weighting.
-        double deviation = 0.0;
-        for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
-            it != neighbors.end(); ++it) {
-          deviation += std::pow(get(sdf_values, it->first) - current_sdf_value, 2);
-        }
-        deviation = std::sqrt(deviation / neighbors.size());
-        if(deviation == 0.0) {
-          deviation = std::numeric_limits<double>::epsilon();  //this might happen
-        }
-        for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin();
-            it != neighbors.end(); ++it) {
-          double spatial_weight = gaussian_function(it->second,
-                                  window_size / 2.0); // window_size => 2*sigma
-          double domain_weight = gaussian_function(get(sdf_values,
-                                 it->first) -  current_sdf_value, 1.5 * deviation);
-          double weight = spatial_weight * domain_weight;
 
-          total_sdf_value += get(sdf_values, it->first) * weight;
+    Filter()(mesh, get_window_size(), sdf_input, sdf_output);
-          total_weight += weight;
+    sdf_values = smoothed_sdf;
-        }
-        get(smoothed_sdf_values, facet_it) = total_sdf_value / total_weight;
-      }
-      sdf_values = smoothed_sdf_values;
-    }
   }
 
   /**
@@ -628,66 +514,6 @@ protected:
     return static_cast<int>(facet_sqrt) + 1;
   }
 
-  /**
-   * Breadth-first search on facets, by treating facets, which share a common edge, are 1-level neighbors.
-   * @param facet root facet
-   * @param max_level maximum distance (number of edges) between root facet and visited facet
-   * @param[out] neighbors visited facets and their distances to root facet
-   */
-  void get_neighbors_by_edge(Facet_const_handle& facet, int max_level,
-                             std::map<Facet_const_handle, int>& neighbors) {
-    typedef std::pair<Facet_const_handle, int> facet_level_pair;
-    std::queue<facet_level_pair> facet_queue;
-    facet_queue.push(facet_level_pair(facet, 0));
-    while(!facet_queue.empty()) {
-      const facet_level_pair& pair = facet_queue.front();
-      bool inserted = neighbors.insert(pair).second;
-      if(inserted && pair.second < max_level) {
-        typename Facet::Halfedge_around_facet_const_circulator facet_circulator =
-          pair.first->facet_begin();
-        do {
-          if(!facet_circulator->opposite()->is_border()) {
-            facet_queue.push(facet_level_pair(facet_circulator->opposite()->facet(),
-                                              pair.second + 1));
-          }
-        } while(++facet_circulator != pair.first->facet_begin());
-      }
-      facet_queue.pop();
-    }
-  }
-
-  /**
-   * Breadth-first search on facets, by treating facets, which share a common vertex, are 1-level neighbors.
-   * @param facet root facet
-   * @param max_level maximum distance (number of edges) between root facet and visited facet
-   * @param[out] neighbors visited facets and their distances to root facet
-   */
-  void get_neighbors_by_vertex(Facet_const_handle& facet,
-                               std::map<Facet_const_handle, int>& neighbors, int max_level) {
-    typedef std::pair<Facet_const_handle, int> facet_level_pair;
-    std::queue<facet_level_pair> facet_queue;
-    facet_queue.push(facet_level_pair(facet, 0));
-    while(!facet_queue.empty()) {
-      const facet_level_pair& pair = facet_queue.front();
-      bool inserted = neighbors.insert(pair).second;
-      if(inserted && pair.second < max_level) {
-        Facet_const_handle facet = pair.first;
-        Halfedge_const_iterator edge = facet->halfedge();
-        do { // loop on three vertices of the facet
-          Vertex_const_iterator vertex = edge->vertex();
-          typename Facet::Halfedge_around_vertex_const_circulator vertex_circulator =
-            vertex->vertex_begin();
-          do { // for each vertex loop on incoming edges (through those edges loop on neighbor facets which includes the vertex)
-            if(!vertex_circulator->is_border()) {
-              facet_queue.push(facet_level_pair(vertex_circulator->facet(), pair.second + 1));
-            }
-          } while(++vertex_circulator != vertex->vertex_begin());
-        } while((edge = edge->next()) != facet->halfedge());
-      }
-      facet_queue.pop();
-    }
-  }
-
   /**
    * Finds facets which have zero sdf values.
    * Sdf values on these facets are assigned to average sdf value of its neighbors.
@@ -832,12 +658,7 @@ protected:
     return data[ facet_index_map[facet] ];
   }
 
-  /**
+
-   * Gauissian function for weighting.
-   */
-  double gaussian_function(double value, double deviation) {
-    return exp(-0.5 * (std::pow(value / deviation, 2)));
-  }
 
 public:
   /**
@@ -1006,6 +827,7 @@ public:
 #endif
   }
 };
+}//namespace internal
 } //namespace CGAL
 
 #undef CGAL_NORMALIZATION_ALPHA

Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h

@@ -4,12 +4,12 @@
  * @file mesh_segmentation.h
  * The API which contains free template functions for SDF computation and mesh segmentation.
  */
-#include <CGAL/Surface_mesh_segmentation.h>
+#include <CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h>
 
 #define CGAL_DEFAULT_CONE_ANGLE (2.0 / 3.0) * CGAL_PI /**< Default opening angle for cone */
 #define CGAL_DEFAULT_NUMBER_OF_RAYS 25                /**< Default number of rays picked from cone for each facet */
 #define CGAL_DEFAULT_NUMBER_OF_LEVELS 5               /**< Default number of clusters for soft clustering */
-#define CGAL_DEFAULT_SMOOTHING_LAMBDA 0.46            /**< Default factor which indicates importance of surface features in energy minization*/
+#define CGAL_DEFAULT_SMOOTHING_LAMBDA 0.23            /**< Default factor which indicates importance of surface features in energy minization */
 
 /** CGAL */
 namespace CGAL
@@ -22,7 +22,8 @@ namespace CGAL
  * If still there is any facet which has no SDF value, minimum SDF value greater than zero is assigned to it.
  *  - Smoothed with bilateral filtering.
  *  - Linearly normalized between [0-1].
- * @param polyhedron `CGAL Polyhedron` on which SDF is computed
+ *
+ * @param polyhedron `CGAL Polyhedron` on which SDF values are computed
  * @param[out] sdf_values `WritablePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type
  * @param cone_angle opening angle for cone, expressed in radians
  * @param number_of_rays number of rays picked from cone for each facet