21 #ifndef mia_internal_watershed_hh
22 #define mia_internal_watershed_hh
35 class TWatershed :
public watershed_traits<dim>::Handler::Product
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;
50 template <
template <
typename>
class Image,
typename T>
53 struct PixelWithLocation {
60 template <
template <
typename>
class Image,
typename T>
61 bool grow(
const PixelWithLocation& p, Image<unsigned int>& labels,
const Image<T>& data)
const;
63 friend bool operator < (
const PixelWithLocation& lhs,
const PixelWithLocation& rhs)
65 mia::less_then<Position> l;
66 return lhs.value > rhs.value ||
67 ( lhs.value == rhs.value && l(rhs.pos, lhs.pos));
70 std::vector<MPosition> m_neighborhood;
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;
99 m_with_borders(with_borders),
103 m_neighborhood.reserve(neighborhood->size() - 1);
105 for (
auto i = neighborhood->begin(); i != neighborhood->end(); ++i)
106 if ( *i != MPosition::_0 )
107 m_neighborhood.push_back(*i);
110 m_togradnorm = Handler::instance().produce(
"gradnorm");
113 static const unsigned int boundary_label = std::numeric_limits<unsigned int>::max();
116 template <
template <
typename>
class Image,
typename T>
117 bool TWatershed<dim>::grow(
const PixelWithLocation& p, Image<unsigned int>& labels,
const Image<T>& data)
const
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);
129 for (
auto i = m_neighborhood.begin(); i != m_neighborhood.end(); ++i) {
130 Position new_pos( p.pos + *i);
132 if (new_pos < size && !labels(new_pos) && data(new_pos) <= value) {
133 locations.push(new_pos);
134 backtrack.push_back(new_pos);
138 while (!locations.empty()) {
140 auto loc = locations.top();
142 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);
159 if (!new_pos_label) {
160 new_positions.push_back(new_pos);
170 first_label = new_pos_label;
171 }
else if (first_label != new_pos_label) {
173 loc_is_boundary =
true;
178 if (!loc_is_boundary) {
179 labels(loc) = first_label;
180 backtrack.push_back(loc);
182 if (first_label != label) {
187 if (!has_backtracked) {
188 for_each(backtrack.begin(), backtrack.end(),
189 [&first_label, &labels](
const Position & p) {
190 labels(p) = first_label;
193 has_backtracked =
true;
201 backtrack.push_back(loc);
205 for_each(new_positions.begin(), new_positions.end(),
206 [&locations](
const Position & p) {
212 while (!locations.empty() && locations.top() == loc)
216 return has_backtracked;
220 template <
template <
typename>
class Image,
typename T>
223 auto sizeND = data.get_size();
224 Image<unsigned int> labels(data.get_size());
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;
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();
237 p.value = *i > thresh ? *i : thresh;
239 if (p.value <= thresh) {
243 if (!grow(p, labels, data))
253 while (!pixels.empty()) {
254 auto pixel = pixels.top();
258 if (labels(pixel.pos)) {
262 unsigned int first_label = 0;
263 bool is_boundary =
false;
266 for (
auto i = m_neighborhood.begin(); i != m_neighborhood.end() && !is_boundary; ++i) {
269 if (new_pos < sizeND) {
270 auto l = labels(new_pos);
275 else if (first_label != l)
276 is_boundary = m_with_borders;
283 labels(pixel.pos) = first_label;
287 cvdebug() <<
" set " << pixel.pos <<
" with " << data(pixel.pos) <<
" to " << labels(pixel.pos) <<
"\n";
294 labels(pixel.pos) = next_label;
296 if (!grow(pixel, labels, data))
303 cvmsg() <<
"Got " << next_label <<
" distinct bassins\n";
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; });
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(); });
317 Image<unsigned int> *result =
new Image<unsigned int>(data.get_size(), data);
318 copy(labels.begin(), labels.end(), result->begin());
328 return m_togradnorm ?
mia::filter(*
this, *m_togradnorm->filter(image)) :
334 watershed_traits<dim>::Handler::Interface(
"ws"),
335 m_with_borders(false),
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"));
349 typename watershed_traits<dim>::Handler::Product *
352 return new TWatershed<dim>(m_neighborhood, m_with_borders, m_thresh, m_eval_grad);
358 return "basic watershead segmentation.";
plugin for the templated version of the standard watershed algorithm
templated version of the standard watershed algorithm
watershed_traits< dim >::PNeighbourhood PNeighbourhood
TWatershed< dim >::result_type operator()(const Image< T > &data) const
watershed_traits< dim >::Handler Handler
TWatershed(PNeighbourhood neighborhood, bool with_borders, float treash, bool eval_grad)
PNeighbourhood::element_type::value_type MPosition
CImage::dimsize_type Position
friend bool operator<(const PixelWithLocation &lhs, const PixelWithLocation &rhs)
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
#define NS_MIA_END
conveniance define to end the mia namespace
static F::result_type filter(const F &f, const B &b)
vstream & cvmsg()
send messages to this stream adapter
#define cverb
define a shortcut to the raw output stream
CParameter * make_coi_param(T &value, S1 lower_bound, S2 upper_bound, bool required, const char *descr)
CParameter * make_param(T &value, bool required, const char *descr)
CTParameter< bool > CBoolParameter
boolean parameter
static const unsigned int boundary_label