My Project
watershed.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef mia_internal_watershed_hh
22 #define mia_internal_watershed_hh
23 
24 #include <mia/core/filter.hh>
25 #include <queue>
26 
28 
34 template <int dim>
35 class TWatershed : public watershed_traits<dim>::Handler::Product
36 {
37 public:
38  typedef typename watershed_traits<dim>::PNeighbourhood PNeighbourhood;
39  typedef typename PNeighbourhood::element_type::value_type MPosition;
40  typedef typename watershed_traits<dim>::Handler Handler;
41  typedef typename Handler::Product CFilter;
42  typedef typename CFilter::Pointer PFilter;
43  typedef typename CFilter::Image CImage;
44  typedef typename CImage::Pointer PImage;
45  typedef typename CImage::dimsize_type Position;
46 
47 
48  TWatershed(PNeighbourhood neighborhood, bool with_borders, float treash, bool eval_grad);
49 
50  template <template <typename> class Image, typename T>
51  typename TWatershed<dim>::result_type operator () (const Image<T>& data) const ;
52 private:
53  struct PixelWithLocation {
54  Position pos;
55  float value;
56  };
57 
58  typename TWatershed<dim>::result_type do_filter(const CImage& image) const;
59 
60  template <template <typename> class Image, typename T>
61  bool grow(const PixelWithLocation& p, Image<unsigned int>& labels, const Image<T>& data) const;
62 
63  friend bool operator < (const PixelWithLocation& lhs, const PixelWithLocation& rhs)
64  {
65  mia::less_then<Position> l;
66  return lhs.value > rhs.value ||
67  ( lhs.value == rhs.value && l(rhs.pos, lhs.pos));
68  }
69 
70  std::vector<MPosition> m_neighborhood;
71  PFilter m_togradnorm;
72  bool m_with_borders;
73  float m_thresh;
74 };
75 
76 
82 template <int dim>
83 class TWatershedFilterPlugin: public watershed_traits<dim>::Handler::Interface
84 {
85 public:
87 private:
88  virtual typename watershed_traits<dim>::Handler::Product *do_create()const;
89  virtual const std::string do_get_descr()const;
90  typename watershed_traits<dim>::PNeighbourhood m_neighborhood;
91  bool m_with_borders;
92  float m_thresh;
93  bool m_eval_grad;
94 };
95 
96 
97 template <int dim>
98 TWatershed<dim>::TWatershed(PNeighbourhood neighborhood, bool with_borders, float thresh, bool eval_grad):
99  m_with_borders(with_borders),
100  m_thresh(thresh)
101 
102 {
103  m_neighborhood.reserve(neighborhood->size() - 1);
104 
105  for (auto i = neighborhood->begin(); i != neighborhood->end(); ++i)
106  if ( *i != MPosition::_0 )
107  m_neighborhood.push_back(*i);
108 
109  if (eval_grad)
110  m_togradnorm = Handler::instance().produce("gradnorm");
111 }
112 
113 static const unsigned int boundary_label = std::numeric_limits<unsigned int>::max();
114 
115 template <int dim>
116 template <template <typename> class Image, typename T>
117 bool TWatershed<dim>::grow(const PixelWithLocation& p, Image<unsigned int>& labels, const Image<T>& data) const
118 {
119  const auto size = data.get_size();
120  std::vector<Position> backtrack;
121  std::priority_queue<Position, std::vector<Position>, mia::less_then<Position>> locations;
122  bool has_backtracked = false;
123  backtrack.push_back(p.pos);
124  std::vector<Position> new_positions;
125  new_positions.reserve(m_neighborhood.size());
126  float value = p.value;
127  unsigned int label = labels(p.pos);
128 
129  for (auto i = m_neighborhood.begin(); i != m_neighborhood.end(); ++i) {
130  Position new_pos( p.pos + *i);
131 
132  if (new_pos < size && !labels(new_pos) && data(new_pos) <= value) {
133  locations.push(new_pos);
134  backtrack.push_back(new_pos);
135  }
136  }
137 
138  while (!locations.empty()) {
139  // incoming locations are always un-labelled, and the gradient value is equal or below the target value
140  auto loc = locations.top();
141  locations.pop();
142  new_positions.clear();
143  unsigned int first_label = 0;
144  bool loc_is_boundary = false;
145 
146  for (auto i = m_neighborhood.begin(); i != m_neighborhood.end() && !loc_is_boundary; ++i) {
147  Position new_pos( loc + *i);
148 
149  if (! (new_pos < size) )
150  continue;
151 
152  cverb << data(new_pos);
153 
154  if (data(new_pos) > value)
155  continue;
156 
157  unsigned int new_pos_label = labels(new_pos);
158 
159  if (!new_pos_label) {
160  new_positions.push_back(new_pos);
161  continue;
162  }
163 
164  // already visited?
165  if (new_pos_label == label || new_pos_label == boundary_label)
166  continue;
167 
168  // first label hit
169  if (!first_label) {
170  first_label = new_pos_label;
171  } else if (first_label != new_pos_label) {
172  // hit two different labels (apart from the original one)
173  loc_is_boundary = true;
174  }
175  }
176 
177  if (first_label) {
178  if (!loc_is_boundary) {
179  labels(loc) = first_label;
180  backtrack.push_back(loc);
181 
182  if (first_label != label) {
183  // we hit a single label from a lower gradient value, this means
184  // we are connected to an already labeled basin ->
185  // first time = backtrack
186  // later = boundary
187  if (!has_backtracked) {
188  for_each(backtrack.begin(), backtrack.end(),
189  [&first_label, &labels](const Position & p) {
190  labels(p) = first_label;
191  });
192  label = first_label;
193  has_backtracked = true;
194  } else
195  labels(loc) = m_with_borders ? boundary_label : label;
196  }
197  } else
198  labels(loc) = m_with_borders ? boundary_label : label;
199  } else {
200  labels(loc) = label;
201  backtrack.push_back(loc);
202  }
203 
204  if (labels(loc) != boundary_label) {
205  for_each(new_positions.begin(), new_positions.end(),
206  [&locations](const Position & p) {
207  locations.push(p);
208  });
209  }
210 
211  // is there a queue that doesn't repeat values?
212  while (!locations.empty() && locations.top() == loc)
213  locations.pop();
214  }
215 
216  return has_backtracked;
217 }
218 
219 template <int dim>
220 template <template <typename> class Image, typename T>
222 {
223  auto sizeND = data.get_size();
224  Image<unsigned int> labels(data.get_size());
225  // evaluate the real thresh hold based on the actual gradient range
226  auto gradient_range = std::minmax_element(data.begin(), data.end());
227  float thresh = m_thresh * (*gradient_range.second - *gradient_range.first) + *gradient_range.first;
228  std::priority_queue<PixelWithLocation> pixels;
229  PixelWithLocation p;
230  auto i = data.begin_range(Position::_0, data.get_size());
231  auto e = data.end_range(Position::_0, data.get_size());
232  auto l = labels.begin();
233  long next_label = 1;
234 
235  while (i != e) {
236  p.pos = i.pos();
237  p.value = *i > thresh ? *i : thresh;
238 
239  if (p.value <= thresh) {
240  if (!*l) {
241  *l = next_label;
242 
243  if (!grow(p, labels, data))
244  ++next_label;
245  }
246  } else
247  pixels.push(p);
248 
249  ++i;
250  ++l;
251  }
252 
253  while (!pixels.empty()) {
254  auto pixel = pixels.top();
255  pixels.pop();
256 
257  // this label was set because we grew an initial region
258  if (labels(pixel.pos)) {
259  continue;
260  }
261 
262  unsigned int first_label = 0;
263  bool is_boundary = false;
264  // check if neighborhood is already labeled
265 
266  for (auto i = m_neighborhood.begin(); i != m_neighborhood.end() && !is_boundary; ++i) {
267  Position new_pos( pixel.pos + *i);
268 
269  if (new_pos < sizeND) {
270  auto l = labels(new_pos);
271 
272  if ( l && l != boundary_label) {
273  if (!first_label)
274  first_label = l;
275  else if (first_label != l)
276  is_boundary = m_with_borders;
277  }
278  }
279  }
280 
281  if (first_label) {
282  if (!is_boundary)
283  labels(pixel.pos) = first_label;
284  else
285  labels(pixel.pos) = boundary_label;
286 
287  cvdebug() << " set " << pixel.pos << " with " << data(pixel.pos) << " to " << labels(pixel.pos) << "\n";
288  continue;
289  }
290 
291  // new label to assign
292  // if a new label is assigned, we have to grow the region of equal gradient values
293  // to assure we catch the whole bassin
294  labels(pixel.pos) = next_label;
295 
296  if (!grow(pixel, labels, data))
297  ++next_label;
298  }
299 
300  // convert to smalles possible intensity range and convert the boundary label to highest
301  // intensity value
302  CImage *r = NULL;
303  cvmsg() << "Got " << next_label << " distinct bassins\n";
304 
305  if (next_label < 255) {
306  Image<unsigned char> *result = new Image<unsigned char>(data.get_size(), data);
307  std::transform(labels.begin(), labels.end(), result->begin(),
308  [](unsigned int p)-> unsigned char { return (p != boundary_label) ? static_cast<unsigned char>(p) : 255; });
309  r = result;
310  } else if (next_label < std::numeric_limits<unsigned short>::max()) {
311  Image<unsigned short> *result = new Image<unsigned short>(data.get_size(), data);
312  std::transform(labels.begin(), labels.end(), result->begin(),
313  [](unsigned int p)-> unsigned short { return (p != boundary_label) ? static_cast<unsigned short>(p) :
314  std::numeric_limits<unsigned short>::max(); });
315  r = result;
316  } else {
317  Image<unsigned int> *result = new Image<unsigned int>(data.get_size(), data);
318  copy(labels.begin(), labels.end(), result->begin());
319  r = result;
320  }
321 
322  return PImage(r);
323 }
324 
325 template <int dim>
326 typename TWatershed<dim>::result_type TWatershed<dim>::do_filter(const CImage& image) const
327 {
328  return m_togradnorm ? mia::filter(*this, *m_togradnorm->filter(image)) :
329  mia::filter(*this, image);
330 }
331 
332 template <int dim>
334  watershed_traits<dim>::Handler::Interface("ws"),
335  m_with_borders(false),
336  m_thresh(0.0),
337  m_eval_grad(false)
338 {
339  this->add_parameter("n", make_param(m_neighborhood, "sphere:r=1", false, "Neighborhood for watershead region growing"));
340  this->add_parameter("mark", new mia::CBoolParameter(m_with_borders, false, "Mark the segmented watersheds with a special gray scale value"));
341  this->add_parameter("thresh", make_coi_param(m_thresh, 0, 1.0, false, "Relative gradient norm threshold. The actual value "
342  "threshold value is thresh * (max_grad - min_grad) + min_grad. Bassins "
343  "separated by gradients with a lower norm will be joined"));
344  this->add_parameter("evalgrad", new mia::CBoolParameter(m_eval_grad, false, "Set to 1 if the input image does "
345  "not represent a gradient norm image"));
346 }
347 
348 template <int dim>
349 typename watershed_traits<dim>::Handler::Product *
351 {
352  return new TWatershed<dim>(m_neighborhood, m_with_borders, m_thresh, m_eval_grad);
353 }
354 
355 template <int dim>
356 const std::string TWatershedFilterPlugin<dim>::do_get_descr()const
357 {
358  return "basic watershead segmentation.";
359 }
360 
362 
363 #endif
TWatershedFilterPlugin
plugin for the templated version of the standard watershed algorithm
Definition: watershed.hh:83
TWatershed::MPosition
PNeighbourhood::element_type::value_type MPosition
Definition: watershed.hh:39
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
TWatershed::PFilter
CFilter::Pointer PFilter
Definition: watershed.hh:42
TWatershed::PNeighbourhood
watershed_traits< dim >::PNeighbourhood PNeighbourhood
Definition: watershed.hh:38
TWatershed::CFilter
Handler::Product CFilter
Definition: watershed.hh:41
make_coi_param
CParameter * make_coi_param(T &value, S1 lower_bound, S2 upper_bound, bool required, const char *descr)
Definition: parameter.hh:343
TWatershed
templated version of the standard watershed algorithm
Definition: watershed.hh:35
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
make_param
CParameter * make_param(T &value, bool required, const char *descr)
Definition: parameter.hh:280
filter
static F::result_type filter(const F &f, const B &b)
Definition: core/filter.hh:258
cvdebug
CDebugSink & cvdebug()
Definition: msgstream.hh:226
TWatershedFilterPlugin::TWatershedFilterPlugin
TWatershedFilterPlugin()
Definition: watershed.hh:333
TWatershed::PImage
CImage::Pointer PImage
Definition: watershed.hh:44
TWatershed::operator()
TWatershed< dim >::result_type operator()(const Image< T > &data) const
Definition: watershed.hh:221
boundary_label
static const unsigned int boundary_label
Definition: watershed.hh:113
cvmsg
vstream & cvmsg()
send messages to this stream adapter
Definition: msgstream.hh:331
TWatershed::Position
CImage::dimsize_type Position
Definition: watershed.hh:45
cverb
#define cverb
define a shortcut to the raw output stream
Definition: msgstream.hh:341
TWatershed::Handler
watershed_traits< dim >::Handler Handler
Definition: watershed.hh:40
CBoolParameter
CTParameter< bool > CBoolParameter
boolean parameter
Definition: parameter.hh:561
TWatershed::CImage
CFilter::Image CImage
Definition: watershed.hh:43
filter.hh
TWatershed::TWatershed
TWatershed(PNeighbourhood neighborhood, bool with_borders, float treash, bool eval_grad)
Definition: watershed.hh:98
TWatershed::operator<
friend bool operator<(const PixelWithLocation &lhs, const PixelWithLocation &rhs)
Definition: watershed.hh:63