* 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@lrde.epita.fr

+ Split the interface of the FFT off from its implementation. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain roland@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