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 8e3c3f37d6f..14bec111050 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 @@ -50,21 +50,23 @@ #include #include -#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE -#include -#include -#if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated. -#include -#else -#include -#endif -#else + +#ifndef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE namespace MaxFlow { #include } #endif +#include +#include + +#if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated. +# include +#else +# include +#endif + #include @@ -196,14 +198,15 @@ struct Alpha_expansion_old_API_wrapper_graph // since it is done by exploring all edges to find opposite //////////////////////////////////////////////////////////////////////////////////////// -#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE +} // namespace internal + /** * @brief Implements alpha-expansion graph cut algorithm. * * For representing graph, it uses adjacency_list with OutEdgeList = vecS, VertexList = listS. * Also no pre-allocation is made for vertex-list. */ -class Alpha_expansion_graph_cut_boost +class Alpha_expansion_boost_adjacency_list { private: typedef boost::adjacency_list_traits @@ -225,11 +228,73 @@ private: typedef boost::graph_traits Traits; typedef boost::color_traits ColorTraits; +public: + typedef Traits::vertex_descriptor Vertex_descriptor; typedef Traits::vertex_iterator Vertex_iterator; typedef Traits::edge_descriptor Edge_descriptor; typedef Traits::edge_iterator Edge_iterator; +private: + + Graph graph; + Vertex_descriptor cluster_source; + Vertex_descriptor cluster_sink; + +public: + + void clear_graph() + { + graph.clear(); + cluster_source = boost::add_vertex(graph); + cluster_sink = boost::add_vertex(graph); + } + + Vertex_descriptor add_vertex() + { + return boost::add_vertex(graph); + } + + void add_tweight (Vertex_descriptor& v, double w1, double w2) + { + add_edge (cluster_source, v, w1, 0); + add_edge (v, cluster_sink, w2, 0); + } + + void init_vertices() + { + // initialize vertex indices, it is necessary since we are using VertexList = listS + Vertex_iterator v_begin, v_end; + Traits::vertices_size_type index = 0; + for(boost::tie(v_begin, v_end) = vertices(graph); v_begin != v_end; ++v_begin) { + boost::put(boost::vertex_index, graph, *v_begin, index++); + } + } + + double max_flow() + { +#if BOOST_VERSION >= 104400 + return boost::boykov_kolmogorov_max_flow(graph, cluster_source, + cluster_sink); +#else + return boost::kolmogorov_max_flow(graph, cluster_source, cluster_sink); +#endif + } + + template + void update(VertexLabelMap vertex_label_map, + const std::vector& inserted_vertices, + InputVertexDescriptor vd, + std::size_t vertex_i, + std::size_t alpha) + { + boost::default_color_type color = boost::get(boost::vertex_color, graph, + inserted_vertices[vertex_i]); + if(get (vertex_label_map, vd) != alpha + && color == ColorTraits::white()) //new comers (expansion occurs) + put (vertex_label_map, vd, alpha); + } + /** * Adds two directional edges between @a v1 and @a v2 * @param v1 first vertex @@ -239,9 +304,8 @@ private: * @param graph to be added * @return pair of added edges, first: v1->v2 and second: v2->v1 */ - boost::tuple - add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, - double w2, Graph& graph) const { + void add_edge (Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, double w2) + { Edge_descriptor v1_v2, v2_v1; bool v1_v2_added, v2_v1_added; @@ -257,194 +321,13 @@ private: 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); - } - - - -public: - /** - * Applies alpha-expansion graph-cut for energy minimization. - * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1]) - * @param edge_weights contains weights for each edge in @a edges (correspondence according to order) - * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex] - * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]), - * and as output it returns final labeling of vertices (i.e. assigned cluster-id to each facet) - * @return result of energy function - */ - double operator()(const std::vector >& - edges, - const std::vector& edge_weights, - const std::vector >& probability_matrix, - std::vector& labels) const - { - Alpha_expansion_old_API_wrapper_graph graph (edges, edge_weights, probability_matrix, labels); - - return (*this)(graph, - graph.edge_weight_map(), - graph.vertex_index_map(), - graph.vertex_label_map(), - graph.vertex_label_probability_map()); - } - - template - double operator()(const InputGraph& input_graph, - Edge_weight_map edge_weight_map, - Vertex_index_map vertex_index_map, - Vertex_label_map vertex_label_map, - Vertex_label_probability_map vertex_label_probability_map) const - { - typedef boost::graph_traits GT; - typedef typename GT::edge_descriptor input_edge_descriptor; - typedef typename GT::vertex_descriptor input_vertex_descriptor; - - // TODO: check this hardcoded parameter - const double tolerance = 1e-10; - - double min_cut = (std::numeric_limits::max)(); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - double vertex_creation_time, edge_creation_time, cut_time; - vertex_creation_time = edge_creation_time = cut_time = 0.0; - #endif - - std::vector inserted_vertices; - inserted_vertices.resize(num_vertices (input_graph)); - - std::size_t number_of_labels = get(vertex_label_probability_map, *std::begin(vertices(input_graph))).size(); - - Graph graph; - - bool success; - do { - success = false; - - for (std::size_t alpha = 0; alpha < number_of_labels; ++ alpha) - { - graph.clear(); - - Vertex_descriptor cluster_source = boost::add_vertex(graph); - Vertex_descriptor cluster_sink = boost::add_vertex(graph); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - Timer timer; - timer.start(); - #endif - - // For E-Data - // add every input vertex as a vertex to the graph, put edges to source & sink vertices - for (input_vertex_descriptor vd : vertices(input_graph)) - { - std::size_t vertex_i = get(vertex_index_map, vd); - Vertex_descriptor new_vertex = boost::add_vertex(graph); - inserted_vertices[vertex_i] = new_vertex; - double source_weight = get(vertex_label_probability_map, vd)[alpha]; - // since it is expansion move, current alpha labeled vertices will be assigned to alpha again, - // making sink_weight 'infinity' guarantee this. - double sink_weight = (get(vertex_label_map, vd) == alpha ? - (std::numeric_limits::max)() - : get(vertex_label_probability_map, vd)[get(vertex_label_map, vd)]); - - 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); - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - vertex_creation_time += timer.time(); - timer.reset(); - #endif - - // For E-Smooth - // add edge between every vertex, - for (input_edge_descriptor ed : edges(input_graph)) - { - input_vertex_descriptor vd1 = source(ed, input_graph); - input_vertex_descriptor vd2 = target(ed, input_graph); - std::size_t idx1 = get (vertex_index_map, vd1); - std::size_t idx2 = get (vertex_index_map, vd2); - - double weight = get (edge_weight_map, ed); - - Vertex_descriptor v1 = inserted_vertices[idx1], - v2 = inserted_vertices[idx2]; - - std::size_t label_1 = get (vertex_label_map, vd1); - std::size_t label_2 = get (vertex_label_map, vd2); - if(label_1 == label_2) { - if(label_1 != alpha) { - add_edge_and_reverse(v1, v2, weight, weight, graph); - } - } else { - Vertex_descriptor inbetween = boost::add_vertex(graph); - - double w1 = (label_1 == alpha) ? 0 : weight; - double w2 = (label_2 == alpha) ? 0 : weight; - 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, 0.0, graph); - } - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - edge_creation_time += timer.time(); - #endif - - // initialize vertex indices, it is necessary since we are using VertexList = listS - Vertex_iterator v_begin, v_end; - Traits::vertices_size_type index = 0; - for(boost::tie(v_begin, v_end) = vertices(graph); v_begin != v_end; ++v_begin) { - boost::put(boost::vertex_index, graph, *v_begin, index++); - } - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - timer.reset(); - #endif -#if BOOST_VERSION >= 104400 - double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source, - cluster_sink); -#else - double flow = boost::kolmogorov_max_flow(graph, cluster_source, cluster_sink); -#endif - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - cut_time += timer.time(); - #endif - - if(min_cut - flow <= flow * tolerance) { - continue; - } - min_cut = flow; - success = true; - //update labeling - for (input_vertex_descriptor vd : vertices (input_graph)) - { - std::size_t vertex_i = get (vertex_index_map, vd); - boost::default_color_type color = boost::get(boost::vertex_color, graph, - inserted_vertices[vertex_i]); - if(get (vertex_label_map, vd) != alpha - && color == ColorTraits::white()) { //new comers (expansion occurs) - put (vertex_label_map, vd, alpha); - } - } - } - } while(success); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << - std::endl; - CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; - CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; - #endif - - return min_cut; } }; // another implementation using compressed_sparse_row_graph // for now there is a performance problem while setting reverse edges // if that can be solved, it is faster than Alpha_expansion_graph_cut_boost -class Alpha_expansion_graph_cut_boost_CSR +class Alpha_expansion_boost_compressed_sparse_row { private: // CSR only accepts bundled props @@ -470,217 +353,127 @@ private: typedef boost::graph_traits Traits; typedef boost::color_traits ColorTraits; +public: + typedef Traits::vertex_descriptor Vertex_descriptor; typedef Traits::vertex_iterator Vertex_iterator; typedef Traits::edge_descriptor Edge_descriptor; typedef Traits::edge_iterator Edge_iterator; - void - add_edge_and_reverse(std::size_t v1 , std::size_t v2, double w1, double w2, - std::vector >& edge_map, - std::vector& edge_weights) const { +private: + + Graph graph; + std::size_t nb_vertices; + std::vector > edge_map; + std::vector edge_map_weights; + +public: + void clear_graph() + { + nb_vertices = 2; + edge_map.clear(); + edge_map_weights.clear(); + // edge_map.reserve(labels.size() * + // 8); // there is no way to know exact edge count, it is a heuristic value + // edge_map_weights.reserve(labels.size() * 8); + } + + Vertex_descriptor add_vertex() + { + return (nb_vertices ++); + } + + void add_tweight (Vertex_descriptor& v, double w1, double w2) + { + add_edge (0, v, w1, 0); + add_edge (v, 1, w2, 0); + } + + void init_vertices() + { +#if BOOST_VERSION >= 104000 + graph = Graph(boost::edges_are_unsorted, edge_map.begin(), edge_map.end(), + edge_map_weights.begin(), nb_vertices); +#else + graph= Graph(edge_map.begin(), edge_map.end(), + edge_map_weights.begin(), nb_vertices); +#endif + + // PERFORMANCE PROBLEM + // need to set reverse edge map, I guess there is no way to do that before creating the graph + // since we do not have edge_descs + // however from our edge_map, we know that each (2i, 2i + 1) is reverse pairs, how to facilitate that ? + // will look it back + Graph::edge_iterator ei, ee; + for(boost::tie(ei, ee) = boost::edges(graph); ei != ee; ++ei) { + Graph::vertex_descriptor v1 = boost::source(*ei, graph); + Graph::vertex_descriptor v2 = boost::target(*ei, graph); + std::pair opp_edge = boost::edge(v2, v1, graph); + + CGAL_assertion(opp_edge.second); + graph[opp_edge.first].edge_reverse = + *ei; // and edge_reverse of *ei will be (or already have been) set by the opp_edge + } + } + + double max_flow() + { +#if BOOST_VERSION >= 104400 + // since properties are bundled, defaults does not work need to specify them + return boost::boykov_kolmogorov_max_flow + (graph, + boost::get(&EdgeP::edge_capacity, graph), + boost::get(&EdgeP::edge_residual_capacity, graph), + boost::get(&EdgeP::edge_reverse, graph), + boost::get(&VertexP::vertex_predecessor, graph), + boost::get(&VertexP::vertex_color, graph), + boost::get(&VertexP::vertex_distance_t, graph), + boost::get(boost::vertex_index, + graph), // this is not bundled, get it from graph (CRS provides one) + 0, 1); +#else + return boost::kolmogorov_max_flow + (graph, + boost::get(&EdgeP::edge_capacity, graph), + boost::get(&EdgeP::edge_residual_capacity, graph), + boost::get(&EdgeP::edge_reverse, graph), + boost::get(&VertexP::vertex_predecessor, graph), + boost::get(&VertexP::vertex_color, graph), + boost::get(&VertexP::vertex_distance_t, graph), + boost::get(boost::vertex_index, + graph), // this is not bundled, get it from graph + 0, 1); +#endif + } + + template + void update(VertexLabelMap vertex_label_map, + const std::vector& inserted_vertices, + InputVertexDescriptor vd, + std::size_t vertex_i, + std::size_t alpha) + { + boost::default_color_type color = graph[vertex_i + 2].vertex_color; + if(get(vertex_label_map, vd)!= alpha + && color == ColorTraits::white()) //new comers (expansion occurs) + put(vertex_label_map, vd, alpha); + } + + void add_edge(Vertex_descriptor v1, Vertex_descriptor v2, double w1, double w2) + { edge_map.push_back(std::make_pair(v1, v2)); EdgeP p1; p1.edge_capacity = w1; - edge_weights.push_back(p1); + edge_map_weights.push_back(p1); edge_map.push_back(std::make_pair(v2, v1)); EdgeP p2; p2.edge_capacity = w2; - edge_weights.push_back(p2); + edge_map_weights.push_back(p2); } + -public: - /** - * Applies alpha-expansion graph-cut for energy minimization. - * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1]) - * @param edge_weights contains weights for each edge in @a edges (correspondence according to order) - * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex] - * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]), - * and as output it returns final labeling of vertices (i.e. assigned cluster-id to each facet) - * @return result of energy function - */ - double operator()(const std::vector >& - edges, - const std::vector& edge_weights, - const std::vector >& probability_matrix, - std::vector& labels) const { - const double tolerance = 1e-10; - - double min_cut = (std::numeric_limits::max)(); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - double vertex_creation_time, edge_creation_time, graph_creation_time, - reverse_mapping_time, cut_time; - vertex_creation_time = edge_creation_time = graph_creation_time = - reverse_mapping_time = cut_time = 0.0; - #endif - - Graph graph; - - bool success; - do { - success = false; - std::size_t alpha = 0; - - for(std::vector >::const_iterator it = - probability_matrix.begin(); - it != probability_matrix.end(); ++it, ++alpha) { - std::vector > edge_map; - std::vector edge_map_weights; - edge_map.reserve(labels.size() * - 8); // there is no way to know exact edge count, it is a heuristic value - edge_map_weights.reserve(labels.size() * 8); - - std::size_t cluster_source = 0; - std::size_t cluster_sink = 1; - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - Timer timer; - timer.start(); - #endif - // For E-Data - // add every facet as a vertex to the graph, put edges to source & sink vertices - for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) { - 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, vertex_i + 2, source_weight, 0.0, edge_map, - edge_map_weights); - add_edge_and_reverse(vertex_i + 2, cluster_sink, sink_weight, 0.0, edge_map, - edge_map_weights); - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - vertex_creation_time += timer.time(); - timer.reset(); - #endif - // For E-Smooth - // add edge between every vertex, - std::size_t num_vert = labels.size() + 2; - 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) { - std::size_t v1 = edge_it->first + 2, v2 = edge_it->second + 2; - std::size_t label_1 = labels[edge_it->first], label_2 = labels[edge_it->second]; - if(label_1 == label_2) { - if(label_1 != alpha) { - add_edge_and_reverse(v1, v2, *weight_it, *weight_it, edge_map, - edge_map_weights); - } - } else { - std::size_t inbetween = num_vert++; - - double w1 = (label_1 == alpha) ? 0 : *weight_it; - double w2 = (label_2 == alpha) ? 0 : *weight_it; - add_edge_and_reverse(inbetween, v1, w1, w1, edge_map, edge_map_weights); - add_edge_and_reverse(inbetween, v2, w2, w2, edge_map, edge_map_weights); - add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, edge_map, - edge_map_weights); - } - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - edge_creation_time += timer.time(); - timer.reset(); - #endif -#if BOOST_VERSION >= 104000 - Graph graph(boost::edges_are_unsorted, edge_map.begin(), edge_map.end(), - edge_map_weights.begin(), num_vert); -#else - Graph graph(edge_map.begin(), edge_map.end(), - edge_map_weights.begin(), num_vert); -#endif - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - graph_creation_time += timer.time(); - timer.reset(); - #endif - - // PERFORMANCE PROBLEM - // need to set reverse edge map, I guess there is no way to do that before creating the graph - // since we do not have edge_descs - // however from our edge_map, we know that each (2i, 2i + 1) is reverse pairs, how to facilitate that ? - // will look it back - Graph::edge_iterator ei, ee; - for(boost::tie(ei, ee) = boost::edges(graph); ei != ee; ++ei) { - Graph::vertex_descriptor v1 = boost::source(*ei, graph); - Graph::vertex_descriptor v2 = boost::target(*ei, graph); - std::pair opp_edge = boost::edge(v2, v1, graph); - - CGAL_assertion(opp_edge.second); - graph[opp_edge.first].edge_reverse = - *ei; // and edge_reverse of *ei will be (or already have been) set by the opp_edge - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - reverse_mapping_time += timer.time(); - timer.reset(); - #endif - -#if BOOST_VERSION >= 104400 - // since properties are bundled, defaults does not work need to specify them - double flow = boost::boykov_kolmogorov_max_flow(graph, - boost::get(&EdgeP::edge_capacity, graph), - boost::get(&EdgeP::edge_residual_capacity, graph), - boost::get(&EdgeP::edge_reverse, graph), - boost::get(&VertexP::vertex_predecessor, graph), - boost::get(&VertexP::vertex_color, graph), - boost::get(&VertexP::vertex_distance_t, graph), - boost::get(boost::vertex_index, - graph), // this is not bundled, get it from graph (CRS provides one) - cluster_source, - cluster_sink - ); -#else - double flow = boost::kolmogorov_max_flow(graph, - boost::get(&EdgeP::edge_capacity, graph), - boost::get(&EdgeP::edge_residual_capacity, graph), - boost::get(&EdgeP::edge_reverse, graph), - boost::get(&VertexP::vertex_predecessor, graph), - boost::get(&VertexP::vertex_color, graph), - boost::get(&VertexP::vertex_distance_t, graph), - boost::get(boost::vertex_index, - graph), // this is not bundled, get it from graph - cluster_source, - cluster_sink - ); -#endif - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - cut_time += timer.time(); - #endif - - if(min_cut - flow <= flow * tolerance) { - continue; - } - min_cut = flow; - success = true; - //update labeling - for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) { - boost::default_color_type color = graph[vertex_i + 2].vertex_color; - if(labels[vertex_i] != alpha - && color == ColorTraits::white()) { //new comers (expansion occurs) - labels[vertex_i] = alpha; - } - } - } - } while(success); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << - std::endl; - CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; - CGAL_TRACE_STREAM << "graph creation time: " << graph_creation_time << - std::endl; - CGAL_TRACE_STREAM << "reverse mapping time: " << reverse_mapping_time << - std::endl; - CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; - #endif - return min_cut; - } }; -#endif #ifndef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE /** @@ -689,132 +482,71 @@ public: * For underlying max-flow algorithm, it uses the MAXFLOW software implemented by Boykov & Kolmogorov. * Also no pre-allocation is made. */ -class Alpha_expansion_graph_cut_boykov_kolmogorov +class Alpha_expansion_MaxFlow { public: + + typedef MaxFlow::Graph::node_id Vertex_descriptor; + +private: + + MaxFlow::Graph graph; + +public: + + void clear_graph() + { + graph = MaxFlow::Graph(); + } + + Vertex_descriptor add_vertex() + { + return graph.add_node(); + } + + void add_tweight (Vertex_descriptor& v, double w1, double w2) + { + graph.add_tweights(v, w1, w2); + } + + void init_vertices() + { + } + + double max_flow() + { + return graph.maxflow(); + } + + template + void update(VertexLabelMap vertex_label_map, + const std::vector& inserted_vertices, + InputVertexDescriptor vd, + std::size_t vertex_i, + std::size_t alpha) + { + if(get(vertex_label_map, vd) != alpha + && graph.what_segment(inserted_vertices[vertex_i]) == MaxFlow::Graph::SINK) + put(vertex_label_map, vd, alpha); + } + /** - * Applies alpha-expansion graph-cut for energy minimization. - * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1]) - * @param edge_weights contains weights for each edge in @a edges (correspondence according to order) - * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex] - * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]), - * and as output it returns final labeling of vertices - * @return result of energy function + * Adds two directional edges between @a v1 and @a v2 + * @param v1 first vertex + * @param v2 second vertex + * @param w1 weight for edge from v1 to v2 (v1->v2) + * @param w2 weight for edge from v2 to v1 (v2->v1) + * @param graph to be added + * @return pair of added edges, first: v1->v2 and second: v2->v1 */ - double operator()(const std::vector >& - edges, - const std::vector& edge_weights, - const std::vector >& probability_matrix, - std::vector& labels) const { - const double tolerance = 1e-10; - - double min_cut = (std::numeric_limits::max)(); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - double vertex_creation_time, edge_creation_time, cut_time; - vertex_creation_time = edge_creation_time = cut_time = 0.0; - #endif - - std::vector inserted_vertices; - inserted_vertices.resize(labels.size()); - bool success; - do { - success = false; - std::size_t alpha = 0; - for(std::vector >::const_iterator it = - probability_matrix.begin(); - it != probability_matrix.end(); ++it, ++alpha) { - MaxFlow::Graph graph; - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - Timer timer; - timer.start(); - #endif - // 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 < labels.size(); ++vertex_i) { - MaxFlow::Graph::node_id new_vertex = graph.add_node(); - inserted_vertices[vertex_i] = 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]; - graph.add_tweights(new_vertex, source_weight, sink_weight); - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - vertex_creation_time += timer.time(); - timer.reset(); - #endif - // 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) { - MaxFlow::Graph::node_id v1 = inserted_vertices[edge_it->first]; - MaxFlow::Graph::node_id v2 = inserted_vertices[edge_it->second]; - std::size_t label_1 = labels[edge_it->first], label_2 = labels[edge_it->second]; - if(label_1 == label_2) { - if(label_1 != alpha) { - graph.add_edge(v1, v2, *weight_it, *weight_it); - } - } else { - MaxFlow::Graph::node_id inbetween = graph.add_node(); - - double w1 = (label_1 == alpha) ? 0 : *weight_it; - double w2 = (label_2 == alpha) ? 0 : *weight_it; - graph.add_edge(inbetween, v1, w1, w1); - graph.add_edge(inbetween, v2, w2, w2); - - graph.add_tweights(inbetween, 0.0, *weight_it); - } - } - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - edge_creation_time += timer.time(); - timer.reset(); - #endif - - double flow = graph.maxflow(); - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - cut_time += timer.time(); - #endif - - if(min_cut - flow <= flow * tolerance) { - continue; - } - - min_cut = flow; - success = true; - //update labeling - for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) { - if(labels[vertex_i] != alpha - && graph.what_segment(inserted_vertices[vertex_i]) == MaxFlow::Graph::SINK) { - labels[vertex_i] = alpha; - } - } - } - } while(success); - - #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << - std::endl; - CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; - CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; - #endif - return min_cut; + void add_edge (Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, double w2) + { + graph.add_edge(v1, v2, w1, w2); } }; #endif //CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE -}//namespace internal /// @endcond -struct Alpha_expansion_boost_adjacency_list { }; -struct Alpha_expansion_boost_compressed_sparse_row { }; -struct Alpha_expansion_MaxFlow { }; - - template ::value) - return internal::Alpha_expansion_graph_cut_boost()(input_graph, - edge_weight_map, - vertex_index_map, - vertex_label_map, - vertex_label_probability_map); + typedef boost::graph_traits GT; + typedef typename GT::edge_descriptor input_edge_descriptor; + typedef typename GT::vertex_descriptor input_vertex_descriptor; + + typedef AlphaExpansionImplementation Alpha_expansion; + typedef typename Alpha_expansion::Vertex_descriptor Vertex_descriptor; + + // TODO: check this hardcoded parameter + const double tolerance = 1e-10; + + double min_cut = (std::numeric_limits::max)(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + double vertex_creation_time, edge_creation_time, cut_time; + vertex_creation_time = edge_creation_time = cut_time = 0.0; +#endif + + std::vector inserted_vertices; + inserted_vertices.resize(num_vertices (input_graph)); + + std::size_t number_of_labels = get(vertex_label_probability_map, *std::begin(vertices(input_graph))).size(); + + bool success; + do { + success = false; + + for (std::size_t alpha = 0; alpha < number_of_labels; ++ alpha) + { + alpha_expansion.clear_graph(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + Timer timer; + timer.start(); +#endif + + // For E-Data + // add every input vertex as a vertex to the graph, put edges to source & sink vertices + for (input_vertex_descriptor vd : vertices(input_graph)) + { + std::size_t vertex_i = get(vertex_index_map, vd); + Vertex_descriptor new_vertex = alpha_expansion.add_vertex(); + inserted_vertices[vertex_i] = new_vertex; + double source_weight = get(vertex_label_probability_map, vd)[alpha]; + // since it is expansion move, current alpha labeled vertices will be assigned to alpha again, + // making sink_weight 'infinity' guarantee this. + double sink_weight = (get(vertex_label_map, vd) == alpha ? + (std::numeric_limits::max)() + : get(vertex_label_probability_map, vd)[get(vertex_label_map, vd)]); + + alpha_expansion.add_tweight(new_vertex, source_weight, sink_weight); + } +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + vertex_creation_time += timer.time(); + timer.reset(); +#endif + + // For E-Smooth + // add edge between every vertex, + for (input_edge_descriptor ed : edges(input_graph)) + { + input_vertex_descriptor vd1 = source(ed, input_graph); + input_vertex_descriptor vd2 = target(ed, input_graph); + std::size_t idx1 = get (vertex_index_map, vd1); + std::size_t idx2 = get (vertex_index_map, vd2); + + double weight = get (edge_weight_map, ed); + + Vertex_descriptor v1 = inserted_vertices[idx1], + v2 = inserted_vertices[idx2]; + + std::size_t label_1 = get (vertex_label_map, vd1); + std::size_t label_2 = get (vertex_label_map, vd2); + if(label_1 == label_2) { + if(label_1 != alpha) { + alpha_expansion.add_edge(v1, v2, weight, weight); + } + } else { + Vertex_descriptor inbetween = alpha_expansion.add_vertex(); + + double w1 = (label_1 == alpha) ? 0 : weight; + double w2 = (label_2 == alpha) ? 0 : weight; + alpha_expansion.add_edge(inbetween, v1, w1, w1); + alpha_expansion.add_edge(inbetween, v2, w2, w2); + alpha_expansion.add_tweight(inbetween, 0., weight); + } + } +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + edge_creation_time += timer.time(); +#endif + + alpha_expansion.init_vertices(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + timer.reset(); +#endif + + double flow = alpha_expansion.max_flow(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + cut_time += timer.time(); +#endif + + if(min_cut - flow <= flow * tolerance) { + continue; + } + min_cut = flow; + success = true; + //update labeling + for (input_vertex_descriptor vd : vertices (input_graph)) + { + std::size_t vertex_i = get (vertex_index_map, vd); + alpha_expansion.update(vertex_label_map, inserted_vertices, vd, vertex_i, alpha); + } + } + } while(success); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << + std::endl; + CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; + CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; +#endif + + return min_cut; } template >& edges, + const std::vector& edge_weights, + const std::vector >& probability_matrix, + std::vector& labels) +{ + internal::Alpha_expansion_old_API_wrapper_graph graph (edges, edge_weights, probability_matrix, labels); + + return alpha_expansion_graph_cut(graph, + graph.edge_weight_map(), + graph.vertex_index_map(), + graph.vertex_label_map(), + graph.vertex_label_probability_map()); +} + }//namespace CGAL #endif //CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H