diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h index 0d947a4f75e..da4230153f4 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h @@ -6,7 +6,6 @@ * (https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/BGL) which introduces some features * for face-based traversal / manipulation by FaceGraphs * +) Deciding on which parameters will be taken from user - * +) Make it more readable: calculate_sdf_value_of_facet function. * * About paper (and correctness / efficiency etc.): * +) Weighting ray distances with inverse of their angles: not sure how to weight exactly @@ -28,6 +27,8 @@ #include #include +#include + #include #include #include @@ -70,11 +71,13 @@ protected: Tree; typedef typename Tree::Object_and_primitive_id Object_and_primitive_id; + typedef typename Tree::Primitive_id + Primitive_id; typedef std::map Face_value_map; typedef std::map Face_center_map; /*Sampled points from disk, t1 = coordinate-x, t2 = coordinate-y, t3 = angle with cone-normal (weight). */ - typedef CGAL::Triple Disk_sample; + typedef CGAL::Triple Disk_sample; typedef std::vector > Disk_samples_list; template @@ -111,12 +114,10 @@ public: double calculate_sdf_value_of_facet (const Facet_handle& facet, const Tree& tree) const; - void cast_and_return_minimum (const Ray& ray, const Tree& tree, - const Facet_handle& facet, - bool& is_found, double& min_distance) const; - void cast_and_return_minimum_use_closest (const Ray& ray, const Tree& tree, - const Facet_handle& facet, - bool& is_found, double& min_distance) const; + boost::optional cast_and_return_minimum(const Ray& ray, + const Tree& tree, const Facet_handle& facet) const; + boost::optional cast_and_return_minimum_use_closest(const Ray& ray, + const Tree& tree, const Facet_handle& facet) const; double calculate_sdf_value_from_rays (std::vector& ray_distances, std::vector& ray_weights) const; @@ -157,8 +158,8 @@ inline Surface_mesh_segmentation::Surface_mesh_segmentation( calculate_sdf_values(); SEG_DEBUG(std::cout << t) apply_GMM_fitting(); - //write_sdf_values("sdf_values_sample_dino_2.txt"); - //read_sdf_values("sdf_values_sample_dino.txt"); + //write_sdf_values("sdf_values_sample_dino_ws.txt"); + //read_sdf_values("sdf_values_sample_elephant.txt"); } template @@ -211,55 +212,57 @@ Surface_mesh_segmentation::calculate_sdf_value_of_facet( sample_it != disk_samples.end(); ++sample_it) { Vector disk_vector = v1 * sample_it->first + v2 * sample_it->second; Ray ray(center, normal + disk_vector); - double min_distance; - bool is_found; - cast_and_return_minimum_use_closest(ray, tree, facet, is_found, min_distance); - if(!is_found) { + + //boost::optional min_distance = cast_and_return_minimum(ray, tree, facet); + boost::optional min_distance = cast_and_return_minimum_use_closest(ray, + tree, facet); + if(!min_distance) { continue; } ray_weights.push_back(sample_it->third); - ray_distances.push_back(min_distance); + ray_distances.push_back(*min_distance); } return calculate_sdf_value_from_rays_with_trimmed_mean(ray_distances, ray_weights); } template -void Surface_mesh_segmentation::cast_and_return_minimum( - const Ray& ray, const Tree& tree, const Facet_handle& facet, - bool& is_found, double& min_distance) const +boost::optional +Surface_mesh_segmentation::cast_and_return_minimum( + const Ray& ray, const Tree& tree, const Facet_handle& facet) const { + boost::optional min_distance; std::list intersections; - tree.all_intersection(ray, std::back_inserter(intersections)); + tree.all_intersections(ray, std::back_inserter(intersections)); Vector min_i_ray; - typename Tree::Primitive_id min_id; - is_found = false; + Primitive_id min_id; for(typename std::list::iterator op_it = intersections.begin(); op_it != intersections.end() ; ++op_it) { - CGAL::Object object = op_it->first; - typename Tree::Primitive_id id = op_it->second; - Point i_point; - if(id == facet) { + CGAL::Object object = op_it->first; + Primitive_id id = op_it->second; + if(id == facet) { continue; //Since center is located on related facet, we should skip it if there is an intersection with it. } + + Point i_point; if(!CGAL::assign(i_point, object)) { continue; //What to do here (in case of intersection object is a segment), I am not sure ??? } + Vector i_ray = (ray.source() - i_point); double new_distance = i_ray.squared_length(); - if(!is_found || new_distance < min_distance) { + if(!min_distance || new_distance < min_distance) { min_distance = new_distance; min_id = id; min_i_ray = i_ray; - is_found = true; } } - if(!is_found) { - return; + if(!min_distance) { + return min_distance; } - min_distance = sqrt(min_distance); + Point min_v1 = min_id->halfedge()->vertex()->point(); Point min_v2 = min_id->halfedge()->next()->vertex()->point(); Point min_v3 = min_id->halfedge()->next()->next()->vertex()->point(); @@ -267,43 +270,52 @@ void Surface_mesh_segmentation::cast_and_return_minimum( if(CGAL::angle(CGAL::ORIGIN + min_i_ray, Point(CGAL::ORIGIN), CGAL::ORIGIN + min_normal) != CGAL::ACUTE) { - is_found = false; + return boost::optional(); } - //total_intersection += intersections.size(); - //total_cast++; + return (min_distance = sqrt(*min_distance)); } template -void Surface_mesh_segmentation::cast_and_return_minimum_use_closest ( +boost::optional +Surface_mesh_segmentation::cast_and_return_minimum_use_closest ( const Ray& ray, const Tree& tree, - const Facet_handle& facet, bool& is_found, double& min_distance) const + const Facet_handle& facet) const { - is_found = false; + //static double dist = 0.1; + //boost::optional min_distance_2 = dist++; + //return min_distance_2; + + boost::optional min_distance; boost::optional intersection = tree.closest_intersection(ray); + //boost::optional intersection = tree.any_intersection(ray); if(!intersection) { - return; + return min_distance; } CGAL::Object object = intersection->first; - typename Tree::Primitive_id min_id = intersection->second; + Primitive_id min_id = intersection->second; + + if(min_id == facet) { + CGAL_error(); // There should be no such case, after center-facet arrangments. + } + Point i_point; if(!CGAL::assign(i_point, object)) { - return; + return min_distance; } Vector min_i_ray = ray.source() - i_point; - min_distance = sqrt(min_i_ray.squared_length()); - Point min_v1 = min_id->halfedge()->vertex()->point(); Point min_v2 = min_id->halfedge()->next()->vertex()->point(); Point min_v3 = min_id->halfedge()->next()->next()->vertex()->point(); Vector min_normal = CGAL::normal(min_v1, min_v2, min_v3) * -1.0; if(CGAL::angle(CGAL::ORIGIN + min_i_ray, Point(CGAL::ORIGIN), - CGAL::ORIGIN + min_normal) == CGAL::ACUTE) { - is_found = true; + CGAL::ORIGIN + min_normal) != CGAL::ACUTE) { + return min_distance; } + return (min_distance = sqrt(min_i_ray.squared_length())); } template @@ -629,7 +641,7 @@ inline void Surface_mesh_segmentation::apply_GMM_fitting() sdf_vector.push_back(pair_it->second); } SEG_DEBUG(Timer t) - internal::Expectation_maximization fitter(number_of_centers, sdf_vector); + internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 50); SEG_DEBUG(std::cout << t) std::vector center_memberships; fitter.fill_with_center_ids(center_memberships); diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h index 34b1ff63dc9..b7ac4fc782f 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h @@ -2,18 +2,26 @@ #define CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H /* NEED TO BE DONE */ /* About implementation: - * Calculating probability multiple times, solution: storing matrix(cluster, point) for probability values. - * About safe division: where centers have no members, in graph cut it may happen + * + There are two implementations : one of them is using probability matrix (center x point), the other is not. + * + Also code size is twice large, I did not try to simplify before deciding on whether we are going to make + * both approaches (with matrix and without) available or not. + * + There are a lot of parameters (with default values) and initialization types, + * so I am planning to use whether 'Named Parameter Idiom' or Boost Parameter Library... + * */ + #include #include #include #include #include #include - #include #include + +#define DEF_MAX_ITER 100 +#define DEF_THRESHOLD 1e-4 +#define USE_MATRIX true //#define ACTIVATE_SEGMENTATION_EM_DEBUG #ifdef ACTIVATE_SEGMENTATION_EM_DEBUG #define SEG_DEBUG(x) x; @@ -23,6 +31,8 @@ namespace CGAL { +namespace internal +{ class Gaussian_point; @@ -46,6 +56,9 @@ public: double e_over = -0.5 * pow((x - mean) / deviation, 2); return exp(e_over); } + double probability_with_coef(double x) const { + return probability(x) * mixing_coefficient; + } bool operator < (const Gaussian_center& center) const { return mean < center.mean; } @@ -60,12 +73,13 @@ public: Gaussian_point(double data): data(data), total_membership(0.0) { } - void calculate_total_membership(const std::vector& centers) { + double calculate_total_membership(const std::vector& centers) { total_membership = 0.0; for(std::vector::const_iterator it = centers.begin(); it != centers.end(); ++it) { total_membership += it->probability(data) * it->mixing_coefficient; } + return total_membership; } }; @@ -106,22 +120,52 @@ public: int maximum_iteration; bool is_converged; protected: - int seed; + unsigned int seed; + std::vector > probability_matrix; public: + /* For uniform initialization, with one run */ Expectation_maximization(int number_of_centers, const std::vector& data, - const std::vector& initial_centers = std::vector(), - int number_of_runs = 50, int maximum_iteration = 100) - : points(data.begin(), data.end()), threshold(1e-4), - maximum_iteration(maximum_iteration), is_converged(false), seed(time(NULL)) { - srand(seed); - if(initial_centers.empty()) { - calculate_fitting_with_multiple_run(number_of_centers, number_of_runs); - } else { - initiate_centers(number_of_centers, initial_centers); - calculate_fitting(); + double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER, + bool use_matrix = USE_MATRIX) + : points(data.begin(), data.end()), threshold(threshold), + maximum_iteration(maximum_iteration), is_converged(false), seed(0), + probability_matrix( use_matrix ? std::vector > + (number_of_centers, std::vector(points.size())) : + std::vector >()) { + for(int i = 0; i < 20; ++i) { + initiate_centers_uniformly(number_of_centers); + use_matrix ? calculate_fitting_with_matrix() : calculate_fitting(); } } + /* For initialization with provided center ids per point, with one run */ + Expectation_maximization(int number_of_centers, const std::vector& data, + const std::vector& initial_center_ids, + double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER, + bool use_matrix = USE_MATRIX) + : points(data.begin(), data.end()), threshold(threshold), + maximum_iteration(maximum_iteration), is_converged(false), seed(0), + probability_matrix( use_matrix ? std::vector > + (number_of_centers, std::vector(points.size())) : + std::vector >()) { + initiate_centers_from_memberships(number_of_centers, initial_center_ids); + use_matrix ? calculate_fitting_with_matrix() : calculate_fitting(); + } + /* For initialization with random center selection (Forgy), with multiple run */ + Expectation_maximization(int number_of_centers, const std::vector& data, + int number_of_runs, + double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER, + bool use_matrix = USE_MATRIX) + : points(data.begin(), data.end()), threshold(threshold), + maximum_iteration(maximum_iteration), is_converged(false), + seed(static_cast(time(NULL))), + probability_matrix( use_matrix ? std::vector > + (number_of_centers, std::vector(points.size())) : + std::vector >()) { + srand(seed); + calculate_fitting_with_multiple_run(number_of_centers, number_of_runs, + use_matrix); + } void fill_with_center_ids(std::vector& data_centers) { data_centers.reserve(points.size()); @@ -142,6 +186,7 @@ public: } } protected: + void initiate_centers_randomly(int number_of_centers) { centers.clear(); /* Randomly generate means of centers */ @@ -149,7 +194,7 @@ protected: double initial_mixing_coefficient = 1.0; double initial_deviation = 1.0 / (2.0 * number_of_centers); for(int i = 0; i < number_of_centers; ++i) { - double random_index = rand() % points.size(); + int random_index = rand() % points.size(); //center_indexes.push_back(random_index); Gaussian_point mean_point = points[random_index]; double initial_mean = mean_point.data; @@ -162,7 +207,7 @@ protected: void initiate_centers_uniformly(int number_of_centers) { centers.clear(); - /* Uniformly generate centers */ + /* Uniformly generate means of centers */ double initial_deviation = 1.0 / (2.0 * number_of_centers); double initial_mixing_coefficient = 1.0; for(int i = 0; i < number_of_centers; ++i) { @@ -174,15 +219,14 @@ protected: } void initiate_centers_from_memberships(int number_of_centers, - const std::vector& initial_centers) { - centers.clear(); + const std::vector& initial_center_ids) { /* Calculate mean */ - int number_of_point = initial_centers.size(); + int number_of_points = initial_center_ids.size(); centers = std::vector(number_of_centers); std::vector member_count(number_of_centers, 0); - for(int i = 0; i < number_of_point; ++i) { - int center_id = initial_centers[i]; + for(int i = 0; i < number_of_points; ++i) { + int center_id = initial_center_ids[i]; double data = points[i].data; centers[center_id].mean += data; member_count[center_id] += 1; @@ -191,11 +235,11 @@ protected: for(int i = 0; i < number_of_centers; ++i) { centers[i].mean /= member_count[i]; centers[i].mixing_coefficient = member_count[i] / static_cast - (number_of_point); + (number_of_points); } /* Calculate deviation */ - for(int i = 0; i < number_of_point; ++i) { - int center_id = initial_centers[i]; + for(int i = 0; i < number_of_points; ++i) { + int center_id = initial_center_ids[i]; double data = points[i].data; centers[center_id].deviation += pow(data - centers[center_id].mean, 2); } @@ -205,29 +249,6 @@ protected: sort(centers.begin(), centers.end()); } - void initiate_centers(int number_of_centers, - const std::vector& initial_centers) { - if(initial_centers.empty()) { - initiate_centers_uniformly(number_of_centers); - } else { - initiate_centers_from_memberships(number_of_centers, initial_centers); - } - } - /*Calculates total membership values for a point */ - void calculate_membership() { - for(std::vector::iterator point_it = points.begin(); - point_it != points.end(); ++point_it) { - point_it->calculate_total_membership(centers); - } - } - double calculate_deviation(const Gaussian_point& point) const { - double deviation = 0.0; - for(std::vector::const_iterator point_it = points.begin(); - point_it != points.end(); ++point_it) { - deviation += pow(point_it->data - point.data, 2); - } - return sqrt(deviation / points.size()); - } /*Calculates new parameter values for each cluster */ void calculate_parameters() { for(std::vector::iterator center_it = centers.begin(); @@ -235,35 +256,109 @@ protected: center_it->calculate_parameters(points); } } + + void calculate_parameters_with_matrix() { + for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { + // Calculate new mean + double new_mean = 0.0, total_membership = 0.0; + for(std::size_t point_i = 0; point_i < points.size(); ++point_i) { + double membership = probability_matrix[center_i][point_i]; + new_mean += membership * points[point_i].data; + total_membership += membership; + } + new_mean /= total_membership; + // Calculate new deviation + double new_deviation = 0.0; + for(std::size_t point_i = 0; point_i < points.size(); ++point_i) { + double membership = probability_matrix[center_i][point_i]; + new_deviation += membership * pow(points[point_i].data - new_mean, 2); + } + new_deviation = sqrt(new_deviation/total_membership); + // Assign new parameters + centers[center_i].mixing_coefficient = total_membership; + centers[center_i].deviation = new_deviation; + centers[center_i].mean = new_mean; + } + } + /*Calculates how much this adjustment is likely to represent data*/ double calculate_likelihood() { double likelihood = 0.0; for(std::vector::iterator point_it = points.begin(); point_it != points.end(); ++point_it) { - double point_likelihood = 0.0; - for(std::vector::iterator center_it = centers.begin(); - center_it != centers.end(); ++center_it) { - point_likelihood += center_it->mixing_coefficient * center_it->probability( - point_it->data); - } + double point_likelihood = point_it->calculate_total_membership(centers); likelihood += log(point_likelihood); } return likelihood; } - double iterate() { - calculate_membership(); + double calculate_likelihood_with_matrix() { + /** + * Calculate Log-likelihood + * The trick (merely a trick) is while calculating log-likelihood, we also store probability results into matrix, + * so that in next iteration we do not have to calculate them again. + */ + double likelihood = 0.0; + for(std::size_t point_i = 0; point_i < points.size(); ++point_i) { + double total_membership = 0.0; + for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { + double membership = centers[center_i].probability_with_coef( + points[point_i].data); + total_membership += membership; + probability_matrix[center_i][point_i] = membership; + } + for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { + probability_matrix[center_i][point_i] /= total_membership; + } + likelihood += log(total_membership); + } + return likelihood; + } + + double iterate(bool first_iteration) { + // E-step + if(first_iteration) { + calculate_likelihood(); + } + // M-step calculate_parameters(); + // Likelihood step + // It also stores total_membership values into points, so corresponds to E-step of next iteration. return calculate_likelihood(); } + double iterate_with_matrix(bool first_iteration) { + // E-step + if(first_iteration) { + calculate_likelihood_with_matrix(); + } + // M-step + calculate_parameters_with_matrix(); + // Likelihood step + // It also stores membership values into matrix, so corresponds to E-step of next iteration. + return calculate_likelihood_with_matrix(); + } + double calculate_fitting() { double likelihood = -(std::numeric_limits::max)(), prev_likelihood; int iteration_count = 0; is_converged = false; while(!is_converged && iteration_count++ < maximum_iteration) { prev_likelihood = likelihood; - likelihood = iterate(); + likelihood = iterate(iteration_count == 1); + is_converged = likelihood - prev_likelihood < threshold * fabs(likelihood); + } + SEG_DEBUG(std::cout << likelihood << " " << iteration_count << std::endl) + return likelihood; + } + + double calculate_fitting_with_matrix() { + double likelihood = -(std::numeric_limits::max)(), prev_likelihood; + int iteration_count = 0; + is_converged = false; + while(!is_converged && iteration_count++ < maximum_iteration) { + prev_likelihood = likelihood; + likelihood = iterate_with_matrix(iteration_count == 1); is_converged = likelihood - prev_likelihood < threshold * fabs(likelihood); } SEG_DEBUG(std::cout << likelihood << " " << iteration_count << std::endl) @@ -271,15 +366,16 @@ protected: } void calculate_fitting_with_multiple_run(int number_of_centers, - int number_of_run) { + int number_of_run, bool use_matrix) { double max_likelihood = -(std::numeric_limits::max)(); std::vector max_centers; while(number_of_run-- > 0) { initiate_centers_randomly(number_of_centers); - double likelihood = calculate_fitting(); - write_center_parameters("center_paramters.txt"); + double likelihood = use_matrix ? calculate_fitting_with_matrix() : + calculate_fitting(); + //write_center_parameters("center_paramters.txt"); if(likelihood > max_likelihood) { max_centers = centers; max_likelihood = likelihood; @@ -307,7 +403,11 @@ protected: } } }; +}//namespace internal }//namespace CGAL +#undef DEF_MAX_ITER +#undef DEF_THRESHOLD +#undef USE_MATRIX #ifdef SEG_DEBUG #undef SEG_DEBUG #endif diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h index 77e3f8de740..11c89145407 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h @@ -16,6 +16,8 @@ namespace CGAL { +namespace internal +{ class K_means_point; @@ -81,13 +83,14 @@ public: int maximum_iteration; bool is_converged; protected: - int seed; + unsigned int seed; public: K_means_clustering(int number_of_centers, const std::vector& data, int number_of_run = 50, int maximum_iteration = 100) : points(data.begin(), data.end()), maximum_iteration(maximum_iteration), - is_converged(false), seed(time(NULL)) { + is_converged(false), + seed(static_cast(time(NULL))) { srand(seed); calculate_clustering_with_multiple_run(number_of_centers, number_of_run); } @@ -180,6 +183,7 @@ public: } } }; +}//namespace internal }//namespace CGAL #ifdef SEG_DEBUG #undef SEG_DEBUG