[Olena-patches] oln 10.183: Slow Gaussian and neighborhood3d

Index: olena/ChangeLog from Niels Van Vliet &lt;niels(a)lrde.epita.fr&gt; * olena/oln/core/neighborhood3d.hh: Add some common neighborhoods. * olena/oln/convol/slow_gaussian.hh: Add file. * olena/oln/convol/slow_gaussian.hxx: Add file. Index: olena/oln/core/neighborhood3d.hh --- olena/oln/core/neighborhood3d.hh Mon, 29 Mar 2004 09:53:11 +0200 odou_s (oln/c/35_neighborho 1.19 600) +++ olena/oln/core/neighborhood3d.hh Tue, 11 May 2004 09:50:46 +0200 van-vl_n (oln/c/35_neighborho 1.19 600) @@ -260,6 +260,55 @@ return neighb; } + + /*! + ** \brief Create an ellipsoid neighborhood (3 dimension). + ** \arg zradius radius Z. + ** \arg yradius radius Y. + ** \arg xradius radius X. + ** \pre zradius > 0 + ** \pre yradius > 0 + ** \pre xradius > 0 + ** + ** The ellipsoid formula is : + ** \f$$\frac{x^2}{xradius^2} + \frac{y^2}{yradius^2} + \frac{z^2}{zradius^2} = 1$\f$ + */ + inline neighborhood3d + mk_neighb_ellipsoid(float zradius, float yradius, float xradius) + { + precondition(zradius > 0); + precondition(yradius > 0); + precondition(xradius > 0); + + neighborhood3d neighb; + coord zmax = (coord)roundf(zradius); + float zr2 = zradius * zradius; + float yr2 = yradius * yradius; + for (coord z = -zmax; z <= zmax; ++z) + { + /* + x^2 y^2 z^2 + --------- + --------- + --------- = 1 + xradius^2 yradius^2 zradius^2 + */ + + /* Set x to 0 in the above formula to find ymax. */ + float v = 1 - z * z / zr2; + if (v < 0) v = 0; // Can happen because zmax has been rounded. + coord ymax = (coord)roundf(yradius * sqrtf(v)); + for (coord y = -ymax; y <= ymax; ++y) + { + float w = v - y * y / yr2; + if (w < 0) w = 0; // Can happen because ymax has been rounded. + coord xmax = (coord)roundf(xradius * sqrtf(w)); + for (coord x = z > 0 ? 0 : 1; x <= xmax; ++x) + neighb.add(z, y, x); + } + } + return neighb; + } + + /*! ** \brief Create a cube neighborhood (3 dimension). ** \arg width Number of slice, colunm and row. @@ -285,6 +334,16 @@ return win; } + /*! + ** \brief Create a ball neighborhood (3 dimension). + ** \arg radius The radius. + ** \return The new neighborhood (3d). + */ + inline neighborhood3d + mk_neighb_ball(float radius) + { + return mk_neighb_ellipsoid(radius, radius, radius); + } } // end of oln #endif // OLENA_CORE_NEIGHBORHOOD3D_HH Index: olena/oln/convol/slow_gaussian.hh --- olena/oln/convol/slow_gaussian.hh Wed, 12 May 2004 23:10:01 +0200 van-vl_n () +++ olena/oln/convol/slow_gaussian.hh Wed, 12 May 2004 23:04:19 +0200 van-vl_n (oln/r/7_slow_gauss 644) @@ -0,0 +1,80 @@ +// Copyright (C) 2004 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 59 Temple Place - Suite 330, Boston, +// MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + + +#ifndef SLOW_GAUSSIAN_HH +# define SLOW_GAUSSIAN_HH +# include <oln/core/image.hh> +# include <ntg/float.hh> +# include "convolution.hh" + +namespace oln { + namespace convol { + namespace slow { + + /*! Gaussian kernel. + ** + ** In the specialization, the function kernel must be defined. + */ + template <int Dim> + struct gaussian_kernel + { + }; + + /*! 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. + */ + template <class I> + 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)()) + { + 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)); + 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; + } + }//End namespace slow + }// End namespace convol +}// End namespace oln + +#include "slow_gaussian.hxx" +#endif // end SLOW_GAUSSIAN_HH Index: olena/oln/convol/slow_gaussian.hxx --- olena/oln/convol/slow_gaussian.hxx Wed, 12 May 2004 23:10:01 +0200 van-vl_n () +++ olena/oln/convol/slow_gaussian.hxx Wed, 12 May 2004 23:07:36 +0200 van-vl_n (oln/r/8_slow_gauss 644) @@ -0,0 +1,158 @@ +// Copyright (C) 2004 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 59 Temple Place - Suite 330, Boston, +// MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + + +#ifndef SLOW_GAUSSIAN_HXX +# define SLOW_GAUSSIAN_HXX +# include <oln/core/w_window1d.hh> +# 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 { + namespace internal { + /* FIXME: The hierarchy w_window is a REAL mess. + a w_window1d is not an abstract::w_window, + and w_windownd does not supprot concrete_type, etc + */ + template<class T> + T + normalise(const T &in) + { + T w = in; + double sum = 0; + for (unsigned i = 0; i < in.card(); ++i) + sum += in.w(i); + sum = 1. / sum; + for (unsigned i = 0; i < in.card(); ++i) + w.set(w.dp(i), in.w(i) * sum); + return w; + } + + }// End namespace internal + + //! Kernel of a 1D Gaussian. + template<> + struct gaussian_kernel<1> + { + typedef oln::w_window1d<double> 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.)); + 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)) + ); + return w; + } + }; + + //! Kernel of a 2D Gaussian. + template<> + struct gaussian_kernel<2> + { + typedef w_window2d<double> 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.); + + for (int x = -size; x <= +size; ++x) + for (int y = -size; y <= +size; ++y) + 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)) + ) + / ( 2. * sigma * sigma) + ) + ); + return w; + } + }; + + //! Kernel of a 3D Gaussian. + template<> + struct gaussian_kernel<3> + { + typedef w_window3d<double> ret; + enum {dim = 3}; + + static w_window3d<double> + kernel_norm(double sigma, const point3d & mu = point3d(0, 0, 0)) + { + return internal::normalise(kernel(sigma, mu)); + } + static w_window3d<double> + kernel(double sigma, const point3d & mu = point3d(0, 0, 0)) + { + assert(sigma >= 0); + w_window3d<double> 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.))); + + for (int x = -size; x <= +size; ++x) + for (int y = -size; y <= +size; ++y) + for (int z = -size; z <= +size; ++z) + 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))); + return w; + } + }; + + }//End namespace slow + }// End namespace convol +}// End namespace oln +#endif // end SLOW_GAUSSIAN_HXX