[Olena-patches] olena-2.0-39-g8b1b95c Add: forward and backward difference kenrels for TV method, replaced tv method with primal dual version.

--- milena/mln/inpainting/derivative_kernels.hh | 91 +++++-- milena/mln/inpainting/tv.hh | 372 ++++++++++++++++++++++---- 2 files changed, 381 insertions(+), 82 deletions(-) diff --git a/milena/mln/inpainting/derivative_kernels.hh b/milena/mln/inpainting/derivative_kernels.hh index a7f3c34..04e1535 100644 --- a/milena/mln/inpainting/derivative_kernels.hh +++ b/milena/mln/inpainting/derivative_kernels.hh @@ -1,30 +1,71 @@ +#ifndef DERIVATIVE_KERNELS_HH_ +# define DERIVATIVE_KERNELS_HH_ + namespace mln { namespace inpainting { - namespace - { - const float a = 61; - const float b = 17; - - 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 }; + namespace internal + { + const float a = 2; + const float b = 1; + } + + template <typename T> + const T* make_dx_kernel() + { + static const T ws[] = { internal::b, internal::a, internal::b, + 0, 0, 0, + -internal::b, -internal::a, -internal::b }; + return ws; + } + + template <typename T> + const T* make_dy_kernel() + { + static const T ws[] = { internal::b, 0, -internal::b, + internal::a, 0, -internal::a, + internal::b, 0, -internal::b }; - return ws; - } + return ws; + } + + template <typename T> + const T* make_forward_horizontal_kernel() + { + static const T ws[] = { 0, 0, 0, + 0, -1, 1, + 0, 0, 0 }; + return ws; + } + + template <typename T> + const T* make_forward_vertical_kernel() + { + static const T ws[] = { 0, 0, 0, + 0, -1, 0, + 0, 1, 0 }; + + return ws; + } + + template <typename T> + const T* make_backward_horizontal_kernel() + { + static const T ws[] = { 0, 0, 0, + -1, 1, 0, + 0, 0, 0 }; + return ws; + } + + template <typename T> + const T* make_backward_vertical_kernel() + { + static const T ws[] = { 0, -1, 0, + 0, 1, 0, + 0, 0, 0 }; + + return ws; } template <typename T> @@ -51,10 +92,12 @@ namespace mln const T* make_dxy_kernel() { static const T ws[] = {0.1051, 0, -0.1051, - 0, 0, 0, - -0.1051, 0, 0.1051}; + 0, 0, 0, + -0.1051, 0, 0.1051}; return ws; } } } + +#endif diff --git a/milena/mln/inpainting/tv.hh b/milena/mln/inpainting/tv.hh index d5bed23..44bf24b 100644 --- a/milena/mln/inpainting/tv.hh +++ b/milena/mln/inpainting/tv.hh @@ -1,8 +1,16 @@ +#ifndef TV_HH_ +# define TV_HH_ + #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/l1.hh> #include <mln/norm/l2.hh> +#include <mln/inpainting/derivative_kernels.hh> + +#include <mln/extension/adjust_duplicate.hh> +#include <mln/inpainting/fast_isotropic_diffusion.hh> namespace mln { @@ -10,84 +18,332 @@ namespace mln { struct tv { + public: + tv(unsigned max_iterations = 10000) + : max_iterations_ (max_iterations) + { + } + typedef float float_type; + + template <template <typename> class I, typename T> + void operator()(I<T>& src, + const I<bool>& mask); + + private: + unsigned max_iterations_; }; - template <typename T> - inline - void interval(T& v, - tv::float_type min, - tv::float_type max) + namespace internal { - for (unsigned i = 0; i < mln_dim(T); ++i) - v[i] = std::max(std::min(v[i], max), min); + 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 T> + inline + void normalize_linf(T& u, + T& v) + { + for (unsigned i = 0; i < mln_dim(T); ++i) + { + float norm = std::sqrt(u[i] * u[i] + v[i] * v[i]); + + norm = std::max(1.0f, norm); + + u[i] /= norm; + v[i] /= norm; + } + } + + template <typename I, typename T> + void linear(I& y, + const I& x, + const T& a, + const T& b, + const mln_value(I)& c = literal::zero) + { + mln_pixter(const I) p_x(x); + mln_pixter(I) p_y(y); + + for_all_2 (p_x, p_y) + { + p_y.val() = a * p_y.val() + b * p_x.val() + c; + } + } + + template <typename I> + void convolve(const I& src, + I& ima, + const float* kernel) + { + static window2d full3x3; + static const bool vals[] = { 1, 1, 1, + 1, 1, 1, + 1, 1, 1 }; + + convert::from_to(vals, full3x3); + + mln_pixter(I) o(ima); + mln_pixter(const I) p(src); + mln_qixter(const I, const window2d) q(p, full3x3); + + for_all_2 (p, o) + { + mln_value(I) accu = literal::zero; + + unsigned k = 0; + for_all (q) + { + accu += q.val() * kernel[k]; + ++k; + } + + o.val() = accu; + } + } + + // debug + template <typename I, typename M> + void print_first_value(const I& src, const M& mask) + { + for (unsigned i = 0; i < src.nelements(); ++i) + { + point2d p = src.point_at_index(i); + + if (mask(p)) + { + std::cout << p << std::endl; + std::cout << src(p) << std::endl; + return; + } + } + } + + template <typename I> + float max_diff_norm(const I& u, const I& old) + { + float max = 0; + + mln_pixter(const I) p(u); + mln_pixter(const I) p_old(old); + + for_all_2 (p, p_old) + { + max = std::max(max, norm::l1(p.val() - p_old.val())); + } + + return max * max; // l2 norm + } } - 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) + // TV regularization inpainting with PrimalDual algorithm + template <template <typename> class I, typename T> + void tv::operator()(I<T>& src, + const I<bool>& mask) { - typedef tv::float_type float_type; + typedef float float_type; - extension::adjust_duplicate(u_x, 1); - extension::adjust_duplicate(u_y, 1); + const float_type* dx_bkern = make_backward_horizontal_kernel<float_type>(); + const float_type* dy_bkern = make_backward_vertical_kernel<float_type>(); + const float_type* dx_fkern = make_forward_horizontal_kernel<float_type>(); + const float_type* dy_fkern = make_forward_vertical_kernel<float_type>(); - const float_type* dx_kern = make_dx_kernel<float_type>(); - const float_type* dy_kern = make_dy_kernel<float_type>(); + const float_type sigma = 1.0f / std::sqrt(8.0f); + + I<T> u_x; + initialize(u_x, src); - static window2d full3x3; - static const bool vals[] = { 1, 1, 1, - 1, 1, 1, - 1, 1, 1 }; + I<T> u_y; + initialize(u_y, src); - convert::from_to(vals, full3x3); + I<T> u_x2; + initialize(u_x2, src); + I<T> u_y2; + initialize(u_y2, src); + + I<T> u = duplicate(src); + + I<T> old = duplicate(u); + + I<T> u_aux = duplicate(u); + + internal::convolve(u, u_x, dx_fkern); + internal::convolve(u, u_y, dy_fkern); + + float_type max_diff = std::numeric_limits<float_type>::max(); + + unsigned iterations = 0; + + while ((max_diff > CONVERGENCE_THRESHOLD) + && (iterations < this->max_iterations_)) + { + internal::convolve(u_aux, u_x2, dx_fkern); + internal::convolve(u_aux, u_y2, dy_fkern); + + internal::linear(u_x, u_x2, 1.0f, sigma); + internal::linear(u_y, u_y2, 1.0f, sigma); + + mln_pixter(I<T>) p_x(u_x); + mln_pixter(I<T>) p_y(u_y); + + for_all_2 (p_x, p_y) + { + internal::normalize_linf(p_x.val(), p_y.val()); + } + + // divergence + { + internal::convolve(u_x, u_x2, dx_bkern); + internal::convolve(u_y, u_y2, dy_bkern); + + internal::linear(u, u_x2, 0.f, 1.f); + internal::linear(u, u_y2, 1.f, 1.f); + } - 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); + internal::linear(u, u, sigma, 0.0f); + internal::linear(u, old, 1.0f, 1.0f); - mln_qixter(const I, const window2d) q_x(p_x, full3x3); - mln_qixter(const I, const window2d) q_y(p_y, full3x3); + mln_pixter(I<T>) p(u); + mln_pixter(const I<T>) p_src(src); + mln_pixter(const I<bool>) p_m(mask); - 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()) + for_all_3 (p, p_src, p_m) + { + if (!p_m.val()) + { + p.val() = p_src.val(); + } + else + { + internal::interval(p.val(), 0.0f, 1.0f); + } + } + + internal::linear(u_aux, u, 0.0f, 2.0f); + internal::linear(u_aux, old, 1.0f, -1.0f); + + max_diff = internal::max_diff_norm(u, old); + + data::paste(u, old); + + ++iterations; + } + + std::cout << iterations << "\t"; + + data::paste(u, src); + } + + template <template <typename> class I, typename T> + struct tv_step + { + public: + typedef float float_type; + + tv_step(I<T>& src, + const I<bool>& mask_) + : mask (mask_) + { + this->src = duplicate(src); + + initialize(u_x, src); + initialize(u_y, src); + initialize(u_x2, src); + initialize(u_y2, src); + + u = duplicate(src); + u_aux = duplicate(u); + + dx_bkern = make_backward_horizontal_kernel<float_type>(); + dy_bkern = make_backward_vertical_kernel<float_type>(); + dx_fkern = make_forward_horizontal_kernel<float_type>(); + dy_fkern = make_forward_vertical_kernel<float_type>(); + + sigma = 1.0f / std::sqrt(8.0f); + } + + void operator()(I<T>& ima) + { + internal::convolve(u_aux, u_x2, dx_fkern); + internal::convolve(u_aux, u_y2, dy_fkern); + + internal::linear(u_x, u_x2, 1.0f, sigma); + internal::linear(u_y, u_y2, 1.0f, sigma); + + mln_pixter(I<T>) p_x(u_x); + mln_pixter(I<T>) p_y(u_y); + + for_all_2 (p_x, p_y) + { + internal::normalize_linf(p_x.val(), p_y.val()); + } + + // divergence { - if (p_m.val()) - { - mln_value(I) old = p.val(); + internal::convolve(u_x, u_x2, dx_bkern); + internal::convolve(u_y, u_y2, dy_bkern); + + internal::linear(u, u_x2, 0.f, 1.f); + internal::linear(u, u_y2, 1.f, 1.f); + } - mln_value(I) p_x = literal::zero; - mln_value(I) p_y = literal::zero; + internal::linear(u, u, sigma, 0.0f); + internal::linear(u, ima, 1.0f, 1.0f); - unsigned k = 0; - for_all_2 (q_x, q_y) + mln_pixter(I<T>) p(u); + mln_pixter(const I<T>) p_src(src); + mln_pixter(const I<bool>) p_m(mask); + + for_all_3 (p, p_src, p_m) + { + if (!p_m.val()) { - p_x += q_x.val() * dx_kern[k]; - p_y += q_y.val() * dy_kern[k]; - ++k; + p.val() = p_src.val(); } - - 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; - } - } - } - } + else + { + internal::interval(p.val(), 0.0f, 1.0f); + } + } + + internal::linear(u_aux, u, 0.0f, 2.0f); + internal::linear(u_aux, ima, 1.0f, -1.0f); + + float max_diff = internal::max_diff_norm(u, ima); + + data::paste(u, ima); + } + + private: + I<T> src; + const I<bool>& mask; + + I<T> u_x; + I<T> u_y; + + I<T> u_x2; + I<T> u_y2; + + I<T> u; + I<T> u_aux; + + const float_type* dx_bkern; + const float_type* dy_bkern; + const float_type* dx_fkern; + const float_type* dy_fkern; + + float_type sigma; + }; } /* mln::inpainting */ } /* mln */ + +#endif -- 1.7.2.5