SMIL  1.0.4
Watershed.hpp
1 /* __HEAD__
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * Copyright (c) 2017-2023, Centre de Morphologie Mathematique
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * * Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  * * Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in the
13  * documentation and/or other materials provided with the distribution.
14  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
15  * names of its contributors may be used to endorse or promote products
16  * derived from this software without specific prior written permission.
17  *
18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE
22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
28  * THE POSSIBILITY OF SUCH DAMAGE.
29  *
30  * Description :
31  * Porting GeoCuts module from Morph-M - This is the Watershed submodule
32  *
33  * History :
34  * - 20/03/2019 - by Jose-Marcio Martins da Cruz
35  * Just created it...
36  * - XX/XX/XXXX -
37  *
38  *
39  * __HEAD__ - Stop here !
40  */
41 
42 #ifndef _D_GEOCUTS_WATERSHED_HPP_
43 #define _D_GEOCUTS_WATERSHED_HPP_
44 
45 /*
46  *
47  *
48  *
49  */
50 #ifndef __BOOST_INCLUDED__
51 #define __BOOST_INCLUDED__
52 
53 #include <boost/config.hpp>
54 // for boost::tie
55 #include <boost/utility.hpp>
56 // for boost::graph_traits
57 #include <boost/graph/graph_traits.hpp>
58 #include <boost/graph/adjacency_list.hpp>
59 // JOE #include <boost/graph/graphviz.hpp>
60 // JOE #include <boost/graph/prim_minimum_spanning_tree.hpp>
61 // JOE #include <boost/graph/dijkstra_shortest_paths.hpp>
62 // JOE #include <boost/graph/johnson_all_pairs_shortest.hpp>
63 
64 #include <boost/version.hpp>
65 
66 #if 0
67 // FROM STAWIASKI JAN 2012
68 #include "../boost_ext/kolmogorov_max_flow_min_cost.hpp"
69 //#include "../boost_ext/maximum_spanning_tree.hpp"
70 //STAWIASKI JAN2012 commented, why?
71 //#include "../boost_ext/boost_compare.hpp"
72 #include <boost/graph/connected_components.hpp>
73 #endif
74 
75 #endif // __BOOST_INCLUDED__
76 
77 #include <vector>
78 
79 using namespace geocuts;
80 
81 namespace smil
82 {
83  /*
84  *
85  *
86  *
87  */
88  template <class T1, class T2>
89  RES_T geoCutsWatershed_MinCut(const Image<T1> &imIn,
90  const Image<T2> &imMarker, const double Power,
91  const StrElt &nl, Image<T2> &imOut)
92  {
93  std::cout << "Enter function Geo-Cuts Watershed" << std::endl;
94 
95  ASSERT_ALLOCATED(&imIn, &imMarker, &imOut);
96  ASSERT_SAME_SIZE(&imIn, &imMarker, &imOut);
97 
98  typename Image<T1>::lineType bufIn = imIn.getPixels();
99  typename Image<T2>::lineType bufMarker = imMarker.getPixels();
100  typename Image<T2>::lineType bufOut = imOut.getPixels();
101 
102  double exposant = Power;
103 
104  // needed for max flow: capacit map, rev_capacity map, etc.
105  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
106  boost::directedS>
107  Traits_d;
108 
109  typedef boost::adjacency_list<
110  boost::vecS, boost::vecS, boost::directedS,
111  boost::property<boost::vertex_name_t, std::string>,
112  boost::property<
113  boost::edge_capacity_t, double,
114  boost::property<boost::edge_residual_capacity_t, double,
115  boost::property<boost::edge_reverse_t,
116  Traits_d::edge_descriptor>>>>
117  Graph_d;
118 
119  Graph_d g;
120 
121  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
122  boost::get(boost::edge_capacity, g);
123 
124  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
125  get(boost::edge_reverse, g);
126 
127  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
128  residual_capacity = get(boost::edge_residual_capacity, g);
129 
130  Graph_d::edge_descriptor e1, e2, e3, e4;
131  Graph_d::vertex_descriptor vSource, vSink;
132  T1 numVertex = maxVal(imIn);
133  // T2 numLabels = maxVal(imMarker);
134 
135  std::cout << "build graph vertices" << std::endl;
136  std::cout << "number of vertices : " << numVertex << std::endl;
137 
138  size_t pixelCount = imIn.getPixelCount();
139  for (off_t i = 0; i < (off_t) pixelCount; i++) {
140  boost::add_vertex(g);
141  bufOut[i] = (T2) 1;
142  }
143 
144  vSource = boost::add_vertex(g);
145  vSink = boost::add_vertex(g);
146 
147  int width = imIn.getWidth();
148  int height = imIn.getHeight();
149  int depth = imIn.getDepth();
150 
151  off_t strideY = width;
152  off_t strideZ = width * height;
153  for (int z = 0; z < depth; z++) {
154  off_t p0Z = z * strideZ;
155  for (int y = 0; y < height; y++) {
156  off_t p0Y = y * strideY;
157  for (int x = 0; x < width; x++) {
158  off_t o1 = p0Z + p0Y + x;
159 
160  // iterators on the Structuring Element
161  vector<IntPoint> pts = filterStrElt(nl);
162  vector<IntPoint>::iterator itBegin, itEnd, it;
163  itBegin = pts.begin();
164  itEnd = pts.end();
165 
166  T1 val1 = bufIn[o1];
167  T2 marker = bufMarker[o1];
168 
169  bool hasEdge = false;
170  if (marker == 2) {
171  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o1, g);
172  boost::tie(e3, hasEdge) = boost::add_edge(o1, vSource, g);
173  capacity[e4] = (std::numeric_limits<double>::max)();
174  capacity[e3] = (std::numeric_limits<double>::max)();
175  rev[e4] = e3;
176  rev[e3] = e4;
177  } else if (marker == 3) {
178  boost::tie(e4, hasEdge) = boost::add_edge(o1, vSink, g);
179  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o1, g);
180  capacity[e4] = (std::numeric_limits<double>::max)();
181  capacity[e3] = (std::numeric_limits<double>::max)();
182  rev[e4] = e3;
183  rev[e3] = e4;
184  }
185 
186  for (it = itBegin; it != itEnd; it++) {
187  if (x + it->x > width - 1 || x + it->x < 0)
188  continue;
189  if (y + it->y > height - 1 || y + it->y < 0)
190  continue;
191  if (z + it->z > depth - 1 || z + it->z < 0)
192  continue;
193 
194  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
195  if (o2 <= o1)
196  continue;
197 
198  T2 val2 = bufIn[o2];
199  double valeur = (255.0 / (std::abs((double) (val1 - val2)) + 1));
200  double cost = std::pow(valeur, exposant);
201 
202  bool hasEdge = false;
203  boost::tie(e4, hasEdge) = boost::add_edge(o1, o2, g);
204  boost::tie(e3, hasEdge) = boost::add_edge(o2, o1, g);
205  capacity[e4] = cost;
206  capacity[e3] = cost;
207  rev[e4] = e3;
208  rev[e3] = e4;
209  }
210  }
211  }
212  }
213 
214  std::cout << "Compute Max flow" << std::endl;
215  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
216  boost::get(boost::vertex_index, g);
217  std::vector<boost::default_color_type> color(boost::num_vertices(g));
218 #if BOOST_VERSION >= 104700
219  double flow =
220  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
221  &color[0], indexmap, vSource, vSink);
222 #else
223  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
224  &color[0], indexmap, vSource, vSink);
225 #endif
226  std::cout << "c The total flow:" << std::endl;
227  std::cout << "s " << flow << std::endl << std::endl;
228 
229  // for all pixels in imIn create a vertex and an edge
230  for (off_t o1 = 0; o1 < (off_t) pixelCount; o1++) {
231  if (color[o1] == color[vSource])
232  bufOut[o1] = (T2) 2;
233  else if (color[o1] == color[vSink])
234  bufOut[o1] = (T2) 3;
235  else if (color[o1] == 1)
236  bufOut[o1] = (T2) 4;
237  }
238 
239  return RES_OK;
240  }
241 
242  /*
243  *
244  *
245  *
246  */
247  template <class T1, class T2>
248  RES_T geoCutsMultiway_Watershed(const Image<T1> &imIn,
249  const Image<T2> &imMarker, const double Power,
250  const StrElt &nl, Image<T2> &imOut)
251  {
252  std::cout << "Enter function Multi way watershed" << std::endl;
253  ASSERT_ALLOCATED(&imIn, &imMarker, &imOut);
254  ASSERT_SAME_SIZE(&imIn, &imMarker, &imOut);
255 
256  typename Image<T1>::lineType bufIn = imIn.getPixels();
257  typename Image<T2>::lineType bufMarker = imMarker.getPixels();
258  typename Image<T2>::lineType bufOut = imOut.getPixels();
259 
260  double exposant = Power;
261 
262  // needed for max flow: capacit map, rev_capacity map, etc.
263  typedef boost::adjacency_list_traits<boost::vecS, boost::vecS,
264  boost::directedS>
265  Traits_d;
266 
267  typedef boost::adjacency_list<
268  boost::vecS, boost::vecS, boost::directedS,
269  boost::property<boost::vertex_name_t, std::string>,
270  boost::property<
271  boost::edge_capacity_t, double,
272  boost::property<boost::edge_residual_capacity_t, double,
273  boost::property<boost::edge_reverse_t,
274  Traits_d::edge_descriptor>>>>
275  Graph_d;
276 
277  Graph_d g;
278 
279  boost::property_map<Graph_d, boost::edge_capacity_t>::type capacity =
280  boost::get(boost::edge_capacity, g);
281 
282  boost::property_map<Graph_d, boost::edge_reverse_t>::type rev =
283  get(boost::edge_reverse, g);
284 
285  boost::property_map<Graph_d, boost::edge_residual_capacity_t>::type
286  residual_capacity = get(boost::edge_residual_capacity, g);
287 
288  Graph_d::edge_descriptor e1, e2, e3, e4;
289  Graph_d::vertex_descriptor vSource, vSink;
290  T1 numVertex = maxVal(imIn);
291  T2 numLabels = maxVal(imMarker);
292 
293  std::cout << "build graph vertices" << std::endl;
294 
295  size_t pixelCount = imIn.getPixelCount();
296  for (off_t i = 0; i < (off_t) pixelCount; i++) {
297  boost::add_vertex(g);
298  bufOut[i] = (T2) 1;
299  }
300 
301  std::cout << "number of Labels: " << numLabels << std::endl;
302  std::cout << "number of vertices: " << numVertex << std::endl;
303 
304  vSource = boost::add_vertex(g);
305  vSink = boost::add_vertex(g);
306 
307  std::cout << "build graph edges" << std::endl;
308 
309  int width = imIn.getWidth();
310  int height = imIn.getHeight();
311  int depth = imIn.getDepth();
312 
313  off_t strideY = width;
314  off_t strideZ = width * height;
315  for (int z = 0; z < depth; z++) {
316  off_t p0Z = z * strideZ;
317  for (int y = 0; y < height; y++) {
318  off_t p0Y = y * strideY;
319  for (int x = 0; x < width; x++) {
320  off_t o1 = p0Z + p0Y + x;
321 
322  // iterators on the Structuring Element
323  vector<IntPoint> pts = filterStrElt(nl);
324  vector<IntPoint>::iterator itBegin, itEnd, it;
325  itBegin = pts.begin();
326  itEnd = pts.end();
327 
328  T1 val1 = bufIn[o1];
329 
330  for (it = itBegin; it != itEnd; it++) {
331  if (x + it->x > width - 1 || x + it->x < 0)
332  continue;
333  if (y + it->y > height - 1 || y + it->y < 0)
334  continue;
335  if (z + it->z > depth - 1 || z + it->z < 0)
336  continue;
337 
338  off_t o2 = o1 + it->z * strideZ + it->y * strideY + it->x;
339  if (o2 <= o1)
340  continue;
341 
342  T2 val2 = bufIn[o2];
343  double valeur = (255.0 / (std::abs((double) (val1 - val2)) + 1));
344  double cost = std::pow(valeur, exposant);
345 
346  bool hasEdge = false;
347  boost::tie(e4, hasEdge) = boost::add_edge(o1, o2, g);
348  boost::tie(e3, hasEdge) = boost::add_edge(o2, o1, g);
349  capacity[e4] = cost;
350  capacity[e3] = cost;
351  rev[e4] = e3;
352  rev[e3] = e4;
353  }
354  }
355  }
356  }
357 
358  for (T2 nbk = 2; nbk <= numLabels; nbk++) {
359  // for all pixels in imIn create a vertex
360  for (off_t o0 = 0; o0 < (off_t) pixelCount; o0++) {
361  T2 val1 = bufMarker[o0];
362  T1 val2 = bufIn[o0];
363 
364  double cost = std::numeric_limits<double>::max();
365  bool hasEdge;
366  if (val1 == nbk) {
367  boost::tie(e4, hasEdge) = boost::edge(vSource, o0, g);
368  if (!hasEdge) {
369  boost::tie(e4, hasEdge) = boost::add_edge(vSource, o0, g);
370  boost::tie(e3, hasEdge) = boost::add_edge(o0, vSource, g);
371  capacity[e4] = cost;
372  capacity[e3] = cost;
373  rev[e4] = e3;
374  rev[e3] = e4;
375  }
376  } else if (val1 > 1 && val1 != nbk) {
377  boost::tie(e4, hasEdge) = boost::edge(val2, vSink, g);
378  if (!hasEdge) {
379  boost::tie(e4, hasEdge) = boost::add_edge(o0, vSink, g);
380  boost::tie(e3, hasEdge) = boost::add_edge(vSink, o0, g);
381  capacity[e4] = cost;
382  capacity[e3] = cost;
383  rev[e4] = e3;
384  rev[e3] = e4;
385  }
386  }
387  }
388 
389  std::cout << "Compute Max flow" << std::endl;
390  boost::property_map<Graph_d, boost::vertex_index_t>::type indexmap =
391  boost::get(boost::vertex_index, g);
392  std::vector<boost::default_color_type> color(boost::num_vertices(g));
393 #if BOOST_VERSION >= 104700
394  double flow =
395  boykov_kolmogorov_max_flow(g, capacity, residual_capacity, rev,
396  &color[0], indexmap, vSource, vSink);
397 #else
398  double flow = kolmogorov_max_flow(g, capacity, residual_capacity, rev,
399  &color[0], indexmap, vSource, vSink);
400 #endif
401  std::cout << "c The total flow:" << std::endl;
402  std::cout << "s " << flow << std::endl << std::endl;
403 
404  // for all pixels in imIn create a vertex and an edge
405  for (off_t o1 = 0; o1 < (off_t) pixelCount; o1++) {
406  T1 val1 = (int) bufIn[o1];
407  T2 val2 = (int) bufOut[o1];
408  T2 val3 = (int) bufMarker[o1];
409 
410  if (val2 == 1) {
411  if (color[val1] == color[vSource])
412  bufOut[o1] = (T2) nbk;
413  }
414 
415  bool hasEdge;
416  if (val3 == nbk) {
417  boost::tie(e4, hasEdge) = boost::edge(vSource, o1, g);
418  if (hasEdge) {
419  boost::remove_edge(vSource, o1, g);
420  boost::remove_edge(o1, vSource, g);
421  }
422  } else if (val3 > 1) {
423  boost::tie(e4, hasEdge) = boost::edge(o1, vSink, g);
424  if (hasEdge) {
425  boost::remove_edge(o1, vSink, g);
426  boost::remove_edge(vSink, o1, g);
427  }
428  }
429  }
430  }
431  return RES_OK;
432  }
433 
434 } // namespace smil
435 
436 #endif // _D_GEOCUTS_WATERSHED_HPP_
smil::BaseImage::getDepth
size_t getDepth() const
Get image depth (Z)
Definition: DBaseImage.h:90
smil::BaseImage::getWidth
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:80
smil::BaseImage::getPixelCount
size_t getPixelCount() const
Get the number of pixels.
Definition: DBaseImage.h:160
smil::BaseImage::getHeight
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:85
smil::Image::getPixels
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:105
smil::StrElt
Base structuring element.
Definition: DStructuringElement.h:68
smil::pow
RES_T pow(const Image< T1 > &imIn, Image< T2 > &imOut, double exponent=2)
pow() - power of an image
Definition: DImageArith.hpp:905
smil::maxVal
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
maxVal() - Max value of an image
Definition: DMeasures.hpp:348
smil::geoCutsMultiway_Watershed
RES_T geoCutsMultiway_Watershed(const Image< T1 > &imIn, const Image< T2 > &imMarker, const double Power, const StrElt &nl, Image< T2 > &imOut)
Watershed as a mutli_terminal cut (multi label)
Definition: Watershed.hpp:248
smil::geoCutsWatershed_MinCut
RES_T geoCutsWatershed_MinCut(const Image< T1 > &imIn, const Image< T2 > &imMarker, const double Power, const StrElt &nl, Image< T2 > &imOut)
Watershed as a Minimum Cut (2 labels)
Definition: Watershed.hpp:89