Index: doc/ChangeLog
from Simon Odou <simon(a)lrde.epita.fr>
* ref/out/exdoc.config.in: Add -DHAVE_CONFIG_H and libs (for fft).
Index: olena/ChangeLog
from Simon Odou <simon(a)lrde.epita.fr>
* oln/convol/convolution.hh: Move to ...
* oln/convol/slow_convolution.hh: ... Here. Convolve two images
(convolution between image and w_window does not make sense).
* oln/transforms/fft.hh: Shift the DFT.
* tests/convol/Makefile.am: Add lib for fft test.
* tests/convol/tests/slow_sum_2d_1: Fix include.
* tests/convol/tests/slow_sum_2d_2: Likewise.
* tests/convol/tests/sobel_gradient: Likewise.
* oln/makefile.src: Add slow_convolution.hh and fast_convolution.hh.
* oln/convol/slow_gaussian.hh: Adapt to new kernel (image).
* oln/convol/slow_gaussian.hxx: Use image to compute convolution
(not w_window).
* oln/convol/fast_convolution.hh: New. Use fft to convolute.
* tests/convol/tests/fast_convol: New.
Index: olena/oln/convol/slow_convolution.hh
--- olena/oln/convol/convolution.hh Mon, 15 Mar 2004 17:40:54 +0100 van-vl_n (oln/f/39_convolutio 1.4.1.4.1.6.1.3 600)
+++ olena/oln/convol/slow_convolution.hh Mon, 14 Jun 2004 06:05:43 +0200 odou_s (oln/f/39_convolutio 1.4.1.4.1.6.1.3 600)
@@ -26,8 +26,8 @@
// Public License.
-#ifndef OLENA_CONVOL_CONVOLUTION_HH__
-# define OLENA_CONVOL_CONVOLUTION_HH__
+#ifndef OLENA_CONVOL_SLOW_CONVOLUTION_HH__
+# define OLENA_CONVOL_SLOW_CONVOLUTION_HH__
# include <oln/basics.hh>
# include <oln/basics2d.hh>
@@ -50,11 +50,13 @@
**
** \param DestValue Data type of the output image you want.
** \param I Exact type of the input image.
- ** \param Win Exact type of the window.
**
** \arg input The image to process.
** \arg win The window to convolve with.
**
+ ** \deprecated This function is obsolete, prefer :
+ ** convolve(const abstract::non_vectorial_image< I >& input,
+ ** const abstract::image< J >& k)
** \todo FIXME: we must always specify DestValue.
*/
template<class DestValue, class I, class Win>
@@ -80,6 +82,58 @@
}
/*!
+ ** \brief Perform a convolution of two images.
+ **
+ ** \param DestValue Data type of the output image you want.
+ ** \param I Exact type of the input image.
+ ** \param J Exact type of the input image (convolution kernel).
+ **
+ ** \arg input The image to process.
+ ** \arg k The kernel to convolve with.
+ **
+ ** \todo FIXME: we must always specify DestValue.
+ */
+ template<class DestValue, class I, class J>
+ typename mute<I, DestValue>::ret
+ convolve(const abstract::non_vectorial_image< I >& input,
+ const abstract::image< J >& k)
+ {
+ mlc::eq<I::dim, J::dim>::ensure();
+
+ typename mute<I, DestValue>::ret output(input.size());
+
+ // Compute Delta.
+ // \todo FIXME: should be in the image hierarchy.
+ coord delta = 0;
+ for (unsigned i = 0; i < J::dim; i++)
+ if (k.size().nth(i) > delta)
+ delta = k.size().nth(i);
+ input.border_adapt_copy(delta);
+
+ // Computer center of the kernel.
+ // \todo FIXME: should be in the image hierarchy.
+ oln_iter_type(I) p_im(input);
+ oln_iter_type(I) p_k(k);
+ oln_point_type(I) center;
+ unsigned i_center = 0;
+ unsigned real_center = k.npoints() % 2 ? k.npoints() / 2 + 1 :
+ k.npoints() / 2;
+ for_all(p_k)
+ if (++i_center == real_center)
+ center = p_k;
+
+ for_all(p_im)
+ {
+ DestValue sum = ntg_zero_val(DestValue);
+ for_all(p_k)
+ sum += static_cast<DestValue> (k[p_k]) *
+ static_cast<DestValue> (input[p_im - (center - p_k)]);
+ output[p_im] = sum;
+ }
+ return output;
+ }
+
+ /*!
** \brief Perform a convolution of an image and a window.
**
** \param DestValue Data type of the output image you want.
@@ -107,4 +161,4 @@
} // end namespace oln
-#endif // OLENA_CONVOL_CONVOLUTION_HH__
+#endif // OLENA_CONVOL_SLOW_CONVOLUTION_HH__
Index: olena/oln/transforms/fft.hh
--- olena/oln/transforms/fft.hh Fri, 19 Mar 2004 14:01:37 +0100 palma_g (oln/i/36_fft.hh 1.4.1.10 600)
+++ olena/oln/transforms/fft.hh Mon, 14 Jun 2004 05:38:16 +0200 odou_s (oln/i/36_fft.hh 1.4.1.10 600)
@@ -41,6 +41,7 @@
# include <ntg/basics.hh>
# include <ntg/all.hh>
# include <oln/basics2d.hh>
+# include <oln/morpher/piece_morpher.hh>
namespace oln {
@@ -485,6 +486,84 @@
}
/*!
+ ** \brief Shift zero-frequency component of discrete Fourier transform
+ ** to center of spectrum.
+ **
+ ** \param T1 The data type of the image returned.
+ **
+ ** The zero-frequency component of discrete Fourier transform are moved
+ ** to center of the image :
+ **
+ ** \htmlonly
+ ** <table>
+ ** <tr><td>1</td><td>2</td></tr>
+ ** <tr><td>3</td><td>4</td></tr>
+ ** </table>
+ ** becomes
+ ** <table>
+ ** <tr><td>4</td><td>3</td></tr>
+ ** <tr><td>2</td><td>1</td></tr>
+ ** </table>
+ ** \endhtmlonly
+ **
+ */
+ template <class T1>
+ image2d<T1> shift_transform_inv()
+ {
+ image2d<T1> t = transform_inv<T1>();
+ image2d<T1> st(t.size());
+
+ // We have to exchange t_1 with t_1_dest and not directly t_3 because
+ // they have not he same size.
+ typedef morpher::piece_morpher< image2d<T1> > piece_t;
+ piece_t t_1(t, dpoint2d(0, 0),
+ image2d_size((t.size().nrows() - 1) / 2,
+ (t.size().ncols() - 1) / 2,
+ t.border()));
+ piece_t t_1_dest(st, dpoint2d(t.nrows() - t_1.nrows(),
+ t.ncols() - t_1.ncols()),
+ image2d_size(t_1.nrows(), t_1.ncols(),
+ t.border()));
+ piece_t t_2(t, dpoint2d(0, (t.size().ncols() - 1) / 2),
+ image2d_size((t.size().nrows() - 1) / 2,
+ t.size().ncols() - (t.size().ncols() - 1) / 2,
+ t.border()));
+ piece_t t_2_dest(st, dpoint2d(t.nrows() - t_2.nrows(), 0),
+ image2d_size(t_2.nrows(), t_2.ncols(),
+ t.border()));
+ piece_t t_3(t, dpoint2d((t.size().nrows() - 1) / 2, 0),
+ image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2,
+ (t.size().ncols() - 1) / 2,
+ t.border()));
+ piece_t t_3_dest(st, dpoint2d(0, t.ncols() - t_3.ncols()),
+ image2d_size(t_3.nrows(), t_3.ncols(),
+ t.border()));
+ piece_t t_4(t, dpoint2d((t.size().nrows() - 1) / 2,
+ (t.size().ncols() - 1) / 2),
+ image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2,
+ t.size().ncols() - (t.size().ncols() - 1) / 2,
+ t.border()));
+ piece_t t_4_dest(st, dpoint2d(0, 0),
+ image2d_size(t_4.nrows(), t_4.ncols(),
+ t.border()));
+
+ oln_iter_type(piece_t) i1(t_1);
+ for_all(i1)
+ t_1_dest[i1] = t_1[i1];
+ oln_iter_type(piece_t) i2(t_2);
+ for_all(i2)
+ t_2_dest[i2] = t_2[i2];
+ oln_iter_type(piece_t) i3(t_3);
+ for_all(i3)
+ t_3_dest[i3] = t_3[i3];
+ oln_iter_type(piece_t) i4(t_4);
+ for_all(i4)
+ t_4_dest[i4] = t_4[i4];
+
+ return st;
+ }
+
+ /*!
** \brief Compute and return the invert transform.
*/
image2d<T> transform_inv()
@@ -492,6 +571,15 @@
return transform_inv<T>();
}
+ /*!
+ ** \brief Shift zero-frequency component of discrete Fourier transform
+ ** to center of spectrum.
+ */
+ image2d<T> shift_transform_inv()
+ {
+ return shift_transform_inv<T>();
+ }
+
};
/*!
Index: olena/tests/convol/Makefile.am
--- olena/tests/convol/Makefile.am Tue, 05 Aug 2003 21:39:33 +0200 burrus_n (oln/f/38_Makefile.a 1.13 600)
+++ olena/tests/convol/Makefile.am Mon, 14 Jun 2004 07:24:23 +0200 odou_s (oln/f/38_Makefile.a 1.13 600)
@@ -1,4 +1,7 @@
include ../check/Makefile.runtests
+LDADD_RUNTESTS += $(FFTW_LDFLAGS)
+CXXFLAGS_RUNTESTS += $(FFTW_CFLAGS)
+
EXTRA_DIST += sum_on_random.pgm
CLEANFILES += lena-sobel.pgm
Index: olena/tests/convol/tests/slow_sum_2d_1
--- olena/tests/convol/tests/slow_sum_2d_1 Tue, 20 Apr 2004 15:16:21 +0200 van-vl_n (oln/f/41_slow_sum_2 1.2.1.6 600)
+++ olena/tests/convol/tests/slow_sum_2d_1 Mon, 14 Jun 2004 08:14:27 +0200 odou_s (oln/f/41_slow_sum_2 1.2.1.6 600)
@@ -28,7 +28,7 @@
#include <oln/basics2d.hh>
-#include <oln/convol/convolution.hh>
+#include <oln/convol/slow_convolution.hh>
#include <oln/level/compare.hh>
Index: olena/tests/convol/tests/slow_sum_2d_2
--- olena/tests/convol/tests/slow_sum_2d_2 Tue, 20 Apr 2004 15:16:21 +0200 van-vl_n (oln/f/40_slow_sum_2 1.2.1.5 600)
+++ olena/tests/convol/tests/slow_sum_2d_2 Mon, 14 Jun 2004 08:14:37 +0200 odou_s (oln/f/40_slow_sum_2 1.2.1.5 600)
@@ -28,7 +28,7 @@
#include <oln/basics2d.hh>
-#include <oln/convol/convolution.hh>
+#include <oln/convol/slow_convolution.hh>
#include <oln/level/compare.hh>
Index: olena/tests/convol/tests/sobel_gradient
--- olena/tests/convol/tests/sobel_gradient Tue, 20 Apr 2004 15:16:21 +0200 van-vl_n (oln/f/45_sobel_grad 1.3.1.7 600)
+++ olena/tests/convol/tests/sobel_gradient Mon, 14 Jun 2004 08:15:02 +0200 odou_s (oln/f/45_sobel_grad 1.3.1.7 600)
@@ -28,7 +28,7 @@
#include <oln/basics2d.hh>
-#include <oln/convol/convolution.hh>
+#include <oln/convol/slow_convolution.hh>
#include <oln/convert/basics.hh>
#include "check.hh"
Index: doc/ref/out/exdoc.config.in
--- doc/ref/out/exdoc.config.in Fri, 09 Apr 2004 16:53:44 +0200 thivol_d (oln/k/7_exdoc.conf 1.7 600)
+++ doc/ref/out/exdoc.config.in Mon, 14 Jun 2004 07:03:20 +0200 odou_s (oln/k/7_exdoc.conf 1.7 600)
@@ -4,7 +4,8 @@
TAG_CLOSE = endcode
CAPTIONS = cxx # We want to run cxx on the extracted files (see line below)
ALIAS cxx = @CXX@ # Here, cxx means g++ but you can choose other compilers
- OPTIONS = @DOC_CPPFLAGS@ @CXXFLAGS_OPTIMIZE@ @CXXFLAGS_STRICT_ERRORS@ -I@top_srcdir@/integre -I@top_builddir@/olena -I@top_srcdir@/olena -I@top_srcdir@/metalic $1 -o $2 -DIMG_OUT=\"../img/\" -DIMG_IN=\"@top_srcdir@/olena/img/\" # tell how to use the soft, i.e. where to put input and output arguments (default if not overriden) ($1: input, $2: output) FIXME: $* should have explicit name, chek flags
+ # FIXME: we should write the compilation line in the source file (for libs).
+ OPTIONS = -lz -lfftw -lrfftw @DOC_CPPFLAGS@ @CXXFLAGS_OPTIMIZE@ @CXXFLAGS_STRICT_ERRORS@ -I@top_srcdir@/integre -I@top_builddir@/olena -I@top_srcdir@/olena -I@top_srcdir@/metalic -I@top_builddir@ $1 -o $2 -DIMG_OUT=\"../img/\" -DIMG_IN=\"@top_srcdir@/olena/img/\" -DHAVE_CONFIG_H # tell how to use the soft, i.e. where to put input and output arguments (default if not overriden) ($1: input, $2: output) FIXME: $* should have explicit name, chek flags
OUT = out # FIXME: should be obsolete
EXT = cc # Extension of generated file
STD_OUT_EXT = std # Extension of generated file standard output
Index: olena/oln/makefile.src
--- olena/oln/makefile.src Sun, 06 Jun 2004 23:17:15 +0200 thivol_d (oln/r/4_makefile.s 1.2 600)
+++ olena/oln/makefile.src Fri, 11 Jun 2004 03:14:27 +0200 odou_s (oln/r/4_makefile.s 1.2 600)
@@ -31,7 +31,8 @@
convert/rgbxyz.hh \
convert/stretch.hh \
convert/value_to_point.hh \
- convol/convolution.hh \
+ convol/slow_convolution.hh \
+ convol/fast_convolution.hh \
convol/fast_gaussian.hh \
convol/fast_gaussian.hxx \
convol/fast_gaussian_coefficient.hh \
Index: olena/oln/convol/slow_gaussian.hh
--- olena/oln/convol/slow_gaussian.hh Wed, 09 Jun 2004 20:13:57 +0200 van-vl_n (oln/r/7_slow_gauss 1.2 644)
+++ olena/oln/convol/slow_gaussian.hh Mon, 14 Jun 2004 08:31:10 +0200 odou_s (oln/r/7_slow_gauss 1.2 644)
@@ -30,7 +30,7 @@
# define SLOW_GAUSSIAN_HH
# include <oln/core/image.hh>
# include <ntg/float.hh>
-# include "convolution.hh"
+# include <oln/convol/slow_convolution.hh>
# include <oln/core/behavior.hh>
namespace oln {
@@ -87,7 +87,14 @@
{
const typename gaussian_kernel<I::dim>::ret k
= gaussian_kernel<I::dim>::kernel(sigma, radius_factor);
- behavior.adapt_border(in, coord(k.delta()));
+
+ // Compute Delta.
+ // \todo FIXME: should be in the image hierarchy.
+ coord delta = 0;
+ for (unsigned i = 0; i < I::dim; i++)
+ if (in.size().nth(i) > delta)
+ delta = in.size().nth(i);
+ behavior.adapt_border(in, delta);
typename mute<I, ntg::float_d>::ret im =
oln::convol::slow::convolve<ntg::float_d> (in, k);
Index: olena/oln/convol/slow_gaussian.hxx
--- olena/oln/convol/slow_gaussian.hxx Wed, 09 Jun 2004 20:13:57 +0200 van-vl_n (oln/r/8_slow_gauss 1.2 644)
+++ olena/oln/convol/slow_gaussian.hxx Mon, 14 Jun 2004 05:42:37 +0200 odou_s (oln/r/8_slow_gauss 1.2 644)
@@ -28,9 +28,9 @@
#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>
+# include <oln/core/image1d.hh>
+# include <oln/core/image2d.hh>
+# include <oln/core/image3d.hh>
namespace oln {
namespace convol {
@@ -44,29 +44,31 @@
inline T
normalise(const T &in)
{
- T w = in;
+ T w(in.size());
ntg::float_d sum = 0.;
- for (unsigned i = 0; i < in.card(); ++i)
- sum += in.w(i);
+ oln_iter_type(T) i(in);
+ for_all(i)
+ sum += in[i];
sum = 1. / sum;
assert(finite(sum));
- for (unsigned i = 0; i < in.card(); ++i)
- w.set(w.dp(i), in.w(i) * sum);
+ for_all(i)
+ w[i] = in[i] * sum;
return w;
}
+
}// End namespace internal
//! Kernel of a 1D Gaussian.
template<>
struct gaussian_kernel<1>
{
- typedef oln::w_window1d<ntg::float_d> ret;
+ typedef oln::image1d<ntg::float_d> ret;
enum {dim = 1};
//! Return a discrete Gaussian kernel, that is equal to zero
// at a distance of radius_factor * sigma.
- static w_window1d<ntg::float_d>
+ static image1d<ntg::float_d>
kernel(ntg::float_d sigma,
ntg::float_d radius_factor)
{
@@ -79,20 +81,20 @@
//
// 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>
+ static image1d<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);
+ image1d<ntg::float_d> w(size * 2 + 1);
+
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 * x /(2 * sigma * sigma)));
+ for (int x = -size; x <= size; ++x)
+ w(x + size) = inv_sigma_sqrt_2pi * exp(-x * x /(2 * sigma * sigma));
return w;
}
};
@@ -101,12 +103,12 @@
template<>
struct gaussian_kernel<2>
{
- typedef w_window2d<ntg::float_d> ret;
+ typedef image2d<ntg::float_d> ret;
enum {dim = 2};
//! Return a discrete Gaussian kernel, that is equal to zero at a
// distance of radius_factor * sigma.
- static w_window2d<ntg::float_d>
+ static image2d<ntg::float_d>
kernel(ntg::float_d sigma,
ntg::float_d radius_factor)
{
@@ -120,26 +122,27 @@
//
// 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>
+ static image2d<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);
+ image2d<ntg::float_d> w = image2d<ntg::float_d>(size * 2 + 1,
+ size * 2 + 1);
+
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)
+ 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
+ w(x + size, y + size) = inv_sigma_sigma_pi_2
* exp(- (x * x + y * y)
/ ( 2. * sigma * sigma)
- )
);
+ else
+ w(x + size, y + size) = 0;
return w;
}
};
@@ -148,12 +151,12 @@
template<>
struct gaussian_kernel<3>
{
- typedef w_window3d<ntg::float_d> ret;
+ typedef image3d<ntg::float_d> ret;
enum {dim = 3};
//! Return a discrete Gaussian kernel, that is equal to zero at a
// distance of radius_factor * sigma.
- static w_window3d<ntg::float_d>
+ static image3d<ntg::float_d>
kernel(ntg::float_d sigma,
ntg::float_d radius_factor)
{
@@ -168,16 +171,18 @@
//
// 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>
+ static image3d<ntg::float_d>
kernel_values(ntg::float_d sigma,
ntg::float_d radius_factor)
{
precondition(sigma > 0.);
precondition(radius_factor >= 0.);
- w_window3d<ntg::float_d> w;
-
const int size = int(sigma * radius_factor);
+ image3d<ntg::float_d> w(size * 2 + 1,
+ size * 2 + 1,
+ size * 2 + 1);
+
const ntg::float_d k = 1. / (sigma * sigma * sigma *
sqrt((M_PI * 2.) * (M_PI * 2.) *
(M_PI * 2.)));
@@ -186,9 +191,11 @@
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 *
+ w(x + size, y + size, z + size) = k *
exp(-(x * x + y * y + z * z)
- ) / (2. * sigma * sigma));
+ ) / (2. * sigma * sigma);
+ else
+ w(x + size, y + size, z + size) = 0;
return w;
}
};
Index: olena/oln/convol/fast_convolution.hh
--- olena/oln/convol/fast_convolution.hh Mon, 14 Jun 2004 09:04:45 +0200 odou_s ()
+++ olena/oln/convol/fast_convolution.hh Mon, 14 Jun 2004 07:40:15 +0200 odou_s (oln/r/12_fast_convo 644)
@@ -0,0 +1,173 @@
+// 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 OLENA_CONVOL_FAST_CONVOLUTION_HH__
+# define OLENA_CONVOL_FAST_CONVOLUTION_HH__
+
+# include <oln/basics.hh>
+# include <oln/basics2d.hh>
+# include <ntg/all.hh>
+# include <mlc/cmp.hh>
+# include <mlc/array/all.hh>
+# include <oln/transforms/fft.hh>
+# include <oln/morpher/piece_morpher.hh>
+
+namespace oln {
+
+ /*! \brief Algorithms related to convolution.
+ */
+ namespace convol {
+
+ /*! \brief Convolution algorithms.
+ */
+ namespace fast {
+
+ /*!
+ ** \brief Perform a convolution of two images.
+ **
+ ** \param DestValue Data type of the output image you want.
+ ** \param I Exact type of the input image.
+ ** \param J Exact type of the input image (convolution kernel).
+ **
+ ** \arg input The image to process.
+ ** \arg k The kernel to convolve with.
+ **
+ ** \todo FIXME: the kernel could be bigger than the image.
+ ** \todo FIXME: supports only 2 dimensions (like the fft).
+ **
+ ** The algorithm is just the convolution theorem application.
+ **
+ ** This example compute a fast gaussian convolution :
+ ** \code
+ ** #include <oln/basics2d.hh>
+ ** #include <oln/convol/slow_gaussian.hh>
+ ** #include <oln/convol/fast_convolution.hh>
+ **
+ ** using namespace oln;
+ ** using namespace ntg;
+ **
+ ** int main() {
+ ** image2d<int_u8> src(IMG_IN "lena.pgm");
+ ** float_d sigma = 2.5;
+ ** float_d radius = 3;
+ ** image2d<float_d> k = convol::slow::gaussian_kernel<2>::kernel(sigma, radius);
+ ** image2d<float_d> f_src(src.size());
+ ** oln_iter_type_(image2d<float_d>) i(src);
+ ** for_all(i)
+ ** f_src[i] = float_d(src[i]);
+ ** image2d<float_d> tmp = convol::fast::convolve<float_d>(f_src, k);
+ ** for_all(i)
+ ** src[i] = int_u8(tmp[i]);
+ ** save(src, IMG_OUT "oln_convol_fast_convolve.pgm");
+ ** }
+ ** \endcode
+ ** \image html lena_pgm.png
+ ** \image latex lena_pgm.png
+ ** =>
+ ** \image html oln_convol_fast_convolve.pgm
+ ** \image latex oln_convol_fast_convole.pgm
+ **
+ */
+ template<class DestValue, class I, class J>
+ typename mute<I, DestValue>::ret
+ convolve(const abstract::image< I >& input,
+ const abstract::image< J >& k)
+ {
+ mlc::eq<I::dim, J::dim>::ensure();
+ mlc::eq<I::dim, 2>::ensure();
+ assert(input.npoints() > k.npoints());
+
+ // We compute k with a size of input (k is centered in big_k).
+ image2d<oln_value_type(J)> big_k(input.size());
+#define CENTER_DST(I) \
+((I) % 2 ? (I) / 2 + 1 : (I) / 2)
+#define CENTER_SRC(I) \
+((I) % 2 ? (I) / 2 : ((I) - 1) / 2)
+ typedef morpher::piece_morpher< image2d<oln_value_type(J)> > piece_t;
+ piece_t piece_k(big_k,
+ dpoint2d(CENTER_DST(big_k.size().nrows()) -
+ CENTER_SRC(k.size().nrows()) - 1,
+ CENTER_DST(big_k.size().ncols()) -
+ CENTER_SRC(k.size().ncols()) - 1),
+ oln::image2d_size(k.size().nrows(),
+ k.size().ncols(),
+ big_k.border()));
+ oln_iter_type(piece_t) i_k(piece_k);
+ for_all(i_k)
+ piece_k[i_k] = k[i_k];
+
+ transforms::fft<oln_value_type(I), ntg::rect> tr_input(input.exact());
+ transforms::fft<oln_value_type(J), ntg::rect> tr_k(big_k.exact());
+
+ tr_input.transform();
+ typename mute<I, ntg::cplx<ntg::rect, ntg::float_d> >::ret&
+ Input = tr_input.transformed_image();
+ const typename mute<J, ntg::cplx<ntg::rect, ntg::float_d> >::ret
+ K = tr_k.transform();
+
+ oln_iter_type(I) i_input(Input);
+ for_all(i_input) {
+ Input[i_input] *= K[i_input];
+ // Scale.
+ Input[i_input] *= Input.size().nrows() * Input.size().ncols();
+ }
+
+ typename mute<I, DestValue>::ret output = tr_input.shift_transform_inv();
+
+ return output;
+ }
+
+ /*!
+ ** \brief Perform a convolution of an image and a window.
+ **
+ ** \param DestValue Data type of the output image you want.
+ ** \param I Exact type of the input image.
+ ** \param Info Informations about the array.
+ ** \param Win Data type of the array.
+ **
+ ** \arg input The image to process.
+ ** \arg arr The array to convolve with.
+ **
+ ** \todo FIXME: don't use array1d, ..., arraynd.
+ */
+ template<class DestValue, class I, class Info, class Win>
+ typename mute<I, DestValue>::ret
+ convolve(const abstract::image < I >& input,
+ const mlc::array2d<Info, Win >& arr)
+ {
+ return convolve<DestValue>(input, static_cast< image2d<Win> >(arr));
+ // FIXME: Should be abstract::w_window<T_arr>. Adjust #include once done.
+ }
+
+ } // end namespace fast
+
+ } // end namespace convol
+
+} // end namespace oln
+
+#endif // OLENA_CONVOL_FAST_CONVOLUTION_HH__
Index: olena/tests/convol/tests/fast_convol
--- olena/tests/convol/tests/fast_convol Mon, 14 Jun 2004 09:04:45 +0200 odou_s ()
+++ olena/tests/convol/tests/fast_convol Mon, 14 Jun 2004 07:59:18 +0200 odou_s (oln/r/13_fast_convo 644)
@@ -0,0 +1,33 @@
+#include "data.hh"
+#include "check.hh"
+
+#include <oln/basics2d.hh>
+#include <oln/convol/slow_gaussian.hh>
+#include <oln/convol/fast_convolution.hh>
+#include <oln/utils/md5.hh>
+
+using namespace oln;
+using namespace ntg;
+
+int main() {
+ oln::utils::key::value_type data_key[] =
+ {0x7a, 0x42, 0xc0, 0xa5, 0xec, 0xc7, 0x14, 0xa, 0xc7,
+ 0x21, 0xd9, 0x1f, 0xdb, 0xfa, 0x54, 0xb9};
+ oln::utils::key key(data_key);
+
+ image2d<int_u8> src(rdata("lena.pgm"));
+ float_d sigma = 2.5;
+ float_d radius = 3;
+ image2d<float_d> k = convol::slow::gaussian_kernel<2>::kernel(sigma, radius);
+ image2d<float_d> f_src(src.size());
+ oln_iter_type_(image2d<float_d>) i(src);
+ for_all(i)
+ f_src[i] = float_d(src[i]);
+
+ image2d<float_d> tmp = convol::fast::convolve<float_d>(f_src, k);
+
+ for_all(i)
+ src[i] = tmp[i];
+ if (oln::utils::md5(src) != key)
+ return false;
+}
--
Simon Odou
simon(a)lrde.epita.fr
