https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud@lrde.epita.fr> Start cleaning up linear gaussian filters. * mln/linear/gaussian: New directory. * mln/linear/gaussian/internal: New directory. * mln/linear/gaussian.hh: Split code into... * mln/linear/gaussian/impl.hh, * mln/linear/gaussian/filter.hh, * mln/linear/gaussian/internal/coefficients.hh: ...those new files. * tests/linear/gaussian: New directory. * tests/linear/Makefile.am: Update. * tests/linear/gaussian/Makefile.am: New. * tests/linear/gaussian/filter.cc: New. Misc. * mln/labeling/compute.hh: . * mln/labeling/pack.hh (relabel_inplace): Fix missing renaming into... (pack_inplace): ...this. mln/labeling/compute.hh | 2 mln/labeling/pack.hh | 3 mln/linear/gaussian/filter.hh | 735 ------------------------- mln/linear/gaussian/impl.hh | 427 --------------- mln/linear/gaussian/internal/coefficients.hh | 768 ++++----------------------- tests/linear/Makefile.am | 4 tests/linear/gaussian/Makefile.am | 4 tests/linear/gaussian/filter.cc | 27 8 files changed, 185 insertions(+), 1785 deletions(-) Index: mln/linear/gaussian/impl.hh --- mln/linear/gaussian/impl.hh (revision 0) +++ mln/linear/gaussian/impl.hh (working copy) @@ -26,12 +26,12 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_LINEAR_GAUSSIAN_HH -# define MLN_LINEAR_GAUSSIAN_HH +#ifndef MLN_LINEAR_GAUSSIAN_IMPL_HH +# define MLN_LINEAR_GAUSSIAN_IMPL_HH -/// \file mln/linear/gaussian.hh +/// \file mln/linear/gaussian/impl.hh /// -/// Gaussian filter. +/// Gaussian filter implementation. /// /// \todo Add a clean reference Rachid Deriche /// Recursively implementing the gaussian and its derivatives (1993) @@ -52,6 +52,8 @@ # include <mln/level/stretch.hh> # include <mln/algebra/vec.hh> +# include <mln/linear/gaussian/internal/coefficients.hh> + namespace mln { @@ -59,117 +61,20 @@ namespace linear { - /// Gaussian filter of an image \p input - /// - /// \pre output.domain = input.domain - /// - template <typename I> - mln_concrete(I) - gaussian(const Image<I>& input, float sigma); - - - template <typename I> - mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir); - - -# ifndef MLN_INCLUDE_ONLY - - namespace impl - { - - typedef float norm_fun(float, float, - float, float, - float, float, - float, float, - float, float, - int&); - - struct recursivefilter_coef_ + namespace gaussian { - /// - /// Constructor. - /// - recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm); - std::vector<float> n, d, nm, dm; - }; - - inline - recursivefilter_coef_::recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm) - { - n.reserve(5); - d.reserve(5); - nm.reserve(5); - dm.reserve(5); - - b0 /= s; - b1 /= s; - w0 /= s; - w1 /= s; - - float sin0 = std::sin(w0); - float sin1 = std::sin(w1); - float cos0 = std::cos(w0); - float cos1 = std::cos(w1); - - 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] = - std::exp(-b1 - 2*b0) * (c1 * sin1 - cos1 * c0) + - std::exp(-b0 - 2*b1) * (a1 * sin0 - cos0 * a0); - n[2] = - 2 * std::exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 - - cos1 * a1 * sin0 - - cos0 * c1 * sin1) + - c0 * std::exp(-2*b0) + a0 * std::exp(-2*b1); - n[1] = - std::exp(-b1) * (c1 * sin1 - (c0 + 2 * a0) * cos1) + - std::exp(-b0) * (a1 * sin0 - (2 * c0 + a0) * cos0); - n[0] = - a0 + c0; - - d[4] = - std::exp(-2 * b0 - 2 * b1); - d[3] = - -2 * cos0 * std::exp(-b0 - 2*b1) - - 2 * cos1 * std::exp(-b1 - 2*b0); - d[2] = - 4 * cos1 * cos0 * std::exp(-b0 - b1) + - std::exp(-2*b1) + std::exp(-2*b0); - d[1] = - -2 * std::exp(-b1) * cos1 - 2 * std::exp(-b0) * cos0; - - for (unsigned i = 1; i <= 3; ++i) + namespace impl { - dm[i] = d[i]; - nm[i] = float(sign) * (n[i] - d[i] * n[0]); - } - dm[4] = d[4]; - nm[4] = float(sign) * (-d[4] * n[0]); - } +# ifndef MLN_INCLUDE_ONLY template <class WorkType, class I> inline void recursivefilter_(I& ima, - const recursivefilter_coef_& c, + const internal::coefficients& c, const mln_psite(I)& start, const mln_psite(I)& finish, int len, @@ -264,104 +169,6 @@ } - 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); - - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; - float scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); - - float expb1 = std::exp(b1); - float exp2b1 = std::exp(2.f * 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * 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 = std::exp(b1); - float exp2b1 = std::exp(2.f * 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); - - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; - scale0 *= scale0 * scale0; - - 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 = std::exp(b1); - float exp2b1 = std::exp(2.f * b1); - - float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; - scale1 *= scale1 * scale1; - - 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 @@ -493,10 +300,10 @@ generic_filter_common_(trait::value::nature::floating, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out) { - mln_ch_value(O, float) work_img(exact(in).domain()); + mln_ch_value(O, double) work_img(exact(in).domain()); data::paste(in, work_img); extension::adjust_fill(work_img, 4, 0); @@ -517,11 +324,11 @@ generic_filter_common_(trait::value::nature::floating, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out, int dir) { - mln_ch_value(O, float) work_img(exact(in).domain()); + mln_ch_value(O, double) work_img(exact(in).domain()); data::paste(in, work_img); extension::adjust_fill(work_img, 4, 0); @@ -542,10 +349,10 @@ generic_filter_common_(trait::value::nature::scalar, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out) { - mln_ch_value(O, float) work_img(exact(in).domain()); + mln_ch_value(O, double) work_img(exact(in).domain()); data::paste(in, work_img); extension::adjust_fill(work_img, 4, 0); @@ -566,11 +373,11 @@ generic_filter_common_(trait::value::nature::scalar, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out, int dir) { - mln_ch_value(O, float) work_img(exact(in).domain()); + mln_ch_value(O, double) work_img(exact(in).domain()); data::paste(in, work_img); extension::adjust_fill(work_img, 4, 0); @@ -592,10 +399,10 @@ generic_filter_common_(trait::value::nature::vectorial, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out) { - // typedef algebra::vec<3, float> vec3f; + // typedef algebra::vec<3, double> vec3f; // mln_ch_value(O, vec3f) work_img(exact(in).domain()); // FIXME : paste does not work (rgb8 -> vec3f). data::paste(in, out); @@ -614,11 +421,11 @@ generic_filter_common_(trait::value::nature::vectorial, const Image<I>& in, const F& coef, - float sigma, + double sigma, Image<O>& out, int dir) { - // typedef algebra::vec<3, float> vec3f; + // typedef algebra::vec<3, double> vec3f; // mln_ch_value(O, vec3f) work_img(exact(in).domain()); // FIXME : paste does not work (rgb8 -> vec3f). data::paste(in, out); @@ -630,199 +437,17 @@ out, coef, dir); } - } // end of namespace mln::linear::impl - - // Facade. - /*! Apply an approximated gaussian filter of \p sigma on \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir) - { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - 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_norm_coef_); - - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output, dir); - return output; - } - - - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma, int dir) - { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - mln_concrete(I) output; - initialize(output, input); - - 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, dir); - return output; - } - - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma, int dir) - { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - mln_concrete(I) output; - initialize(output, input); - - 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, dir); - return output; - } - - - - - - /*! Apply an approximated gaussian filter of \p sigma on \p input. - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - 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_norm_coef_); - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output); - - return output; - } - - - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - mln_concrete(I) output; - initialize(output, input); - - 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; - } - - - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); +# endif // ! MLN_INCLUDE_ONLY - mln_concrete(I) output; - initialize(output, input); - 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; - } + } // end of namespace mln::linear::gaussian::impl -# endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln::linear::gaussian } // end of namespace mln::linear } // end of namespace mln -#endif // ! MLN_LINEAR_GAUSSIAN_HH +#endif // ! MLN_LINEAR_GAUSSIAN_IMPL_HH Property changes on: mln/linear/gaussian/impl.hh ___________________________________________________________________ Added: svn:mergeinfo Index: mln/linear/gaussian/filter.hh --- mln/linear/gaussian/filter.hh (revision 0) +++ mln/linear/gaussian/filter.hh (working copy) @@ -26,31 +26,19 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_LINEAR_GAUSSIAN_HH -# define MLN_LINEAR_GAUSSIAN_HH +#ifndef MLN_LINEAR_GAUSSIAN_FILTER_HH +# define MLN_LINEAR_GAUSSIAN_FILTER_HH -/// \file mln/linear/gaussian.hh +/// \file mln/linear/gaussian/filter.hh /// /// Gaussian filter. /// /// \todo Add a clean reference Rachid Deriche /// Recursively implementing the gaussian and its derivatives (1993) -# include <vector> -# include <cmath> - # include <mln/core/concept/image.hh> -# include <mln/core/alias/point2d.hh> -# include <mln/core/alias/dpoint1d.hh> -# include <mln/core/alias/dpoint3d.hh> -# include <mln/extension/adjust_fill.hh> -# include <mln/geom/ncols.hh> -# include <mln/geom/nrows.hh> -# include <mln/geom/ninds.hh> -# include <mln/geom/nslis.hh> -# include <mln/data/paste.hh> -# include <mln/level/stretch.hh> -# include <mln/algebra/vec.hh> +# include <mln/linear/gaussian/impl.hh> +# include <mln/linear/gaussian/internal/coefficients.hh> namespace mln @@ -59,580 +47,28 @@ namespace linear { + namespace gaussian + { + /// Gaussian filter of an image \p input /// /// \pre output.domain = input.domain /// template <typename I> mln_concrete(I) - gaussian(const Image<I>& input, float sigma); + filter(const Image<I>& input, double sigma); template <typename I> mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir); - - -# ifndef MLN_INCLUDE_ONLY - - namespace impl - { - - typedef float norm_fun(float, float, - float, float, - float, float, - float, float, - float, float, - int&); - - struct recursivefilter_coef_ - { - - /// - /// Constructor. - /// - recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm); - std::vector<float> n, d, nm, dm; - }; - - inline - recursivefilter_coef_::recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm) - { - n.reserve(5); - d.reserve(5); - nm.reserve(5); - dm.reserve(5); - - b0 /= s; - b1 /= s; - w0 /= s; - w1 /= s; - - float sin0 = std::sin(w0); - float sin1 = std::sin(w1); - float cos0 = std::cos(w0); - float cos1 = std::cos(w1); - - 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] = - std::exp(-b1 - 2*b0) * (c1 * sin1 - cos1 * c0) + - std::exp(-b0 - 2*b1) * (a1 * sin0 - cos0 * a0); - n[2] = - 2 * std::exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 - - cos1 * a1 * sin0 - - cos0 * c1 * sin1) + - c0 * std::exp(-2*b0) + a0 * std::exp(-2*b1); - n[1] = - std::exp(-b1) * (c1 * sin1 - (c0 + 2 * a0) * cos1) + - std::exp(-b0) * (a1 * sin0 - (2 * c0 + a0) * cos0); - n[0] = - a0 + c0; - - d[4] = - std::exp(-2 * b0 - 2 * b1); - d[3] = - -2 * cos0 * std::exp(-b0 - 2*b1) - - 2 * cos1 * std::exp(-b1 - 2*b0); - d[2] = - 4 * cos1 * cos0 * std::exp(-b0 - b1) + - std::exp(-2*b1) + std::exp(-2*b0); - d[1] = - -2 * std::exp(-b1) * cos1 - 2 * std::exp(-b0) * cos0; - - for (unsigned i = 1; i <= 3; ++i) - { - dm[i] = d[i]; - nm[i] = float(sign) * (n[i] - d[i] * n[0]); - } - dm[4] = d[4]; - nm[4] = float(sign) * (-d[4] * n[0]); - } - - - - template <class WorkType, class I> - inline - void - recursivefilter_(I& ima, - const recursivefilter_coef_& c, - const mln_psite(I)& start, - const mln_psite(I)& finish, - int len, - const mln_deduce(I, psite, delta)& d) - { - std::vector<WorkType> tmp1(len); - std::vector<WorkType> tmp2(len); - - // The fourth degree approximation implies to have a special - // look on the four first points we consider that there is - // no signal before 0 (to be discussed) - - // -- - // Causal part - - tmp1[0] = - c.n[0] * ima(start); - - tmp1[1] = - c.n[0] * ima(start + d) - + c.n[1] * ima(start) - - c.d[1] * tmp1[0]; - - tmp1[2] = - c.n[0] * ima(start + d + d) - + c.n[1] * ima(start + d) - + c.n[2] * ima(start) - - c.d[1] * tmp1[1] - - c.d[2] * tmp1[0]; - - tmp1[3] = - c.n[0] * ima(start + d + d + d) - + c.n[1] * ima(start + d + d) - + c.n[2] * ima(start + d) - + c.n[3] * ima(start) - - c.d[1] * tmp1[2] - c.d[2] * tmp1[1] - - c.d[3] * tmp1[0]; - - mln_psite(I) current(start + d + d + d + d); - for (mln_deduce(I, site, coord) i = 4; i < len; ++i) - { - tmp1[i] = - c.n[0] * ima(current) - + c.n[1] * ima(current - d) - + c.n[2] * ima(current - d - d) - + c.n[3] * ima(current - d - d - d) - - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2] - - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4]; - current = current + d; - } - - // Non causal part - - tmp2[len - 1] = WorkType(); // FIXME : = 0, literal::zero ...? - - tmp2[len - 2] = - c.nm[1] * ima(finish); - - tmp2[len - 3] = - c.nm[1] * ima(finish - d) - + c.nm[2] * ima(finish) - - c.dm[1] * tmp2[len - 2]; - - tmp2[len - 4] = - c.nm[1] * ima(finish - d - d) - + c.nm[2] * ima(finish - d) - + c.nm[3] * ima(finish) - - c.dm[1] * tmp2[len - 3] - - c.dm[2] * tmp2[len - 2]; - - current = finish - d - d - d ; - - for (int i = len - 5; i >= 0; --i) - { - tmp2[i] = - c.nm[1] * ima(current) - + c.nm[2] * ima(current + d) - + c.nm[3] * ima(current + d + d) - + c.nm[4] * ima(current + d + d + d) - - c.dm[1] * tmp2[i + 1] - c.dm[2] * tmp2[i + 2] - - c.dm[3] * tmp2[i + 3] - c.dm[4] * tmp2[i + 4]; - current = current - d; - } - - // Combine results from causal and non-causal parts. - current = start; - for (int i = 0; i < len; ++i) - { - ima(current) = tmp1[i] + tmp2[i]; - current = current + d; - } - } - - - 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); - - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; - float scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); - - float expb1 = std::exp(b1); - float exp2b1 = std::exp(2.f * 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * 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)); + filter(const Image<I>& input, double sigma, int dir); - float expb1 = std::exp(b1); - float exp2b1 = std::exp(2.f * 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 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); - - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; - scale0 *= scale0 * scale0; - - 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 = std::exp(b1); - float exp2b1 = std::exp(2.f * b1); - - float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; - scale1 *= scale1 * scale1; - - 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 - generic_filter_(trait::image::dimension::one_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - mln_precondition(dir < I::site::dim); - - - recursivefilter_<mln_value(I)>(img, coef, - point1d(static_cast<def::coord>(-img.border())), - point1d(static_cast<def::coord>(geom::ninds(img) - 1 + - img.border())), - geom::ninds(img) + 2 * img.border(), - dpoint1d(1)); - } - - template <class I, class F> - inline - void - generic_filter_(trait::image::dimension::two_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - - mln_precondition(dir < I::site::dim); - - - if (dir == 0) - { - // Apply on rows. - for (unsigned j = 0; j < geom::ncols(img); ++j) - recursivefilter_<mln_value(I)>(img, coef, - point2d(static_cast<def::coord>(-img.border()), - static_cast<def::coord>(j)), - point2d(static_cast<def::coord>(geom::nrows(img) - 1 + - img.border()), - static_cast<def::coord>(j)), - geom::nrows(img) + 2 * img.border(), - dpoint2d(1, 0)); - } - - if (dir == 1) - { - // Apply on columns. - for (unsigned i = 0; i < geom::nrows(img); ++i) - recursivefilter_<mln_value(I)>(img, coef, - point2d(static_cast<def::coord>(i), - static_cast<def::coord>(-img.border())), - point2d(static_cast<def::coord>(i), - static_cast<def::coord>(geom::ncols(img) - 1 + - img.border())), - geom::ncols(img) + 2 * img.border(), - dpoint2d(0, 1)); - } - } - - template <class I, class F> - inline - void - generic_filter_(trait::image::dimension::three_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - mln_precondition(dir < I::site::dim); - - if (dir == 0) - { - // Apply on slices. - for (unsigned j = 0; j < geom::nrows(img); ++j) - for (unsigned k = 0; k < geom::ncols(img); ++k) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(-img.border()), - static_cast<def::coord>(j), - static_cast<def::coord>(k)), - point3d(static_cast<def::coord>(geom::nslis(img) - 1 + - img.border()), - static_cast<def::coord>(j), - static_cast<def::coord>(k)), - geom::nslis(img) + 2 * - img.border(), - dpoint3d(1, 0, 0)); - } - - - if (dir == 1) - { - // Apply on rows. - for (unsigned i = 0; i < geom::nslis(img); ++i) - for (unsigned k = 0; k < geom::ncols(img); ++k) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(-img.border()), - static_cast<def::coord>(k)), - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(geom::nrows(img) - 1 + - img.border()), - static_cast<def::coord>(k)), - geom::nrows(img) + 2 * - img.border(), - dpoint3d(0, 1, 0)); - } - - if (dir == 2) - { - // Apply on columns. - for (unsigned i = 0; i < geom::nslis(img); ++i) - for (unsigned j = 0; j < geom::nrows(img); ++i) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(j), - static_cast<def::coord>(-img.border())), - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(j), - static_cast<def::coord>(geom::ncols(img) - - 1 + img.border())), - geom::ncols(img) + 2 * - img.border(), - dpoint3d(0, 0, 1)); - } - } - - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::floating, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, i); - - // We don't need to convert work_img - data::paste(work_img, out); - } - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::floating, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, dir); - - // We don't need to convert work_img - data::paste(work_img, out); - } - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::scalar, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, i); - - // Convert work_img into result type - data::paste(level::stretch(mln_value(I)(), work_img), out); - } - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::scalar, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, dir); - - // Convert work_img into result type - data::paste(level::stretch(mln_value(I)(), work_img), out); - } - - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::vectorial, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - // typedef algebra::vec<3, float> vec3f; - // mln_ch_value(O, vec3f) work_img(exact(in).domain()); - // FIXME : paste does not work (rgb8 -> vec3f). - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - out, coef, i); - } - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::vectorial, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - // typedef algebra::vec<3, float> vec3f; - // mln_ch_value(O, vec3f) work_img(exact(in).domain()); - // FIXME : paste does not work (rgb8 -> vec3f). - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - out, coef, dir); - } +# ifndef MLN_INCLUDE_ONLY - } // end of namespace mln::linear::impl - // Facade. + // Facades. /*! Apply an approximated gaussian filter of \p sigma on \p input. * on a specific direction \p dir @@ -646,95 +82,28 @@ template <typename I> inline mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir) + filter(const Image<I>& input, double sigma, int dir) { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); + trace::entering("linear::gaussian::filter"); - 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_norm_coef_); - - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output, dir); - return output; - } - - - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma, int dir) - { mln_precondition(exact(input).is_valid()); mln_precondition(dir < I::site::dim); mln_concrete(I) output; initialize(output, input); - 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_); + internal::coefficients coef = internal::coefficients_not_derivative(sigma); impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), input, coef, sigma, output, dir); - return output; - } - - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma, int dir) - { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - mln_concrete(I) output; - initialize(output, input); - - 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, dir); + trace::exiting("linear::gaussian::filter"); return output; } - /*! Apply an approximated gaussian filter of \p sigma on \p input. * This filter is applied in all the input image direction. * @@ -743,86 +112,32 @@ template <typename I> inline mln_concrete(I) - gaussian(const Image<I>& input, float sigma) + filter(const Image<I>& input, double sigma) { - mln_precondition(exact(input).is_valid()); - - 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_norm_coef_); - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output); + trace::entering("linear::gaussian::filter"); - return output; - } - - - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma) - { mln_precondition(exact(input).is_valid()); mln_concrete(I) output; initialize(output, input); - 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; - } - - - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - mln_concrete(I) output; - initialize(output, input); + internal::coefficients coef = internal::coefficients_not_derivative(sigma); - 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); + + trace::exiting("linear::gaussian::filter"); return output; } + # endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln::linear::gaussian + } // end of namespace mln::linear } // end of namespace mln -#endif // ! MLN_LINEAR_GAUSSIAN_HH +#endif // ! MLN_LINEAR_GAUSSIAN_FILTER_HH Property changes on: mln/linear/gaussian/filter.hh ___________________________________________________________________ Added: svn:mergeinfo Index: mln/linear/gaussian/internal/coefficients.hh --- mln/linear/gaussian/internal/coefficients.hh (revision 0) +++ mln/linear/gaussian/internal/coefficients.hh (working copy) @@ -26,12 +26,12 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_LINEAR_GAUSSIAN_HH -# define MLN_LINEAR_GAUSSIAN_HH +#ifndef MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH +# define MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH -/// \file mln/linear/gaussian.hh +/// \file mln/linear/gaussian/internal/coefficients.hh /// -/// Gaussian filter. +/// Compute coefficients for recursive Gaussian filtering. /// /// \todo Add a clean reference Rachid Deriche /// Recursively implementing the gaussian and its derivatives (1993) @@ -39,18 +39,6 @@ # include <vector> # include <cmath> -# include <mln/core/concept/image.hh> -# include <mln/core/alias/point2d.hh> -# include <mln/core/alias/dpoint1d.hh> -# include <mln/core/alias/dpoint3d.hh> -# include <mln/extension/adjust_fill.hh> -# include <mln/geom/ncols.hh> -# include <mln/geom/nrows.hh> -# include <mln/geom/ninds.hh> -# include <mln/geom/nslis.hh> -# include <mln/data/paste.hh> -# include <mln/level/stretch.hh> -# include <mln/algebra/vec.hh> namespace mln @@ -59,52 +47,48 @@ namespace linear { - /// Gaussian filter of an image \p input - /// - /// \pre output.domain = input.domain - /// - template <typename I> - mln_concrete(I) - gaussian(const Image<I>& input, float sigma); - - - template <typename I> - mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir); - - -# ifndef MLN_INCLUDE_ONLY + namespace gaussian + { - namespace impl + namespace internal { - typedef float norm_fun(float, float, - float, float, - float, float, - float, float, - float, float, + + typedef double norm_fun(double, double, + double, double, + double, double, + double, double, + double, double, int&); - struct recursivefilter_coef_ + struct coefficients { - - /// /// Constructor. - /// - recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm); - std::vector<float> n, d, nm, dm; + // + coefficients(double a0, double a1, + double b0, double b1, + double c0, double c1, + double w0, double w1, + double s, norm_fun norm); + std::vector<double> n, d, nm, dm; }; + + coefficients coefficients_not_derivative(double sigma); + coefficients coefficients_1st_derivative(double sigma); + coefficients coefficients_2nd_derivative(double sigma); + + + +# ifndef MLN_INCLUDE_ONLY + + inline - recursivefilter_coef_::recursivefilter_coef_(float a0, float a1, - float b0, float b1, - float c0, float c1, - float w0, float w1, - float s, norm_fun norm) + coefficients::coefficients(double a0, double a1, + double b0, double b1, + double c0, double c1, + double w0, double w1, + double s, norm_fun norm) { n.reserve(5); d.reserve(5); @@ -116,13 +100,13 @@ w0 /= s; w1 /= s; - float sin0 = std::sin(w0); - float sin1 = std::sin(w1); - float cos0 = std::cos(w0); - float cos1 = std::cos(w1); + double sin0 = std::sin(w0); + double sin1 = std::sin(w1); + double cos0 = std::cos(w0); + double cos1 = std::cos(w1); int sign = 1; - float n_ = norm(a0, a1, b0, b1, c0, c1, cos0, sin0, cos1, sin1, sign); + double n_ = norm(a0, a1, b0, b1, c0, c1, cos0, sin0, cos1, sin1, sign); a0 /= n_; a1 /= n_; @@ -157,672 +141,160 @@ for (unsigned i = 1; i <= 3; ++i) { dm[i] = d[i]; - nm[i] = float(sign) * (n[i] - d[i] * n[0]); + nm[i] = double(sign) * (n[i] - d[i] * n[0]); } dm[4] = d[4]; - nm[4] = float(sign) * (-d[4] * n[0]); - } - - - - template <class WorkType, class I> - inline - void - recursivefilter_(I& ima, - const recursivefilter_coef_& c, - const mln_psite(I)& start, - const mln_psite(I)& finish, - int len, - const mln_deduce(I, psite, delta)& d) - { - std::vector<WorkType> tmp1(len); - std::vector<WorkType> tmp2(len); - - // The fourth degree approximation implies to have a special - // look on the four first points we consider that there is - // no signal before 0 (to be discussed) - - // -- - // Causal part - - tmp1[0] = - c.n[0] * ima(start); - - tmp1[1] = - c.n[0] * ima(start + d) - + c.n[1] * ima(start) - - c.d[1] * tmp1[0]; - - tmp1[2] = - c.n[0] * ima(start + d + d) - + c.n[1] * ima(start + d) - + c.n[2] * ima(start) - - c.d[1] * tmp1[1] - - c.d[2] * tmp1[0]; - - tmp1[3] = - c.n[0] * ima(start + d + d + d) - + c.n[1] * ima(start + d + d) - + c.n[2] * ima(start + d) - + c.n[3] * ima(start) - - c.d[1] * tmp1[2] - c.d[2] * tmp1[1] - - c.d[3] * tmp1[0]; - - mln_psite(I) current(start + d + d + d + d); - for (mln_deduce(I, site, coord) i = 4; i < len; ++i) - { - tmp1[i] = - c.n[0] * ima(current) - + c.n[1] * ima(current - d) - + c.n[2] * ima(current - d - d) - + c.n[3] * ima(current - d - d - d) - - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2] - - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4]; - current = current + d; - } - - // Non causal part - - tmp2[len - 1] = WorkType(); // FIXME : = 0, literal::zero ...? - - tmp2[len - 2] = - c.nm[1] * ima(finish); - - tmp2[len - 3] = - c.nm[1] * ima(finish - d) - + c.nm[2] * ima(finish) - - c.dm[1] * tmp2[len - 2]; - - tmp2[len - 4] = - c.nm[1] * ima(finish - d - d) - + c.nm[2] * ima(finish - d) - + c.nm[3] * ima(finish) - - c.dm[1] * tmp2[len - 3] - - c.dm[2] * tmp2[len - 2]; - - current = finish - d - d - d ; - - for (int i = len - 5; i >= 0; --i) - { - tmp2[i] = - c.nm[1] * ima(current) - + c.nm[2] * ima(current + d) - + c.nm[3] * ima(current + d + d) - + c.nm[4] * ima(current + d + d + d) - - c.dm[1] * tmp2[i + 1] - c.dm[2] * tmp2[i + 2] - - c.dm[3] * tmp2[i + 3] - c.dm[4] * tmp2[i + 4]; - current = current - d; + nm[4] = double(sign) * (-d[4] * n[0]); } - // Combine results from causal and non-causal parts. - current = start; - for (int i = 0; i < len; ++i) - { - ima(current) = tmp1[i] + tmp2[i]; - current = current + d; - } - } 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, + double norm_not_derivative(double a0, double a1, + double b0, double b1, + double c0, double c1, + double cos0, double sin0, + double cos1, double sin1, int& sign) { - float expb0 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); + double expb0 = std::exp(b0); + double exp2b0 = std::exp(2.f * b0); - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; - float scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); + double scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + double scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); - float expb1 = std::exp(b1); - float exp2b1 = std::exp(2.f * b1); + double expb1 = std::exp(b1); + double exp2b1 = std::exp(2.f * b1); - float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; - float scaleC = 2 * c1 * sin1 * expb1 - c0 * (1 - exp2b1); + double scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + double scaleC = 2 * c1 * sin1 * expb1 - c0 * (1 - exp2b1); - float sumA = scaleA / scale0; - float sumC = scaleC / scale1; + double sumA = scaleA / scale0; + double sumC = scaleC / scale1; sign = 1; - return (sumA + sumC); + 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, + double norm_1st_derivative(double a0, double a1, + double b0, double b1, + double c0, double c1, + double cos0, double sin0, + double cos1, double sin1, int& sign) { - float expb0 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); + double expb0 = std::exp(b0); + double exp2b0 = std::exp(2.f * b0); - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + double scale0 = 1 + exp2b0 - 2 * cos0 * expb0; scale0 *= scale0; - float scaleA = - 2 * a1 * sin0 * expb0 * (1 - exp2b0) + + double scaleA = - 2 * a1 * sin0 * expb0 * (1 - exp2b0) + 2 * a0 * expb0 * (2 * expb0 - cos0 * (1 + exp2b0)); - float expb1 = std::exp(b1); - float exp2b1 = std::exp(2.f * b1); + double expb1 = std::exp(b1); + double exp2b1 = std::exp(2.f * b1); - float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + double scale1 = 1 + exp2b1 - 2 * cos1 * expb1; scale1 *= scale1; - float scaleC = - 2 * c1 * sin1 * expb1 * (1 - exp2b1) + + double scaleC = - 2 * c1 * sin1 * expb1 * (1 - exp2b1) + 2 * c0 * expb1 * (2 * expb1 - cos1 * (1 + exp2b1)); - float sumA = scaleA / scale0; - float sumC = scaleC / scale1; + double sumA = scaleA / scale0; + double sumC = scaleC / scale1; sign = -1; - return (sumA + sumC); + 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, + double norm_2nd_derivative(double a0, double a1, + double b0, double b1, + double c0, double c1, + double cos0, double sin0, + double cos1, double sin1, int& sign) { - float expb0 = std::exp(b0); - float exp2b0 = std::exp(2.f * b0); + double expb0 = std::exp(b0); + double exp2b0 = std::exp(2.f * b0); - float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; + double scale0 = 1 + exp2b0 - 2 * cos0 * expb0; scale0 *= scale0 * scale0; - float scaleA = a1 * sin0 * expb0 * + double 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 = std::exp(b1); - float exp2b1 = std::exp(2.f * b1); + double expb1 = std::exp(b1); + double exp2b1 = std::exp(2.f * b1); - float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; + double scale1 = 1 + exp2b1 - 2 * cos1 * expb1; scale1 *= scale1 * scale1; - float scaleC = c1 * sin1 * expb1 * + double 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; + double sumA = scaleA / scale0; + double sumC = scaleC / scale1; sign = 1; - return (sumA + sumC); + return sumA + sumC; } - template <class I, class F> inline - void - generic_filter_(trait::image::dimension::one_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - mln_precondition(dir < I::site::dim); - - - recursivefilter_<mln_value(I)>(img, coef, - point1d(static_cast<def::coord>(-img.border())), - point1d(static_cast<def::coord>(geom::ninds(img) - 1 + - img.border())), - geom::ninds(img) + 2 * img.border(), - dpoint1d(1)); - } - - template <class I, class F> - inline - void - generic_filter_(trait::image::dimension::two_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - - mln_precondition(dir < I::site::dim); - - - if (dir == 0) - { - // Apply on rows. - for (unsigned j = 0; j < geom::ncols(img); ++j) - recursivefilter_<mln_value(I)>(img, coef, - point2d(static_cast<def::coord>(-img.border()), - static_cast<def::coord>(j)), - point2d(static_cast<def::coord>(geom::nrows(img) - 1 + - img.border()), - static_cast<def::coord>(j)), - geom::nrows(img) + 2 * img.border(), - dpoint2d(1, 0)); - } - - if (dir == 1) + coefficients coefficients_not_derivative(double sigma) { - // Apply on columns. - for (unsigned i = 0; i < geom::nrows(img); ++i) - recursivefilter_<mln_value(I)>(img, coef, - point2d(static_cast<def::coord>(i), - static_cast<def::coord>(-img.border())), - point2d(static_cast<def::coord>(i), - static_cast<def::coord>(geom::ncols(img) - 1 + - img.border())), - geom::ncols(img) + 2 * img.border(), - dpoint2d(0, 1)); - } + coefficients tmp(+1.6800, +3.7350, + +1.7830, +1.7230, + -0.6803, -0.2598, + +0.6318, +1.9970, + sigma, + norm_not_derivative); + return tmp; } - template <class I, class F> - inline - void - generic_filter_(trait::image::dimension::three_d, - Image<I>& img_, const F& coef, int dir) - { - I& img = exact(img_); - mln_precondition(dir < I::site::dim); - - if (dir == 0) - { - // Apply on slices. - for (unsigned j = 0; j < geom::nrows(img); ++j) - for (unsigned k = 0; k < geom::ncols(img); ++k) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(-img.border()), - static_cast<def::coord>(j), - static_cast<def::coord>(k)), - point3d(static_cast<def::coord>(geom::nslis(img) - 1 + - img.border()), - static_cast<def::coord>(j), - static_cast<def::coord>(k)), - geom::nslis(img) + 2 * - img.border(), - dpoint3d(1, 0, 0)); - } - - - if (dir == 1) - { - // Apply on rows. - for (unsigned i = 0; i < geom::nslis(img); ++i) - for (unsigned k = 0; k < geom::ncols(img); ++k) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(-img.border()), - static_cast<def::coord>(k)), - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(geom::nrows(img) - 1 + - img.border()), - static_cast<def::coord>(k)), - geom::nrows(img) + 2 * - img.border(), - dpoint3d(0, 1, 0)); - } - - if (dir == 2) - { - // Apply on columns. - for (unsigned i = 0; i < geom::nslis(img); ++i) - for (unsigned j = 0; j < geom::nrows(img); ++i) - recursivefilter_<mln_value(I)>(img, coef, - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(j), - static_cast<def::coord>(-img.border())), - point3d(static_cast<def::coord>(i), - static_cast<def::coord>(j), - static_cast<def::coord>(geom::ncols(img) - - 1 + img.border())), - geom::ncols(img) + 2 * - img.border(), - dpoint3d(0, 0, 1)); - } - } - - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::floating, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, i); - - // We don't need to convert work_img - data::paste(work_img, out); - } - template <class I, class F, class O> inline - void - generic_filter_common_(trait::value::nature::floating, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, dir); - - // We don't need to convert work_img - data::paste(work_img, out); - } - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::scalar, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, i); - - // Convert work_img into result type - data::paste(level::stretch(mln_value(I)(), work_img), out); - } - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::scalar, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - mln_ch_value(O, float) work_img(exact(in).domain()); - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - work_img, coef, dir); - - // Convert work_img into result type - data::paste(level::stretch(mln_value(I)(), work_img), out); - } - - - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::vectorial, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out) - { - // typedef algebra::vec<3, float> vec3f; - // mln_ch_value(O, vec3f) work_img(exact(in).domain()); - // FIXME : paste does not work (rgb8 -> vec3f). - data::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) - for (int i = 0; i < I::site::dim; ++i) - generic_filter_(mln_trait_image_dimension(I)(), - out, coef, i); - } - - template <class I, class F, class O> - inline - void - generic_filter_common_(trait::value::nature::vectorial, - const Image<I>& in, - const F& coef, - float sigma, - Image<O>& out, - int dir) - { - // typedef algebra::vec<3, float> vec3f; - // mln_ch_value(O, vec3f) work_img(exact(in).domain()); - // FIXME : paste does not work (rgb8 -> vec3f). - data::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) - generic_filter_(mln_trait_image_dimension(I)(), - out, coef, dir); - } - - } // end of namespace mln::linear::impl - - // Facade. - - /*! Apply an approximated gaussian filter of \p sigma on \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian(const Image<I>& input, float sigma, int dir) + coefficients coefficients_1st_derivative(double sigma) { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - 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_norm_coef_); - - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output, dir); - return output; + coefficients tmp(-0.6472, -4.5310, + +1.5270, +1.5160, + +0.6494, +0.9557, + +0.6719, +2.0720, + sigma, + norm_1st_derivative); + return tmp; } - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma, int dir) + coefficients coefficients_2nd_derivative(double sigma) { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - mln_concrete(I) output; - initialize(output, input); - - 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, dir); - return output; + coefficients tmp(-1.3310, +3.661, + +1.2400, +1.314, + +0.3225, -1.738, + +0.7480, +2.166, + sigma, + norm_2nd_derivative); } - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input. - * on a specific direction \p dir - * if \p dir = 0, the filter is applied on the first image dimension. - * if \p dir = 1, the filter is applied on the second image dimension. - * And so on... - * - * \pre input.is_valid - * \pre dir < dimension(input) - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma, int dir) - { - mln_precondition(exact(input).is_valid()); - mln_precondition(dir < I::site::dim); - - mln_concrete(I) output; - initialize(output, input); - - 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, dir); - return output; - } - - - - - - /*! Apply an approximated gaussian filter of \p sigma on \p input. - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - 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_norm_coef_); - impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, output); - - return output; - } - - - /*! Apply an approximated first derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_1st_derivative(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - mln_concrete(I) output; - initialize(output, input); - - 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; - } - - - /*! Apply an approximated second derivative gaussian filter of \p sigma on - * \p input - * This filter is applied in all the input image direction. - * - * \pre input.is_valid - */ - template <typename I> - inline - mln_concrete(I) - gaussian_2nd_derivative(const Image<I>& input, float sigma) - { - mln_precondition(exact(input).is_valid()); - - mln_concrete(I) output; - initialize(output, input); +# endif // ! MLN_INCLUDE_ONLY - 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; - } + } // end of namespace mln::linear::gaussian::internal -# endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln::linear::gaussian } // end of namespace mln::linear } // end of namespace mln -#endif // ! MLN_LINEAR_GAUSSIAN_HH +#endif // ! MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH Property changes on: mln/linear/gaussian/internal/coefficients.hh ___________________________________________________________________ Added: svn:mergeinfo Index: mln/labeling/compute.hh --- mln/labeling/compute.hh (revision 3725) +++ mln/labeling/compute.hh (working copy) @@ -128,7 +128,7 @@ const mln_value(L)& nlabels) { mln_precondition(exact(label).is_valid()); - mlc_is_a(mln_value(L), mln::value::Symbolic)::check(); + // mlc_is_a(mln_value(L), mln::value::Symbolic)::check(); (void) a; (void) label; (void) nlabels; Index: mln/labeling/pack.hh --- mln/labeling/pack.hh (revision 3725) +++ mln/labeling/pack.hh (working copy) @@ -34,7 +34,6 @@ /// Remove components and relabel a labeled image in order to have a /// contiguous labeling. - # include <mln/core/concept/image.hh> # include <mln/make/relabelfun.hh> @@ -71,7 +70,7 @@ // template <typename I> void - relabel_inplace(Image<I>& label, mln_value(I)& new_nlabels); + pack_inplace(Image<I>& label, mln_value(I)& new_nlabels); Index: tests/linear/Makefile.am --- tests/linear/Makefile.am (revision 3725) +++ tests/linear/Makefile.am (working copy) @@ -2,7 +2,9 @@ include $(top_srcdir)/milena/tests/tests.mk -SUBDIRS = local +SUBDIRS = \ + gaussian \ + local check_PROGRAMS = \ convolve \ Index: tests/linear/gaussian/Makefile.am --- tests/linear/gaussian/Makefile.am (revision 0) +++ tests/linear/gaussian/Makefile.am (working copy) @@ -3,8 +3,8 @@ include $(top_srcdir)/milena/tests/tests.mk check_PROGRAMS = \ - convolve + filter -convolve_SOURCES = convolve.cc +filter_SOURCES = filter.cc TESTS = $(check_PROGRAMS) Property changes on: tests/linear/gaussian/Makefile.am ___________________________________________________________________ Added: svn:mergeinfo Index: tests/linear/gaussian/filter.cc --- tests/linear/gaussian/filter.cc (revision 0) +++ tests/linear/gaussian/filter.cc (working copy) @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 EPITA Research and Development +// Laboratory (LRDE) // // This file is part of the Olena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -26,27 +26,20 @@ // reasons why the executable file might be covered by the GNU General // Public License. -/// \file tests/linear/gaussian.cc +/// \file tests/linear/gaussian/filter.cc /// -/// Test on mln::linear::gaussian. +/// Test on mln::linear::gaussian::filter. #include <mln/core/image/image2d.hh> -#include <mln/value/int_u8.hh> -#include <mln/value/int_u_sat.hh> +#include <mln/value/int_u8.hh> #include <mln/io/pgm/load.hh> #include <mln/io/pgm/save.hh> -#include <mln/level/transform.hh> -#include <mln/data/paste.hh> -#include <mln/math/round.hh> - -#include <mln/linear/gaussian.hh> +#include <mln/linear/gaussian/filter.hh> #include "tests/data.hh" -#include <mln/trace/all.hh> - int main() @@ -56,12 +49,6 @@ image2d<value::int_u8> lena; io::pgm::load(lena, MLN_IMG_DIR "/lena.pgm"); - image2d<value::int_u8> out = linear::gaussian(lena, 5.1f); + image2d<value::int_u8> out = linear::gaussian::filter(lena, 5.1f); io::pgm::save(out, "out.pgm"); - - - image2d<float> lenaf(lena.domain()); - data::paste(lena, lenaf); - - image2d<float> outf = linear::gaussian(lenaf, 5.1f); } Property changes on: tests/linear/gaussian/filter.cc ___________________________________________________________________ Added: svn:mergeinfo