Merge all three versions with one unique function + separate calls

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h

@@ -50,21 +50,23 @@
 #include <CGAL/boost/graph/named_function_params.h>
 
 #include <boost/version.hpp>
-#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
+
-#include <boost/graph/adjacency_list.hpp>
+#ifndef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
-#include <boost/graph/compressed_sparse_row_graph.hpp>
-#if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated.
-#include <boost/graph/boykov_kolmogorov_max_flow.hpp>
-#else
-#include <boost/graph/kolmogorov_max_flow.hpp>
-#endif
-#else
 namespace MaxFlow
 {
 #include <CGAL/internal/auxiliary/graph.h>
 }
 #endif
 
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/compressed_sparse_row_graph.hpp>
+
+#if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated.
+#  include <boost/graph/boykov_kolmogorov_max_flow.hpp>
+#else
+#  include <boost/graph/kolmogorov_max_flow.hpp>
+#endif
+
 #include <vector>
 
 
@@ -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<boost::vecS, boost::listS, boost::directedS>
@@ -225,11 +228,73 @@ private:
   typedef boost::graph_traits<Graph> Traits;
   typedef boost::color_traits<boost::default_color_type> 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 <typename VertexLabelMap, typename InputVertexDescriptor>
+  void update(VertexLabelMap vertex_label_map,
+              const std::vector<Vertex_descriptor>& 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<Edge_descriptor, Edge_descriptor>
+  void add_edge (Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, double w2)
-  add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1,
+  {
-                       double w2, Graph& graph) const {
     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<std::pair<std::size_t, std::size_t> >&
-                    edges,
-                    const std::vector<double>& edge_weights,
-                    const std::vector<std::vector<double> >& probability_matrix,
-                    std::vector<std::size_t>& 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 <typename InputGraph,
-            typename Edge_weight_map,
-            typename Vertex_index_map,
-            typename Vertex_label_map,
-            typename Vertex_label_probability_map>
-  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<InputGraph> 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<double>::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<Vertex_descriptor> 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<double>::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<Graph> Traits;
   typedef boost::color_traits<boost::default_color_type> 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
+private:
-  add_edge_and_reverse(std::size_t v1 , std::size_t v2, double w1, double w2,
+
-                       std::vector<std::pair<std::size_t, std::size_t> >& edge_map,
+  Graph graph;
-                       std::vector<EdgeP>& edge_weights) const {
+  std::size_t nb_vertices;
+  std::vector<std::pair<std::size_t, std::size_t> > edge_map;
+  std::vector<EdgeP>                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<Graph::edge_descriptor, bool> 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 <typename VertexLabelMap, typename InputVertexDescriptor>
+  void update(VertexLabelMap vertex_label_map,
+              const std::vector<Vertex_descriptor>& 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<std::pair<std::size_t, std::size_t> >&
-                    edges,
-                    const std::vector<double>& edge_weights,
-                    const std::vector<std::vector<double> >& probability_matrix,
-                    std::vector<std::size_t>& labels) const {
-    const double tolerance = 1e-10;
 
-    double min_cut = (std::numeric_limits<double>::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<std::vector<double> >::const_iterator it =
-            probability_matrix.begin();
-          it != probability_matrix.end(); ++it, ++alpha) {
-        std::vector<std::pair<std::size_t, std::size_t> > edge_map;
-        std::vector<EdgeP>                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<double>::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<double>::const_iterator weight_it = edge_weights.begin();
-        for(std::vector<std::pair<std::size_t, std::size_t> >::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<Graph::edge_descriptor, bool> 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 <typename VertexLabelMap, typename InputVertexDescriptor>
+  void update(VertexLabelMap vertex_label_map,
+              const std::vector<Vertex_descriptor>& 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.
+   * Adds two directional edges between @a v1 and @a v2
-   * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1])
+   * @param v1 first vertex
-   * @param edge_weights contains weights for each edge in @a edges (correspondence according to order)
+   * @param v2 second vertex
-   * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex]
+   * @param w1 weight for edge from v1 to v2 (v1->v2)
-   * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]),
+   * @param w2 weight for edge from v2 to v1 (v2->v1)
-   * and as output it returns final labeling of vertices
+   * @param graph to be added
-   * @return result of energy function
+   * @return pair of added edges, first: v1->v2 and second: v2->v1
    */
-  double operator()(const std::vector<std::pair<std::size_t, std::size_t> >&
+  void add_edge (Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, double w2)
-                    edges,
+  {
-                    const std::vector<double>& edge_weights,
+    graph.add_edge(v1, v2, w1, w2);
-                    const std::vector<std::vector<double> >& probability_matrix,
-                    std::vector<std::size_t>& labels) const {
-    const double tolerance = 1e-10;
-
-    double min_cut = (std::numeric_limits<double>::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<MaxFlow::Graph::node_id> inserted_vertices;
-    inserted_vertices.resize(labels.size());
-    bool success;
-    do {
-      success = false;
-      std::size_t alpha = 0;
-      for(std::vector<std::vector<double> >::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<double>::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<double>::const_iterator weight_it = edge_weights.begin();
-        for(std::vector<std::pair<std::size_t, std::size_t> >::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;
   }
 };
 #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 <typename InputGraph,
           typename Edge_weight_map,
           typename Vertex_index_map,
@@ -826,14 +558,132 @@ double alpha_expansion_graph_cut (const InputGraph& input_graph,
                                   Vertex_index_map vertex_index_map,
                                   Vertex_label_map vertex_label_map,
                                   Vertex_label_probability_map vertex_label_probability_map,
-                                  const AlphaExpansionImplementation& = AlphaExpansionImplementation())
+                                  AlphaExpansionImplementation alpha_expansion = AlphaExpansionImplementation())
 {
-  if (std::is_same<AlphaExpansionImplementation,  Alpha_expansion_boost_adjacency_list>::value)
+  typedef boost::graph_traits<InputGraph> GT;
-    return internal::Alpha_expansion_graph_cut_boost()(input_graph,
+  typedef typename GT::edge_descriptor input_edge_descriptor;
-                                                       edge_weight_map,
+  typedef typename GT::vertex_descriptor input_vertex_descriptor;
-                                                       vertex_index_map,
+
-                                                       vertex_label_map,
+  typedef AlphaExpansionImplementation Alpha_expansion;
-                                                       vertex_label_probability_map);
+  typedef typename Alpha_expansion::Vertex_descriptor Vertex_descriptor;
+
+  // TODO: check this hardcoded parameter
+  const double tolerance = 1e-10;
+
+  double min_cut = (std::numeric_limits<double>::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<Vertex_descriptor> 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<double>::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 <typename InputGraph,
@@ -856,5 +706,20 @@ double alpha_expansion_graph_cut (const InputGraph& input_graph,
     (input_graph, edge_weight_map, vertex_index_map, vertex_label_map, vertex_label_probability_map);
 }
 
+// Old API
+inline double alpha_expansion_graph_cut (const std::vector<std::pair<std::size_t, std::size_t> >& edges,
+                                         const std::vector<double>& edge_weights,
+                                         const std::vector<std::vector<double> >& probability_matrix,
+                                         std::vector<std::size_t>& 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