Commit 243fdbc6 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Migrate morphological gradients

parent ae1dc3b3
...@@ -14,6 +14,7 @@ Structural Morphological Operation ...@@ -14,6 +14,7 @@ Structural Morphological Operation
morpho/hit_or_miss morpho/hit_or_miss
morpho/rank_filter morpho/rank_filter
morpho/median_filter morpho/median_filter
morpho/gradient
morpho/opening_by_reconstruction morpho/opening_by_reconstruction
......
Morphological Gradients
=======================
Include :file:`<mln/morpho/gradient.hpp>`
.. cpp:namespace:: mln::morpho
.. cpp:function:: \
Image{I} concrete_t<I> gradient(I f, StructuringElement se)
Image{I} concrete_t<I> external_gradient(I f, StructuringElement se)
Image{I} concrete_t<I> internal_gradient(I f, StructuringElement se)
Compute the morphological gradients. Morphological gradients are
operators enhancing variations of pixel intensity in a neighborhood
determined by a structuring element. Three combinations are currently used:
#. arithmetic difference between the dilation and the erosion; this is basic
morphological gradient, also called **Beucher** gradient.
:math:`\rho_B = \delta_B - \varepsilon_B`
#. arithmetic difference between the dilation and the original image; also called **external gradient**:
:math:`\rho_B^{+} = \delta_B - \mathrm{id}`
#. arithmetic difference between the original image and its erosion; also called the **internal_gradient**:
:math:`\rho_B^{-} = \mathrm{id} - \varepsilon_B`
If input is not integral, the marginal ordering is considered and the 𝑙₂ of the vector
difference is returned:
:math:`\rho_B = \left| \delta_B - \varepsilon_B \right|`
:param f: Input image 𝑓
:param se: Elementary structuring element.
:return: An image whose type is deduced from the input image.
:precondition:
:exception: N/A
Example 1 : Gradient by a square on a gray-level image
------------------------------------------------------
.. code-block:: cpp
#include <mln/morpho/gradient.hpp>
#include <mln/core/se/rect2d.hpp>
// Define a square SE of size 7x7
auto input = ...;
auto rect = mln::se::rect2d(7,7);
auto grad1 = mln::morpho::gradient(input, rect);
auto grad2 = mln::morpho::internal_gradient(input, rect);
auto grad3 = mln::morpho::external_gradient(input, rect);
.. list-table::
* - .. figure:: /images/lena_gray.jpg
- .. figure:: /images/morpho_gradient_1.png
Beucher (thick) gradient
* - .. figure:: /images/morpho_int_gradient_1.png
Internal gradient
- .. figure:: /images/morpho_ext_gradient_1.png
External Gradient
\ No newline at end of file
...@@ -27,6 +27,10 @@ add_image("erosion-cli;dilation;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho ...@@ -27,6 +27,10 @@ add_image("erosion-cli;dilation;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho
add_image("erosion-cli;opening;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_opening_1.png) add_image("erosion-cli;opening;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_opening_1.png)
add_image("erosion-cli;closing;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_closing_1.png) add_image("erosion-cli;closing;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_closing_1.png)
add_image("erosion-cli;median;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_median_1.png) add_image("erosion-cli;median;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_median_1.png)
add_image("erosion-cli;gradient;square;7" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_gradient_1.png)
add_image("erosion-cli;ext_gradient;square;7" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_ext_gradient_1.png)
add_image("erosion-cli;int_gradient;square;7" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_int_gradient_1.png)
add_image("cdt;1" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-4.png) add_image("cdt;1" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-4.png)
add_image("cdt;2" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-8.png) add_image("cdt;2" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-8.png)
add_image("cdt;3" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-2-3.png) add_image("cdt;3" "${DOCUMENTATION_IMAGE_DIR}/F.png" F-2-3.png)
......
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
#include <mln/morpho/experimental/erosion.hpp> #include <mln/morpho/experimental/erosion.hpp>
#include <mln/morpho/experimental/opening.hpp> #include <mln/morpho/experimental/opening.hpp>
#include <mln/morpho/experimental/median_filter.hpp> #include <mln/morpho/experimental/median_filter.hpp>
#include <mln/morpho/experimental/gradient.hpp>
#include <mln/io/experimental/imread.hpp> #include <mln/io/experimental/imread.hpp>
#include <mln/io/experimental/imsave.hpp> #include <mln/io/experimental/imsave.hpp>
...@@ -56,6 +57,9 @@ enum morpho_op_type ...@@ -56,6 +57,9 @@ enum morpho_op_type
kOpening, kOpening,
kClosing, kClosing,
kMedian, kMedian,
kGradient,
kInternalGradient,
kExternalGradient
}; };
std::istream& operator>>(std::istream& in, morpho_op_type& se) std::istream& operator>>(std::istream& in, morpho_op_type& se)
...@@ -73,6 +77,12 @@ std::istream& operator>>(std::istream& in, morpho_op_type& se) ...@@ -73,6 +77,12 @@ std::istream& operator>>(std::istream& in, morpho_op_type& se)
se = kClosing; se = kClosing;
else if (token == "median") else if (token == "median")
se = kMedian; se = kMedian;
else if (token == "gradient")
se = kGradient;
else if (token == "ext_gradient")
se = kExternalGradient;
else if (token == "int_gradient")
se = kInternalGradient;
else else
throw po::invalid_option_value("Invalid Operator"); throw po::invalid_option_value("Invalid Operator");
return in; return in;
...@@ -115,7 +125,7 @@ int main(int argc, char** argv) ...@@ -115,7 +125,7 @@ int main(int argc, char** argv)
{ {
std::cout << "Usage: " << argv[0] std::cout << "Usage: " << argv[0]
<< " [-h,--help] OPERATOR SE size input output\n" << " [-h,--help] OPERATOR SE size input output\n"
"OPERATOR\t Morphological operation to perform [erosion | dilation]\n" "OPERATOR\t Morphological operation to perform [erosion | dilation | opening | closing | median | gradient | ext_gradient | int_gradient]\n"
"SE\t Structuring element to use [square | disc | diamond]\n" "SE\t Structuring element to use [square | disc | diamond]\n"
"size\t Size of the SE\n" "size\t Size of the SE\n"
"input\t Input image (u8)\n" "input\t Input image (u8)\n"
...@@ -146,6 +156,15 @@ int main(int argc, char** argv) ...@@ -146,6 +156,15 @@ int main(int argc, char** argv)
case kMedian: case kMedian:
output = mln::morpho::experimental::median_filter(input, nbh, mln::extension::bm::mirror{}); output = mln::morpho::experimental::median_filter(input, nbh, mln::extension::bm::mirror{});
break; break;
case kGradient:
output = mln::morpho::experimental::gradient(input, nbh);
break;
case kExternalGradient:
output = mln::morpho::experimental::external_gradient(input, nbh);
break;
case kInternalGradient:
output = mln::morpho::experimental::internal_gradient(input, nbh);
break;
} }
mln::io::experimental::imsave(output, vm["output"].as<std::string>()); mln::io::experimental::imsave(output, vm["output"].as<std::string>());
......
...@@ -32,10 +32,14 @@ namespace mln ...@@ -32,10 +32,14 @@ namespace mln
template <class InputImage, class OutputImage, class UnaryFunction> template <class InputImage, class OutputImage, class UnaryFunction>
void transform(InputImage in, const Image<OutputImage>& out, UnaryFunction f); void transform(InputImage in, const Image<OutputImage>& out, UnaryFunction f);
template <class InputImage1, class InputImage2, class OutputImage, class BinaryFunction>
void transform(InputImage1 in1, InputImage2 in2, OutputImage out, BinaryFunction f);
template <class InputImage, class UnaryFunction> template <class InputImage, class UnaryFunction>
image_ch_value_t<InputImage, std::decay_t<std::invoke_result_t<UnaryFunction, image_reference_t<InputImage>>>> image_ch_value_t<InputImage, std::decay_t<std::invoke_result_t<UnaryFunction, image_reference_t<InputImage>>>>
transform(InputImage in, UnaryFunction f); transform(InputImage in, UnaryFunction f);
/******************************************/ /******************************************/
...@@ -62,9 +66,36 @@ namespace mln ...@@ -62,9 +66,36 @@ namespace mln
::ranges::transform(r1, ::ranges::begin(r2), f); ::ranges::transform(r1, ::ranges::begin(r2), f);
} }
template <class InputImage1, class InputImage2, class OutputImage, class BinaryFunction>
void transform(InputImage1 in1, InputImage2 in2, OutputImage out, BinaryFunction f)
{
static_assert(mln::is_a<InputImage1, experimental::Image>());
static_assert(mln::is_a<InputImage2, experimental::Image>());
static_assert(mln::is_a<OutputImage, experimental::Image>());
static_assert(
::ranges::Invocable<BinaryFunction, image_reference_t<InputImage1>, image_reference_t<InputImage2>>());
static_assert(
std::is_convertible_v<
std::invoke_result_t<BinaryFunction, image_reference_t<InputImage1>, image_reference_t<InputImage2>>,
image_value_t<OutputImage>>,
"The result of the function is not implicitely convertible to the output image value type");
mln_entering("mln::transform");
// FIXME: disabled becaused some domain (e.g. morphed) do not implement comparison for now
// mln_precondition(in2.domain() == out.domain());
// mln_precondition(in1.domain() == out.domain());
auto&& ivals1 = in1.new_values();
auto&& ivals2 = in2.new_values();
auto&& ovals = out.new_values();
for (auto [r1, r2, r3] : ranges::view::zip(ranges::rows(ivals1), ranges::rows(ivals2), ranges::rows(ovals)))
::ranges::transform(r1, ::ranges::begin(r2), ::ranges::begin(r3), f);
}
template <class InputImage, class UnaryFunction> template <class InputImage, class UnaryFunction>
image_ch_value_t<InputImage, std::decay_t<std::invoke_result_t<UnaryFunction, image_reference_t<InputImage>>>> image_ch_value_t<InputImage, std::decay_t<std::invoke_result_t<UnaryFunction, image_reference_t<InputImage>>>>
transform(InputImage in, UnaryFunction f) transform(InputImage in, UnaryFunction f)
{ {
static_assert(mln::is_a<InputImage, experimental::Image>()); static_assert(mln::is_a<InputImage, experimental::Image>());
static_assert(::ranges::Invocable<UnaryFunction, image_reference_t<InputImage>>()); static_assert(::ranges::Invocable<UnaryFunction, image_reference_t<InputImage>>());
......
#pragma once
#include <mln/core/algorithm/transform.hpp>
#include <mln/core/image/image.hpp>
#include <mln/core/trace.hpp>
#include <mln/core/concept/new/structuring_elements.hpp>
#include <mln/morpho/experimental/dilation.hpp>
#include <mln/morpho/experimental/erosion.hpp>
#include <mln/core/private/maths_ops.hpp>
/// \file
namespace mln::morpho::experimental
{
/// \ingroup morpho
/// \brief Compute the Beucher gradient.
///
/// The beucher gradient is defined as:
/// \f[
/// |\nabla u(p)| = \left| \bigvee_{q \in \mathcal{N}(p)} f(q) -
/// \bigwedge_{q \in \mathcal{N}(p)} f(q) \right|
/// \f]
///
/// + If the optional \p out image is provided, it must be wide enough to store
/// the results (the function does not perform any resizing).
///
/// + If a optional \p cmp function is provided, the algorithm
/// will internally do an unqualified call to `sup(x, y, cmp)`
/// and `inf(x, y,cmp)`. The default is the product-order so
/// that it works for vectorial type as well.
///
/// \param[in] input Input image.
/// \param[in] se Neighborhood/SE/Window to look around.
/// \param[in] cmp (optional) Comparaison function. The method internally does an
/// unqualified call to `inf(x, y, cmp)` and `sup(x, y, cmp)`. Default
/// is the product-order.
/// \param[in] norm (optional) Norm function used in \f$|x - y|\f$. The default is
/// the ℓ₂-norm.
/// \param[out] out (optional) Output image to write in.
///
/// FIXME: This can be done in a single pass using kernels but you
/// we have to consider:
/// * an incrementable accumlator
/// * an extension that is mirrorizable
namespace details
{
template <class T, class enable = void>
struct grad_op
{
using result_type = decltype(mln::l2norm(std::declval<T>()));
result_type operator()(T a, T b) const noexcept { return mln::l2norm(T(a - b)); }
};
template <class T>
struct grad_op<T, std::enable_if_t<std::is_arithmetic_v<T>>>
{
using result_type = T;
T operator()(T a, T b) const noexcept { return a - b; }
};
template <class I>
using gradient_result_t = image_ch_value_t<I, typename grad_op<image_value_t<I>>::result_type>;
}
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> gradient(InputImage&& input, const SE& se);
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> external_gradient(InputImage&& input, const SE& se);
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> internal_gradient(InputImage&& input, const SE& se);
/*************************/
/*** Implementation ***/
/*************************/
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> gradient(InputImage&& ima, const SE& se)
{
using I = std::remove_reference_t<InputImage>;
static_assert(mln::is_a<I, mln::experimental::Image>());
static_assert(mln::is_a<SE, mln::experimental::StructuringElement>());
mln_entering("mln::morpho::gradient");
auto d = morpho::experimental::dilation(ima, se);
auto e = morpho::experimental::erosion(ima, se);
using R = typename details::grad_op<image_value_t<I>>::result_type;
details::gradient_result_t<I> out = imchvalue<R>(ima);
mln::transform(d, e, out, details::grad_op<image_value_t<I>>());
return out;
}
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> external_gradient(InputImage&& ima, const SE& se)
{
using I = std::remove_reference_t<InputImage>;
static_assert(mln::is_a<I, mln::experimental::Image>());
static_assert(mln::is_a<SE, mln::experimental::StructuringElement>());
mln_entering("mln::morpho::external_gradient");
auto d = morpho::experimental::dilation(ima, se);
using R = typename details::grad_op<image_value_t<I>>::result_type;
details::gradient_result_t<I> out = imchvalue<R>(ima);
mln::transform(d, ima, out, details::grad_op<image_value_t<I>>());
return out;
}
template <class InputImage, class SE>
details::gradient_result_t<std::remove_reference_t<InputImage>> internal_gradient(InputImage&& ima, const SE& se)
{
using I = std::remove_reference_t<InputImage>;
static_assert(mln::is_a<I, mln::experimental::Image>());
static_assert(mln::is_a<SE, mln::experimental::StructuringElement>());
mln_entering("mln::morpho::internal_gradient");
auto e = morpho::experimental::erosion(ima, se);
using R = typename details::grad_op<image_value_t<I>>::result_type;
details::gradient_result_t<I> out = imchvalue<R>(ima);
mln::transform(ima, e, out, details::grad_op<image_value_t<I>>());
return out;
}
} // namespace mln::morpho::experimental
#include <mln/core/algorithm/all_of.hpp> #include <mln/morpho/experimental/gradient.hpp>
#include <mln/core/algorithm/iota.hpp> #include <mln/core/algorithm/iota.hpp>
#include <mln/core/grays.hpp> #include <mln/core/algorithm/paste.hpp>
#include <mln/core/image/image2d.hpp> #include <mln/core/algorithm/all_of.hpp>
#include <mln/core/image/experimental/ndimage.hpp>
#include <mln/core/image/view/operators.hpp> #include <mln/core/image/view/operators.hpp>
#include <mln/core/neighb2d.hpp> #include <mln/core/image/view/filter.hpp>
#include <mln/core/se/rect2d.hpp> #include <mln/core/se/rect2d.hpp>
#include <mln/io/imread.hpp> #include <mln/core/colors.hpp>
#include <mln/io/imsave.hpp>
#include <mln/morpho/structural/gradient.hpp> #include <functional>
#include <fixtures/ImagePath/image_path.hpp> #include <fixtures/ImagePath/image_path.hpp>
#include <fixtures/ImageCompare/image_compare.hpp>
#include <mln/io/experimental/imread.hpp>
#include <mln/io/experimental/imprint.hpp>
#include <gtest/gtest.h> #include <gtest/gtest.h>
using namespace mln; using namespace mln;
TEST(Morpho, gradient_gradient_0) TEST(Morpho, gradient_binary)
{ {
using namespace mln::view::ops; mln::experimental::image2d<bool> ima = {{0, 0, 0, 0, 1, 0, 0}, //
{0, 0, 0, 0, 1, 0, 0}, //
{0, 1, 1, 1, 1, 1, 0}, //
{0, 1, 1, 1, 1, 1, 1}, //
{0, 1, 1, 1, 1, 1, 0}}; //
image2d<uint8> ima(10, 10);
iota(ima, 10);
{ // Fast: border wide enough mln::experimental::se::rect2d win(3, 3);
mln::se::rect2d win(3, 1);
auto out = morpho::structural::gradient(ima, win);
static_assert(std::is_same<decltype(out)::value_type, int>::value, "Error integral promotion should give int.");
ASSERT_TRUE(all_of(out == 1 || out == 2)); {
mln::experimental::image2d<bool> ref = {{0, 0, 0, 1, 1, 1, 0}, //
{1, 1, 1, 1, 1, 1, 1}, //
{1, 1, 1, 1, 1, 1, 1}, //
{1, 1, 0, 0, 0, 1, 1}, //
{1, 1, 0, 0, 0, 1, 1}}; //
auto out = mln::morpho::experimental::gradient(ima, win);
ASSERT_IMAGES_EQ_EXP(out, ref);
} }
}
TEST(Morpho, gradient_gradient_1) {
{ mln::experimental::image2d<bool> ref = {{0, 0, 0, 1, 0, 1, 0}, //
image2d<uint8> ima; {1, 1, 1, 1, 0, 1, 1}, //
io::imread(fixtures::ImagePath::concat_with_filename("small.pgm"), ima); {1, 0, 0, 0, 0, 0, 1}, //
{1, 0, 0, 0, 0, 0, 0}, //
{1, 0, 0, 0, 0, 0, 1}}; //
{ // Fast: border wide enough
mln::se::rect2d win(7, 7); auto out = mln::morpho::experimental::external_gradient(ima, win);
auto out = morpho::structural::gradient(ima, win); ASSERT_IMAGES_EQ_EXP(out, ref);
} }
{
mln::experimental::image2d<bool> ref = {{0, 0, 0, 0, 1, 0, 0}, //
{0, 0, 0, 0, 1, 0, 0}, //
{0, 1, 1, 1, 1, 1, 0}, //
{0, 1, 0, 0, 0, 1, 1}, //
{0, 1, 0, 0, 0, 1, 0}}; //
auto out = mln::morpho::experimental::internal_gradient(ima, win);
ASSERT_IMAGES_EQ_EXP(out, ref);
}
} }
// Border is not wide enough => call dilate + erode TEST(Morpho, gradient_grayscale)
TEST(Morpho, gradient_gradient_2)
{ {
using namespace mln::view::ops; mln::experimental::image2d<uint8_t> ima = {{0, 0, 0, 0, 9, 0, 0}, //
{0, 0, 0, 0, 9, 0, 0}, //
{0, 9, 9, 9, 9, 9, 0}, //
{0, 9, 9, 9, 9, 9, 9}, //
{0, 9, 9, 9, 9, 9, 0}}; //
image2d<uint8> ima(0);
image2d<uint8> ima2;
io::imread(fixtures::ImagePath::concat_with_filename("small.pgm"), ima);
io::imread(fixtures::ImagePath::concat_with_filename("small.pgm"), ima2);
mln::se::rect2d win(3, 3); mln::experimental::se::rect2d win(3, 3);
auto out1 = morpho::structural::gradient(ima, win);
auto out2 = morpho::structural::gradient(ima2, win);
ASSERT_TRUE(all_of(out1 == out2)); {
mln::experimental::image2d<uint8_t> ref = {{0, 0, 0, 9, 9, 9, 0}, //
{9, 9, 9, 9, 9, 9, 9}, //
{9, 9, 9, 9, 9, 9, 9}, //
{9, 9, 0, 0, 0, 9, 9}, //
{9, 9, 0, 0, 0, 9, 9}}; //
auto out = mln::morpho::experimental::gradient(ima, win);
ASSERT_IMAGES_EQ_EXP(out, ref);
}
{
mln::experimental::image2d<uint8_t> ref = {{0, 0, 0, 9, 0, 9, 0}, //
{9, 9, 9, 9, 0, 9, 9}, //
{9, 0, 0, 0, 0, 0, 9}, //
{9, 0, 0, 0, 0, 0, 0}, //
{9, 0, 0, 0, 0, 0, 9}}; //
auto out = mln::morpho::experimental::external_gradient(ima, win);
ASSERT_IMAGES_EQ_EXP(out, ref);
}
{
mln::experimental::image2d<uint8_t> ref = {{0, 0, 0, 0, 9, 0, 0}, //
{0, 0, 0, 0, 9, 0, 0}, //
{0, 9, 9, 9, 9, 9, 0}, //
{0, 9, 0, 0, 0, 9, 9}, //
{0, 9, 0, 0, 0, 9, 0}}; //
auto out = mln::morpho::experimental::internal_gradient(ima, win);
ASSERT_IMAGES_EQ_EXP(out, ref);
}
} }
// Dilation on a with a vmorph / binary case
TEST(Morpho, gradient_gradient_3) // This test is disabled because of a lifetime problem on the filter view (the domain lifetime depends on the image which might be a temporary).
TEST(Morpho, DISABLED_gradient_on_view)
{ {
image2d<uint8> ima; using namespace mln::view::ops;
io::imread(fixtures::ImagePath::concat_with_filename("small.pgm"), ima); mln::experimental::image2d<uint8> ima(10, 5);
mln::iota(ima, 0);
// Morpher has no extension // Morpher has no extension
mln::se::rect2d win(3, 3); mln::experimental::se::rect2d win(3, 3);
auto out = morpho::structural::gradient(ima > 128, win);
}
// Dilation on a with a vmorph / binary case constexpr uint8 _X = 0;
TEST(Morpho, gradient_gradient_4) mln::experimental::image2d<uint8_t> ref = {{_X, _X, _X, _X, _X, _X, _X, _X, _X, _X},
{ {_X, _X, _X, _X, _X, _X, _X, _X, _X, _X},
image2d<uint8> ima; {_X, _X, _X, _X, _X, _X, 11, 12, 12, 11},
io::imread(fixtures::ImagePath::concat_with_filename("small.pgm"), ima); {11, 12, 12, 12, 12, 20, 21, 22, 22, 21},
{11, 12, 12, 12, 12, 12, 12, 12, 12, 11}};
mln::se::rect2d win(3, 3);
image2d<uint8> out;