---
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