SMIL  1.0.4
Mosaic_GeoCuts_impl.full.hpp
1 #ifndef D_MOSAIC_GEOCUTS_IMPL_HPP
2 #define D_MOSAIC_GEOCUTS_IMPL_HPP
3 
4 #include <time.h>
5 
6 #include <boost/config.hpp>
7 // for boost::tie
8 #include <boost/utility.hpp>
9 // for boost::graph_traits
10 #include <boost/graph/graph_traits.hpp>
11 #include <boost/graph/adjacency_list.hpp>
12 
13 #include <boost/version.hpp>
14 #if BOOST_VERSION >= 104700
15 #include <boost/graph/boykov_kolmogorov_max_flow.hpp>
16 #elif BOOST_VERSION >= 103500
17 #include <boost/graph/kolmogorov_max_flow.hpp>
18 #endif
19 
20 #include <vector>
21 
22 static int debugOn = false;
23 
24 #define SMIL_ENTER_FUNCTION(a) \
25  { \
26  if (debugOn) \
27  cout << "Entering function " << __func__ << " " << (a) << endl; \
28  }
29 
30 #define SMIL_REGISTER_ERROR(a)
31 
32 typedef off_t offset_t;
33 
34 namespace smil
35 {
36  /*
37  * Some useful tools
38  */
39  /* This function selects points from an Structuring element
40  * generating a positive offset
41  */
42  inline vector<IntPoint> filterStrElt(StrElt se)
43  {
44  vector<IntPoint> pts(0);
45 
46  vector<IntPoint>::iterator it, itStart, itEnd;
47  itStart = se.points.begin();
48  itEnd = se.points.end();
49 
50  for (it = itStart; it != itEnd; it++) {
51  bool ok = (4 * it->z + 2 * it->y + it->x) > 0;
52  if (ok)
53  pts.push_back(*it);
54 
55  if (debugOn) {
56  std::cout << (ok ? "GOOD " : "BAD ") << std::right << " "
57  << std::setw(6) << it->x << " " << std::setw(6) << it->y
58  << " " << std::setw(6) << it->z << endl;
59  }
60  }
61  return pts;
62  }
63 
64  /*
65  *
66  */
67  using namespace boost;
68 
69  // needed for max flow: capacit map, rev_capacity map, etc.
70  typedef adjacency_list_traits<vecS, vecS, directedS> Traits_T;
71 
72  typedef adjacency_list<
73  vecS, vecS, directedS, property<vertex_name_t, std::string>,
74  property<edge_capacity_t, double,
75  property<edge_residual_capacity_t, double,
76  property<edge_reverse_t, Traits_T::edge_descriptor>>>>
77  Graph_T;
78 
79  // edge capacity
80  typedef property_map<Graph_T, edge_capacity_t>::type EdgeCap_T;
81  // edge reverse
82  typedef property_map<Graph_T, edge_reverse_t>::type EdgeRevCap_T;
83  // edge residual capacity
84  typedef property_map<Graph_T, edge_residual_capacity_t>::type EdgeResCap_T;
85  //
86  typedef property_map<Graph_T, vertex_index_t>::type VertexIndex_T;
87 
88  /*
89  *
90  *
91  *
92  *
93  */
94 #if 1
95  template <class T>
96  RES_T geoCutsMinSurfaces(const Image<T> &imIn, const Image<T> &imGrad,
97  const Image<T> &imMarker, const StrElt &nl,
98  Image<T> &imOut)
99  {
100  SMIL_ENTER_FUNCTION("");
101 
102  ASSERT_ALLOCATED(&imIn, &imGrad, &imMarker, &imOut);
103  ASSERT_SAME_SIZE(&imIn, &imGrad, &imMarker, &imOut);
104 
105  offset_t o0;
106  offset_t o1;
107 
108  Graph_T g;
109  EdgeCap_T capacity = boost::get(boost::edge_capacity, g);
110  EdgeRevCap_T rev = boost::get(boost::edge_reverse, g);
111  EdgeResCap_T residual_capacity =
112  boost::get(boost::edge_residual_capacity, g);
113 
114  bool in1;
115  Graph_T::edge_descriptor e1, e2, e3, e4, e5;
116  Graph_T::vertex_descriptor vSource, vSink;
117  int numVert = 0;
118  int numEdges = 0;
119  int max, not_used;
120 
121  std::cout << "build Region Adjacency Graph" << std::endl;
122  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
123 
124  clock_t t1 = clock();
125 
126  max = maxVal(imIn);
127 
128  numVert = max;
129  std::cout << "number of Vertices : " << numVert << std::endl;
130 
131  for (int i = 1; i <= numVert; i++) {
132  boost::add_vertex(g);
133  }
134 
135  vSource = boost::add_vertex(g);
136  vSink = boost::add_vertex(g);
137 
138  clock_t tt_marker2 = 0, tt_marker3 = 0, tt_new_edge = 0, tt_old_edge = 0;
139  clock_t t2 = clock();
140  std::cout << "Nodes creation time : " << double(t2 - t1) / CLOCKS_PER_SEC
141  << " seconds\n";
142 
143  std::cout << "Building Region Adjacency Graph Edges" << std::endl;
144  t1 = clock();
145 
146  typename Image<T>::lineType bufIn = imIn.getPixels();
147  typename Image<T>::lineType bufOut = imOut.getPixels();
148  typename Image<T>::lineType bufMarker = imMarker.getPixels();
149  typename Image<T>::lineType bufGrad = imGrad.getPixels();
150 
151  int pixelCount = imIn.getPixelCount();
152  for (o1 = 0; o1 < pixelCount; o1++) {
153  int val = (T) bufIn[o1];
154  int marker = (T) bufMarker[o1];
155  int val_prec = 0, marker_prec = 0;
156 
157  if (val <= 0)
158  continue;
159 
160  if (val > 0) {
161  if (marker == 2 && marker_prec != marker && val_prec != val) {
162  clock_t temps_marker2 = clock();
163  boost::tie(e4, in1) = boost::edge(vSource, val, g);
164 
165  if (in1 == 0) {
166  // std::cout<<"Add new edge marker 2"<<std::endl;
167  boost::tie(e4, in1) = boost::add_edge(vSource, val, g);
168  boost::tie(e3, in1) = boost::add_edge(val, vSource, g);
169  capacity[e4] = (std::numeric_limits<double>::max)();
170  capacity[e3] = (std::numeric_limits<double>::max)();
171  rev[e4] = e3;
172  rev[e3] = e4;
173  }
174  tt_marker2 += clock() - temps_marker2;
175  } else {
176  if (marker == 3 && marker_prec != marker && val_prec != val) {
177  clock_t temps_marker3 = clock();
178  boost::tie(e3, in1) = boost::edge(vSink, val, g);
179  if (in1 == 0) {
180  // std::cout<<"Add new edge marker 3"<<std::endl;
181  boost::tie(e4, in1) = boost::add_edge(val, vSink, g);
182  boost::tie(e3, in1) = boost::add_edge(vSink, val, g);
183  capacity[e4] = (std::numeric_limits<double>::max)();
184  capacity[e3] = (std::numeric_limits<double>::max)();
185  rev[e4] = e3;
186  rev[e3] = e4;
187  }
188  tt_marker3 += clock() - temps_marker3;
189  }
190  }
191 
192  double val_grad_o1 = imGrad.getPixel(o1);
193  // val de val2 precedente
194  int val2_prec = val;
195 
196 #if 0
197  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
198  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
199 
200  neighb.setCenter(o1);
201 
202  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
203  const offset_t o2 = nit.getOffset();
204 
205  if (o2 <= o1)
206  continue;
207  if (o2 > o1) {
208  int val2 = imIn.pixelFromOffset(o2);
209  if (val != val2) {
210  double val_grad_o2 = imGrad.pixelFromOffset(o2);
211  double maxi = std::max(val_grad_o1, val_grad_o2);
212  double cost = 10000.0 / (1.0 + std::pow(maxi, 4));
213 
214  if (val2_prec == val2) {
215  // same val2 means same edge (thus, keep e3 and e4)
216  capacity[e4] = capacity[e4] + cost;
217  capacity[e3] = capacity[e3] + cost;
218  } else {
219  boost::tie(e5, in1) = boost::edge(val, val2, g);
220 
221  if (in1 == 0) {
222  clock_t temps_new_edge = clock();
223  // std::cout<<"Add new edge "<< val<<" --
224  // "<<val2<<std::endl;
225  numEdges++;
226  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
227  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
228  capacity[e4] = cost;
229  capacity[e3] = cost;
230  rev[e4] = e3;
231  rev[e3] = e4;
232  tt_new_edge += clock() - temps_new_edge;
233 
234  } else {
235  clock_t temps_old_edge = clock();
236  // std::cout<<"existing edge"<<std::endl;
237  boost::tie(e4, in1) = boost::edge(val, val2, g);
238  boost::tie(e3, in1) = boost::edge(val2, val, g);
239  capacity[e4] = capacity[e4] + cost;
240  capacity[e3] = capacity[e3] + cost;
241  tt_old_edge += clock() - temps_old_edge;
242  }
243  val2_prec = val2;
244  }
245  }
246  }
247  }
248  val_prec = val;
249  marker_prec = marker;
250 #endif
251  }
252  }
253 
254  std::cout << "Number of edges : " << numEdges << std::endl;
255  t2 = clock();
256  std::cout << "Edges creation time : " << double(t2 - t1) / CLOCKS_PER_SEC
257  << " seconds\n";
258  std::cout << "Marker2 : " << double(tt_marker2) / CLOCKS_PER_SEC
259  << " seconds\n";
260  std::cout << "Marker3 : " << double(tt_marker3) / CLOCKS_PER_SEC
261  << " seconds\n";
262  std::cout << "New edges : " << double(tt_new_edge) / CLOCKS_PER_SEC
263  << " seconds\n";
264  std::cout << "Old edges : " << double(tt_old_edge) / CLOCKS_PER_SEC
265  << " seconds\n";
266 
267  boost::property_map<Graph_T, boost::vertex_index_t>::type indexmap =
268  boost::get(boost::vertex_index, g);
269  std::vector<boost::default_color_type> color(boost::num_vertices(g));
270 
271  std::cout << "Compute Max flow" << std::endl;
272  t1 = clock();
273 #if BOOST_VERSION >= 104700
274  double flow =
275  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
276  &color[0], indexmap, vSource, vSink);
277 #else
278  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
279  &color[0], indexmap, vSource, vSink);
280 #endif
281  std::cout << "c The total flow:" << std::endl;
282  std::cout << "s " << flow << std::endl << std::endl;
283  t2 = clock();
284  std::cout << "Flow computation time : " << double(t2 - t1) / CLOCKS_PER_SEC
285  << " seconds\n";
286 
287  t1 = clock();
288  pixelCount = imIn.getPixelCount();
289  for (o1 = 0; o1 < pixelCount; o1++) {
290  int val = (T) bufIn[o1];
291 
292  if (val == 0) {
293  bufOut[o1] = 0;
294  } else {
295  if (color[val] == color[vSource])
296  bufOut[o1] = 2;
297  else if (color[val] == color[vSink])
298  bufOut[o1] = 3;
299  else
300  bufOut[o1] = 4;
301  }
302  }
303 
304  t2 = clock();
305  std::cout << "Computing imOut took : " << double(t2 - t1) / CLOCKS_PER_SEC
306  << " seconds\n";
307 
308  return RES_OK;
309  }
310 #endif
311 
312 #if 0
313  /*
314  *
315  *
316  *
317  *
318  */
319  template <class ImageIn, class ImageGrad, class ImageMarker, class SE,
320  class ImageOut>
321  RES_T geoCutsMinSurfaces_with_Line(
322  const ImageIn &imIn, const ImageGrad &imGrad, const ImageMarker &imMarker,
323  const SE &nl, ImageOut &imOut)
324  {
325  SMIL_ENTER_FUNCTION("");
326 
327  std::cout << "Enter function t_geoCutsMinSurfaces_With_Line" << std::endl;
328 
329  if ((!imOut.isAllocated())) {
330  SMIL_REGISTER_ERROR("Not allocated");
331  return RES_NOT_ALLOCATED;
332  }
333 
334  if ((!imIn.isAllocated())) {
335  SMIL_REGISTER_ERROR("Not allocated");
336  return RES_NOT_ALLOCATED;
337  }
338 
339  if ((!imGrad.isAllocated())) {
340  SMIL_REGISTER_ERROR("Not allocated");
341  return RES_NOT_ALLOCATED;
342  }
343 
344  if ((!imMarker.isAllocated())) {
345  SMIL_REGISTER_ERROR("Not allocated");
346  return RES_NOT_ALLOCATED;
347  }
348 
349  // common image iterator
350  typename ImageIn::const_iterator it, iend;
351  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
352  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
353  offset_t o0;
354  offset_t o1;
355 
356  // needed for max flow: capacit map, rev_capacity map, etc.
357  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
358  boost::directedS>
359  Traits;
360  typedef boost::adjacency_list<
361  boost::listS, boost::vecS, boost::directedS,
362  boost::property<boost::vertex_name_t, std::string>,
363  boost::property<
364  boost::edge_capacity_t, double,
365  boost::property<boost::edge_residual_capacity_t, double,
366  boost::property<boost::edge_reverse_t,
367  Traits::edge_descriptor>>>>
368  Graph_d;
369 
370  Graph_d g;
371 
372  double sigma = 1.0;
373 
374  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
375  boost::get(boost::edge_capacity, g);
376 
377  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
378  get(boost::edge_reverse, g);
379 
380  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
381  residual_capacity = get(boost::edge_residual_capacity, g);
382 
383  bool in1;
384  Graph_d::edge_descriptor e1, e2, e3, e4;
385  Graph_d::vertex_descriptor vSource, vSink;
386  int numVert = 0;
387 
388  std::cout << "build Region Adjacency Graph" << std::endl;
389 
390  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
391  // for all pixels in imIn create a vertex
392  o0 = it.getOffset();
393  int val = imIn.pixelFromOffset(o0);
394 
395  if (val > numVert) {
396  numVert = val;
397  }
398  }
399 
400  std::cout << "build Region Adjacency Graph Vertices :" << numVert
401  << std::endl;
402 
403  for (int i = 0; i <= numVert; i++) {
404  boost::add_vertex(g);Mosaic_GeoCuts.hpp
405  }
406 
407  vSource = boost::add_vertex(g);
408  vSink = boost::add_vertex(g);
409 
410  std::cout << "build Region Adjacency Graph Edges" << std::endl;
411 
412  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
413  o1 = it.getOffset();
414  // std::cout<<o1<<std::endl;
415 
416  int val = imIn.pixelFromOffset(o1);
417  int marker = imMarker.pixelFromOffset(o1);
418 
419  if (marker == 2) {
420  boost::tie(e4, in1) = boost::edge(vSource, val, g);
421  if (in1 == 0) {
422  // std::cout<<"Add new edge marker 2"<<std::endl;
423  boost::tie(e4, in1) = boost::add_edge(vSource, val, g);
424  boost::tie(e3, in1) = boost::add_edge(val, vSource, g);
425  capacity[e4] = (std::numeric_limits<double>::max)();
426  capacity[e3] = (std::numeric_limits<double>::max)();
427  rev[e4] = e3;
428  rev[e3] = e4;
429  }
430  } else if (marker == 3) {
431  boost::tie(e3, in1) = boost::edge(vSink, val, g);
432  if (in1 == 0) {
433  // std::cout<<"Add new edge marker 3"<<std::endl;
434  boost::tie(e4, in1) = boost::add_edge(val, vSink, g);
435  boost::tie(e3, in1) = boost::add_edge(vSink, val, g);
436  capacity[e4] = (std::numeric_limits<double>::max)();
437  capacity[e3] = (std::numeric_limits<double>::max)();
438  rev[e4] = e3;
439  rev[e3] = e4;
440  }
441  }
442 
443  neighb.setCenter(o1);
444 
445  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
446  const offset_t o2 = nit.getOffset();
447  int val2 = imIn.pixelFromOffset(o2);
448 
449  if (o2 > o1) {
450  if (val != val2) {
451  boost::tie(e3, in1) = boost::edge(val, val2, g);
452  // std::cout<<in1<<std::endl;
453  // std::cout<<"Compute Gradient"<<std::endl;
454  double val3 = imGrad.pixelFromOffset(o1);
455  double val4 = imGrad.pixelFromOffset(o2);
456  double maxi = std::max(val3, val4);
457  double cost = 10000000.0 / (1.0 + std::pow(maxi, 5));
458 
459  if (in1 == 0) {
460  // std::cout<<"Add new edge"<<std::endl;
461  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
462  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
463  capacity[e4] = cost;
464  capacity[e3] = cost;
465  rev[e4] = e3;
466  rev[e3] = e4;
467  } else {
468  // std::cout<<"existing edge"<<std::endl;
469  boost::tie(e4, in1) = boost::edge(val, val2, g);
470  boost::tie(e3, in1) = boost::edge(val2, val, g);
471  capacity[e4] = capacity[e4] + cost;
472  capacity[e3] = capacity[e3] + cost;
473  }
474  }
475  }
476  }
477  }
478 
479  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
480  boost::get(boost::vertex_index, g);
481  std::vector<boost::default_color_type> color(boost::num_vertices(g));
482 
483  std::cout << "Compute Max flow" << std::endl;
484 #if BOOST_VERSION >= 104700
485  double flow =
486  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
487  &color[0], indexmap, vSource, vSink);
488 #else
489  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
490  &color[0], indexmap, vSource, vSink);
491 #endif
492  std::cout << "c The total flow:" << std::endl;
493  std::cout << "s " << flow << std::endl << std::endl;
494 
495  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
496  // for all pixels in imIn create a vertex and an edge
497  o1 = it.getOffset();
498  int val = imIn.pixelFromOffset(o1);
499 
500  if (color[val] == color[vSource])
501  imOut.setPixel(o1, 2);
502  else if (color[val] == color[vSink])
503  imOut.setPixel(o1, 3);
504  else
505  imOut.setPixel(o1, 4);
506  }
507 
508  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
509  // for all pixels in imIn create a vertex and an edge
510  o1 = it.getOffset();
511  int valeur_c = imOut.pixelFromOffset(o1);
512  neighb.setCenter(o1);
513 
514  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
515  const offset_t o2 = nit.getOffset();
516  int valeur_n = imOut.pixelFromOffset(o2);
517 
518  if (valeur_n != valeur_c && valeur_n != 0 && valeur_c != 0) {
519  imOut.setPixel(o1, 0);
520  }
521  }
522  }
523 
524  return RES_OK;
525  }
526 
527  /*
528  *
529  *
530  *
531  *
532  */
533  // ImageLabel and ImageMarker should be unsigned integers
534  template <class ImageLabel, class ImageVal, class ImageMarker, class SE,
535  class ImageOut>
536  RES_T geoCutsMinSurfaces_with_steps(
537  const ImageLabel &imLabel, const ImageVal &imVal,
538  const ImageMarker &imMarker, const SE &nl, F_SIMPLE step_x,
539  F_SIMPLE step_y, F_SIMPLE step_z, ImageOut &imOut)
540  {
541  SMIL_ENTER_FUNCTION("(multi-valued version)");
542 
543  // std::cout << "Enter function t_geoCutsMinSurfaces_with_steps
544  // (multi_valued version)" << std::endl;
545 
546  if (!imOut.isAllocated() || !imLabel.isAllocated() ||
547  !imVal.isAllocated() || !imMarker.isAllocated()) {
548  SMIL_REGISTER_ERROR("Image not allocated");
549  return RES_NOT_ALLOCATED;
550  }
551 
552  // common image iterator
553  typename ImageLabel::const_iterator it, iend;
554  morphee::selement::Neighborhood<SE, ImageLabel> neighb(imLabel, nl);
555  typename morphee::selement::Neighborhood<SE, ImageLabel>::iterator nit,
556  nend;
557  offset_t o0;
558  offset_t o1;
559 
560  // needed for max flow: capacit map, rev_capacity map, etc.
561  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
562  boost::directedS>
563  Traits;
564  typedef boost::adjacency_list<
565  boost::vecS, boost::vecS, boost::directedS,
566  boost::property<boost::vertex_name_t, std::string>,
567  boost::property<
568  boost::edge_capacity_t, F_DOUBLE,
569  boost::property<boost::edge_residual_capacity_t, F_DOUBLE,
570  boost::property<boost::edge_reverse_t,
571  Traits::edge_descriptor>>>>
572  Graph_d;
573 
574  // if we had computed the number of vertices before, we could directly
575  // initialize g with it
576  Graph_d g;
577 
578  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
579  boost::get(boost::edge_capacity, g);
580  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
581  get(boost::edge_reverse, g);
582  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
583  residual_capacity = get(boost::edge_residual_capacity, g);
584 
585  bool in1;
586  Graph_d::edge_descriptor e1, e2, e3, e4, e5;
587  Graph_d::vertex_descriptor vSource, vSink;
588  int numVert = 0;
589  int numEdges = 0;
590 
591  // Vertices creation
592  // std::cout<<"build Region Adjacency Graph Vertices"<<std::endl;
593  clock_t t1 = clock();
594  int max, not_used;
595  morphee::stats::t_measMinMax(imLabel, not_used, max);
596  numVert = max;
597  // std::cout<<"number of Vertices (without source and sink):
598  // "<<numVert<<std::endl;
599  // Warning : numVert+1 nodes created, but node 0 is not used (in order to
600  // simplify correspondance between labels and nodes)
601  for (int i = 0; i <= numVert; i++) {
602  boost::add_vertex(g);
603  }
604  vSource = boost::add_vertex(g);
605  vSink = boost::add_vertex(g);
606 
607  clock_t tt_marker2 = 0, tt_marker3 = 0, tt_new_edge = 0, tt_old_edge = 0;
608  clock_t t2 = clock();
609  // std::cout << "Nodes creation time : " << F_DOUBLE(t2-t1) /
610  // CLOCKS_PER_SEC << " seconds\n" ;
611 
612  // Edges creation
613  // std::cout<<"Building Region Adjacency Graph Edges"<<std::endl;
614  t1 = clock();
615  for (it = imLabel.begin(), iend = imLabel.end(); it != iend; ++it) {
616  o1 = it.getOffset();
617  typename ImageLabel::value_type label = imLabel.pixelFromOffset(o1),
618  label_prec = 0;
619  typename ImageMarker::value_type marker = imMarker.pixelFromOffset(o1),
620  marker_prec = 0;
621 
622  if (label > 0) {
623  if (marker == 2 && marker_prec != marker &&
624  label_prec != label) // add edge to Source
625  {
626  clock_t temps_marker2 = clock();
627  boost::tie(e4, in1) = boost::edge(vSource, label, g);
628  if (in1 == 0) // if in1 == 0 : edge had not yet been added
629  {
630  boost::tie(e4, in1) = boost::add_edge(vSource, label, g);
631  boost::tie(e3, in1) = boost::add_edge(label, vSource, g);
632  capacity[e4] = (std::numeric_limits<F_DOUBLE>::max)();
633  capacity[e3] = (std::numeric_limits<F_DOUBLE>::max)();
634  rev[e4] = e3;
635  rev[e3] = e4;
636  }
637  tt_marker2 += clock() - temps_marker2;
638  } else if (marker == 3 && marker_prec != marker &&
639  label_prec != label) // add edge to Sink
640  {
641  clock_t temps_marker3 = clock();
642  boost::tie(e3, in1) = boost::edge(vSink, label, g);
643  if (in1 == 0) // if in1 == 0 : edge had not yet been added
644  {
645  // std::cout<<"Add new edge marker 3"<<std::endl;
646  boost::tie(e4, in1) = boost::add_edge(label, vSink, g);
647  boost::tie(e3, in1) = boost::add_edge(vSink, label, g);
648  capacity[e4] = (std::numeric_limits<F_DOUBLE>::max)();
649  capacity[e3] = (std::numeric_limits<F_DOUBLE>::max)();
650  rev[e4] = e3;
651  rev[e3] = e4;
652  }
653  tt_marker3 += clock() - temps_marker3;
654  }
655 
656  neighb.setCenter(o1);
657  typename ImageVal::value_type val_o1 = imVal.pixelFromOffset(o1);
658  typename ImageLabel::value_type label2_prec =
659  label; // val de label2 precedente
660 
661  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
662  const offset_t o2 = nit.getOffset();
663 
664  if (o2 > o1) {
665  typename ImageLabel::value_type label2 =
666  imLabel.pixelFromOffset(o2);
667  if (label != label2) {
668  typename ImageVal::value_type val_o2 = imVal.pixelFromOffset(o2);
669  // F_SIMPLE max_diff = 255; // EDF TODO : use data traits
670  F_SIMPLE diff = t_Distance_L1(val_o1, val_o2);
671  // F_SIMPLE leak = std::sqrt(max_diff - diff);
672 
673  INT8 dx, dy, dz;
674  dx = std::abs(it.getX() - nit.getX());
675  dy = std::abs(it.getY() - nit.getY());
676  dz = std::abs(it.getZ() - nit.getZ());
677 
678  // F_SIMPLE dist ;
679  // dist = std::sqrt(std::pow(step_x*dx,2) +
680  // std::pow(step_y*dy,2) + std::pow(step_z*dz,2));
681  // if(dist == 0)
682  // {
683  // std::cout <<
684  // "ERROR : Distance between pixels equal to zero! " <<
685  // " Setting it to 1.\n" ;
686  // dist = 1 ;
687  // }
688 
689  F_SIMPLE surf = 1; // TODO : this only works with 4-connexity
690  // (in 2d) or 6-connexity (in 3d)
691  if (dx == 0)
692  surf *= step_x;
693  if (dy == 0)
694  surf *= step_y;
695  if (dz == 0)
696  surf *= step_z;
697  if (surf == 0) {
698  // std::cout << "ERROR : Surface between pixels equal to zero!
699  // Setting it to 1.\n" ;
700  surf = 1;
701  }
702 
703  // if (o2%1000 == 0)
704  // {
705  // std::cout << " surf : " << surf ;
706  // }
707 
708  // F_SIMPLE grad = diff/dist ;
709  // F_SIMPLE weighted_surf = grad * surf ;
710 
711  // Cette fonction devrait être remplacée par une fonction
712  // paramètre
713  // F_DOUBLE cost = 10000.0/(1.0+std::pow(diff/dist,4));
714  // F_DOUBLE cost = dist/(1+diff);
715  // F_DOUBLE cost = leak * surf ;
716  F_DOUBLE cost = surf / (1 + diff);
717  // if (o2%1000 == 0)
718  // {
719  // std::cout << " cost : " << cost << "\n";
720  // }
721 
722  // std::cout << " dx: " << (double)dx <<
723  // " dy: " << (double)dy <<
724  // " dz: " << (double)dz <<
725  // " dist: " << (double)dist <<
726  // " surf: " << (double)surf <<
727  // " grad: " << (double)grad <<
728  // " w_s: " << (double)weighted_surf <<
729  // " cost: " << (double)cost << "\n";
730 
731  if (label2_prec == label2) // same label2 means same edge (thus,
732  // keep e3 and e4)
733  {
734  capacity[e4] = capacity[e4] + cost;
735  capacity[e3] = capacity[e3] + cost;
736  } else {
737  boost::tie(e5, in1) = boost::edge(label, label2, g);
738  if (in1 == 0) {
739  clock_t temps_new_edge = clock();
740  // std::cout<<"Add new edge "<< label<<" --
741  // "<<label2<<std::endl;
742  numEdges++;
743  boost::tie(e4, in1) = boost::add_edge(label, label2, g);
744  boost::tie(e3, in1) = boost::add_edge(label2, label, g);
745  capacity[e4] = cost;
746  capacity[e3] = cost;
747  rev[e4] = e3;
748  rev[e3] = e4;
749  tt_new_edge += clock() - temps_new_edge;
750 
751  } else {
752  clock_t temps_old_edge = clock();
753  // std::cout<<"existing edge"<<std::endl;
754  boost::tie(e4, in1) = boost::edge(label, label2, g);
755  boost::tie(e3, in1) = boost::edge(label2, label, g);
756  capacity[e4] = capacity[e4] + cost;
757  capacity[e3] = capacity[e3] + cost;
758  tt_old_edge += clock() - temps_old_edge;
759  }
760  label2_prec = label2;
761  }
762  }
763  }
764  }
765  label_prec = label;
766  marker_prec = marker;
767  }
768  }
769 
770  t2 = clock();
771  // std::cout << "Number of initial edges : " << numEdges << std::endl;
772  // std::cout << "Edges creation time : " << F_DOUBLE(t2-t1) /
773  // CLOCKS_PER_SEC << " seconds\n" ; std::cout << "Marker2 : " <<
774  // F_DOUBLE(tt_marker2) / CLOCKS_PER_SEC << " seconds\n"; std::cout
775  // << "Marker3 : " << F_DOUBLE(tt_marker3) / CLOCKS_PER_SEC << "
776  // seconds\n"; std::cout << "New edges : " << F_DOUBLE(tt_new_edge) /
777  // CLOCKS_PER_SEC << " seconds\n"; std::cout << "Old edges : " <<
778  // F_DOUBLE(tt_old_edge) / CLOCKS_PER_SEC << " seconds\n";
779 
780  // We should test that the same region node is not connected
781  // simultaneously to source and sink :
782  // * iterate on sink neighbors ;
783  // * if neigbhbor is linked to source, remove edges to source and sink (in
784  // fact, set its capacity to 0, to avoid the modification of the graph
785  // structure)
786  //
787  t1 = clock();
788  Graph_d::vertex_descriptor sink_neighbour;
789  typename boost::graph_traits<Graph_d>::adjacency_iterator ai, ai_end;
790  UINT32 rem_edges = 0;
791  for (boost::tie(ai, ai_end) = adjacent_vertices(vSink, g); ai != ai_end;
792  ++ai) {
793  sink_neighbour = *ai;
794  tie(e1, in1) = edge(sink_neighbour, vSource, g);
795  if (in1) {
796  // remove_edge(vSource, sink_neighbour, g);
797  // remove_edge(sink_neighbour, vSource, g);
798  capacity[e1] = 0;
799  capacity[rev[e1]] = 0;
800  rem_edges++;
801  tie(e2, in1) = edge(vSink, sink_neighbour, g);
802  capacity[e2] = 0;
803  capacity[rev[e2]] = 0;
804  rem_edges++;
805  }
806  }
807  t2 = clock();
808  // std::cout << "Graph post-processing : Removal of " << rem_edges << "
809  // edges in : " << F_DOUBLE(t2-t1) / CLOCKS_PER_SEC << " seconds\n" ;
810 
811  // Prepare to run the max-flow algorithm
812  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
813  boost::get(boost::vertex_index, g);
814  std::vector<boost::default_color_type> color(boost::num_vertices(g));
815  // std::cout << "Compute Max flow" << std::endl;
816  t1 = clock();
817 
818 #if BOOST_VERSION >= 104700
819  F_DOUBLE flow =
820  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
821  &color[0], indexmap, vSource, vSink);
822 #else
823  F_DOUBLE flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
824  &color[0], indexmap, vSource, vSink);
825 #endif
826  // std::cout << "c The total flow:" << std::endl;
827  // std::cout << "s " << flow << std::endl;
828  t2 = clock();
829  // std::cout << "Flow computation time : " << F_DOUBLE(t2-t1) /
830  // CLOCKS_PER_SEC << " seconds\n" ;
831 
832  t1 = clock();
833  for (it = imLabel.begin(), iend = imLabel.end(); it != iend; ++it) {
834  o1 = it.getOffset();
835  typename ImageLabel::value_type label = imLabel.pixelFromOffset(o1);
836 
837  if (label == 0) {
838  imOut.setPixel(o1, 0);
839  } else // if color[label] == 0 : source node ; else, sink node (accord
840  // to boost graph doc)
841  {
842  if (color[label] == color[vSource])
843  imOut.setPixel(o1, 2);
844  else
845  imOut.setPixel(o1, 3);
846  }
847  }
848  t2 = clock();
849  // std::cout << "Computing imOut took : " << F_DOUBLE(t2-t1) /
850  // CLOCKS_PER_SEC << " seconds\n" ;
851 
852  // std::cout << "\n";
853 
854  return RES_OK;
855  }
856 
857  /*
858  *
859  *
860  *
861  *
862  */
863  // ImageLabel and ImageMarker should be unsigned integers
864  template <class ImageLabel, class ImageVal, class ImageMarker, class SE,
865  class ImageOut>
866  RES_T geoCutsMinSurfaces_with_steps_vGradient(
867  const ImageLabel &imLabel, const ImageVal &imVal,
868  const ImageMarker &imMarker, const SE &nl, F_SIMPLE step_x,
869  F_SIMPLE step_y, F_SIMPLE step_z, ImageOut &imOut)
870  {
871  SMIL_ENTER_FUNCTION("(multi-valued version)");
872 
873  // std::cout << "Enter function t_geoCutsMinSurfaces_with_steps
874  // (multi_valued version)" << std::endl;
875 
876  if (!imOut.isAllocated() || !imLabel.isAllocated() ||
877  !imVal.isAllocated() || !imMarker.isAllocated()) {
878  SMIL_REGISTER_ERROR("Image not allocated");
879  return RES_NOT_ALLOCATED;
880  }
881 
882  // common image iterator
883  typename ImageLabel::const_iterator it, iend;
884  morphee::selement::Neighborhood<SE, ImageLabel> neighb(imLabel, nl);
885  typename morphee::selement::Neighborhood<SE, ImageLabel>::iterator nit,
886  nend;
887  offset_t o0;
888  offset_t o1;
889 
890  // needed for max flow: capacit map, rev_capacity map, etc.
891  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
892  boost::directedS>
893  Traits;
894  typedef boost::adjacency_list<
895  boost::vecS, boost::vecS, boost::directedS,
896  boost::property<boost::vertex_name_t, std::string>,
897  boost::property<
898  boost::edge_capacity_t, F_DOUBLE,
899  boost::property<boost::edge_residual_capacity_t, F_DOUBLE,
900  boost::property<boost::edge_reverse_t,
901  Traits::edge_descriptor>>>>
902  Graph_d;
903 
904  // if we had computed the number of vertices before, we could directly
905  // initialize g with it
906  Graph_d g;
907 
908  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
909  boost::get(boost::edge_capacity, g);
910  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
911  get(boost::edge_reverse, g);
912  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
913  residual_capacity = get(boost::edge_residual_capacity, g);
914 
915  bool in1;
916  Graph_d::edge_descriptor e1, e2, e3, e4, e5;
917  Graph_d::vertex_descriptor vSource, vSink;
918  int numVert = 0;
919  int numEdges = 0;
920 
921  // Vertices creation
922  // std::cout<<"build Region Adjacency Graph Vertices"<<std::endl;
923  clock_t t1 = clock();
924  int max, not_used;
925  morphee::stats::t_measMinMax(imLabel, not_used, max);
926  numVert = max;
927  // std::cout<<"number of Vertices (without source and sink):
928  // "<<numVert<<std::endl;
929  // Warning : numVert+1 nodes created, but node 0 is not used (in order to
930  // simplify correspondance between labels and nodes)
931  for (int i = 0; i <= numVert; i++) {
932  boost::add_vertex(g);
933  }
934  vSource = boost::add_vertex(g);
935  vSink = boost::add_vertex(g);
936 
937  clock_t tt_marker2 = 0, tt_marker3 = 0, tt_new_edge = 0, tt_old_edge = 0;
938  clock_t t2 = clock();
939  // std::cout << "Nodes creation time : " << F_DOUBLE(t2-t1) /
940  // CLOCKS_PER_SEC << " seconds\n" ;
941 
942  // Edges creation
943  // std::cout<<"Building Region Adjacency Graph Edges"<<std::endl;
944  t1 = clock();
945  for (it = imLabel.begin(), iend = imLabel.end(); it != iend; ++it) {
946  o1 = it.getOffset();
947  typename ImageLabel::value_type label = imLabel.pixelFromOffset(o1),
948  label_prec = 0;
949  typename ImageMarker::value_type marker = imMarker.pixelFromOffset(o1),
950  marker_prec = 0;
951 
952  if (label > 0) {
953  if (marker == 2 && marker_prec != marker &&
954  label_prec != label) // add edge to Source
955  {
956  clock_t temps_marker2 = clock();
957  boost::tie(e4, in1) = boost::edge(vSource, label, g);
958  if (in1 == 0) // if in1 == 0 : edge had not yet been added
959  {
960  boost::tie(e4, in1) = boost::add_edge(vSource, label, g);
961  boost::tie(e3, in1) = boost::add_edge(label, vSource, g);
962  capacity[e4] = (std::numeric_limits<F_DOUBLE>::max)();
963  capacity[e3] = (std::numeric_limits<F_DOUBLE>::max)();
964  rev[e4] = e3;
965  rev[e3] = e4;
966  }
967  tt_marker2 += clock() - temps_marker2;
968  } else if (marker == 3 && marker_prec != marker &&
969  label_prec != label) // add edge to Sink
970  {
971  clock_t temps_marker3 = clock();
972  boost::tie(e3, in1) = boost::edge(vSink, label, g);
973  if (in1 == 0) // if in1 == 0 : edge had not yet been added
974  {
975  // std::cout<<"Add new edge marker 3"<<std::endl;
976  boost::tie(e4, in1) = boost::add_edge(label, vSink, g);
977  boost::tie(e3, in1) = boost::add_edge(vSink, label, g);
978  capacity[e4] = (std::numeric_limits<F_DOUBLE>::max)();
979  capacity[e3] = (std::numeric_limits<F_DOUBLE>::max)();
980  rev[e4] = e3;
981  rev[e3] = e4;
982  }
983  tt_marker3 += clock() - temps_marker3;
984  }
985 
986  neighb.setCenter(o1);
987  typename ImageVal::value_type val_o1 = imVal.pixelFromOffset(o1);
988  typename ImageLabel::value_type label2_prec =
989  label; // val de label2 precedente
990 
991  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
992  const offset_t o2 = nit.getOffset();
993 
994  if (o2 > o1) {
995  typename ImageLabel::value_type label2 =
996  imLabel.pixelFromOffset(o2);
997  if (label != label2) {
998  typename ImageVal::value_type val_o2 = imVal.pixelFromOffset(o2);
999  F_SIMPLE diff = std::max(val_o1, val_o2);
1000 
1001  INT8 dx, dy, dz;
1002  dx = std::abs(it.getX() - nit.getX());
1003  dy = std::abs(it.getY() - nit.getY());
1004  dz = std::abs(it.getZ() - nit.getZ());
1005 
1006  // F_SIMPLE dist ;
1007  // dist = std::sqrt(std::pow(step_x*dx,2) +
1008  // std::pow(step_y*dy,2) + std::pow(step_z*dz,2)); if(dist == 0)
1009  // {
1010  // std::cout << "ERROR : Distance between pixels equal to
1011  // zero! Setting it to 1.\n" ; dist = 1 ;
1012  // }
1013 
1014  F_SIMPLE surf = 1; // TODO : this only works with 4-connexity
1015  // (in 2d) or 6-connexity (in 3d)
1016  if (dx == 0)
1017  surf *= step_x;
1018  if (dy == 0)
1019  surf *= step_y;
1020  if (dz == 0)
1021  surf *= step_z;
1022  if (surf == 0) {
1023  // std::cout << "ERROR : Surface between pixels equal to zero!
1024  // Setting it to 1.\n" ;
1025  surf = 1;
1026  }
1027 
1028  // if (o2%1000 == 0)
1029  // {
1030  // std::cout << " surf : " << surf ;
1031  // }
1032 
1033  // F_SIMPLE grad = diff/dist ;
1034  // F_SIMPLE weighted_surf = grad * surf ;
1035 
1036  // Cette fonction devrait être remplacée par une fonction
1037  // paramètre
1038  // F_DOUBLE cost = 10000.0/(1.0+std::pow(diff/dist,4));
1039  // F_DOUBLE cost = dist/(1+diff);
1040  // F_DOUBLE cost = leak * surf ;
1041  F_DOUBLE cost = surf / (1 + diff);
1042  // if (o2%1000 == 0)
1043  // {
1044  // std::cout << " cost : " << cost << "\n";
1045  // }
1046 
1047  // std::cout << "dx: " << (double)dx << " dy: " << (double)dy
1048  // << " dz: " << (double)dz << " dist: " << (double)dist << "
1049  // surf: " << (double)surf << " grad: " << (double)grad << "
1050  // w_s: " << (double)weighted_surf << " cost: " << (double)cost
1051  // << "\n";
1052 
1053  if (label2_prec == label2) // same label2 means same edge (thus,
1054  // keep e3 and e4)
1055  {
1056  capacity[e4] = capacity[e4] + cost;
1057  capacity[e3] = capacity[e3] + cost;
1058  } else {
1059  boost::tie(e5, in1) = boost::edge(label, label2, g);
1060  if (in1 == 0) {
1061  clock_t temps_new_edge = clock();
1062  // std::cout<<"Add new edge "<< label<<" --
1063  // "<<label2<<std::endl;
1064  numEdges++;
1065  boost::tie(e4, in1) = boost::add_edge(label, label2, g);
1066  boost::tie(e3, in1) = boost::add_edge(label2, label, g);
1067  capacity[e4] = cost;
1068  capacity[e3] = cost;
1069  rev[e4] = e3;
1070  rev[e3] = e4;
1071  tt_new_edge += clock() - temps_new_edge;
1072 
1073  } else {
1074  clock_t temps_old_edge = clock();
1075  // std::cout<<"existing edge"<<std::endl;
1076  boost::tie(e4, in1) = boost::edge(label, label2, g);
1077  boost::tie(e3, in1) = boost::edge(label2, label, g);
1078  capacity[e4] = capacity[e4] + cost;
1079  capacity[e3] = capacity[e3] + cost;
1080  tt_old_edge += clock() - temps_old_edge;
1081  }
1082  label2_prec = label2;
1083  }
1084  }
1085  }
1086  }
1087  label_prec = label;
1088  marker_prec = marker;
1089  }
1090  }
1091 
1092  t2 = clock();
1093  // std::cout << "Number of initial edges : " << numEdges << std::endl;
1094  // std::cout << "Edges creation time : " << F_DOUBLE(t2-t1) /
1095  // CLOCKS_PER_SEC << " seconds\n" ; std::cout << "Marker2 : " <<
1096  // F_DOUBLE(tt_marker2) / CLOCKS_PER_SEC << " seconds\n"; std::cout
1097  // << "Marker3 : " << F_DOUBLE(tt_marker3) / CLOCKS_PER_SEC << "
1098  // seconds\n"; std::cout << "New edges : " << F_DOUBLE(tt_new_edge) /
1099  // CLOCKS_PER_SEC << " seconds\n"; std::cout << "Old edges : " <<
1100  // F_DOUBLE(tt_old_edge) / CLOCKS_PER_SEC << " seconds\n";
1101 
1102  // We should test that the same region node is not connected
1103  // simultaneously to source and sink :
1104  // * iterate on sink neighbors ;
1105  // * if neigbhbor is linked to source, remove edges to source and sink (in
1106  // fact, set its capacity to 0, to avoid the modification of the graph
1107  // structure)
1108  //
1109  t1 = clock();
1110  Graph_d::vertex_descriptor sink_neighbour;
1111  typename boost::graph_traits<Graph_d>::adjacency_iterator ai, ai_end;
1112  UINT32 rem_edges = 0;
1113  for (boost::tie(ai, ai_end) = adjacent_vertices(vSink, g); ai != ai_end;
1114  ++ai) {
1115  sink_neighbour = *ai;
1116  tie(e1, in1) = edge(sink_neighbour, vSource, g);
1117  if (in1) {
1118  // remove_edge(vSource, sink_neighbour, g);
1119  // remove_edge(sink_neighbour, vSource, g);
1120  capacity[e1] = 0;
1121  capacity[rev[e1]] = 0;
1122  rem_edges++;
1123  tie(e2, in1) = edge(vSink, sink_neighbour, g);
1124  capacity[e2] = 0;
1125  capacity[rev[e2]] = 0;
1126  rem_edges++;
1127  }
1128  }
1129  t2 = clock();
1130  // std::cout << "Graph post-processing : Removal of " << rem_edges << "
1131  // edges in : " << F_DOUBLE(t2-t1) / CLOCKS_PER_SEC << " seconds\n" ;
1132 
1133  // Prepare to run the max-flow algorithm
1134  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
1135  boost::get(boost::vertex_index, g);
1136  std::vector<boost::default_color_type> color(boost::num_vertices(g));
1137  // std::cout << "Compute Max flow" << std::endl;
1138  t1 = clock();
1139 #if BOOST_VERSION >= 104700
1140  F_DOUBLE flow =
1141  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1142  &color[0], indexmap, vSource, vSink);
1143 #else
1144  F_DOUBLE flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1145  &color[0], indexmap, vSource, vSink);
1146 #endif
1147  // std::cout << "c The total flow:" << std::endl;
1148  // std::cout << "s " << flow << std::endl;
1149  t2 = clock();
1150  // std::cout << "Flow computation time : " << F_DOUBLE(t2-t1) /
1151  // CLOCKS_PER_SEC << " seconds\n" ;
1152 
1153  t1 = clock();
1154  for (it = imLabel.begin(), iend = imLabel.end(); it != iend; ++it) {
1155  o1 = it.getOffset();
1156  typename ImageLabel::value_type label = imLabel.pixelFromOffset(o1);
1157 
1158  if (label == 0) {
1159  imOut.setPixel(o1, 0);
1160  } else // if color[label] == 0 : source node ; else, sink node (accord
1161  // to boost graph doc)
1162  {
1163  if (color[label] == color[vSource])
1164  imOut.setPixel(o1, 2);
1165  else
1166  imOut.setPixel(o1, 3);
1167  }
1168  }
1169  t2 = clock();
1170  // std::cout << "Computing imOut took : " << F_DOUBLE(t2-t1) /
1171  // CLOCKS_PER_SEC << " seconds\n" ;
1172 
1173  // std::cout << "\n";
1174 
1175  return RES_OK;
1176  }
1177 
1178  /*
1179  *
1180  *
1181  *
1182  *
1183  */
1184  template <class ImageIn, class ImageGrad, class ImageMarker, class SE,
1185  class ImageOut>
1186  RES_T geoCutsMinSurfaces_with_steps_old(
1187  const ImageIn &imIn, const ImageGrad &imGrad, const ImageMarker &imMarker,
1188  const SE &nl, F_SIMPLE step_x, F_SIMPLE step_y, F_SIMPLE step_z,
1189  ImageOut &imOut)
1190  {
1191  SMIL_ENTER_FUNCTION("");
1192 
1193  std::cout << "Enter function t_geoCutsMinSurfaces" << std::endl;
1194 
1195  if ((!imOut.isAllocated())) {
1196  SMIL_REGISTER_ERROR("Not allocated");
1197  return RES_NOT_ALLOCATED;
1198  }
1199 
1200  if ((!imIn.isAllocated())) {
1201  SMIL_REGISTER_ERROR("Not allocated");
1202  return RES_NOT_ALLOCATED;
1203  }
1204 
1205  if ((!imGrad.isAllocated())) {
1206  SMIL_REGISTER_ERROR("Not allocated");
1207  return RES_NOT_ALLOCATED;
1208  }
1209 
1210  if ((!imMarker.isAllocated())) {
1211  SMIL_REGISTER_ERROR("Not allocated");
1212  return RES_NOT_ALLOCATED;
1213  }
1214 
1215  // common image iterator
1216  typename ImageIn::const_iterator it, iend;
1217  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
1218  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
1219  offset_t o0;
1220  offset_t o1;
1221 
1222  // needed for max flow: capacit map, rev_capacity map, etc.
1223  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
1224  boost::directedS>
1225  Traits;
1226  typedef boost::adjacency_list<
1227  boost::vecS, boost::vecS, boost::directedS,
1228  boost::property<boost::vertex_name_t, std::string>,
1229  boost::property<
1230  boost::edge_capacity_t, double,
1231  boost::property<boost::edge_residual_capacity_t, double,
1232  boost::property<boost::edge_reverse_t,
1233  Traits::edge_descriptor>>>>
1234  Graph_d;
1235 
1236  Graph_d g;
1237 
1238  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
1239  boost::get(boost::edge_capacity, g);
1240 
1241  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
1242  get(boost::edge_reverse, g);
1243 
1244  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
1245  residual_capacity = get(boost::edge_residual_capacity, g);
1246 
1247  bool in1;
1248  Graph_d::edge_descriptor e1, e2, e3, e4, e5;
1249  Graph_d::vertex_descriptor vSource, vSink;
1250  int numVert = 0;
1251  int numEdges = 0;
1252  int max, not_used;
1253 
1254  std::cout << "build Region Adjacency Graph" << std::endl;
1255 
1256  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
1257 
1258  clock_t t1 = clock();
1259 
1260  morphee::stats::t_measMinMax(imIn, not_used, max);
1261  numVert = max;
1262  std::cout << "number of Vertices : " << numVert << std::endl;
1263 
1264  for (int i = 1; i <= numVert; i++) {
1265  boost::add_vertex(g);
1266  }
1267 
1268  vSource = boost::add_vertex(g);
1269  vSink = boost::add_vertex(g);
1270 
1271  clock_t tt_marker2 = 0, tt_marker3 = 0, tt_new_edge = 0, tt_old_edge = 0;
1272  clock_t t2 = clock();
1273  std::cout << "Nodes creation time : " << double(t2 - t1) / CLOCKS_PER_SEC
1274  << " seconds\n";
1275 
1276  std::cout << "Building Region Adjacency Graph Edges" << std::endl;
1277  t1 = clock();
1278 
1279  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1280  o1 = it.getOffset();
1281  // pas de raison que ce soit int ; il faut rendre ceci générique
1282  int val = imIn.pixelFromOffset(o1);
1283  int marker = imMarker.pixelFromOffset(o1);
1284  int val_prec = 0, marker_prec = 0;
1285 
1286  if (val > 0) {
1287  if (marker == 2 && marker_prec != marker && val_prec != val) {
1288  clock_t temps_marker2 = clock();
1289  boost::tie(e4, in1) = boost::edge(vSource, val, g);
1290 
1291  if (in1 == 0) {
1292  // std::cout<<"Add new edge marker 2"<<std::endl;
1293  boost::tie(e4, in1) = boost::add_edge(vSource, val, g);
1294  boost::tie(e3, in1) = boost::add_edge(val, vSource, g);
1295  capacity[e4] = (std::numeric_limits<double>::max)();
1296  capacity[e3] = (std::numeric_limits<double>::max)();
1297  rev[e4] = e3;
1298  rev[e3] = e4;
1299  }
1300  tt_marker2 += clock() - temps_marker2;
1301  } else if (marker == 3 && marker_prec != marker && val_prec != val) {
1302  clock_t temps_marker3 = clock();
1303  boost::tie(e3, in1) = boost::edge(vSink, val, g);
1304  if (in1 == 0) {
1305  // std::cout<<"Add new edge marker 3"<<std::endl;
1306  boost::tie(e4, in1) = boost::add_edge(val, vSink, g);
1307  boost::tie(e3, in1) = boost::add_edge(vSink, val, g);
1308  capacity[e4] = (std::numeric_limits<double>::max)();
1309  capacity[e3] = (std::numeric_limits<double>::max)();
1310  rev[e4] = e3;
1311  rev[e3] = e4;
1312  }
1313  tt_marker3 += clock() - temps_marker3;
1314  }
1315 
1316  neighb.setCenter(o1);
1317  // Enlever double et int ; prendre types génériques
1318  double val_grad_o1 = imGrad.pixelFromOffset(o1);
1319  // typename mageGrad::value_type val_grad_o1 =
1320  // imGrad.pixelFromOffset(o1);
1321  int val2_prec = val; // val de val2 precedente
1322 
1323  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
1324  const offset_t o2 = nit.getOffset();
1325 
1326  if (o2 > o1) {
1327  // enlever int ; prendre type générique
1328  int val2 = imIn.pixelFromOffset(o2);
1329  if (val != val2) {
1330  // enlever double ; prendre type générique
1331  double val_grad_o2 = imGrad.pixelFromOffset(o2);
1332  // prendre la distance L1, remplacer double et F_SIMPLE par
1333  // valeur générique
1334  double diff = t_Distance_L1(val_grad_o1, val_grad_o2);
1335  // double diff = std::abs(val_grad_o1 - val_grad_o2);
1336 
1337  F_SIMPLE dist =
1338  std::sqrt(std::pow(step_x * (it.getX() - nit.getX()), 2) +
1339  std::pow(step_y * (it.getY() - nit.getY()), 2) +
1340  std::pow(step_z * (it.getZ() - nit.getZ()), 2));
1341 
1342  // Cette fonction devrait être remplacée par une fonction
1343  // paramètre
1344  // double cost = 10000.0/(1.0+std::pow(diff/dist,4));
1345  double cost = dist / (1 + diff);
1346 
1347  if (val2_prec ==
1348  val2) // same val2 means same edge (thus, keep e3 and e4)
1349  {
1350  capacity[e4] = capacity[e4] + cost;
1351  capacity[e3] = capacity[e3] + cost;
1352  } else {
1353  boost::tie(e5, in1) = boost::edge(val, val2, g);
1354 
1355  if (in1 == 0) {
1356  clock_t temps_new_edge = clock();
1357  // std::cout<<"Add new edge "<< val<<" --
1358  // "<<val2<<std::endl;
1359  numEdges++;
1360  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
1361  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
1362  capacity[e4] = cost;
1363  capacity[e3] = cost;
1364  rev[e4] = e3;
1365  rev[e3] = e4;
1366  tt_new_edge += clock() - temps_new_edge;
1367 
1368  } else {
1369  clock_t temps_old_edge = clock();
1370  // std::cout<<"existing edge"<<std::endl;
1371  boost::tie(e4, in1) = boost::edge(val, val2, g);
1372  boost::tie(e3, in1) = boost::edge(val2, val, g);
1373  capacity[e4] = capacity[e4] + cost;
1374  capacity[e3] = capacity[e3] + cost;
1375  tt_old_edge += clock() - temps_old_edge;
1376  }
1377  val2_prec = val2;
1378  }
1379  }
1380  }
1381  }
1382  val_prec = val;
1383  marker_prec = marker;
1384  }
1385  }
1386 
1387  std::cout << "Number of edges : " << numEdges << std::endl;
1388  t2 = clock();
1389  std::cout << "Edges creation time : " << double(t2 - t1) / CLOCKS_PER_SEC
1390  << " seconds\n";
1391  std::cout << "Marker2 : " << double(tt_marker2) / CLOCKS_PER_SEC
1392  << " seconds\n";
1393  std::cout << "Marker3 : " << double(tt_marker3) / CLOCKS_PER_SEC
1394  << " seconds\n";
1395  std::cout << "New edges : " << double(tt_new_edge) / CLOCKS_PER_SEC
1396  << " seconds\n";
1397  std::cout << "Old edges : " << double(tt_old_edge) / CLOCKS_PER_SEC
1398  << " seconds\n";
1399 
1400  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
1401  boost::get(boost::vertex_index, g);
1402  std::vector<boost::default_color_type> color(boost::num_vertices(g));
1403 
1404  std::cout << "Compute Max flow" << std::endl;
1405  t1 = clock();
1406 #if BOOST_VERSION >= 104700
1407  double flow =
1408  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1409  &color[0], indexmap, vSource, vSink);
1410 #else
1411  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1412  &color[0], indexmap, vSource, vSink);
1413 #endif
1414  std::cout << "c The total flow:" << std::endl;
1415  std::cout << "s " << flow << std::endl;
1416  t2 = clock();
1417  std::cout << "Flow computation time : " << double(t2 - t1) / CLOCKS_PER_SEC
1418  << " seconds\n";
1419 
1420  t1 = clock();
1421  int miss = 0;
1422  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1423  o1 = it.getOffset();
1424  // enlever int ; prendre type générique
1425  int val = imIn.pixelFromOffset(o1);
1426 
1427  if (val == 0) {
1428  imOut.setPixel(o1, 0);
1429  } else {
1430  if (color[val] == color[vSource])
1431  imOut.setPixel(o1, 2);
1432  else if (color[val] == color[vSink])
1433  imOut.setPixel(o1, 3);
1434  else {
1435  imOut.setPixel(o1, 20);
1436  miss++;
1437  }
1438  }
1439  }
1440  t2 = clock();
1441  std::cout << "Computing imOut took : " << double(t2 - t1) / CLOCKS_PER_SEC
1442  << " seconds\n";
1443 
1444  if (miss > 0)
1445  std::cout << "WARNING : Missclassified nodes : " << miss << "\n";
1446 
1447  return RES_OK;
1448  }
1449 #endif
1450 
1451  /*
1452  *
1453  *
1454  *
1455  *
1456  */
1457  template <class T>
1458  RES_T geoCutsMultiWay_MinSurfaces(const Image<T> &imIn,
1459  const Image<T> &imGrad,
1460  const Image<T> &imMarker, const StrElt &nl,
1461  Image<T> &imOut)
1462  {
1463  SMIL_ENTER_FUNCTION("");
1464 
1465  ASSERT_ALLOCATED(&imIn, &imGrad, &imMarker, &imOut);
1466  ASSERT_SAME_SIZE(&imIn, &imGrad, &imMarker, &imOut);
1467 
1468  offset_t o0, o1;
1469 
1470  typename Image<T>::lineType bufIn = imIn.getPixels();
1471  typename Image<T>::lineType bufOut = imOut.getPixels();
1472  typename Image<T>::lineType bufMarker = imMarker.getPixels();
1473  typename Image<T>::lineType bufGrad = imGrad.getPixels();
1474 
1475  std::cout << "build Region Adjacency Graph" << std::endl;
1476 
1477  Graph_T g;
1478  EdgeCap_T capacity = boost::get(boost::edge_capacity, g);
1479  EdgeRevCap_T rev = boost::get(boost::edge_reverse, g);
1480  EdgeResCap_T residual_capacity =
1481  boost::get(boost::edge_residual_capacity, g);
1482 
1483  bool in1;
1484  Graph_T::edge_descriptor e1, e2, e3, e4;
1485  Graph_T::vertex_descriptor vSource, vSink;
1486  int numVert = 0;
1487  int numLabels = 0;
1488 
1489  int pixelCount = imIn.getPixelCount();
1490  for (o0 = 0; o0 < pixelCount; o0++) {
1491  int val = (int) bufIn[o0];
1492  int val2 = (int) bufMarker[o0];
1493 
1494  bufOut[o0] = 1;
1495  if (val2 > numLabels) {
1496  numLabels = val2;
1497  }
1498 
1499  if (val > numVert) {
1500  numVert = val;
1501  }
1502  }
1503 
1504  std::cout << "number of labels :" << numLabels << std::endl;
1505 
1506  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
1507 
1508  for (int i = 0; i <= numVert; i++) {
1509  boost::add_vertex(g);
1510  }
1511 
1512  vSource = boost::add_vertex(g);
1513  vSink = boost::add_vertex(g);
1514 
1515  std::vector<boost::default_color_type> color(boost::num_vertices(g));
1516  VertexIndex_T indexmap = boost::get(boost::vertex_index, g);
1517 
1518  std::cout << "build Region Adjacency Graph Edges" << std::endl;
1519 
1520  vector<IntPoint> pts = filterStrElt(nl);
1521  vector<IntPoint>::iterator itBegin, itEnd;
1522  itBegin = pts.begin();
1523  itEnd = pts.end();
1524 
1525  int width = imIn.getWidth();
1526  int height = imIn.getHeight();
1527  int depth = imIn.getDepth();
1528 
1529  int strideY = width;
1530  int strideZ = width * height;
1531  for (int z = 0; z < depth; z++) {
1532  off_t p0 = z * strideZ;
1533  for (int y = 0; y < height; y++) {
1534  p0 += y * strideY;
1535  for (int x = 0; x < width; x++) {
1536  int val = (int) bufIn[o1];
1537  // int marker = (int) bufMarker[o1]; // XXX unused ???
1538 
1539  o1 = p0 + x;
1540 
1541  vector<IntPoint>::iterator it;
1542  for (it = itBegin; it != itEnd; it++) {
1543  if (x + it->x > width - 1 || x + it->x < 0)
1544  continue;
1545  if (y + it->y > height - 1 || y + it->y < 0)
1546  continue;
1547  if (z + it->z > depth - 1 || z + it->z < 0)
1548  continue;
1549  offset_t o2 = o1 + it->z *strideZ + it->y *strideY + it->x;
1550  if (o2 <= o1)
1551  continue;
1552 
1553  int val2 = bufIn[o2];
1554  if (val == val2)
1555  continue;
1556 
1557  boost::tie(e3, in1) = boost::edge(val, val2, g);
1558  // std::cout<<in1<<std::endl;
1559  // std::cout<<"Compute Gradient"<<std::endl;
1560  double val3 = (double) bufGrad[o1];
1561  double val4 = (double) bufGrad[o2];
1562  double maxi = std::max(val3, val4);
1563  double cost = 10000.0 / (1.0 + std::pow(maxi, 4));
1564 
1565  if (in1 == 0) {
1566  // std::cout<<"Add new edge"<<std::endl;
1567  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
1568  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
1569  capacity[e4] = cost;
1570  capacity[e3] = cost;
1571  rev[e4] = e3;
1572  rev[e3] = e4;
1573  } else {
1574  // std::cout<<"existing edge"<<std::endl;
1575  boost::tie(e4, in1) = boost::edge(val, val2, g);
1576  boost::tie(e3, in1) = boost::edge(val2, val, g);
1577  capacity[e4] = capacity[e4] + cost;
1578  capacity[e3] = capacity[e3] + cost;
1579  }
1580  }
1581  }
1582  }
1583  }
1584 
1585  for (int nbk = 2; nbk <= numLabels; nbk++) {
1586  // for all pixels in imIn create a vertex
1587  for (int o0 = 0; o0 < pixelCount; o0++) {
1588  int val = (int) bufMarker[o0];
1589  int val2 = (int) bufIn[o0];
1590 
1591  if (val == nbk) {
1592  boost::tie(e4, in1) = boost::edge(vSource, val2, g);
1593  if (in1 == 0) {
1594  boost::tie(e4, in1) = boost::add_edge(vSource, val2, g);
1595  boost::tie(e3, in1) = boost::add_edge(val2, vSource, g);
1596  capacity[e4] = (std::numeric_limits<double>::max)();
1597  capacity[e3] = (std::numeric_limits<double>::max)();
1598  rev[e4] = e3;
1599  rev[e3] = e4;
1600  }
1601  } else if (val > 1 && val != nbk) {
1602  boost::tie(e4, in1) = boost::edge(val2, vSink, g);
1603  if (in1 == 0) {
1604  boost::tie(e4, in1) = boost::add_edge(val2, vSink, g);
1605  boost::tie(e3, in1) = boost::add_edge(vSink, val2, g);
1606  capacity[e4] = (std::numeric_limits<double>::max)();
1607  capacity[e3] = (std::numeric_limits<double>::max)();
1608  rev[e4] = e3;
1609  rev[e3] = e4;
1610  }
1611  }
1612  }
1613 
1614  std::cout << "Compute Max flow" << nbk << std::endl;
1615 #if BOOST_VERSION >= 104700
1616  double flow =
1617  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1618  &color[0], indexmap, vSource, vSink);
1619 #else
1620  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1621  &color[0], indexmap, vSource, vSink);
1622 #endif
1623 
1624  std::cout << "c The total flow:" << std::endl;
1625  std::cout << "s " << flow << std::endl << std::endl;
1626 
1627  // for all pixels in imIn create a vertex and an edge
1628  for (int o1 = 0; o1 < pixelCount; o1++) {
1629  int val = (int) bufIn[o1];
1630  int val2 = (int) bufOut[o1];
1631  int val3 = (int) bufMarker[o1];
1632 
1633  if (val2 == 1) {
1634  if (color[val] == color[vSource])
1635  bufOut[o1] = (T) nbk;
1636  }
1637 
1638  if (val3 == nbk) {
1639  boost::tie(e4, in1) = boost::edge(vSource, val, g);
1640  if (in1 == 1) {
1641  boost::remove_edge(vSource, val, g);
1642  boost::remove_edge(val, vSource, g);
1643  }
1644  } else if (val3 > 1 && val3 != nbk) {
1645  boost::tie(e4, in1) = boost::edge(val, vSink, g);
1646  if (in1 == 1) {
1647  boost::remove_edge(val, vSink, g);
1648  boost::remove_edge(vSink, val, g);
1649  }
1650  }
1651  }
1652  }
1653 
1654  return RES_OK;
1655  }
1656 
1657 #if 0
1658  /*
1659  *
1660  *
1661  *
1662  *
1663  */
1664  template <class ImageIn, class ImageGrad, class ImageMosaic,
1665  class ImageMarker, class SE, class ImageOut>
1666  RES_T geoCutsOptimize_Mosaic(
1667  const ImageIn &imIn, const ImageGrad &imGrad, const ImageMosaic &imMosaic,
1668  const ImageMarker &imMarker, const SE &nl, ImageOut &imOut)
1669  {
1670  SMIL_ENTER_FUNCTION("");
1671 
1672  std::cout << "Enter function t_geoCutsOptimize_Mosaic" << std::endl;
1673 
1674  if ((!imOut.isAllocated())) {
1675  SMIL_REGISTER_ERROR("Not allocated");
1676  return RES_NOT_ALLOCATED;
1677  }
1678 
1679  if ((!imIn.isAllocated())) {
1680  SMIL_REGISTER_ERROR("Not allocated");
1681  return RES_NOT_ALLOCATED;
1682  }
1683 
1684  if ((!imGrad.isAllocated())) {
1685  SMIL_REGISTER_ERROR("Not allocated");
1686  return RES_NOT_ALLOCATED;
1687  }
1688 
1689  if ((!imMosaic.isAllocated())) {
1690  SMIL_REGISTER_ERROR("Not allocated");
1691  return RES_NOT_ALLOCATED;
1692  }
1693 
1694  if ((!imMarker.isAllocated())) {
1695  SMIL_REGISTER_ERROR("Not allocated");
1696  return RES_NOT_ALLOCATED;
1697  }
1698 
1699  // common image iterator
1700  typename ImageIn::const_iterator it, iend;
1701  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
1702  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
1703  offset_t o0;
1704  offset_t o1;
1705 
1706  // needed for max flow: capacit map, rev_capacity map, etc.
1707  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
1708  boost::directedS>
1709  Traits;
1710  typedef boost::adjacency_list<
1711  boost::listS, boost::vecS, boost::directedS,
1712  boost::property<boost::vertex_name_t, std::string>,
1713  boost::property<
1714  boost::edge_capacity_t, double,
1715  boost::property<boost::edge_residual_capacity_t, double,
1716  boost::property<boost::edge_reverse_t,
1717  Traits::edge_descriptor>>>>
1718  Graph_d;
1719 
1720  Graph_d g;
1721 
1722  double sigma = 1.0;
1723 
1724  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
1725  boost::get(boost::edge_capacity, g);
1726 
1727  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
1728  get(boost::edge_reverse, g);
1729 
1730  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
1731  residual_capacity = get(boost::edge_residual_capacity, g);
1732 
1733  bool in1;
1734  Graph_d::edge_descriptor e1, e2, e3, e4;
1735  Graph_d::vertex_descriptor vSource, vSink;
1736  int numVert = 0;
1737  int numLabels = 0;
1738 
1739  std::cout << "build Region Adjacency Graph" << std::endl;
1740 
1741  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1742  // for all pixels in imIn create a vertex
1743  o0 = it.getOffset();
1744  int val = imIn.pixelFromOffset(o0);
1745  int val2 = imMarker.pixelFromOffset(o0);
1746 
1747  imOut.setPixel(o0, 1);
1748 
1749  if (val2 > numLabels) {
1750  numLabels = val2;
1751  }
1752 
1753  if (val > numVert) {
1754  numVert = val;
1755  }
1756  }
1757  std::cout << "numlabels = " << numLabels << std::endl;
1758 
1759  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
1760 
1761  for (int i = 0; i <= numVert; i++) {
1762  boost::add_vertex(g);
1763  }
1764 
1765  vSource = boost::add_vertex(g);
1766  vSink = boost::add_vertex(g);
1767 
1768  std::vector<boost::default_color_type> color(boost::num_vertices(g));
1769  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
1770  boost::get(boost::vertex_index, g);
1771 
1772  std::cout << "build Region Adjacency Graph Edges" << std::endl;
1773 
1774  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1775  // for all pixels in imIn create a vertex and an edge
1776  o1 = it.getOffset();
1777  int val = imIn.pixelFromOffset(o1);
1778  int marker = imMarker.pixelFromOffset(o1);
1779  neighb.setCenter(o1);
1780 
1781  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
1782  const offset_t o2 = nit.getOffset();
1783  int val2 = imIn.pixelFromOffset(o2);
1784 
1785  if (o2 > o1) {
1786  if (val != val2) {
1787  boost::tie(e3, in1) = boost::edge(val, val2, g);
1788  // std::cout<<in1<<std::endl;
1789  // std::cout<<"Compute Gradient"<<std::endl;
1790  double val3 = imGrad.pixelFromOffset(o1);
1791  double val4 = imGrad.pixelFromOffset(o2);
1792  double maxi = std::max(val3, val4);
1793  double cost = 1000 / (1 + 0.1 * (maxi) * (maxi));
1794 
1795  if (in1 == 0) {
1796  // std::cout<<"Add new edge"<<std::endl;
1797  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
1798  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
1799  capacity[e4] = cost;
1800  capacity[e3] = cost;
1801  rev[e4] = e3;
1802  rev[e3] = e4;
1803  } else {
1804  // std::cout<<"existing edge"<<std::endl;
1805  boost::tie(e4, in1) = boost::edge(val, val2, g);
1806  boost::tie(e3, in1) = boost::edge(val2, val, g);
1807  capacity[e4] = capacity[e4] + cost;
1808  capacity[e3] = capacity[e3] + cost;
1809  }
1810  }
1811  }
1812  }
1813  }
1814 
1815  for (int label1 = 1; label1 < numLabels; label1++) {
1816  for (int label2 = label1 + 1; label2 <= numLabels; label2++) {
1817  if (label1 != label2) {
1818  std::cout << "Optimize Pair of labels: " << label1 << " " << label2
1819  << std::endl;
1820 
1821  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1822  // for all pixels in imIn create a vertex
1823  o0 = it.getOffset();
1824  int val = imMarker.pixelFromOffset(o0);
1825  int val2 = imIn.pixelFromOffset(o0);
1826 
1827  if (val == label1) {
1828  boost::tie(e4, in1) = boost::edge(vSource, val2, g);
1829  if (in1 == 0) {
1830  boost::tie(e4, in1) = boost::add_edge(vSource, val2, g);
1831  boost::tie(e3, in1) = boost::add_edge(val2, vSource, g);
1832  capacity[e4] = (std::numeric_limits<double>::max)();
1833  capacity[e3] = (std::numeric_limits<double>::max)();
1834  rev[e4] = e3;
1835  rev[e3] = e4;
1836  }
1837  } else if (val == label2) {
1838  boost::tie(e4, in1) = boost::edge(val2, vSink, g);
1839  if (in1 == 0) {
1840  boost::tie(e4, in1) = boost::add_edge(val2, vSink, g);
1841  boost::tie(e3, in1) = boost::add_edge(vSink, val2, g);
1842  capacity[e4] = (std::numeric_limits<double>::max)();
1843  capacity[e3] = (std::numeric_limits<double>::max)();
1844  rev[e4] = e3;
1845  rev[e3] = e4;
1846  }
1847  }
1848  }
1849 
1850  std::cout << "Compute Max flow" << std::endl;
1851 #if BOOST_VERSION >= 104700
1852  double flow =
1853  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1854  &color[0], indexmap, vSource, vSink);
1855 #else
1856  double flow =
1857  kolmogorov_max_flow(g, capacity, residual_capacity, rev,
1858  &color[0], indexmap, vSource, vSink);
1859 #endif
1860  std::cout << "c The total flow:" << std::endl;
1861  std::cout << "s " << flow << std::endl << std::endl;
1862 
1863  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1864  // for all pixels in imIn create a vertex and an edge
1865  o1 = it.getOffset();
1866  int val = imIn.pixelFromOffset(o1);
1867  int val2 = imMosaic.pixelFromOffset(o1);
1868  int val3 = imMarker.pixelFromOffset(o1);
1869 
1870  if (val2 == label1 || val2 == label2) {
1871  if (color[val] == color[vSource])
1872  imOut.setPixel(o1, label1);
1873  else if (color[val] == color[vSink])
1874  imOut.setPixel(o1, label2);
1875  else
1876  imOut.setPixel(o1, label2);
1877  }
1878 
1879  if (val3 == label1) {
1880  boost::tie(e4, in1) = boost::edge(vSource, val, g);
1881  if (in1 == 1) {
1882  boost::remove_edge(vSource, val, g);
1883  boost::remove_edge(val, vSource, g);
1884  }
1885  } else if (val3 == label2) {
1886  boost::tie(e4, in1) = boost::edge(val, vSink, g);
1887  if (in1 == 1) {
1888  boost::remove_edge(val, vSink, g);
1889  boost::remove_edge(vSink, val, g);
1890  }
1891  }
1892  }
1893  }
1894  }
1895  }
1896 
1897  return RES_OK;
1898  }
1899 
1900  /*
1901  *
1902  *
1903  *
1904  *
1905  */
1906  template <class ImageIn, class ImageGrad, class ImageCurvature,
1907  class ImageMarker, typename _Beta, class SE, class ImageOut>
1908  RES_T geoCutsRegularized_MinSurfaces(
1909  const ImageIn &imIn, const ImageGrad &imGrad,
1910  const ImageCurvature &imCurvature, const ImageMarker &imMarker,
1911  const _Beta Beta, const SE &nl, ImageOut &imOut)
1912  {
1913  SMIL_ENTER_FUNCTION("");
1914 
1915  std::cout << "Enter function t_geoCutsRegularized_MinSurfaces"
1916  << std::endl;
1917 
1918  if ((!imOut.isAllocated())) {
1919  SMIL_REGISTER_ERROR("Not allocated");
1920  return RES_NOT_ALLOCATED;
1921  }
1922 
1923  if ((!imIn.isAllocated())) {
1924  SMIL_REGISTER_ERROR("Not allocated");
1925  return RES_NOT_ALLOCATED;
1926  }
1927 
1928  if ((!imGrad.isAllocated())) {
1929  SMIL_REGISTER_ERROR("Not allocated");
1930  return RES_NOT_ALLOCATED;
1931  }
1932 
1933  if ((!imCurvature.isAllocated())) {
1934  SMIL_REGISTER_ERROR("Not allocated");
1935  return RES_NOT_ALLOCATED;
1936  }
1937 
1938  if ((!imMarker.isAllocated())) {
1939  SMIL_REGISTER_ERROR("Not allocated");
1940  return RES_NOT_ALLOCATED;
1941  }
1942 
1943  // common image iterator
1944  typename ImageIn::const_iterator it, iend;
1945  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
1946  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
1947  offset_t o0;
1948  offset_t o1;
1949 
1950  // needed for max flow: capacit map, rev_capacity map, etc.
1951  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
1952  boost::directedS>
1953  Traits;
1954  typedef boost::adjacency_list<
1955  boost::listS, boost::vecS, boost::directedS,
1956  boost::property<boost::vertex_name_t, std::string>,
1957  boost::property<
1958  boost::edge_capacity_t, double,
1959  boost::property<boost::edge_residual_capacity_t, double,
1960  boost::property<boost::edge_reverse_t,
1961  Traits::edge_descriptor>>>>
1962  Graph_d;
1963 
1964  Graph_d g;
1965 
1966  double beta = Beta;
1967 
1968  std::cout << "Reg. :" << beta << std::endl;
1969 
1970  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
1971  boost::get(boost::edge_capacity, g);
1972 
1973  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
1974  get(boost::edge_reverse, g);
1975 
1976  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
1977  residual_capacity = get(boost::edge_residual_capacity, g);
1978 
1979  bool in1;
1980  Graph_d::edge_descriptor e1, e2, e3, e4;
1981  Graph_d::vertex_descriptor vSource, vSink;
1982  int numVert = 0;
1983 
1984  std::cout << "build Region Adjacency Graph" << std::endl;
1985 
1986  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
1987  // for all pixels in imIn create a vertex
1988  o0 = it.getOffset();
1989  int val = imIn.pixelFromOffset(o0);
1990  if (val > numVert) {
1991  numVert = val;
1992  }
1993  }
1994 
1995  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
1996 
1997  for (int i = 1; i <= numVert; i++) {
1998  boost::add_vertex(g);
1999  }
2000 
2001  vSource = boost::add_vertex(g);
2002  vSink = boost::add_vertex(g);
2003 
2004  std::cout << "build Region Adjacency Graph Edges" << std::endl;
2005 
2006  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2007  // for all pixels in imIn create a vertex and an edge
2008  o1 = it.getOffset();
2009  int val = imIn.pixelFromOffset(o1);
2010  int marker = imMarker.pixelFromOffset(o1);
2011 
2012  if (val > 0) {
2013  if (marker == 2) {
2014  boost::tie(e4, in1) = boost::edge(vSource, val, g);
2015  if (in1 == 0) {
2016  boost::tie(e4, in1) = boost::add_edge(vSource, val, g);
2017  boost::tie(e3, in1) = boost::add_edge(val, vSource, g);
2018  capacity[e4] = (std::numeric_limits<double>::max)();
2019  capacity[e3] = (std::numeric_limits<double>::max)();
2020  rev[e4] = e3;
2021  rev[e3] = e4;
2022  }
2023  } else if (marker == 3) {
2024  boost::tie(e3, in1) = boost::edge(vSink, val, g);
2025  if (in1 == 0) {
2026  boost::tie(e4, in1) = boost::add_edge(val, vSink, g);
2027  boost::tie(e3, in1) = boost::add_edge(vSink, val, g);
2028  capacity[e4] = (std::numeric_limits<double>::max)();
2029  capacity[e3] = (std::numeric_limits<double>::max)();
2030  rev[e4] = e3;
2031  rev[e3] = e4;
2032  }
2033  }
2034 
2035  neighb.setCenter(o1);
2036 
2037  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
2038  const offset_t o2 = nit.getOffset();
2039  int val2 = imIn.pixelFromOffset(o2);
2040 
2041  if (o2 > o1) {
2042  if (val != val2) {
2043  boost::tie(e3, in1) = boost::edge(val, val2, g);
2044 
2045  double val3 = imGrad.pixelFromOffset(o1);
2046  double val4 = imGrad.pixelFromOffset(o2);
2047 
2048  double val5 = imCurvature.pixelFromOffset(o1);
2049  double val6 = imCurvature.pixelFromOffset(o2);
2050 
2051  double maxigrad = std::max(val3, val4);
2052 
2053  double maxicurv = std::max(val5, val6);
2054 
2055  // if (maxicurv >0) std::cout<<maxicurv<<std::endl;
2056 
2057  double costcurvature = (beta) *maxicurv;
2058 
2059  double costgradient = 10000.0 / (1 + std::pow(maxigrad, 2));
2060 
2061  double cost = costgradient + costcurvature;
2062 
2063  if (in1 == 0) {
2064  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
2065  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
2066  capacity[e4] = cost;
2067  capacity[e3] = cost;
2068  rev[e4] = e3;
2069  rev[e3] = e4;
2070  } else {
2071  boost::tie(e4, in1) = boost::edge(val, val2, g);
2072  boost::tie(e3, in1) = boost::edge(val2, val, g);
2073  capacity[e4] = capacity[e4] + cost;
2074  capacity[e3] = capacity[e3] + cost;
2075  }
2076  }
2077  }
2078  }
2079  }
2080  }
2081 
2082  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
2083  boost::get(boost::vertex_index, g);
2084  std::vector<boost::default_color_type> color(boost::num_vertices(g));
2085 
2086  std::cout << "Compute Max flow" << std::endl;
2087 #if BOOST_VERSION >= 104700
2088  double flow =
2089  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2090  &color[0], indexmap, vSource, vSink);
2091 #else
2092  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2093  &color[0], indexmap, vSource, vSink);
2094 #endif
2095  std::cout << "c The total flow:" << std::endl;
2096  std::cout << "s " << flow << std::endl << std::endl;
2097 
2098  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2099  // for all pixels in imIn create a vertex and an edge
2100  o1 = it.getOffset();
2101  int val = imIn.pixelFromOffset(o1);
2102 
2103  if (val == 0) {
2104  imOut.setPixel(o1, 0);
2105  } else {
2106  if (color[val] == color[vSource])
2107  imOut.setPixel(o1, 3);
2108  if (color[val] == 1)
2109  imOut.setPixel(o1, 4);
2110  if (color[val] == color[vSink])
2111  imOut.setPixel(o1, 2);
2112  }
2113  }
2114 
2115  return RES_OK;
2116  }
2117 
2118  /*
2119  *
2120  *
2121  *
2122  *
2123  */
2124  template <class ImageIn, class ImageMosaic, class ImageMarker, class SE,
2125  class ImageOut>
2126  RES_T geoCutsSegment_Graph(const ImageIn &imIn, const ImageMosaic &imMosaic,
2127  const ImageMarker &imMarker, const SE &nl,
2128  ImageOut &imOut)
2129  {
2130  SMIL_ENTER_FUNCTION("");
2131 
2132  std::cout << "Enter function optimize mosaic t_geoCutsSegment_Graph"
2133  << std::endl;
2134 
2135  if ((!imOut.isAllocated())) {
2136  SMIL_REGISTER_ERROR("Not allocated");
2137  return RES_NOT_ALLOCATED;
2138  }
2139 
2140  if ((!imIn.isAllocated())) {
2141  SMIL_REGISTER_ERROR("Not allocated");
2142  return RES_NOT_ALLOCATED;
2143  }
2144 
2145  if ((!imMosaic.isAllocated())) {
2146  SMIL_REGISTER_ERROR("Not allocated");
2147  return RES_NOT_ALLOCATED;
2148  }
2149 
2150  if ((!imMarker.isAllocated())) {
2151  SMIL_REGISTER_ERROR("Not allocated");
2152  return RES_NOT_ALLOCATED;
2153  }
2154 
2155  // common image iterator
2156  typename ImageIn::const_iterator it, iend;
2157  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
2158  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
2159  offset_t o0;
2160  offset_t o1;
2161 
2162  // needed for max flow: capacit map, rev_capacity map, etc.
2163  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
2164  boost::directedS>
2165  Traits;
2166  typedef boost::adjacency_list<
2167  boost::listS, boost::vecS, boost::directedS,
2168  boost::property<boost::vertex_name_t, std::string>,
2169  boost::property<
2170  boost::edge_capacity_t, double,
2171  boost::property<boost::edge_residual_capacity_t, double,
2172  boost::property<boost::edge_reverse_t,
2173  Traits::edge_descriptor>>>>
2174  Graph_d;
2175 
2176  Graph_d g;
2177 
2178  double sigma = 1.0;
2179 
2180  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
2181  boost::get(boost::edge_capacity, g);
2182 
2183  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
2184  get(boost::edge_reverse, g);
2185 
2186  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
2187  residual_capacity = get(boost::edge_residual_capacity, g);
2188 
2189  bool in1;
2190  Graph_d::edge_descriptor e1, e2, e3, e4;
2191  Graph_d::vertex_descriptor vSource, vSink;
2192  int numVert = 0;
2193  double meanclasse1 = 0.0;
2194  double meanclasse12 = 0.0;
2195 
2196  double meanclasse2 = 0.5;
2197 
2198  double sigma1 = 0.25;
2199  double sigma2 = 0.5;
2200 
2201  double max_value = 0.0;
2202  double max_longueur = 0.0;
2203 
2204  std::cout << "build Region Adjacency Graph" << std::endl;
2205 
2206  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2207  // for all pixels in imIn create a vertex
2208  o0 = it.getOffset();
2209  int val = imMosaic.pixelFromOffset(o0);
2210  int val2 = imIn.pixelFromOffset(o0);
2211  int val3 = imMarker.pixelFromOffset(o0);
2212 
2213  if (val > numVert) {
2214  numVert = val;
2215  }
2216  if (val2 > max_value) {
2217  max_value = val2;
2218  }
2219  if (val3 > max_longueur) {
2220  max_longueur = val3;
2221  }
2222  }
2223 
2224  std::cout << "build Region Adjacency Graph Vertices" << std::endl;
2225 
2226  std::cout << "Number of Vertices : " << numVert << std::endl;
2227 
2228  std::cout << "Max value : " << max_value << std::endl;
2229 
2230  for (int i = 0; i <= numVert; i++) {
2231  boost::add_vertex(g);
2232  }
2233 
2234  vSource = boost::add_vertex(g);
2235  vSink = boost::add_vertex(g);
2236 
2237  std::cout << "build Region Adjacency Graph Edges" << std::endl;
2238 
2239  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2240  // for all pixels in imIn create a vertex and an edge
2241  o1 = it.getOffset();
2242  int val = imMosaic.pixelFromOffset(o1);
2243 
2244  if (val > 0) {
2245  boost::tie(e4, in1) = boost::edge(vSource, val, g);
2246 
2247  if (in1 == 0) {
2248  // std::cout<<"Add new edge marker 2"<<std::endl;
2249  double valee = (double) imIn.pixelFromOffset(o1) / max_value;
2250  double longueur =
2251  (double) imMarker.pixelFromOffset(o1) / max_longueur;
2252 
2253  double cost1 = 4 * (1 - longueur) + (valee - meanclasse1) *
2254  (valee - meanclasse1) /
2255  (2 * sigma1 * sigma1);
2256  double cost12 = 4 * (1 - longueur) + (valee - meanclasse12) *
2257  (valee - meanclasse12) /
2258  (2 * sigma1 * sigma1);
2259  double cost2 = 4 * (1 - 0.17) + (valee - meanclasse2) *
2260  (valee - meanclasse2) /
2261  (2 * sigma2 * sigma2);
2262 
2263  /*
2264  double cost1 =
2265  (valee-meanclasse1)*(valee-meanclasse1)/(2*sigma1*sigma1); double
2266  cost12 =
2267  (valee-meanclasse12)*(valee-meanclasse12)/(2*sigma1*sigma1);
2268  double cost2 =
2269  (valee-meanclasse2)*(valee-meanclasse2)/(2*sigma2*sigma2);
2270  */
2271 
2272  boost::tie(e4, in1) = boost::add_edge(vSource, val, g);
2273  boost::tie(e3, in1) = boost::add_edge(val, vSource, g);
2274  capacity[e4] = std::min(cost1, cost12);
2275  capacity[e3] = std::min(cost1, cost12);
2276  rev[e4] = e3;
2277  rev[e3] = e4;
2278 
2279  boost::tie(e4, in1) = boost::add_edge(vSink, val, g);
2280  boost::tie(e3, in1) = boost::add_edge(val, vSink, g);
2281  capacity[e4] = cost2;
2282  capacity[e3] = cost2;
2283  rev[e4] = e3;
2284  rev[e3] = e4;
2285  }
2286 
2287  neighb.setCenter(o1);
2288 
2289  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
2290  const offset_t o2 = nit.getOffset();
2291  int val2 = imMosaic.pixelFromOffset(o2);
2292 
2293  if (o2 > o1) {
2294  if (val != val2 && val2 > 0) {
2295  boost::tie(e3, in1) = boost::edge(val, val2, g);
2296 
2297  double valee1 = (double) imIn.pixelFromOffset(o1) / max_value;
2298  double valee2 = (double) imIn.pixelFromOffset(o2) / max_value;
2299 
2300  double longueur1 =
2301  (double) imMarker.pixelFromOffset(o1) / max_longueur;
2302  double longueur2 =
2303  (double) imMarker.pixelFromOffset(o2) / max_longueur;
2304 
2305  double cost_diff = 0.01 * std::exp(-0.01 * (valee1 - valee2) *
2306  (valee1 - valee2));
2307  double cost_longueur = 0.1;
2308 
2309  if (in1 == 0) {
2310  // std::cout<<"Add new edge"<<std::endl;
2311  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
2312  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
2313  capacity[e4] = cost_longueur;
2314  capacity[e3] = cost_longueur;
2315  rev[e4] = e3;
2316  rev[e3] = e4;
2317  }
2318  }
2319  }
2320  }
2321  }
2322  }
2323 
2324  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
2325  boost::get(boost::vertex_index, g);
2326  std::vector<boost::default_color_type> color(boost::num_vertices(g));
2327 
2328  std::cout << "Compute Max flow" << std::endl;
2329 #if BOOST_VERSION >= 104700
2330  double flow =
2331  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2332  &color[0], indexmap, vSource, vSink);
2333 #else
2334  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2335  &color[0], indexmap, vSource, vSink);
2336 #endif
2337  std::cout << "c The total flow:" << std::endl;
2338  std::cout << "s " << flow << std::endl << std::endl;
2339 
2340  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2341  // for all pixels in imIn create a vertex and an edge
2342  o1 = it.getOffset();
2343  int val = imMosaic.pixelFromOffset(o1);
2344 
2345  if (color[val] == color[vSource])
2346  imOut.setPixel(o1, 10);
2347  if (color[val] == 1)
2348  imOut.setPixel(o1, 0);
2349  if (color[val] == color[vSink])
2350  imOut.setPixel(o1, 30);
2351  }
2352 
2353  return RES_OK;
2354  }
2355 
2356  /*
2357  *
2358  *
2359  *
2360  *
2361  */
2362  template <class ImageIn, class ImageMosaic, class ImageMarker, typename _Beta,
2363  typename _Sigma, class SE, class ImageOut>
2364  RES_T MAP_MRF_edge_preserving(
2365  const ImageIn &imIn, const ImageMosaic &imMosaic,
2366  const ImageMarker &imMarker, const _Beta Beta, const _Sigma Sigma,
2367  const SE &nl, ImageOut &imOut)
2368  {
2369  SMIL_ENTER_FUNCTION("");
2370 
2371  std::cout << "Enter function t_MAP_MRF_edge_preserving" << std::endl;
2372 
2373  if ((!imOut.isAllocated())) {
2374  SMIL_REGISTER_ERROR("Not allocated");
2375  return RES_NOT_ALLOCATED;
2376  }
2377 
2378  if ((!imIn.isAllocated())) {
2379  SMIL_REGISTER_ERROR("Not allocated");
2380  return RES_NOT_ALLOCATED;
2381  }
2382 
2383  if ((!imMosaic.isAllocated())) {
2384  SMIL_REGISTER_ERROR("Not allocated");
2385  return RES_NOT_ALLOCATED;
2386  }
2387 
2388  if ((!imMarker.isAllocated())) {
2389  SMIL_REGISTER_ERROR("Not allocated");
2390  return RES_NOT_ALLOCATED;
2391  }
2392 
2393  // common image iterator
2394  typename ImageIn::const_iterator it, iend;
2395  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
2396  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
2397  offset_t o0;
2398  offset_t o1;
2399 
2400  // needed for max flow: capacit map, rev_capacity map, etc.
2401  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
2402  boost::directedS>
2403  Traits;
2404  typedef boost::adjacency_list<
2405  boost::listS, boost::vecS, boost::directedS,
2406  boost::property<boost::vertex_name_t, std::string>,
2407  boost::property<
2408  boost::edge_capacity_t, double,
2409  boost::property<boost::edge_residual_capacity_t, double,
2410  boost::property<boost::edge_reverse_t,
2411  Traits::edge_descriptor>>>>
2412  Graph_d;
2413 
2414  Graph_d g;
2415 
2416  double sigma = Sigma;
2417 
2418  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
2419  boost::get(boost::edge_capacity, g);
2420 
2421  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
2422  get(boost::edge_reverse, g);
2423 
2424  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
2425  residual_capacity = get(boost::edge_residual_capacity, g);
2426 
2427  bool in1;
2428  Graph_d::edge_descriptor e1, e2, e3, e4;
2429  Graph_d::vertex_descriptor vSource, vSink;
2430  int numVert = 1;
2431 
2432  std::cout << "build Region Adjacency Graph" << std::endl;
2433 
2434  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2435  // for all pixels in imIn create a vertex
2436  o0 = it.getOffset();
2437  int val = imMosaic.pixelFromOffset(o0);
2438  imOut.setPixel(o0, 0);
2439 
2440  if (val > numVert) {
2441  numVert = val;
2442  }
2443  }
2444 
2445  double *mean = new double[numVert + 1];
2446  int *nb_val = new int[numVert + 1];
2447  int *marker = new int[numVert + 1];
2448 
2449  for (int i = 0; i <= numVert; i++) {
2450  boost::add_vertex(g);
2451  mean[i] = 0;
2452  nb_val[i] = 0;
2453  marker[i] = 0;
2454  }
2455 
2456  vSource = boost::add_vertex(g);
2457  vSink = boost::add_vertex(g);
2458 
2459  double meanforeground = 0;
2460  double meanbackground = 0;
2461  double nb_foreground = 0;
2462  double nb_background = 0;
2463 
2464  std::cout << "Compute Mean Value in Regions" << std::endl;
2465 
2466  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2467  // for all pixels in imIn create a vertex
2468  o0 = it.getOffset();
2469  int val = imIn.pixelFromOffset(o0);
2470  int val2 = imMosaic.pixelFromOffset(o0);
2471  int val3 = imMarker.pixelFromOffset(o0);
2472  double valeur = (double) val;
2473 
2474  mean[val2] = mean[val2] + (valeur / 255.0);
2475  nb_val[val2] = nb_val[val2] + 1;
2476 
2477  if (val3 == 2) {
2478  meanforeground = meanforeground + (valeur / 255.0);
2479  nb_foreground++;
2480  marker[val2] = 2;
2481  } else if (val3 == 3) {
2482  meanbackground = meanbackground + (valeur / 255.0);
2483  nb_background++;
2484  marker[val2] = 3;
2485  }
2486  }
2487 
2488  meanforeground = meanforeground / nb_foreground;
2489  meanbackground = meanbackground / nb_background;
2490 
2491  std::cout << "Mean Foreground " << meanforeground << std::endl;
2492  std::cout << "Mean Background " << meanbackground << std::endl;
2493 
2494  std::cout << "Compute terminal links" << std::endl;
2495 
2496  double sigmab = 0.2;
2497 
2498  for (int i = 0; i <= numVert; i++) {
2499  mean[i] = mean[i] / (nb_val[i]);
2500 
2501  if (marker[i] == 2 && nb_val[i] > 0) {
2502  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
2503  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
2504  capacity[e4] = (std::numeric_limits<double>::max)();
2505  capacity[e3] = (std::numeric_limits<double>::max)();
2506  rev[e4] = e3;
2507  rev[e3] = e4;
2508  } else if (marker[i] == 3 && nb_val[i] > 0) {
2509  boost::tie(e4, in1) = boost::add_edge(vSink, i, g);
2510  boost::tie(e3, in1) = boost::add_edge(i, vSink, g);
2511  capacity[e4] = (std::numeric_limits<double>::max)();
2512  capacity[e3] = (std::numeric_limits<double>::max)();
2513  rev[e4] = e3;
2514  rev[e3] = e4;
2515 
2516  }
2517 
2518  else if (nb_val[i] > 0) {
2519  double valee = mean[i];
2520  double sigma = Sigma;
2521  double val2 = (valee - meanforeground) * (valee - meanforeground) /
2522  (2 * sigma * sigma);
2523  double val1 = (valee - meanbackground) * (valee - meanbackground) /
2524  (2 * sigmab * sigmab);
2525 
2526  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
2527  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
2528  capacity[e4] = nb_val[i] * val1;
2529  capacity[e3] = nb_val[i] * val1;
2530  rev[e4] = e3;
2531  rev[e3] = e4;
2532 
2533  boost::tie(e4, in1) = boost::add_edge(i, vSink, g);
2534  boost::tie(e3, in1) = boost::add_edge(vSink, i, g);
2535  capacity[e4] = nb_val[i] * val2;
2536  capacity[e3] = nb_val[i] * val2;
2537  rev[e4] = e3;
2538  rev[e3] = e4;
2539  }
2540  }
2541 
2542  std::cout << "build Region Adjacency Graph Edges" << std::endl;
2543 
2544  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2545  // for all pixels in imIn create a vertex and an edge
2546  o1 = it.getOffset();
2547  int val = imMosaic.pixelFromOffset(o1);
2548  neighb.setCenter(o1);
2549 
2550  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
2551  const offset_t o2 = nit.getOffset();
2552  int val2 = imMosaic.pixelFromOffset(o2);
2553 
2554  if (o2 > o1) {
2555  if (val != val2) {
2556  boost::tie(e3, in1) = boost::edge(val, val2, g);
2557  double cost = (double) Beta -
2558  ((double) Beta) *
2559  std::pow(std::abs(mean[val] - mean[val2]), 6.0);
2560  if (in1 == 0) {
2561  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
2562  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
2563  capacity[e4] = cost;
2564  capacity[e3] = cost;
2565  rev[e4] = e3;
2566  rev[e3] = e4;
2567  } else {
2568  boost::tie(e4, in1) = boost::edge(val, val2, g);
2569  boost::tie(e3, in1) = boost::edge(val2, val, g);
2570  capacity[e4] = capacity[e4] + cost;
2571  capacity[e3] = capacity[e3] + cost;
2572  }
2573  }
2574  }
2575  }
2576  }
2577 
2578  delete[] nb_val;
2579  delete[] mean;
2580  delete[] marker;
2581 
2582  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
2583  boost::get(boost::vertex_index, g);
2584  std::vector<boost::default_color_type> color(boost::num_vertices(g));
2585 
2586  std::cout << "Compute Max flow" << std::endl;
2587 
2588 #if BOOST_VERSION >= 104700
2589  double flow =
2590  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2591  &color[0], indexmap, vSource, vSink);
2592 #else
2593  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2594  &color[0], indexmap, vSource, vSink);
2595 #endif
2596  std::cout << "c The total flow:" << std::endl;
2597  std::cout << "s " << flow << std::endl << std::endl;
2598 
2599  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2600  // for all pixels in imIn create a vertex and an edge
2601  o1 = it.getOffset();
2602  int val = imMosaic.pixelFromOffset(o1);
2603 
2604  if (color[val] == color[vSource])
2605  imOut.setPixel(o1, 2);
2606  else if (color[val] == color[vSink])
2607  imOut.setPixel(o1, 3);
2608  else
2609  imOut.setPixel(o1, 4);
2610  }
2611 
2612  return RES_OK;
2613  }
2614 
2615  /*
2616  *
2617  *
2618  *
2619  *
2620  */
2621  template <class ImageIn, class ImageMosaic, class ImageMarker, typename _Beta,
2622  typename _Sigma, class SE, class ImageOut>
2623  RES_T MAP_MRF_Ising(const ImageIn &imIn, const ImageMosaic &imMosaic,
2624  const ImageMarker &imMarker, const _Beta Beta,
2625  const _Sigma Sigma, const SE &nl, ImageOut &imOut)
2626  {
2627  SMIL_ENTER_FUNCTION("");
2628 
2629  std::cout << "Enter function t_MAP_MRF_Ising" << std::endl;
2630 
2631  if ((!imOut.isAllocated())) {
2632  SMIL_REGISTER_ERROR("Not allocated");
2633  return RES_NOT_ALLOCATED;
2634  }
2635 
2636  if ((!imIn.isAllocated())) {
2637  SMIL_REGISTER_ERROR("Not allocated");
2638  return RES_NOT_ALLOCATED;
2639  }
2640 
2641  if ((!imMosaic.isAllocated())) {
2642  SMIL_REGISTER_ERROR("Not allocated");
2643  return RES_NOT_ALLOCATED;
2644  }
2645 
2646  if ((!imMarker.isAllocated())) {
2647  SMIL_REGISTER_ERROR("Not allocated");
2648  return RES_NOT_ALLOCATED;
2649  }
2650 
2651  // common image iterator
2652 
2653  typename ImageIn::const_iterator it, iend;
2654  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
2655  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
2656  offset_t o0;
2657  offset_t o1;
2658 
2659  // needed for max flow: capacit map, rev_capacity map, etc.
2660  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
2661  boost::directedS>
2662  Traits;
2663  typedef boost::adjacency_list<
2664  boost::listS, boost::vecS, boost::directedS,
2665  boost::property<boost::vertex_name_t, std::string>,
2666  boost::property<
2667  boost::edge_capacity_t, double,
2668  boost::property<boost::edge_residual_capacity_t, double,
2669  boost::property<boost::edge_reverse_t,
2670  Traits::edge_descriptor>>>>
2671  Graph_d;
2672 
2673  Graph_d g;
2674 
2675  double sigma = Sigma;
2676 
2677  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
2678  boost::get(boost::edge_capacity, g);
2679 
2680  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
2681  get(boost::edge_reverse, g);
2682 
2683  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
2684  residual_capacity = get(boost::edge_residual_capacity, g);
2685 
2686  bool in1;
2687  Graph_d::edge_descriptor e1, e2, e3, e4;
2688  Graph_d::vertex_descriptor vSource, vSink;
2689  int numVert = 1;
2690 
2691  std::cout << "build Region Adjacency Graph" << std::endl;
2692 
2693  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2694  // for all pixels in imIn create a vertex
2695  o0 = it.getOffset();
2696  int val = imMosaic.pixelFromOffset(o0);
2697  imOut.setPixel(o0, 0);
2698 
2699  if (val > numVert) {
2700  numVert = val;
2701  }
2702  }
2703 
2704  double *mean = new double[numVert + 1];
2705  int *nb_val = new int[numVert + 1];
2706  int *marker = new int[numVert + 1];
2707 
2708  for (int i = 0; i <= numVert; i++) {
2709  boost::add_vertex(g);
2710  mean[i] = 0;
2711  nb_val[i] = 0;
2712  marker[i] = 0;
2713  }
2714 
2715  vSource = boost::add_vertex(g);
2716  vSink = boost::add_vertex(g);
2717 
2718  double meanforeground = 0;
2719  double meanbackground = 0;
2720  double nb_foreground = 0;
2721  double nb_background = 0;
2722 
2723  std::cout << "Compute Mean Value in Regions" << std::endl;
2724 
2725  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2726  // for all pixels in imIn create a vertex
2727  o0 = it.getOffset();
2728  int val = imIn.pixelFromOffset(o0);
2729  int val2 = imMosaic.pixelFromOffset(o0);
2730  int val3 = imMarker.pixelFromOffset(o0);
2731  double valeur = (double) val;
2732 
2733  mean[val2] = mean[val2] + (valeur / 255.0);
2734  nb_val[val2] = nb_val[val2] + 1;
2735 
2736  if (val3 == 2) {
2737  meanforeground = meanforeground + (valeur / 255.0);
2738  nb_foreground++;
2739  marker[val2] = 2;
2740  } else if (val3 == 3) {
2741  meanbackground = meanbackground + (valeur / 255.0);
2742  nb_background++;
2743  marker[val2] = 3;
2744  }
2745  }
2746 
2747  meanforeground = meanforeground / nb_foreground;
2748  meanbackground = meanbackground / nb_background;
2749 
2750  std::cout << "Mean Foreground " << meanforeground << std::endl;
2751  std::cout << "Mean Background " << meanbackground << std::endl;
2752 
2753  std::cout << "Compute terminal links" << std::endl;
2754 
2755  double sigmab = 0.2;
2756  sigmab = Sigma;
2757 
2758  for (int i = 0; i <= numVert; i++) {
2759  mean[i] = mean[i] / (nb_val[i]);
2760 
2761  if (marker[i] == 2 && nb_val[i] > 0) {
2762  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
2763  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
2764  capacity[e4] = (std::numeric_limits<double>::max)();
2765  capacity[e3] = (std::numeric_limits<double>::max)();
2766  rev[e4] = e3;
2767  rev[e3] = e4;
2768  } else if (marker[i] == 3 && nb_val[i] > 0) {
2769  boost::tie(e4, in1) = boost::add_edge(vSink, i, g);
2770  boost::tie(e3, in1) = boost::add_edge(i, vSink, g);
2771  capacity[e4] = (std::numeric_limits<double>::max)();
2772  capacity[e3] = (std::numeric_limits<double>::max)();
2773  rev[e4] = e3;
2774  rev[e3] = e4;
2775 
2776  } else if (nb_val[i] > 0) {
2777  double valee = mean[i];
2778  double sigma = Sigma;
2779  double val2 = (valee - meanforeground) * (valee - meanforeground) /
2780  (2 * sigma * sigma);
2781  double val1 = (valee - meanbackground) * (valee - meanbackground) /
2782  (2 * sigmab * sigmab);
2783 
2784  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
2785  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
2786  capacity[e4] = nb_val[i] * val1;
2787  capacity[e3] = nb_val[i] * val1;
2788  rev[e4] = e3;
2789  rev[e3] = e4;
2790 
2791  boost::tie(e4, in1) = boost::add_edge(i, vSink, g);
2792  boost::tie(e3, in1) = boost::add_edge(vSink, i, g);
2793  capacity[e4] = nb_val[i] * val2;
2794  capacity[e3] = nb_val[i] * val2;
2795  rev[e4] = e3;
2796  rev[e3] = e4;
2797  }
2798  }
2799 
2800  delete[] nb_val;
2801  delete[] mean;
2802  int numEdges = 0;
2803 
2804  std::cout << "build Region Adjacency Graph Edges" << std::endl;
2805 
2806  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2807  // for all pixels in imIn create a vertex and an edge
2808  o1 = it.getOffset();
2809  int val = imMosaic.pixelFromOffset(o1);
2810  neighb.setCenter(o1);
2811 
2812  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
2813  const offset_t o2 = nit.getOffset();
2814  int val2 = imMosaic.pixelFromOffset(o2);
2815 
2816  if (o2 > o1) {
2817  if (val != val2) {
2818  boost::tie(e3, in1) = boost::edge(val, val2, g);
2819  double cost = (double) Beta;
2820  if (in1 == 0) {
2821  numEdges++;
2822  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
2823  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
2824  capacity[e4] = cost;
2825  capacity[e3] = cost;
2826  rev[e4] = e3;
2827  rev[e3] = e4;
2828  } else {
2829  boost::tie(e4, in1) = boost::edge(val, val2, g);
2830  boost::tie(e3, in1) = boost::edge(val2, val, g);
2831  capacity[e4] = capacity[e4] + cost;
2832  capacity[e3] = capacity[e3] + cost;
2833  }
2834  }
2835  }
2836  }
2837  }
2838 
2839  std::cout << "Number of vertices " << numVert << std::endl;
2840  std::cout << "Number of Edges " << numEdges << std::endl;
2841 
2842  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
2843  boost::get(boost::vertex_index, g);
2844  std::vector<boost::default_color_type> color(boost::num_vertices(g));
2845 
2846  std::cout << "Compute Max flow" << std::endl;
2847 #if BOOST_VERSION >= 104700
2848  double flow =
2849  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2850  &color[0], indexmap, vSource, vSink);
2851 #else
2852  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
2853  &color[0], indexmap, vSource, vSink);
2854 #endif
2855  std::cout << "c The total flow:" << std::endl;
2856  std::cout << "s " << flow << std::endl << std::endl;
2857 
2858  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2859  // for all pixels in imIn create a vertex and an edge
2860  o1 = it.getOffset();
2861  int val = imMosaic.pixelFromOffset(o1);
2862 
2863  if (color[val] == color[vSource])
2864  imOut.setPixel(o1, 2);
2865  else if (color[val] == color[vSink])
2866  imOut.setPixel(o1, 3);
2867  else
2868  imOut.setPixel(o1, 4);
2869  }
2870 
2871  return RES_OK;
2872  }
2873 
2874  /*
2875  *
2876  *
2877  *
2878  *
2879  */
2880  template <class ImageIn, class ImageMosaic, class ImageMarker, typename _Beta,
2881  typename _Sigma, class SE, class ImageOut>
2882  RES_T MAP_MRF_Potts(const ImageIn &imIn, const ImageMosaic &imMosaic,
2883  const ImageMarker &imMarker, const _Beta Beta,
2884  const _Sigma Sigma, const SE &nl, ImageOut &imOut)
2885  {
2886  SMIL_ENTER_FUNCTION("");
2887 
2888  std::cout << "Enter function t_MAP_MRF_Potts" << std::endl;
2889 
2890  if ((!imOut.isAllocated())) {
2891  SMIL_REGISTER_ERROR("Not allocated");
2892  return RES_NOT_ALLOCATED;
2893  }
2894 
2895  if ((!imIn.isAllocated())) {
2896  SMIL_REGISTER_ERROR("Not allocated");
2897  return RES_NOT_ALLOCATED;
2898  }
2899 
2900  if ((!imMosaic.isAllocated())) {
2901  SMIL_REGISTER_ERROR("Not allocated");
2902  return RES_NOT_ALLOCATED;
2903  }
2904 
2905  if ((!imMarker.isAllocated())) {
2906  SMIL_REGISTER_ERROR("Not allocated");
2907  return RES_NOT_ALLOCATED;
2908  }
2909 
2910  // common image iterator
2911  typename ImageIn::const_iterator it, iend;
2912  morphee::selement::Neighborhood<SE, ImageIn> neighb(imIn, nl);
2913  typename morphee::selement::Neighborhood<SE, ImageIn>::iterator nit, nend;
2914  offset_t o0;
2915  offset_t o1;
2916 
2917  // needed for max flow: capacit map, rev_capacity map, etc.
2918  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
2919  boost::directedS>
2920  Traits;
2921  typedef boost::adjacency_list<
2922  boost::listS, boost::vecS, boost::directedS,
2923  boost::property<boost::vertex_name_t, std::string>,
2924  boost::property<
2925  boost::edge_capacity_t, double,
2926  boost::property<boost::edge_residual_capacity_t, double,
2927  boost::property<boost::edge_reverse_t,
2928  Traits::edge_descriptor>>>>
2929  Graph_d;
2930 
2931  Graph_d g;
2932 
2933  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
2934  boost::get(boost::edge_capacity, g);
2935 
2936  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
2937  get(boost::edge_reverse, g);
2938 
2939  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
2940  residual_capacity = get(boost::edge_residual_capacity, g);
2941 
2942  bool in1;
2943  Graph_d::edge_descriptor e1, e2, e3, e4;
2944  Graph_d::vertex_descriptor vSource, vSink, vSource2, vSink2;
2945 
2946  int numVert = 0;
2947 
2948  double meanlabel1 = 0;
2949  double meanlabel2 = 0;
2950  double meanlabel3 = 0;
2951 
2952  double nb_label1 = 0;
2953  double nb_label2 = 0;
2954  double nb_label3 = 0;
2955 
2956  double Sigmal = (double) Sigma;
2957 
2958  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2959  // for all pixels in imIn create a vertex
2960  o0 = it.getOffset();
2961  int val = imMosaic.pixelFromOffset(o0);
2962 
2963  imOut.setPixel(o0, 0);
2964 
2965  if (val > numVert) {
2966  numVert = val;
2967  }
2968  }
2969 
2970  double *mean = new double[numVert + 1];
2971  int *nb_val = new int[numVert + 1];
2972  int *marker = new int[numVert + 1];
2973 
2974  std::cout << "number of regions : " << numVert << std::endl;
2975 
2976  for (int i = 0; i <= numVert; i++) {
2977  mean[i] = 0;
2978  nb_val[i] = 0;
2979  marker[i] = 0;
2980  }
2981 
2982  std::cout << "build Region Adjacency Graph" << std::endl;
2983 
2984  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
2985  // for all pixels in imIn create a vertex
2986  o0 = it.getOffset();
2987  imOut.setPixel(o0, 0);
2988 
2989  int val = imIn.pixelFromOffset(o0);
2990  int val2 = imMosaic.pixelFromOffset(o0);
2991  int val3 = imMarker.pixelFromOffset(o0);
2992 
2993  double valeur = (double) val;
2994 
2995  mean[val2] = mean[val2] + (valeur / 255.0);
2996  nb_val[val2] = nb_val[val2]++;
2997 
2998  if (val3 == 2) {
2999  meanlabel1 = meanlabel1 + (valeur / 255.0);
3000  nb_label1++;
3001  marker[val2] = 2;
3002  } else if (val3 == 3) {
3003  meanlabel2 = meanlabel2 + (valeur / 255.0);
3004  nb_label2++;
3005  marker[val2] = 3;
3006  } else if (val3 == 4) {
3007  meanlabel3 = meanlabel3 + (valeur / 255.0);
3008  nb_label3++;
3009  marker[val2] = 4;
3010  }
3011  }
3012 
3013  meanlabel1 = meanlabel1 / nb_label1;
3014  meanlabel2 = meanlabel2 / nb_label2;
3015  meanlabel3 = meanlabel3 / nb_label3;
3016 
3017  std::cout << "Mean Label 1 " << meanlabel1 << std::endl;
3018  std::cout << "Mean Label 2 " << meanlabel2 << std::endl;
3019  std::cout << "Mean Label 3 " << meanlabel3 << std::endl;
3020 
3021  std::cout << "Compute terminal links" << std::endl;
3022 
3023  for (int i = 0; i <= numVert; i++) {
3024  boost::add_vertex(g);
3025  mean[i] = mean[i] / (nb_val[i]);
3026  }
3027 
3028  vSource = boost::add_vertex(g);
3029  vSink = boost::add_vertex(g);
3030 
3031  for (int i = 0; i <= numVert; i++) {
3032  if (marker[i] == 4 && nb_val[i] > 0) {
3033  boost::tie(e4, in1) = boost::add_edge(vSink, i, g);
3034  boost::tie(e3, in1) = boost::add_edge(i, vSink, g);
3035  capacity[e4] = (std::numeric_limits<double>::max)();
3036  capacity[e3] = (std::numeric_limits<double>::max)();
3037  rev[e4] = e3;
3038  rev[e3] = e4;
3039  } else if (marker[i] > 0 && nb_val[i] > 0) {
3040  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
3041  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
3042  capacity[e4] = (std::numeric_limits<double>::max)();
3043  capacity[e3] = (std::numeric_limits<double>::max)();
3044  rev[e4] = e3;
3045  rev[e3] = e4;
3046  } else if (nb_val[i] > 0) {
3047  double valee = mean[i];
3048 
3049  double val1 =
3050  (valee - meanlabel3) * (valee - meanlabel3) / (2 * 0.05 * 0.05);
3051 
3052  double val2 = 0;
3053  double val3 = 0;
3054 
3055  val2 =
3056  (valee - meanlabel1) * (valee - meanlabel1) / (2 * Sigmal * Sigmal);
3057  val3 = (valee - meanlabel2) * (valee - meanlabel2) / (2 * 0.2 * 0.2);
3058 
3059  boost::tie(e4, in1) = boost::add_edge(vSource, i, g);
3060  boost::tie(e3, in1) = boost::add_edge(i, vSource, g);
3061  capacity[e4] = val1;
3062  capacity[e3] = val1;
3063  rev[e4] = e3;
3064  rev[e3] = e4;
3065 
3066  boost::tie(e4, in1) = boost::add_edge(i, vSink, g);
3067  boost::tie(e3, in1) = boost::add_edge(vSink, i, g);
3068  capacity[e4] = std::min(val2, val3);
3069  capacity[e3] = std::min(val2, val3);
3070  rev[e4] = e3;
3071  rev[e3] = e4;
3072  }
3073  }
3074 
3075  std::cout << "build Region Adjacency Graph Edges" << std::endl;
3076 
3077  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
3078  // for all pixels in imIn create a vertex and an edge
3079  o1 = it.getOffset();
3080  int val = imMosaic.pixelFromOffset(o1);
3081  neighb.setCenter(o1);
3082 
3083  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
3084  const offset_t o2 = nit.getOffset();
3085  int val2 = imMosaic.pixelFromOffset(o2);
3086 
3087  if (o2 > o1) {
3088  if (val != val2) {
3089  boost::tie(e3, in1) = boost::edge(val, val2, g);
3090  double cost = (double) Beta;
3091  if (in1 == 0) {
3092  boost::tie(e4, in1) = boost::add_edge(val, val2, g);
3093  boost::tie(e3, in1) = boost::add_edge(val2, val, g);
3094  capacity[e4] = cost;
3095  capacity[e3] = cost;
3096  rev[e4] = e3;
3097  rev[e3] = e4;
3098  } else {
3099  boost::tie(e4, in1) = boost::edge(val, val2, g);
3100  boost::tie(e3, in1) = boost::edge(val2, val, g);
3101  capacity[e4] = capacity[e4] + cost;
3102  capacity[e3] = capacity[e3] + cost;
3103  }
3104  }
3105  }
3106  }
3107  }
3108 
3109  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
3110  boost::get(boost::vertex_index, g);
3111  std::vector<boost::default_color_type> color(boost::num_vertices(g));
3112 
3113  std::cout << "Compute Max flow" << std::endl;
3114 #if BOOST_VERSION >= 104700
3115  double flow =
3116  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
3117  &color[0], indexmap, vSource, vSink);
3118 #else
3119  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
3120  &color[0], indexmap, vSource, vSink);
3121 #endif
3122  std::cout << "c The total flow:" << std::endl;
3123  std::cout << "s " << flow << std::endl << std::endl;
3124 
3125  std::cout << "Source Label:" << color[vSource] << std::endl;
3126  std::cout << "Sink Label:" << color[vSink] << std::endl;
3127 
3128  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
3129  // for all pixels in imIn create a vertex and an edge
3130  o1 = it.getOffset();
3131  int valimout = imOut.pixelFromOffset(o1);
3132  int val = imMosaic.pixelFromOffset(o1);
3133 
3134  if (nb_val[val] > 0) {
3135  if (valimout == 0) {
3136  if (color[val] == 4)
3137  imOut.setPixel(o1, 4);
3138  }
3139  }
3140  }
3141 
3142  Graph_d g2;
3143 
3144  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity2 =
3145  boost::get(boost::edge_capacity, g2);
3146 
3147  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev2 =
3148  get(boost::edge_reverse, g2);
3149 
3150  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
3151  residual_capacity2 = get(boost::edge_residual_capacity, g2);
3152 
3153  for (int i = 0; i <= numVert; i++) {
3154  boost::add_vertex(g2);
3155  }
3156 
3157  vSource2 = boost::add_vertex(g2);
3158  vSink2 = boost::add_vertex(g2);
3159 
3160  for (int i = 0; i <= numVert; i++) {
3161  if (marker[i] == 3 && nb_val[i] > 0) {
3162  boost::tie(e4, in1) = boost::add_edge(vSink2, i, g2);
3163  boost::tie(e3, in1) = boost::add_edge(i, vSink2, g2);
3164  capacity2[e4] = (std::numeric_limits<double>::max)();
3165  capacity2[e3] = (std::numeric_limits<double>::max)();
3166  rev2[e4] = e3;
3167  rev2[e3] = e4;
3168  } else if (marker[i] == 2 && nb_val[i] > 0) {
3169  boost::tie(e4, in1) = boost::add_edge(vSource2, i, g2);
3170  boost::tie(e3, in1) = boost::add_edge(i, vSource2, g2);
3171  capacity2[e4] = (std::numeric_limits<double>::max)();
3172  capacity2[e3] = (std::numeric_limits<double>::max)();
3173  rev2[e4] = e3;
3174  rev2[e3] = e4;
3175  } else if (nb_val[i] > 0) {
3176  double valee = mean[i];
3177  double val1 =
3178  (valee - meanlabel2) * (valee - meanlabel2) / (2 * 0.2 * 0.2);
3179  double val2 =
3180  (valee - meanlabel1) * (valee - meanlabel1) / (2 * Sigmal * Sigmal);
3181 
3182  boost::tie(e4, in1) = boost::add_edge(vSource2, i, g2);
3183  boost::tie(e3, in1) = boost::add_edge(i, vSource2, g2);
3184  capacity2[e4] = val1;
3185  capacity2[e3] = val1;
3186  rev2[e4] = e3;
3187  rev2[e3] = e4;
3188 
3189  boost::tie(e4, in1) = boost::add_edge(i, vSink2, g2);
3190  boost::tie(e3, in1) = boost::add_edge(vSink2, i, g2);
3191  capacity2[e4] = val2;
3192  capacity2[e3] = val2;
3193  rev2[e4] = e3;
3194  rev2[e3] = e4;
3195  }
3196  }
3197 
3198  std::cout << "build Region Adjacency Graph Edges" << std::endl;
3199 
3200  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
3201  // for all pixels in imIn create a vertex and an edge
3202  o1 = it.getOffset();
3203  int val = imMosaic.pixelFromOffset(o1);
3204  neighb.setCenter(o1);
3205 
3206  for (nit = neighb.begin(), nend = neighb.end(); nit != nend; ++nit) {
3207  const offset_t o2 = nit.getOffset();
3208  int val2 = imMosaic.pixelFromOffset(o2);
3209 
3210  if (o2 > o1) {
3211  if (val != val2) {
3212  boost::tie(e3, in1) = boost::edge(val, val2, g2);
3213  double cost = (double) Beta;
3214  if (in1 == 0) {
3215  boost::tie(e4, in1) = boost::add_edge(val, val2, g2);
3216  boost::tie(e3, in1) = boost::add_edge(val2, val, g2);
3217  capacity2[e4] = cost;
3218  capacity2[e3] = cost;
3219  rev2[e4] = e3;
3220  rev2[e3] = e4;
3221  } else {
3222  boost::tie(e4, in1) = boost::edge(val, val2, g2);
3223  boost::tie(e3, in1) = boost::edge(val2, val, g2);
3224  capacity2[e4] = capacity2[e4] + cost;
3225  capacity2[e3] = capacity2[e3] + cost;
3226  }
3227  }
3228  }
3229  }
3230  }
3231 
3232  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap2 =
3233  boost::get(boost::vertex_index, g2);
3234  std::vector<boost::default_color_type> color2(boost::num_vertices(g2));
3235 
3236  std::cout << "Compute Max flow" << std::endl;
3237 #if BOOST_VERSION >= 104700
3238  flow = boykov_kolmogorov_max_flow(g2, capacity2, residual_capacity2, rev2,
3239  &color2[0], indexmap2, vSource2, vSink2);
3240 #else
3241  flow = kolmogorov_max_flow(g2, capacity2, residual_capacity2, rev2,
3242  &color2[0], indexmap2, vSource2, vSink2);
3243 #endif
3244  std::cout << "c The total flow:" << std::endl;
3245  std::cout << "s " << flow << std::endl << std::endl;
3246 
3247  std::cout << "Source Label:" << color[vSource2] << std::endl;
3248  std::cout << "Sink Label:" << color[vSink2] << std::endl;
3249 
3250  for (it = imIn.begin(), iend = imIn.end(); it != iend; ++it) {
3251  // for all pixels in imIn create a vertex and an edge
3252  o1 = it.getOffset();
3253  int valimout = imOut.pixelFromOffset(o1);
3254  int val = imMosaic.pixelFromOffset(o1);
3255 
3256  if (nb_val[val] > 0) {
3257  if (valimout == 0) {
3258  if (color2[val] == 4)
3259  imOut.setPixel(o1, 3);
3260  }
3261  }
3262  }
3263 
3264  delete[] nb_val;
3265  delete[] mean;
3266 
3267  return RES_OK;
3268  }
3269 #endif
3270 
3271 /*
3272  *
3273  *
3274  *
3275  *
3276  */
3277 #define indent " "
3278 
3279  inline void printSE(StrElt ss)
3280  {
3281  cout << indent << "Structuring Element" << endl;
3282  cout << indent << "Type: " << ss.seT << endl;
3283  cout << indent << "Size: " << ss.size << endl;
3284  size_t ptNbr = ss.points.size();
3285  cout << indent << "Point Nbr: " << ptNbr << endl;
3286  if (ptNbr == 0)
3287  return;
3288 
3289  vector<IntPoint>::iterator itStart = ss.points.begin();
3290  vector<IntPoint>::iterator it, itEnd;
3291  itEnd = ss.points.end();
3292 
3293  vector<IntPoint> z(0);
3294  z = filterStrElt(ss);
3295  itStart = z.begin();
3296  itEnd = z.end();
3297  for (it = itStart; it != itEnd; it++) {
3298  // cout << "z " << "* " << it->x << " " << it->y << " " << it->z <<
3299  // endl;
3300  }
3301  cout << " =================== " << endl;
3302  }
3303 
3304  void testHandleSE(StrElt &se)
3305  {
3306  se = CubeSE(2);
3307  printSE(se);
3308  se = se.noCenter();
3309  printSE(se);
3310  int sz = se.getSize();
3311  se = se.homothety(sz);
3312  printSE(se);
3313 #if 0
3314  imIn.printSelf();
3315  cout << "* width " << imIn.getWidth() << endl;
3316  cout << "* height " << imIn.getHeight() << endl;
3317  cout << "* depth " << imIn.getDepth() << endl;
3318 #endif
3319  }
3320 
3321 } // namespace smil
3322 
3323 #endif // D_MOSAIC_GEOCUTS_IMPL_HPP
smil::sqrt
RES_T sqrt(const Image< T1 > &imIn, Image< T2 > &imOut)
sqrt() - square root of an image
Definition: DImageArith.hpp:926
smil::diff
RES_T diff(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
diff() - Difference between two images.
Definition: DImageArith.hpp:448
smil::exp
RES_T exp(const Image< T1 > &imIn, Image< T2 > &imOut, int base=0)
exp() - exponential of an image
Definition: DImageArith.hpp:881
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::mean
RES_T mean(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Mean filter.
Definition: DMorphoFilter.hpp:258
smil::label
size_t label(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
label() - Image labelization
Definition: DMorphoLabel.hpp:564
smil::maxVal
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
maxVal() - Max value of an image
Definition: DMeasures.hpp:348
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