--- milena/mln/inpainting/derivative_kernels.hh | 34 +++-- milena/mln/inpainting/fast_isotropic_diffusion.hh | 2 - milena/mln/inpainting/lpde.hh | 64 +++++++--- milena/mln/inpainting/mean_diffusion.hh | 133 +++++++++++++++++++++ milena/mln/inpainting/tv.hh | 93 ++++++++++++++ 5 files changed, 292 insertions(+), 34 deletions(-) create mode 100644 milena/mln/inpainting/mean_diffusion.hh create mode 100644 milena/mln/inpainting/tv.hh diff --git a/milena/mln/inpainting/derivative_kernels.hh b/milena/mln/inpainting/derivative_kernels.hh index 5b642ee..a7f3c34 100644 --- a/milena/mln/inpainting/derivative_kernels.hh +++ b/milena/mln/inpainting/derivative_kernels.hh @@ -2,23 +2,29 @@ namespace mln { namespace inpainting { - template <typename T> - const T* make_dx_kernel() + namespace { - static const T ws[] = { 0.1, 0.2, 0.1, - 0, 0, 0, - -0.1, -0.2, -0.1 }; - return ws; - } + const float a = 61; + const float b = 17; - template <typename T> - const T* make_dy_kernel() - { - static const T ws[] = { 0.1, 0, -0.1, - 0.2, 0, -0.2, - 0.1, 0, -0.1 }; + template <typename T> + const T* make_dx_kernel() + { + static const T ws[] = { b, a, b, + 0, 0, 0, + -b, -a, -b }; + return ws; + } + + template <typename T> + const T* make_dy_kernel() + { + static const T ws[] = { b, 0, -b, + a, 0, -a, + b, 0, -b }; - return ws; + return ws; + } } template <typename T> diff --git a/milena/mln/inpainting/fast_isotropic_diffusion.hh b/milena/mln/inpainting/fast_isotropic_diffusion.hh index d63f9d4..44a8989 100644 --- a/milena/mln/inpainting/fast_isotropic_diffusion.hh +++ b/milena/mln/inpainting/fast_isotropic_diffusion.hh @@ -76,8 +76,6 @@ namespace mln std::swap(u_, u_next_); } - - std::cout << iterations << "\t"; } } /* mln::inpainting */ diff --git a/milena/mln/inpainting/lpde.hh b/milena/mln/inpainting/lpde.hh index 5f34d83..f38b7a3 100644 --- a/milena/mln/inpainting/lpde.hh +++ b/milena/mln/inpainting/lpde.hh @@ -24,7 +24,10 @@ #include <mln/core/alias/vec2d.hh> #include <mln/inpainting/derivative_kernels.hh> + #include <mln/inpainting/fast_isotropic_diffusion.hh> +#include <mln/inpainting/mean_diffusion.hh> +#include <mln/inpainting/tv.hh> #include <mln/core/concept/accumulator.hh> @@ -42,7 +45,7 @@ namespace mln struct lpde { public: - typedef double float_type; + typedef float float_type; lpde(float dt = 1.0f, float dx = 1.0f, float t_f = 20.0); @@ -96,19 +99,30 @@ namespace mln return components_product(mln_trait_value_nature(T)(), a, b); } + template <typename T> + T normalize(T v) + { + lpde::float_type magnitude = norm::l2(v); + + magnitude = std::max((lpde::float_type)1.0, magnitude); + + return v / magnitude; + } + template <typename I, typename M> lpde::float_type diffuse(const I& u, I& u_next, - //const I& f, const std::vector<lpde::float_type>& w, const M& mask, - lpde::float_type& max_diff_norm) + lpde::float_type& max_diff_norm, + I& u_x, + I& u_y) { typedef lpde::float_type float_type; typedef mln_value(I) T; const float_type activation_threshold = 0; - + const unsigned derivative_kernels_size = 5; const float_type* derivative_kernel[derivative_kernels_size]; @@ -127,13 +141,20 @@ namespace mln mln_pixter(const I) p(u); mln_pixter(I) p_next(u_next); - //mln_pixter(const I) p_f(f); mln_pixter(const M) p_m(mask); + + mln_pixter(I) pu_x(u_x); + mln_pixter(I) pu_y(u_y); + mln_qixter(const I, const window2d) q(p, full3x3); float_type sse = 0; - - for_all_3 (p, p_next, p_m) + + for (p.start(), p_next.start(), p_m.start(), + pu_x.start(), pu_y.start(); + p.is_valid(); + p.next(), p_next.next(), p_m.next(), + pu_x.next(), pu_y.next()) { if (p_m.val()) { @@ -142,9 +163,6 @@ namespace mln if (fabs(w[1]) > activation_threshold) p_next.val() = p.val() * w[1]; - // if (fabs(w[0]) > activation_threshold) - // p_next.val() += w[0] * p_f.val(); - T p_derivative[derivative_kernels_size]; for (unsigned i = 0; i < derivative_kernels_size; ++i) @@ -172,6 +190,9 @@ namespace mln const T p_yy2 = components_product(p_yy, p_yy); const T p_xy2 = components_product(p_xy, p_xy); + pu_x.val() = normalize(p_x); + pu_y.val() = normalize(p_y); + //||∇u||^2 if (fabs(w[2]) > activation_threshold) p_next.val() += w[2] * (p_x2 + p_y2); @@ -200,6 +221,8 @@ namespace mln } } + inpainting::tv(u_next, u_x, u_y, mask, max_diff_norm); + return sse; } @@ -207,10 +230,11 @@ namespace mln void lpde::operator()(I<T>& src, const I<bool>& mask) { - //fast_isotropic_diffusion fast_iso; - //fast_iso(src, mask); - - //const I<T> f(src); + // mean_diffusion mean_diff; + // mean_diff(src, mask); + + // fast_isotropic_diffusion fast_iso; + // fast_iso(src, mask); I<T>& u_ = src; I<T> u_next_(src); @@ -229,17 +253,21 @@ namespace mln float_type min_max_diff_norm = std::numeric_limits<float_type>::max(); + I<T> u_x(u_.domain()); + I<T> u_y(u_.domain()); + float t = 0; while (t < this->t_f && max_diff_norm > CONVERGENCE_THRESHOLD - //&& sse < sse_old - && !std::isnan(sse)) + && !std::isnan(max_diff_norm) + ) { + max_diff_norm = 0; sse_old = sse; - sse = diffuse(*u, *u_next, /* f, */ w, mask, max_diff_norm); + sse = diffuse(*u, *u_next, w, mask, max_diff_norm, u_x, u_y); if (max_diff_norm < min_max_diff_norm) min_max_diff_norm = max_diff_norm; @@ -258,4 +286,4 @@ namespace mln } } /* mln::inpainting */ -} /* mln*/ +} /* mln */ diff --git a/milena/mln/inpainting/mean_diffusion.hh b/milena/mln/inpainting/mean_diffusion.hh new file mode 100644 index 0000000..e0e3f48 --- /dev/null +++ b/milena/mln/inpainting/mean_diffusion.hh @@ -0,0 +1,133 @@ +#include <mln/core/macros.hh> +#include <mln/core/concept/image.hh> +#include <mln/core/alias/point2d.hh> +#include <mln/pw/value.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/core/image/dmorph/image_if.hh> + +#include <mln/core/concept/function.hh> + +#include <mln/data/transform.hh> + +#include <mln/literal/colors.hh> +#include <mln/literal/zero.hh> +#include <mln/pw/image.hh> +#include <mln/convert/to_p_array.hh> +#include <mln/core/var.hh> + +namespace mln +{ + + namespace inpainting + { + + struct mean_diffusion + { + public: + typedef float float_type; + + template <template <typename> class I, typename T> + void operator()(I<T>& src, + const I<bool>& mask_); + }; + + namespace + { + template <typename I, typename N> + struct contour : public Function_v2b< contour<I,N> > + { + typedef bool result; + + contour(const I& ima, const N& neighborhood) + : ima_ (ima), + neighborhood_ (neighborhood) + { + } + + bool operator()(const mln_psite(I)& p) const + { + if (ima_(p)) + { + mln_niter(N) n(neighborhood_, p); + + for_all(n) + { + if (!ima_(n)) + return true; + } + } + + return false; + } + + private: + const I& ima_; + const N& neighborhood_; + }; + + template <template <typename> class I, typename T> + void diffuse(I<T>& src, + const p_array<mln_psite(I<bool>)>& ctr, + I<bool>& mask, + const neighb2d& neighborhood) + { + + mln_piter(p_array<mln_psite(I<bool>)>) p(ctr); + mln_niter_(neighb2d) n(neighborhood, p); + + for_all(p) + { + T accu = literal::zero; + + unsigned nb = 0; + + for_all(n) + { + if (!mask(n)) + { + accu += src(n); + ++nb; + } + } + + src(p) = accu / nb; + } + } + + } /* !mln::inpainting::{anonymous} */ + + template <template <typename> class I, typename T> + void mean_diffusion::operator()(I<T>& src, + const I<bool>& mask_) + { + const neighb2d neighborhood = c8(); + + I<bool> mask(mask_.domain()); + data::paste(mask_, mask); + + unsigned contour_nsites = std::numeric_limits<unsigned>::max(); + + while (contour_nsites > 0) + { + contour<I<bool>, neighb2d> ctr = contour<I<bool>, neighb2d>(mask, neighborhood); + p_array<mln_psite(I<bool>)> ctr_set = convert::to_p_array(mask | ctr); + contour_nsites = ctr_set.nsites(); + + if (contour_nsites > 0) + { + diffuse(src, ctr_set, mask, neighborhood); + + mln_piter(p_array<mln_psite(I<bool>)>) p(ctr_set); + for_all (p) + { + mask(p) = false; + } + } + } + } + + } /* mln::inpainting */ + +} /* mln */ + diff --git a/milena/mln/inpainting/tv.hh b/milena/mln/inpainting/tv.hh new file mode 100644 index 0000000..d5bed23 --- /dev/null +++ b/milena/mln/inpainting/tv.hh @@ -0,0 +1,93 @@ +#include <mln/core/alias/window2d.hh> +#include <mln/core/concept/iterator.hh> +#include <mln/algebra/vec.hh> +#include <mln/convert/from_to.hh> +#include <mln/norm/l2.hh> + +namespace mln +{ + namespace inpainting + { + struct tv + { + typedef float float_type; + }; + + template <typename T> + inline + void interval(T& v, + tv::float_type min, + tv::float_type max) + { + for (unsigned i = 0; i < mln_dim(T); ++i) + v[i] = std::max(std::min(v[i], max), min); + } + + template <typename I, typename M> + void tv(I& u, + const I& u_x, + const I& u_y, + const M& mask, + tv::float_type& max_diff_norm) + { + typedef tv::float_type float_type; + + extension::adjust_duplicate(u_x, 1); + extension::adjust_duplicate(u_y, 1); + + const float_type* dx_kern = make_dx_kernel<float_type>(); + const float_type* dy_kern = make_dy_kernel<float_type>(); + + static window2d full3x3; + static const bool vals[] = { 1, 1, 1, + 1, 1, 1, + 1, 1, 1 }; + + convert::from_to(vals, full3x3); + + mln_pixter(I) p(u); + mln_pixter(const M) p_m(mask); + mln_pixter(const I) p_x(u_x); + mln_pixter(const I) p_y(u_y); + + mln_qixter(const I, const window2d) q_x(p_x, full3x3); + mln_qixter(const I, const window2d) q_y(p_y, full3x3); + + for (p.start(), p_m.start(), p_x.start(), p_y.start(); + p.is_valid(); + p.next(), p_m.next(), p_x.next(), p_y.next()) + { + if (p_m.val()) + { + mln_value(I) old = p.val(); + + mln_value(I) p_x = literal::zero; + mln_value(I) p_y = literal::zero; + + unsigned k = 0; + for_all_2 (q_x, q_y) + { + p_x += q_x.val() * dx_kern[k]; + p_y += q_y.val() * dy_kern[k]; + ++k; + } + + mln_value(I) dv = p_x + p_y; + + p.val() += dv; + + interval(p.val(), 0.0, 255.0); + + float_type diff = norm::l2(old - p.val()); + + if (diff > max_diff_norm) + { + max_diff_norm = diff; + } + } + } + } + + } /* mln::inpainting */ + +} /* mln */ -- 1.7.2.5