Index: olena/ChangeLog from Niels Van Vliet <niels@lrde.epita.fr> * olena/oln/core/abstract/behavior.hh (adapt_border): Change constness of the image modified to const. * olena/oln/core/behavior.hh: Likewise. * olena/oln/convol/slow_gaussian.hh: Add behavior. Remove parameter mu. Add doc and example. Add radius_factor. * olena/oln/convol/slow_gaussian.hxx: Clean code. Remove parameter mu. Add radius_factor. * olena/oln/morpho/gradient.inc: Fix doc. Index: olena/oln/morpho/gradient.inc --- olena/oln/morpho/gradient.inc Sun, 14 Mar 2004 19:03:34 +0100 van-vl_n (oln/43_gradient.i 1.16 600) +++ olena/oln/morpho/gradient.inc Wed, 09 Jun 2004 10:48:12 +0200 van-vl_n (oln/43_gradient.i 1.16 600) @@ -71,7 +71,7 @@ ** IMG_OUT "oln_morpho_beucher_gradient.pbm"); ** return 0; ** } -** \encode +** \endcode ** ** \image html lena128_pgm.png ** \image latex lena128_pgm.png Index: olena/oln/core/abstract/behavior.hh --- olena/oln/core/abstract/behavior.hh Thu, 11 Mar 2004 17:12:19 +0100 thivol_d (oln/j/46_behavior.h 1.3 600) +++ olena/oln/core/abstract/behavior.hh Wed, 09 Jun 2004 19:01:52 +0200 van-vl_n (oln/j/46_behavior.h 1.3 600) @@ -58,7 +58,7 @@ ** Adapt the border of an image regarding the kind of behavior wanted. */ template <class I> - void adapt_border(oln::abstract::image<I> &im, coord border_size) const + void adapt_border(const oln::abstract::image<I> &im, coord border_size) const { mlc_dispatch(adapt_border)(im, border_size); }; Index: olena/oln/core/behavior.hh --- olena/oln/core/behavior.hh Fri, 12 Mar 2004 20:17:58 +0100 thivol_d (oln/j/47_behavior.h 1.3 600) +++ olena/oln/core/behavior.hh Wed, 09 Jun 2004 19:03:10 +0200 van-vl_n (oln/j/47_behavior.h 1.3 600) @@ -47,7 +47,8 @@ typedef mlc_exact_vt_type(self_type, Exact) exact_type; template <class I> - void adapt_border_impl(oln::abstract::image<I> &im, coord border_size) const + void adapt_border_impl(const oln::abstract::image<I> &im, + coord border_size) const { im.border_adapt_mirror(border_size); }; @@ -75,9 +76,11 @@ }; template <class I> - void adapt_border_impl(abstract::image<I> &im, coord border_size) const + void adapt_border_impl(const abstract::image<I> &im, + coord border_size) const { - im.border_adapt_assign(border_size, ntg::cast::force<oln_value_type(I)>(value_)); + im.border_adapt_assign(border_size, + ntg::cast::force<oln_value_type(I)>(value_)); }; protected: @@ -101,7 +104,8 @@ typedef mlc_exact_vt_type(self_type, Exact) exact_type; template <class I> - void adapt_border_impl(abstract::image<I> &im, coord border_size) const + void adapt_border_impl(const abstract::image<I> &im, + coord border_size) const { im.border_adapt_copy(border_size); }; Index: olena/oln/convol/slow_gaussian.hh --- olena/oln/convol/slow_gaussian.hh Wed, 12 May 2004 23:10:18 +0200 van-vl_n (oln/r/7_slow_gauss 1.1 600) +++ olena/oln/convol/slow_gaussian.hh Wed, 09 Jun 2004 19:14:55 +0200 van-vl_n (oln/r/7_slow_gauss 1.1 600) @@ -31,6 +31,7 @@ # include <oln/core/image.hh> # include <ntg/float.hh> # include "convolution.hh" +# include <oln/core/behavior.hh> namespace oln { namespace convol { @@ -45,33 +46,67 @@ { }; - /*! Gaussian filter. + /*! Slow Gaussian filter. ** - ** \todo FIXME: this class has been done on hurry. ** \arg in input image. ** \arg sigma sigma. - ** \arg mu point corresponding to the center. + ** \arg radius_factor the kernel of the convolution is done with a + ** kernel null at a distance of sigma * radius_factor. + ** \arg behavior Object to know how to work on borders. By default + ** mirror is used. + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/convol/slow_gaussian.hh> + ** #include <ntg/int.hh> + ** + ** using namespace oln; + ** using namespace ntg; + ** + ** int main() + ** { + ** image2d<int_u8> lena(IMG_IN "lena.pgm"); + ** + ** float_d sigma = 2.5; + ** lena = convol::slow::gaussian(lena, sigma, 3); + ** save(lena, IMG_OUT "oln_convol_slow_gaussian.pgm"); + ** } + ** \endcode + ** \image html lena_pgm.png + ** \image latex lena_pgm.png + ** => + ** \image html oln_convol_slow_gaussian.pgm + ** \image latex oln_convol_slow_gaussian.pgm */ - template <class I> + + template <class I, class BE> oln_concrete_type(I) gaussian(const oln::abstract::non_vectorial_image<I> & in, const ntg::float_d sigma, - const oln_point_type(I) & mu = oln_point_type(I)()) + const ntg::float_d radius_factor, + const oln::abstract::behavior<BE> &behavior) { - static const typename gaussian_kernel<I::dim>::ret k - = gaussian_kernel<I::dim>::kernel_norm(sigma, mu); - in.border_adapt_mirror(k.delta()); - //FIXME: we should use fast convol. - //FIXME: This is hugely. - typename mute<I, double>::ret im = oln::convol::slow::convolve<double> - (in, - gaussian_kernel<I::dim>::kernel(sigma, mu)); + const typename gaussian_kernel<I::dim>::ret k + = gaussian_kernel<I::dim>::kernel(sigma, radius_factor); + behavior.adapt_border(in, coord(k.delta())); + + typename mute<I, ntg::float_d>::ret im = + oln::convol::slow::convolve<ntg::float_d> (in, k); oln_concrete_type(I) out(in.size()); oln_iter_type(I) it(out); for_all(it) out[it] = oln_value_type(I)(im[it]); return out; } + + //! Gaussian filter, with default borders. + template <class I> + oln_concrete_type(I) + gaussian(const oln::abstract::non_vectorial_image<I> & in, + const ntg::float_d sigma, + const ntg::float_d radius_factor) + { + return gaussian(in, sigma, radius_factor, mirror_behavior<>()); + } }//End namespace slow }// End namespace convol }// End namespace oln Index: olena/oln/convol/slow_gaussian.hxx --- olena/oln/convol/slow_gaussian.hxx Wed, 12 May 2004 23:10:18 +0200 van-vl_n (oln/r/8_slow_gauss 1.1 600) +++ olena/oln/convol/slow_gaussian.hxx Wed, 09 Jun 2004 18:31:54 +0200 van-vl_n (oln/r/8_slow_gauss 1.1 600) @@ -32,9 +32,6 @@ # include <oln/core/w_window2d.hh> # include <oln/core/w_window3d.hh> -//! SIGMA_COEF times sigma is the size of the kernel. -# define SIGMA_COEF 3 - namespace oln { namespace convol { namespace slow { @@ -44,14 +41,15 @@ and w_windownd does not supprot concrete_type, etc */ template<class T> - T + inline T normalise(const T &in) { T w = in; - double sum = 0; + ntg::float_d sum = 0.; for (unsigned i = 0; i < in.card(); ++i) sum += in.w(i); sum = 1. / sum; + assert(finite(sum)); for (unsigned i = 0; i < in.card(); ++i) w.set(w.dp(i), in.w(i) * sum); return w; @@ -63,26 +61,38 @@ template<> struct gaussian_kernel<1> { - typedef oln::w_window1d<double> ret; + typedef oln::w_window1d<ntg::float_d> ret; enum {dim = 1}; - static w_window1d<double> - kernel_norm(double sigma, const point1d & mu = point1d(0)) - { - return internal::normalise(kernel(sigma, mu)); - } - - w_window1d<double> - static kernel(double sigma, const point1d & mu = point1d(0)) - { - assert(sigma >= 0); - w_window1d<double> w; - - const int size = int(sigma * SIGMA_COEF); - const double inv_sigma_sqrt_2pi = 1. / (sigma * sqrt(M_PI * 2.)); + //! Return a discrete Gaussian kernel, that is equal to zero + // at a distance of radius_factor * sigma. + static w_window1d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + return internal::normalise(kernel_values(sigma, radius_factor)); + } + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window1d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + w_window1d<ntg::float_d> w; + + const int size = int(sigma * radius_factor); + const ntg::float_d inv_sigma_sqrt_2pi + = 1. / (sigma * sqrt(M_PI * 2.)); for (int x = -size; x <= +size; ++x) - w.add(x, inv_sigma_sqrt_2pi * exp( -(x - mu.nth(0)) * (x - mu.nth(0)) / (2 * sigma * sigma)) - ); + w.add(x, inv_sigma_sqrt_2pi * exp(-x * x /(2 * sigma * sigma))); return w; } }; @@ -91,27 +101,42 @@ template<> struct gaussian_kernel<2> { - typedef w_window2d<double> ret; + typedef w_window2d<ntg::float_d> ret; enum {dim = 2}; - static w_window2d<double> - kernel_norm(double sigma, const point2d & mu = point2d(0, 0)) - { - return internal::normalise(kernel(sigma, mu)); - } - w_window2d<double> - static kernel(double sigma, const point2d & mu = point2d(0, 0)) - { - assert(sigma >= 0); - w_window2d<double> w; - const int size = int(sigma * SIGMA_COEF); - const double inv_sigma_sigma_pi_2 = 1. / (sigma * sigma * M_PI * 2.); + //! Return a discrete Gaussian kernel, that is equal to zero at a + // distance of radius_factor * sigma. + static w_window2d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + return internal::normalise(kernel_values(sigma, radius_factor)); + } + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window2d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + w_window2d<ntg::float_d> w; + + const int size = int(sigma * radius_factor); + const ntg::float_d inv_sigma_sigma_pi_2 = 1. / (sigma * sigma * + M_PI * 2.); for (int x = -size; x <= +size; ++x) for (int y = -size; y <= +size; ++y) + if (x * x + y * y <= size * size) w.add(x, y, inv_sigma_sigma_pi_2 - * exp(- ((x - mu.nth(0)) * (x - mu.nth(0)) + (y - mu.nth(1)) * (y - mu.nth(1)) - ) + * exp(- (x * x + y * y) / ( 2. * sigma * sigma) ) ); @@ -123,31 +148,47 @@ template<> struct gaussian_kernel<3> { - typedef w_window3d<double> ret; + typedef w_window3d<ntg::float_d> ret; enum {dim = 3}; - static w_window3d<double> - kernel_norm(double sigma, const point3d & mu = point3d(0, 0, 0)) + //! Return a discrete Gaussian kernel, that is equal to zero at a + // distance of radius_factor * sigma. + static w_window3d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) { - return internal::normalise(kernel(sigma, mu)); + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + return internal::normalise(kernel_values(sigma, radius_factor)); } - static w_window3d<double> - kernel(double sigma, const point3d & mu = point3d(0, 0, 0)) + + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window3d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) { - assert(sigma >= 0); - w_window3d<double> w; + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + w_window3d<ntg::float_d> w; - const int size = int(sigma * SIGMA_COEF); - const double k = 1. / (sigma * sigma * sigma * sqrt((M_PI * 2.) * (M_PI * 2.) * (M_PI * 2.))); + const int size = int(sigma * radius_factor); + const ntg::float_d k = 1. / (sigma * sigma * sigma * + sqrt((M_PI * 2.) * (M_PI * 2.) * + (M_PI * 2.))); for (int x = -size; x <= +size; ++x) for (int y = -size; y <= +size; ++y) for (int z = -size; z <= +size; ++z) + if (x * x + y * y + z * z <= size) w.add(x, y, z, k * - exp(-((x - mu.nth(0)) * (x - mu.nth(0)) + - (y - mu.nth(1)) * (y - mu.nth(1)) + - (z - mu.nth(2)) * (z - mu.nth(2)) - ) / (2. * sigma * sigma))); + exp(-(x * x + y * y + z * z) + ) / (2. * sigma * sigma)); return w; } };