Added test files, and code review.

.gitattributes

@@ -4109,6 +4109,13 @@ 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/package_info/Surface_mesh_segmentation/maintainer -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/Disk_samplers_test.cpp -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/Expectation_maximization_test.cpp -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/Filters_test.cpp -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/K_means_clustering_test.cpp -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/Utils.h -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/data/cactus.off -text
+Surface_mesh_segmentation/test/Surface_mesh_segmentation/mesh_segmentation_test.cpp -text
 Surface_mesh_simplification/doc_tex/OLD[!!-~]Surface_mesh_simplification.tex -text
 Surface_mesh_simplification/doc_tex/Surface_mesh_simplification/fig/Illustration-Simplification-ALL.jpg -text
 Surface_mesh_simplification/doc_tex/Surface_mesh_simplification/fig/Illustration-Simplification-ALL.pdf -text svneol=unset#application/pdf

Surface_mesh_segmentation/doc/Concepts/SegmentationGeomTraits.h

@@ -17,7 +17,7 @@ public:
 /// \name Number type
 /// @{
 /*! 
-A number type model of FieldWithSqrt (double or float is recommanded)
+A number type model of FieldWithSqrt (double or float is recommended)
 */ 
 typedef Hidden_type NT;
 /// @}

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

@@ -5,6 +5,7 @@
  * @brief This file contains 3 sampling methods, which can be used as a template parameter for CGAL::internal::SDF_calculation.
  */
 #include <cmath>
+#include <CGAL/number_type_basic.h>
 
 #define CGAL_ANGLE_ST_DEV_DIVIDER 3.0
 
@@ -40,10 +41,10 @@ namespace internal
 /**
  * @brief Uses Vogel's method to sample points from unit-disk.
  *
- * Template parameter @a Tuple should have a constructor which takes 3 double parameters.
+ * @tparam Tuple should have a constructor which takes 3 double parameters.
  * @see Disk_samplers.h, SDF_calculation
  */
-template<class Tuple>
+template<class Tuple, bool uniform = false>
 class Vogel_disk_sampling
 {
 public:
@@ -60,8 +61,7 @@ public:
   template<class OutputIterator>
   void operator()(int number_of_points,
                   double cone_angle,
-                  OutputIterator out_it,
-                  bool uniform = false) const {
+                  OutputIterator out_it) const {
     const double golden_ratio = 3.0 - std::sqrt(5.0);
 
     if(uniform) {
@@ -116,7 +116,7 @@ public:
 /**
  * @brief Uses polar mapping to sample points from unit-disk.
  *
- * Template parameter @a Tuple should have a constructor which takes 3 double parameters.
+ * @tparam Tuple should have a constructor which takes 3 double parameters.
  */
 template<class Tuple>
 class Polar_disk_sampling
@@ -183,7 +183,7 @@ public:
 /**
  * @brief Uses concentric mapping to sample points from unit-disk.
  *
- * Template parameter @a Tuple should have a constructor which takes 3 double parameters.
+ * @tparam Tuple should have a constructor which takes 3 double parameters.
  */
 template<class Tuple>
 class Concentric_disk_sampling
@@ -208,7 +208,8 @@ public:
     const int number_of_points_sqrt = static_cast<int>(std::sqrt(
                                         static_cast<double>(number_of_points)));
     const double length_of_normal = 1.0 / tan(cone_angle / 2.0);
-    const double fraction = 2.0 / (number_of_points_sqrt -1);
+    const double fraction = (number_of_points_sqrt == 1) ? 0.0
+                            : 2.0 / (number_of_points_sqrt -1);
     // use cone_angle / 3 as one standard deviation while weighting.
     const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER;
 

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

@@ -7,10 +7,12 @@
 #include <limits>
 
 #include <CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h>
+#include <CGAL/assertions.h>
 
-#define CGAL_DEFAULT_MAXIMUM_ITERATION 15
+#define CGAL_DEFAULT_MAXIMUM_ITERATION 10
+#define CGAL_DEFAULT_NUMBER_OF_RUN 15
 #define CGAL_DEFAULT_THRESHOLD 1e-3
-#define CGAL_DEFAULT_NUMBER_OF_RUN 20
+
 #define CGAL_DEFAULT_SEED 1340818006
 
 namespace CGAL

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

@@ -142,6 +142,8 @@ public:
     data_centers.reserve(points.size());
     for(std::vector<K_means_point>::iterator point_it = points.begin();
         point_it != points.end(); ++point_it) {
+      point_it->calculate_new_center(
+        centers); // just refind closest center (incase order of centers are changed etc)
       data_centers.push_back(point_it->center_id);
     }
   }

Surface_mesh_segmentation/test/Surface_mesh_segmentation/Disk_samplers_test.cpp

@@ -0,0 +1,63 @@
+#include <CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h>
+
+#include <boost/tuple/tuple.hpp>
+#include <vector>
+
+typedef boost::tuple<double, double, double> boost_tuple;
+
+void print(const std::vector<boost_tuple>& samples)
+{
+	const int map_size = 31;
+    const int map_size_2 = 45;
+	std::vector<std::vector<bool> > sample_map(map_size, std::vector<bool>(map_size_2, false));
+	
+	for(std::vector<boost_tuple>::const_iterator sample_it = samples.begin(); 
+		sample_it != samples.end(); ++sample_it)
+	{
+		double x = (sample_it->get<0>() +1)/2;
+		double y = (sample_it->get<1>() +1)/2;
+		x *= (map_size-1);
+		y *= (map_size_2-1);
+		int x_c  = static_cast<int>(x + 0.49);
+		int y_c  = static_cast<int>(y + 0.49);
+		sample_map[x_c][y_c] = true;
+	}
+	for(int i = 0; i < map_size; ++i)
+	{
+		for(int j = 0; j < map_size_2; ++j)
+		{
+			if(sample_map[i][j]){ std::cout << "*"; }
+			else                { std::cout << " "; }
+		}
+		std::cout << std::endl;
+	}
+	std::cout << std::endl;
+}
+/**
+ * Uses disk sampling functors to sample points from unit-disk.
+ * It also prints sampled points for visual debugging.
+ *
+ * Note that it always return EXIT_SUCCESS
+ */
+int main(void)
+{
+	CGAL::internal::Vogel_disk_sampling<boost_tuple> sampling_1;
+	CGAL::internal::Vogel_disk_sampling<boost_tuple, true> sampling_2;
+	CGAL::internal::Polar_disk_sampling<boost_tuple> sampling_3;
+	CGAL::internal::Concentric_disk_sampling<boost_tuple> sampling_4;
+
+	std::vector<boost_tuple> samples_1;
+	std::vector<boost_tuple> samples_2;
+	std::vector<boost_tuple> samples_3;
+    std::vector<boost_tuple> samples_4;
+    
+	sampling_1(64, 2.0 / 3.0 * CGAL_PI, std::back_inserter(samples_1));   
+	sampling_2(64, 2.0 / 3.0 * CGAL_PI, std::back_inserter(samples_2));  
+	sampling_3(64, 2.0 / 3.0 * CGAL_PI, std::back_inserter(samples_3));
+	sampling_4(64, 2.0 / 3.0 * CGAL_PI, std::back_inserter(samples_4));
+
+    print(samples_1);
+	print(samples_2);
+	print(samples_3);
+	print(samples_4);	
+}	

Surface_mesh_segmentation/test/Surface_mesh_segmentation/Expectation_maximization_test.cpp

@@ -0,0 +1,84 @@
+
+#include <CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h>
+
+#include <boost/random.hpp>
+#include <boost/random/normal_distribution.hpp>
+/**
+ * Generates sample points using a few gauissians.
+ * Then applies gmm fitting on these generated points.
+ * Provides a heuristic score for each gmm fitting result.
+ *
+ * Note that it always return EXIT_SUCCESS
+ */
+int main(void)
+{
+    boost::mt19937 engine;
+	engine.seed(1340818006);
+
+    // generate random data using gauissians below
+    std::vector< boost::normal_distribution<double> > distributions;
+    distributions.push_back(boost::normal_distribution<double>(0.1, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(0.4, 0.1));
+    distributions.push_back(boost::normal_distribution<double>(0.55, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(0.7, 0.1));
+    distributions.push_back(boost::normal_distribution<double>(0.9, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(1.0, 0.05));
+    
+    std::vector<double> data;    
+    for(std::vector< boost::normal_distribution<double> >::iterator it = distributions.begin();
+      it != distributions.end(); ++it)
+    {
+        boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> > var_nor(engine, *it);
+        
+        for(int i = 0; i < 300; ++i) { data.push_back(var_nor()); }
+    }
+    
+    // calculate closest center (using above gauissians) for each generated points
+    // we will compare it with gmm fitting results
+	// also we might want to compute mixing coef for each center and select centers according to mixing_coef * prob(data)
+    std::vector<int> data_centers;
+    for(std::vector<double>::iterator it = data.begin(); it != data.end(); ++it)
+    {
+        int center_id = -1, center_counter = 0;;
+        double min_distance = (std::numeric_limits<double>::max)();        
+        for(std::vector< boost::normal_distribution<double> >::iterator dis_it = distributions.begin();
+          dis_it != distributions.end(); ++dis_it, ++center_counter)
+        {
+            double distance = std::abs(*it - dis_it->mean());            
+            if(min_distance > distance) 
+            {
+                min_distance = distance;
+                center_id = center_counter;
+            }            
+        }
+        data_centers.push_back(center_id);
+    }
+    
+    // apply gmm fitting clustering
+    typedef CGAL::internal::Expectation_maximization E_M;
+    std::vector<E_M> gmm_fitters;
+    gmm_fitters.push_back(E_M(distributions.size(), data, E_M::PLUS_INITIALIZATION));
+    gmm_fitters.push_back(E_M(distributions.size(), data, E_M::RANDOM_INITIALIZATION));
+    gmm_fitters.push_back(E_M(distributions.size(), data, E_M::K_MEANS_INITIALIZATION));
+
+    std::vector< std::vector<int> > calculated_centers(gmm_fitters.size());
+    std::vector< std::vector<int> >::iterator calc_centers_it = calculated_centers.begin();
+    for(std::vector<E_M>::iterator it = gmm_fitters.begin(); it != gmm_fitters.end(); ++it, ++calc_centers_it)
+    {
+        it->fill_with_center_ids(*calc_centers_it);
+    }
+    
+    std::cout << "Compare results of EM with 'expected' (but be aware, it is not optimal result in terms of likelihood)" << std::endl;
+    std::cout << "Another words a clustering which has higher likelihood can result in worse score in here" << std::endl;
+    for(std::vector< std::vector<int> >::iterator calc_centers_it = calculated_centers.begin();
+        calc_centers_it != calculated_centers.end(); ++calc_centers_it)
+    {
+        int true_count = 0;
+        std::vector<int>::iterator calculated_it = calc_centers_it->begin();
+        for(std::vector<int>::iterator it = data_centers.begin(); it != data_centers.end(); ++it, ++calculated_it)
+        {
+            if( (*it) == (*calculated_it) ) { ++true_count; }
+        }
+        std::cout << "[0,1]: " << static_cast<double>(true_count) / data_centers.size() << std::endl;
+    }
+}	

Surface_mesh_segmentation/test/Surface_mesh_segmentation/Filters_test.cpp

@@ -0,0 +1,69 @@
+#include <CGAL/internal/Surface_mesh_segmentation/Filters.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+
+#include <cmath>
+#include <map>
+#include <boost/property_map/property_map.hpp>
+
+#include "Utils.h"
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+
+/**
+ * Uses bilateral and median filtering to smooth associated values with facets.
+ * It also prints average differences between them (smoothed values with two methods).
+ * Note that it always return EXIT_SUCCESS if .off file is read successfully.
+ */
+int main(void)
+{	
+	Polyhedron mesh;
+	if( !read_to_polyhedron("./data/cactus.off", mesh) ) { return 1; }
+		
+	typedef std::map< typename Polyhedron::Facet_const_handle, double> Facet_double_map;
+    Facet_double_map internal_1;
+    boost::associative_property_map<Facet_double_map> value_pmap_1(internal_1);
+	
+	Facet_double_map internal_2;
+    boost::associative_property_map<Facet_double_map> value_pmap_2(internal_2);
+	
+	double counter = 0.0;
+	for(Polyhedron::Facet_const_iterator facet_it = mesh.facets_begin(); 
+        facet_it != mesh.facets_end(); ++facet_it, ++counter)   
+    {
+        value_pmap_1[facet_it] = value_pmap_2[facet_it] = counter;                               
+    }
+	const double average_value = counter / 2.0;
+	
+	CGAL::internal::Bilateral_filtering<Polyhedron, CGAL::internal::Neighbor_selector_by_edge<Polyhedron> > filter_1;
+	CGAL::internal::Bilateral_filtering<Polyhedron, CGAL::internal::Neighbor_selector_by_vertex<Polyhedron> > filter_2;
+	CGAL::internal::Median_filtering<Polyhedron, CGAL::internal::Neighbor_selector_by_edge<Polyhedron> > filter_3;
+	CGAL::internal::Median_filtering<Polyhedron, CGAL::internal::Neighbor_selector_by_vertex<Polyhedron> > filter_4;
+	
+	filter_1(mesh, 2, value_pmap_1);
+	filter_3(mesh, 2, value_pmap_2);
+	
+	double average_dif = 0.0;
+	for(Polyhedron::Facet_const_iterator facet_it = mesh.facets_begin(); 
+        facet_it != mesh.facets_end(); ++facet_it)   
+    {
+		average_dif += std::abs(value_pmap_1[facet_it] - value_pmap_2[facet_it]);
+	}		
+	average_dif = (average_dif / mesh.size_of_facets()) / average_value;
+	std::cout << "average differences between bilateral and median filters: " << average_dif << std::endl;
+	
+	filter_2(mesh, 2, value_pmap_1);
+	filter_4(mesh, 2, value_pmap_2);
+	
+	average_dif = 0.0;
+	for(Polyhedron::Facet_const_iterator facet_it = mesh.facets_begin(); 
+        facet_it != mesh.facets_end(); ++facet_it)   
+    {
+		average_dif += std::abs(value_pmap_1[facet_it] - value_pmap_2[facet_it]);
+	}		
+	average_dif = (average_dif / mesh.size_of_facets()) / average_value;
+	std::cout << "average differences between bilateral and median filters (2) :" << average_dif << std::endl;	
+}

Surface_mesh_segmentation/test/Surface_mesh_segmentation/K_means_clustering_test.cpp

@@ -0,0 +1,84 @@
+
+#include <CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h>
+
+#include <boost/random.hpp>
+#include <boost/random/normal_distribution.hpp>
+/**
+ * Generates sample points using a few gauissians.
+ * Then applies k-means on these generated points.
+ * Provides a heuristic score for each k-means clustering result.
+ *
+ * Note that it always return EXIT_SUCCESS
+ */
+int main(void)
+{
+    boost::mt19937 engine;
+	engine.seed(1340818006);
+
+    // generate random data using gauissians below
+    std::vector< boost::normal_distribution<double> > distributions;
+    distributions.push_back(boost::normal_distribution<double>(0.1, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(0.4, 0.1));
+    distributions.push_back(boost::normal_distribution<double>(0.55, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(0.7, 0.1));
+    distributions.push_back(boost::normal_distribution<double>(0.9, 0.05));
+    distributions.push_back(boost::normal_distribution<double>(1.0, 0.05));
+    
+    std::vector<double> data;    
+    for(std::vector< boost::normal_distribution<double> >::iterator it = distributions.begin();
+      it != distributions.end(); ++it)
+    {
+        boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> > var_nor(engine, *it);
+        
+        for(int i = 0; i < 300; ++i) { data.push_back(var_nor()); }
+    }
+    
+    // calculate closest center (using above gauissians) for each generated points
+    // we will compare it with k-means results
+    std::vector<int> data_centers;
+    for(std::vector<double>::iterator it = data.begin(); it != data.end(); ++it)
+    {
+        int center_id = -1, center_counter = 0;;
+        double min_distance = (std::numeric_limits<double>::max)();        
+        for(std::vector< boost::normal_distribution<double> >::iterator dis_it = distributions.begin();
+          dis_it != distributions.end(); ++dis_it, ++center_counter)
+        {
+            double distance = std::abs(*it - dis_it->mean());            
+            if(min_distance > distance) 
+            {
+                min_distance = distance;
+                center_id = center_counter;
+            }            
+        }
+        data_centers.push_back(center_id);
+    }
+    
+    // apply k-means clustering
+    typedef CGAL::internal::K_means_clustering K_means;
+    std::vector<K_means> k_means;
+    k_means.push_back(K_means(distributions.size(), data, K_means::PLUS_INITIALIZATION));
+    k_means.push_back(K_means(distributions.size(), data, K_means::PLUS_INITIALIZATION, 2, 5));
+    k_means.push_back(K_means(distributions.size(), data, K_means::RANDOM_INITIALIZATION));
+    k_means.push_back(K_means(distributions.size(), data, K_means::RANDOM_INITIALIZATION, 2, 5));
+
+    std::vector< std::vector<int> > calculated_centers(k_means.size());
+    std::vector< std::vector<int> >::iterator calc_centers_it = calculated_centers.begin();
+    for(std::vector<K_means>::iterator it = k_means.begin(); it != k_means.end(); ++it, ++calc_centers_it)
+    {
+        it->fill_with_center_ids(*calc_centers_it);
+    }
+    
+    std::cout << "Compare results of k-means with 'expected' (but be aware, it is not optimal result in terms of within-cluster error)" << std::endl;
+    std::cout << "Another words a clustering which has smaller within-cluster error can result in worse score in here" << std::endl;
+    for(std::vector< std::vector<int> >::iterator calc_centers_it = calculated_centers.begin();
+        calc_centers_it != calculated_centers.end(); ++calc_centers_it)
+    {
+        int true_count = 0;
+        std::vector<int>::iterator calculated_it = calc_centers_it->begin();
+        for(std::vector<int>::iterator it = data_centers.begin(); it != data_centers.end(); ++it, ++calculated_it)
+        {
+            if( (*it) == (*calculated_it) ) { ++true_count; }
+        }
+        std::cout << "[0,1]: " << static_cast<double>(true_count) / data_centers.size() << std::endl;
+    }
+}	

Surface_mesh_segmentation/test/Surface_mesh_segmentation/Utils.h

@@ -0,0 +1,14 @@
+#include <fstream>
+#include <iostream>
+
+template<class Polyhedron>
+bool read_to_polyhedron(const char* file_name, Polyhedron& mesh)
+{
+	std::ifstream input(file_name);
+    
+    if ( !input || !(input >> mesh) || mesh.empty() ){
+		std::cerr << "Problem occured while reading off file";
+		return false;
+    }
+	return true;
+}

Surface_mesh_segmentation/test/Surface_mesh_segmentation/mesh_segmentation_test.cpp

@@ -0,0 +1,51 @@
+#include <CGAL/mesh_segmentation.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+
+#include <cmath>
+#include <map>
+#include <boost/property_map/property_map.hpp>
+
+#include "Utils.h"
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+
+/**
+ * Note that it always return EXIT_SUCCESS if .off file is read successfully.
+ */
+int main(void)
+{	
+	Polyhedron mesh;
+	if( !read_to_polyhedron("./data/cactus.off", mesh) ) { return 1; }
+	
+	typedef std::map<Polyhedron::Facet_const_handle, double> Facet_double_map;
+    Facet_double_map internal_map;
+    boost::associative_property_map<Facet_double_map> sdf_property_map(internal_map);	
+	
+	std::pair<double, double> min_max_sdf = CGAL::sdf_values_computation(mesh, sdf_property_map);
+	std::cout << "minimum sdf: " << min_max_sdf.first << " maximum sdf: " << min_max_sdf.second << std::endl;
+	
+	
+	typedef std::map<Polyhedron::Facet_const_handle, int> Facet_int_map;
+		Facet_int_map internal_segment_map;
+    boost::associative_property_map<Facet_int_map> segment_property_map(internal_segment_map);
+	
+    int nb_segments = CGAL::surface_mesh_segmentation_from_sdf_values(
+		mesh, sdf_property_map, segment_property_map);
+		
+	if(nb_segments != 3)
+	{
+		std::cout << "Number of segments should be 3 for cactus model (since it is pretty easy model to segment)" << std::endl;
+	}
+
+	int nb_segments_2 = CGAL::surface_mesh_segmentation(mesh, segment_property_map);
+	
+	if(nb_segments_2 != nb_segments) 
+	{
+		std::cout << "Inconsistency between 'surface_mesh_segmentation' and 'surface_mesh_segmentation_from_sdf_values'" 
+			<< std::endl;
+	}
+}