October 2012 - Olena-patches

olena-2.0-44-g027f767 Add functions.
by Guillaume Lazzara 18 Oct '12

18 Oct '12

* mln/fun/vv2v/span.hh, * mln/fun/vvvv2v/max.hh, * mln/fun/vvvv2v/min.hh, * mln/fun/vvvv2v/span.hh: New. * mln/math/max.hh, * mln/math/min.hh: New overloads. --- milena/ChangeLog | 12 ++++++++++ milena/mln/fun/vv2v/{mean.hh => span.hh} | 33 ++++++++++++++------------- milena/mln/fun/vvvv2v/{mean.hh => max.hh} | 20 ++++++++-------- milena/mln/fun/vvvv2v/mean.hh | 2 +- milena/mln/fun/vvvv2v/{mean.hh => min.hh} | 20 ++++++++-------- milena/mln/fun/vvvv2v/{mean.hh => span.hh} | 34 ++++++++++++++++----------- milena/mln/math/max.hh | 17 ++++++++++++++ milena/mln/math/min.hh | 17 ++++++++++++++ 8 files changed, 104 insertions(+), 51 deletions(-) copy milena/mln/fun/vv2v/{mean.hh => span.hh} (70%) copy milena/mln/fun/vvvv2v/{mean.hh => max.hh} (79%) copy milena/mln/fun/vvvv2v/{mean.hh => min.hh} (79%) copy milena/mln/fun/vvvv2v/{mean.hh => span.hh} (67%) diff --git a/milena/ChangeLog b/milena/ChangeLog index 7493137..5d2bea5 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,15 @@ +2012-10-18 Guillaume Lazzara <z(a)lrde.epita.fr> + + Add functions. + + * mln/fun/vv2v/span.hh, + * mln/fun/vvvv2v/max.hh, + * mln/fun/vvvv2v/min.hh, + * mln/fun/vvvv2v/span.hh: New. + + * mln/math/max.hh, + * mln/math/min.hh: New overloads. + 2012-10-05 Guillaume Lazzara <z(a)lrde.epita.fr> Improve accu::median_interval. diff --git a/milena/mln/fun/vv2v/mean.hh b/milena/mln/fun/vv2v/span.hh similarity index 70% copy from milena/mln/fun/vv2v/mean.hh copy to milena/mln/fun/vv2v/span.hh index 510b9a6..e9cd002 100644 --- a/milena/mln/fun/vv2v/mean.hh +++ b/milena/mln/fun/vv2v/span.hh @@ -23,15 +23,15 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_FUN_VV2V_MEAN_HH -# define MLN_FUN_VV2V_MEAN_HH +#ifndef MLN_FUN_VV2V_SPAN_HH +# define MLN_FUN_VV2V_SPAN_HH /// \file /// -/// Functor that computes the mean of two values. +/// Functor that computes the spanimum of two values. # include <mln/core/concept/function.hh> -# include <mln/math/mean.hh> +# include <mln/value/interval.hh> namespace mln @@ -45,24 +45,25 @@ namespace mln // FIXME: Doc. - /// \brief A functor computing the mean of two values. - template <typename L, typename R = L> - struct mean : public Function_vv2v< mean<L,R> >, - private mlc_converts_to(R,L)::check_t + /// \brief A functor computing the span of two interval values. + template <typename T> + struct span : public Function_vv2v< span<T> > { - typedef R result; - R operator()(const L& v1, const L& v2) const; + typedef value::interval<T> result; + + value::interval<T> operator()(const value::interval<T>& v1, + const value::interval<T>& v2) const; }; # ifndef MLN_INCLUDE_ONLY - template <typename L, typename R> - inline - R - mean<L,R>::operator()(const L& v1, const L& v2) const + template <typename T> + value::interval<T> + span<T>::operator()(const value::interval<T>& v1, + const value::interval<T>& v2) const { - return R(mln::math::mean(v1, v2); + return value::span(v1, v2); } # endif // ! MLN_INCLUDE_ONLY @@ -74,4 +75,4 @@ namespace mln } // end of namespace mln -#endif // ! MLN_FUN_VV2V_MEAN_HH +#endif // ! MLN_FUN_VV2V_SPAN_HH diff --git a/milena/mln/fun/vvvv2v/mean.hh b/milena/mln/fun/vvvv2v/max.hh similarity index 79% copy from milena/mln/fun/vvvv2v/mean.hh copy to milena/mln/fun/vvvv2v/max.hh index 0dba050..311bcfe 100644 --- a/milena/mln/fun/vvvv2v/mean.hh +++ b/milena/mln/fun/vvvv2v/max.hh @@ -23,15 +23,15 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_FUN_VVVV2V_MEAN_HH -# define MLN_FUN_VVVV2V_MEAN_HH +#ifndef MLN_FUN_VVVV2V_MAX_HH +# define MLN_FUN_VVVV2V_MAX_HH /// \file /// -/// Functor that computes the mean of two values. +/// Functor that computes the max of two values. # include <mln/core/concept/function.hh> -# include <mln/math/mean.hh> +# include <mln/math/max.hh> namespace mln @@ -45,9 +45,9 @@ namespace mln // FIXME: Doc. - /// \brief A functor computing the mean of two values. + /// \brief A functor computing the max of two values. template <typename T, typename R=T> - struct mean : public Function_vvvv2v< mean<T> > + struct max : public Function_vvvv2v< max<T> > { typedef R result; R operator()(const T& v1, const T& v2, const T& v3, const T& v4) const; @@ -56,12 +56,12 @@ namespace mln # ifndef MLN_INCLUDE_ONLY - template <typename T> + template <typename T, typename R> inline R - mean<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const + max<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const { - return R(mln::math::mean(v1, v2, v3, v4)); + return R(mln::math::max(v1, v2, v3, v4)); } # endif // ! MLN_INCLUDE_ONLY @@ -73,4 +73,4 @@ namespace mln } // end of namespace mln -#endif // ! MLN_FUN_VVVV2V_MEAN_HH +#endif // ! MLN_FUN_VVVV2V_MAX_HH diff --git a/milena/mln/fun/vvvv2v/mean.hh b/milena/mln/fun/vvvv2v/mean.hh index 0dba050..d9810ab 100644 --- a/milena/mln/fun/vvvv2v/mean.hh +++ b/milena/mln/fun/vvvv2v/mean.hh @@ -56,7 +56,7 @@ namespace mln # ifndef MLN_INCLUDE_ONLY - template <typename T> + template <typename T, typename R> inline R mean<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const diff --git a/milena/mln/fun/vvvv2v/mean.hh b/milena/mln/fun/vvvv2v/min.hh similarity index 79% copy from milena/mln/fun/vvvv2v/mean.hh copy to milena/mln/fun/vvvv2v/min.hh index 0dba050..f0c4b63 100644 --- a/milena/mln/fun/vvvv2v/mean.hh +++ b/milena/mln/fun/vvvv2v/min.hh @@ -23,15 +23,15 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_FUN_VVVV2V_MEAN_HH -# define MLN_FUN_VVVV2V_MEAN_HH +#ifndef MLN_FUN_VVVV2V_MIN_HH +# define MLN_FUN_VVVV2V_MIN_HH /// \file /// -/// Functor that computes the mean of two values. +/// Functor that computes the min of two values. # include <mln/core/concept/function.hh> -# include <mln/math/mean.hh> +# include <mln/math/min.hh> namespace mln @@ -45,9 +45,9 @@ namespace mln // FIXME: Doc. - /// \brief A functor computing the mean of two values. + /// \brief A functor computing the min of two values. template <typename T, typename R=T> - struct mean : public Function_vvvv2v< mean<T> > + struct min : public Function_vvvv2v< min<T> > { typedef R result; R operator()(const T& v1, const T& v2, const T& v3, const T& v4) const; @@ -56,12 +56,12 @@ namespace mln # ifndef MLN_INCLUDE_ONLY - template <typename T> + template <typename T, typename R> inline R - mean<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const + min<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const { - return R(mln::math::mean(v1, v2, v3, v4)); + return R(mln::math::min(v1, v2, v3, v4)); } # endif // ! MLN_INCLUDE_ONLY @@ -73,4 +73,4 @@ namespace mln } // end of namespace mln -#endif // ! MLN_FUN_VVVV2V_MEAN_HH +#endif // ! MLN_FUN_VVVV2V_MIN_HH diff --git a/milena/mln/fun/vvvv2v/mean.hh b/milena/mln/fun/vvvv2v/span.hh similarity index 67% copy from milena/mln/fun/vvvv2v/mean.hh copy to milena/mln/fun/vvvv2v/span.hh index 0dba050..8c99bd1 100644 --- a/milena/mln/fun/vvvv2v/mean.hh +++ b/milena/mln/fun/vvvv2v/span.hh @@ -23,15 +23,15 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_FUN_VVVV2V_MEAN_HH -# define MLN_FUN_VVVV2V_MEAN_HH +#ifndef MLN_FUN_VVVV2V_SPAN_HH +# define MLN_FUN_VVVV2V_SPAN_HH /// \file /// -/// Functor that computes the mean of two values. +/// Functor that computes the spanimum of two values. # include <mln/core/concept/function.hh> -# include <mln/math/mean.hh> +# include <mln/value/interval.hh> namespace mln @@ -45,23 +45,29 @@ namespace mln // FIXME: Doc. - /// \brief A functor computing the mean of two values. - template <typename T, typename R=T> - struct mean : public Function_vvvv2v< mean<T> > + /// \brief A functor computing the span of two interval values. + template <typename T> + struct span : public Function_vvvv2v< span<T> > { - typedef R result; - R operator()(const T& v1, const T& v2, const T& v3, const T& v4) const; + typedef value::interval<T> result; + + value::interval<T> operator()(const value::interval<T>& v1, + const value::interval<T>& v2, + const value::interval<T>& v3, + const value::interval<T>& v4) const; }; # ifndef MLN_INCLUDE_ONLY template <typename T> - inline - R - mean<T,R>::operator()(const T& v1, const T& v2, const T& v3, const T& v4) const + value::interval<T> + span<T>::operator()(const value::interval<T>& v1, + const value::interval<T>& v2, + const value::interval<T>& v3, + const value::interval<T>& v4) const { - return R(mln::math::mean(v1, v2, v3, v4)); + return value::span(v1, v2, v3, v4); } # endif // ! MLN_INCLUDE_ONLY @@ -73,4 +79,4 @@ namespace mln } // end of namespace mln -#endif // ! MLN_FUN_VVVV2V_MEAN_HH +#endif // ! MLN_FUN_VVVV2V_SPAN_HH diff --git a/milena/mln/math/max.hh b/milena/mln/math/max.hh index 0fe8ea4..786aa2e 100644 --- a/milena/mln/math/max.hh +++ b/milena/mln/math/max.hh @@ -43,6 +43,9 @@ namespace mln template <typename T> T max(const T& v1, const T& v2); + /// Return the maximum of four values. + template <typename T> + T max(const T& v1, const T& v2, const T& v3, const T& v4); } // end of namespace mln::math @@ -55,6 +58,14 @@ namespace mln return v1 > v2 ? v1 : v2; } + /// \internal Generic implementation of the maximum function. + template <typename T> + T max_(const T& v1, const T& v2, const T& v3, const T& v4) + { + return max_(max_(exact(v1), exact(v2)), max_(exact(v3), exact(v4))); + } + + namespace math { @@ -64,6 +75,12 @@ namespace mln return max_(exact(v1), exact(v2)); } + template <typename T> + T max(const T& v1, const T& v2, const T& v3, const T& v4) + { + return max_(exact(v1), exact(v2), exact(v3), exact(v4)); + } + } // end of namespace mln::math # endif // ! MLN_INCLUDE_ONLY diff --git a/milena/mln/math/min.hh b/milena/mln/math/min.hh index d2a5503..706c103 100644 --- a/milena/mln/math/min.hh +++ b/milena/mln/math/min.hh @@ -43,6 +43,9 @@ namespace mln template <typename T> T min(const T& v1, const T& v2); + /// Return the minimum of four values. + template <typename T> + T min(const T& v1, const T& v2, const T& v3, const T& v4); } // end of namespace mln::math @@ -55,6 +58,14 @@ namespace mln return v1 < v2 ? v1 : v2; } + /// \internal Generic implementation of the minimum function. + template <typename T> + T min_(const T& v1, const T& v2, const T& v3, const T& v4) + { + return min_(min_(v1, v2), min_(v3, v4)); + } + + namespace math { @@ -64,6 +75,12 @@ namespace mln return min_(exact(v1), exact(v2)); } + template <typename T> + T min(const T& v1, const T& v2, const T& v3, const T& v4) + { + return min_(exact(v1), exact(v2), exact(v3), exact(v4)); + } + } // end of namespace mln::math # endif // ! MLN_INCLUDE_ONLY -- 1.7.2.5

1 0

branch fft updated: olena-2.0-80-g70839e9
by git＠lrde.epita.fr 10 Oct '12

10 Oct '12

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "Olena, a generic and efficient image processing platform". The branch fft has been updated via 70839e9dd38681e1a878aef7704f9fe5cea714df (commit) via 6468cd2dfbfd4fd467726fa1ff21a6c8b3f7c6e7 (commit) via 6fc43639af2664bee4bf8c15cd2bc2b064f1df2c (commit) via d031727ef459696d4c41eea15c5718c8673144e6 (commit) from 9891bf6decdbf37fcf373c0a569e2972943b0eb2 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- 70839e9 Split the interface of the FFT off from its implementation. 6468cd2 Stylistic changes in the Fast Fourier Transform. 6fc4363 Remove dead code in the Fast Fourier Transform. d031727 Minor corrections in the Fast Fourier Transform. ----------------------------------------------------------------------- Summary of changes: milena/ChangeLog | 31 ++ milena/mln/transform/fft.hh | 1205 ++++++++++++++++++++++--------------------- 2 files changed, 640 insertions(+), 596 deletions(-) hooks/post-receive -- Olena, a generic and efficient image processing platform

1 0

olena-2.0-80-g70839e9 Split the interface of the FFT off from its implementation.
by Roland Levillain 10 Oct '12

10 Oct '12

* 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 faa5f7b..4fc30ca 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.2.5

1 0

olena-2.0-79-g6468cd2 Stylistic changes in the Fast Fourier Transform.
by Roland Levillain 10 Oct '12

10 Oct '12

* 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 55c1dcc..faa5f7b 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.2.5

1 0

olena-2.0-78-g6fc4363 Remove dead code in the Fast Fourier Transform.
by Roland Levillain 10 Oct '12

10 Oct '12

* 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 9a9e6c4..55c1dcc 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.2.5

1 0

olena-2.0-77-gd031727 Minor corrections in the Fast Fourier Transform.
by Roland Levillain 10 Oct '12

10 Oct '12

* 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 bd06a16..9a9e6c4 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.2.5

1 0

branch fft created: olena-2.0-76-g9891bf6
by git＠lrde.epita.fr 09 Oct '12

09 Oct '12

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "Olena, a generic and efficient image processing platform". The branch fft has been created at 9891bf6decdbf37fcf373c0a569e2972943b0eb2 (commit) - Log ----------------------------------------------------------------- 9891bf6 Improve the test of the Fast Fourier Transform. a8e615d Update copyright headers. a5f7553 headers.mk, tests/unit_test/unit-tests.mk: Regen. 9be8364 Exercise the Fast Fourier Transform. 81f0816 Fix the Fast Fourier Transform. 118fb8a Ask configure to try to find FFTW (version 3). 2d7f714 Import the FFTW-based Fast Fourier Transform. ----------------------------------------------------------------------- hooks/post-receive -- Olena, a generic and efficient image processing platform

1 0

olena-2.0-76-g9891bf6 Improve the test of the Fast Fourier Transform.
by Roland Levillain 09 Oct '12

09 Oct '12

* 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 037656f..bd06a16 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 f458a92..2eaad1f 100644 --- a/milena/tests/transform/Makefile.am +++ b/milena/tests/transform/Makefile.am @@ -54,4 +54,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.2.5

1 0

olena-2.0-75-ga8e615d Update copyright headers.
by Roland Levillain 09 Oct '12

09 Oct '12

* 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 29ead4f..037656f 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.2.5

1 0

olena-2.0-74-ga5f7553 headers.mk, tests/unit_test/unit-tests.mk: Regen.
by Roland Levillain 09 Oct '12

09 Oct '12

--- milena/ChangeLog | 4 ++++ milena/headers.mk | 1 + milena/tests/unit_test/unit-tests.mk | 2 ++ 3 files changed, 7 insertions(+), 0 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 88f2402..29ead4f 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,9 @@ 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> + Exercise the Fast Fourier Transform. * tests/transform/fft.cc: Fix the path to the input image. diff --git a/milena/headers.mk b/milena/headers.mk index dbd26e9..2397af6 100644 --- a/milena/headers.mk +++ b/milena/headers.mk @@ -1134,6 +1134,7 @@ mln/transform/distance_and_influence_zone_geodesic.hh \ mln/transform/distance_front.hh \ mln/transform/distance_geodesic.hh \ mln/transform/essential.hh \ +mln/transform/fft.hh \ mln/transform/hough.hh \ mln/transform/influence_zone_front.hh \ mln/transform/influence_zone_geodesic.hh \ diff --git a/milena/tests/unit_test/unit-tests.mk b/milena/tests/unit_test/unit-tests.mk index 4fc2dad..f6c772a 100644 --- a/milena/tests/unit_test/unit-tests.mk +++ b/milena/tests/unit_test/unit-tests.mk @@ -1176,6 +1176,7 @@ mln_transform_distance_and_influence_zone_geodesic \ mln_transform_distance_front \ mln_transform_distance_geodesic \ mln_transform_essential \ +mln_transform_fft \ mln_transform_hough \ mln_transform_influence_zone_front \ mln_transform_influence_zone_geodesic \ @@ -2440,6 +2441,7 @@ mln_transform_distance_and_influence_zone_geodesic_SOURCES = mln_transform_dista mln_transform_distance_front_SOURCES = mln_transform_distance_front.cc mln_transform_distance_geodesic_SOURCES = mln_transform_distance_geodesic.cc mln_transform_essential_SOURCES = mln_transform_essential.cc +mln_transform_fft_SOURCES = mln_transform_fft.cc mln_transform_hough_SOURCES = mln_transform_hough.cc mln_transform_influence_zone_front_SOURCES = mln_transform_influence_zone_front.cc mln_transform_influence_zone_geodesic_SOURCES = mln_transform_influence_zone_geodesic.cc -- 1.7.2.5

1 0