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>
26
28
30
31template <int dim>
32class TSeededWS : public watershed_traits<dim>::Handler::Product
33{
34public:
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;
52private:
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
61template <int dim>
62class TSeededWSFilterPlugin: public watershed_traits<dim>::Handler::Interface
63{
64public:
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();
71private:
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
80template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim, bool supported>
81struct seeded_ws {
82 static R apply(const Image<T>& image, const Image<S>& seed, N n, bool with_borders);
83};
84
85template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim>
86struct 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
95template <template <typename> class Image, typename T, typename N, typename R, int dim>
96struct 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
116template <typename L, typename Position>
117struct PixelWithLocation {
118 Position pos;
119 float value;
120 L label;
121};
122
123template <typename L, typename Position>
124bool 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
131template <template <typename> class Image, typename T, typename S, int dim>
132class TRunSeededWatershed
133{
134public:
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();
147private:
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
161template <template <typename> class Image, typename T, typename S, int dim>
162TRunSeededWatershed<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
176template <template <typename> class Image, typename T, typename S, int dim>
177void 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
208template <template <typename> class Image, typename T, typename S, int dim>
209typename TRunSeededWatershed<Image, T, S, dim>::PImage
210TRunSeededWatershed<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
249template <template <typename> class Image, typename T, typename S, typename N, typename R, int dim, bool supported>
250R 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
257template <int dim>
258TSeededWS<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
268template <int dim>
269template <template <typename> class Image, typename T>
270typename 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
292template <int dim>
293typename 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
306template <int dim>
307TSeededWSFilterPlugin<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
319template <int dim>
320typename 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
326template <int dim>
327const 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
CTParameter< bool > CBoolParameter
boolean parameter
Definition: parameter.hh:561
CParameter * make_param(T &value, bool required, const char *descr)
Definition: parameter.hh:280
base class for all filer type functors.
Definition: core/filter.hh:70