21 #ifndef mia_internal_watershed_hh
22 #define mia_internal_watershed_hh
35 class TWatershed :
public watershed_traits<dim>::Handler::Product {
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;
49 template <
template <
typename>
class Image,
typename T>
52 struct PixelWithLocation {
59 template <
template <
typename>
class Image,
typename T>
60 bool grow(
const PixelWithLocation& p, Image<unsigned int>& labels,
const Image<T>& data)
const;
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));
68 std::vector<MPosition> m_neighborhood;
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;
96 m_with_borders(with_borders),
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);
106 m_togradnorm = Handler::instance().produce(
"gradnorm");
109 static const unsigned int boundary_label = std::numeric_limits<unsigned int>::max();
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
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;
120 backtrack.push_back(p.pos);
122 std::vector<Position> new_positions;
123 new_positions.reserve(m_neighborhood.size());
125 float value = p.value;
126 unsigned int label = labels(p.pos);
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);
136 while (!locations.empty()) {
138 auto loc = locations.top();
141 new_positions.clear();
143 unsigned int first_label = 0;
144 bool loc_is_boundary =
false;
146 for (
auto i = m_neighborhood.begin(); i != m_neighborhood.end() && !loc_is_boundary; ++i) {
147 Position new_pos( loc + *i);
149 if (! (new_pos < size) )
152 cverb << data(new_pos);
154 if (data(new_pos) > value)
157 unsigned int new_pos_label = labels(new_pos);
158 if (!new_pos_label) {
159 new_positions.push_back(new_pos);
169 first_label = new_pos_label;
170 }
else if (first_label != new_pos_label) {
172 loc_is_boundary =
true;
176 if (!loc_is_boundary) {
177 labels(loc) = first_label;
178 backtrack.push_back(loc);
179 if (first_label != label) {
185 if (!has_backtracked) {
186 for_each(backtrack.begin(), backtrack.end(),
187 [&first_label, &labels](
const Position& p){labels(p) = first_label;});
189 has_backtracked =
true;
197 backtrack.push_back(loc);
201 for_each(new_positions.begin(), new_positions.end(),
202 [&locations](
const Position& p){locations.push(p);});
206 while (!locations.empty() && locations.top() == loc)
210 return has_backtracked;
214 template <
template <
typename>
class Image,
typename T>
217 auto sizeND = data.get_size();
218 Image<unsigned int> labels(data.get_size());
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;
224 std::priority_queue<PixelWithLocation> pixels;
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();
232 p.value = *i > thresh ? *i : thresh;
233 if (p.value <= thresh) {
236 if (!grow(p, labels, data))
245 while (!pixels.empty()) {
246 auto pixel = pixels.top();
249 if (labels(pixel.pos)) {
253 unsigned int first_label = 0;
254 bool is_boundary =
false;
257 for (
auto i = m_neighborhood.begin(); i != m_neighborhood.end() && !is_boundary; ++i) {
259 if (new_pos < sizeND) {
260 auto l = labels(new_pos);
265 if (first_label != l)
266 is_boundary = m_with_borders;
272 labels(pixel.pos) = first_label;
275 cvdebug() <<
" set " << pixel.pos <<
" with " << data(pixel.pos) <<
" to "<< labels(pixel.pos) <<
"\n";
281 labels(pixel.pos) = next_label;
282 if (!grow(pixel, labels, data))
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; });
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(); });
301 Image<unsigned int> * result =
new Image<unsigned int>(data.get_size(), data);
302 copy(labels.begin(), labels.end(), result->begin());
312 return m_togradnorm ?
mia::filter(*
this, *m_togradnorm->filter(image)):
318 watershed_traits<dim>::Handler::Interface(
"ws"),
319 m_with_borders(false),
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",
new mia::CFloatParameter(m_thresh, 0, 1.0,
false,
"Relative gradient norm threshold. The actual value threshhold value "
326 "is thresh * (max_grad - min_grad) + min_grad. Bassins separated by gradients "
327 "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 not represent a gradient norm image"));
332 typename watershed_traits<dim>::Handler::Product *
335 return new TWatershed<dim>(m_neighborhood, m_with_borders, m_thresh, m_eval_grad);
341 return "basic watershead segmentation.";