// Copyright (c) 2014 GeometryFactory (France). All rights reserved. // // This file is part of CGAL (www.cgal.org) // // $URL: https://github.com/CGAL/cgal/blob/v5.2/BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h $ // $Id: alpha_expansion_graphcut.h e893ac1 2020-08-18T10:06:51+02:00 Sébastien Loriot // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Ilker O. Yaz, Simon Giraudot #ifndef CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H #define CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H #include #include #include #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT #include #endif #include #include #include #include #include #include #if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated. # include #else # include #endif #include namespace CGAL { /// \cond SKIP_IN_MANUAL namespace internal { struct Alpha_expansion_old_API_wrapper_graph { typedef std::size_t vertex_descriptor; typedef std::size_t edge_descriptor; typedef boost::directed_tag directed_category; typedef boost::disallow_parallel_edge_tag edge_parallel_category; typedef boost::edge_list_graph_tag traversal_category; typedef boost::counting_iterator counting_iterator; typedef CGAL::Iterator_range counting_range; typedef CGAL::Identity_property_map Vertex_index_map; typedef CGAL::Pointer_property_map::type Vertex_label_map; struct Vertex_label_cost_map { typedef std::size_t key_type; typedef std::vector value_type; typedef value_type reference; typedef boost::readable_property_map_tag category; const std::vector >* cost_matrix; Vertex_label_cost_map (const std::vector >* cost_matrix) : cost_matrix (cost_matrix) { } friend reference get (const Vertex_label_cost_map& pmap, key_type idx) { std::vector out; out.reserve (pmap.cost_matrix->size()); for (std::size_t i = 0; i < pmap.cost_matrix->size(); ++ i) out.push_back ((*pmap.cost_matrix)[i][idx]); return out; } }; typedef CGAL::Pointer_property_map::const_type Edge_cost_map; const std::vector >& edges; const std::vector& edge_costs; const std::vector >& cost_matrix; std::vector& labels; Alpha_expansion_old_API_wrapper_graph (const std::vector >& edges, const std::vector& edge_costs, const std::vector >& cost_matrix, std::vector& labels) : edges (edges), edge_costs (edge_costs), cost_matrix (cost_matrix), labels (labels) { } friend counting_range vertices (const Alpha_expansion_old_API_wrapper_graph& graph) { return CGAL::make_range (boost::counting_iterator(0), boost::counting_iterator(graph.labels.size())); } friend std::size_t num_vertices (const Alpha_expansion_old_API_wrapper_graph& graph) { return graph.labels.size(); } friend counting_range edges (const Alpha_expansion_old_API_wrapper_graph& graph) { return CGAL::make_range (boost::counting_iterator(0), boost::counting_iterator(graph.edges.size())); } friend vertex_descriptor source (edge_descriptor ed, const Alpha_expansion_old_API_wrapper_graph& graph) { return graph.edges[ed].first; } friend vertex_descriptor target (edge_descriptor ed, const Alpha_expansion_old_API_wrapper_graph& graph) { return graph.edges[ed].second; } Vertex_index_map vertex_index_map() const { return Vertex_index_map(); } Vertex_label_map vertex_label_map() { return CGAL::make_property_map(labels); } Vertex_label_cost_map vertex_label_cost_map() const { return Vertex_label_cost_map(&cost_matrix); } Edge_cost_map edge_cost_map() const { return CGAL::make_property_map(edge_costs); } }; //////////////////////////////////////////////////////////////////////////////////////// // Comments about performance: // // 1) With BGL: // * Using adjacency_list: // ** Without pre-allocating vertex-list // | OutEdgeList | VertexList | Performance | // | listS | listS | 25.2 | // | vecS | listS | 22.7 | // | listS | vecS | 30.7 | // | vecS | vecS | 26.1 | // // ** With pre-allocating vertex-list with max-node size // (Note: exact number of vertices are not certain at the beginning) // | OutEdgeList | VertexList | Performance | // | listS | vecS | 25.2 | // | vecS | vecS | 23.4 | // // * Didn't try adjacency_matrix since our graph is sparse // ( Also one can check BGL book, performance section ) // // Decision: // * Alpha_expansion_graph_cut_boost: use adjacency_list without // pre-allocating vertex-list. // // 2) With Boykov-Kolmogorov MAXFLOW software: // (http://pub.ist.ac.at/~vnk/software/maxflow-v2.21.src.tar.gz) // | Performance | // | 3.1 | // * Alpha_expansion_graph_cut_boykov_kolmogorov provides an implementation. // MAXFLOW does not provide any option for pre-allocation (It is possible with v_3.02 though). // // Typical Benchmark result provided by Ilker // | construction of vertices | construction of edges | graph cut | Total // ----------------------------------------------------------------------------------------------------------- // boost with an adjacency list | 1.53 | 1.51 | 3.00 | 6.04 // boost with CSR | 0.11 (gather in a vector) | 0.15 (gather in a vector) | 2.67 | 2.93 // MaxFlow | 0.042 | 0.076 | 1.043 | 1.161 // // The main issue for now with CSR is the construction of the opposite edge map that is too costly, // since it is done by exploring all edges to find opposite //////////////////////////////////////////////////////////////////////////////////////// } // 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_boost_adjacency_list_impl { private: typedef boost::adjacency_list_traits Adjacency_list_traits; typedef boost::adjacency_list > > >, // 3 edge properties boost::property > > > Graph; 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(std::size_t(get (vertex_label_map, vd)) != alpha && color == ColorTraits::white()) //new comers (expansion occurs) put (vertex_label_map, vd, static_cast::value_type>(alpha)); } 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; boost::tie(v1_v2, v1_v2_added) = boost::add_edge(v1, v2, graph); boost::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); } }; // 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_boost_compressed_sparse_row_impl { private: // CSR only accepts bundled props struct VertexP { boost::default_color_type vertex_color; double vertex_distance_t; // ? do not now there is another way to take it, I think since edge_descriptor does not rely on properties // this should be fine... boost::compressed_sparse_row_graph::edge_descriptor vertex_predecessor; }; struct EdgeP { double edge_capacity; double edge_residual_capacity; boost::compressed_sparse_row_graph::edge_descriptor edge_reverse; }; typedef boost::compressed_sparse_row_graph Graph; 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; 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&, 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_map_weights.push_back(p1); edge_map.push_back(std::make_pair(v2, v1)); EdgeP p2; p2.edge_capacity = w2; edge_map_weights.push_back(p2); } }; // tags struct Alpha_expansion_boost_adjacency_list_tag { }; struct Alpha_expansion_boost_compressed_sparse_row_tag { }; struct Alpha_expansion_MaxFlow_tag { }; // forward declaration class Alpha_expansion_MaxFlow_impl; /// \endcond // NOTE: latest performances check (2019-07-22) // // Using a random graph with 50000 vertices, 100000 edges and 30 labels: // // METHOD TIMING MEMORY // Boost Adjacency list 49s 122MiB // Boost CSR 187s 77MiB // MaxFlow 12s 717MiB /** \ingroup PkgBGLPartition regularizes a partition of a graph into `n` labels using the alpha expansion algorithm \cgalCite{Boykov2001FastApproximate}. For a graph \f$(V,E)\f$, this function computes a partition `f` that minimizes the following cost function: \f[ \mathrm{C}(f) = \sum_{\{v0,v1\} \in E} C_E(v0,v1) + \sum_{v \in V} C_V(f_v) \f] where \f$C_E(v0,v1)\f$ is the edge cost of assigning a different label to \f$v0\f$ and \f$v1\f$, and \f$C_V(f_v)\f$ is the vertex cost of assigning the label \f$f\f$ to the vertex \f$v\f$. \tparam InputGraph a model of `VertexAndEdgeListGraph` \tparam EdgeCostMap a model of `ReadablePropertyMap` with `boost::graph_traits::%edge_descriptor` as key and `double` as value \tparam VertexLabelCostMap a model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` as key and `std::vector` as value \tparam VertexLabelMap a model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` as key and `std::size_t` as value \tparam NamedParameters a sequence of named parameters \param input_graph the input graph. \param edge_cost_map a property map providing the weight of each edge. \param vertex_label_map a property map providing the label of each vertex. This map will be updated by the algorithm with the regularized version of the partition. \param vertex_label_cost_map a property map providing, for each vertex, an `std::vector` containing the cost of this vertex to belong to each label. Each `std::vector` should have the same size `n` (which is the number of labels), each label being indexed from `0` to `n-1`. For example, `get(vertex_label_cost_map, vd)[label_idx]` returns the cost of vertex `vd` to belong to the label `label_idx`. \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below \cgalNamedParamsBegin \cgalParamNBegin{vertex_index_map} \cgalParamDescription{a property map associating to each vertex of `input_graph` a unique index between `0` and `num_vertices(input_graph) - 1`} \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` as key type and `std::size_t` as value type} \cgalParamDefault{an automatically indexed internal map} \cgalParamExtra{If this parameter is not passed, internal machinery will create and initialize a face index property map, either using the internal property map if it exists or using an external map. The latter might result in - slightly - worsened performance in case of non-constant complexity for index access.} \cgalParamNEnd \cgalParamNBegin{implementation_tag} \cgalParamDescription{a tag used to select which implementation of the alpha expansion should be used. Available implementation tags are: - `CGAL::Alpha_expansion_boost_adjacency_list` - `CGAL::Alpha_expansion_boost_compressed_sparse_row_tag` - `CGAL::Alpha_expansion_MaxFlow_tag`} \cgalParamDefault{`CGAL::Alpha_expansion_boost_adjacency_list`} \cgalParamNEnd \cgalNamedParamsEnd \note The `MaxFlow` implementation is provided by the \ref PkgSurfaceMeshSegmentationRef under GPL license. The header `` must be included if users want to use this implementation. */ template double alpha_expansion_graphcut (const InputGraph& input_graph, EdgeCostMap edge_cost_map, VertexLabelCostMap vertex_label_cost_map, VertexLabelMap vertex_label_map, const NamedParameters& np) { using parameters::choose_parameter; using parameters::get_parameter; typedef boost::graph_traits GT; typedef typename GT::edge_descriptor input_edge_descriptor; typedef typename GT::vertex_descriptor input_vertex_descriptor; typedef typename GetInitializedVertexIndexMap::type VertexIndexMap; VertexIndexMap vertex_index_map = CGAL::get_initialized_vertex_index_map(input_graph, np); typedef typename GetImplementationTag::type Impl_tag; // select implementation typedef typename std::conditional ::value, Alpha_expansion_boost_adjacency_list_impl, typename std::conditional ::value, Alpha_expansion_boost_compressed_sparse_row_impl, Alpha_expansion_MaxFlow_impl>::type>::type Alpha_expansion; typedef typename Alpha_expansion::Vertex_descriptor Vertex_descriptor; Alpha_expansion alpha_expansion; // 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_cost_map, *(vertices(input_graph).first)).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 : CGAL::make_range(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_cost_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 = (std::size_t(get(vertex_label_map, vd)) == alpha ? (std::numeric_limits::max)() : get(vertex_label_cost_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 : CGAL::make_range(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_cost_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 : CGAL::make_range(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; } /// \cond SKIP_IN_MANUAL // variant with default NP template double alpha_expansion_graphcut (const InputGraph& input_graph, EdgeCostMap edge_cost_map, VertexLabelCostMap vertex_label_cost_map, VertexLabelMap vertex_label_map) { return alpha_expansion_graphcut (input_graph, edge_cost_map, vertex_label_cost_map, vertex_label_map, CGAL::parameters::all_default()); } // Old API inline double alpha_expansion_graphcut (const std::vector >& edges, const std::vector& edge_costs, const std::vector >& cost_matrix, std::vector& labels) { internal::Alpha_expansion_old_API_wrapper_graph graph (edges, edge_costs, cost_matrix, labels); return alpha_expansion_graphcut(graph, graph.edge_cost_map(), graph.vertex_label_cost_map(), graph.vertex_label_map(), CGAL::parameters::vertex_index_map (graph.vertex_index_map())); } template double alpha_expansion_graphcut (const std::vector >& edges, const std::vector& edge_costs, const std::vector >& cost_matrix, std::vector& labels, const AlphaExpansionImplementationTag&) { internal::Alpha_expansion_old_API_wrapper_graph graph (edges, edge_costs, cost_matrix, labels); return alpha_expansion_graphcut(graph, graph.edge_cost_map(), graph.vertex_label_cost_map(), graph.vertex_label_map(), CGAL::parameters::vertex_index_map (graph.vertex_index_map()). implementation_tag (AlphaExpansionImplementationTag())); } /// \endcond }//namespace CGAL namespace boost { template <> struct property_map { typedef CGAL::internal::Alpha_expansion_old_API_wrapper_graph::Vertex_index_map type; typedef CGAL::internal::Alpha_expansion_old_API_wrapper_graph::Vertex_index_map const_type; }; } #endif //CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H