From 8e0538d2c7cb3099bf79c6839fbb6580347310cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=8Dlker=20Yaz?= Date: Tue, 19 Jun 2012 18:17:19 +0000 Subject: [PATCH] First version of alpha expansion graph cut. --- .../include/CGAL/Surface_mesh_segmentation.h | 31 +-- .../Alpha_expansion_graph_cut.h | 192 ++++++++++++++++-- .../Expectation_maximization.h | 3 +- 3 files changed, 190 insertions(+), 36 deletions(-) diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h index 64e2fce0692..e6fd8159ce4 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h @@ -126,6 +126,8 @@ public: bool use_minimum_segment; double multiplier_for_segment; static long miss_counter; + + internal::Expectation_maximization fitter; //member functions public: Surface_mesh_segmentation(Polyhedron* mesh, @@ -196,8 +198,8 @@ inline Surface_mesh_segmentation::Surface_mesh_segmentation( #endif // //write_sdf_values("sdf_values_sample_dino_ws.txt"); - //read_sdf_values("sdf_values_sample_elephant.txt"); - //apply_GMM_fitting(); + //read_sdf_values("sdf_values_sample_camel.txt"); + apply_GMM_fitting(); apply_graph_cut(); } @@ -292,15 +294,13 @@ Surface_mesh_segmentation::calculate_sdf_value_of_facet( segment_distance = min_distance; //first assignment of the segment_distance } else { // use segment_distance to limit rays as segments ray_direction = ray_direction / sqrt(ray_direction.squared_length()); - ray_direction = ray_direction * (*segment_distance * multiplier_for_segment); + Segment segment(center, CGAL::operator+(center, ray_direction)); boost::tie(is_intersected, intersection_is_acute, min_distance) = cast_and_return_minimum(segment, tree, facet); - if(!is_intersected) { - //continue; // for utopia case - just continue on miss + if(!is_intersected) { //no intersection is found ++miss_counter; - //Ray ray(segment.target(), ray_direction); Ray ray(segment.target(), ray_direction); boost::tie(is_intersected, intersection_is_acute, min_distance) = cast_and_return_minimum(ray, tree, facet); @@ -308,8 +308,9 @@ Surface_mesh_segmentation::calculate_sdf_value_of_facet( continue; } min_distance += CGAL::sqrt( - segment.squared_length()); // since we our source is segment.target() - } else if(!intersection_is_acute) { + segment.squared_length()); // since our source is segment.target() + } else if( + !intersection_is_acute) { // intersection is found, but it is not acceptable (so, do not continue ray-segment casting) continue; } @@ -897,8 +898,8 @@ inline void Surface_mesh_segmentation::apply_GMM_fitting() } SEG_DEBUG(CGAL::Timer t) SEG_DEBUG(t.start()) - // apply em with 5 runs, number of runs might become a parameter. - internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10); + //internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10); + fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 10); SEG_DEBUG(std::cout << t.time() << std::endl) std::vector center_memberships; fitter.fill_with_center_ids(center_memberships); @@ -927,7 +928,6 @@ inline void Surface_mesh_segmentation::apply_K_means_clustering() pair_it != sdf_values.end(); ++pair_it, ++center_it) { centers.insert(std::pair(pair_it->first, (*center_it))); } - //center_memberships_temp = center_memberships; //remove } template inline void @@ -944,8 +944,9 @@ Surface_mesh_segmentation::apply_GMM_fitting_with_K_means_init() std::vector center_memberships; clusterer.fill_with_center_ids(center_memberships); //std::vector center_memberships = center_memberships_temp; - internal::Expectation_maximization fitter(number_of_centers, sdf_vector, - center_memberships); + //internal::Expectation_maximization fitter(number_of_centers, sdf_vector, center_memberships); + fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, + center_memberships); center_memberships.clear(); fitter.fill_with_center_ids(center_memberships); std::vector::iterator center_it = center_memberships.begin(); @@ -976,7 +977,7 @@ 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 *= 15; //lambda, will be variable. + angle *= 10; //lambda, will be variable. // we may also want to consider edge lengths, also penalize convex angles. edge_weights.push_back(angle); } @@ -987,7 +988,7 @@ void Surface_mesh_segmentation::apply_graph_cut() pair_it != sdf_values.end(); ++pair_it) { sdf_vector.push_back(pair_it->second); } - internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 50); + //internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 3); //fill probability matrix. std::vector > probability_matrix; fitter.fill_with_minus_log_probabilities(probability_matrix); 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 8cbee29b8a9..b77a57b9288 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 @@ -43,10 +43,173 @@ public: typedef Traits::vertex_iterator Vertex_iterator; typedef Traits::edge_iterator Edge_iterator; - Alpha_expansion_graph_cut(std::vector >& edges, - std::vector& edge_weights, std::vector& labels, - std::vector >& probability_matrix, + Alpha_expansion_graph_cut(const std::vector >& edges, + const std::vector& edge_weights, std::vector& labels, + const std::vector >& probability_matrix, std::vector& center_ids) { + apply_alpha_expansion_2(edges, edge_weights, probability_matrix, labels); + center_ids = labels; + } + + boost::tuple + add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, + double w2, Graph& graph) { + Edge_descriptor v1_v2, v2_v1; + bool v1_v2_added, v2_v1_added; + + tie(v1_v2, v1_v2_added) = boost::add_edge(v1, v2, graph); + tie(v2_v1, v2_v1_added) = boost::add_edge(v2, v1, graph); + + CGAL_assertion(v1_v2_added && v2_v1_added); + //put edge capacities + boost::put(boost::edge_reverse, graph, v1_v2, v2_v1); + boost::put(boost::edge_reverse, graph, v2_v1, v1_v2); + + //map reverse edges + boost::put(boost::edge_capacity, graph, v1_v2, w1); + boost::put(boost::edge_capacity, graph, v2_v1, w2); + + return boost::make_tuple(v1_v2, v2_v1); + } + + 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) { + 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); + } + + 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_alpha_expansion_2(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) { + Vertex_descriptor v1 = inserted_vertices[edge_it->first], + v2 = inserted_vertices[edge_it->second]; + 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; + 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; + + add_edge_and_reverse(inbetween, v1, w1, w1, graph); + add_edge_and_reverse(inbetween, v2, w2, w2, graph); + } + } + + 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); @@ -86,8 +249,8 @@ public: boost::put(boost::edge_reverse, graph, sink_e, sink_e_rev); boost::put(boost::edge_reverse, graph, sink_e_rev, sink_e); } - std::vector::iterator weight_it = edge_weights.begin(); - for(std::vector >::iterator edge_it = edges.begin(); + 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]; @@ -101,18 +264,14 @@ public: CGAL_assertion(edge_added && edge_rev_added); //put edge capacities - if(labels[edge_it->first] == labels[edge_it->second]) { - boost::put(boost::edge_capacity, graph, edge, *weight_it); - boost::put(boost::edge_capacity, graph, edge_rev, *weight_it); - } else { - boost::put(boost::edge_capacity, graph, edge, *weight_it); - boost::put(boost::edge_capacity, graph, edge_rev, *weight_it); - } + 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); @@ -121,19 +280,12 @@ public: 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 + if(color == ColorTraits::black()) { // in sink center_ids.push_back(1); } else { center_ids.push_back(0); } } - std::cout << flow << std::endl; - //for(std::pair pair_v = boost::vertices(graph); - // pair_v.first != pair_v.second; ++(pair_v.first)) - //{ - // if(cluster_source == (*pair_v.first) || cluster_sink == (*pair_v.first)) { continue; } - // - //} } }; }//namespace internal 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 86f58f98965..fb96abf7c24 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 @@ -70,6 +70,7 @@ protected: std::vector > membership_matrix; public: + Expectation_maximization() { } /* 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) @@ -136,7 +137,7 @@ public: sum += probability; probabilities[center_i][point_i] = probability; } -#if 0 +#if 1 // pdf values scaled so that their sum will equal to 1. for(std::size_t point_i = 0; point_i < points.size(); ++point_i) { double probability = probabilities[center_i][point_i] / sum;