From: Florian Lesaint <florian.lesaint@lrde.epita.fr> To: olena-patches@lrde.epita.fr Subject: milena r2731: Add gaussian_1st_derivative and gaussian_2nd_derivative URL: https://svn.lrde.epita.fr/svn/oln/branches/cleanup-2008/milena ChangeLog: 2008-10-30 Florian Lesaint <florian.lesaint@lrde.epita.fr> Add gaussian_1st_derivative and gaussian_2nd_derivative. * mln/extension/adjust_fill.hh: Add a adjust_fill version which takes a delta. * mln/level/stretch.hh: Fix code. * mln/level/transform.hh: Revert work done by Nicolas Ballas to make gaussian work. * mln/linear/gaussian.hh: Add gaussian_1st_derivative and gaussian_2nd_derivative. * tests/linear/gaussian.cc: Fix test. --- extension/adjust_fill.hh | 27 ++++++ level/stretch.hh | 4 level/transform.hh | 11 -- linear/gaussian.hh | 199 +++++++++++++++++++++++++++++++++++++++-------- 4 files changed, 202 insertions(+), 39 deletions(-) Index: branches/cleanup-2008/milena/tests/linear/gaussian.cc =================================================================== Index: branches/cleanup-2008/milena/mln/level/transform.hh =================================================================== --- branches/cleanup-2008/milena/mln/level/transform.hh (revision 2730) +++ branches/cleanup-2008/milena/mln/level/transform.hh (revision 2731) @@ -146,26 +146,23 @@ } // end of namespace mln::level::impl::generic + } // end of namespace mln::level::impl // Facade. - namespace internal - { template <typename I, typename F, typename O> inline - void transform_(const Image<I>& input, const Function_v2v<F>& f, + void transform(const Image<I>& input, const Function_v2v<F>& f, Image<O>& output) { trace::entering("level::transform"); mln_precondition(exact(output).domain() >= exact(input).domain()); - transform_dispatch_(exact(input), exact(f), exact(output)); + impl::internal::transform_dispatch_(exact(input), exact(f), exact(output)); trace::exiting("level::transform"); } - } - } // end of namespace mln::level::impl template <typename I, typename F> @@ -178,7 +175,7 @@ mln_precondition(exact(input).has_data()); mln_ch_value(I, mln_result(F)) output; initialize(output, input); - impl::internal::transform_(input, f, output); + transform(input, f, output); trace::exiting("level::transform"); return output; Index: branches/cleanup-2008/milena/mln/level/stretch.hh =================================================================== --- branches/cleanup-2008/milena/mln/level/stretch.hh (revision 2730) +++ branches/cleanup-2008/milena/mln/level/stretch.hh (revision 2731) @@ -70,8 +70,6 @@ { trace::entering("level::impl::stretch"); - initialize(output, input); - mln_value(I) min_, max_; estim::min_max(input, min_, max_); if (max_ != min_) @@ -97,6 +95,8 @@ { trace::entering("level::stretch"); + initialize(output, input); + mln_precondition(exact(output).domain() == exact(input).domain()); impl::stretch(mln_value(O)(), input, output); Index: branches/cleanup-2008/milena/mln/linear/gaussian.hh =================================================================== --- branches/cleanup-2008/milena/mln/linear/gaussian.hh (revision 2730) +++ branches/cleanup-2008/milena/mln/linear/gaussian.hh (revision 2731) @@ -32,6 +32,9 @@ /*! \file mln/linear/gaussian.hh * * \brief Gaussian filter. + * + * \todo Add a clean reference David Deriche + * Recursively implementing the gaussian and its derivatives (1993) */ # include <mln/core/concept/image.hh> @@ -40,6 +43,9 @@ # include <mln/geom/ncols.hh> # include <mln/geom/nrows.hh> +# include <mln/extension/adjust_fill.hh> +# include <mln/level/stretch.hh> + # include <mln/algebra/vec.hh> # include <vector> @@ -64,6 +70,13 @@ namespace impl { + typedef float norm_fun(float, float, + float, float, + float, float, + float, float, + float, float, + int&); + struct recursivefilter_coef_ { @@ -74,9 +87,8 @@ float b0, float b1, float c0, float c1, float w0, float w1, - float s); + float s, norm_fun norm); std::vector<float> n, d, nm, dm; - float sumA, sumC; }; inline @@ -84,7 +96,7 @@ float b0, float b1, float c0, float c1, float w0, float w1, - float s) + float s, norm_fun norm) { n.reserve(5); d.reserve(5); @@ -101,20 +113,13 @@ float cos0 = cos(w0); float cos1 = cos(w1); - sumA = - (2.0 * a1 * exp( b0 ) * cos0 * cos0 - a0 * sin0 * exp( 2.0 * b0 ) - + a0 * sin0 - 2.0 * a1 * exp( b0 )) / - (( 2.0 * cos0 * exp( b0 ) - exp( 2.0 * b0 ) - 1 ) * sin0); - - sumC = - (2.0 * c1 * exp( b1 ) * cos1 * cos1 - c0 * sin1 * exp( 2.0 * b1 ) - + c0 * sin1 - 2.0 * c1 * exp( b1 )) - / (( 2.0 * cos1 * exp( b1 ) - exp( 2.0 * b1 ) - 1 ) * sin1); - - a0 /= (sumA + sumC); - a1 /= (sumA + sumC); - c0 /= (sumA + sumC); - c1 /= (sumA + sumC); + int sign = 1; + float n_ = norm(a0, a1, b0, b1, c0, c1, cos0, sin0, cos1, sin1, sign); + + a0 /= n_; + a1 /= n_; + c0 /= n_; + c1 /= n_; n[3] = exp( -b1 - 2*b0 ) * (c1 * sin1 - cos1 * c0) + @@ -130,7 +135,8 @@ n[0] = a0 + c0; - d[4] = exp(-2 * b0 - 2 * b1); + d[4] = + exp(-2 * b0 - 2 * b1); d[3] = -2 * cos0 * exp(-b0 - 2*b1) - 2 * cos1 * exp(-b1 - 2*b0); @@ -143,10 +149,10 @@ for (unsigned i = 1; i <= 3; ++i) { dm[i] = d[i]; - nm[i] = n[i] - d[i] * n[0]; + nm[i] = sign * (n[i] - d[i] * n[0]); } dm[4] = d[4]; - nm[4] = -d[4] * n[0]; + nm[4] = sign * (-d[4] * n[0]); } @@ -251,10 +257,101 @@ } + inline + float gaussian_norm_coef_(float a0, float a1, + float b0, float b1, + float c0, float c1, + float cos0, float sin0, + float cos1, float sin1, + int& sign) + { + float expb0 = exp(b0); + float exp2b0 = exp(2.0 * b0); + + float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + float scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); + + float expb1 = exp(b1); + float exp2b1 = exp(2.0 * b1); + + float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + float scaleC = 2 * c1 * sin1 * expb1 - c0 * (1 - exp2b1); + + float sumA = scaleA / scale0; + float sumC = scaleC / scale1; + + sign = 1; + + return (sumA + sumC); + } + + inline + float gaussian_1st_deriv_coef_norm_(float a0, float a1, + float b0, float b1, + float c0, float c1, + float cos0, float sin0, + float cos1, float sin1, + int& sign) + { + float expb0 = exp(b0); + float exp2b0 = exp(2.0 * b0); + + float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + scale0 *= scale0; + float scaleA = - 2 * a1 * sin0 * expb0 * (1 - exp2b0) + + 2 * a0 * expb0 * (2 * expb0 - cos0 * (1 + exp2b0)); + + float expb1 = exp(b1); + float exp2b1 = exp(2.0 * b1); + + float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + scale1 *= scale1; + float scaleC = - 2 * c1 * sin1 * expb1 * (1 - exp2b1) + + 2 * c0 * expb1 * (2 * expb1 - cos1 * (1 + exp2b1)); + + float sumA = scaleA / scale0; + float sumC = scaleC / scale1; + + sign = -1; + + return (sumA + sumC); + } + + inline + float gaussian_2nd_deriv_coef_norm_(float a0, float a1, + float b0, float b1, + float c0, float c1, + float cos0, float sin0, + float cos1, float sin1, + int& sign) + { + float expb0 = exp(b0); + float exp2b0 = exp(2.0 * b0); + + float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + + float scaleA = a1 * sin0 * expb0 * (1 + expb0 * (2 * cos0 * (1 + exp2b0) + exp2b0 - 6)) + + a0 * expb0 * (2 * expb0 * (2 - cos0 * cos0) * (1 - exp2b0) - cos0 * (1 - exp2b0 * exp2b0)); + + float expb1 = exp(b1); + float exp2b1 = exp(2.0 * b1); + + float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + + float scaleC = c1 * sin1 * expb1 * (1 + expb1 * (2 * cos1 * (1 + exp2b1) + exp2b1 - 6)) + + c0 * expb1 * (2 * expb1 * (2 - cos1 * cos1) * (1 - exp2b1) - cos1 * (1 - exp2b1 * exp2b1)); + + float sumA = scaleA / scale0; + float sumC = scaleC / scale1; + sign = 1; + + return (sumA + sumC); + } + template <class I, class F> inline void - gaussian_(Image<I>& img_, const F& coef) + generic_filter_(Image<I>& img_, const F& coef) { I& img = exact(img_); @@ -273,14 +370,13 @@ point2d(i, geom::ncols(img) - 1 + img.border()), geom::ncols(img) + 2 * img.border(), dpoint2d(0, 1)); - } template <class I, class F, class O> inline void - gaussian_common_(trait::value::nature::scalar, + generic_filter_common_(trait::value::nature::scalar, const Image<I>& in, const F& coef, float sigma, @@ -288,22 +384,24 @@ { mln_ch_value(O, float) work_img(exact(in).domain()); level::paste(in, work_img); + extension::adjust_fill(work_img, 4, 0); // On tiny sigma, Derich algorithm doesn't work. // It is the same thing that to convolve with a Dirac. if (sigma > 0.006) - gaussian_(work_img, coef); + generic_filter_(work_img, coef); /* Convert the result image to the user-requested datatype. FIXME: We are making an unnecessary copy in case the user expects a ntg::float_s image. */ - level::paste(work_img, out); + level::stretch(work_img, out); + //level::paste(work_img, out); } template <class I, class F, class O> inline void - gaussian_common_(trait::value::nature::vectorial, + generic_filter_common_(trait::value::nature::vectorial, const Image<I>& in, const F& coef, float sigma, @@ -314,8 +412,10 @@ // FIXME : paste does not work (rgb8 -> vec3f). level::paste(in, out); + // On tiny sigma, Derich algorithm doesn't work. + // It is the same thing that to convolve with a Dirac. if (sigma > 0.006) - gaussian_(out, coef); + generic_filter_(out, coef); } } // end of namespace mln::linear::impl @@ -330,20 +430,59 @@ mln_precondition(exact(input).has_data()); mln_concrete(I) output; - initialize(output, input); impl::recursivefilter_coef_ coef(1.68f, 3.735f, 1.783f, 1.723f, -0.6803f, -0.2598f, 0.6318f, 1.997f, - sigma); - impl::gaussian_common_(mln_trait_value_nature(mln_value(I))(), + sigma, impl::gaussian_norm_coef_); + impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), input, coef, sigma, output); return output; } + template <class I> + inline + mln_concrete(I) + gaussian_1st_derivative(const Image<I>& input, float sigma) + { + mln_precondition(exact(input).has_data()); + + mln_concrete(I) output; + + impl::recursivefilter_coef_ + coef(-0.6472f, -4.531f, + 1.527f, 1.516f, + 0.6494f, 0.9557f, + 0.6719f, 2.072f, + sigma, impl::gaussian_1st_deriv_coef_norm_); + impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), + input, coef, sigma, output); + return output; + } + + template <class I> + inline + mln_concrete(I) + gaussian_2nd_derivative(const Image<I>& input, float sigma) + { + mln_precondition(exact(input).has_data()); + + mln_concrete(I) output; + + impl::recursivefilter_coef_ + coef(-1.331f, 3.661f, + 1.24f, 1.314f, + 0.3225f, -1.738f, + 0.748f, 2.166f, + sigma, impl::gaussian_2nd_deriv_coef_norm_); + impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), + input, coef, sigma, output); + return output; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::linear Index: branches/cleanup-2008/milena/mln/extension/adjust_fill.hh =================================================================== --- branches/cleanup-2008/milena/mln/extension/adjust_fill.hh (revision 2730) +++ branches/cleanup-2008/milena/mln/extension/adjust_fill.hh (revision 2731) @@ -71,12 +71,28 @@ const Neighborhood<N>& nbh, const mln_value(I)& val); + template <typename I> + void adjust_fill(const Image<I>& ima, + unsigned delta, + const mln_value(I)& val); # ifndef MLN_INCLUDE_ONLY namespace impl { + template <typename I, typename V> + void do_adjust_fill(const I& ima, + unsigned delta, + const V& val) + { + mln_precondition(exact(ima).has_data()); + // mln_precondition(exact(win_like).is_valid()); + + border::adjust(ima, delta); + extension::fill(ima, val); + } + template <typename I, typename W, typename V> void do_adjust_fill(const I& ima, const W& win_like, @@ -124,6 +140,17 @@ trace::exiting("extension::adjust_fill"); } + template <typename I> + void adjust_fill(const Image<I>& ima, + unsigned delta, + const mln_value(I)& val) + { + trace::entering("extension::adjust_fill"); + impl::do_adjust_fill(ima, delta, val); + trace::exiting("extension::adjust_fill"); + } + + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::extension