Commit 582356c0 authored by Michaël Roynard's avatar Michaël Roynard

Merge branch 'development/wst' into 'dev'

Meyer watershed segmentation.

See merge request !12
parents 5cd131c0 6e1e6aaf
Labeling Module
===============
.. toctree::
:maxdepth: 1
labeling/local_extrema
Local extrema labelization
==========================
Include :file:`<mln/labeling/local_extrema.hpp>`
.. cpp:namespace:: mln::labeling
#. .. cpp:function:: \
template <class Label_t, class I, class N> \
ch_value_t<I, Label_t> local_minima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel)
#. .. cpp:function:: \
template <class Label_t, class I, class N> \
ch_value_t<I, Label_t> local_maxima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel)
#. .. cpp:function:: \
template <class Label_t, class I, class N, class Compare> \
ch_value_t<I, Label_t> local_minima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel, Compare cmp)
Labelize the local minima (1) or maxima (2) of an image with an optional
comparaison function ``cmp`` (3).
:tparam Label_t: The type of label (must be :cpp:concept:`Integral` and
not boolean)
:param input: Input image
:param nbh: Neighborhood considered
:param nlabel (out): Number of extrema detected
:param cmp (optional): Comparison function
:return: A labelized image
:exception: N/A
Notes
-----
Complexity
----------
The algorithm is quasi-linear and requires :math:`n` extra-memory space.
......@@ -14,3 +14,16 @@ Structural Morphological Operation
morpho/hit_or_miss
morpho/opening_by_reconstruction
morpho/closing_by_reconstruction
Segmentation
************
.. toctree::
:maxdepth: 1
morpho/watershed
Closing
=======
Include :file:`<mln/morpho/structural/closing.hpp>`
.. cpp:namespace:: mln::morpho::structural
#. .. cpp:function:: \
template <class InputImage, class StructuringElement> \
concrete_t<InputImage> closing(const InputImage& ima, const StructuringElement& se)
......
Closing by reconstruction
=========================
Include :file:`<mln/morpho/closing_by_reconstruction.hpp>`
.. cpp:namespace:: mln::morpho
.. cpp:function:: \
template <class InputImage1, class InputImage2, class Neighborhood> \
concrete_t<InputImage1> closing_by_reconstruction(const InputImage1& f, const InputImage2& markers, const Neighborhood& nbh)
......
.. _dilation:
Dilation
========
Include :file:`<mln/morpho/structural/dilate.hpp>`
.. cpp:namespace:: mln::morpho::structural
#. .. cpp:function:: \
template <class InputImage, class StructuringElement> \
concrete_t<InputImage> dilate(const InputImage& ima, const StructuringElement& se)
......
Erosion
=======
Include :file:`<mln/morpho/structural/erode.hpp>`
#. .. cpp:function:: \
template <class InputImage, class StructuringElement> \
concrete_t<InputImage> erode(const InputImage& ima, const StructuringElement& se)
......
Opening
=======
Include :file:`<mln/morpho/structural/opening.hpp>`
.. cpp:namespace:: mln::morpho::structural
#. .. cpp:function:: \
template <class InputImage, class StructuringElement> \
concrete_t<InputImage> opening(const InputImage& ima, const StructuringElement& se)
......
......@@ -3,6 +3,10 @@
Opening by reconstruction
=========================
Include :file:`<mln/morpho/opening_by_reconstruction.hpp>`
.. cpp:namespace:: mln::morpho
.. cpp:function:: \
template <class InputImage1, class InputImage2, class Neighborhood> \
concrete_t<InputImage1> opening_by_reconstruction(const InputImage1& f, const InputImage2& markers, const Neighborhood& nbh)
......
Watershed
=========
Include :file:`<mln/morpho/watershed.hpp>`
.. cpp:namespace:: mln::morpho
.. cpp:function:: template <class Label_t, class I, class N> \
ch_value_t<I, Label_t> watershed(const Image<I>& ima, const Neighborhood<N>& nbh, int& nlabel)
Watershed by immersion as defined in [BM92]_. The catchment basins are
labeled from 1 to n, and the special label 0 is used for watershed lines.
:tparam Label_t: The type of label (must be *signed* :cpp:concept:`Integral`)
:param input: Input image
:param nbh: Neighborhood considered
:param nlabel (out): Number of catchment basins
:return: A labelized image
:exception: N/A
Notes
-----
Complexity
----------
The algorithm is quasi-linear and requires :math:`n` extra-memory space.
References
----------
.. [BM92] Beucher, S., & Meyer, F. (1992). The morphological approach to
segmentation: the watershed transformation. Optical Engineering-New York-Marcel
Dekker Incorporated-, 34, 433-433.
Example 1: Cell segmentation
----------------------------
1. The distance transform is performed. Maxima correspond to cell centers.
A :ref:`dilation <dilation>` by a small disc removes the non-meaningfull maxima.
2. Invertion of the distance image so that maxima become minima
3. Watershed segmentation
4. Input labelization w.r.t in segmentation labels
.. literalinclude:: /snipsets/blobs_watershed.cpp
:start-after: BEGIN
:end-before: END
.. list-table::
* - .. image:: /images/blobs_binary.png
:width: 100%
- .. image:: /images/blobs_distance_transform.png
:width: 100%
* - Input image (note that blobs may touch)
- Distance transform and dilation (after heat LUT)
.. figure:: /images/blobs_segmentation.png
:width: 100%
:align: center
Segmented blobs and watershed lines (labels displayed in false colors).
......@@ -3,9 +3,12 @@ Reference
.. toctree::
:maxdepth: 2
concept
images
core
algorithms
labeling
morpho
......@@ -28,6 +28,7 @@ add_image("erosion-cli;opening;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_
add_image("erosion-cli;closing;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_closing_1.png)
add_image("staff_lines" "${DOCUMENTATION_IMAGE_DIR}/staff_lines.pbm"
morpho_hitormiss_1.png morpho_hitormiss_2.png staff_lines_markers.png morpho_reconstruction_1.png morpho_reconstruction_2.png)
add_image("blobs_watershed" "${DOCUMENTATION_IMAGE_DIR}/blobs_binary.png" blobs_distance_transform.png blobs_segmentation.png)
add_custom_target(build-images
......@@ -38,4 +39,5 @@ link_libraries(${Boost_PROGRAM_OPTIONS_LIBRARY_RELEASE} ${FreeImage_LIBRARY})
add_executable(erosion-cli erosion-cli.cpp)
add_executable(staff_lines staff_lines.cpp)
add_executable(component_tree_1 component_tree_1.cpp)
add_executable(blobs_watershed blobs_watershed.cpp)
#include <mln/core/image/image2d.hpp>
#include <mln/core/neighb2d.hpp>
#include <mln/core/se/disc.hpp>
#include <mln/core/algorithm/transform.hpp>
#include <mln/core/neighborhood/dyn_wneighborhood.hpp>
#include <mln/data/stretch.hpp>
#include <mln/transform/chamfer_distance_transform.hpp>
#include <mln/morpho/watershed.hpp>
#include <mln/morpho/structural/dilate.hpp>
#include <mln/io/imread.hpp>
#include <mln/io/imsave.hpp>
class c23_t : public mln::dyn_wneighborhood_base<std::array<mln::point2d, 8>,
std::array<int, 8>,
mln::constant_neighborhood_tag, c23_t>
{
public:
static constexpr std::array<mln::point2d, 8> dpoints = {
{{-1, -1}, {-1, +0}, {-1, 1}, //
{+0, -1}, {+0, 1}, //
{+1, -1}, {+1, +0}, {+1, 1}} //
};
static constexpr std::array<int, 8> weights = {
{3, 2, 3,
2, 2,
3, 2, 3}
};
};
constexpr std::array<mln::point2d, 8> c23_t::dpoints;
constexpr std::array<int, 8> c23_t::weights;
static const std::vector<mln::rgb8> regions_table =
{
{0, 0, 0}, {255, 255, 255}, {0, 127, 255}, {127, 255, 0}, {255, 0, 127}, {148, 148, 148},
{0, 255, 128}, {128, 0, 255}, {255, 128, 0}, {1, 0, 159}, {0, 159, 1}, {159, 1, 0},
{96, 254, 255}, {255, 96, 254}, {254, 255, 96}, {21, 125, 126}, {126, 21, 125}, {125, 126, 21},
{116, 116, 255}, {116, 255, 116}, {255, 116, 116}, {0, 227, 228}, {228, 0, 227}, {227, 228, 0},
{28, 27, 255}, {27, 255, 28}, {255, 28, 27}, {59, 59, 59}, {176, 195, 234}, {255, 196, 175},
{68, 194, 171}, {171, 68, 194}, {194, 171, 68}, {71, 184, 72}, {72, 71, 184}, {184, 72, 71},
{188, 255, 169}, {63, 179, 252}, {179, 252, 63}, {252, 63, 179}, {0, 9, 80}, {9, 80, 0},
{80, 0, 9}, {250, 175, 255}, {213, 134, 199}, {95, 100, 115}, {0, 163, 188}, {163, 188, 0},
{188, 0, 163}, {0, 73, 203}, {73, 203, 0}, {203, 0, 73}, {0, 189, 94}, {94, 0, 189},
{189, 94, 0}, {119, 243, 187}, {32, 125, 55}, {55, 32, 125}, {125, 55, 32}, {185, 102, 255},
{255, 185, 102}, {168, 209, 120}, {119, 166, 208}, {192, 96, 135}, {41, 255, 182}, {130, 153, 83},
{55, 88, 247}, {55, 247, 89}, {247, 55, 88}, {0, 75, 87}, {75, 87, 0}, {87, 0, 75},
{59, 135, 200}, {127, 213, 51}, {162, 255, 255}, {182, 37, 255}, {255, 182, 37}, {117, 57, 228},
{210, 163, 142}, {228, 117, 57}, {246, 255, 193}, {123, 107, 188}, {107, 194, 123}, {5, 59, 145},
{59, 145, 5}, {145, 5, 59}, {198, 39, 119}, {23, 197, 40}, {40, 23, 197}, {197, 40, 23},
{158, 199, 178}, {121, 201, 255}, {223, 223, 134}, {84, 253, 39}, {15, 203, 149}, {149, 15, 203},
{203, 149, 15}, {90, 144, 152}, {139, 75, 143}, {132, 97, 71}, {219, 65, 224}, {224, 219, 65},
{40, 255, 255}, {69, 223, 218}, {0, 241, 74}, {74, 0, 241}, {241, 74, 0}, {51, 171, 122},
{227, 211, 220}, {87, 127, 61}, {176, 124, 90}, {13, 39, 36}, {255, 142, 165}, {255, 38, 255},
{255, 255, 38}, {107, 50, 83}, {165, 142, 224}, {9, 181, 255}, {181, 255, 9}, {255, 9, 181},
{70, 238, 140}, {5, 74, 255}, {255, 5, 74}, {51, 84, 138}, {101, 172, 31}, {17, 115, 177},
{0, 0, 221}, {0, 221, 0}, {221, 0, 0}, {200, 255, 220}, {50, 41, 0}, {205, 150, 255},
{116, 45, 178}, {189, 255, 113}, {44, 0, 47}, {171, 119, 40}, {255, 107, 205}, {172, 115, 177},
{236, 73, 133}, {168, 0, 109}, {207, 46, 168}, {203, 181, 188}, {35, 188, 212}, {52, 97, 90},
{184, 209, 39}, {152, 164, 41}, {70, 46, 227}, {227, 70, 46}, {255, 156, 211}, {222, 146, 98},
{95, 56, 136}, {152, 54, 102}, {0, 142, 86}, {86, 0, 142}, {142, 86, 0}, {96, 223, 86},
{46, 135, 246}, {120, 208, 4}, {158, 233, 212}, {214, 92, 177}, {88, 147, 104}, {147, 240, 149},
{148, 93, 227}, {133, 255, 72}, {194, 27, 209}, {255, 255, 147}, {0, 93, 44}, {158, 36, 160},
{0, 233, 182}, {217, 94, 96}, {88, 103, 218}, {38, 154, 163}, {139, 114, 118}, {43, 0, 94},
{174, 164, 113}, {114, 188, 168}, {119, 23, 0}, {93, 86, 42}, {202, 226, 255}, {155, 191, 80},
{136, 158, 255}, {62, 247, 0}, {88, 146, 234}, {229, 183, 0}, {36, 212, 110}, {161, 143, 0},
{210, 191, 105}, {0, 164, 133}, {89, 30, 41}, {132, 0, 164}, {42, 89, 30}, {217, 222, 178},
{11, 22, 121}, {22, 107, 221}, {255, 151, 69}, {3, 158, 45}, {45, 3, 158}, {158, 45, 3},
{29, 42, 86}, {22, 122, 9}, {110, 209, 213}, {57, 221, 53}, {91, 101, 159}, {45, 140, 93},
{37, 213, 247}, {0, 34, 185}, {34, 185, 0}, {185, 0, 34}, {172, 0, 236}, {78, 180, 210},
{221, 107, 231}, {43, 49, 162}, {49, 162, 43}, {162, 43, 49}, {213, 248, 36}, {214, 0, 114},
{248, 36, 213}, {243, 34, 149}, {167, 158, 185}, {224, 122, 144}, {149, 245, 34}, {98, 31, 255},
{255, 98, 31}, {193, 200, 152}, {95, 80, 255}, {63, 123, 128}, {72, 62, 102}, {148, 62, 255},
{108, 226, 151}, {255, 99, 159}, {126, 255, 226}, {136, 223, 98}, {255, 95, 80}, {15, 153, 225},
{211, 41, 73}, {41, 71, 212}, {187, 217, 83}, {79, 235, 180}, {127, 166, 0}, {243, 135, 251},
{0, 41, 229}, {229, 0, 41}, {216, 255, 82}, {249, 174, 141}, {255, 215, 249}, {79, 31, 167},
{167, 79, 31}, {185, 102, 213}, {83, 215, 255}, {40, 2, 4}, {220, 171, 224}, {4, 0, 41},
{90, 50, 6}, {113, 15, 221}, {221, 113, 15}, {115, 0, 33},
};
mln::rgb8 regions_lut(int x)
{
static int sz = static_cast<int>(regions_table.size());
return regions_table[x % sz];
}
mln::rgb8 heat_lut(float x)
{
assert(0 <= x && x <= 1);
float x0 = 1.f / 4.f;
float x1 = 2.f / 4.f;
float x2 = 3.f / 4.f;
if (x < x0)
{
auto g = static_cast<mln::uint8>(x / x0 * 255);
return mln::rgb8{0, g, 255};
}
else if (x < x1)
{
auto b = static_cast<mln::uint8>((x1 - x) / x0 * 255);
return mln::rgb8{0, 255, b};
}
else if (x < x2)
{
auto r = static_cast<mln::uint8>((x - x1) / x0 * 255);
return mln::rgb8{r, 255, 0};
}
else
{
auto b = static_cast<mln::uint8>((1.f - x) / x0 * 255);
return mln::rgb8{255, b, 0};
}
}
int main(int argc, char** argv)
{
if (argc < 4)
{
std::cout << "Usage:" << argv[0] << " input distance output\n";
std::abort();
}
mln::image2d<mln::uint8> input;
mln::io::imread(argv[1], input);
// BEGIN
// (1) Compute the distance transform
auto d = mln::transforms::chamfer_distance_transform<mln::uint8>(input, c23_t());
// Remove non-meaninfull extrema
d = mln::morpho::structural::dilate(d, mln::se::disc(5));
// (2) Inverse the distance
auto dinv = mln::transform(d, [](mln::uint8 x) -> mln::uint8 { return mln::value_traits<mln::uint8>::max() - x; });
// (3) Run the watershed segmentation
int nlabel;
auto ws = mln::morpho::watershed<mln::int16>(dinv, mln::c8, nlabel);
// (4) Labelize input
auto output = mln::where(ws == 0, 1, mln::where(input, ws, 0));
// END
auto d_stretched = mln::data::stretch<float>(d);
mln::io::imsave(mln::transform(d_stretched, heat_lut), argv[2]);
mln::io::imsave(mln::transform(output, regions_lut), argv[3]);
}
......@@ -93,6 +93,8 @@ namespace mln
extended_by_value_image_pixel() = default;
extended_by_value_image_pixel(const extended_by_value_image_pixel&) = default;
extended_by_value_image_pixel(Morpher* ima) : m_ima(ima) {}
extended_by_value_image_pixel(Morpher& ima, const Pix& pix) : m_ima(&ima), m_pix(pix) {}
// Conversion non-const -> const pixel
......
......@@ -100,7 +100,7 @@ namespace mln
ndimage_pixel() = default;
ndimage_pixel(image_type* ima, T* ptr) : details::ndimage_pixel_base<T, dim>(ptr), m_ima(ima) {}
ndimage_pixel(image_type* ima, T* ptr = nullptr) : details::ndimage_pixel_base<T, dim>(ptr), m_ima(ima) {}
/// \brief Copy / copy conversion constructor
template <class U, class J, typename = std::enable_if_t<std::is_convertible<U*, T*>::value>>
......
#ifndef MLN_LABELING_LOCAL_MINIMA_HPP
# define MLN_LABELING_LOCAL_MINIMA_HPP
# include <mln/core/trace.hpp>
# include <mln/core/image/image.hpp>
# include <mln/core/neighborhood/neighborhood.hpp>
# include <mln/core/extension/extension.hpp>
# include <mln/core/extension/fill.hpp>
namespace mln
{
namespace labeling
{
/// Compute and labelize the local mimima of an image
///
/// \param[in] input Input image
/// \param[in] nbh neighborhood
/// \param[out] nlabel Number of minima
/// \param[in] cmp (optional) Strict weak ordering comparison function for values
template <class Label_t, class I, class N, class Compare = std::less<mln_value(I)>>
mln_ch_value(I, Label_t) local_minima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel, Compare cmp = Compare());
/// Compute and labelize the local maxima of an image
///
/// \param[in] input Input image
/// \param[in] nbh neighborhood
/// \param[out] nlabel Number of maxima
template <class Label_t, class I, class N>
mln_ch_value(I, Label_t) local_maxima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel);
/******************************************/
/**** Implementation ****/
/******************************************/
namespace details
{
template <class Parent, class Point>
Point zfindroot(Parent& parent, Point& p)
{
auto& par = parent(p);
if (par != p)
p = zfindroot(parent, par);
return p;
}
// Output must initialized to Label_t(-1)
template <class I, class O, class N, class Compare>
int local_minima(const I& input, const N& nbh, O& output,
Compare cmp)
{
using Label_t = mln_value(O);
constexpr Label_t kUnseen = static_cast<Label_t>(-1);
mln_ch_value(I, mln_point(I)) parent = imchvalue<mln_point(I)>(input);
mln_pixter(pxIn, input);
mln_pixter(pxOut, output);
mln_pixter(pxPar, parent);
mln_iter(nxIn, nbh(pxIn));
mln_iter(nxOut, nbh(pxOut));
mln_iter(nxPar, nbh(pxPar));
{
// Forward scan
mln_forall(pxIn, pxOut, pxPar)
{
auto p = pxPar->point();
auto p_input_value = pxIn->val();
// Make set
auto root = p;
bool is_possible_minimum = true;
mln_forall(nxIn, nxOut, nxPar)
{
if (nxOut->val() == kUnseen)
continue;
auto n_input_value = nxIn->val();
if (cmp(n_input_value, p_input_value)) // input(n) < input(p)
{
is_possible_minimum = false;
}
else if (cmp(p_input_value, n_input_value)) // input(p) < input(n) -> n is not a local minimum
{
const auto r = zfindroot(parent, nxPar->val());
output(r) = false;
}
else // Flat zone linking
{
const auto n_root = zfindroot(parent, nxPar->val());
if (n_root != root)
{
is_possible_minimum &= output(n_root);
auto roots = std::minmax(root, n_root);
parent(roots.second) = roots.first;
root = roots.first;
}
}
}
pxPar->val() = root;
pxOut->val() = is_possible_minimum;
output(root) = is_possible_minimum;
}
}
int n_label = 0;
// Labelisation
{
mln_forall(pxOut, pxPar)
{
if (pxOut->val())
{
if (pxPar->point() == pxPar->val()) // New minimum
pxOut->val() = ++n_label;
else
{
auto r = zfindroot(parent, pxPar->val());
pxOut->val() = output(r);
}
}
}
}
mln_assertion((n_label <= value_traits<Label_t>::max()) &&
"Overflow detected (the number of label has exceed the label type capacity");
return n_label;
}
}
template <class Label_t, class I, class N, class Compare>
mln_ch_value(I, Label_t)
local_minima(const Image<I>& input_, const Neighborhood<N>& nbh_, int& nlabel, Compare cmp)
{
mln_entering("mln::labeling::local_minima");
static_assert(std::is_integral<Label_t>::value && !std::is_same<Label_t, bool>::value,
"Label type must a non-bool integral type");
constexpr Label_t kUnseen = static_cast<Label_t>(-1);
const I& input = exact(input_);
const N& nbh = exact(nbh_);
mln_ch_value(I, Label_t) output = imchvalue<Label_t>(input).adjust(nbh).init(kUnseen);
if (extension::need_adjust(output, nbh))
{
auto out = extension::add_value_extension(output, kUnseen);
nlabel = details::local_minima(input, nbh, out, cmp);
}
else
{
mln::extension::fill(output, kUnseen);
nlabel = details::local_minima(input, nbh, output, cmp);
}
return output;
}
template <class Label_t, class I, class N>
mln_ch_value(I, Label_t)
local_maxima(const Image<I>& input, const Neighborhood<N>& nbh, int& nlabel)
{
return local_minima<Label_t>(input, nbh, nlabel, std::greater<mln_value(I)> ());
}
} // end of namespace mln::labeling
} // end of namespace mln
#endif //!MLN_LABELING_LOCAL_MINIMA_HPP
#ifndef MLN_MORPHO_PRIVATE_PQUEUE_HPP
# define MLN_MORPHO_PRIVATE_PQUEUE_HPP
# include <mln/morpho/private/pqueue_hqueue_fifo.hpp>
namespace mln
{
namespace morpho
{
namespace details
{
/// Provides a priority queue of points with a the fifo property
template <class I>
class pqueue_fifo
{
using key_type = mln_value(I);
using value_type = mln_point(I);
public:
template <class J>
pqueue_fifo(const Image<J>& f);
void push(const key_type& priority, const value_type& element);
void pop();
std::pair<key_type, value_type> top() const;
bool empty() const;
private:
static_assert(value_traits<key_type>::quant <= 16,
"Only low quantized type supported.");
pqueue_hqueue_fifo<mln_concrete(I)> m_delegate;
};
/******************************************/
/**** Implementation ****/
/******************************************/
template <class I>
template <class J>
pqueue_fifo<I>::pqueue_fifo(const Image<J>& f)
: m_delegate(f)
{
}
template <class I>
void pqueue_fifo<I>::push(const key_type& k, const value_type& v)
{
m_delegate.push(k, v);
}
template <class I>
void pqueue_fifo<I>::pop()
{
m_delegate.pop();
}
template <class I>
bool pqueue_fifo<I>::empty() const
{
return m_delegate.empty();