--- milena/mln/inpainting/lpde.hh | 67 +++++++++++++++++++++++++---------------- 1 files changed, 41 insertions(+), 26 deletions(-) diff --git a/milena/mln/inpainting/lpde.hh b/milena/mln/inpainting/lpde.hh index 04a903a..44d036e 100644 --- a/milena/mln/inpainting/lpde.hh +++ b/milena/mln/inpainting/lpde.hh @@ -94,8 +94,9 @@ namespace mln } template <typename I, typename M> - float diffuse(I& u, - const I& f, + float diffuse(const I& u, + I& u_next, + //const I& f, const std::vector<float>& w, const M& mask) { @@ -105,7 +106,7 @@ namespace mln const unsigned derivative_kernels_size = 5; const float* derivative_kernel[derivative_kernels_size]; - + derivative_kernel[0] = make_dx_kernel(); derivative_kernel[1] = make_dy_kernel(); derivative_kernel[2] = make_dxx_kernel(); @@ -119,25 +120,26 @@ namespace mln convert::from_to(vals, full3x3); - mln_pixter(I) p(u); + 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_qixter(const I, window2d) q(p, full3x3); + mln_qixter(const I, const window2d) q(p, full3x3); float sse = 0; - + unsigned car = 0; - for_all_2 (p, p_m) - { - //if (p_m.val()) + for_all_3 (p, p_next, p_m) + { + //if (p_m.val()) { - T old = p.val(); + const T old = p.val(); if (fabs(w[1]) > activation_threshold) - p.val() *= w[1]; + p_next.val() = p.val() * w[1]; // if (fabs(w[0]) > activation_threshold) - // p.val() += w[0] * p_f.val(); + // p_next.val() += w[0] * p_f.val(); T p_derivative[derivative_kernels_size]; @@ -168,25 +170,25 @@ namespace mln //||∇u||^2 if (fabs(w[2]) > activation_threshold) - p.val() += w[2] * (p_x2 + p_y2); + p_next.val() += w[2] * (p_x2 + p_y2); // laplacian if (fabs(w[3]) > activation_threshold) - p.val() += w[3] * (p_xx + p_yy); + p_next.val() += w[3] * (p_xx + p_yy); // tr(H^2(u)) if (fabs(w[4]) > activation_threshold) - p.val() += w[4] * (p_xx2 + 2 * p_xy2 + p_yy2); + p_next.val() += w[4] * (p_xx2 + 2 * p_xy2 + p_yy2); // (∇u)^T H_u (∇u) if (fabs(w[5]) > activation_threshold) - p.val() += w[5] * (components_product(p_x2, p_xx) - + 2 * components_product(components_product(p_x, p_y), p_xy) - + components_product(p_y2, p_yy)); + p_next.val() += w[5] * (components_product(p_x2, p_xx) + + 2 * components_product(components_product(p_x, p_y), p_xy) + + components_product(p_y2, p_yy)); - sse += norm::l2(old - p.val()); + sse += norm::l2(old - p_next.val()); } - } + } return sse; } @@ -195,18 +197,31 @@ namespace mln void lpde::operator()(I<T>& src, I<bool>& mask) { - extension::adjust_duplicate(src, 1); + //const I<T> f(src); + + I<T>& u_ = src; + I<T> u_next_(src); + + extension::adjust_duplicate(u_, 1); + extension::adjust_duplicate(u_next_, 1); - const I<T> f(src); + I<T>* u = &u_; + I<T>* u_next = &u_next_; - I<T>& u = src; + float sse = std::numeric_limits<float>::max() / 2.0; - float sse = std::numeric_limits<float>::max(); + float sse_old = std::numeric_limits<float>::max(); float t = 0; - while (t < this->t_f && sse > 0.01 && !std::isnan(sse)) + while (t < this->t_f && sse < sse_old && !std::isnan(sse)) { - sse = diffuse(u, f, w, mask); + sse_old = sse; + + sse = diffuse(*u, *u_next, /* f, */ w, mask); + + I<T>* swap = u; + u = u_next; + u_next = swap; std::cout << "sse = " << sse << std::endl; -- 1.7.2.5