[Olena] Re: oln 10.192: convolution

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; +}