SMIL  1.0.4
MinSurfaces.hpp
1 /* __HEAD__
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * Copyright (c) 2017-2023, Centre de Morphologie Mathematique
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * * Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  * * Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in the
13  * documentation and/or other materials provided with the distribution.
14  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
15  * names of its contributors may be used to endorse or promote products
16  * derived from this software without specific prior written permission.
17  *
18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE
22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
28  * THE POSSIBILITY OF SUCH DAMAGE.
29  *
30  * Description :
31  * Porting GeoCuts module from Morph-M - This is the MinSurfaces part
32  *
33  * History :
34  * - 20/03/2019 - by Jose-Marcio Martins da Cruz
35  * Just created it...
36  * - XX/XX/XXXX -
37  *
38  *
39  * __HEAD__ - Stop here !
40  */
41 
42 #ifndef _D_GEOCUTS_MINSURFACE_HPP_
43 #define _D_GEOCUTS_MINSURFACE_HPP_
44 
45 #if 0
46 #include <morphee/selement/include/selementNeighborList.hpp>
47 #include <morphee/selement/include/private/selementNeighborhood_T.hpp>
48 #include <morphee/image/include/private/image_T.hpp>
49 #include <morphee/image/include/imageUtils.hpp>
50 #endif
51 
52 /*
53  *
54  *
55  *
56  */
57 #ifndef __BOOST_INCLUDED__
58 #define __BOOST_INCLUDED__
59 
60 #include <boost/config.hpp>
61 // for boost::tie
62 #include <boost/utility.hpp>
63 // for boost::graph_traits
64 #include <boost/graph/graph_traits.hpp>
65 #include <boost/graph/adjacency_list.hpp>
66 // JOE #include <boost/graph/graphviz.hpp>
67 // JOE #include <boost/graph/prim_minimum_spanning_tree.hpp>
68 // JOE #include <boost/graph/dijkstra_shortest_paths.hpp>
69 // JOE #include <boost/graph/johnson_all_pairs_shortest.hpp>
70 
71 #include <boost/version.hpp>
72 
73 #if 0
74 // FROM STAWIASKI JAN 2012
75 #include "../boost_ext/kolmogorov_max_flow_min_cost.hpp"
76 //#include "../boost_ext/maximum_spanning_tree.hpp"
77 //STAWIASKI JAN2012 commented, why?
78 //#include "../boost_ext/boost_compare.hpp"
79 #include <boost/graph/connected_components.hpp>
80 #endif
81 
82 #endif // __BOOST_INCLUDED__
83 
84 #include <vector>
85 
86 typedef struct {
87  float x;
88  float y;
89  float p;
90 } morceau;
91 
92 typedef std::list<morceau> affine_par_morceaux;
93 
94 using namespace geocuts;
95 
96 namespace smil
97 {
98  /*
99  *
100  *
101  *
102  */
103  template <class T1, class T2>
104  RES_T geoCuts(const Image<T1> &imIn, const Image<T1> &imGradx,
105  const Image<T1> &imGrady, const Image<T2> &imMarker,
106  const StrElt &nl, Image<T2> &imOut)
107  {
108  std::cout << "Enter function GeoCuts" << std::endl;
109 
110  ASSERT_ALLOCATED(&imIn, &imGradx, &imGrady, &imMarker, &imOut);
111  ASSERT_SAME_SIZE(&imIn, &imGradx, &imGrady, &imMarker, &imOut);
112 
113  typename Image<T1>::lineType bufIn = imIn.getPixels();
114  typename Image<T1>::lineType bufGradx = imGradx.getPixels();
115  typename Image<T1>::lineType bufGrady = imGrady.getPixels();
116  typename Image<T2>::lineType bufMarker = imMarker.getPixels();
117  typename Image<T2>::lineType bufOut = imOut.getPixels();
118 
119  // needed for max flow: capacit map, rev_capacity map, etc.
120  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
121  boost::directedS>
122  Traits_d;
123 
124  typedef boost::adjacency_list<
125  boost::listS, boost::vecS, boost::directedS,
126  boost::property<boost::vertex_name_t, std::string>,
127  boost::property<
128  boost::edge_capacity_t, double,
129  boost::property<boost::edge_residual_capacity_t, double,
130  boost::property<boost::edge_reverse_t,
131  Traits_d::edge_descriptor>>>>
132  Graph_d;
133 
134  Graph_d g;
135 
136  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
137  boost::get(boost::edge_capacity, g);
138 
139  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
140  get(boost::edge_reverse, g);
141 
142  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
143  residual_capacity = get(boost::edge_residual_capacity, g);
144 
145  Graph_d::edge_descriptor e1, e2, e3, e4;
146  Graph_d::vertex_descriptor vSource, vSink;
147  T1 numVertex = maxVal(imIn);
148  // T2 numLabels = maxVal(imMarker);
149 
150  std::cout << "build graph vertices" << std::endl;
151  std::cout << "number of vertices : " << numVertex << std::endl;
152 
153  size_t pixelCount = imIn.getPixelCount();
154  for (off_t i = 0; i < (off_t) pixelCount; i++) {
155  boost::add_vertex(g);
156  bufOut[i] = (T2) 1;
157  }
158 
159  vSource = boost::add_vertex(g);
160  vSink = boost::add_vertex(g);
161 
162  int width = imIn.getWidth();
163  int height = imIn.getHeight();
164  int depth = imIn.getDepth();
165 
166  off_t strideY = width;
167  off_t strideZ = width * height;
168  for (int z = 0; z < depth; z++) {
169  off_t p0Z = z * strideZ;
170  for (int y = 0; y < height; y++) {
171  off_t p0Y = y * strideY;
172  for (int x = 0; x < width; x++) {
173  off_t o1 = p0Z + p0Y + x;
174 
175  // iterators on the Structuring Element
176  vector<IntPoint> pts = filterStrElt(nl);
177  vector<IntPoint>::iterator itBegin, itEnd, it;
178  itBegin = pts.begin();
179  itEnd = pts.end();
180 
181  T1 val1 = bufIn[o1];
182  T2 marker = bufMarker[o1];
183 
184  bool hasEdge = false;
185  if (marker == 2) {
186  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o1, g);
187  boost::tie(e3, hasEdge) = boost::add_edge(o1, vSource, g);
188  capacity[e4] = (std::numeric_limits<double>::max)();
189  capacity[e3] = (std::numeric_limits<double>::max)();
190  rev[e4] = e3;
191  rev[e3] = e4;
192  } else if (marker == 3) {
193  boost::tie(e4, hasEdge) = boost::add_edge(o1, vSink, g);
194  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o1, g);
195  capacity[e4] = (std::numeric_limits<double>::max)();
196  capacity[e3] = (std::numeric_limits<double>::max)();
197  rev[e4] = e3;
198  rev[e3] = e4;
199  }
200 
201  for (it = itBegin; it != itEnd; it++) {
202  if (x + it->x > width - 1 || x + it->x < 0)
203  continue;
204  if (y + it->y > height - 1 || y + it->y < 0)
205  continue;
206  if (z + it->z > depth - 1 || z + it->z < 0)
207  continue;
208 
209  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
210  if (o2 <= o1)
211  continue;
212 
213  T2 val2 = bufIn[o2];
214  double cost = 1000 / (1 + 1.5 * (val1 - val2) * (val1 - val2));
215 
216  bool hasEdge = false;
217  boost::tie(e4, hasEdge) = boost::add_edge(o1, o2, g);
218  boost::tie(e3, hasEdge) = boost::add_edge(o2, o1, g);
219  capacity[e4] = cost;
220  capacity[e3] = cost;
221  rev[e4] = e3;
222  rev[e3] = e4;
223  }
224 
225  T1 valRight = 0;
226  T1 valLeft = 0;
227  T1 valUp = 0;
228  T1 valDown = 0;
229  for (it = itBegin; it != itEnd; it++) {
230  if (x + it->x > width - 1 || x + it->x < 0)
231  continue;
232  if (y + it->y > height - 1 || y + it->y < 0)
233  continue;
234  if (z + it->z > depth - 1 || z + it->z < 0)
235  continue;
236 
237  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
238  if (o2 == o1)
239  continue;
240 
241  // JOE - this is quite different from the original Stawiaski code
242  // JOE - doubt : should consider the Gradient on o2 or o1 ???
243  if (o2 > o1 + strideY) {
244  valDown = (T1) bufGrady[o2];
245  continue;
246  }
247  if (o2 < o1 + strideY) {
248  valUp = (T1) bufGrady[o2];
249  continue;
250  }
251  if (o2 > o1)
252  valRight = (T1) bufGradx[o2];
253  else
254  valLeft = (T1) bufGradx[o2];
255  }
256 
257  double divergence = 10. * ((valLeft - valRight) + (valUp - valDown));
258 
259  if (divergence < 0 && marker == 0) {
260  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o1, g);
261  boost::tie(e3, hasEdge) = boost::add_edge(o1, vSource, g);
262  capacity[e4] = std::abs(divergence);
263  capacity[e3] = std::abs(divergence);
264  rev[e4] = e3;
265  rev[e3] = e4;
266  } else if (divergence > 0 && marker == 0) {
267  boost::tie(e4, hasEdge) = boost::add_edge(o1, vSink, g);
268  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o1, g);
269  capacity[e4] = divergence;
270  capacity[e3] = divergence;
271  rev[e4] = e3;
272  rev[e3] = e4;
273  }
274  }
275  }
276  }
277 
278  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
279  boost::get(boost::vertex_index, g);
280  std::vector<boost::default_color_type> color(boost::num_vertices(g));
281 
282  std::cout << "Compute Max flow" << std::endl;
283 #if BOOST_VERSION >= 104700
284  double flow =
285  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
286  &color[0], indexmap, vSource, vSink);
287 #else
288  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
289  &color[0], indexmap, vSource, vSink);
290 #endif
291 
292  std::cout << "c The total flow:" << std::endl;
293  std::cout << "s " << flow << std::endl << std::endl;
294 
295  // for all pixels in imIn create a vertex and an edge
296  for (off_t o1 = 0; o1 < (off_t) pixelCount; o1++) {
297  if (color[o1] == color[vSource])
298  bufOut[o1] = (T2) 2;
299  else if (color[o1] == 1)
300  bufOut[o1] = (T2) 4;
301  else if (color[o1] == color[vSink])
302  bufOut[o1] = (T2) 3;
303  }
304 
305  return RES_OK;
306  }
307 
308  /*
309  *
310  *
311  *
312  */
313  template <class T1, class T2>
314  RES_T geoCutsMinSurfaces(const Image<T1> &imIn, const Image<T2> &imMarker,
315  const StrElt &nl, Image<T2> &imOut)
316  {
317  std::cout << "Enter function Geo-Cuts Watershed" << std::endl;
318 
319  ASSERT_ALLOCATED(&imIn, &imMarker, &imOut);
320  ASSERT_SAME_SIZE(&imIn, &imMarker, &imOut);
321 
322  typename Image<T1>::lineType bufIn = imIn.getPixels();
323  typename Image<T2>::lineType bufMarker = imMarker.getPixels();
324  typename Image<T2>::lineType bufOut = imOut.getPixels();
325 
326  // needed for max flow: capacit map, rev_capacity map, etc.
327  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
328  boost::directedS>
329  Traits_d;
330 
331  typedef boost::adjacency_list<
332  boost::vecS, boost::vecS, boost::directedS,
333  boost::property<boost::vertex_name_t, std::string>,
334  boost::property<
335  boost::edge_capacity_t, double,
336  boost::property<boost::edge_residual_capacity_t, double,
337  boost::property<boost::edge_reverse_t,
338  Traits_d::edge_descriptor>>>>
339  Graph_d;
340 
341  Graph_d g;
342 
343  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
344  boost::get(boost::edge_capacity, g);
345 
346  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
347  get(boost::edge_reverse, g);
348 
349  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
350  residual_capacity = get(boost::edge_residual_capacity, g);
351 
352  Graph_d::edge_descriptor e1, e2, e3, e4;
353  Graph_d::vertex_descriptor vSource, vSink;
354  T1 numVertex = maxVal(imIn);
355  // T2 numLabels = maxVal(imMarker);
356  T2 numEdges = 0;
357 
358  std::cout << "build graph vertices" << std::endl;
359  std::cout << "number of vertices : " << numVertex << std::endl;
360 
361  std::cout << "build graph vertices" << std::endl;
362  size_t pixelCount = imIn.getPixelCount();
363  for (off_t i = 0; i < (off_t) pixelCount; i++) {
364  boost::add_vertex(g);
365  bufOut[i] = (T2) 1;
366  }
367 
368  vSource = boost::add_vertex(g);
369  vSink = boost::add_vertex(g);
370 
371  // OK here
372 
373  int width = imIn.getWidth();
374  int height = imIn.getHeight();
375  int depth = imIn.getDepth();
376 
377  off_t strideY = width;
378  off_t strideZ = width * height;
379  for (int z = 0; z < depth; z++) {
380  off_t p0Z = z * strideZ;
381  for (int y = 0; y < height; y++) {
382  off_t p0Y = y * strideY;
383  for (int x = 0; x < width; x++) {
384  off_t o1 = p0Z + p0Y + x;
385 
386  // iterators on the Structuring Element
387  vector<IntPoint> pts = filterStrElt(nl);
388  vector<IntPoint>::iterator itBegin, itEnd, it;
389  itBegin = pts.begin();
390  itEnd = pts.end();
391 
392  T1 val1 = bufIn[o1];
393  T2 marker = bufMarker[o1];
394 
395  bool hasEdge = false;
396  if (marker == 2) {
397  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o1, g);
398  boost::tie(e3, hasEdge) = boost::add_edge(o1, vSource, g);
399  capacity[e4] = (std::numeric_limits<double>::max)();
400  capacity[e3] = (std::numeric_limits<double>::max)();
401  rev[e4] = e3;
402  rev[e3] = e4;
403  } else if (marker == 3) {
404  boost::tie(e4, hasEdge) = boost::add_edge(o1, vSink, g);
405  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o1, g);
406  capacity[e4] = (std::numeric_limits<double>::max)();
407  capacity[e3] = (std::numeric_limits<double>::max)();
408  rev[e4] = e3;
409  rev[e3] = e4;
410  }
411 
412  for (it = itBegin; it != itEnd; it++) {
413  if (x + it->x > width - 1 || x + it->x < 0)
414  continue;
415  if (y + it->y > height - 1 || y + it->y < 0)
416  continue;
417  if (z + it->z > depth - 1 || z + it->z < 0)
418  continue;
419 
420  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
421  if (o2 <= o1)
422  continue;
423 
424  T2 val2 = bufIn[o2];
425  double valeur = std::abs((double) (val2 - val1));
426  double cost = 256 / (1 + std::pow(valeur, 2));
427 
428  numEdges++;
429  bool hasEdge = false;
430  boost::tie(e4, hasEdge) = boost::add_edge(o1, o2, g);
431  boost::tie(e3, hasEdge) = boost::add_edge(o2, o1, g);
432  capacity[e4] = cost;
433  capacity[e3] = cost;
434  rev[e4] = e3;
435  rev[e3] = e4;
436  }
437  }
438  }
439  }
440  std::cout << "number of Edges : " << numEdges << std::endl;
441 
442  std::cout << "Compute Max flow" << std::endl;
443  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
444  boost::get(boost::vertex_index, g);
445  std::vector<boost::default_color_type> color(boost::num_vertices(g));
446 #if BOOST_VERSION >= 104700
447  double flow =
448  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
449  &color[0], indexmap, vSource, vSink);
450 #else
451  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
452  &color[0], indexmap, vSource, vSink);
453 #endif
454  std::cout << "c The total flow:" << std::endl;
455  std::cout << "s " << flow << std::endl << std::endl;
456 
457  // for all pixels in imIn create a vertex and an edge
458  for (off_t o1 = 0; o1 < (off_t) pixelCount; o1++) {
459  if (color[o1] == color[vSource])
460  bufOut[o1] = (T2) 2;
461  else if (color[o1] == color[vSink])
462  bufOut[o1] = (T2) 3;
463  else if (color[o1] == 1)
464  bufOut[o1] = (T2) 4;
465  }
466 
467  return RES_OK;
468  }
469 
470  /*
471  *
472  *
473  *
474  */
475  template <class T1, class T2>
476  RES_T geoCutsMultiway_MinSurfaces(const Image<T1> &imIn,
477  const Image<T2> &imMarker, const StrElt &nl,
478  Image<T2> &imOut)
479  {
480  ASSERT_ALLOCATED(&imIn, &imMarker, &imOut);
481  ASSERT_SAME_SIZE(&imIn, &imMarker, &imOut);
482 
483  typename Image<T1>::lineType bufIn = imIn.getPixels();
484  typename Image<T2>::lineType bufMarker = imMarker.getPixels();
485  typename Image<T2>::lineType bufOut = imOut.getPixels();
486 
487  // needed for max flow: capacit map, rev_capacity map, etc.
488  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
489  boost::directedS>
490  Traits_d;
491 
492  typedef boost::adjacency_list<
493  boost::listS, boost::vecS, boost::directedS,
494  boost::property<boost::vertex_name_t, std::string>,
495  boost::property<
496  boost::edge_capacity_t, double,
497  boost::property<boost::edge_residual_capacity_t, double,
498  boost::property<boost::edge_reverse_t,
499  Traits_d::edge_descriptor>>>>
500  Graph_d;
501 
502  Graph_d g;
503 
504  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
505  boost::get(boost::edge_capacity, g);
506 
507  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
508  get(boost::edge_reverse, g);
509 
510  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
511  residual_capacity = get(boost::edge_residual_capacity, g);
512 
513  Graph_d::edge_descriptor e1, e2, e3, e4;
514  Graph_d::vertex_descriptor vSource, vSink;
515 
516  std::cout << "build graph vertices" << std::endl;
517  size_t pixelCount = imIn.getPixelCount();
518 
519  size_t numVertex = pixelCount;
520  size_t numLabels = maxVal(imMarker);
521 
522  for (off_t i = 0; i < (off_t) pixelCount; i++) {
523  boost::add_vertex(g);
524  bufOut[i] = 1;
525  }
526 
527  std::cout << "number of Labels: " << numLabels << std::endl;
528  std::cout << "number of vertices: " << numVertex << std::endl;
529 
530  vSource = boost::add_vertex(g);
531  vSink = boost::add_vertex(g);
532 
533  std::cout << "build graph edges" << std::endl;
534 
535  int width = imIn.getWidth();
536  int height = imIn.getHeight();
537  int depth = imIn.getDepth();
538 
539  off_t strideY = width;
540  off_t strideZ = width * height;
541  for (int z = 0; z < depth; z++) {
542  off_t p0Z = z * strideZ;
543  for (int y = 0; y < height; y++) {
544  off_t p0Y = y * strideY;
545  for (int x = 0; x < width; x++) {
546  off_t o1 = p0Z + p0Y + x;
547 
548  // iterators on the Structuring Element
549  vector<IntPoint> pts = filterStrElt(nl);
550  vector<IntPoint>::iterator itBegin, itEnd, it;
551  itBegin = pts.begin();
552  itEnd = pts.end();
553 
554  T1 val1 = bufIn[o1];
555 
556  for (it = itBegin; it != itEnd; it++) {
557  if (x + it->x > width - 1 || x + it->x < 0)
558  continue;
559  if (y + it->y > height - 1 || y + it->y < 0)
560  continue;
561  if (z + it->z > depth - 1 || z + it->z < 0)
562  continue;
563 
564  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
565  if (o2 <= o1)
566  continue;
567 
568  T2 val2 = bufIn[o2];
569 
570  double cost = 10000.0 / (1.0 + 1.5 * (val1 - val2) * (val1 - val2));
571 
572  bool hasEdge = false;
573  boost::tie(e4, hasEdge) = boost::add_edge(o1, o2, g);
574  boost::tie(e3, hasEdge) = boost::add_edge(o2, o1, g);
575  capacity[e4] = cost + 1;
576  capacity[e3] = cost + 1;
577  rev[e4] = e3;
578  rev[e3] = e4;
579  }
580  }
581  }
582  }
583 
584  for (T2 nbk = 2; nbk <= numLabels; nbk++) {
585  // for all pixels in imIn create a vertex
586  for (off_t o0 = 0; o0 < (off_t) pixelCount; o0++) {
587  T2 val1 = bufMarker[o0];
588  T1 val2 = bufIn[o0];
589 
590  double cost = std::numeric_limits<double>::max();
591  bool hasEdge;
592  if (val1 == nbk) {
593  boost::tie(e4, hasEdge) = boost::edge(vSource, o0, g);
594  if (!hasEdge) {
595  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o0, g);
596  boost::tie(e3, hasEdge) = boost::add_edge(o0, vSource, g);
597  capacity[e4] = cost;
598  capacity[e3] = cost;
599  rev[e4] = e3;
600  rev[e3] = e4;
601  }
602  } else if (val1 > 1 && val1 != nbk) {
603  boost::tie(e4, hasEdge) = boost::edge(val2, vSink, g);
604  if (!hasEdge) {
605  boost::tie(e4, hasEdge) = boost::add_edge(o0, vSink, g);
606  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o0, g);
607  capacity[e4] = cost;
608  capacity[e3] = cost;
609  rev[e4] = e3;
610  rev[e3] = e4;
611  }
612  }
613  }
614 
615  std::cout << "Compute Max flow" << nbk << std::endl;
616  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
617  boost::get(boost::vertex_index, g);
618  std::vector<boost::default_color_type> color(boost::num_vertices(g));
619 
620 #if BOOST_VERSION >= 104700
621  double flow =
622  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
623  &color[0], indexmap, vSource, vSink);
624 #else
625  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
626  &color[0], indexmap, vSource, vSink);
627 #endif
628 
629  std::cout << "c The total flow:" << std::endl;
630  std::cout << "s " << flow << std::endl << std::endl;
631 
632  // for all pixels in imIn create a vertex and an edge
633  for (off_t o1 = 0; o1 < (off_t) pixelCount; o1++) {
634  T1 val1 = (int) bufIn[o1];
635  T2 val2 = (int) bufOut[o1];
636  T2 val3 = (int) bufMarker[o1];
637 
638  if (val2 == 1) {
639  if (color[val1] == color[vSource])
640  bufOut[o1] = (T2) nbk;
641  }
642 
643  bool hasEdge;
644  if (val3 == nbk) {
645  boost::tie(e4, hasEdge) = boost::edge(vSource, o1, g);
646  if (hasEdge) {
647  boost::remove_edge(vSource, o1, g);
648  boost::remove_edge(o1, vSource, g);
649  }
650  } else if (val3 > 1) {
651  boost::tie(e4, hasEdge) = boost::edge(o1, vSink, g);
652  if (hasEdge) {
653  boost::remove_edge(o1, vSink, g);
654  boost::remove_edge(vSink, o1, g);
655  }
656  }
657  }
658  }
659 
660  return RES_OK;
661  }
662 
663 } // namespace smil
664 
665 #endif // _D_GEOCUTS_MINSURFACE_HPP_
smil::BaseImage::getDepth
size_t getDepth() const
Get image depth (Z)
Definition: DBaseImage.h:90
smil::BaseImage::getWidth
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:80
smil::BaseImage::getPixelCount
size_t getPixelCount() const
Get the number of pixels.
Definition: DBaseImage.h:160
smil::BaseImage::getHeight
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:85
smil::Image::getPixels
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:105
smil::StrElt
Base structuring element.
Definition: DStructuringElement.h:68
smil::pow
RES_T pow(const Image< T1 > &imIn, Image< T2 > &imOut, double exponent=2)
pow() - power of an image
Definition: DImageArith.hpp:905
smil::maxVal
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
maxVal() - Max value of an image
Definition: DMeasures.hpp:348
smil::geoCuts
RES_T geoCuts(const Image< T1 > &imIn, const Image< T1 > &imGradx, const Image< T1 > &imGrady, const Image< T2 > &imMarker, const StrElt &nl, Image< T2 > &imOut)
Returns Geo Cuts algorithm.
Definition: MinSurfaces.hpp:104
smil::geoCutsMultiway_MinSurfaces
RES_T geoCutsMultiway_MinSurfaces(const Image< T1 > &imIn, const Image< T2 > &imMarker, const StrElt &nl, Image< T2 > &imOut)
Multiple object segmentation, Geo Cuts algorithm on a pixel adjacency graph, ImMarker is composed of ...
Definition: MinSurfaces.hpp:476
smil::geoCutsMinSurfaces
RES_T geoCutsMinSurfaces(const Image< T1 > &imIn, const Image< T2 > &imMarker, const StrElt &nl, Image< T2 > &imOut)
Geo Cuts algorithm on a pixel adjacency graph, ImMarker is composed of three values 0 for unmarked pi...
Definition: MinSurfaces.hpp:314
morceau
Definition: GeoCutsAlgo_impl_T.hpp:57