Index: olena/ChangeLog
from Niels Van Vliet <niels(a)lrde.epita.fr>
* 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
