Olena-patches

olena-patches@lrde.epita.fr

9625 discussions

olena: olena-2.0-632-g1e9b9bd Update copyright headers.

by Roland Levillain

* mln/transform/fft.hh, * tests/transform/fft.cc: Here. --- milena/ChangeLog | 8 ++++++++ milena/mln/transform/fft.hh | 28 +++++++++++++--------------- milena/tests/transform/fft.cc | 29 +++++++++++++---------------- 3 files changed, 34 insertions(+), 31 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 681452a..eb7148d 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,13 @@ 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + Update copyright headers. + + * mln/transform/fft.hh, + * tests/transform/fft.cc: + Here. + +2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + * headers.mk, tests/unit_test/unit-tests.mk: Regen. 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index cf0f1bb..a010e47 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -1,30 +1,28 @@ // Copyright (C) 2007, 2008, 2009, 2012 EPITA Research and Development // Laboratory // -// This file is part of the Milena 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 file is part of Olena. // -// This library is distributed in the hope that it will be useful, +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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, 51 Franklin Street, Fifth Floor, -// Boston, MA 02111-1307, USA. +// along with Olena. If not, see <http://www.gnu.org/licenses/>. // // As a special exception, you may use this file as part of a free -// software library without restriction. Specifically, if other files +// software project 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. +// 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 MLN_TRANSFORM_FFT_HH # define MLN_TRANSFORM_FFT_HH diff --git a/milena/tests/transform/fft.cc b/milena/tests/transform/fft.cc index bf7b8b9..c5c5bec 100644 --- a/milena/tests/transform/fft.cc +++ b/milena/tests/transform/fft.cc @@ -1,30 +1,27 @@ -// -*- c++ -*- // Copyright (C) 2004, 2012 EPITA Research and Development Laboratory // -// This file is part of the Milena 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 file is part of Olena. // -// This library is distributed in the hope that it will be useful, +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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, Inc., 51 Franklin Street, Fifth Floor, -// Boston, MA 02110-1301, USA. +// along with Olena. If not, see <http://www.gnu.org/licenses/>. // // As a special exception, you may use this file as part of a free -// software library without restriction. Specifically, if other files +// software project 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. +// 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. #include <mln/transform/fft.hh> #include <mln/value/int_u8.hh> -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-633-g870ac7a Improve the test of the Fast Fourier Transform.

by Roland Levillain

* tests/transform/fft.cc: Ensure the inverse transform of the transform is pretty much the same as the input image. Save the cropped transform. Remove dead code. * tests/transform/Makefile.am (MOSTLYCLEANFILES): Add fft_trans_cropped.pgm. --- milena/ChangeLog | 11 +++++++++++ milena/tests/transform/Makefile.am | 1 + milena/tests/transform/fft.cc | 18 ++++++++++++++++-- 3 files changed, 28 insertions(+), 2 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index eb7148d..9cfd3a6 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,16 @@ 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + Improve the test of the Fast Fourier Transform. + + * tests/transform/fft.cc: Ensure the inverse transform of the + transform is pretty much the same as the input image. + Save the cropped transform. + Remove dead code. + * tests/transform/Makefile.am (MOSTLYCLEANFILES): + Add fft_trans_cropped.pgm. + +2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + Update copyright headers. * mln/transform/fft.hh, diff --git a/milena/tests/transform/Makefile.am b/milena/tests/transform/Makefile.am index b266926..a49ae65 100644 --- a/milena/tests/transform/Makefile.am +++ b/milena/tests/transform/Makefile.am @@ -46,4 +46,5 @@ MOSTLYCLEANFILES = \ fft_copy.pgm \ fft_low_pass.pgm \ fft_trans_clipped.pgm \ + fft_trans_cropped.pgm \ fft_trans_log.pgm diff --git a/milena/tests/transform/fft.cc b/milena/tests/transform/fft.cc index c5c5bec..1dce9df 100644 --- a/milena/tests/transform/fft.cc +++ b/milena/tests/transform/fft.cc @@ -30,6 +30,12 @@ #include <mln/opt/at.hh> #include <mln/debug/println.hh> +#include <mln/core/image/flat_image.hh> +#include <mln/fun/vv2b/le.hh> +#include <mln/fun/vv2v/diff_abs.hh> +#include <mln/data/transform.hh> +#include <mln/test/predicate.hh> + #include "tests/data.hh" #define CHECK(Condition) \ @@ -61,11 +67,16 @@ int main () io::pgm::save(im3, "fft_copy.pgm"); image2d<int_u8> fft = fourier.transformed_image_log_magn<int_u8>(true); - // debug::println(fourier.transformed_image()); io::pgm::save(fft, "fft.pgm"); std::cout << "Test: Image == F-1(F(Image)) ... " << std::flush; - // CHECK (im1 == im3); + /* FIXME: Milena lacks some feature to make write the following in a + shorter fashion (fun-morpher accepting binary (vv2v) functions or + function composition in pw-functions. */ + CHECK(test::predicate(data::transform(im1, im3, + fun::vv2v::diff_abs<int_u8>()), + flat_image<int_u8, box2d>(1, im1.domain()), + fun::vv2b::le<int_u8>())); image2d<int_u8> out = fourier.transformed_image_clipped_magn<int_u8>(0.01); @@ -83,6 +94,9 @@ int main () for (int col = 40; col < im2.ncols() - 40; ++col) opt::at(im2, row, col) = 0; + fft = fourier.transformed_image_log_magn<int_u8>(true); + io::pgm::save(fft, "fft_trans_cropped.pgm"); + out = fourier.transform_inv<int_u8>(); io::pgm::save(out, "fft_low_pass.pgm"); -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-634-g6919bfe Minor corrections in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh (internal::_fft<T>::transformed_image() const): Return a const reference instead of a copy. (internal::_fft<T>::transformed_image_clipped_magn(double, bool)): Drop the const before the first argument. Properly initialize `max'. (internal::_fft<T>::transformed_image_log_magn(double, double, bool)): Properly initialize `max'. --- milena/ChangeLog | 13 +++++++++++++ milena/mln/transform/fft.hh | 10 +++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 9cfd3a6..e163188 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,16 @@ +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + + Minor corrections in the Fast Fourier Transform. + + * mln/transform/fft.hh + (internal::_fft<T>::transformed_image() const): Return a const + reference instead of a copy. + (internal::_fft<T>::transformed_image_clipped_magn(double, bool)): + Drop the const before the first argument. + Properly initialize `max'. + (internal::_fft<T>::transformed_image_log_magn(double, double, bool)): + Properly initialize `max'. + 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> Improve the test of the Fast Fourier Transform. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index a010e47..ed1e87f 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -81,7 +81,7 @@ namespace mln { ** ** Const version. */ - const image2d< std::complex<T> > transformed_image() const + const image2d< std::complex<T> >& transformed_image() const { return trans_im; } @@ -156,14 +156,14 @@ namespace mln { */ template <class R> - image2d<R> transformed_image_clipped_magn(const double clip, + image2d<R> transformed_image_clipped_magn(double clip, bool ordered = true) const { // Check that R is real image2d<R> new_im(trans_im.domain()); // check that clip is >=0 and <=1 ? - double max; + double max = mln_min(double); mln_piter(image2d<T>) it(trans_im.domain()); for_all(it) @@ -209,7 +209,7 @@ namespace mln { ** \param clip Value used for clipping. ** \param ordered Kind of traversal. */ - image2d<T> transformed_image_clipped_magn(const double clip, + image2d<T> transformed_image_clipped_magn(double clip, bool ordered = true) const { return transformed_image_clipped_magn<T>(clip, ordered); @@ -271,7 +271,7 @@ namespace mln { image2d<R> new_im(trans_im.domain()); - double max = 0; + double max = mln_min(double); mln_piter(image2d<R>) it(trans_im.domain()); for_all(it) -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-635-g089d5ff Remove dead code in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 +++ milena/mln/transform/fft.hh | 93 ------------------------------------------- 2 files changed, 6 insertions(+), 93 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index e163188..6a5f5d7 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Remove dead code in the Fast Fourier Transform. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Minor corrections in the Fast Fourier Transform. * mln/transform/fft.hh diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index ed1e87f..dfbe1d7 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -479,88 +479,6 @@ namespace mln { } /*! - ** \brief Shift zero-frequency component of discrete Fourier transform - ** to center of spectrum. - ** - ** \param R 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 R> - image2d<R> shift_transform_inv() - { - image2d<R> t = transform_inv<R>(); - image2d<R> 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<R> > 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() @@ -568,17 +486,6 @@ namespace mln { 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>(); - } - */ - }; /*! -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-636-g588a6d4 Stylistic changes in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 + milena/mln/transform/fft.hh | 630 ++++++++++++++++++++----------------------- 2 files changed, 305 insertions(+), 331 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 6a5f5d7..8f96157 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Stylistic changes in the Fast Fourier Transform. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Remove dead code in the Fast Fourier Transform. * mln/transform/fft.hh: Here. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index dfbe1d7..9b70b4f 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -27,17 +27,20 @@ #ifndef MLN_TRANSFORM_FFT_HH # define MLN_TRANSFORM_FFT_HH -# include <mln/core/image/image2d.hh> -# include <mln/estim/min_max.hh> -# include <mln/opt/at.hh> +# include <complex> -# include <complex> +# include <fftw3.h> -# include <fftw3.h> +# include <mln/core/image/image2d.hh> +# include <mln/estim/min_max.hh> +# include <mln/opt/at.hh> -namespace mln { - namespace internal { +namespace mln +{ + + namespace internal + { /// Dispatch traits for fftw enum fft_dispatch { fft_cplx, fft_real }; @@ -45,316 +48,279 @@ namespace mln { template <typename T> struct fft_trait; - /*! - ** \brief fft trait. - */ + /// \brief FFT traits for `double'. template <> struct fft_trait<double> { - static const fft_dispatch which = fft_real; ///< Real dispatch. - typedef double fftw_input; ///< Type of input. + /// Type of dispatch (real numbers). + static const fft_dispatch which = fft_real; + /// Type of input (double-precision floating point numbers). + typedef double fftw_input; }; - /*! - ** \brief fft trait. - ** - ** \specialization for ntg::cplx<R, T> - */ + /// \brief FFT traits for `std::complex<T>'. template <typename T> struct fft_trait< std::complex<T> > { - static const fft_dispatch which = fft_cplx; ///< Complex dispatch. - typedef std::complex <T> fftw_input; ///< Type of input. + /// Type of dispatch (complex numbers). + static const fft_dispatch which = fft_cplx; + /// Type of input (complex numbers of the standard C++ library). + typedef std::complex<T> fftw_input; }; - /*! - ** _fft<ntg::cplx_representation, T> - ** - ** \param T Data type. - */ + /// \brief Internal structure containining packing data and + /// instructions for the (forward) and inverse (backward) + /// transforms. + /// + /// \tparam T Data type. template <class T> - class _fft + class fft { public: - /*! - ** \brief Accessor to transformed image. - ** - ** Const version. - */ + /// \brief Accessor to transformed image (const version). const image2d< std::complex<T> >& transformed_image() const { return trans_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** Non const version. - */ + /// \brief Accessor to transformed image (non const version). image2d< std::complex<T> >& transformed_image() { return trans_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** |T[p]|. - ** - ** \tparam R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (magnitude). + /// + /// For each point \p p of the transformed image \p T, + /// an image containing <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_magn(bool ordered = true) const { - // FIXME : check that R is real + // FIXME: Ensure R is real. + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); image2d<R> new_im(trans_im.domain()); - if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) opt::at(new_im, row, col) = - std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % trans_im.ncols())); + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)); else { mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) new_im(it) = std::norm(trans_im(it)); } - return new_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** |T[p]|. - ** - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (magnitude). + /// + /// For each point \p p of the transformed image \p T, + /// an image containing <tt>|T[p]|</tt> is returned. + /// + /// \param ordered Kind of traversal. image2d<T> transformed_image_magn(bool ordered = true) const { return transformed_image_magn<T>(ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param clip Value used for clipping. - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param clip Value used for clipping. + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(double clip, bool ordered = true) const { - // Check that R is real + // FIXME: Ensure R is real. + // FIXME: Ensure `clip' is between 0 and 1? - image2d<R> new_im(trans_im.domain()); - // check that clip is >=0 and <=1 ? double max = mln_min(double); mln_piter(image2d<T>) it(trans_im.domain()); - for_all(it) if (std::norm(trans_im(it)) > max) max = std::norm(trans_im(it)); + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) - { - if (std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % - trans_im.ncols())) >= max * clip) - opt::at(new_im, row, col) = mln_max(R); - else - opt::at(new_im, row, col) = - (double) mln_max(R) * - std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % - trans_im.ncols())) / (max * clip); - } - else - { - for_all(it) - { - if (std::norm(trans_im(it)) >= max * clip) - new_im(it) = mln_max(R); + if (std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) + >= max * clip) + opt::at(new_im, row, col) = mln_max(R); else - new_im(it) = (double) mln_max(R) * - std::norm(trans_im(it)) / (max * clip); - } - } - + opt::at(new_im, row, col) = + (double) mln_max(R) * + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) / (max * clip); + else + for_all(it) + if (std::norm(trans_im(it)) >= max * clip) + new_im(it) = mln_max(R); + else + new_im(it) = + (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); return new_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param clip Value used for clipping. - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \param clip Value used for clipping. + /// \param ordered Kind of traversal. image2d<T> transformed_image_clipped_magn(double clip, bool ordered = true) const { return transformed_image_clipped_magn<T>(clip, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - - image2d<T> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<T>(1, ordered); - } - - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(bool ordered = true) const { return transformed_image_clipped_magn<R>(1, ordered); } - // FIXME: Find a more elegant way to fix range interval on a and b. - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [a; b].\n - ** - ** \param a Lower bound. - ** \param b Upper bound. - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \param ordered Kind of traversal. + image2d<T> transformed_image_clipped_magn(bool ordered = true) const + { + return transformed_image_clipped_magn<T>(1, ordered); + } + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [a, b]. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param a Lower bound. + /// \param b Upper bound. + /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_log_magn(double a, - double b, + image2d<R> transformed_image_log_magn(double a, double b, bool ordered = true) const { - // Check that R is real - // 0 <= a <= 1000 - // 0 <= b <= 1000 + /* FIXME: Find a more elegant way to fix range interval on a + and b (Note from Roland: what does it mean?). */ + + // FIXME: Ensure R is real. + // FIXME: Ensure 0 <= a <= 1000. + // FIXME: Ensure 0 <= b <= 1000. image2d<R> new_im(trans_im.domain()); double max = mln_min(double); mln_piter(image2d<R>) it(trans_im.domain()); - for_all(it) if (std::norm(trans_im(it)) > max) max = std::norm(trans_im(it)); - + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) opt::at(new_im, row, col) = - log(a + b * std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % trans_im.ncols()))) / - log (a + b * max) * mln_max(R); + log(a + b * std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols))) + / log (a + b * max) * mln_max(R); else { mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = log(a + b * std::norm(trans_im(it))) / - log (a + b * max) * mln_max(R); + new_im(it) = + log(a + b * std::norm(trans_im(it))) + / log (a + b * max) * mln_max(R); } - return new_im; } - // FIXME: Find a more elegant way to fix boundaries of a and b. - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [a; b].\n - ** - ** \param a Lower bound. - ** \param b Upper bound. - ** \param ordered Kind of traversal. - */ - - image2d<T> transformed_image_log_magn(double a, - double b, - bool ordered = true) const + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [a, b]. + /// + /// \param a Lower bound. + /// \param b Upper bound. + /// \param ordered Kind of traversal. + image2d<T> transformed_image_log_magn(double a, double b, + bool ordered = true) const { return transformed_image_log_magn<T>(a, b, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [1; 100].\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [1, 100]. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_log_magn(bool ordered = true) const { return transformed_image_log_magn<R>(1, 100, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [1; 100].\n - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [1, 100]. + /// + /// \param ordered Kind of traversal. image2d<T> transformed_image_log_magn(bool ordered = true) const { return transformed_image_log_magn<T>(1, 100, ordered); } - - /*! - ** \brief Destructor. - ** - ** Let memory free for other big images !!! - */ - ~_fft() + ~fft() { fftw_free(in); fftw_free(out); @@ -363,218 +329,222 @@ namespace mln { } protected: - - typename fft_trait<T>::fftw_input *in; ///< Input image. - std::complex<T> *out; ///< Complex image. - fftw_plan p; ///< Plan. - fftw_plan p_inv; ///< inverted plan. - image2d< std::complex<T> > trans_im; ///< Transformed image. - + /// Input image (used internally by FFTW). + typename fft_trait<T>::fftw_input* in; + /// Output (complex) image used internally by FFTW). + std::complex<T>* out; + /// (Forward) plan (command to be executed by FFTW). + fftw_plan p; + /// Inverse (backward) plan (command to be executed by FFTW). + fftw_plan p_inv; + /// Transformed image (``exported'' Milena image accessible + /// through accessors). + image2d< std::complex<T> > trans_im; }; - } // end of namespace internal + } // end of namespace mln::internal - namespace transform { - /*! - ** \brief fft class declaration. - ** - ** \param T Data type. - ** \param which Dispatch. - */ + namespace transform + { + + /// \brief FFT engine, producing a Discrete Fourier Transform + /// using FFTW (version 3). + /// + /// \tparam T Data type. + /// + /// \param which Dispatch traits. template <class T, internal::fft_dispatch which = internal::fft_trait<T>::which > class fft; - /*! - ** \brief fft specialization for fft real dispatch. - ** - ** \param T Data type. - */ + + /// \brief FFT engine (specialization for images of real values). + /// + /// \tparam T Data type. template <class T> - class fft<T, internal::fft_real> : public internal::_fft<T> + class fft<T, internal::fft_real> : public internal::fft<T> { - public: - - /*! - ** \brief Constructor. - ** - ** Initialization of data in order to compute the fft. - ** - ** \param original_im Image to process. - */ + /// \brief Constructor. + /// + /// Initialization of data for the computation of the FFT. + /// + /// \param original_im Image to process. template <typename D> fft(const image2d<D>& original_im) { - this->in = (T*) fftw_malloc(original_im.nrows() * original_im.ncols() * sizeof(T)); - this->out = (std::complex<T>*) - fftw_malloc(original_im.nrows() * (original_im.ncols() / 2 + 1) * sizeof(std::complex<T>)); - - for (unsigned row = 0; row < original_im.nrows(); ++row) - for (unsigned col = 0; col < original_im.ncols(); ++col) - this->in[row * original_im.ncols() + col] = opt::at(original_im, row, col); - - this->p = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), - this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_c2r_2d (original_im.nrows(), original_im.ncols(), - reinterpret_cast<fftw_complex*>(this->out), this->in, FFTW_ESTIMATE); - + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); + this->out = (std::complex<T>*) fftw_malloc(nrows + * (ncols / 2 + 1) + * sizeof(std::complex<T>)); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + this->in[row * ncols + col] = + opt::at(original_im, row, col); + + this->p = + fftw_plan_dft_r2c_2d (nrows, ncols, + this->in, + reinterpret_cast<fftw_complex*>(this->out), + FFTW_ESTIMATE); + this->p_inv = + fftw_plan_dft_c2r_2d (nrows, ncols, + reinterpret_cast<fftw_complex*>(this->out), + this->in, + FFTW_ESTIMATE); this->trans_im = image2d< std::complex<T> >(original_im.domain()); } - /*! - ** \brief Compute and return the transform. - */ + /// \brief Compute and return the transform (as a complex image). image2d< std::complex<T> > transform() { fftw_execute(this->p); - unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) { - this->out[i] = std::complex<T> (this->out[i].real() / denom, this->out[i].imag() / denom); + this->out[i] = std::complex<T> (this->out[i].real() / denom, + this->out[i].imag() / denom); opt::at(this->trans_im, row, col) = this->out[i]; ++i; } - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = this->trans_im.ncols() - 1; col > this->trans_im.ncols() / 2; --col) - opt::at(this->trans_im, row, col) = opt::at(this->trans_im, this->trans_im.nrows() - row - 1, - this->trans_im.ncols() - col - 1); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = ncols - 1; col > ncols / 2; --col) + opt::at(this->trans_im, row, col) = + opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); return this->trans_im; } - /*! - ** \brief Compute and return the invert transform. - ** - ** \param R Data type of output image. - */ + /// \brief Compute and return the inverse transform (as a real + /// image) of the FFT. + /// + /// \tparam R Value type of output image. template <class R> image2d<R> transform_inv() { - // FIXME : Check that R is real + // FIXME: Ensure R is real. - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) - this->out[row * (this->trans_im.ncols() / 2 + 1) + col] = + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + this->out[row * (ncols / 2 + 1) + col] = opt::at(this->trans_im, row, col); fftw_execute(this->p_inv); image2d<R> new_im(this->trans_im.domain()); int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) ? - (this->in[i] <= mln_max(R) ? - (R)this->in [i] : - mln_min(R)) : - mln_max(R)); + opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) + ? (this->in[i] <= mln_max(R) + ? (R) this->in [i] + : mln_min(R)) + : mln_max(R)); ++i; } return new_im; } - /*! - ** \brief Compute and return the invert transform. - */ + /// \brief Compute and return the inverse transform of the FFT. image2d<T> transform_inv() { return transform_inv<T>(); } - }; - /*! - ** \brief fft specialization for fft complex dispatch. - ** - ** \param T Data type. - ** \param R Complex representation kind. - */ + + /// \brief FFT engine (specialization for images of complex values). + /// + /// \tparam T Data type. template <class T> - class fft<T, internal::fft_cplx> : public internal::_fft<T> + class fft<T, internal::fft_cplx> : public internal::fft<T> { - public: - - /*! - ** \brief Constructor. - ** - ** Initialization of data in order to compute the fft. - ** - ** \param original_im Image to process. - */ + /// \brief Constructor. + /// + /// Initialization of data for the computation of the FFT. + /// + /// \param original_im Image to process. fft(const image2d< std::complex<T> >& original_im) { - this->in = fftw_malloc(sizeof(std::complex<T>) * - original_im.nrows() * original_im.ncols()); - this->out = fftw_malloc(sizeof(std::complex<T>) * - original_im.nrows() * original_im.ncols()); + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - for (unsigned row = 0; row < original_im.nrows(); ++row) - for (unsigned col = 0; col < original_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - this->in[row * original_im.ncols() + col].re = - original_im(row, col).real(); - this->in[row * original_im.ncols() + col].im = - original_im(row, col).imag(); + this->in[row * ncols + col].re = original_im(row, col).real(); + this->in[row * ncols + col].im = original_im(row, col).imag(); } - this->p = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), - this->in, this->out, - FFTW_FORWARD, FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), - this->out, this->in, + this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, + FFTW_FORWARD, FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, FFTW_BACKWARD, FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); } - /*! - ** \brief Compute and return the transform. - */ + /// \brief Compute and return the transform (as a complex image). image2d< std::complex<T> > transform() { fftw_execute(this->p); - unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { this->out[i].re /= denom; this->out[i].im /= denom; - this->trans_im(row, col) = std::complex<double>(this->out[i].re, - this->out[i].im); + this->trans_im(row, col) = + std::complex<double>(this->out[i].re, this->out[i].im); ++i; } return this->trans_im; } - /*! - ** \brief Compute and return the invert transform. - ** - ** \param R Data type of output image. - */ + /// \brief Compute and return the inverse transform (as a + /// complex image). + /// + /// \tparam R Component type of the value type of the output + /// image. template <class R> image2d< std::complex<R> > transform_inv() { - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - this->out[row * this->trans_im.ncols() + col].re = this->trans_im(row, col).real(); - this->out[row * this->trans_im.ncols() + col].im = this->trans_im(row, col).imag(); + this->out[row * ncols + col].re = + this->trans_im(row, col).real(); + this->out[row * ncols + col].im = + this->trans_im(row, col).imag(); } fftw_execute(this->p_inv); image2d< std::complex<R> > new_im(this->trans_im.size()); int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { new_im(row, col) = this->in[i]; ++i; @@ -582,9 +552,7 @@ namespace mln { return new_im; } - /*! - ** \brief Compute and return the invert transform. - */ + /// \brief Compute and return the inverse transform of the FFT. image2d< std::complex<T> > transform_inv() { return transform_inv<T>(); @@ -592,8 +560,8 @@ namespace mln { }; - } // end of namespace transform + } // end of namespace mln::transform -} // end of namespace oln +} // end of namespace mln #endif // ! MLN_TRANSFORM_FFT_HH -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-637-g8b31df6 Split the interface of the FFT off from its implementation.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 + milena/mln/transform/fft.hh | 722 ++++++++++++++++++++++++++----------------- 2 files changed, 436 insertions(+), 292 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 8f96157..72bfb1a 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Split the interface of the FFT off from its implementation. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Stylistic changes in the Fast Fourier Transform. * mln/transform/fft.hh: Here. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index 9b70b4f..36e4b80 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -42,9 +42,11 @@ namespace mln namespace internal { - /// Dispatch traits for fftw + /// Dispatch tags for FFT traits. enum fft_dispatch { fft_cplx, fft_real }; + + /// \brief FFT traits. template <typename T> struct fft_trait; @@ -68,6 +70,7 @@ namespace mln typedef std::complex<T> fftw_input; }; + /// \brief Internal structure containining packing data and /// instructions for the (forward) and inverse (backward) /// transforms. @@ -78,16 +81,10 @@ namespace mln { public: /// \brief Accessor to transformed image (const version). - const image2d< std::complex<T> >& transformed_image() const - { - return trans_im; - } + const image2d< std::complex<T> >& transformed_image() const; /// \brief Accessor to transformed image (non const version). - image2d< std::complex<T> >& transformed_image() - { - return trans_im; - } + image2d< std::complex<T> >& transformed_image(); /// \brief Accessor to the transformed image (magnitude). /// @@ -98,28 +95,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_magn(bool ordered = true) const - { - // FIXME: Ensure R is real. - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - image2d<R> new_im(trans_im.domain()); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - opt::at(new_im, row, col) = - std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)); - else - { - mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = std::norm(trans_im(it)); - } - return new_im; - } + image2d<R> transformed_image_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (magnitude). /// @@ -127,10 +103,7 @@ namespace mln /// an image containing <tt>|T[p]|</tt> is returned. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_magn(bool ordered = true) const - { - return transformed_image_magn<T>(ordered); - } + image2d<T> transformed_image_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -144,43 +117,7 @@ namespace mln /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(double clip, - bool ordered = true) const - { - // FIXME: Ensure R is real. - // FIXME: Ensure `clip' is between 0 and 1? - - double max = mln_min(double); - mln_piter(image2d<T>) it(trans_im.domain()); - for_all(it) - if (std::norm(trans_im(it)) > max) - max = std::norm(trans_im(it)); - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - image2d<R> new_im(trans_im.domain()); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - if (std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)) - >= max * clip) - opt::at(new_im, row, col) = mln_max(R); - else - opt::at(new_im, row, col) = - (double) mln_max(R) * - std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)) / (max * clip); - else - for_all(it) - if (std::norm(trans_im(it)) >= max * clip) - new_im(it) = mln_max(R); - else - new_im(it) = - (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); - return new_im; - } + bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -191,10 +128,7 @@ namespace mln /// \param clip Value used for clipping. /// \param ordered Kind of traversal. image2d<T> transformed_image_clipped_magn(double clip, - bool ordered = true) const - { - return transformed_image_clipped_magn<T>(clip, ordered); - } + bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -206,10 +140,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<R>(1, ordered); - } + image2d<R> transformed_image_clipped_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -218,10 +149,7 @@ namespace mln /// containing a clipped value of <tt>|T[p]|</tt> is returned. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<T>(1, ordered); - } + image2d<T> transformed_image_clipped_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -237,43 +165,7 @@ namespace mln /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_log_magn(double a, double b, - bool ordered = true) const - { - /* FIXME: Find a more elegant way to fix range interval on a - and b (Note from Roland: what does it mean?). */ - - // FIXME: Ensure R is real. - // FIXME: Ensure 0 <= a <= 1000. - // FIXME: Ensure 0 <= b <= 1000. - - image2d<R> new_im(trans_im.domain()); - - double max = mln_min(double); - mln_piter(image2d<R>) it(trans_im.domain()); - for_all(it) - if (std::norm(trans_im(it)) > max) - max = std::norm(trans_im(it)); - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - opt::at(new_im, row, col) = - log(a + b * std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols))) - / log (a + b * max) * mln_max(R); - else - { - mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = - log(a + b * std::norm(trans_im(it))) - / log (a + b * max) * mln_max(R); - } - return new_im; - } + bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -286,10 +178,7 @@ namespace mln /// \param b Upper bound. /// \param ordered Kind of traversal. image2d<T> transformed_image_log_magn(double a, double b, - bool ordered = true) const - { - return transformed_image_log_magn<T>(a, b, ordered); - } + bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -302,10 +191,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_log_magn(bool ordered = true) const - { - return transformed_image_log_magn<R>(1, 100, ordered); - } + image2d<R> transformed_image_log_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -315,18 +201,9 @@ namespace mln /// interval [1, 100]. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_log_magn(bool ordered = true) const - { - return transformed_image_log_magn<T>(1, 100, ordered); - } - - ~fft() - { - fftw_free(in); - fftw_free(out); - fftw_destroy_plan(p); - fftw_destroy_plan(p_inv); - } + image2d<T> transformed_image_log_magn(bool ordered = true) const; + + ~fft(); protected: /// Input image (used internally by FFTW). @@ -359,7 +236,7 @@ namespace mln class fft; - /// \brief FFT engine (specialization for images of real values). + /// \brief oFFT engine (specialization for images of real values). /// /// \tparam T Data type. template <class T> @@ -372,96 +249,20 @@ namespace mln /// /// \param original_im Image to process. template <typename D> - fft(const image2d<D>& original_im) - { - const unsigned nrows = original_im.nrows(); - const unsigned ncols = original_im.ncols(); - - this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); - this->out = (std::complex<T>*) fftw_malloc(nrows - * (ncols / 2 + 1) - * sizeof(std::complex<T>)); - - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - this->in[row * ncols + col] = - opt::at(original_im, row, col); - - this->p = - fftw_plan_dft_r2c_2d (nrows, ncols, - this->in, - reinterpret_cast<fftw_complex*>(this->out), - FFTW_ESTIMATE); - this->p_inv = - fftw_plan_dft_c2r_2d (nrows, ncols, - reinterpret_cast<fftw_complex*>(this->out), - this->in, - FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); - } + fft(const image2d<D>& original_im); /// \brief Compute and return the transform (as a complex image). - image2d< std::complex<T> > transform() - { - fftw_execute(this->p); - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - unsigned denom = nrows * ncols; - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col <= ncols / 2; ++col) - { - this->out[i] = std::complex<T> (this->out[i].real() / denom, - this->out[i].imag() / denom); - opt::at(this->trans_im, row, col) = this->out[i]; - ++i; - } - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = ncols - 1; col > ncols / 2; --col) - opt::at(this->trans_im, row, col) = - opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); - return this->trans_im; - } + image2d< std::complex<T> > transform(); /// \brief Compute and return the inverse transform (as a real /// image) of the FFT. /// /// \tparam R Value type of output image. template <class R> - image2d<R> transform_inv() - { - // FIXME: Ensure R is real. - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col <= ncols / 2; ++col) - this->out[row * (ncols / 2 + 1) + col] = - opt::at(this->trans_im, row, col); - - fftw_execute(this->p_inv); - - image2d<R> new_im(this->trans_im.domain()); - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) - ? (this->in[i] <= mln_max(R) - ? (R) this->in [i] - : mln_min(R)) - : mln_max(R)); - ++i; - } - return new_im; - } + image2d<R> transform_inv(); /// \brief Compute and return the inverse transform of the FFT. - image2d<T> transform_inv() - { - return transform_inv<T>(); - } + image2d<T> transform_inv(); }; @@ -477,48 +278,10 @@ namespace mln /// Initialization of data for the computation of the FFT. /// /// \param original_im Image to process. - fft(const image2d< std::complex<T> >& original_im) - { - const unsigned nrows = original_im.nrows(); - const unsigned ncols = original_im.ncols(); - - this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->in[row * ncols + col].re = original_im(row, col).real(); - this->in[row * ncols + col].im = original_im(row, col).imag(); - } - - this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, - FFTW_FORWARD, FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, - FFTW_BACKWARD, FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); - } + fft(const image2d< std::complex<T> >& original_im); /// \brief Compute and return the transform (as a complex image). - image2d< std::complex<T> > transform() - { - fftw_execute(this->p); - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - unsigned denom = nrows * ncols; - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->out[i].re /= denom; - this->out[i].im /= denom; - this->trans_im(row, col) = - std::complex<double>(this->out[i].re, this->out[i].im); - ++i; - } - return this->trans_im; - } + image2d< std::complex<T> > transform(); /// \brief Compute and return the inverse transform (as a /// complex image). @@ -526,42 +289,417 @@ namespace mln /// \tparam R Component type of the value type of the output /// image. template <class R> - image2d< std::complex<R> > transform_inv() - { - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->out[row * ncols + col].re = - this->trans_im(row, col).real(); - this->out[row * ncols + col].im = - this->trans_im(row, col).imag(); - } - - fftw_execute(this->p_inv); - - image2d< std::complex<R> > new_im(this->trans_im.size()); - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - new_im(row, col) = this->in[i]; - ++i; - } - return new_im; - } + image2d< std::complex<R> > transform_inv(); /// \brief Compute and return the inverse transform of the FFT. - image2d< std::complex<T> > transform_inv() - { - return transform_inv<T>(); - } - + image2d< std::complex<T> > transform_inv(); }; } // end of namespace mln::transform + + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <class T> + inline + const image2d< std::complex<T> >& + fft<T>::transformed_image() const + { + return trans_im; + } + + template <class T> + inline + image2d< std::complex<T> >& + fft<T>::transformed_image() + { + return trans_im; + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_magn(bool ordered) const + { + // FIXME: Ensure R is real. + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + for_all(it) + new_im(it) = std::norm(trans_im(it)); + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_magn(bool ordered) const + { + return transformed_image_magn<T>(ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_clipped_magn(double clip, bool ordered) const + { + // FIXME: Ensure R is real. + // FIXME: Ensure `clip' is between 0 and 1? + + double max = mln_min(double); + mln_piter(image2d<T>) it(trans_im.domain()); + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + if (std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) + >= max * clip) + opt::at(new_im, row, col) = mln_max(R); + else + opt::at(new_im, row, col) = + (double) mln_max(R) * + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) / (max * clip); + else + for_all(it) + if (std::norm(trans_im(it)) >= max * clip) + new_im(it) = mln_max(R); + else + new_im(it) = + (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_clipped_magn(double clip, bool ordered) const + { + return transformed_image_clipped_magn<T>(clip, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_clipped_magn(bool ordered) const + { + return transformed_image_clipped_magn<R>(1, ordered); + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_clipped_magn(bool ordered) const + { + return transformed_image_clipped_magn<T>(1, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_log_magn(double a, double b, bool ordered) const + { + /* FIXME: Find a more elegant way to fix range interval on a + and b (Note from Roland: what does it mean?). */ + + // FIXME: Ensure R is real. + // FIXME: Ensure 0 <= a <= 1000. + // FIXME: Ensure 0 <= b <= 1000. + + image2d<R> new_im(trans_im.domain()); + + double max = mln_min(double); + mln_piter(image2d<R>) it(trans_im.domain()); + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + log(a + b * std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols))) + / log (a + b * max) * mln_max(R); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + for_all(it) + new_im(it) = + log(a + b * std::norm(trans_im(it))) + / log (a + b * max) * mln_max(R); + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_log_magn(double a, double b, bool ordered) const + { + return transformed_image_log_magn<T>(a, b, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_log_magn(bool ordered) const + { + return transformed_image_log_magn<R>(1, 100, ordered); + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_log_magn(bool ordered) const + { + return transformed_image_log_magn<T>(1, 100, ordered); + } + + template <class T> + inline + fft<T>::~fft() + { + fftw_free(in); + fftw_free(out); + fftw_destroy_plan(p); + fftw_destroy_plan(p_inv); + } + + } // end of namespace mln::internal + + + namespace transform + { + + /*--------------------------------------------------------. + | FFT engine (specialization for images of real values). | + `--------------------------------------------------------*/ + + template <class T> + template <typename D> + inline + fft<T, internal::fft_real>::fft(const image2d<D>& original_im) + { + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); + this->out = (std::complex<T>*) fftw_malloc(nrows + * (ncols / 2 + 1) + * sizeof(std::complex<T>)); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + this->in[row * ncols + col] = + opt::at(original_im, row, col); + + this->p = + fftw_plan_dft_r2c_2d (nrows, ncols, + this->in, + reinterpret_cast<fftw_complex*>(this->out), + FFTW_ESTIMATE); + this->p_inv = + fftw_plan_dft_c2r_2d (nrows, ncols, + reinterpret_cast<fftw_complex*>(this->out), + this->in, + FFTW_ESTIMATE); + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_real>::transform() + { + fftw_execute(this->p); + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + { + this->out[i] = std::complex<T> (this->out[i].real() / denom, + this->out[i].imag() / denom); + opt::at(this->trans_im, row, col) = this->out[i]; + ++i; + } + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = ncols - 1; col > ncols / 2; --col) + opt::at(this->trans_im, row, col) = + opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); + return this->trans_im; + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T, internal::fft_real>::transform_inv() + { + // FIXME: Ensure R is real. + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + this->out[row * (ncols / 2 + 1) + col] = + opt::at(this->trans_im, row, col); + + fftw_execute(this->p_inv); + + image2d<R> new_im(this->trans_im.domain()); + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) + ? (this->in[i] <= mln_max(R) + ? (R) this->in [i] + : mln_min(R)) + : mln_max(R)); + ++i; + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T, internal::fft_real>::transform_inv() + { + return transform_inv<T>(); + } + + + /*-----------------------------------------------------------. + | FFT engine (specialization for images of complex values). | + `-----------------------------------------------------------*/ + + template <class T> + inline + fft<T, internal::fft_cplx>::fft(const image2d< std::complex<T> >& original_im) + { + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->in[row * ncols + col].re = original_im(row, col).real(); + this->in[row * ncols + col].im = original_im(row, col).imag(); + } + + this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, + FFTW_FORWARD, FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, + FFTW_BACKWARD, FFTW_ESTIMATE); + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_cplx>::transform() + { + fftw_execute(this->p); + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->out[i].re /= denom; + this->out[i].im /= denom; + this->trans_im(row, col) = + std::complex<double>(this->out[i].re, this->out[i].im); + ++i; + } + return this->trans_im; + } + + template <class T> + template <class R> + inline + image2d< std::complex<R> > + fft<T, internal::fft_cplx>::transform_inv() + { + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->out[row * ncols + col].re = + this->trans_im(row, col).real(); + this->out[row * ncols + col].im = + this->trans_im(row, col).imag(); + } + + fftw_execute(this->p_inv); + + image2d< std::complex<R> > new_im(this->trans_im.size()); + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + new_im(row, col) = this->in[i]; + ++i; + } + return new_im; + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_cplx>::transform_inv() + { + return transform_inv<T>(); + } + + } // end of namespace mln::transform + +# endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln #endif // ! MLN_TRANSFORM_FFT_HH -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-627-g5f1d602 Import the FFTW-based Fast Fourier Transform.

by Alexandre Abraham

* mln/transform/fft.hh, * tests/transform/fft.cc: New. Imported from the old directory sanbox/abraham/. Signed-off-by: Roland Levillain <roland(a)lrde.epita.fr> --- milena/ChangeLog | 9 + milena/mln/transform/fft.hh | 693 +++++++++++++++++++++++++++++++++++++++++ milena/tests/transform/fft.cc | 92 ++++++ 3 files changed, 794 insertions(+) create mode 100644 milena/mln/transform/fft.hh create mode 100644 milena/tests/transform/fft.cc diff --git a/milena/ChangeLog b/milena/ChangeLog index 2e0ed94..83a2c76 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,12 @@ +2012-10-09 Alexandre Abraham <abraham(a)lrde.epita.fr> + + Import the FFTW-based Fast Fourier Transform. + + * mln/transform/fft.hh, + * tests/transform/fft.cc: + New. + Imported from the old directory sanbox/abraham/. + 2013-08-23 Roland Levillain <roland(a)lrde.epita.fr> Fix a couple of Milena tests. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh new file mode 100644 index 0000000..0dafa7b --- /dev/null +++ b/milena/mln/transform/fft.hh @@ -0,0 +1,693 @@ +// Copyright (C) 2007, 2008, 2009 EPITA Research and Development Laboratory +// +// This file is part of the Milena 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, 51 Franklin Street, Fifth Floor, +// 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 MLN_TRANSFORM_FFT_HH +# define MLN_TRANSFORM_FFT_HH + +# include <mln/core/image/image2d.hh> +# include <mln/estim/min_max.hh> +# include <mln/opt/at.hh> + +# include <complex> + +# include <fftw3.h> + +namespace mln { + + namespace internal { + + /// Dispatch traits for fftw + enum fft_dispatch { fft_cplx, fft_real }; + + template <typename T> + struct fft_trait; + + /*! + ** \brief fft trait. + */ + template <> + struct fft_trait<double> + { + static const fft_dispatch which = fft_real; ///< Real dispatch. + typedef double fftw_input; ///< Type of input. + }; + + /*! + ** \brief fft trait. + ** + ** \specialization for ntg::cplx<R, T> + */ + template <typename T> + struct fft_trait< std::complex<T> > + { + static const fft_dispatch which = fft_cplx; ///< Complex dispatch. + typedef std::complex <T> fftw_input; ///< Type of input. + }; + + /*! + ** _fft<ntg::cplx_representation, T> + ** + ** \param T Data type. + */ + template <class T> + class _fft + { + public: + /*! + ** \brief Accessor to transformed image. + ** + ** Const version. + */ + const image2d< std::complex<T> > transformed_image() const + { + return trans_im; + } + + /*! + ** \brief Accessor to transformed image. + ** + ** Non const version. + */ + image2d< std::complex<T> >& transformed_image() + { + return trans_im; + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** |T[p]|. + ** + ** \tparam R Data type of the resulting image. + ** + ** \param ordered Kind of traversal. + */ + template <class R> + image2d<R> transformed_image_magn(bool ordered = true) const + { + // FIXME : check that R is real + + image2d<R> new_im(trans_im.domain()); + + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), + (col + trans_im.ncols() / 2) % trans_im.ncols())); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + + for_all(it) + new_im(it) = std::norm(trans_im(it)); + } + + return new_im; + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** |T[p]|. + ** + ** \param ordered Kind of traversal. + */ + image2d<T> transformed_image_magn(bool ordered = true) const + { + return transformed_image_magn<T>(ordered); + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a clipped value of |T[p]|.\n + ** + ** \param R Data type of the resulting image. + ** + ** \param clip Value used for clipping. + ** \param ordered Kind of traversal. + */ + + template <class R> + image2d<R> transformed_image_clipped_magn(const double clip, + bool ordered = true) const + { + // Check that R is real + + image2d<R> new_im(trans_im.domain()); + // check that clip is >=0 and <=1 ? + double max; + mln_piter(image2d<T>) it(trans_im.domain()); + + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + { + if (std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), + (col + trans_im.ncols() / 2) % + trans_im.ncols())) >= max * clip) + opt::at(new_im, row, col) = mln_max(R); + else + opt::at(new_im, row, col) = + (double) mln_max(R) * + std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), + (col + trans_im.ncols() / 2) % + trans_im.ncols())) / (max * clip); + } + else + { + for_all(it) + { + if (std::norm(trans_im(it)) >= max * clip) + new_im(it) = mln_max(R); + else + new_im(it) = (double) mln_max(R) * + std::norm(trans_im(it)) / (max * clip); + } + } + + return new_im; + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a clipped value of |T[p]|.\n + ** + ** \param clip Value used for clipping. + ** \param ordered Kind of traversal. + */ + image2d<T> transformed_image_clipped_magn(const double clip, + bool ordered = true) const + { + return transformed_image_clipped_magn<T>(clip, ordered); + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a clipped value of |T[p]|.\n + ** + ** \param R Data type of the resulting image. + ** + ** \param ordered Kind of traversal. + */ + + image2d<T> transformed_image_clipped_magn(bool ordered = true) const + { + return transformed_image_clipped_magn<T>(1, ordered); + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a clipped value of |T[p]|.\n + ** + ** \param R Data type of the resulting image. + ** + ** \param ordered Kind of traversal. + */ + + template <class R> + image2d<R> transformed_image_clipped_magn(bool ordered = true) const + { + return transformed_image_clipped_magn<R>(1, ordered); + } + + // FIXME: Find a more elegant way to fix range interval on a and b. + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a log translated value of |T[p]| on interval [a; b].\n + ** + ** \param a Lower bound. + ** \param b Upper bound. + ** \param ordered Kind of traversal. + */ + + template <class R> + image2d<R> transformed_image_log_magn(double a, + double b, + bool ordered = true) const + { + // Check that R is real + // 0 <= a <= 1000 + // 0 <= b <= 1000 + + image2d<R> new_im(trans_im.domain()); + + double max = 0; + mln_piter(image2d<R>) it(trans_im.domain()); + + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + log(a + b * std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), + (col + trans_im.ncols() / 2) % trans_im.ncols()))) / + log (a + b * max) * mln_max(R); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + + for_all(it) + new_im(it) = log(a + b * std::norm(trans_im(it))) / + log (a + b * max) * mln_max(R); + } + + return new_im; + } + + // FIXME: Find a more elegant way to fix boundaries of a and b. + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a log translated value of |T[p]| on interval [a; b].\n + ** + ** \param a Lower bound. + ** \param b Upper bound. + ** \param ordered Kind of traversal. + */ + + image2d<T> transformed_image_log_magn(double a, + double b, + bool ordered = true) const + { + return transformed_image_log_magn<T>(a, b, ordered); + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a log translated value of |T[p]| on interval [1; 100].\n + ** + ** \param R Data type of the resulting image. + ** + ** \param ordered Kind of traversal. + */ + + template <class R> + image2d<R> transformed_image_log_magn(bool ordered = true) const + { + return transformed_image_log_magn<R>(1, 100, ordered); + } + + /*! + ** \brief Accessor to transformed image. + ** + ** For each point p of the transformed image T, you get + ** a log translated value of |T[p]| on interval [1; 100].\n + ** + ** \param ordered Kind of traversal. + */ + + image2d<T> transformed_image_log_magn(bool ordered = true) const + { + return transformed_image_log_magn<T>(1, 100, ordered); + } + + + /*! + ** \brief Destructor. + ** + ** Let memory free for other big images !!! + */ + ~_fft() + { + fftw_free(in); + fftw_free(out); + fftw_destroy_plan(p); + fftw_destroy_plan(p_inv); + } + + protected: + + typename fft_trait<T>::fftw_input *in; ///< Input image. + std::complex<T> *out; ///< Complex image. + fftw_plan p; ///< Plan. + fftw_plan p_inv; ///< inverted plan. + image2d< std::complex<T> > trans_im; ///< Transformed image. + + }; + + } // end of namespace internal + + namespace transform { + + /*! + ** \brief fft class declaration. + ** + ** \param T Data type. + ** \param which Dispatch. + */ + template <class T, + internal::fft_dispatch which = internal::fft_trait<T>::which > + class fft; + + /*! + ** \brief fft specialization for fft real dispatch. + ** + ** \param T Data type. + */ + template <class T> + class fft<T, internal::fft_real> : public internal::_fft<T> + { + + public: + + /*! + ** \brief Constructor. + ** + ** Initialization of data in order to compute the fft. + ** + ** \param original_im Image to process. + */ + template <typename D> + fft(const image2d<D>& original_im) + { + this->in = (T*) fftw_malloc(original_im.nrows() * original_im.ncols() * sizeof(T)); + this->out = (std::complex<T>*) + fftw_malloc(original_im.nrows() * (original_im.ncols() / 2 + 1) * sizeof(std::complex<T>)); + + for (unsigned row = 0; row < original_im.nrows(); ++row) + for (unsigned col = 0; col < original_im.ncols(); ++col) + this->in[row * original_im.ncols() + col] = opt::at(original_im, row, col); + + this->p = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), + this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), + this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); + + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + /*! + ** \brief Compute and return the transform. + */ + image2d< std::complex<T> > transform() + { + fftw_execute(this->p); + + unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + int i = 0; + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) + { + this->out[i] = std::complex<T> (this->out[i].real() / denom, this->out[i].imag() / denom); + opt::at(this->trans_im, row, col) = this->out[i]; + ++i; + } + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = this->trans_im.ncols() - 1; col > this->trans_im.ncols() / 2; --col) + opt::at(this->trans_im, row, col) = opt::at(this->trans_im, this->trans_im.nrows() - row - 1, + this->trans_im.ncols() - col - 1); + return this->trans_im; + } + + /*! + ** \brief Compute and return the invert transform. + ** + ** \param R Data type of output image. + */ + template <class R> + image2d<R> transform_inv() + { + // FIXME : Check that R is real + + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) + this->out[row * (this->trans_im.ncols() / 2 + 1) + col] = + opt::at(this->trans_im, row, col).real(); + + fftw_execute(this->p_inv); + + image2d<R> new_im(this->trans_im.domain()); + int i = 0; + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + { + opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) ? + (this->in[i] <= mln_max(R) ? + (R)this->in [i] : + mln_min(R)) : + mln_max(R)); + ++i; + } + return new_im; + } + + /*! + ** \brief Shift zero-frequency component of discrete Fourier transform + ** to center of spectrum. + ** + ** \param R 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 R> + image2d<R> shift_transform_inv() + { + image2d<R> t = transform_inv<R>(); + image2d<R> 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<R> > 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() + { + 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>(); + } + */ + + }; + + /*! + ** \brief fft specialization for fft complex dispatch. + ** + ** \param T Data type. + ** \param R Complex representation kind. + */ + template <class T> + class fft<T, internal::fft_cplx> : public internal::_fft<T> + { + + public: + + /*! + ** \brief Constructor. + ** + ** Initialization of data in order to compute the fft. + ** + ** \param original_im Image to process. + */ + fft(const image2d< std::complex<T> >& original_im) + { + this->in = fftw_malloc(sizeof(std::complex<T>) * + original_im.nrows() * original_im.ncols()); + this->out = fftw_malloc(sizeof(std::complex<T>) * + original_im.nrows() * original_im.ncols()); + + for (unsigned row = 0; row < original_im.nrows(); ++row) + for (unsigned col = 0; col < original_im.ncols(); ++col) + { + this->in[row * original_im.ncols() + col].re = + original_im(row, col).real(); + this->in[row * original_im.ncols() + col].im = + original_im(row, col).imag(); + } + + this->p = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), + this->in, this->out, + FFTW_FORWARD, FFTW_ESTIMATE); + this->p_inv = fftw2d_plan_dft_2d(original_im.nrows(), original_im.ncols(), + this->in, this->out, + FFTW_BACKWARD, FFTW_ESTIMATE); + + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + /*! + ** \brief Compute and return the transform. + */ + image2d< std::complex<T> > transform() + { + fftw_execute(this->p); + + unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + int i = 0; + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + { + this->out[i].re /= denom; + this->out[i].im /= denom; + this->trans_im(row, col) = std::complex<double>(this->out[i].re, + this->out[i].im); + ++i; + } + return this->trans_im; + } + + /*! + ** \brief Compute and return the invert transform. + ** + ** \param R Data type of output image. + */ + template <class R> + image2d< std::complex<R> > transform_inv() + { + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + { + this->out[row * this->trans_im.ncols() + col].re = this->trans_im(row, col).real(); + this->out[row * this->trans_im.ncols() + col].im = this->trans_im(row, col).imag(); + } + + fftw_execute(this->p_inv); + + image2d< std::complex<R> > new_im(this->trans_im.size()); + int i = 0; + for (unsigned row = 0; row < this->trans_im.nrows(); ++row) + for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + { + new_im(row, col) = this->in[i]; + ++i; + } + return new_im; + } + + /*! + ** \brief Compute and return the invert transform. + */ + image2d< std::complex<T> > transform_inv() + { + return transform_inv<T>(); + } + + }; + + } // end of namespace transform + +} // end of namespace oln + +#endif // ! MLN_TRANSFORM_FFT_HH diff --git a/milena/tests/transform/fft.cc b/milena/tests/transform/fft.cc new file mode 100644 index 0000000..d7f7e63 --- /dev/null +++ b/milena/tests/transform/fft.cc @@ -0,0 +1,92 @@ +// -*- c++ -*- +// Copyright (C) 2004 EPITA Research and Development Laboratory +// +// This file is part of the Milena 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, Inc., 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301, 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. + +#include <mln/transform/fft.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/opt/at.hh> +#include <mln/debug/println.hh> + +#define CHECK(Condition) \ + if (Condition) \ + std::cout << "OK" << std::endl; \ + else \ + { \ + std::cout << "FAIL" << std::endl; \ + status = 1; \ + } + +using namespace mln; +using namespace transform; +using namespace value; + +int main () +{ + int status = 0; + + image2d<int_u8> im1; + io::pgm::load(im1, "lena.pgm"); + + fft<double> fourier(im1); + + image2d< std::complex<double> > im2 = fourier.transform(); + + image2d<int_u8> im3 = fourier.transform_inv<int_u8>(); + + io::pgm::save(im3, "fft_copy.pgm"); + + image2d<int_u8> fft = fourier.transformed_image_log_magn<int_u8>(true); + // debug::println(fourier.transformed_image()); + io::pgm::save(fft, "fft.pgm"); + + std::cout << "Test: Image == F-1(F(Image)) ... " << std::flush; + // CHECK (im1 == im3); + + image2d<int_u8> out = fourier.transformed_image_clipped_magn<int_u8>(0.01); + + io::pgm::save(out, "fft_trans_clipped.pgm"); + + out = fourier.transformed_image_log_magn<int_u8>(1, 100); + + io::pgm::save(out, "fft_trans_log.pgm"); + + for (int row = 40; row < im2.nrows() - 40; ++row) + for (int col = 0; col < im2.ncols(); ++col) + opt::at(im2, row, col) = 0; + + for (int row = 0; row < im2.nrows(); ++row) + for (int col = 40; col < im2.ncols() - 40; ++col) + opt::at(im2, row, col) = 0; + + out = fourier.transform_inv<int_u8>(); + + io::pgm::save(out, "fft_low_pass.pgm"); + + return status; +} -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-628-g31f1d70 Ask configure to try to find FFTW (version 3).

by Roland Levillain

* configure.ac: Here. --- ChangeLog | 6 ++++++ configure.ac | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/ChangeLog b/ChangeLog index affb146..dd4134d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + + Ask configure to try to find FFTW (version 3). + + * configure.ac: Here. + 2013-08-22 Roland Levillain <roland(a)lrde.epita.fr> * AUTHORS: Update. diff --git a/configure.ac b/configure.ac index c8a46c8..095e81a 100644 --- a/configure.ac +++ b/configure.ac @@ -315,6 +315,12 @@ OLN_WITH_LIB([GDCM], [gdcm-2.0/gdcmReader.h], [gdcmCommon], [gdcm], [ -lgdcmexpat -lgdcmjpeg12 -lgdcmjpeg16 -lgdcmjpeg8], [ -lgdcmopenjpeg -lgdcmuuid -lgdcmzlib])) +# --------------- # +# FFTW3 library. # +# --------------- # + +OLN_WITH_LIB([FFTW3], [fftw3.h], [fftw3], [fftw3], [FFTW3], [-lm]) + ## --------- ## ## Swilena. ## -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-629-gbee3fa8 Fix the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh (mln::transform::fft<T, internal::fft_real>::fft): Properly initialize the plan of the inverse transform using fftw_plan_dft_c2r_2d. (mln::transform::fft<T, internal::fft_real>::transform_inv): Use the entire complex values (not just their real parts) to compute the inverse transform. (mln::transform::class fft<T, internal::fft_cplx>::fft): Call the proper function (fftw_plan_dft_2d) to create the plan of the inverse transform. Switch the input and the output so that the inverse transform uses the right images. --- milena/ChangeLog | 17 +++++++++++++++++ milena/mln/transform/fft.hh | 15 ++++++++------- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 83a2c76..8bc21f5 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,20 @@ +2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + + Fix the Fast Fourier Transform. + + * mln/transform/fft.hh + (mln::transform::fft<T, internal::fft_real>::fft): + Properly initialize the plan of the inverse transform using + fftw_plan_dft_c2r_2d. + (mln::transform::fft<T, internal::fft_real>::transform_inv): + Use the entire complex values (not just their real parts) to + compute the inverse transform. + (mln::transform::class fft<T, internal::fft_cplx>::fft): + Call the proper function (fftw_plan_dft_2d) to create the plan of + the inverse transform. + Switch the input and the output so that the inverse transform uses + the right images. + 2012-10-09 Alexandre Abraham <abraham(a)lrde.epita.fr> Import the FFTW-based Fast Fourier Transform. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index 0dafa7b..cf0f1bb 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2007, 2008, 2009 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008, 2009, 2012 EPITA Research and Development +// Laboratory // // This file is part of the Milena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -418,8 +419,8 @@ namespace mln { this->p = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), - this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_c2r_2d (original_im.nrows(), original_im.ncols(), + reinterpret_cast<fftw_complex*>(this->out), this->in, FFTW_ESTIMATE); this->trans_im = image2d< std::complex<T> >(original_im.domain()); } @@ -460,7 +461,7 @@ namespace mln { for (unsigned row = 0; row < this->trans_im.nrows(); ++row) for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) this->out[row * (this->trans_im.ncols() / 2 + 1) + col] = - opt::at(this->trans_im, row, col).real(); + opt::at(this->trans_im, row, col); fftw_execute(this->p_inv); @@ -620,9 +621,9 @@ namespace mln { this->p = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), this->in, this->out, FFTW_FORWARD, FFTW_ESTIMATE); - this->p_inv = fftw2d_plan_dft_2d(original_im.nrows(), original_im.ncols(), - this->in, this->out, - FFTW_BACKWARD, FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), + this->out, this->in, + FFTW_BACKWARD, FFTW_ESTIMATE); this->trans_im = image2d< std::complex<T> >(original_im.domain()); } -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-630-g02720ff Exercise the Fast Fourier Transform.

by Roland Levillain

* tests/transform/fft.cc: Fix the path to the input image. * tests/transform/Makefile.am (check_PROGRAMS) [HAVE_FFTW3]: Add fft. (fft_SOURCES, fft_CPPFLAGS, fft_LDFLAGS) [HAVE_FFTW3]: New. (MOSTLYCLEANFILES): New. --- milena/ChangeLog | 10 ++++++++++ milena/tests/transform/Makefile.am | 22 ++++++++++++++++++++-- milena/tests/transform/fft.cc | 6 ++++-- 3 files changed, 34 insertions(+), 4 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 8bc21f5..fc060d8 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,15 @@ 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + Exercise the Fast Fourier Transform. + + * tests/transform/fft.cc: Fix the path to the input image. + * tests/transform/Makefile.am + (check_PROGRAMS) [HAVE_FFTW3]: Add fft. + (fft_SOURCES, fft_CPPFLAGS, fft_LDFLAGS) [HAVE_FFTW3]: New. + (MOSTLYCLEANFILES): New. + +2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> + Fix the Fast Fourier Transform. * mln/transform/fft.hh diff --git a/milena/tests/transform/Makefile.am b/milena/tests/transform/Makefile.am index 40de7f1..b266926 100644 --- a/milena/tests/transform/Makefile.am +++ b/milena/tests/transform/Makefile.am @@ -1,5 +1,5 @@ -# Copyright (C) 2008, 2009, 2010 EPITA Research and Development -# Laboratory (LRDE) +# Copyright (C) 2008, 2009, 2010, 2012 EPITA Research and Development +# Laboratory (LRDE). # # This file is part of Olena. # @@ -28,4 +28,22 @@ check_PROGRAMS = \ influence_zone_geodesic_saturated + +if HAVE_FFTW3 + +check_PROGRAMS += fft +fft_SOURCES = fft.cc +fft_CPPFLAGS = $(AM_CPPFLAGS) $(FFTW3_CPPFLAGS) +fft_LDFLAGS = $(AM_LDFLAGS) $(FFTW3_LDFLAGS) + +endif HAVE_FFTW3 + + TESTS = $(check_PROGRAMS) + +MOSTLYCLEANFILES = \ + fft.pgm \ + fft_copy.pgm \ + fft_low_pass.pgm \ + fft_trans_clipped.pgm \ + fft_trans_log.pgm diff --git a/milena/tests/transform/fft.cc b/milena/tests/transform/fft.cc index d7f7e63..bf7b8b9 100644 --- a/milena/tests/transform/fft.cc +++ b/milena/tests/transform/fft.cc @@ -1,5 +1,5 @@ // -*- c++ -*- -// Copyright (C) 2004 EPITA Research and Development Laboratory +// Copyright (C) 2004, 2012 EPITA Research and Development Laboratory // // This file is part of the Milena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -33,6 +33,8 @@ #include <mln/opt/at.hh> #include <mln/debug/println.hh> +#include "tests/data.hh" + #define CHECK(Condition) \ if (Condition) \ std::cout << "OK" << std::endl; \ @@ -51,7 +53,7 @@ int main () int status = 0; image2d<int_u8> im1; - io::pgm::load(im1, "lena.pgm"); + io::pgm::load(im1, MLN_IMG_DIR "/lena.pgm"); fft<double> fourier(im1); -- 1.7.10.4

11 years, 4 months

Jump to page:

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

Olena-patches