My Project
seededwatershed.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_seededwatershed_hh
22 #define mia_internal_seededwatershed_hh
23 
24 #include <queue>
25 #include <mia/template/dimtrait.hh>
26 
28 
30 
31 template <int dim>
32 class TSeededWS : public watershed_traits<dim>::Handler::Product
33 {
34 public:
35  typedef typename watershed_traits<dim>::PNeighbourhood PNeighbourhood;
36  typedef typename PNeighbourhood::element_type::value_type MPosition;
37  typedef typename watershed_traits<dim>::Handler Handler;
38  typedef typename watershed_traits<dim>::FileHandler FileHandler;
39  typedef typename Handler::Product CFilter;
40  typedef typename FileHandler::Instance::DataKey DataKey;
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  TSeededWS(const DataKey& mask_image, PNeighbourhood neighborhood,
48  bool with_borders, bool input_is_gradient);
49 
50  template <template <typename> class Image, typename T>
51  typename TSeededWS<dim>::result_type operator () (const Image<T>& data) const;
52 private:
53  virtual PImage do_filter(const CImage& image) const;
54 
55  DataKey m_label_image_key;
56  PNeighbourhood m_neighborhood;
57  PFilter m_togradnorm;
58  bool m_with_borders;
59 };
60 
61 template <int dim>
62 class TSeededWSFilterPlugin: public watershed_traits<dim>::Handler::Interface
63 {
64 public:
65  typedef typename watershed_traits<dim>::PNeighbourhood PNeighbourhood;
66  typedef typename watershed_traits<dim>::Handler Handler;
67  typedef typename watershed_traits<dim>::FileHandler FileHandler;
68  typedef typename Handler::Product CFilter;
69 
70  TSeededWSFilterPlugin();
71 private:
72  virtual CFilter *do_create()const;
73  virtual const std::string do_get_descr()const;
74  std::string m_seed_image_file;
75  PNeighbourhood m_neighborhood;
76  bool m_with_borders;
77  bool m_input_is_gradient;
78 };
79 
80 template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim, bool supported>
81 struct seeded_ws {
82  static R apply(const Image<T>& image, const Image<S>& seed, N n, bool with_borders);
83 };
84 
85 template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim>
86 struct seeded_ws<Image, T, S, N, R, dim, false> {
87  static R apply(const Image<T>& /*image*/, const Image<S>& /*seed*/, N /*n*/, bool /*with_borders*/)
88  {
89  throw create_exception<std::invalid_argument>("C2DRunSeededWS: seed data type '", __type_descr<S>::value, "' not supported");
90  }
91 };
92 
93 
94 // helper to dispatch based on the pixel type of the seed image
95 template <template <typename> class Image, typename T, typename N, typename R, int dim>
96 struct dispatch_RunSeededWS : public TFilter<R> {
97 
98  dispatch_RunSeededWS(N neighborhood, const Image<T>& image, bool with_borders):
99  m_neighborhood(neighborhood),
100  m_image(image),
101  m_with_borders(with_borders)
102  {}
103 
104  template <typename S>
105  R operator () (const Image<S>& seed) const
106  {
107  const bool supported = std::is_integral<S>::value && !std::is_same<S, bool>::value;
108  return seeded_ws<Image, T, S, N, R, dim, supported>::apply(m_image, seed, m_neighborhood, m_with_borders);
109  }
110  N m_neighborhood;
111  const Image<T>& m_image;
112  bool m_with_borders;
113 };
114 
115 
116 template <typename L, typename Position>
117 struct PixelWithLocation {
118  Position pos;
119  float value;
120  L label;
121 };
122 
123 template <typename L, typename Position>
124 bool operator < (const PixelWithLocation<L, Position>& lhs, const PixelWithLocation<L, Position>& rhs)
125 {
126  mia::less_then<Position> l;
127  return lhs.value > rhs.value ||
128  ( lhs.value == rhs.value && l(rhs.pos, lhs.pos));
129 }
130 
131 template <template <typename> class Image, typename T, typename S, int dim>
132 class TRunSeededWatershed
133 {
134 public:
135  typedef typename watershed_traits<dim>::PNeighbourhood PNeighbourhood;
136  typedef typename PNeighbourhood::element_type::value_type MPosition;
137  typedef typename watershed_traits<dim>::Handler Handler;
138  typedef typename Handler::Product CFilter;
139  typedef typename CFilter::Pointer PFilter;
140  typedef typename CFilter::Image CImage;
141  typedef typename CImage::Pointer PImage;
142  typedef typename CImage::dimsize_type Position;
143 
144 
145  TRunSeededWatershed(const Image<T>& image, const Image<S>& seed, PNeighbourhood neighborhood, bool with_borders);
146  PImage run();
147 private:
148  void add_neighborhood(const PixelWithLocation<S, Position>& pixel);
149  const Image<T>& m_image;
150  const Image<S>& m_seed;
151  PNeighbourhood m_neighborhood;
152  Image<unsigned char> m_visited;
153  Image<unsigned char> m_stored;
154  Image<S> *m_result;
155  PImage m_presult;
156  std::priority_queue<PixelWithLocation<S, Position>> m_seeds;
157  S m_watershed;
158  bool m_with_borders;
159 };
160 
161 template <template <typename> class Image, typename T, typename S, int dim>
162 TRunSeededWatershed<Image, T, S, dim>::TRunSeededWatershed(const Image<T>& image, const Image<S>& seed,
163  PNeighbourhood neighborhood, bool with_borders):
164  m_image(image),
165  m_seed(seed),
166  m_neighborhood(neighborhood),
167  m_visited(seed.get_size()),
168  m_stored(seed.get_size()),
169  m_watershed(std::numeric_limits<S>::max()),
170  m_with_borders(with_borders)
171 {
172  m_result = new Image<S>(seed.get_size(), image);
173  m_presult.reset(m_result);
174 }
175 
176 template <template <typename> class Image, typename T, typename S, int dim>
177 void TRunSeededWatershed<Image, T, S, dim>::add_neighborhood(const PixelWithLocation<S, Position>& pixel)
178 {
179  PixelWithLocation<S, Position> new_pixel;
180  new_pixel.label = pixel.label;
181  bool hit_boundary = false;
182 
183  for (auto i = m_neighborhood->begin(); i != m_neighborhood->end(); ++i) {
184  Position new_pos( pixel.pos + *i);
185 
186  if (new_pos != pixel.pos && new_pos < m_seed.get_size()) {
187  if (!m_visited(new_pos)) {
188  if (!m_stored(new_pos) ) {
189  new_pixel.value = m_image(new_pos);
190  new_pixel.pos = new_pos;
191  m_seeds.push(new_pixel);
192  m_stored(new_pos) = 1;
193  }
194  } else {
195  hit_boundary |= (*m_result)(new_pos) != pixel.label &&
196  (*m_result)(new_pos) != m_watershed;
197  }
198  }
199  }
200 
201  // set pixel to new label
202  if (!m_visited(pixel.pos)) {
203  m_visited(pixel.pos) = true;
204  (*m_result)(pixel.pos) = (m_with_borders && hit_boundary) ? m_watershed : pixel.label;
205  }
206 }
207 
208 template <template <typename> class Image, typename T, typename S, int dim>
209 typename TRunSeededWatershed<Image, T, S, dim>::PImage
210 TRunSeededWatershed<Image, T, S, dim>::run()
211 {
212  // copy seed and read initial pixels
213  auto iv = m_visited.begin();
214  auto is = m_seed.begin_range(Position::_0, m_seed.get_size());
215  auto es = m_seed.end_range(Position::_0, m_seed.get_size());
216  auto ir = m_result->begin();
217  auto ii = m_image.begin();
218  auto ist = m_stored.begin();
219  PixelWithLocation<S, Position> pixel;
220 
221  while (is != es) {
222  *ir = *is;
223 
224  if (*is) {
225  *iv = 1;
226  *ist;
227  pixel.pos = is.pos();
228  pixel.value = *ii;
229  pixel.label = *is;
230  m_seeds.push(pixel);
231  }
232 
233  ++iv;
234  ++is;
235  ++ir;
236  ++ist;
237  ++ii;
238  }
239 
240  while (!m_seeds.empty()) {
241  auto p = m_seeds.top();
242  m_seeds.pop();
243  add_neighborhood(p);
244  }
245 
246  return m_presult;
247 }
248 
249 template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim, bool supported>
250 R seeded_ws<Image, T, S, N, R, dim, supported>::apply(const Image<T>& image,
251  const Image<S>& seed, N neighborhood, bool with_borders)
252 {
253  TRunSeededWatershed<Image, T, S, dim> ws(image, seed, neighborhood, with_borders);
254  return ws.run();
255 }
256 
257 template <int dim>
258 TSeededWS<dim>::TSeededWS(const DataKey& label_image_key, PNeighbourhood neighborhood, bool with_borders, bool input_is_gradient):
259  m_label_image_key(label_image_key),
260  m_neighborhood(neighborhood),
261  m_with_borders(with_borders)
262 
263 {
264  if (!input_is_gradient)
265  m_togradnorm = Handler::instance().produce("gradnorm");
266 }
267 
268 template <int dim>
269 template <template <typename> class Image, typename T>
270 typename TSeededWS<dim>::result_type TSeededWS<dim>::operator () (const Image<T>& data) const
271 {
272  // read start label image
273  auto in_image_list = m_label_image_key.get();
274 
275  if (!in_image_list || in_image_list->empty())
276  throw create_exception<std::runtime_error>( "C2DSeededWS: no seed image could be loaded");
277 
278  if (in_image_list->size() > 1)
279  cvwarn() << "C2DSeededWS:got more than one seed image. Ignoring all but first";
280 
281  auto seed = (*in_image_list)[0];
282 
283  if (seed->get_size() != data.get_size()) {
284  throw create_exception<std::invalid_argument>( "C2DSeededWS: seed and input differ in size: seed ", seed->get_size()
285  , ", input ", data.get_size());
286  }
287 
288  dispatch_RunSeededWS<Image, T, PNeighbourhood, PImage, dim> ws(m_neighborhood, data, m_with_borders);
289  return mia::filter(ws, *seed);
290 }
291 
292 template <int dim>
293 typename TSeededWS<dim>::PImage TSeededWS<dim>::do_filter(const CImage& image) const
294 {
295  PImage result;
296 
297  if (m_togradnorm) {
298  auto grad = m_togradnorm->filter(image);
299  result = mia::filter(*this, *grad);
300  } else
301  result = mia::filter(*this, image);
302 
303  return result;
304 }
305 
306 template <int dim>
307 TSeededWSFilterPlugin<dim>::TSeededWSFilterPlugin():
308  Handler::Interface("sws"),
309  m_with_borders(false),
310  m_input_is_gradient(false)
311 {
312  this->add_parameter("seed", new CStringParameter(m_seed_image_file, CCmdOptionFlags::required_input,
313  "seed input image containing the lables for the initial regions"));
314  this->add_parameter("n", make_param(m_neighborhood, "sphere:r=1", false, "Neighborhood for watershead region growing"));
315  this->add_parameter("mark", new CBoolParameter(m_with_borders, false, "Mark the segmented watersheds with a special gray scale value"));
316  this->add_parameter("grad", new CBoolParameter(m_input_is_gradient, false, "Interpret the input image as gradient. "));
317 }
318 
319 template <int dim>
320 typename TSeededWSFilterPlugin<dim>::CFilter *TSeededWSFilterPlugin<dim>::do_create()const
321 {
322  auto seed = FileHandler::instance().load_to_pool(m_seed_image_file);
323  return new TSeededWS<dim>(seed, m_neighborhood, m_with_borders, m_input_is_gradient);
324 }
325 
326 template <int dim>
327 const std::string TSeededWSFilterPlugin<dim>::do_get_descr()const
328 {
329  return "seeded watershead. The algorithm extracts exactly so "
330  "many reagions as initial labels are given in the seed image.";
331 }
332 
334 
336 #endif
bool operator<(const T2DVector< T > &a, const T2DVector< S > &b)
Definition: 2d/vector.hh:495
an string parameter
Definition: parameter.hh:537
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
static F::result_type filter(const F &f, const B &b)
Definition: core/filter.hh:258
vstream & cvwarn()
send warnings to this stream adapter
Definition: msgstream.hh:321
CParameter * make_param(T &value, bool required, const char *descr)
Definition: parameter.hh:280
CTParameter< bool > CBoolParameter
boolean parameter
Definition: parameter.hh:561
base class for all filer type functors.
Definition: core/filter.hh:70