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