~ubuntu-branches/debian/stretch/cgal/stretch

« back to all changes in this revision

Viewing changes to include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h

  • Committer: Package Import Robot
  • Author(s): Joachim Reichel
  • Date: 2014-04-05 10:56:43 UTC
  • mfrom: (1.2.4)
  • Revision ID: package-import@ubuntu.com-20140405105643-jgnrpu2thtx23zfs
Tags: 4.4-1
* New upstream release.
* Remove patches do-not-link-example-with-qt4-support-library.patch and
  fix_jet_fitting_3.patch (applied upstream).

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#ifndef CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H
 
2
// Copyright (c) 2014  GeometryFactory Sarl (France).
 
3
// All rights reserved.
 
4
//
 
5
// This file is part of CGAL (www.cgal.org).
 
6
// You can redistribute it and/or modify it under the terms of the GNU
 
7
// General Public License as published by the Free Software Foundation,
 
8
// either version 3 of the License, or (at your option) any later version.
 
9
//
 
10
// Licensees holding a valid commercial license may use this file in
 
11
// accordance with the commercial license agreement provided with the software.
 
12
//
 
13
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 
14
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 
15
 
 
16
 
 
17
#define CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H
 
18
 
 
19
/// @cond CGAL_DOCUMENT_INTERNAL
 
20
 
 
21
/**
 
22
 * @file Alpha_expansion_graph_cut.h
 
23
 * @brief This file contains 3 graph-cut algorithms, which can be used as a template parameter for CGAL::internal::Surface_mesh_segmentation.
 
24
 *
 
25
 * Main differences between implementations are underlying max-flow algorithm and graph type (i.e. results are the same, performance differs).
 
26
 *
 
27
 * By default, we use MAXFLOW and the class Alpha_expansion_graph_cut_boykov_kolmogorov.
 
28
 * For deactivating MAXFLOW software and using boost implementation instead, define CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE.
 
29
 * It deactivates Alpha_expansion_graph_cut_boykov_kolmogorov, activate boost versions
 
30
 * and makes CGAL::internal::Surface_mesh_segmentation using Alpha_expansion_graph_cut_boost
 
31
 * as default implementation for the graph-cut.
 
32
 *
 
33
 * Also algorithms can be used by their-own for applying alpha-expansion graph-cut on any graph.
 
34
 *
 
35
 */
 
36
#include <CGAL/assertions.h>
 
37
#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
38
#include <CGAL/Timer.h>
 
39
#endif
 
40
#include <CGAL/trace.h>
 
41
 
 
42
#include <boost/version.hpp>
 
43
#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
 
44
#include <boost/graph/adjacency_list.hpp>
 
45
#include <boost/graph/compressed_sparse_row_graph.hpp>
 
46
#if BOOST_VERSION >= 104400 // at this version kolmogorov_max_flow become depricated.
 
47
#include <boost/graph/boykov_kolmogorov_max_flow.hpp>
 
48
#else
 
49
#include <boost/graph/kolmogorov_max_flow.hpp>
 
50
#endif
 
51
#else
 
52
namespace MaxFlow
 
53
{
 
54
#include <CGAL/internal/auxiliary/graph.h>
 
55
}
 
56
#endif
 
57
 
 
58
#include <vector>
 
59
 
 
60
 
 
61
 
 
62
 
 
63
namespace CGAL
 
64
{
 
65
namespace internal
 
66
{
 
67
 
 
68
////////////////////////////////////////////////////////////////////////////////////////
 
69
//   Comments about performance:
 
70
//
 
71
// 1) With BGL:
 
72
//     * Using adjacency_list:
 
73
//     ** Without pre-allocating vertex-list
 
74
//       | OutEdgeList | VertexList | Performance |
 
75
//       |    listS    |   listS    |     25.2    |
 
76
//       |    vecS     |   listS    |     22.7    |
 
77
//       |    listS    |   vecS     |     30.7    |
 
78
//       |    vecS     |   vecS     |     26.1    |
 
79
//
 
80
//     ** With pre-allocating vertex-list with max-node size
 
81
//        (Note: exact number of vertices are not certain at the beginning)
 
82
//       | OutEdgeList | VertexList | Performance |
 
83
//       |    listS    |   vecS     |     25.2    |
 
84
//       |    vecS     |   vecS     |     23.4    |
 
85
//
 
86
//     * Didn't try adjacency_matrix since our graph is sparse
 
87
//     ( Also one can check BGL book, performance section )
 
88
//
 
89
//    Decision:
 
90
//     * Alpha_expansion_graph_cut_boost: use adjacency_list<vecS, listS> without
 
91
//       pre-allocating vertex-list.
 
92
//
 
93
// 2) With Boykov-Kolmogorov MAXFLOW software:
 
94
//   (http://pub.ist.ac.at/~vnk/software/maxflow-v2.21.src.tar.gz)
 
95
//                                  | Performance |
 
96
//                                  |     3.1     |
 
97
//     * Alpha_expansion_graph_cut_boykov_kolmogorov provides an implementation.
 
98
//       MAXFLOW does not provide any option for pre-allocation (It is possible with v_3.02 though).
 
99
//
 
100
// Typical Benchmark result provided by Ilker
 
101
//                                 | construction of vertices  |  construction of edges    | graph cut  | Total
 
102
//   -----------------------------------------------------------------------------------------------------------
 
103
//   boost with an adjacency list  |         1.53              | 1.51                      |  3.00      | 6.04
 
104
//   boost with CSR                | 0.11 (gather in a vector) | 0.15 (gather in a vector) |  2.67      | 2.93
 
105
//   MaxFlow                       |       0.042               | 0.076                     |  1.043     | 1.161
 
106
//
 
107
// The main issue for now with CSR is the construction of the opposite edge map that is too costly,
 
108
// since it is done by exploring all edges to find opposite
 
109
////////////////////////////////////////////////////////////////////////////////////////
 
110
 
 
111
#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
 
112
/**
 
113
 * @brief Implements alpha-expansion graph cut algorithm.
 
114
 *
 
115
 * For representing graph, it uses adjacency_list with OutEdgeList = vecS, VertexList = listS.
 
116
 * Also no pre-allocation is made for vertex-list.
 
117
 */
 
118
class Alpha_expansion_graph_cut_boost
 
119
{
 
120
private:
 
121
  typedef boost::adjacency_list_traits<boost::vecS, boost::listS, boost::directedS>
 
122
  Adjacency_list_traits;
 
123
 
 
124
  typedef boost::adjacency_list<boost::vecS, boost::listS, boost::directedS,
 
125
          // 4 vertex properties
 
126
          boost::property<boost::vertex_index_t, std::size_t,
 
127
          boost::property<boost::vertex_color_t, boost::default_color_type,
 
128
          boost::property<boost::vertex_distance_t, double,
 
129
          boost::property<boost::vertex_predecessor_t, Adjacency_list_traits::edge_descriptor >
 
130
          > > >,
 
131
          // 3 edge properties
 
132
          boost::property<boost::edge_capacity_t, double,
 
133
          boost::property<boost::edge_residual_capacity_t, double,
 
134
          boost::property<boost::edge_reverse_t, Adjacency_list_traits::edge_descriptor> >
 
135
          > > Graph;
 
136
 
 
137
  typedef boost::graph_traits<Graph> Traits;
 
138
  typedef boost::color_traits<boost::default_color_type> ColorTraits;
 
139
 
 
140
  typedef Traits::vertex_descriptor Vertex_descriptor;
 
141
  typedef Traits::vertex_iterator   Vertex_iterator;
 
142
  typedef Traits::edge_descriptor   Edge_descriptor;
 
143
  typedef Traits::edge_iterator     Edge_iterator;
 
144
 
 
145
  /**
 
146
   * Adds two directional edges between @a v1 and @a v2
 
147
   * @param v1 first vertex
 
148
   * @param v2 second vertex
 
149
   * @param w1 weight for edge from v1 to v2 (v1->v2)
 
150
   * @param w2 weight for edge from v2 to v1 (v2->v1)
 
151
   * @param graph to be added
 
152
   * @return pair of added edges, first: v1->v2 and second: v2->v1
 
153
   */
 
154
  boost::tuple<Edge_descriptor, Edge_descriptor>
 
155
  add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1,
 
156
                       double w2, Graph& graph) const {
 
157
    Edge_descriptor v1_v2, v2_v1;
 
158
    bool v1_v2_added, v2_v1_added;
 
159
 
 
160
    boost::tie(v1_v2, v1_v2_added) = boost::add_edge(v1, v2, graph);
 
161
    boost::tie(v2_v1, v2_v1_added) = boost::add_edge(v2, v1, graph);
 
162
 
 
163
    CGAL_assertion(v1_v2_added && v2_v1_added);
 
164
    //put edge capacities
 
165
    boost::put(boost::edge_reverse, graph, v1_v2, v2_v1);
 
166
    boost::put(boost::edge_reverse, graph, v2_v1, v1_v2);
 
167
 
 
168
    //map reverse edges
 
169
    boost::put(boost::edge_capacity, graph, v1_v2, w1);
 
170
    boost::put(boost::edge_capacity, graph, v2_v1, w2);
 
171
 
 
172
    return boost::make_tuple(v1_v2, v2_v1);
 
173
  }
 
174
 
 
175
public:
 
176
  /**
 
177
   * Applies alpha-expansion graph-cut for energy minimization.
 
178
   * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1])
 
179
   * @param edge_weights contains weights for each edge in @a edges (correspondence according to order)
 
180
   * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex]
 
181
   * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]),
 
182
   * and as output it returns final labeling of vertices (i.e. assigned cluster-id to each facet)
 
183
   * @return result of energy function
 
184
   */
 
185
  double operator()(const std::vector<std::pair<std::size_t, std::size_t> >&
 
186
                    edges,
 
187
                    const std::vector<double>& edge_weights,
 
188
                    const std::vector<std::vector<double> >& probability_matrix,
 
189
                    std::vector<std::size_t>& labels) const {
 
190
    const double tolerance = 1e-10;
 
191
 
 
192
    double min_cut = (std::numeric_limits<double>::max)();
 
193
 
 
194
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
195
    double vertex_creation_time, edge_creation_time, cut_time;
 
196
    vertex_creation_time = edge_creation_time = cut_time = 0.0;
 
197
    #endif
 
198
 
 
199
    std::vector<Vertex_descriptor> inserted_vertices;
 
200
    inserted_vertices.resize(labels.size());
 
201
 
 
202
    Graph graph;
 
203
 
 
204
    bool success;
 
205
    do {
 
206
      success = false;
 
207
      std::size_t alpha = 0;
 
208
 
 
209
      for(std::vector<std::vector<double> >::const_iterator it =
 
210
            probability_matrix.begin();
 
211
          it != probability_matrix.end(); ++it, ++alpha) {
 
212
        graph.clear();
 
213
 
 
214
        Vertex_descriptor cluster_source = boost::add_vertex(graph);
 
215
        Vertex_descriptor cluster_sink = boost::add_vertex(graph);
 
216
 
 
217
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
218
        Timer timer;
 
219
        timer.start();
 
220
        #endif
 
221
 
 
222
        // For E-Data
 
223
        // add every facet as a vertex to the graph, put edges to source & sink vertices
 
224
        for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) {
 
225
          Vertex_descriptor new_vertex = boost::add_vertex(graph);
 
226
          inserted_vertices[vertex_i] = new_vertex;
 
227
          double source_weight = probability_matrix[alpha][vertex_i];
 
228
          // since it is expansion move, current alpha labeled vertices will be assigned to alpha again,
 
229
          // making sink_weight 'infinity' guarantee this.
 
230
          double sink_weight = (labels[vertex_i] == alpha) ?
 
231
                               (std::numeric_limits<double>::max)()
 
232
                               : probability_matrix[labels[vertex_i]][vertex_i];
 
233
 
 
234
          add_edge_and_reverse(cluster_source, new_vertex, source_weight, 0.0, graph);
 
235
          add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, 0.0, graph);
 
236
        }
 
237
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
238
        vertex_creation_time += timer.time();
 
239
        timer.reset();
 
240
        #endif
 
241
 
 
242
        // For E-Smooth
 
243
        // add edge between every vertex,
 
244
        std::vector<double>::const_iterator weight_it = edge_weights.begin();
 
245
        for(std::vector<std::pair<std::size_t, std::size_t> >::const_iterator edge_it =
 
246
              edges.begin(); edge_it != edges.end();
 
247
            ++edge_it, ++weight_it) {
 
248
          Vertex_descriptor v1 = inserted_vertices[edge_it->first],
 
249
                            v2 = inserted_vertices[edge_it->second];
 
250
          std::size_t label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
 
251
          if(label_1 == label_2) {
 
252
            if(label_1 != alpha) {
 
253
              add_edge_and_reverse(v1, v2, *weight_it, *weight_it, graph);
 
254
            }
 
255
          } else {
 
256
            Vertex_descriptor inbetween = boost::add_vertex(graph);
 
257
 
 
258
            double w1 = (label_1 == alpha) ? 0 : *weight_it;
 
259
            double w2 = (label_2 == alpha) ? 0 : *weight_it;
 
260
            add_edge_and_reverse(inbetween, v1, w1, w1, graph);
 
261
            add_edge_and_reverse(inbetween, v2, w2, w2, graph);
 
262
            add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
 
263
          }
 
264
        }
 
265
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
266
        edge_creation_time += timer.time();
 
267
        #endif
 
268
 
 
269
        // initialize vertex indices, it is neccessary since we are using VertexList = listS
 
270
        Vertex_iterator v_begin, v_end;
 
271
        Traits::vertices_size_type index = 0;
 
272
        for(boost::tie(v_begin, v_end) = vertices(graph); v_begin != v_end; ++v_begin) {
 
273
          boost::put(boost::vertex_index, graph, *v_begin, index++);
 
274
        }
 
275
 
 
276
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
277
        timer.reset();
 
278
        #endif
 
279
#if BOOST_VERSION >= 104400
 
280
        double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source,
 
281
                      cluster_sink);
 
282
#else
 
283
        double flow = boost::kolmogorov_max_flow(graph, cluster_source, cluster_sink);
 
284
#endif
 
285
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
286
        cut_time += timer.time();
 
287
        #endif
 
288
 
 
289
        if(min_cut - flow < flow * tolerance) {
 
290
          continue;
 
291
        }
 
292
        min_cut = flow;
 
293
        success = true;
 
294
        //update labeling
 
295
        for(std::size_t vertex_i = 0; vertex_i < inserted_vertices.size(); ++vertex_i) {
 
296
          boost::default_color_type color = boost::get(boost::vertex_color, graph,
 
297
                                            inserted_vertices[vertex_i]);
 
298
          if(labels[vertex_i] != alpha
 
299
              && color == ColorTraits::white()) { //new comers (expansion occurs)
 
300
            labels[vertex_i] = alpha;
 
301
          }
 
302
        }
 
303
      }
 
304
    } while(success);
 
305
 
 
306
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
307
    CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time <<
 
308
                      std::endl;
 
309
    CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl;
 
310
    CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl;
 
311
    #endif
 
312
 
 
313
    return min_cut;
 
314
  }
 
315
};
 
316
 
 
317
// another implementation using compressed_sparse_row_graph
 
318
// for now there is a performance problem while setting reverse edges
 
319
// if that can be solved, it is faster than Alpha_expansion_graph_cut_boost
 
320
class Alpha_expansion_graph_cut_boost_CSR
 
321
{
 
322
private:
 
323
  // CSR only accepts bundled props
 
324
  struct VertexP {
 
325
    boost::default_color_type vertex_color;
 
326
    double vertex_distance_t;
 
327
    // ? do not now there is another way to take it, I think since edge_descriptor does not rely on properties
 
328
    // this should be fine...
 
329
    boost::compressed_sparse_row_graph<boost::directedS>::edge_descriptor
 
330
    vertex_predecessor;
 
331
  };
 
332
 
 
333
  struct EdgeP {
 
334
    double edge_capacity;
 
335
    double edge_residual_capacity;
 
336
    boost::compressed_sparse_row_graph<boost::directedS>::edge_descriptor
 
337
    edge_reverse;
 
338
  };
 
339
 
 
340
  typedef boost::compressed_sparse_row_graph<boost::directedS,
 
341
          VertexP, EdgeP> Graph;
 
342
 
 
343
  typedef boost::graph_traits<Graph> Traits;
 
344
  typedef boost::color_traits<boost::default_color_type> ColorTraits;
 
345
 
 
346
  typedef Traits::vertex_descriptor Vertex_descriptor;
 
347
  typedef Traits::vertex_iterator   Vertex_iterator;
 
348
  typedef Traits::edge_descriptor   Edge_descriptor;
 
349
  typedef Traits::edge_iterator     Edge_iterator;
 
350
 
 
351
  void
 
352
  add_edge_and_reverse(std::size_t v1 , std::size_t v2, double w1, double w2,
 
353
                       std::vector<std::pair<std::size_t, std::size_t> >& edge_map,
 
354
                       std::vector<EdgeP>& edge_weights) const {
 
355
    edge_map.push_back(std::make_pair(v1, v2));
 
356
    EdgeP p1;
 
357
    p1.edge_capacity = w1;
 
358
    edge_weights.push_back(p1);
 
359
 
 
360
    edge_map.push_back(std::make_pair(v2, v1));
 
361
    EdgeP p2;
 
362
    p2.edge_capacity = w2;
 
363
    edge_weights.push_back(p2);
 
364
  }
 
365
 
 
366
public:
 
367
  /**
 
368
   * Applies alpha-expansion graph-cut for energy minimization.
 
369
   * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1])
 
370
   * @param edge_weights contains weights for each edge in @a edges (correspondence according to order)
 
371
   * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex]
 
372
   * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]),
 
373
   * and as output it returns final labeling of vertices (i.e. assigned cluster-id to each facet)
 
374
   * @return result of energy function
 
375
   */
 
376
  double operator()(const std::vector<std::pair<std::size_t, std::size_t> >&
 
377
                    edges,
 
378
                    const std::vector<double>& edge_weights,
 
379
                    const std::vector<std::vector<double> >& probability_matrix,
 
380
                    std::vector<std::size_t>& labels) const {
 
381
    const double tolerance = 1e-10;
 
382
 
 
383
    double min_cut = (std::numeric_limits<double>::max)();
 
384
 
 
385
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
386
    double vertex_creation_time, edge_creation_time, graph_creation_time,
 
387
           reverse_mapping_time, cut_time;
 
388
    vertex_creation_time = edge_creation_time = graph_creation_time =
 
389
                             reverse_mapping_time = cut_time = 0.0;
 
390
    #endif
 
391
 
 
392
    Graph graph;
 
393
 
 
394
    bool success;
 
395
    do {
 
396
      success = false;
 
397
      std::size_t alpha = 0;
 
398
 
 
399
      for(std::vector<std::vector<double> >::const_iterator it =
 
400
            probability_matrix.begin();
 
401
          it != probability_matrix.end(); ++it, ++alpha) {
 
402
        std::vector<std::pair<std::size_t, std::size_t> > edge_map;
 
403
        std::vector<EdgeP>                edge_map_weights;
 
404
        edge_map.reserve(labels.size() *
 
405
                         8); // there is no way to know exact edge count, it is a heuristic value
 
406
        edge_map_weights.reserve(labels.size() * 8);
 
407
 
 
408
        std::size_t cluster_source = 0;
 
409
        std::size_t cluster_sink = 1;
 
410
 
 
411
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
412
        Timer timer;
 
413
        timer.start();
 
414
        #endif
 
415
        // For E-Data
 
416
        // add every facet as a vertex to the graph, put edges to source & sink vertices
 
417
        for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) {
 
418
          double source_weight = probability_matrix[alpha][vertex_i];
 
419
          // since it is expansion move, current alpha labeled vertices will be assigned to alpha again,
 
420
          // making sink_weight 'infinity' guarantee this.
 
421
          double sink_weight = (labels[vertex_i] == alpha) ?
 
422
                               (std::numeric_limits<double>::max)()
 
423
                               : probability_matrix[labels[vertex_i]][vertex_i];
 
424
 
 
425
          add_edge_and_reverse(cluster_source, vertex_i + 2, source_weight, 0.0, edge_map,
 
426
                               edge_map_weights);
 
427
          add_edge_and_reverse(vertex_i + 2, cluster_sink, sink_weight, 0.0, edge_map,
 
428
                               edge_map_weights);
 
429
        }
 
430
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
431
        vertex_creation_time += timer.time();
 
432
        timer.reset();
 
433
        #endif
 
434
        // For E-Smooth
 
435
        // add edge between every vertex,
 
436
        std::size_t num_vert = labels.size() + 2;
 
437
        std::vector<double>::const_iterator weight_it = edge_weights.begin();
 
438
        for(std::vector<std::pair<std::size_t, std::size_t> >::const_iterator edge_it =
 
439
              edges.begin(); edge_it != edges.end();
 
440
            ++edge_it, ++weight_it) {
 
441
          std::size_t v1 = edge_it->first + 2, v2 = edge_it->second + 2;
 
442
          std::size_t label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
 
443
          if(label_1 == label_2) {
 
444
            if(label_1 != alpha) {
 
445
              add_edge_and_reverse(v1, v2, *weight_it, *weight_it, edge_map,
 
446
                                   edge_map_weights);
 
447
            }
 
448
          } else {
 
449
            std::size_t inbetween = num_vert++;
 
450
 
 
451
            double w1 = (label_1 == alpha) ? 0 : *weight_it;
 
452
            double w2 = (label_2 == alpha) ? 0 : *weight_it;
 
453
            add_edge_and_reverse(inbetween, v1, w1, w1, edge_map, edge_map_weights);
 
454
            add_edge_and_reverse(inbetween, v2, w2, w2, edge_map, edge_map_weights);
 
455
            add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, edge_map,
 
456
                                 edge_map_weights);
 
457
          }
 
458
        }
 
459
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
460
        edge_creation_time += timer.time();
 
461
        timer.reset();
 
462
        #endif
 
463
#if BOOST_VERSION >= 104000
 
464
        Graph graph(boost::edges_are_unsorted, edge_map.begin(), edge_map.end(),
 
465
                    edge_map_weights.begin(), num_vert);
 
466
#else
 
467
        Graph graph(edge_map.begin(), edge_map.end(),
 
468
                    edge_map_weights.begin(), num_vert);
 
469
#endif
 
470
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
471
        graph_creation_time += timer.time();
 
472
        timer.reset();
 
473
        #endif
 
474
 
 
475
        // PERFORMANCE PROBLEM
 
476
        // need to set reverse edge map, I guess there is no way to do that before creating the graph
 
477
        // since we do not have edge_descs
 
478
        // however from our edge_map, we know that each (2i, 2i + 1) is reverse pairs, how to facilitate that ?
 
479
        // will look it back
 
480
        Graph::edge_iterator ei, ee;
 
481
        for(boost::tie(ei, ee) = boost::edges(graph); ei != ee; ++ei) {
 
482
          Graph::vertex_descriptor v1 = boost::source(*ei, graph);
 
483
          Graph::vertex_descriptor v2 = boost::target(*ei, graph);
 
484
          std::pair<Graph::edge_descriptor, bool> opp_edge = boost::edge(v2, v1, graph);
 
485
 
 
486
          CGAL_assertion(opp_edge.second);
 
487
          graph[opp_edge.first].edge_reverse =
 
488
            *ei; // and edge_reverse of *ei will be (or already have been) set by the opp_edge
 
489
        }
 
490
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
491
        reverse_mapping_time += timer.time();
 
492
        timer.reset();
 
493
        #endif
 
494
 
 
495
#if BOOST_VERSION >= 104400
 
496
        // since properties are bundled, defaults does not work need to specify them
 
497
        double flow = boost::boykov_kolmogorov_max_flow(graph,
 
498
                      boost::get(&EdgeP::edge_capacity, graph),
 
499
                      boost::get(&EdgeP::edge_residual_capacity, graph),
 
500
                      boost::get(&EdgeP::edge_reverse, graph),
 
501
                      boost::get(&VertexP::vertex_predecessor, graph),
 
502
                      boost::get(&VertexP::vertex_color, graph),
 
503
                      boost::get(&VertexP::vertex_distance_t, graph),
 
504
                      boost::get(boost::vertex_index,
 
505
                                 graph), // this is not bundled, get it from graph (CRS provides one)
 
506
                      cluster_source,
 
507
                      cluster_sink
 
508
                                                       );
 
509
#else
 
510
        double flow = boost::kolmogorov_max_flow(graph,
 
511
                      boost::get(&EdgeP::edge_capacity, graph),
 
512
                      boost::get(&EdgeP::edge_residual_capacity, graph),
 
513
                      boost::get(&EdgeP::edge_reverse, graph),
 
514
                      boost::get(&VertexP::vertex_predecessor, graph),
 
515
                      boost::get(&VertexP::vertex_color, graph),
 
516
                      boost::get(&VertexP::vertex_distance_t, graph),
 
517
                      boost::get(boost::vertex_index,
 
518
                                 graph), // this is not bundled, get it from graph
 
519
                      cluster_source,
 
520
                      cluster_sink
 
521
                                                );
 
522
#endif
 
523
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
524
        cut_time += timer.time();
 
525
        #endif
 
526
 
 
527
        if(min_cut - flow < flow * tolerance) {
 
528
          continue;
 
529
        }
 
530
        min_cut = flow;
 
531
        success = true;
 
532
        //update labeling
 
533
        for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) {
 
534
          boost::default_color_type color =  graph[vertex_i + 2].vertex_color;
 
535
          if(labels[vertex_i] != alpha
 
536
              && color == ColorTraits::white()) { //new comers (expansion occurs)
 
537
            labels[vertex_i] = alpha;
 
538
          }
 
539
        }
 
540
      }
 
541
    } while(success);
 
542
 
 
543
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
544
    CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time <<
 
545
                      std::endl;
 
546
    CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl;
 
547
    CGAL_TRACE_STREAM << "graph creation time: " << graph_creation_time <<
 
548
                      std::endl;
 
549
    CGAL_TRACE_STREAM << "reverse mapping time: " << reverse_mapping_time <<
 
550
                      std::endl;
 
551
    CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl;
 
552
    #endif
 
553
    return min_cut;
 
554
  }
 
555
};
 
556
#endif
 
557
 
 
558
#ifndef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
 
559
/**
 
560
 * @brief Implements alpha-expansion graph cut algorithm.
 
561
 *
 
562
 * For underlying max-flow algorithm, it uses the MAXFLOW software implemented by Boykov & Kolmogorov.
 
563
 *  Also no pre-allocation is made.
 
564
 */
 
565
class Alpha_expansion_graph_cut_boykov_kolmogorov
 
566
{
 
567
public:
 
568
  /**
 
569
   * Applies alpha-expansion graph-cut for energy minimization.
 
570
   * @param edges contains incident vertex-id pairs for each edge (vertex-ids should be between [0, number of vertices -1])
 
571
   * @param edge_weights contains weights for each edge in @a edges (correspondence according to order)
 
572
   * @param probability_matrix contains responsibility of the center on the vertex probability[center][vertex]
 
573
   * @param[in, out] labels as input it contains initial labeling of vertices (i.e. a center-id between [0, number of centers -1]),
 
574
   * and as output it returns final labeling of vertices
 
575
   * @return result of energy function
 
576
   */
 
577
  double operator()(const std::vector<std::pair<std::size_t, std::size_t> >&
 
578
                    edges,
 
579
                    const std::vector<double>& edge_weights,
 
580
                    const std::vector<std::vector<double> >& probability_matrix,
 
581
                    std::vector<std::size_t>& labels) const {
 
582
    const double tolerance = 1e-10;
 
583
 
 
584
    double min_cut = (std::numeric_limits<double>::max)();
 
585
 
 
586
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
587
    double vertex_creation_time, edge_creation_time, cut_time;
 
588
    vertex_creation_time = edge_creation_time = cut_time = 0.0;
 
589
    #endif
 
590
 
 
591
    std::vector<MaxFlow::Graph::node_id> inserted_vertices;
 
592
    inserted_vertices.resize(labels.size());
 
593
    bool success;
 
594
    do {
 
595
      success = false;
 
596
      std::size_t alpha = 0;
 
597
      for(std::vector<std::vector<double> >::const_iterator it =
 
598
            probability_matrix.begin();
 
599
          it != probability_matrix.end(); ++it, ++alpha) {
 
600
        MaxFlow::Graph graph;
 
601
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
602
        Timer timer;
 
603
        timer.start();
 
604
        #endif
 
605
        // For E-Data
 
606
        // add every facet as a vertex to graph, put edges to source-sink vertices
 
607
        for(std::size_t vertex_i = 0; vertex_i <  labels.size(); ++vertex_i) {
 
608
          MaxFlow::Graph::node_id new_vertex = graph.add_node();
 
609
          inserted_vertices[vertex_i] = new_vertex;
 
610
 
 
611
          double source_weight = probability_matrix[alpha][vertex_i];
 
612
          // since it is expansion move, current alpha labeled vertices will be assigned to alpha again,
 
613
          // making sink_weight 'infinity' guarantee this.
 
614
          double sink_weight = (labels[vertex_i] == alpha) ?
 
615
                               (std::numeric_limits<double>::max)()
 
616
                               : probability_matrix[labels[vertex_i]][vertex_i];
 
617
          graph.add_tweights(new_vertex, source_weight, sink_weight);
 
618
        }
 
619
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
620
        vertex_creation_time += timer.time();
 
621
        timer.reset();
 
622
        #endif
 
623
        // For E-Smooth
 
624
        // add edge between every vertex,
 
625
        std::vector<double>::const_iterator weight_it = edge_weights.begin();
 
626
        for(std::vector<std::pair<std::size_t, std::size_t> >::const_iterator edge_it =
 
627
              edges.begin(); edge_it != edges.end();
 
628
            ++edge_it, ++weight_it) {
 
629
          MaxFlow::Graph::node_id v1 = inserted_vertices[edge_it->first];
 
630
          MaxFlow::Graph::node_id v2 = inserted_vertices[edge_it->second];
 
631
          std::size_t label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
 
632
          if(label_1 == label_2) {
 
633
            if(label_1 != alpha) {
 
634
              graph.add_edge(v1, v2, *weight_it, *weight_it);
 
635
            }
 
636
          } else {
 
637
            MaxFlow::Graph::node_id inbetween = graph.add_node();
 
638
 
 
639
            double w1 = (label_1 == alpha) ? 0 : *weight_it;
 
640
            double w2 = (label_2 == alpha) ? 0 : *weight_it;
 
641
            graph.add_edge(inbetween, v1, w1, w1);
 
642
            graph.add_edge(inbetween, v2, w2, w2);
 
643
 
 
644
            graph.add_tweights(inbetween, 0.0, *weight_it);
 
645
          }
 
646
        }
 
647
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
648
        edge_creation_time += timer.time();
 
649
        timer.reset();
 
650
        #endif
 
651
 
 
652
        double flow = graph.maxflow();
 
653
        #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
654
        cut_time += timer.time();
 
655
        #endif
 
656
 
 
657
        if(min_cut - flow < flow * tolerance) {
 
658
          continue;
 
659
        }
 
660
 
 
661
        min_cut = flow;
 
662
        success = true;
 
663
        //update labeling
 
664
        for(std::size_t vertex_i = 0; vertex_i < labels.size(); ++vertex_i) {
 
665
          if(labels[vertex_i] != alpha
 
666
              && graph.what_segment(inserted_vertices[vertex_i]) == MaxFlow::Graph::SINK) {
 
667
            labels[vertex_i] = alpha;
 
668
          }
 
669
        }
 
670
      }
 
671
    } while(success);
 
672
 
 
673
    #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
 
674
    CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time <<
 
675
                      std::endl;
 
676
    CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl;
 
677
    CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl;
 
678
    #endif
 
679
    return min_cut;
 
680
  }
 
681
};
 
682
#endif //CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
 
683
}//namespace internal
 
684
/// @endcond
 
685
}//namespace CGAL
 
686
#endif //CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H