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
34template <int dim>
35class TWatershed : public watershed_traits<dim>::Handler::Product
36{
37public:
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 ;
52private:
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
82template <int dim>
83class TWatershedFilterPlugin: public watershed_traits<dim>::Handler::Interface
84{
85public:
87private:
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
97template <int dim>
98TWatershed<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
113static const unsigned int boundary_label = std::numeric_limits<unsigned int>::max();
114
115template <int dim>
116template <template <typename> class Image, typename T>
117bool 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
219template <int dim>
220template <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
325template <int dim>
326typename 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
332template <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
348template <int dim>
349typename watershed_traits<dim>::Handler::Product *
351{
352 return new TWatershed<dim>(m_neighborhood, m_with_borders, m_thresh, m_eval_grad);
353}
354
355template <int dim>
356const std::string TWatershedFilterPlugin<dim>::do_get_descr()const
357{
358 return "basic watershead segmentation.";
359}
360
362
363#endif
plugin for the templated version of the standard watershed algorithm
Definition watershed.hh:84
templated version of the standard watershed algorithm
Definition watershed.hh:36
watershed_traits< dim >::PNeighbourhood PNeighbourhood
Definition watershed.hh:38
TWatershed< dim >::result_type operator()(const Image< T > &data) const
Definition watershed.hh:221
CImage::Pointer PImage
Definition watershed.hh:44
watershed_traits< dim >::Handler Handler
Definition watershed.hh:40
TWatershed(PNeighbourhood neighborhood, bool with_borders, float treash, bool eval_grad)
Definition watershed.hh:98
PNeighbourhood::element_type::value_type MPosition
Definition watershed.hh:39
CFilter::Image CImage
Definition watershed.hh:43
Handler::Product CFilter
Definition watershed.hh:41
CImage::dimsize_type Position
Definition watershed.hh:45
CFilter::Pointer PFilter
Definition watershed.hh:42
friend bool operator<(const PixelWithLocation &lhs, const PixelWithLocation &rhs)
Definition watershed.hh:63
#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)
#define cverb
define a shortcut to the raw output stream
Definition msgstream.hh:341
vstream & cvmsg()
send messages to this stream adapter
Definition msgstream.hh:331
CDebugSink & cvdebug()
Definition msgstream.hh:226
CParameter * make_coi_param(T &value, S1 lower_bound, S2 upper_bound, bool required, const char *descr)
Definition parameter.hh:343
CParameter * make_param(T &value, bool required, const char *descr)
Definition parameter.hh:280
static const unsigned int boundary_label
Definition watershed.hh:113