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