diff --git a/Polyhedron/demo/Polyhedron/Mesh_segmentation_dialog.ui b/Polyhedron/demo/Polyhedron/Mesh_segmentation_dialog.ui index 65100129821..2e200a1a02a 100644 --- a/Polyhedron/demo/Polyhedron/Mesh_segmentation_dialog.ui +++ b/Polyhedron/demo/Polyhedron/Mesh_segmentation_dialog.ui @@ -5,8 +5,8 @@ 0 0 - 328 - 149 + 352 + 151 @@ -15,30 +15,6 @@ - - - - 7 - - - - - - - 2 - - - - - - - Number of clusters: - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - @@ -49,6 +25,13 @@ + + + + 5 + + + @@ -70,19 +53,53 @@ - + - Show SDF values + Number of clusters: - - true + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - + + + 5 + + + + + - Show clusters + Smoothness: + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + 0.500000000000000 + + + 15.500000000000000 + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal @@ -104,9 +121,16 @@ - + - Apply + Caluclate SDF + + + + + + + Partition diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_mesh_segmentation_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_mesh_segmentation_plugin.cpp index 48afb54796f..14b141e59fe 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_mesh_segmentation_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_mesh_segmentation_plugin.cpp @@ -37,10 +37,11 @@ public: } } - void colorize(CGAL::Surface_mesh_segmentation& segmentation, std::vector& color_vector, bool sdf, bool cluster); + void colorize(CGAL::Surface_mesh_segmentation& segmentation, std::vector& color_vector, bool sdf); public slots: void on_actionSegmentation_triggered(); - void on_Apply_button_clicked(); + void on_Partition_button_clicked(); + void on_SDF_button_clicked(); private: QAction* actionSegmentation; @@ -63,17 +64,16 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_actionSegmentation_triggered() ui_dialog = new Ui::Mesh_segmentation_dialog(); ui_dialog->setupUi(dialog); - connect(ui_dialog->Apply_button, SIGNAL(clicked()), this, SLOT(on_Apply_button_clicked())); + connect(ui_dialog->Partition_button, SIGNAL(clicked()), this, SLOT(on_Partition_button_clicked())); + connect(ui_dialog->SDF_button, SIGNAL(clicked()), this, SLOT(on_SDF_button_clicked())); dialog->show(); } //void calculata_sdf_values(CGAL::Surface_mesh_segmentation* segmentation) //{ // segmentation->calculate_sdf_values(); //} - -void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked() +void Polyhedron_demo_mesh_segmentation_plugin::on_SDF_button_clicked() { - const Scene_interface::Item_id index = scene->mainSelectionIndex(); Scene_polyhedron_item* item = qobject_cast(scene->item(index)); Scene_polyhedron_with_color_item* item_colored = qobject_cast(scene->item(index)); @@ -84,8 +84,6 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked() int number_of_rays_sqrt = ui_dialog->Number_of_rays_spin_box->value(); double cone_angle = (ui_dialog->Cone_angle_spin_box->value() / 180.0) * CGAL_PI; int number_of_clusters = ui_dialog->Number_of_clusters_spin_box->value(); - bool show_sdf_values = ui_dialog->Show_sdf_values_check_box->isChecked(); - bool show_clusters = ui_dialog->Show_clusters_check_box->isChecked(); if(!item_colored) { @@ -99,7 +97,7 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked() CGAL::Surface_mesh_segmentation* segmentation = new CGAL::Surface_mesh_segmentation(new_item->polyhedron(), number_of_rays_sqrt, cone_angle, number_of_clusters); - colorize(*segmentation, color_vector, show_sdf_values, show_clusters); + colorize(*segmentation, color_vector, true); new_item->set_color_vector(color_vector); new_item->segmentation = segmentation; @@ -122,40 +120,73 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked() // ui_dialog->Sdf_value_calculation_bar->setValue(static_cast(segmentation->get_process_of_sdf_calculation() * 100)); segmentation->calculate_sdf_values(); //segmentation->apply_GMM_fitting_with_K_means_init(); - } - segmentation->number_of_centers = number_of_clusters; - segmentation->apply_GMM_fitting_with_K_means_init(); + } - colorize(*segmentation, color_vector, show_sdf_values, show_clusters); + colorize(*segmentation, color_vector, true); item_colored->set_color_vector(color_vector); scene->itemChanged(index); scene->setSelectedItem(index); } QApplication::restoreOverrideCursor(); +} + +void Polyhedron_demo_mesh_segmentation_plugin::on_Partition_button_clicked() +{ + + const Scene_interface::Item_id index = scene->mainSelectionIndex(); + Scene_polyhedron_with_color_item* item_colored = qobject_cast(scene->item(index)); + + if(!item_colored) { return; } + + QApplication::setOverrideCursor(Qt::WaitCursor); + int number_of_rays_sqrt = ui_dialog->Number_of_rays_spin_box->value(); + double cone_angle = (ui_dialog->Cone_angle_spin_box->value() / 180.0) * CGAL_PI; + int number_of_clusters = ui_dialog->Number_of_clusters_spin_box->value(); + double smoothness = ui_dialog->Smoothness_spin_box->value(); + + std::vector color_vector; + CGAL::Surface_mesh_segmentation* segmentation = item_colored->segmentation; + + segmentation->number_of_centers = number_of_clusters; + segmentation->apply_GMM_fitting_and_K_means(); + + segmentation->smoothing_lambda = smoothness; + segmentation->apply_graph_cut(); + segmentation->assign_segments(); + + colorize(*segmentation, color_vector, false); + item_colored->set_color_vector(color_vector); + scene->itemChanged(index); + scene->setSelectedItem(index); + + QApplication::restoreOverrideCursor(); //qApp->processEvents(); } void Polyhedron_demo_mesh_segmentation_plugin::colorize(CGAL::Surface_mesh_segmentation& segmentation, - std::vector& color_vector, bool sdf, bool cluster) + std::vector& color_vector, bool sdf) { + + QColor predefined_list[] = { Qt::red, Qt::darkRed, Qt::green, Qt::darkGreen, Qt::blue, Qt::darkBlue + , Qt::cyan, Qt::darkCyan, Qt::magenta, Qt::darkMagenta, Qt::yellow, Qt::darkYellow, Qt::gray, Qt::darkGray, Qt::lightGray, Qt::white, Qt::black}; + color_vector.reserve(segmentation.sdf_values.size()); int findex = 0; for(CGAL::Surface_mesh_segmentation::Facet_iterator facet_it = segmentation.mesh->facets_begin(); facet_it != segmentation.mesh->facets_end(); ++facet_it, ++findex) { facet_it->set_patch_id(findex); - int color = 0; - int sdf_color = (int) (255 * segmentation.sdf_values[facet_it]); - int cluster_color = (int) (255.0 / segmentation.number_of_centers) * segmentation.centers[facet_it]; - if(sdf && cluster) + //int cluster_color = (int) (255.0 / segmentation.number_of_centers) * segmentation.centers[facet_it]; + if(sdf) { - color = (sdf_color + cluster_color) / 2; + int sdf_color = (int) (255 * segmentation.sdf_values[facet_it]); + color_vector.push_back(QColor(sdf_color,sdf_color,sdf_color)); } else { - color = sdf ? sdf_color : cluster_color; + color_vector.push_back(predefined_list[segmentation.segments[facet_it] % 17]); } - color_vector.push_back(QColor(color,color,color)); + } } diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h index 490a91bcdd7..3dc7c36ea86 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h @@ -73,8 +73,8 @@ public: typedef typename Polyhedron::Edge_iterator Edge_iterator; typedef typename Polyhedron::Vertex_handle Vertex_handle; protected: - typedef typename Kernel::Ray_3 Ray; - typedef typename Kernel::Plane_3 Plane; + typedef typename Kernel::Ray_3 Ray; + typedef typename Kernel::Plane_3 Plane; typedef typename Kernel::Segment_3 Segment; typedef typename CGAL::AABB_polyhedron_triangle_primitive @@ -90,9 +90,10 @@ protected: typedef std::map Face_center_map; typedef std::map Edge_angle_map; typedef std::map Face_segment_map; + /*Sampled points from disk, t1 = coordinate-x, t2 = coordinate-y, t3 = angle with cone-normal (weight). */ - typedef CGAL::Triple Disk_sample; - typedef std::vector > Disk_samples_list; + typedef CGAL::Triple Disk_sample; + typedef std::vector Disk_samples_list; template struct Compare_second_element { @@ -151,12 +152,12 @@ public: boost::tuple cast_and_return_minimum_use_closest( const Query& ray, const Tree& tree, const Facet_handle& facet) const; - double calculate_sdf_value_from_rays (std::vector& ray_distances, - std::vector& ray_weights) const; - double calculate_sdf_value_from_rays_with_mean (std::vector& + double calculate_sdf_value_from_rays(std::vector& ray_distances, + std::vector& ray_weights) const; + double calculate_sdf_value_from_rays_with_mean(std::vector& + ray_distances, std::vector& ray_weights) const; + double calculate_sdf_value_from_rays_with_trimmed_mean(std::vector& ray_distances, std::vector& ray_weights) const; - double calculate_sdf_value_from_rays_with_trimmed_mean ( - std::vector& ray_distances, std::vector& ray_weights) const; void arrange_center_orientation(const Plane& plane, const Vector& unit_normal, Point& center) const; @@ -174,12 +175,13 @@ public: void apply_GMM_fitting(); void apply_K_means_clustering(); void apply_GMM_fitting_with_K_means_init(); + void apply_GMM_fitting_and_K_means(); + void apply_graph_cut(); void apply_graph_cut_with_EM(); void apply_graph_cut_multiple_run(int number_of_run = 5); void assign_segments(); - void dfs(Facet_handle facet, int segment_id); void write_sdf_values(const char* file_name); @@ -212,10 +214,10 @@ inline Surface_mesh_segmentation::Surface_mesh_segmentation( //write_sdf_values("sdf_values_sample_teddy.txt"); //write_sdf_values("18_3.txt"); - - //read_sdf_values("sdf_values_sample_camel.txt"); - apply_GMM_fitting_with_K_means_init(); - apply_graph_cut(); + //read_sdf_values("C:/Users/Ilker/Documents/Visual Studio 2008/Projects/SurfaceMeshSegmentation/SurfaceMeshSegmentation/models/sdf_values_sample_elephant.txt"); + //read_sdf_values("D:/GSoC/MeshsegBenchmark-1.0-full/MeshsegBenchmark-1.0/data/off/173.txt"); + //apply_GMM_fitting_with_K_means_init(); + //apply_graph_cut(); } template @@ -230,13 +232,6 @@ inline void Surface_mesh_segmentation::calculate_sdf_values() double sdf = calculate_sdf_value_of_facet(facet_it, tree); sdf_values.insert(std::pair(facet_it, sdf)); } - - //std::cout << Listing_intersection_traits_ray_or_segment_triangle - // < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list > >::inter_counter << std::endl; - //std::cout << Listing_intersection_traits_ray_or_segment_triangle - // < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list > >::true_inter_counter << std::endl; - //std::cout << Listing_intersection_traits_ray_or_segment_triangle - // < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list > >::do_inter_counter << std::endl; normalize_sdf_values(); smooth_sdf_values(); } @@ -348,16 +343,17 @@ Surface_mesh_segmentation::calculate_sdf_value_of_facet( ray_weights); } -// just for Ray and Segment template -template +template //Query can be templated for just Ray and Segment types. boost::tuple Surface_mesh_segmentation::cast_and_return_minimum( const Query& query, const Tree& tree, const Facet_handle& facet) const { - boost::tuple min_distance = boost::make_tuple(false, false, - 0); + // get<0> : if any intersection is found then true + // get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true + // get<2> : distance between ray/segment origin and intersection point. + boost::tuple min_distance(false, false, 0); std::list intersections; #if 1 //SL: the difference with all_intersections is that in the traversal traits, we do do_intersect before calling intersection. @@ -381,13 +377,9 @@ Surface_mesh_segmentation::cast_and_return_minimum( } const Point* i_point; - //Point i_point; - //AF: Use object_cast as it is faster than assign - //IOY: ok, also using pointer-returning version to get rid of copying if(!(i_point = CGAL::object_cast(&object))) { - continue; + continue; // Continue in case of segment. } - //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 = (query.source() - *i_point); double new_distance = i_ray.squared_length(); @@ -412,7 +404,6 @@ Surface_mesh_segmentation::cast_and_return_minimum( return min_distance; } min_distance.get<1>() = true; // founded intersection is acceptable. - CGAL_assertion(min_distance.get<2>() > 0); min_distance.get<2>() = sqrt(min_distance.get<2>()); return min_distance; } @@ -424,12 +415,10 @@ Surface_mesh_segmentation::cast_and_return_minimum_use_closest ( const Query& ray, const Tree& tree, const Facet_handle& facet) const { - //static double dist = 0.1; - //boost::optional min_distance_2 = dist++; - //return min_distance_2; - - boost::tuple min_distance = boost::make_tuple(false, false, - 0); + // get<0> : if any intersection is found then true + // get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true + // get<2> : distance between ray/segment origin and intersection point. + boost::tuple min_distance(false, false, 0); #if 1 Closest_intersection_traits traversal_traits; tree.traversal(ray, traversal_traits); @@ -733,7 +722,7 @@ double Surface_mesh_segmentation::calculate_dihedral_angle_of_edge_2 double n_angle = CGAL::Mesh_3::dihedral_angle(a, b, c, d) / 180.0; bool n_concave = n_angle > 0; double folded_angle = 1 + ((n_concave ? -1 : +1) * n_angle); - folded_angle = (CGAL::max)(folded_angle, epsilon); + folded_angle = (std::max)(folded_angle, epsilon); if(!n_concave) { return epsilon; // we may want to also penalize convex angles as well... @@ -905,28 +894,65 @@ inline void Surface_mesh_segmentation::smooth_sdf_values() } sdf_values = smoothed_sdf_values; } - template inline void Surface_mesh_segmentation::apply_GMM_fitting() { centers.clear(); std::vector sdf_vector; sdf_vector.reserve(sdf_values.size()); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - sdf_vector.push_back(pair_it->second); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it) { + sdf_vector.push_back(sdf_values[facet_it]); } SEG_DEBUG(CGAL::Timer t) SEG_DEBUG(t.start()) //internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10); - fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 15); + fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 20); SEG_DEBUG(std::cout << "GMM fitting time: " << t.time() << std::endl) std::vector center_memberships; fitter.fill_with_center_ids(center_memberships); std::vector::iterator center_it = center_memberships.begin(); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it, ++center_it) { - centers.insert(std::pair(pair_it->first, (*center_it))); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it, ++center_it) { + centers.insert(std::pair(facet_it, (*center_it))); + } +} + +template +inline void +Surface_mesh_segmentation::apply_GMM_fitting_and_K_means() +{ + centers.clear(); + std::vector sdf_vector; + sdf_vector.reserve(sdf_values.size()); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it) { + sdf_vector.push_back(sdf_values[facet_it]); + } + internal::Expectation_maximization gmm_random_init(number_of_centers, + sdf_vector, 20); + + internal::K_means_clustering k_means(number_of_centers, sdf_vector); + std::vector center_memberships; + k_means.fill_with_center_ids(center_memberships); + internal::Expectation_maximization gmm_k_means_init(number_of_centers, + sdf_vector, center_memberships); + + if(gmm_k_means_init.final_likelihood > gmm_random_init.final_likelihood) { + fitter = gmm_k_means_init; + } else { + fitter = gmm_random_init; + } + center_memberships.clear(); + fitter.fill_with_center_ids(center_memberships); + std::vector::iterator center_it = center_memberships.begin(); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it, ++center_it) { + centers.insert(std::pair(facet_it, (*center_it))); } } @@ -936,17 +962,19 @@ inline void Surface_mesh_segmentation::apply_K_means_clustering() centers.clear(); std::vector sdf_vector; sdf_vector.reserve(sdf_values.size()); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - sdf_vector.push_back(pair_it->second); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it) { + sdf_vector.push_back(sdf_values[facet_it]); } internal::K_means_clustering clusterer(number_of_centers, sdf_vector); std::vector center_memberships; clusterer.fill_with_center_ids(center_memberships); std::vector::iterator center_it = center_memberships.begin(); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it, ++center_it) { - centers.insert(std::pair(pair_it->first, (*center_it))); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it, ++center_it) { + centers.insert(std::pair(facet_it, (*center_it))); } //center_memberships_temp = center_memberships; //remove } @@ -957,9 +985,10 @@ Surface_mesh_segmentation::apply_GMM_fitting_with_K_means_init() centers.clear(); std::vector sdf_vector; sdf_vector.reserve(sdf_values.size()); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - sdf_vector.push_back(pair_it->second); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it) { + sdf_vector.push_back(sdf_values[facet_it]); } internal::K_means_clustering clusterer(number_of_centers, sdf_vector); std::vector center_memberships; @@ -971,9 +1000,10 @@ Surface_mesh_segmentation::apply_GMM_fitting_with_K_means_init() center_memberships.clear(); fitter.fill_with_center_ids(center_memberships); std::vector::iterator center_it = center_memberships.begin(); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it, ++center_it) { - centers.insert(std::pair(pair_it->first, (*center_it))); + for(Facet_iterator facet_it = mesh->facets_begin(); + facet_it != mesh->facets_end(); + ++facet_it, ++center_it) { + centers.insert(std::pair(facet_it, (*center_it))); } } @@ -998,11 +1028,12 @@ void Surface_mesh_segmentation::apply_graph_cut() int index_f2 = facet_indices[edge_it->opposite()->facet()]; edges.push_back(std::pair(index_f1, index_f2)); angle = -log(angle); - angle = (CGAL::max)(angle, 1e-5); + angle = (std::max)(angle, 1e-5); angle *= smoothing_lambda; //lambda, will be variable. // we may also want to consider edge lengths, also penalize convex angles. edge_weights.push_back(angle); } + //apply gmm fitting std::vector sdf_vector; sdf_vector.reserve(sdf_values.size()); diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h index 9865022329d..fa1670e1271 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h @@ -23,7 +23,6 @@ public: typedef boost::adjacency_list& edge_weights, const std::vector >& probability_matrix, std::vector& labels, double* result = NULL) { - double min_cut = apply_alpha_expansion_2(edges, edge_weights, - probability_matrix, labels); + double min_cut = apply_alpha_expansion(edges, edge_weights, probability_matrix, + labels); if(result != NULL) { *result = min_cut; } @@ -75,10 +74,10 @@ public: return boost::make_tuple(v1_v2, v2_v1); } - double apply_alpha_expansion_2(const std::vector >& edges, - const std::vector& edge_weights, - const std::vector >& probability_matrix, - std::vector& labels) { + double apply_alpha_expansion(const std::vector >& edges, + const std::vector& edge_weights, + const std::vector >& probability_matrix, + std::vector& labels) { int number_of_clusters = probability_matrix.size(); double min_cut = (std::numeric_limits::max)(); bool success; @@ -97,13 +96,14 @@ public: inserted_vertices.push_back(new_vertex); double source_weight = probability_matrix[alpha][vertex_i]; + // since it is expansion move, current alpha labeled vertices will be assigned to alpha again, + // making sink_weight 'infinity' guarantee this. double sink_weight = (labels[vertex_i] == alpha) ? (std::numeric_limits::max)() : probability_matrix[labels[vertex_i]][vertex_i]; - add_edge_and_reverse(cluster_source, new_vertex, source_weight, source_weight, - graph); - add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, sink_weight, graph); + 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); } // For E-Smooth // add edge between every vertex, @@ -115,17 +115,15 @@ 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 + // (their edges between sink is infitinity) double w1 = (label_1 == alpha) ? 0 : *weight_it; - //double w1 = *weight_it; add_edge_and_reverse(v1, v2, w1, w1, graph); } else { Vertex_descriptor inbetween = boost::add_vertex(graph); - add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); double w1 = (label_1 == alpha) ? 0 : *weight_it; double w2 = (label_2 == alpha) ? 0 : *weight_it; - //double w1 = *weight_it; - //double w2 = *weight_it; add_edge_and_reverse(inbetween, v1, w1, w1, graph); add_edge_and_reverse(inbetween, v2, w2, w2, graph); add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); @@ -144,7 +142,8 @@ public: for(std::size_t vertex_i = 0; vertex_i < inserted_vertices.size(); ++vertex_i) { boost::default_color_type color = boost::get(boost::vertex_color, graph, inserted_vertices[vertex_i]); - if(labels[vertex_i] != alpha && color == ColorTraits::white()) { //new comers + if(labels[vertex_i] != alpha + && color == ColorTraits::white()) { //new comers (expansion occurs) labels[vertex_i] = alpha; } } @@ -153,171 +152,6 @@ public: } while(success); return min_cut; } - -#if 0 - void apply_alpha_expansion(const std::vector >& edges, - const std::vector& edge_weights, - const std::vector >& probability_matrix, - std::vector& labels) { - int number_of_clusters = probability_matrix.size(); - double min_cut = (std::numeric_limits::max)(); - bool success; - do { - success = false; - for(int alpha = 0; alpha < number_of_clusters; ++alpha) { - Graph graph; - Vertex_descriptor cluster_source = boost::add_vertex(graph); - Vertex_descriptor cluster_sink = boost::add_vertex(graph); - std::vector inserted_vertices; - // For E-Data - // add every facet as a vertex to graph, put edges to source-sink vertices - for(std::size_t vertex_i = 0; vertex_i < probability_matrix[0].size(); - ++vertex_i) { - Vertex_descriptor new_vertex = boost::add_vertex(graph); - inserted_vertices.push_back(new_vertex); - - double source_weight = probability_matrix[alpha][vertex_i]; - double sink_weight = (labels[vertex_i] == alpha) ? - (std::numeric_limits::max)() : - probability_matrix[labels[vertex_i]][vertex_i]; - - 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); - } - // For E-Smooth - // add edge between every vertex, - std::vector::const_iterator weight_it = edge_weights.begin(); - for(std::vector >::const_iterator edge_it = edges.begin(); - edge_it != edges.end(); - ++edge_it, ++weight_it) { -#if 1 - Vertex_descriptor v1 = inserted_vertices[edge_it->first]; - Vertex_descriptor v2 = inserted_vertices[edge_it->second]; - add_edge_and_reverse(v1, v2, *weight_it, *weight_it, graph); -#else - int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second]; - if(label_1 == label_2) { - //double w1 = label_1 == alpha ? 0 : *weight_it; - double w1 = *weight_it; - add_edge_and_reverse(v1, v2, w1, w1, graph); - } else { - Vertex_descriptor inbetween = boost::add_vertex(graph); - add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); - - //double w1 = label_1 == alpha ? 0 : *weight_it; - //double w2 = label_2 == alpha ? 0 : *weight_it; - double w1 = *weight_it; - double w2 = *weight_it; - add_edge_and_reverse(inbetween, v1, w1, w1, graph); - add_edge_and_reverse(inbetween, v2, w2, w2, graph); - //add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); - } -#endif - } - - double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source, - cluster_sink); - if(min_cut - flow < flow * 1e-3) { - continue; - } - std::cout << "prev flow: " << min_cut << "new flow: " << flow << std::endl; - min_cut = flow; - success = true; - //update labeling - for(std::size_t vertex_i = 0; vertex_i < inserted_vertices.size(); ++vertex_i) { - boost::default_color_type color = boost::get(boost::vertex_color, graph, - inserted_vertices[vertex_i]); - if(labels[vertex_i] != alpha && color == ColorTraits::white()) { //new comers - labels[vertex_i] = alpha; - } - } - - } - } while(success); - } - - void apply_simple_cut(const std::vector >& edges, - const std::vector& edge_weights, std::vector& labels, - const std::vector >& probability_matrix, - std::vector& center_ids) { - //Graph graph(edges.begin(), edges.end(), vertex_weights.size(), edge_weights.size()); - Graph graph; - Vertex_descriptor cluster_source = boost::add_vertex(graph); - Vertex_descriptor cluster_sink = boost::add_vertex(graph); - std::vector inserted_vertices; - //inserted_vertices.reserve(vertex_weights.size()); - for(std::size_t vertex_i = 0; vertex_i < probability_matrix[0].size(); - ++vertex_i) { - Vertex_descriptor new_vertex = boost::add_vertex(graph); - inserted_vertices.push_back(new_vertex); - - Edge_descriptor source_e, source_e_rev, sink_e, sink_e_rev; - bool source_e_added, source_e_rev_added, sink_e_added, sink_e_rev_added; - - tie(source_e, source_e_added) = boost::add_edge(cluster_source, new_vertex, - graph); - tie(source_e_rev, source_e_rev_added) = boost::add_edge(new_vertex, - cluster_source, graph); - - tie(sink_e, sink_e_added) = boost::add_edge(new_vertex, cluster_sink, graph); - tie(sink_e_rev, sink_e_rev_added) = boost::add_edge(cluster_sink, new_vertex, - graph); - - CGAL_assertion(source_e_added && source_e_rev_added && sink_e_added - && sink_e_rev_added); - - //put vertex capacities (to edges between itself and source & sink) - boost::put(boost::edge_capacity, graph, source_e, - probability_matrix[0][vertex_i]); - boost::put(boost::edge_capacity, graph, source_e_rev, 0); - boost::put(boost::edge_capacity, graph, sink_e, - probability_matrix[1][vertex_i]); - boost::put(boost::edge_capacity, graph, sink_e_rev, 0); - //map reverse edges - boost::put(boost::edge_reverse, graph, source_e, source_e_rev); - boost::put(boost::edge_reverse, graph, source_e_rev, source_e); - boost::put(boost::edge_reverse, graph, sink_e, sink_e_rev); - boost::put(boost::edge_reverse, graph, sink_e_rev, sink_e); - } - std::vector::const_iterator weight_it = edge_weights.begin(); - for(std::vector >::const_iterator edge_it = edges.begin(); - edge_it != edges.end(); - ++edge_it, ++weight_it) { - Vertex_descriptor v1 = inserted_vertices[edge_it->first]; - Vertex_descriptor v2 = inserted_vertices[edge_it->second]; - - bool edge_added, edge_rev_added; - Edge_descriptor edge, edge_rev; - tie(edge, edge_added) = boost::add_edge(v1, v2, graph); - tie(edge_rev, edge_rev_added) = boost::add_edge(v2, v1, graph); - - CGAL_assertion(edge_added && edge_rev_added); - - //put edge capacities - boost::put(boost::edge_capacity, graph, edge, *weight_it); - boost::put(boost::edge_capacity, graph, edge_rev, *weight_it); - - //map reverse edges - boost::put(boost::edge_reverse, graph, edge, edge_rev); - boost::put(boost::edge_reverse, graph, edge_rev, edge); - } - - double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source, - cluster_sink); - - for(std::vector::iterator vertex_it = - inserted_vertices.begin(); - vertex_it != inserted_vertices.end(); ++vertex_it) { - boost::default_color_type color = boost::get(boost::vertex_color, graph, - *vertex_it); - if(color == ColorTraits::black()) { // in sink - center_ids.push_back(1); - } else { - center_ids.push_back(0); - } - } - } -#endif }; }//namespace internal }//namespace CGAL diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut_with_EM.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut_with_EM.h index c42804d38e3..5f23fd79701 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut_with_EM.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut_with_EM.h @@ -122,16 +122,12 @@ public: int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second]; if(label_1 == label_2) { double w1 = (label_1 == alpha) ? 0 : *weight_it; - //double w1 = *weight_it; add_edge_and_reverse(v1, v2, w1, w1, graph); } else { Vertex_descriptor inbetween = boost::add_vertex(graph); - add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); double w1 = (label_1 == alpha) ? 0 : *weight_it; double w2 = (label_2 == alpha) ? 0 : *weight_it; - //double w1 = *weight_it; - //double w2 = *weight_it; add_edge_and_reverse(inbetween, v1, w1, w1, graph); add_edge_and_reverse(inbetween, v2, w2, w2, graph); add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph); @@ -156,13 +152,13 @@ public: } } labels_result = labels; - } - Expectation_maximization em(probability_matrix.size(), data, labels); - labels.clear(); - probability_matrix.clear(); - em.fill_with_center_ids(labels); - em.fill_with_minus_log_probabilities(probability_matrix); + + } + int matrix_size = probability_matrix.size(); + probability_matrix.clear(); + Expectation_maximization(data).refresh_parameters_and_probabilities(matrix_size, + labels, probability_matrix); } while(success); // recalculate clustering labels = labels_result; 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 f4dbb17da00..81d5298e4d8 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 @@ -19,6 +19,7 @@ #define DEF_MAX_ITER 10 #define DEF_THRESHOLD 1e-4 #define USE_MATRIX true +#define DEF_SEED 1340818006 #define ACTIVATE_SEGMENTATION_EM_DEBUG #ifdef ACTIVATE_SEGMENTATION_EM_DEBUG @@ -65,19 +66,22 @@ public: double threshold; int maximum_iteration; bool is_converged; + double final_likelihood; protected: unsigned int seed; std::vector > membership_matrix; public: Expectation_maximization() { } + Expectation_maximization(const std::vector& data): points(data) { } /* For uniform initialization, with one run */ Expectation_maximization(int number_of_centers, const std::vector& data, double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER) : points(data), threshold(threshold), maximum_iteration(maximum_iteration), is_converged(false), seed(0), membership_matrix(std::vector >(number_of_centers, - std::vector(points.size()))) { + std::vector(points.size()))), + final_likelihood(-(std::numeric_limits::max)()) { initiate_centers_uniformly(number_of_centers); fit(); } @@ -86,7 +90,8 @@ public: std::vector& initial_center_ids, double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER) : points(data), threshold(threshold), maximum_iteration(maximum_iteration), - is_converged(false), seed(0) { + is_converged(false), seed(0), + final_likelihood(-(std::numeric_limits::max)()) { number_of_centers = initiate_centers_from_memberships(number_of_centers, initial_center_ids); membership_matrix = std::vector >(number_of_centers, @@ -101,8 +106,13 @@ public: is_converged(false), seed(static_cast(time(NULL))), membership_matrix(std::vector >(number_of_centers, - std::vector(points.size()))) { + std::vector(points.size()))), + final_likelihood(-(std::numeric_limits::max)()) { +#if 0 srand(seed); +#else + srand(DEF_SEED); +#endif fit_with_multiple_run(number_of_centers, number_of_runs); } @@ -139,9 +149,9 @@ public: } for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { double probability = probabilities[center_i][point_i] / total_probability; - probability = (CGAL::max)(probability, epsilon); + probability = (std::max)(probability, epsilon); probability = -log(probability); - probabilities[center_i][point_i] = (CGAL::max)(probability, + probabilities[center_i][point_i] = (std::max)(probability, 1e-5); // this is another epsilon, edge can not hold 0 to source in graph-cut. } } @@ -173,12 +183,82 @@ public: // { // double probability = probabilities[center_i][point_i]; // probability = (probability - min_value) / max_min_dif; - // probability = (CGAL::max)(probability, epsilon); + // probability = (std::max)(probability, epsilon); // probabilities[center_i][point_i] = -log(probability); // } // #endif //} } + + // experimental, going to be removed most prob ! + void refresh_parameters_and_probabilities(int number_of_centers, + std::vector& initial_center_ids, + std::vector >& probabilities) { + std::vector cluster_exist(number_of_centers, false); + std::vector cluster_shift(number_of_centers, 0); + for(std::vector::iterator id_it = initial_center_ids.begin(); + id_it != initial_center_ids.end(); ++id_it) { + cluster_exist[*id_it] = true; + } + + /* Calculate mean */ + 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_points; ++i) { + int center_id = initial_center_ids[i]; + double data = points[i]; + centers[center_id].mean += data; + member_count[center_id] += 1; + } + /* Assign mean, and mixing coef */ + for(int i = 0; i < number_of_centers; ++i) { + if(!cluster_exist[i]) { + continue; + } + centers[i].mean /= member_count[i]; + centers[i].mixing_coefficient = member_count[i] / static_cast + (number_of_points); + } + /* Calculate deviation */ + for(int i = 0; i < number_of_points; ++i) { + int center_id = initial_center_ids[i]; + double data = points[i]; + centers[center_id].deviation += pow(data - centers[center_id].mean, 2); + } + for(int i = 0; i < number_of_centers; ++i) { + if(!cluster_exist[i]) { + continue; + } + centers[i].deviation = sqrt(centers[i].deviation / member_count[i]); + } + + double epsilon = 1e + -5; // this epsilon should be consistent with epsilon in calculate_dihedral_angle_of_edge! + probabilities = std::vector > + (number_of_centers, std::vector(points.size())); + for(std::size_t point_i = 0; point_i < points.size(); ++point_i) { + double total_probability = 0.0; + for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { + if(!cluster_exist[center_i]) { + probabilities[center_i][point_i] = 0; + continue; + } + double probability = centers[center_i].probability(points[point_i]); + total_probability += probability; + probabilities[center_i][point_i] = probability; + } + for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) { + double probability = probabilities[center_i][point_i] / total_probability; + probability = (std::max)(probability, epsilon); + probability = -log(probability); + probabilities[center_i][point_i] = (std::max)(probability, + 1e-5); // this is another epsilon, edge can not hold 0 to source in graph-cut. + } + } + //sort(centers.begin(), centers.end()); + } protected: // Initialization functions for centers. @@ -213,6 +293,7 @@ protected: sort(centers.begin(), centers.end()); } + int initiate_centers_from_memberships(int number_of_centers, std::vector& initial_center_ids) { /* For handling clusters that have 0 members (in initial_center_ids) */ @@ -336,29 +417,28 @@ protected: while(!is_converged && iteration_count++ < maximum_iteration) { prev_likelihood = likelihood; likelihood = iterate(iteration_count == 1); - is_converged = likelihood - prev_likelihood < threshold * fabs(likelihood); + is_converged = likelihood - prev_likelihood < threshold * std::fabs(likelihood); } SEG_DEBUG(std::cout << "likelihood: " << likelihood << "iteration: " << iteration_count << std::endl) + if(final_likelihood < likelihood) { + final_likelihood = likelihood; + } return likelihood; } // Calls fit() number_of_run times, and stores best found centers void fit_with_multiple_run(int number_of_centers, int number_of_run) { - double max_likelihood = -(std::numeric_limits::max)(); std::vector max_centers; while(number_of_run-- > 0) { initiate_centers_randomly(number_of_centers); - double likelihood = fit(); - //write_center_parameters("center_paramters.txt"); - if(likelihood > max_likelihood) { + if(likelihood == final_likelihood) { max_centers = centers; - max_likelihood = likelihood; } } - SEG_DEBUG(std::cout << "max likelihood: " << max_likelihood << std::endl) + SEG_DEBUG(std::cout << "max likelihood: " << final_likelihood << std::endl) centers = max_centers; } @@ -383,6 +463,7 @@ protected: }; }//namespace internal }//namespace CGAL +#undef DEF_SEED #undef DEF_MAX_ITER #undef DEF_THRESHOLD #undef USE_MATRIX 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 fa821e38e79..2fa0c0ac481 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 @@ -6,7 +6,9 @@ #include #include #include +#include +#define DEF_SEED 1340818006 //#define ACTIVATE_SEGMENTATION_K_MEANS_DEBUG #ifdef ACTIVATE_SEGMENTATION_K_MEANS_DEBUG #define SEG_DEBUG(x) x; @@ -59,9 +61,9 @@ public: * Returns true if center is changed.*/ bool calculate_new_center(std::vector& centers) { int new_center_id = 0; - double min_distance = fabs(centers[0].mean - data); + double min_distance = std::fabs(centers[0].mean - data); for(std::size_t i = 1; i < centers.size(); ++i) { - double new_distance = fabs(centers[i].mean - data); + double new_distance = std::fabs(centers[i].mean - data); if(new_distance < min_distance) { new_center_id = i; min_distance = new_distance; @@ -87,11 +89,15 @@ protected: public: K_means_clustering(int number_of_centers, const std::vector& data, - int number_of_run = 100, int maximum_iteration = 100) + int number_of_run = 20, int maximum_iteration = 100) : points(data.begin(), data.end()), maximum_iteration(maximum_iteration), is_converged(false), seed(static_cast(time(NULL))) { +#if 0 srand(seed); +#else + srand(DEF_SEED); +#endif calculate_clustering_with_multiple_run(number_of_centers, number_of_run); } /* For each data point, data_center is filled by its center's id. */ @@ -121,6 +127,35 @@ public: sort(centers.begin(), centers.end()); } + void initiate_centers_k_means_plus_plus(int number_of_centers) { + centers.clear(); + std::vector distance_square_cumulative(points.size()); + std::vector distance_square(points.size(), + (std::numeric_limits::max)()); + double initial_mean = points[rand() % points.size()].data; + centers.push_back(K_means_center(initial_mean)); + + for(int i = 1; i < number_of_centers; ++i) { + double cumulative_distance_square = 0.0; + for(std::size_t j = 0; j < points.size(); ++j) { + double new_distance = pow(centers.back().mean - points[j].data, 2); + if(new_distance < distance_square[j]) { + distance_square[j] = new_distance; + } + cumulative_distance_square += distance_square[j]; + distance_square_cumulative[j] = cumulative_distance_square; + } + + double random_ds = (rand() / static_cast(RAND_MAX)) * + (distance_square_cumulative.back()); + int selection_index = lower_bound(distance_square_cumulative.begin(), + distance_square_cumulative.end(), random_ds) + - distance_square_cumulative.begin(); + double initial_mean = points[selection_index].data; + centers.push_back(K_means_center(initial_mean)); + } + } + void calculate_clustering() { int iteration_count = 0; bool any_center_changed = true; @@ -149,8 +184,12 @@ public: double error = (std::numeric_limits::max)(); while(number_of_run-- > 0) { clear_center_ids(); - +#if 0 initiate_centers_randomly(number_of_centers); +#else + initiate_centers_k_means_plus_plus(number_of_centers); +#endif + calculate_clustering(); #if 0 double new_error = sum_of_squares(); @@ -198,6 +237,7 @@ public: }; }//namespace internal }//namespace CGAL +#undef DEF_SEED #ifdef SEG_DEBUG #undef SEG_DEBUG #endif